1 /**************************************************************************
\r
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
\r
4 * Author: The ALICE Off-line Project. *
\r
5 * Contributors are mentioned in the code where appropriate. *
\r
7 * Permission to use, copy, modify and distribute this software and its *
\r
8 * documentation strictly for non-commercial purposes is hereby granted *
\r
9 * without fee, provided that the above copyright notice appears in all *
\r
10 * copies and that both the copyright notice and this permission notice *
\r
11 * appear in the supporting documentation. The authors make no claims *
\r
12 * about the suitability of this software for any purpose. It is *
\r
13 * provided "as is" without express or implied warranty. *
\r
14 **************************************************************************/
\r
16 /* $Id: AliTPCclustererKr.cxx,v 1.7 2008/02/06 17:24:53 matyja Exp $ */
\r
18 //-----------------------------------------------------------------
\r
19 // Implementation of the TPC Kr cluster class
\r
21 // Origin: Adam Matyja, INP PAN, adam.matyja@ifj.edu.pl
\r
22 //-----------------------------------------------------------------
\r
25 Instruction - how to use that:
\r
26 There are two macros prepared. One is for preparing clusters from MC
\r
27 samples: FindKrClusters.C. The output is kept in TPC.RecPoints.root.
\r
28 The other macro is prepared for data analysis: FindKrClustersRaw.C.
\r
29 The output is created for each processed file in root file named adc.root.
\r
30 For each data subsample the same named file is created. So be careful
\r
31 do not overwrite them.
\r
33 Additional selection criteria to select the GOLD cluster
\r
35 // open file with clusters
\r
36 TFile f("Krypton.root");
\r
37 TTree * tree = (TTree*)f.Get("Kr")
\r
38 TCut cutR0("cutR0","fADCcluster/fSize<100"); // adjust it according v seetings -
\r
39 TCut cutR1("cutR1","fADCcluster/fSize>7"); // cosmic tracks and noise removal
\r
40 TCut cutR2("cutR2","fMax.fAdc/fADCcluster<0.2"); // digital noise removal
\r
41 TCut cutR3("cutR3","fMax.fAdc/fADCcluster>0.01"); // noise removal
\r
42 TCut cutS1("cutS1","fSize<200"); // adjust it according v seetings - cosmic tracks
\r
43 TCut cutAll = cutR0+cutR1+cutR2+cutR3+cutS1;
\r
44 This values are typical values to be applied in selectors
\r
51 To run clusterizaton for MC type:
\r
54 If you don't want to use the standard selection criteria then you
\r
55 have to do following:
\r
57 // load the standard setup
\r
58 AliRunLoader* rl = AliRunLoader::Open("galice.root");
\r
59 AliTPCLoader *tpcl = (AliTPCLoader*)rl->GetLoader("TPCLoader");
\r
62 gAlice=rl->GetAliRun();
\r
63 TDirectory *cwd = gDirectory;
\r
64 AliTPCv4 *tpc = (AliTPCv4*)gAlice->GetDetector("TPC");
\r
65 Int_t ver = tpc->IsVersion();
\r
67 AliTPCParam *param=(AliTPCParamSR *)gDirectory->Get("75x40_100x60_150x60");
\r
68 AliTPCDigitsArray *digarr=new AliTPCDigitsArray;
\r
69 digarr->Setup(param);
\r
73 Int_t nevmax=rl->GetNumberOfEvents();
\r
74 for(Int_t nev=0;nev<nevmax ;nev++){
\r
76 TTree* input_tree= tpcl->TreeD();//load tree with digits
\r
77 digarr->ConnectTree(input_tree);
\r
78 TTree *output_tree =tpcl->TreeR();//load output tree
\r
80 AliTPCclustererKr *clusters = new AliTPCclustererKr();
\r
81 clusters->SetParam(param);
\r
82 clusters->SetInput(input_tree);
\r
83 clusters->SetOutput(output_tree);
\r
84 clusters->SetDigArr(digarr);
\r
86 //If you want to change the cluster finder parameters for MC there are
\r
89 //1. signal threshold (everything below the given number is treated as 0)
\r
90 clusters->SetMinAdc(3);
\r
92 //2. number of neighbouring timebins to be considered
\r
93 clusters->SetMinTimeBins(2);
\r
95 //3. distance of the cluster center to the center of a pad in pad-padrow plane
\r
96 //(in cm). Remenber that this is still quantified by pad size.
\r
97 clusters->SetMaxPadRangeCm(2.5);
\r
99 //4. distance of the cluster center to the center of a padrow in pad-padrow
\r
100 //plane (in cm). Remenber that this is still quantified by pad size.
\r
101 clusters->SetMaxRowRangeCm(3.5);
\r
103 //5. distance of the cluster center to the max time bin on a pad (in tackts)
\r
104 //ie. fabs(centerT - time)<7
\r
105 clusters->SetMaxTimeRange(7);
\r
107 //6. cut reduce peak at 0. There are noises which appear mostly as two
\r
108 //timebins on one pad.
\r
109 clusters->SetValueToSize(3.5);
\r
112 clusters->finderIO();
\r
113 tpcl->WriteRecPoints("OVERWRITE");
\r
115 delete rl;//cleans everything
\r
118 ********* DATA *********
\r
121 To run clusterizaton for DATA for file named raw_data.root type:
\r
122 .x FindKrClustersRaw.C("raw_data.root")
\r
124 If you want to change some criteria do the following:
\r
127 // remove Altro warnings
\r
129 AliLog::SetClassDebugLevel("AliTPCRawStream",-5);
\r
130 AliLog::SetClassDebugLevel("AliRawReaderDate",-5);
\r
131 AliLog::SetClassDebugLevel("AliTPCAltroMapping",-5);
\r
132 AliLog::SetModuleDebugLevel("RAW",-5);
\r
135 // Get database with noises
\r
137 // char *ocdbpath = gSystem->Getenv("OCDB_PATH");
\r
138 char *ocdbpath ="local:///afs/cern.ch/alice/tpctest/OCDB";
\r
140 ocdbpath="alien://folder=/alice/data/2007/LHC07w/OCDB/";
\r
142 AliCDBManager * man = AliCDBManager::Instance();
\r
143 man->SetDefaultStorage(ocdbpath);
\r
145 AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
\r
146 AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping();
\r
149 TFile *hfile=new TFile("adc.root","RECREATE","ADC file");
\r
150 // Create a ROOT Tree
\r
151 TTree *mytree = new TTree("Kr","Krypton cluster tree");
\r
153 //define infput file
\r
154 const char *fileName="data.root";
\r
155 AliRawReader *reader = new AliRawReaderRoot(fileName);
\r
156 //AliRawReader *reader = new AliRawReaderDate(fileName);
\r
158 AliAltroRawStreamFast* stream = new AliAltroRawStreamFast(reader);
\r
159 stream->SelectRawData("TPC");
\r
161 //one general output
\r
162 AliTPCclustererKr *clusters = new AliTPCclustererKr();
\r
163 clusters->SetOutput(mytree);
\r
164 clusters->SetRecoParam(0);//standard reco parameters
\r
165 AliTPCParamSR *param=new AliTPCParamSR();
\r
166 clusters->SetParam(param);//TPC parameters(sectors, timebins, etc.)
\r
168 //set cluster finder parameters (from data):
\r
169 //1. zero suppression parameter
\r
170 clusters->SetZeroSup(param->GetZeroSup());
\r
173 clusters->SetFirstBin(60);
\r
176 clusters->SetLastBin(950);
\r
179 clusters->SetMaxNoiseAbs(2);
\r
181 //5. maximal amount of sigma of noise
\r
182 clusters->SetMaxNoiseSigma(3);
\r
184 //The remaining parameters are the same paramters as for MC (see MC section
\r
186 clusters->SetMinAdc(3);
\r
187 clusters->SetMinTimeBins(2);
\r
188 clusters->SetMaxPadRangeCm(2.5);
\r
189 clusters->SetMaxRowRangeCm(3.5);
\r
190 clusters->SetMaxTimeRange(7);
\r
191 clusters->SetValueToSize(3.5);
\r
193 while (reader->NextEvent()) {
\r
194 clusters->FinderIO(reader);
\r
204 #include "AliTPCclustererKr.h"
\r
205 #include "AliTPCclusterKr.h"
\r
206 //#include <vector>
\r
208 #include "TObject.h"
\r
209 #include "AliPadMax.h"
\r
210 #include "AliSimDigits.h"
\r
211 #include "AliTPCv4.h"
\r
212 #include "AliTPCParam.h"
\r
213 #include "AliTPCDigitsArray.h"
\r
214 #include "AliTPCvtpr.h"
\r
215 #include "AliTPCClustersRow.h"
\r
219 #include "TTreeStream.h"
\r
221 #include "AliTPCTransform.h"
\r
223 //used in raw data finder
\r
224 #include "AliTPCROC.h"
\r
225 #include "AliTPCCalPad.h"
\r
226 #include "AliTPCAltroMapping.h"
\r
227 #include "AliTPCcalibDB.h"
\r
228 #include "AliTPCRawStream.h"
\r
229 #include "AliTPCRawStreamV3.h"
\r
230 #include "AliTPCRecoParam.h"
\r
231 #include "AliTPCReconstructor.h"
\r
232 #include "AliRawReader.h"
\r
233 #include "AliTPCCalROC.h"
\r
234 #include "AliRawEventHeaderBase.h"
\r
236 ClassImp(AliTPCclustererKr)
\r
239 AliTPCclustererKr::AliTPCclustererKr()
\r
254 // fMaxPadRange(4),
\r
255 // fMaxRowRange(3),
\r
258 fMaxPadRangeCm(2.5),
\r
259 fMaxRowRangeCm(3.5),
\r
270 // default constructor
\r
274 AliTPCclustererKr::AliTPCclustererKr(const AliTPCclustererKr ¶m)
\r
289 // fMaxPadRange(4),
\r
290 // fMaxRowRange(3),
\r
293 fMaxPadRangeCm(2.5),
\r
294 fMaxRowRangeCm(3.5),
\r
305 // copy constructor
\r
307 fParam = param.fParam;
\r
308 fRecoParam = param.fRecoParam;
\r
309 fRawData = param.fRawData;
\r
310 fInput = param.fInput ;
\r
311 fOutput = param.fOutput;
\r
312 fDigarr = param.fDigarr;
\r
313 fZeroSup = param.fZeroSup ;
\r
314 fFirstBin = param.fFirstBin ;
\r
315 fLastBin = param.fLastBin ;
\r
316 fMaxNoiseAbs = param.fMaxNoiseAbs ;
\r
317 fMaxNoiseSigma = param.fMaxNoiseSigma ;
\r
318 fMinAdc = param.fMinAdc;
\r
319 fMinTimeBins = param.fMinTimeBins;
\r
320 // fMaxPadRange = param.fMaxPadRange ;
\r
321 // fMaxRowRange = param.fMaxRowRange ;
\r
322 fMaxTimeRange = param.fMaxTimeRange;
\r
323 fValueToSize = param.fValueToSize;
\r
324 fMaxPadRangeCm = param.fMaxPadRangeCm;
\r
325 fMaxRowRangeCm = param.fMaxRowRangeCm;
\r
326 fIsolCut = param.fIsolCut;
\r
327 fDebugLevel = param.fDebugLevel;
\r
328 fHistoRow = param.fHistoRow ;
\r
329 fHistoPad = param.fHistoPad ;
\r
330 fHistoTime = param.fHistoTime;
\r
331 fHistoRowPad = param.fHistoRowPad;
\r
332 fTimeStamp = param.fTimeStamp;
\r
337 AliTPCclustererKr & AliTPCclustererKr::operator = (const AliTPCclustererKr & param)
\r
340 // assignment operator
\r
342 fParam = param.fParam;
\r
343 fRecoParam = param.fRecoParam;
\r
344 fRawData = param.fRawData;
\r
345 fInput = param.fInput ;
\r
346 fOutput = param.fOutput;
\r
347 fDigarr = param.fDigarr;
\r
348 fZeroSup = param.fZeroSup ;
\r
349 fFirstBin = param.fFirstBin ;
\r
350 fLastBin = param.fLastBin ;
\r
351 fMaxNoiseAbs = param.fMaxNoiseAbs ;
\r
352 fMaxNoiseSigma = param.fMaxNoiseSigma ;
\r
353 fMinAdc = param.fMinAdc;
\r
354 fMinTimeBins = param.fMinTimeBins;
\r
355 // fMaxPadRange = param.fMaxPadRange ;
\r
356 // fMaxRowRange = param.fMaxRowRange ;
\r
357 fMaxTimeRange = param.fMaxTimeRange;
\r
358 fValueToSize = param.fValueToSize;
\r
359 fMaxPadRangeCm = param.fMaxPadRangeCm;
\r
360 fMaxRowRangeCm = param.fMaxRowRangeCm;
\r
361 fIsolCut = param.fIsolCut;
\r
362 fDebugLevel = param.fDebugLevel;
\r
363 fHistoRow = param.fHistoRow ;
\r
364 fHistoPad = param.fHistoPad ;
\r
365 fHistoTime = param.fHistoTime;
\r
366 fHistoRowPad = param.fHistoRowPad;
\r
367 fTimeStamp = param.fTimeStamp;
\r
372 AliTPCclustererKr::~AliTPCclustererKr()
\r
380 void AliTPCclustererKr::SetRecoParam(AliTPCRecoParam *recoParam)
\r
383 // set reconstruction parameters
\r
386 fRecoParam = recoParam;
\r
388 //set default parameters if not specified
\r
389 fRecoParam = AliTPCReconstructor::GetRecoParam();
\r
390 if (!fRecoParam) fRecoParam = AliTPCRecoParam::GetLowFluxParam();
\r
396 ////____________________________________________________________________________
\r
398 void AliTPCclustererKr::SetInput(TTree * tree)
\r
401 // set input tree with digits
\r
404 if (!fInput->GetBranch("Segment")){
\r
405 cerr<<"AliTPCclusterKr::FindClusterKr(): no proper input tree !\n";
\r
411 void AliTPCclustererKr::SetOutput(TTree * /*tree*/)
\r
416 fOutput = new TTreeSRedirector("Krypton.root");
\r
419 ////____________________________________________________________________________
\r
421 Int_t AliTPCclustererKr::FinderIO()
\r
423 // Krypton cluster finder for simulated events from MC
\r
426 Error("Digits2Clusters", "input tree not initialised");
\r
431 Error("Digits2Clusters", "output tree not initialised");
\r
441 Int_t AliTPCclustererKr::FinderIO(AliRawReader* rawReader)
\r
443 // Krypton cluster finder for the TPC raw data
\r
444 // this method is unsing AliAltroRawStreamV3
\r
445 // fParam must be defined before
\r
446 if (!rawReader) return 1;
\r
448 fRawData=kTRUE; //set flag to data
\r
451 Error("Digits2Clusters", "output tree not initialised");
\r
455 fParam->SetMaxTBin(fRecoParam->GetLastBin());//set number of timebins from reco -> param
\r
456 // used later for memory allocation
\r
458 AliRawEventHeaderBase* eventHeader = (AliRawEventHeaderBase*)rawReader->GetEventHeader();
\r
460 fTimeStamp = eventHeader->Get("Timestamp");
\r
461 fRun = rawReader->GetRunNumber();
\r
465 Bool_t isAltro=kFALSE;
\r
467 AliTPCROC * roc = AliTPCROC::Instance();
\r
468 AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
\r
469 AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping();
\r
471 AliTPCRawStreamV3 input(rawReader,(AliAltroMapping**)mapping);
\r
473 const Int_t kNIS = fParam->GetNInnerSector();//number of inner sectors
\r
474 const Int_t kNOS = fParam->GetNOuterSector();//number of outer sectors
\r
475 const Int_t kNS = kNIS + kNOS;//all sectors
\r
479 AliTPCDigitsArray *digarr=new AliTPCDigitsArray(kFALSE);//data not sim
\r
480 digarr->Setup(fParam);//as usually parameters
\r
482 for(Int_t iSec = 0; iSec < kNS; iSec++) {
\r
483 AliTPCCalROC * noiseROC;
\r
484 AliTPCCalROC noiseDummy(iSec);
\r
486 noiseROC = &noiseDummy;//noise=0
\r
488 noiseROC = noiseTPC->GetCalROC(iSec); // noise per given sector
\r
490 Int_t nRows = 0; //number of rows in sector
\r
491 Int_t nDDLs = 0; //number of DDLs
\r
492 Int_t indexDDL = 0; //DDL index
\r
494 nRows = fParam->GetNRowLow();
\r
496 indexDDL = iSec * 2;
\r
498 nRows = fParam->GetNRowUp();
\r
500 indexDDL = (iSec-kNIS) * 4 + kNIS * 2;
\r
504 // Load the raw data for corresponding DDLs
\r
506 rawReader->Reset();
\r
507 rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1);
\r
510 while (input.NextDDL()){
\r
511 // Allocate memory for rows in sector (pads(depends on row) x timebins)
\r
512 if (!digarr->GetRow(iSec,0)){
\r
513 for(Int_t iRow = 0; iRow < nRows; iRow++) {
\r
514 digarr->CreateRow(iSec,iRow);
\r
515 }//end loop over rows
\r
518 while ( input.NextChannel() ) {
\r
519 Int_t iRow = input.GetRow();
\r
520 Int_t iPad = input.GetPad();
\r
521 //check row consistency
\r
522 if (iRow < 0 ) continue;
\r
523 if (iRow < 0 || iRow >= nRows){
\r
524 AliError(Form("Pad-row index (%d) outside the range (%d -> %d) !",
\r
525 iRow, 0, nRows -1));
\r
529 //check pad consistency
\r
530 if (iPad < 0 || iPad >= (Int_t)(roc->GetNPads(iSec,iRow))) {
\r
531 AliError(Form("Pad index (%d) outside the range (%d -> %d) !",
\r
532 iPad, 0, roc->GetNPads(iSec,iRow) ));
\r
536 //loop over bunches
\r
537 while ( input.NextBunch() ){
\r
538 Int_t startTbin = (Int_t)input.GetStartTimeBin();
\r
539 Int_t bunchlength = (Int_t)input.GetBunchLength();
\r
540 const UShort_t *sig = input.GetSignals();
\r
542 for (Int_t iTime = 0; iTime<bunchlength; iTime++){
\r
543 Int_t iTimeBin=startTbin-iTime;
\r
545 if(fDebugLevel==72){
\r
546 fHistoRow->Fill(iRow);
\r
547 fHistoPad->Fill(iPad);
\r
548 fHistoTime->Fill(iTimeBin);
\r
549 fHistoRowPad->Fill(iPad,iRow);
\r
550 }else if(fDebugLevel>=0&&fDebugLevel<72){
\r
551 if(iSec==fDebugLevel){
\r
552 fHistoRow->Fill(iRow);
\r
553 fHistoPad->Fill(iPad);
\r
554 fHistoTime->Fill(iTimeBin);
\r
555 fHistoRowPad->Fill(iPad,iRow);
\r
557 }else if(fDebugLevel==73){
\r
559 fHistoRow->Fill(iRow);
\r
560 fHistoPad->Fill(iPad);
\r
561 fHistoTime->Fill(iTimeBin);
\r
562 fHistoRowPad->Fill(iPad,iRow);
\r
564 }else if(fDebugLevel==74){
\r
566 fHistoRow->Fill(iRow);
\r
567 fHistoPad->Fill(iPad);
\r
568 fHistoTime->Fill(iTimeBin);
\r
569 fHistoRowPad->Fill(iPad,iRow);
\r
573 //check time consistency
\r
574 if ( iTimeBin < fRecoParam->GetFirstBin() || iTimeBin >= fRecoParam->GetLastBin()){
\r
575 //cout<<iTimeBin<<endl;
\r
577 AliFatal(Form("Timebin index (%d) outside the range (%d -> %d) !",
\r
578 iTimeBin, 0, fRecoParam->GetLastBin() -1));
\r
581 Float_t signal=(Float_t)sig[iTime];
\r
582 if (signal <= fZeroSup ||
\r
583 iTimeBin < fFirstBin ||
\r
584 iTimeBin > fLastBin
\r
586 digarr->GetRow(iSec,iRow)->SetDigitFast(0,iTimeBin,iPad);
\r
589 if (!noiseROC) continue;
\r
590 Double_t noiseOnPad = noiseROC->GetValue(iRow,iPad);//noise on given pad and row in sector
\r
591 if (noiseOnPad > fMaxNoiseAbs){
\r
592 digarr->GetRow(iSec,iRow)->SetDigitFast(0,iTimeBin,iPad);
\r
593 continue; // consider noisy pad as dead
\r
595 if(signal <= fMaxNoiseSigma * noiseOnPad){
\r
596 digarr->GetRow(iSec,iRow)->SetDigitFast(0,iTimeBin,iPad);
\r
599 digarr->GetRow(iSec,iRow)->SetDigitFast(TMath::Nint(signal),iTimeBin,iPad);
\r
600 }// end loop signals in bunch
\r
601 }// end loop bunches
\r
604 }// end sector loop
\r
606 if(isAltro) FindClusterKrIO();
\r
616 Int_t AliTPCclustererKr::FinderIOold(AliRawReader* rawReader)
\r
618 // Krypton cluster finder for the TPC raw data
\r
620 // fParam must be defined before
\r
621 if (!rawReader) return 1;
\r
623 if(rawReader)fRawData=kTRUE; //set flag to data
\r
626 Error("Digits2Clusters", "output tree not initialised");
\r
630 fParam->SetMaxTBin(fRecoParam->GetLastBin());//set number of timebins from reco -> param
\r
631 // used later for memory allocation
\r
633 AliRawEventHeaderBase* eventHeader = (AliRawEventHeaderBase*)rawReader->GetEventHeader();
\r
635 fTimeStamp = eventHeader->Get("Timestamp");
\r
636 fRun = rawReader->GetRunNumber();
\r
639 Bool_t isAltro=kFALSE;
\r
641 AliTPCROC * roc = AliTPCROC::Instance();
\r
642 AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
\r
643 AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping();
\r
645 AliTPCRawStream input(rawReader,(AliAltroMapping**)mapping);
\r
647 const Int_t kNIS = fParam->GetNInnerSector();//number of inner sectors
\r
648 const Int_t kNOS = fParam->GetNOuterSector();//number of outer sectors
\r
649 const Int_t kNS = kNIS + kNOS;//all sectors
\r
653 AliTPCDigitsArray *digarr=new AliTPCDigitsArray(kFALSE);//data not sim
\r
654 digarr->Setup(fParam);//as usually parameters
\r
657 // Loop over sectors
\r
659 for(Int_t iSec = 0; iSec < kNS; iSec++) {
\r
660 AliTPCCalROC * noiseROC;
\r
662 noiseROC =new AliTPCCalROC(iSec);//noise=0
\r
665 noiseROC = noiseTPC->GetCalROC(iSec); // noise per given sector
\r
667 Int_t nRows = 0; //number of rows in sector
\r
668 Int_t nDDLs = 0; //number of DDLs
\r
669 Int_t indexDDL = 0; //DDL index
\r
671 nRows = fParam->GetNRowLow();
\r
673 indexDDL = iSec * 2;
\r
675 nRows = fParam->GetNRowUp();
\r
677 indexDDL = (iSec-kNIS) * 4 + kNIS * 2;
\r
681 // Load the raw data for corresponding DDLs
\r
683 rawReader->Reset();
\r
684 rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1);
\r
688 // Allocate memory for rows in sector (pads(depends on row) x timebins)
\r
689 for(Int_t iRow = 0; iRow < nRows; iRow++) {
\r
690 digarr->CreateRow(iSec,iRow);
\r
691 }//end loop over rows
\r
693 rawReader->Reset();
\r
694 rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1);
\r
697 // Begin loop over altro data
\r
699 while (input.Next()) {
\r
701 //check sector consistency
\r
702 if (input.GetSector() != iSec)
\r
703 AliFatal(Form("Sector index mismatch ! Expected (%d), but got (%d) !",iSec,input.GetSector()));
\r
705 Int_t iRow = input.GetRow();
\r
706 Int_t iPad = input.GetPad();
\r
707 Int_t iTimeBin = input.GetTime();
\r
710 if(fDebugLevel==72){
\r
711 fHistoRow->Fill(iRow);
\r
712 fHistoPad->Fill(iPad);
\r
713 fHistoTime->Fill(iTimeBin);
\r
714 fHistoRowPad->Fill(iPad,iRow);
\r
715 }else if(fDebugLevel>=0&&fDebugLevel<72){
\r
716 if(iSec==fDebugLevel){
\r
717 fHistoRow->Fill(iRow);
\r
718 fHistoPad->Fill(iPad);
\r
719 fHistoTime->Fill(iTimeBin);
\r
720 fHistoRowPad->Fill(iPad,iRow);
\r
722 }else if(fDebugLevel==73){
\r
724 fHistoRow->Fill(iRow);
\r
725 fHistoPad->Fill(iPad);
\r
726 fHistoTime->Fill(iTimeBin);
\r
727 fHistoRowPad->Fill(iPad,iRow);
\r
729 }else if(fDebugLevel==74){
\r
731 fHistoRow->Fill(iRow);
\r
732 fHistoPad->Fill(iPad);
\r
733 fHistoTime->Fill(iTimeBin);
\r
734 fHistoRowPad->Fill(iPad,iRow);
\r
738 //check row consistency
\r
739 if (iRow < 0 ) continue;
\r
740 if (iRow < 0 || iRow >= nRows){
\r
741 AliError(Form("Pad-row index (%d) outside the range (%d -> %d) !",
\r
742 iRow, 0, nRows -1));
\r
746 //check pad consistency
\r
747 if (iPad < 0 || iPad >= (Int_t)(roc->GetNPads(iSec,iRow))) {
\r
748 AliError(Form("Pad index (%d) outside the range (%d -> %d) !",
\r
749 iPad, 0, roc->GetNPads(iSec,iRow) ));
\r
753 //check time consistency
\r
754 if ( iTimeBin < fRecoParam->GetFirstBin() || iTimeBin >= fRecoParam->GetLastBin()){
\r
755 //cout<<iTimeBin<<endl;
\r
757 AliFatal(Form("Timebin index (%d) outside the range (%d -> %d) !",
\r
758 iTimeBin, 0, fRecoParam->GetLastBin() -1));
\r
762 Int_t signal = input.GetSignal();
\r
763 if (signal <= fZeroSup ||
\r
764 iTimeBin < fFirstBin ||
\r
765 iTimeBin > fLastBin
\r
767 digarr->GetRow(iSec,iRow)->SetDigitFast(0,iTimeBin,iPad);
\r
770 if (!noiseROC) continue;
\r
771 Double_t noiseOnPad = noiseROC->GetValue(iRow,iPad);//noise on given pad and row in sector
\r
772 if (noiseOnPad > fMaxNoiseAbs){
\r
773 digarr->GetRow(iSec,iRow)->SetDigitFast(0,iTimeBin,iPad);
\r
774 continue; // consider noisy pad as dead
\r
776 if(signal <= fMaxNoiseSigma * noiseOnPad){
\r
777 digarr->GetRow(iSec,iRow)->SetDigitFast(0,iTimeBin,iPad);
\r
780 digarr->GetRow(iSec,iRow)->SetDigitFast(signal,iTimeBin,iPad);
\r
781 }//end of loop over altro data
\r
782 }//end of loop over sectors
\r
785 if(isAltro) FindClusterKrIO();
\r
791 void AliTPCclustererKr::CleanSector(Int_t sector){
\r
793 // clean isolated digits
\r
795 const Int_t kNRows=fParam->GetNRow(sector);//number of rows in sector
\r
796 for(Int_t iRow=0; iRow<kNRows; ++iRow){
\r
797 AliSimDigits *digrow;
\r
799 digrow = (AliSimDigits*)fDigarr->GetRow(sector,iRow);//real data
\r
801 digrow = (AliSimDigits*)fDigarr->LoadRow(sector,iRow);//MC
\r
803 if(!digrow) continue;
\r
804 digrow->ExpandBuffer(); //decrunch
\r
805 const Int_t kNPads = digrow->GetNCols(); // number of pads
\r
806 const Int_t kNTime = digrow->GetNRows(); // number of timebins
\r
807 for(Int_t iPad=1;iPad<kNPads-1;iPad++){
\r
808 Short_t* val = digrow->GetDigitsColumn(iPad);
\r
810 for(Int_t iTimeBin=1;iTimeBin<kNTime-1;iTimeBin++){
\r
811 if (val[iTimeBin]<=0) continue;
\r
812 if (val[iTimeBin-1]+val[iTimeBin+1]<fIsolCut) {val[iTimeBin]=0; continue;}
\r
813 if (val[iTimeBin-kNTime]+val[iTimeBin+kNTime]<fIsolCut) {val[iTimeBin]=0; continue;}
\r
815 if (val[iTimeBin-1-kNTime]+val[iTimeBin+1+kNTime]<fIsolCut) {val[iTimeBin]=0; continue;}
\r
816 if (val[iTimeBin+1-kNTime]+val[iTimeBin-1+kNTime]<fIsolCut) {val[iTimeBin]=0; continue;}
\r
824 ////____________________________________________________________________________
\r
825 Int_t AliTPCclustererKr::FindClusterKrIO()
\r
829 //fParam and fDigarr must be set to run this method
\r
832 Int_t clusterCounter=0;
\r
833 const Int_t nTotalSector=fParam->GetNSector();//number of sectors
\r
834 for(Int_t iSec=0; iSec<nTotalSector; ++iSec){
\r
837 //vector of maxima for each sector
\r
838 //std::vector<AliPadMax*> maximaInSector;
\r
839 TObjArray *maximaInSector=new TObjArray();//to store AliPadMax*
\r
842 // looking for the maxima on the pad
\r
845 const Int_t kNRows=fParam->GetNRow(iSec);//number of rows in sector
\r
846 for(Int_t iRow=0; iRow<kNRows; ++iRow){
\r
847 AliSimDigits *digrow;
\r
849 digrow = (AliSimDigits*)fDigarr->GetRow(iSec,iRow);//real data
\r
851 digrow = (AliSimDigits*)fDigarr->LoadRow(iSec,iRow);//MC
\r
853 if(digrow){//if pointer exist
\r
854 digrow->ExpandBuffer(); //decrunch
\r
855 const Int_t kNPads = digrow->GetNCols(); // number of pads
\r
856 const Int_t kNTime = digrow->GetNRows(); // number of timebins
\r
857 for(Int_t iPad=0;iPad<kNPads;iPad++){
\r
859 Int_t timeBinMax=-1;//timebin of maximum
\r
860 Int_t valueMaximum=-1;//value of maximum in adc
\r
861 Int_t increaseBegin=-1;//timebin when increase starts
\r
862 Int_t sumAdc=0;//sum of adc on the pad in maximum surrounding
\r
863 bool ifIncreaseBegin=true;//flag - check if increasing started
\r
864 bool ifMaximum=false;//flag - check if it could be maximum
\r
865 Short_t* val = digrow->GetDigitsColumn(iPad);
\r
866 for(Int_t iTimeBin=1;iTimeBin<kNTime-1;iTimeBin++){
\r
868 if (val[iTimeBin]==-1) break; // 0 until the end
\r
869 for( ; iTimeBin<kNTime-2&&val[iTimeBin]<fMinAdc ;iTimeBin++) {}
\r
872 Short_t adc = val[iTimeBin];
\r
874 if(adc<fMinAdc){//standard was 3 for fMinAdc
\r
876 if(iTimeBin-increaseBegin<fMinTimeBins){//at least 2 time bins
\r
881 ifIncreaseBegin=true;
\r
885 //insert maximum, default values and set flags
\r
886 //Double_t xCord,yCord;
\r
887 //GetXY(iSec,iRow,iPad,xCord,yCord);
\r
888 Double_t x[]={iRow,iPad,iTimeBin};
\r
890 AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
\r
892 transform->Transform(x,i,0,1);
\r
894 AliPadMax *oneMaximum = new AliPadMax(AliTPCvtpr(valueMaximum,
\r
900 x[2]/*timeBinMax*/),
\r
904 maximaInSector->AddLast(oneMaximum);
\r
910 ifIncreaseBegin=true;
\r
921 if(ifIncreaseBegin){
\r
922 ifIncreaseBegin=false;
\r
923 increaseBegin=iTimeBin;
\r
926 if(adc>valueMaximum){
\r
927 timeBinMax=iTimeBin;
\r
932 if(iTimeBin==kNTime-1 && ifMaximum && kNTime-increaseBegin>fMinTimeBins){//on the edge
\r
933 //at least 3 timebins
\r
934 //insert maximum, default values and set flags
\r
935 //Double_t xCord,yCord;
\r
936 //GetXY(iSec,iRow,iPad,xCord,yCord);
\r
937 Double_t x[]={iRow,iPad,iTimeBin};
\r
939 //AliTPCTransform trafo;
\r
940 //trafo.Transform(x,i,0,1);
\r
942 AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
\r
944 transform->Transform(x,i,0,1);
\r
946 AliPadMax *oneMaximum = new AliPadMax(AliTPCvtpr(valueMaximum,
\r
952 x[2]/*timeBinMax*/),
\r
956 maximaInSector->AddLast(oneMaximum);
\r
962 ifIncreaseBegin=true;
\r
967 }//end loop over timebins
\r
968 }//end loop over pads
\r
970 // cout<<"Pointer does not exist!!"<<endl;
\r
971 }//end if poiner exists
\r
972 }//end loop over rows
\r
974 MakeClusters(maximaInSector,iSec,clusterCounter);
\r
976 maximaInSector->SetOwner(kTRUE);
\r
977 maximaInSector->Delete();
\r
978 delete maximaInSector;
\r
980 cout<<"Number of clusters in event: "<<clusterCounter<<endl;
\r
984 void AliTPCclustererKr::MakeClusters(TObjArray * maximaInSector, Int_t iSec, Int_t &clusterCounter){
\r
991 Int_t maxTimeBin=0;
\r
997 Int_t entriesArr = maximaInSector->GetEntriesFast();
\r
998 for(Int_t it1 = 0; it1 < entriesArr; ++it1 ) {
\r
1000 AliPadMax *mp1=(AliPadMax *)maximaInSector->UncheckedAt(it1);
\r
1001 if (!mp1) continue;
\r
1002 AliTPCclusterKr clusterKr;
\r
1004 Int_t nUsedPads=1;
\r
1005 Int_t clusterValue=0;
\r
1006 clusterValue+=(mp1)->GetSum();
\r
1007 list<Int_t> nUsedRows;
\r
1008 nUsedRows.push_back((mp1)->GetRow());
\r
1010 maxDig =(mp1)->GetAdc() ;
\r
1011 maxSumAdc =(mp1)->GetSum() ;
\r
1012 maxTimeBin =(mp1)->GetTime();
\r
1013 maxPad =(mp1)->GetPad() ;
\r
1014 maxRow =(mp1)->GetRow() ;
\r
1015 maxX =(mp1)->GetX();
\r
1016 maxY =(mp1)->GetY();
\r
1017 maxT =(mp1)->GetT();
\r
1019 AliSimDigits *digrowTmp;
\r
1021 digrowTmp = (AliSimDigits*)fDigarr->GetRow(iSec,(mp1)->GetRow());
\r
1023 digrowTmp = (AliSimDigits*)fDigarr->LoadRow(iSec,(mp1)->GetRow());
\r
1026 digrowTmp->ExpandBuffer(); //decrunch
\r
1028 for(Int_t itb=(mp1)->GetBegin(); itb<((mp1)->GetEnd())+1; itb++){
\r
1029 Int_t adcTmp = digrowTmp->GetDigitUnchecked(itb,(mp1)->GetPad());
\r
1030 AliTPCvtpr *vtpr=new AliTPCvtpr(adcTmp,itb,(mp1)->GetPad(),(mp1)->GetRow(),(mp1)->GetX(),(mp1)->GetY(),(mp1)->GetT());
\r
1031 clusterKr.AddDigitToCluster(vtpr);
\r
1033 clusterKr.SetCenter();//set centr of the cluster
\r
1035 for(Int_t it2 = it1+1; it2 < entriesArr; ++it2 ) {
\r
1036 AliPadMax *mp2=(AliPadMax *)maximaInSector->UncheckedAt(it2);
\r
1037 if (!mp2) continue;
\r
1038 if (TMath::Abs(clusterKr.GetCenterX() - (mp2)->GetX()) > fMaxPadRangeCm) continue;
\r
1039 if (TMath::Abs(clusterKr.GetCenterY() - (mp2)->GetY()) > fMaxRowRangeCm) continue;
\r
1040 if (TMath::Abs(clusterKr.GetCenterT() - (mp2)->GetT()) > fMaxTimeRange) continue;
\r
1043 clusterValue+=(mp2)->GetSum();
\r
1046 nUsedRows.push_back((mp2)->GetRow());
\r
1048 AliSimDigits *digrowTmp1;
\r
1050 digrowTmp1 = (AliSimDigits*)fDigarr->GetRow(iSec,(mp2)->GetRow());
\r
1052 digrowTmp1 = (AliSimDigits*)fDigarr->LoadRow(iSec,(mp2)->GetRow());
\r
1054 digrowTmp1->ExpandBuffer(); //decrunch
\r
1056 for(Int_t itb=(mp2)->GetBegin(); itb<(mp2)->GetEnd()+1; itb++){
\r
1057 Int_t adcTmp = digrowTmp1->GetDigitUnchecked(itb,(mp2)->GetPad());
\r
1058 AliTPCvtpr *vtpr=new AliTPCvtpr(adcTmp,itb,(mp2)->GetPad(),(mp2)->GetRow(),(mp2)->GetX(),(mp2)->GetY(),(mp2)->GetT());
\r
1059 clusterKr.AddDigitToCluster(vtpr);
\r
1062 clusterKr.SetCenter();//set center of the cluster
\r
1064 //which one is bigger
\r
1065 if( (mp2)->GetAdc() > maxDig ){
\r
1066 maxDig =(mp2)->GetAdc() ;
\r
1067 maxSumAdc =(mp2)->GetSum() ;
\r
1068 maxTimeBin =(mp2)->GetTime();
\r
1069 maxPad =(mp2)->GetPad() ;
\r
1070 maxRow =(mp2)->GetRow() ;
\r
1071 maxX =(mp2)->GetX() ;
\r
1072 maxY =(mp2)->GetY() ;
\r
1073 maxT =(mp2)->GetT() ;
\r
1074 } else if ( (mp2)->GetAdc() == maxDig ){
\r
1075 if( (mp2)->GetSum() > maxSumAdc){
\r
1076 maxDig =(mp2)->GetAdc() ;
\r
1077 maxSumAdc =(mp2)->GetSum() ;
\r
1078 maxTimeBin =(mp2)->GetTime();
\r
1079 maxPad =(mp2)->GetPad() ;
\r
1080 maxRow =(mp2)->GetRow() ;
\r
1081 maxX =(mp2)->GetX() ;
\r
1082 maxY =(mp2)->GetY() ;
\r
1083 maxT =(mp2)->GetT() ;
\r
1086 delete maximaInSector->RemoveAt(it2);
\r
1089 delete maximaInSector->RemoveAt(it1);
\r
1090 clusterKr.SetSize();
\r
1091 //through out clusters on the edge and noise
\r
1092 //if(clusterValue/clusterKr.fCluster.size()<fValueToSize)continue;
\r
1093 if(clusterValue/(clusterKr.GetSize())<fValueToSize)continue;
\r
1095 clusterKr.SetADCcluster(clusterValue);
\r
1096 clusterKr.SetNPads(nUsedPads);
\r
1097 clusterKr.SetMax(AliTPCvtpr(maxDig,maxTimeBin,maxPad,maxRow,maxX,maxY,maxT));
\r
1098 clusterKr.SetSec(iSec);
\r
1099 clusterKr.SetSize();
\r
1102 nUsedRows.unique();
\r
1103 clusterKr.SetNRows(nUsedRows.size());
\r
1104 clusterKr.SetCenter();
\r
1106 clusterKr.SetRMS();//Set pad,row,timebin RMS
\r
1107 clusterKr.Set1D();//Set size in pads and timebins
\r
1109 clusterKr.SetTimeStamp(fTimeStamp);
\r
1110 clusterKr.SetRun(fRun);
\r
1115 //save each cluster into file
\r
1117 (*fOutput)<<"Kr"<<
\r
1118 "Cl.="<<&clusterKr<<
\r
1121 //end of save each cluster into file adc.root
\r
1127 ////____________________________________________________________________________
\r
1130 void AliTPCclustererKr::GetXY(Int_t sec,Int_t row,Int_t pad,Double_t& xGlob,Double_t& yGlob){
\r
1132 //gives global XY coordinate of the pad
\r
1135 Double_t yLocal = fParam->GetPadRowRadii(sec,row);//radius of row in sector in cm
\r
1137 Int_t padmax = fParam->GetNPads(sec,row);//number of pads in a given row
\r
1139 if(sec<fParam->GetNInnerSector())padXSize=0.4;
\r
1140 else padXSize=0.6;
\r
1141 Double_t xLocal=(pad+0.5-padmax/2.)*padXSize;//x-value of the center of pad
\r
1144 fParam->AdjustCosSin((Int_t)sec,cos,sin);//return sinus and cosinus of the sector
\r
1146 Double_t xGlob1 = xLocal * cos + yLocal * sin;
\r
1147 Double_t yGlob1 = -xLocal * sin + yLocal * cos;
\r
1151 rot=TMath::Pi()/2.;
\r
1153 xGlob = xGlob1 * TMath::Cos(rot) + yGlob1 * TMath::Sin(rot);
\r
1154 yGlob = -xGlob1 * TMath::Sin(rot) + yGlob1 * TMath::Cos(rot);
\r
1157 if(sec<18||(sec>=36&&sec<54)) xGlob =-1*xGlob;
\r