Line data Source code
1 : /**************************************************************************
2 : * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 : * *
4 : * Author: The ALICE Off-line Project. *
5 : * Contributors are mentioned in the code where appropriate. *
6 : * *
7 : * Permission to use, copy, modify and distribute this software and its *
8 : * documentation strictly for non-commercial purposes is hereby granted *
9 : * without fee, provided that the above copyright notice appears in all *
10 : * copies and that both the copyright notice and this permission notice *
11 : * appear in the supporting documentation. The authors make no claims *
12 : * about the suitability of this software for any purpose. It is *
13 : * provided "as is" without express or implied warranty. *
14 : **************************************************************************/
15 :
16 : /* $Id: AliTPCclustererKr.cxx,v 1.7 2008/02/06 17:24:53 matyja Exp $ */
17 :
18 : //-----------------------------------------------------------------
19 : // Implementation of the TPC Kr cluster class
20 : //
21 : // Origin: Adam Matyja, INP PAN, adam.matyja@ifj.edu.pl
22 : //-----------------------------------------------------------------
23 :
24 : /*
25 : Instruction - how to use that:
26 : There are two macros prepared. One is for preparing clusters from MC
27 : samples: FindKrClusters.C. The output is kept in TPC.RecPoints.root.
28 : The other macro is prepared for data analysis: FindKrClustersRaw.C.
29 : The output is created for each processed file in root file named adc.root.
30 : For each data subsample the same named file is created. So be careful
31 : do not overwrite them.
32 :
33 : Additional selection criteria to select the GOLD cluster
34 : Example:
35 : // open file with clusters
36 : TFile f("Krypton.root");
37 : TTree * tree = (TTree*)f.Get("Kr")
38 : TCut cutR0("cutR0","fADCcluster/fSize<100"); // adjust it according v seetings -
39 : TCut cutR1("cutR1","fADCcluster/fSize>7"); // cosmic tracks and noise removal
40 : TCut cutR2("cutR2","fMax.fAdc/fADCcluster<0.2"); // digital noise removal
41 : TCut cutR3("cutR3","fMax.fAdc/fADCcluster>0.01"); // noise removal
42 : TCut cutS1("cutS1","fSize<200"); // adjust it according v seetings - cosmic tracks
43 : TCut cutAll = cutR0+cutR1+cutR2+cutR3+cutS1;
44 : This values are typical values to be applied in selectors
45 :
46 :
47 : *
48 : **** MC ****
49 : *
50 :
51 : To run clusterizaton for MC type:
52 : .x FindKrClusters.C
53 :
54 : If you don't want to use the standard selection criteria then you
55 : have to do following:
56 :
57 : // load the standard setup
58 : AliRunLoader* rl = AliRunLoader::Open("galice.root");
59 : AliTPCLoader *tpcl = (AliTPCLoader*)rl->GetLoader("TPCLoader");
60 : tpcl->LoadDigits();
61 : rl->LoadgAlice();
62 : gAlice=rl->GetAliRun();
63 : TDirectory *cwd = gDirectory;
64 : AliTPCv4 *tpc = (AliTPCv4*)gAlice->GetDetector("TPC");
65 : Int_t ver = tpc->IsVersion();
66 : rl->CdGAFile();
67 : AliTPCParam *param=(AliTPCParamSR *)gDirectory->Get("75x40_100x60_150x60");
68 : AliTPCDigitsArray *digarr=new AliTPCDigitsArray;
69 : digarr->Setup(param);
70 : cwd->cd();
71 :
72 : //loop over events
73 : Int_t nevmax=rl->GetNumberOfEvents();
74 : for(Int_t nev=0;nev<nevmax ;nev++){
75 : rl->GetEvent(nev);
76 : TTree* input_tree= tpcl->TreeD();//load tree with digits
77 : digarr->ConnectTree(input_tree);
78 : TTree *output_tree =tpcl->TreeR();//load output tree
79 :
80 : AliTPCclustererKr *clusters = new AliTPCclustererKr();
81 : clusters->SetParam(param);
82 : clusters->SetInput(input_tree);
83 : clusters->SetOutput(output_tree);
84 : clusters->SetDigArr(digarr);
85 :
86 : //If you want to change the cluster finder parameters for MC there are
87 : //several of them:
88 :
89 : //1. signal threshold (everything below the given number is treated as 0)
90 : clusters->SetMinAdc(3);
91 :
92 : //2. number of neighbouring timebins to be considered
93 : clusters->SetMinTimeBins(2);
94 :
95 : //3. distance of the cluster center to the center of a pad in pad-padrow plane
96 : //(in cm). Remenber that this is still quantified by pad size.
97 : clusters->SetMaxPadRangeCm(2.5);
98 :
99 : //4. distance of the cluster center to the center of a padrow in pad-padrow
100 : //plane (in cm). Remenber that this is still quantified by pad size.
101 : clusters->SetMaxRowRangeCm(3.5);
102 :
103 : //5. distance of the cluster center to the max time bin on a pad (in tackts)
104 : //ie. fabs(centerT - time)<7
105 : clusters->SetMaxTimeRange(7);
106 :
107 : //6. cut reduce peak at 0. There are noises which appear mostly as two
108 : //timebins on one pad.
109 : clusters->SetValueToSize(3.5);
110 :
111 :
112 : clusters->finderIO();
113 : tpcl->WriteRecPoints("OVERWRITE");
114 : }
115 : delete rl;//cleans everything
116 :
117 : *
118 : ********* DATA *********
119 : *
120 :
121 : To run clusterizaton for DATA for file named raw_data.root type:
122 : .x FindKrClustersRaw.C("raw_data.root")
123 :
124 : If you want to change some criteria do the following:
125 :
126 : //
127 : // remove Altro warnings
128 : //
129 : AliLog::SetClassDebugLevel("AliRawReaderDate",-5);
130 : AliLog::SetClassDebugLevel("AliTPCAltroMapping",-5);
131 : AliLog::SetModuleDebugLevel("RAW",-5);
132 :
133 : //
134 : // Get database with noises
135 : //
136 : // char *ocdbpath = gSystem->Getenv("OCDB_PATH");
137 : char *ocdbpath ="local:///afs/cern.ch/alice/tpctest/OCDB";
138 : if (ocdbpath==0){
139 : ocdbpath="alien://folder=/alice/data/2007/LHC07w/OCDB/";
140 : }
141 : AliCDBManager * man = AliCDBManager::Instance();
142 : man->SetDefaultStorage(ocdbpath);
143 : man->SetRun(0);
144 : AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
145 : AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping();
146 :
147 : //define tree
148 : TFile *hfile=new TFile("adc.root","RECREATE","ADC file");
149 : // Create a ROOT Tree
150 : TTree *mytree = new TTree("Kr","Krypton cluster tree");
151 :
152 : //define infput file
153 : const char *fileName="data.root";
154 : AliRawReader *reader = new AliRawReaderRoot(fileName);
155 : //AliRawReader *reader = new AliRawReaderDate(fileName);
156 : reader->Reset();
157 : AliAltroRawStreamFast* stream = new AliAltroRawStreamFast(reader);
158 : stream->SelectRawData("TPC");
159 :
160 : //one general output
161 : AliTPCclustererKr *clusters = new AliTPCclustererKr();
162 : clusters->SetOutput(mytree);
163 : clusters->SetRecoParam(0);//standard reco parameters
164 : AliTPCParamSR *param=new AliTPCParamSR();
165 : clusters->SetParam(param);//TPC parameters(sectors, timebins, etc.)
166 :
167 : //set cluster finder parameters (from data):
168 : //1. zero suppression parameter
169 : clusters->SetZeroSup(param->GetZeroSup());
170 :
171 : //2. first bin
172 : clusters->SetFirstBin(60);
173 :
174 : //3. last bin
175 : clusters->SetLastBin(950);
176 :
177 : //4. maximal noise
178 : clusters->SetMaxNoiseAbs(2);
179 :
180 : //5. maximal amount of sigma of noise
181 : clusters->SetMaxNoiseSigma(3);
182 :
183 : //The remaining parameters are the same paramters as for MC (see MC section
184 : //points 1-6)
185 : clusters->SetMinAdc(3);
186 : clusters->SetMinTimeBins(2);
187 : clusters->SetMaxPadRangeCm(2.5);
188 : clusters->SetMaxRowRangeCm(3.5);
189 : clusters->SetMaxTimeRange(7);
190 : clusters->SetValueToSize(3.5);
191 :
192 : while (reader->NextEvent()) {
193 : clusters->FinderIO(reader);
194 : }
195 :
196 : hfile->Write();
197 : hfile->Close();
198 : delete stream;
199 :
200 :
201 : */
202 :
203 : #include "AliTPCclustererKr.h"
204 : #include "AliTPCclusterKr.h"
205 : //#include <vector>
206 : #include <list>
207 : #include "TObject.h"
208 : #include "AliPadMax.h"
209 : #include "AliSimDigits.h"
210 : //#include "AliTPCv4.h"
211 : #include "AliTPCParam.h"
212 : #include "AliTPCDigitsArray.h"
213 : #include "AliTPCvtpr.h"
214 : #include "AliTPCClustersRow.h"
215 : #include "TTree.h"
216 : #include "TH1F.h"
217 : #include "TH2F.h"
218 : #include "TTreeStream.h"
219 :
220 : #include "AliTPCTransform.h"
221 :
222 : //used in raw data finder
223 : #include "AliTPCROC.h"
224 : #include "AliTPCCalPad.h"
225 : #include "AliTPCAltroMapping.h"
226 : #include "AliTPCcalibDB.h"
227 : #include "AliTPCRawStreamV3.h"
228 : #include "AliTPCRecoParam.h"
229 : #include "AliTPCReconstructor.h"
230 : #include "AliRawReader.h"
231 : #include "AliTPCCalROC.h"
232 : #include "AliRawEventHeaderBase.h"
233 :
234 : using std::cerr;
235 : using std::cout;
236 : using std::endl;
237 : using std::list;
238 16 : ClassImp(AliTPCclustererKr)
239 :
240 :
241 : AliTPCclustererKr::AliTPCclustererKr()
242 0 : :TObject(),
243 0 : fRawData(kFALSE),
244 0 : fInput(0),
245 0 : fOutput(0),
246 0 : fParam(0),
247 0 : fDigarr(0),
248 0 : fRecoParam(0),
249 0 : fZeroSup(2),
250 0 : fFirstBin(60),
251 0 : fLastBin(950),
252 0 : fMaxNoiseAbs(2),
253 0 : fMaxNoiseSigma(3),
254 0 : fMinAdc(3),
255 0 : fMinTimeBins(2),
256 : // fMaxPadRange(4),
257 : // fMaxRowRange(3),
258 0 : fMaxTimeRange(7),
259 0 : fValueToSize(3.5),
260 0 : fMaxPadRangeCm(2.5),
261 0 : fMaxRowRangeCm(3.5),
262 0 : fIsolCut(3),
263 0 : fDebugLevel(-1),
264 0 : fHistoRow(0),
265 0 : fHistoPad(0),
266 0 : fHistoTime(0),
267 0 : fHistoRowPad(0),
268 0 : fTimeStamp(0),
269 0 : fRun(0)
270 0 : {
271 : //
272 : // default constructor
273 : //
274 0 : }
275 :
276 : AliTPCclustererKr::AliTPCclustererKr(const AliTPCclustererKr ¶m)
277 0 : :TObject(),
278 0 : fRawData(kFALSE),
279 0 : fInput(0),
280 0 : fOutput(0),
281 0 : fParam(0),
282 0 : fDigarr(0),
283 0 : fRecoParam(0),
284 0 : fZeroSup(2),
285 0 : fFirstBin(60),
286 0 : fLastBin(950),
287 0 : fMaxNoiseAbs(2),
288 0 : fMaxNoiseSigma(3),
289 0 : fMinAdc(3),
290 0 : fMinTimeBins(2),
291 : // fMaxPadRange(4),
292 : // fMaxRowRange(3),
293 0 : fMaxTimeRange(7),
294 0 : fValueToSize(3.5),
295 0 : fMaxPadRangeCm(2.5),
296 0 : fMaxRowRangeCm(3.5),
297 0 : fIsolCut(3),
298 0 : fDebugLevel(-1),
299 0 : fHistoRow(0),
300 0 : fHistoPad(0),
301 0 : fHistoTime(0),
302 0 : fHistoRowPad(0),
303 0 : fTimeStamp(0),
304 0 : fRun(0)
305 0 : {
306 : //
307 : // copy constructor
308 : //
309 0 : fParam = param.fParam;
310 0 : fRecoParam = param.fRecoParam;
311 0 : fRawData = param.fRawData;
312 0 : fInput = param.fInput ;
313 0 : fOutput = param.fOutput;
314 0 : fDigarr = param.fDigarr;
315 0 : fZeroSup = param.fZeroSup ;
316 0 : fFirstBin = param.fFirstBin ;
317 0 : fLastBin = param.fLastBin ;
318 0 : fMaxNoiseAbs = param.fMaxNoiseAbs ;
319 0 : fMaxNoiseSigma = param.fMaxNoiseSigma ;
320 0 : fMinAdc = param.fMinAdc;
321 0 : fMinTimeBins = param.fMinTimeBins;
322 : // fMaxPadRange = param.fMaxPadRange ;
323 : // fMaxRowRange = param.fMaxRowRange ;
324 0 : fMaxTimeRange = param.fMaxTimeRange;
325 0 : fValueToSize = param.fValueToSize;
326 0 : fMaxPadRangeCm = param.fMaxPadRangeCm;
327 0 : fMaxRowRangeCm = param.fMaxRowRangeCm;
328 0 : fIsolCut = param.fIsolCut;
329 0 : fDebugLevel = param.fDebugLevel;
330 0 : fHistoRow = param.fHistoRow ;
331 0 : fHistoPad = param.fHistoPad ;
332 0 : fHistoTime = param.fHistoTime;
333 0 : fHistoRowPad = param.fHistoRowPad;
334 0 : fTimeStamp = param.fTimeStamp;
335 0 : fRun = param.fRun;
336 :
337 0 : }
338 :
339 : AliTPCclustererKr & AliTPCclustererKr::operator = (const AliTPCclustererKr & param)
340 : {
341 : //
342 : // assignment operator
343 : //
344 0 : if (this == ¶m) return (*this);
345 :
346 0 : fParam = param.fParam;
347 0 : fRecoParam = param.fRecoParam;
348 0 : fRawData = param.fRawData;
349 0 : fInput = param.fInput ;
350 0 : fOutput = param.fOutput;
351 0 : fDigarr = param.fDigarr;
352 0 : fZeroSup = param.fZeroSup ;
353 0 : fFirstBin = param.fFirstBin ;
354 0 : fLastBin = param.fLastBin ;
355 0 : fMaxNoiseAbs = param.fMaxNoiseAbs ;
356 0 : fMaxNoiseSigma = param.fMaxNoiseSigma ;
357 0 : fMinAdc = param.fMinAdc;
358 0 : fMinTimeBins = param.fMinTimeBins;
359 : // fMaxPadRange = param.fMaxPadRange ;
360 : // fMaxRowRange = param.fMaxRowRange ;
361 0 : fMaxTimeRange = param.fMaxTimeRange;
362 0 : fValueToSize = param.fValueToSize;
363 0 : fMaxPadRangeCm = param.fMaxPadRangeCm;
364 0 : fMaxRowRangeCm = param.fMaxRowRangeCm;
365 0 : fIsolCut = param.fIsolCut;
366 0 : fDebugLevel = param.fDebugLevel;
367 0 : fHistoRow = param.fHistoRow ;
368 0 : fHistoPad = param.fHistoPad ;
369 0 : fHistoTime = param.fHistoTime;
370 0 : fHistoRowPad = param.fHistoRowPad;
371 0 : fTimeStamp = param.fTimeStamp;
372 0 : fRun = param.fRun;
373 0 : return (*this);
374 0 : }
375 :
376 : AliTPCclustererKr::~AliTPCclustererKr()
377 0 : {
378 : //
379 : // destructor
380 : //
381 0 : delete fOutput;
382 0 : }
383 :
384 : void AliTPCclustererKr::SetRecoParam(AliTPCRecoParam *recoParam)
385 : {
386 : //
387 : // set reconstruction parameters
388 : //
389 0 : if (recoParam) {
390 0 : fRecoParam = recoParam;
391 0 : }else{
392 : //set default parameters if not specified
393 0 : fRecoParam = AliTPCReconstructor::GetRecoParam();
394 0 : if (!fRecoParam) fRecoParam = AliTPCRecoParam::GetLowFluxParam();
395 : }
396 0 : return;
397 : }
398 :
399 :
400 : ////____________________________________________________________________________
401 : //// I/O
402 : void AliTPCclustererKr::SetInput(TTree * tree)
403 : {
404 : //
405 : // set input tree with digits
406 : //
407 0 : fInput = tree;
408 0 : if (!fInput->GetBranch("Segment")){
409 0 : cerr<<"AliTPCclusterKr::FindClusterKr(): no proper input tree !\n";
410 0 : fInput=0;
411 0 : return;
412 : }
413 0 : }
414 :
415 : void AliTPCclustererKr::SetOutput(TTree * /*tree*/)
416 : {
417 : //
418 : //dummy
419 : //
420 0 : fOutput = new TTreeSRedirector("Krypton.root");
421 0 : }
422 :
423 : ////____________________________________________________________________________
424 : //// with new I/O
425 : Int_t AliTPCclustererKr::FinderIO()
426 : {
427 : // Krypton cluster finder for simulated events from MC
428 :
429 0 : if (!fInput) {
430 0 : Error("Digits2Clusters", "input tree not initialised");
431 0 : return 10;
432 : }
433 :
434 0 : if (!fOutput) {
435 0 : Error("Digits2Clusters", "output tree not initialised");
436 0 : return 11;
437 : }
438 :
439 0 : FindClusterKrIO();
440 0 : return 0;
441 0 : }
442 :
443 :
444 :
445 : Int_t AliTPCclustererKr::FinderIO(AliRawReader* rawReader)
446 : {
447 : // Krypton cluster finder for the TPC raw data
448 : // this method is unsing AliAltroRawStreamV3
449 : // fParam must be defined before
450 0 : if (!rawReader) return 1;
451 : //
452 0 : fRawData=kTRUE; //set flag to data
453 :
454 0 : if (!fOutput) {
455 0 : Error("Digits2Clusters", "output tree not initialised");
456 0 : return 11;
457 : }
458 :
459 0 : fParam->SetMaxTBin(fRecoParam->GetLastBin());//set number of timebins from reco -> param
460 : // used later for memory allocation
461 :
462 0 : AliRawEventHeaderBase* eventHeader = (AliRawEventHeaderBase*)rawReader->GetEventHeader();
463 0 : if (eventHeader){
464 0 : fTimeStamp = eventHeader->Get("Timestamp");
465 0 : fRun = rawReader->GetRunNumber();
466 0 : }
467 :
468 :
469 : Bool_t isAltro=kFALSE;
470 :
471 0 : AliTPCROC * roc = AliTPCROC::Instance();
472 0 : AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
473 0 : AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping();
474 : //
475 0 : AliTPCRawStreamV3 input(rawReader,(AliAltroMapping**)mapping);
476 :
477 0 : const Int_t kNIS = fParam->GetNInnerSector();//number of inner sectors
478 0 : const Int_t kNOS = fParam->GetNOuterSector();//number of outer sectors
479 0 : const Int_t kNS = kNIS + kNOS;//all sectors
480 :
481 :
482 : //crate TPC view
483 0 : AliTPCDigitsArray *digarr=new AliTPCDigitsArray(kFALSE);//data not sim
484 0 : digarr->Setup(fParam);//as usually parameters
485 :
486 0 : for(Int_t iSec = 0; iSec < kNS; iSec++) {
487 : AliTPCCalROC * noiseROC;
488 0 : AliTPCCalROC noiseDummy(iSec);
489 0 : if(noiseTPC==0x0){
490 : noiseROC = &noiseDummy;//noise=0
491 0 : }else{
492 0 : noiseROC = noiseTPC->GetCalROC(iSec); // noise per given sector
493 : }
494 : Int_t nRows = 0; //number of rows in sector
495 : Int_t nDDLs = 0; //number of DDLs
496 : Int_t indexDDL = 0; //DDL index
497 0 : if (iSec < kNIS) {
498 0 : nRows = fParam->GetNRowLow();
499 : nDDLs = 2;
500 0 : indexDDL = iSec * 2;
501 0 : }else {
502 0 : nRows = fParam->GetNRowUp();
503 : nDDLs = 4;
504 0 : indexDDL = (iSec-kNIS) * 4 + kNIS * 2;
505 : }
506 :
507 : //
508 : // Load the raw data for corresponding DDLs
509 : //
510 0 : rawReader->Reset();
511 0 : rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1);
512 :
513 :
514 0 : while (input.NextDDL()){
515 : // Allocate memory for rows in sector (pads(depends on row) x timebins)
516 0 : if (!digarr->GetRow(iSec,0)){
517 0 : for(Int_t iRow = 0; iRow < nRows; iRow++) {
518 0 : digarr->CreateRow(iSec,iRow);
519 : }//end loop over rows
520 0 : }
521 : //loop over pads
522 0 : while ( input.NextChannel() ) {
523 0 : Int_t iRow = input.GetRow();
524 0 : Int_t iPad = input.GetPad();
525 : //check row consistency
526 0 : if (iRow < 0 ) continue;
527 0 : if (iRow < 0 || iRow >= nRows){
528 0 : AliError(Form("Pad-row index (%d) outside the range (%d -> %d) !",
529 : iRow, 0, nRows -1));
530 0 : continue;
531 : }
532 :
533 : //check pad consistency
534 0 : if (iPad < 0 || iPad >= (Int_t)(roc->GetNPads(iSec,iRow))) {
535 0 : AliError(Form("Pad index (%d) outside the range (%d -> %d) !",
536 : iPad, 0, roc->GetNPads(iSec,iRow) ));
537 0 : continue;
538 : }
539 :
540 : //loop over bunches
541 0 : while ( input.NextBunch() ){
542 0 : Int_t startTbin = (Int_t)input.GetStartTimeBin();
543 0 : Int_t bunchlength = (Int_t)input.GetBunchLength();
544 0 : const UShort_t *sig = input.GetSignals();
545 : isAltro=kTRUE;
546 0 : for (Int_t iTime = 0; iTime<bunchlength; iTime++){
547 0 : Int_t iTimeBin=startTbin-iTime;
548 : //
549 0 : if(fDebugLevel==72){
550 0 : fHistoRow->Fill(iRow);
551 0 : fHistoPad->Fill(iPad);
552 0 : fHistoTime->Fill(iTimeBin);
553 0 : fHistoRowPad->Fill(iPad,iRow);
554 0 : }else if(fDebugLevel>=0&&fDebugLevel<72){
555 0 : if(iSec==fDebugLevel){
556 0 : fHistoRow->Fill(iRow);
557 0 : fHistoPad->Fill(iPad);
558 0 : fHistoTime->Fill(iTimeBin);
559 0 : fHistoRowPad->Fill(iPad,iRow);
560 : }
561 0 : }else if(fDebugLevel==73){
562 0 : if(iSec<36){
563 0 : fHistoRow->Fill(iRow);
564 0 : fHistoPad->Fill(iPad);
565 0 : fHistoTime->Fill(iTimeBin);
566 0 : fHistoRowPad->Fill(iPad,iRow);
567 : }
568 0 : }else if(fDebugLevel==74){
569 0 : if(iSec>=36){
570 0 : fHistoRow->Fill(iRow);
571 0 : fHistoPad->Fill(iPad);
572 0 : fHistoTime->Fill(iTimeBin);
573 0 : fHistoRowPad->Fill(iPad,iRow);
574 : }
575 : }
576 :
577 : //check time consistency
578 0 : if ( iTimeBin < fRecoParam->GetFirstBin() || iTimeBin >= fRecoParam->GetLastBin()){
579 : //cout<<iTimeBin<<endl;
580 0 : continue;
581 : AliFatal(Form("Timebin index (%d) outside the range (%d -> %d) !",
582 : iTimeBin, 0, fRecoParam->GetLastBin() -1));
583 : }
584 : //signal
585 0 : Float_t signal=(Float_t)sig[iTime];
586 0 : if (signal <= fZeroSup ||
587 0 : iTimeBin < fFirstBin ||
588 0 : iTimeBin > fLastBin
589 : ) {
590 0 : digarr->GetRow(iSec,iRow)->SetDigitFast(0,iTimeBin,iPad);
591 0 : continue;
592 : }
593 0 : if (!noiseROC) continue;
594 0 : Double_t noiseOnPad = noiseROC->GetValue(iRow,iPad);//noise on given pad and row in sector
595 0 : if (noiseOnPad > fMaxNoiseAbs){
596 0 : digarr->GetRow(iSec,iRow)->SetDigitFast(0,iTimeBin,iPad);
597 0 : continue; // consider noisy pad as dead
598 : }
599 0 : if(signal <= fMaxNoiseSigma * noiseOnPad){
600 0 : digarr->GetRow(iSec,iRow)->SetDigitFast(0,iTimeBin,iPad);
601 0 : continue;
602 : }
603 0 : digarr->GetRow(iSec,iRow)->SetDigitFast(TMath::Nint(signal),iTimeBin,iPad);
604 0 : }// end loop signals in bunch
605 : }// end loop bunches
606 0 : } // end loop pads
607 : }// end ddl loop
608 0 : }// end sector loop
609 0 : SetDigArr(digarr);
610 0 : if(isAltro) FindClusterKrIO();
611 0 : delete digarr;
612 :
613 : return 0;
614 0 : }
615 :
616 : void AliTPCclustererKr::CleanSector(Int_t sector){
617 : //
618 : // clean isolated digits
619 : //
620 0 : const Int_t kNRows=fParam->GetNRow(sector);//number of rows in sector
621 0 : for(Int_t iRow=0; iRow<kNRows; ++iRow){
622 : AliSimDigits *digrow;
623 0 : if(fRawData){
624 0 : digrow = (AliSimDigits*)fDigarr->GetRow(sector,iRow);//real data
625 0 : }else{
626 0 : digrow = (AliSimDigits*)fDigarr->LoadRow(sector,iRow);//MC
627 : }
628 0 : if(!digrow) continue;
629 0 : digrow->ExpandBuffer(); //decrunch
630 0 : const Int_t kNPads = digrow->GetNCols(); // number of pads
631 0 : const Int_t kNTime = digrow->GetNRows(); // number of timebins
632 0 : for(Int_t iPad=1;iPad<kNPads-1;iPad++){
633 0 : Short_t* val = digrow->GetDigitsColumn(iPad);
634 :
635 0 : for(Int_t iTimeBin=1;iTimeBin<kNTime-1;iTimeBin++){
636 0 : if (val[iTimeBin]<=0) continue;
637 0 : if (val[iTimeBin-1]+val[iTimeBin+1]<fIsolCut) {val[iTimeBin]=0; continue;}
638 0 : if (val[iTimeBin-kNTime]+val[iTimeBin+kNTime]<fIsolCut) {val[iTimeBin]=0; continue;}
639 : //
640 0 : if (val[iTimeBin-1-kNTime]+val[iTimeBin+1+kNTime]<fIsolCut) {val[iTimeBin]=0; continue;}
641 0 : if (val[iTimeBin+1-kNTime]+val[iTimeBin-1+kNTime]<fIsolCut) {val[iTimeBin]=0; continue;}
642 :
643 : }
644 : }
645 0 : }
646 0 : }
647 :
648 :
649 : ////____________________________________________________________________________
650 : Int_t AliTPCclustererKr::FindClusterKrIO()
651 : {
652 :
653 : //
654 : //fParam and fDigarr must be set to run this method
655 : //
656 :
657 0 : Int_t clusterCounter=0;
658 0 : const Int_t nTotalSector=fParam->GetNSector();//number of sectors
659 0 : for(Int_t iSec=0; iSec<nTotalSector; ++iSec){
660 0 : CleanSector(iSec);
661 :
662 : //vector of maxima for each sector
663 : //std::vector<AliPadMax*> maximaInSector;
664 0 : TObjArray *maximaInSector=new TObjArray();//to store AliPadMax*
665 :
666 : //
667 : // looking for the maxima on the pad
668 : //
669 :
670 0 : const Int_t kNRows=fParam->GetNRow(iSec);//number of rows in sector
671 0 : for(Int_t iRow=0; iRow<kNRows; ++iRow){
672 : AliSimDigits *digrow;
673 0 : if(fRawData){
674 0 : digrow = (AliSimDigits*)fDigarr->GetRow(iSec,iRow);//real data
675 0 : }else{
676 0 : digrow = (AliSimDigits*)fDigarr->LoadRow(iSec,iRow);//MC
677 : }
678 0 : if(digrow){//if pointer exist
679 0 : digrow->ExpandBuffer(); //decrunch
680 0 : const Int_t kNPads = digrow->GetNCols(); // number of pads
681 0 : const Int_t kNTime = digrow->GetNRows(); // number of timebins
682 0 : for(Int_t iPad=0;iPad<kNPads;iPad++){
683 :
684 : Int_t timeBinMax=-1;//timebin of maximum
685 : Int_t valueMaximum=-1;//value of maximum in adc
686 : Int_t increaseBegin=-1;//timebin when increase starts
687 : Int_t sumAdc=0;//sum of adc on the pad in maximum surrounding
688 : bool ifIncreaseBegin=true;//flag - check if increasing started
689 : bool ifMaximum=false;//flag - check if it could be maximum
690 0 : Short_t* val = digrow->GetDigitsColumn(iPad);
691 0 : for(Int_t iTimeBin=1;iTimeBin<kNTime-1;iTimeBin++){
692 0 : if (!ifMaximum) {
693 0 : if (val[iTimeBin]==-1) break; // 0 until the end
694 0 : for( ; iTimeBin<kNTime-2&&val[iTimeBin]<fMinAdc ;iTimeBin++) {}
695 : }
696 : //
697 0 : Short_t adc = val[iTimeBin];
698 :
699 0 : if(adc<fMinAdc){//standard was 3 for fMinAdc
700 0 : if(ifMaximum){
701 0 : if(iTimeBin-increaseBegin<fMinTimeBins){//at least 2 time bins
702 : timeBinMax=-1;
703 : valueMaximum=-1;
704 : increaseBegin=-1;
705 : sumAdc=0;
706 : ifIncreaseBegin=true;
707 : ifMaximum=false;
708 0 : continue;
709 : }
710 : //insert maximum, default values and set flags
711 : //Double_t xCord,yCord;
712 : //GetXY(iSec,iRow,iPad,xCord,yCord);
713 0 : Double_t x[]={static_cast<Double_t>(iRow),static_cast<Double_t>(iPad),static_cast<Double_t>(iTimeBin)};
714 0 : Int_t i[]={iSec};
715 0 : AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
716 :
717 0 : transform->Transform(x,i,0,1);
718 :
719 0 : AliPadMax *oneMaximum = new AliPadMax(AliTPCvtpr(valueMaximum,
720 0 : timeBinMax,
721 0 : iPad,
722 0 : iRow,
723 0 : x[0],//xCord,
724 0 : x[1],//yCord,
725 0 : x[2]/*timeBinMax*/),
726 0 : increaseBegin,
727 0 : iTimeBin-1,
728 0 : sumAdc);
729 0 : maximaInSector->AddLast(oneMaximum);
730 :
731 : timeBinMax=-1;
732 : valueMaximum=-1;
733 : increaseBegin=-1;
734 : sumAdc=0;
735 : ifIncreaseBegin=true;
736 : ifMaximum=false;
737 0 : }
738 0 : continue;
739 : }
740 :
741 :
742 :
743 :
744 :
745 :
746 0 : if(ifIncreaseBegin){
747 : ifIncreaseBegin=false;
748 : increaseBegin=iTimeBin;
749 0 : }
750 :
751 0 : if(adc>valueMaximum){
752 : timeBinMax=iTimeBin;
753 : valueMaximum=adc;
754 : ifMaximum=true;
755 0 : }
756 0 : sumAdc+=adc;
757 0 : if(iTimeBin==kNTime-1 && ifMaximum && kNTime-increaseBegin>fMinTimeBins){//on the edge
758 : //at least 3 timebins
759 : //insert maximum, default values and set flags
760 : //Double_t xCord,yCord;
761 : //GetXY(iSec,iRow,iPad,xCord,yCord);
762 0 : Double_t x[]={static_cast<Double_t>(iRow),static_cast<Double_t>(iPad),static_cast<Double_t>(iTimeBin)};
763 0 : Int_t i[]={iSec};
764 : //AliTPCTransform trafo;
765 : //trafo.Transform(x,i,0,1);
766 :
767 0 : AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
768 :
769 0 : transform->Transform(x,i,0,1);
770 :
771 0 : AliPadMax *oneMaximum = new AliPadMax(AliTPCvtpr(valueMaximum,
772 0 : timeBinMax,
773 0 : iPad,
774 0 : iRow,
775 0 : x[0],//xCord,
776 0 : x[1],//yCord,
777 0 : x[2]/*timeBinMax*/),
778 0 : increaseBegin,
779 0 : iTimeBin-1,
780 0 : sumAdc);
781 0 : maximaInSector->AddLast(oneMaximum);
782 :
783 : timeBinMax=-1;
784 : valueMaximum=-1;
785 : increaseBegin=-1;
786 : sumAdc=0;
787 : ifIncreaseBegin=true;
788 : ifMaximum=false;
789 : continue;
790 0 : }
791 :
792 0 : }//end loop over timebins
793 : }//end loop over pads
794 : // }else{
795 : // cout<<"Pointer does not exist!!"<<endl;
796 0 : }//end if poiner exists
797 : }//end loop over rows
798 :
799 0 : MakeClusters(maximaInSector,iSec,clusterCounter);
800 : //
801 0 : maximaInSector->SetOwner(kTRUE);
802 0 : maximaInSector->Delete();
803 0 : delete maximaInSector;
804 : }//end sector for
805 0 : cout<<"Number of clusters in event: "<<clusterCounter<<endl;
806 0 : return 0;
807 0 : }
808 :
809 : void AliTPCclustererKr::MakeClusters(TObjArray * maximaInSector, Int_t iSec, Int_t &clusterCounter){
810 : //
811 : // Make clusters
812 : //
813 :
814 : Int_t maxDig=0;
815 : Int_t maxSumAdc=0;
816 : Int_t maxTimeBin=0;
817 : Int_t maxPad=0;
818 : Int_t maxRow=0;
819 : Double_t maxX=0;
820 : Double_t maxY=0;
821 : Double_t maxT=0;
822 0 : Int_t entriesArr = maximaInSector->GetEntriesFast();
823 0 : for(Int_t it1 = 0; it1 < entriesArr; ++it1 ) {
824 :
825 0 : AliPadMax *mp1=(AliPadMax *)maximaInSector->UncheckedAt(it1);
826 0 : if (!mp1) continue;
827 0 : AliTPCclusterKr clusterKr;
828 :
829 : Int_t nUsedPads=1;
830 : Int_t clusterValue=0;
831 0 : clusterValue+=(mp1)->GetSum();
832 0 : list<Int_t> nUsedRows;
833 0 : nUsedRows.push_back((mp1)->GetRow());
834 :
835 0 : maxDig =(mp1)->GetAdc() ;
836 0 : maxSumAdc =(mp1)->GetSum() ;
837 0 : maxTimeBin =(mp1)->GetTime();
838 0 : maxPad =(mp1)->GetPad() ;
839 0 : maxRow =(mp1)->GetRow() ;
840 0 : maxX =(mp1)->GetX();
841 0 : maxY =(mp1)->GetY();
842 0 : maxT =(mp1)->GetT();
843 :
844 : AliSimDigits *digrowTmp;
845 0 : if(fRawData){
846 0 : digrowTmp = (AliSimDigits*)fDigarr->GetRow(iSec,(mp1)->GetRow());
847 0 : }else{
848 0 : digrowTmp = (AliSimDigits*)fDigarr->LoadRow(iSec,(mp1)->GetRow());
849 : }
850 :
851 0 : digrowTmp->ExpandBuffer(); //decrunch
852 :
853 0 : for(Int_t itb=(mp1)->GetBegin(); itb<((mp1)->GetEnd())+1; itb++){
854 0 : Int_t adcTmp = digrowTmp->GetDigitUnchecked(itb,(mp1)->GetPad());
855 0 : AliTPCvtpr *vtpr=new AliTPCvtpr(adcTmp,itb,(mp1)->GetPad(),(mp1)->GetRow(),(mp1)->GetX(),(mp1)->GetY(),(mp1)->GetT());
856 0 : clusterKr.AddDigitToCluster(vtpr);
857 : }
858 0 : clusterKr.SetCenter();//set centr of the cluster
859 :
860 0 : for(Int_t it2 = it1+1; it2 < entriesArr; ++it2 ) {
861 0 : AliPadMax *mp2=(AliPadMax *)maximaInSector->UncheckedAt(it2);
862 0 : if (!mp2) continue;
863 0 : if (TMath::Abs(clusterKr.GetCenterX() - (mp2)->GetX()) > fMaxPadRangeCm) continue;
864 0 : if (TMath::Abs(clusterKr.GetCenterY() - (mp2)->GetY()) > fMaxRowRangeCm) continue;
865 0 : if (TMath::Abs(clusterKr.GetCenterT() - (mp2)->GetT()) > fMaxTimeRange) continue;
866 :
867 : {
868 0 : clusterValue+=(mp2)->GetSum();
869 :
870 0 : nUsedPads++;
871 0 : nUsedRows.push_back((mp2)->GetRow());
872 :
873 : AliSimDigits *digrowTmp1;
874 0 : if(fRawData){
875 0 : digrowTmp1 = (AliSimDigits*)fDigarr->GetRow(iSec,(mp2)->GetRow());
876 0 : }else{
877 0 : digrowTmp1 = (AliSimDigits*)fDigarr->LoadRow(iSec,(mp2)->GetRow());
878 : }
879 0 : digrowTmp1->ExpandBuffer(); //decrunch
880 :
881 0 : for(Int_t itb=(mp2)->GetBegin(); itb<(mp2)->GetEnd()+1; itb++){
882 0 : Int_t adcTmp = digrowTmp1->GetDigitUnchecked(itb,(mp2)->GetPad());
883 0 : AliTPCvtpr *vtpr=new AliTPCvtpr(adcTmp,itb,(mp2)->GetPad(),(mp2)->GetRow(),(mp2)->GetX(),(mp2)->GetY(),(mp2)->GetT());
884 0 : clusterKr.AddDigitToCluster(vtpr);
885 : }
886 :
887 0 : clusterKr.SetCenter();//set center of the cluster
888 :
889 : //which one is bigger
890 0 : if( (mp2)->GetAdc() > maxDig ){
891 0 : maxDig =(mp2)->GetAdc() ;
892 0 : maxSumAdc =(mp2)->GetSum() ;
893 0 : maxTimeBin =(mp2)->GetTime();
894 0 : maxPad =(mp2)->GetPad() ;
895 0 : maxRow =(mp2)->GetRow() ;
896 0 : maxX =(mp2)->GetX() ;
897 0 : maxY =(mp2)->GetY() ;
898 0 : maxT =(mp2)->GetT() ;
899 0 : } else if ( (mp2)->GetAdc() == maxDig ){
900 0 : if( (mp2)->GetSum() > maxSumAdc){
901 0 : maxDig =(mp2)->GetAdc() ;
902 0 : maxSumAdc =(mp2)->GetSum() ;
903 0 : maxTimeBin =(mp2)->GetTime();
904 0 : maxPad =(mp2)->GetPad() ;
905 0 : maxRow =(mp2)->GetRow() ;
906 0 : maxX =(mp2)->GetX() ;
907 0 : maxY =(mp2)->GetY() ;
908 0 : maxT =(mp2)->GetT() ;
909 0 : }
910 : }
911 0 : delete maximaInSector->RemoveAt(it2);
912 : }
913 0 : }//inside loop
914 0 : delete maximaInSector->RemoveAt(it1);
915 0 : clusterKr.SetSize();
916 : //through out clusters on the edge and noise
917 : //if(clusterValue/clusterKr.fCluster.size()<fValueToSize)continue;
918 0 : if(clusterValue/(clusterKr.GetSize())<fValueToSize)continue;
919 :
920 0 : clusterKr.SetADCcluster(clusterValue);
921 0 : clusterKr.SetNPads(nUsedPads);
922 0 : clusterKr.SetMax(AliTPCvtpr(maxDig,maxTimeBin,maxPad,maxRow,maxX,maxY,maxT));
923 0 : clusterKr.SetSec(iSec);
924 0 : clusterKr.SetSize();
925 :
926 0 : nUsedRows.sort();
927 0 : nUsedRows.unique();
928 0 : clusterKr.SetNRows(nUsedRows.size());
929 0 : clusterKr.SetCenter();
930 :
931 0 : clusterKr.SetRMS();//Set pad,row,timebin RMS
932 0 : clusterKr.Set1D();//Set size in pads and timebins
933 :
934 0 : clusterKr.SetTimeStamp(fTimeStamp);
935 0 : clusterKr.SetRun(fRun);
936 :
937 0 : clusterCounter++;
938 :
939 :
940 : //save each cluster into file
941 0 : if (fOutput){
942 0 : (*fOutput)<<"Kr"<<
943 0 : "Cl.="<<&clusterKr<<
944 : "\n";
945 : }
946 : //end of save each cluster into file adc.root
947 0 : }//outer loop
948 0 : }
949 :
950 :
951 :
952 : ////____________________________________________________________________________
953 :
954 :
955 : void AliTPCclustererKr::GetXY(Int_t sec,Int_t row,Int_t pad,Double_t& xGlob,Double_t& yGlob){
956 : //
957 : //gives global XY coordinate of the pad
958 : //
959 :
960 0 : Double_t yLocal = fParam->GetPadRowRadii(sec,row);//radius of row in sector in cm
961 :
962 0 : Int_t padmax = fParam->GetNPads(sec,row);//number of pads in a given row
963 : Float_t padXSize;
964 0 : if(sec<fParam->GetNInnerSector())padXSize=0.4;
965 : else padXSize=0.6;
966 0 : Double_t xLocal=(pad+0.5-padmax/2.)*padXSize;//x-value of the center of pad
967 :
968 0 : Float_t sin,cos;
969 0 : fParam->AdjustCosSin((Int_t)sec,cos,sin);//return sinus and cosinus of the sector
970 :
971 0 : Double_t xGlob1 = xLocal * cos + yLocal * sin;
972 0 : Double_t yGlob1 = -xLocal * sin + yLocal * cos;
973 :
974 :
975 : Double_t rot=0;
976 0 : rot=TMath::Pi()/2.;
977 :
978 0 : xGlob = xGlob1 * TMath::Cos(rot) + yGlob1 * TMath::Sin(rot);
979 0 : yGlob = -xGlob1 * TMath::Sin(rot) + yGlob1 * TMath::Cos(rot);
980 :
981 0 : yGlob=-1*yGlob;
982 0 : if(sec<18||(sec>=36&&sec<54)) xGlob =-1*xGlob;
983 :
984 :
985 : return;
986 0 : }
|
