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
269 // default constructor
\r
273 AliTPCclustererKr::AliTPCclustererKr(const AliTPCclustererKr ¶m)
\r
288 // fMaxPadRange(4),
\r
289 // fMaxRowRange(3),
\r
292 fMaxPadRangeCm(2.5),
\r
293 fMaxRowRangeCm(3.5),
\r
303 // copy constructor
\r
305 fParam = param.fParam;
\r
306 fRecoParam = param.fRecoParam;
\r
307 fRawData = param.fRawData;
\r
308 fInput = param.fInput ;
\r
309 fOutput = param.fOutput;
\r
310 fDigarr = param.fDigarr;
\r
311 fZeroSup = param.fZeroSup ;
\r
312 fFirstBin = param.fFirstBin ;
\r
313 fLastBin = param.fLastBin ;
\r
314 fMaxNoiseAbs = param.fMaxNoiseAbs ;
\r
315 fMaxNoiseSigma = param.fMaxNoiseSigma ;
\r
316 fMinAdc = param.fMinAdc;
\r
317 fMinTimeBins = param.fMinTimeBins;
\r
318 // fMaxPadRange = param.fMaxPadRange ;
\r
319 // fMaxRowRange = param.fMaxRowRange ;
\r
320 fMaxTimeRange = param.fMaxTimeRange;
\r
321 fValueToSize = param.fValueToSize;
\r
322 fMaxPadRangeCm = param.fMaxPadRangeCm;
\r
323 fMaxRowRangeCm = param.fMaxRowRangeCm;
\r
324 fIsolCut = param.fIsolCut;
\r
325 fDebugLevel = param.fDebugLevel;
\r
326 fHistoRow = param.fHistoRow ;
\r
327 fHistoPad = param.fHistoPad ;
\r
328 fHistoTime = param.fHistoTime;
\r
329 fHistoRowPad = param.fHistoRowPad;
\r
330 fTimeStamp = param.fTimeStamp;
\r
334 AliTPCclustererKr & AliTPCclustererKr::operator = (const AliTPCclustererKr & param)
\r
337 // assignment operator
\r
339 fParam = param.fParam;
\r
340 fRecoParam = param.fRecoParam;
\r
341 fRawData = param.fRawData;
\r
342 fInput = param.fInput ;
\r
343 fOutput = param.fOutput;
\r
344 fDigarr = param.fDigarr;
\r
345 fZeroSup = param.fZeroSup ;
\r
346 fFirstBin = param.fFirstBin ;
\r
347 fLastBin = param.fLastBin ;
\r
348 fMaxNoiseAbs = param.fMaxNoiseAbs ;
\r
349 fMaxNoiseSigma = param.fMaxNoiseSigma ;
\r
350 fMinAdc = param.fMinAdc;
\r
351 fMinTimeBins = param.fMinTimeBins;
\r
352 // fMaxPadRange = param.fMaxPadRange ;
\r
353 // fMaxRowRange = param.fMaxRowRange ;
\r
354 fMaxTimeRange = param.fMaxTimeRange;
\r
355 fValueToSize = param.fValueToSize;
\r
356 fMaxPadRangeCm = param.fMaxPadRangeCm;
\r
357 fMaxRowRangeCm = param.fMaxRowRangeCm;
\r
358 fIsolCut = param.fIsolCut;
\r
359 fDebugLevel = param.fDebugLevel;
\r
360 fHistoRow = param.fHistoRow ;
\r
361 fHistoPad = param.fHistoPad ;
\r
362 fHistoTime = param.fHistoTime;
\r
363 fHistoRowPad = param.fHistoRowPad;
\r
364 fTimeStamp = param.fTimeStamp;
\r
368 AliTPCclustererKr::~AliTPCclustererKr()
\r
376 void AliTPCclustererKr::SetRecoParam(AliTPCRecoParam *recoParam)
\r
379 // set reconstruction parameters
\r
382 fRecoParam = recoParam;
\r
384 //set default parameters if not specified
\r
385 fRecoParam = AliTPCReconstructor::GetRecoParam();
\r
386 if (!fRecoParam) fRecoParam = AliTPCRecoParam::GetLowFluxParam();
\r
392 ////____________________________________________________________________________
\r
394 void AliTPCclustererKr::SetInput(TTree * tree)
\r
397 // set input tree with digits
\r
400 if (!fInput->GetBranch("Segment")){
\r
401 cerr<<"AliTPCclusterKr::FindClusterKr(): no proper input tree !\n";
\r
407 void AliTPCclustererKr::SetOutput(TTree * /*tree*/)
\r
412 fOutput = new TTreeSRedirector("Krypton.root");
\r
415 ////____________________________________________________________________________
\r
417 Int_t AliTPCclustererKr::FinderIO()
\r
419 // Krypton cluster finder for simulated events from MC
\r
422 Error("Digits2Clusters", "input tree not initialised");
\r
427 Error("Digits2Clusters", "output tree not initialised");
\r
437 Int_t AliTPCclustererKr::FinderIO(AliRawReader* rawReader)
\r
439 // Krypton cluster finder for the TPC raw data
\r
440 // this method is unsing AliAltroRawStreamV3
\r
441 // fParam must be defined before
\r
443 if(rawReader)fRawData=kTRUE; //set flag to data
\r
446 Error("Digits2Clusters", "output tree not initialised");
\r
450 fParam->SetMaxTBin(fRecoParam->GetLastBin());//set number of timebins from reco -> param
\r
451 // used later for memory allocation
\r
453 AliRawEventHeaderBase* eventHeader = (AliRawEventHeaderBase*)rawReader->GetEventHeader();
\r
455 fTimeStamp = eventHeader->Get("Timestamp");
\r
459 Bool_t isAltro=kFALSE;
\r
461 AliTPCROC * roc = AliTPCROC::Instance();
\r
462 AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
\r
463 AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping();
\r
465 AliTPCRawStreamV3 input(rawReader,(AliAltroMapping**)mapping);
\r
467 const Int_t kNIS = fParam->GetNInnerSector();//number of inner sectors
\r
468 const Int_t kNOS = fParam->GetNOuterSector();//number of outer sectors
\r
469 const Int_t kNS = kNIS + kNOS;//all sectors
\r
473 AliTPCDigitsArray *digarr=new AliTPCDigitsArray(kFALSE);//data not sim
\r
474 digarr->Setup(fParam);//as usually parameters
\r
476 for(Int_t iSec = 0; iSec < kNS; iSec++) {
\r
477 AliTPCCalROC * noiseROC;
\r
478 AliTPCCalROC noiseDummy(iSec);
\r
480 noiseROC = &noiseDummy;//noise=0
\r
482 noiseROC = noiseTPC->GetCalROC(iSec); // noise per given sector
\r
484 Int_t nRows = 0; //number of rows in sector
\r
485 Int_t nDDLs = 0; //number of DDLs
\r
486 Int_t indexDDL = 0; //DDL index
\r
488 nRows = fParam->GetNRowLow();
\r
490 indexDDL = iSec * 2;
\r
492 nRows = fParam->GetNRowUp();
\r
494 indexDDL = (iSec-kNIS) * 4 + kNIS * 2;
\r
498 // Load the raw data for corresponding DDLs
\r
500 rawReader->Reset();
\r
501 rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1);
\r
504 while (input.NextDDL()){
\r
505 // Allocate memory for rows in sector (pads(depends on row) x timebins)
\r
506 if (!digarr->GetRow(iSec,0)){
\r
507 for(Int_t iRow = 0; iRow < nRows; iRow++) {
\r
508 digarr->CreateRow(iSec,iRow);
\r
509 }//end loop over rows
\r
512 while ( input.NextChannel() ) {
\r
513 Int_t iRow = input.GetRow();
\r
514 Int_t iPad = input.GetPad();
\r
515 //check row consistency
\r
516 if (iRow < 0 ) continue;
\r
517 if (iRow < 0 || iRow >= nRows){
\r
518 AliError(Form("Pad-row index (%d) outside the range (%d -> %d) !",
\r
519 iRow, 0, nRows -1));
\r
523 //check pad consistency
\r
524 if (iPad < 0 || iPad >= (Int_t)(roc->GetNPads(iSec,iRow))) {
\r
525 AliError(Form("Pad index (%d) outside the range (%d -> %d) !",
\r
526 iPad, 0, roc->GetNPads(iSec,iRow) ));
\r
530 //loop over bunches
\r
531 while ( input.NextBunch() ){
\r
532 Int_t startTbin = (Int_t)input.GetStartTimeBin();
\r
533 Int_t bunchlength = (Int_t)input.GetBunchLength();
\r
534 const UShort_t *sig = input.GetSignals();
\r
536 for (Int_t iTime = 0; iTime<bunchlength; iTime++){
\r
537 Int_t iTimeBin=startTbin-iTime;
\r
539 if(fDebugLevel==72){
\r
540 fHistoRow->Fill(iRow);
\r
541 fHistoPad->Fill(iPad);
\r
542 fHistoTime->Fill(iTimeBin);
\r
543 fHistoRowPad->Fill(iPad,iRow);
\r
544 }else if(fDebugLevel>=0&&fDebugLevel<72){
\r
545 if(iSec==fDebugLevel){
\r
546 fHistoRow->Fill(iRow);
\r
547 fHistoPad->Fill(iPad);
\r
548 fHistoTime->Fill(iTimeBin);
\r
549 fHistoRowPad->Fill(iPad,iRow);
\r
551 }else if(fDebugLevel==73){
\r
553 fHistoRow->Fill(iRow);
\r
554 fHistoPad->Fill(iPad);
\r
555 fHistoTime->Fill(iTimeBin);
\r
556 fHistoRowPad->Fill(iPad,iRow);
\r
558 }else if(fDebugLevel==74){
\r
560 fHistoRow->Fill(iRow);
\r
561 fHistoPad->Fill(iPad);
\r
562 fHistoTime->Fill(iTimeBin);
\r
563 fHistoRowPad->Fill(iPad,iRow);
\r
567 //check time consistency
\r
568 if ( iTimeBin < fRecoParam->GetFirstBin() || iTimeBin >= fRecoParam->GetLastBin()){
\r
569 //cout<<iTimeBin<<endl;
\r
571 AliFatal(Form("Timebin index (%d) outside the range (%d -> %d) !",
\r
572 iTimeBin, 0, fRecoParam->GetLastBin() -1));
\r
575 Float_t signal=(Float_t)sig[iTime];
\r
576 if (signal <= fZeroSup ||
\r
577 iTimeBin < fFirstBin ||
\r
578 iTimeBin > fLastBin
\r
580 digarr->GetRow(iSec,iRow)->SetDigitFast(0,iTimeBin,iPad);
\r
583 if (!noiseROC) continue;
\r
584 Double_t noiseOnPad = noiseROC->GetValue(iRow,iPad);//noise on given pad and row in sector
\r
585 if (noiseOnPad > fMaxNoiseAbs){
\r
586 digarr->GetRow(iSec,iRow)->SetDigitFast(0,iTimeBin,iPad);
\r
587 continue; // consider noisy pad as dead
\r
589 if(signal <= fMaxNoiseSigma * noiseOnPad){
\r
590 digarr->GetRow(iSec,iRow)->SetDigitFast(0,iTimeBin,iPad);
\r
593 digarr->GetRow(iSec,iRow)->SetDigitFast(TMath::Nint(signal),iTimeBin,iPad);
\r
594 }// end loop signals in bunch
\r
595 }// end loop bunches
\r
598 }// end sector loop
\r
600 if(isAltro) FindClusterKrIO();
\r
610 Int_t AliTPCclustererKr::FinderIOold(AliRawReader* rawReader)
\r
612 // Krypton cluster finder for the TPC raw data
\r
614 // fParam must be defined before
\r
616 if(rawReader)fRawData=kTRUE; //set flag to data
\r
619 Error("Digits2Clusters", "output tree not initialised");
\r
623 fParam->SetMaxTBin(fRecoParam->GetLastBin());//set number of timebins from reco -> param
\r
624 // used later for memory allocation
\r
626 AliRawEventHeaderBase* eventHeader = (AliRawEventHeaderBase*)rawReader->GetEventHeader();
\r
628 fTimeStamp = eventHeader->Get("Timestamp");
\r
631 Bool_t isAltro=kFALSE;
\r
633 AliTPCROC * roc = AliTPCROC::Instance();
\r
634 AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
\r
635 AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping();
\r
637 AliTPCRawStream input(rawReader,(AliAltroMapping**)mapping);
\r
639 const Int_t kNIS = fParam->GetNInnerSector();//number of inner sectors
\r
640 const Int_t kNOS = fParam->GetNOuterSector();//number of outer sectors
\r
641 const Int_t kNS = kNIS + kNOS;//all sectors
\r
645 AliTPCDigitsArray *digarr=new AliTPCDigitsArray(kFALSE);//data not sim
\r
646 digarr->Setup(fParam);//as usually parameters
\r
649 // Loop over sectors
\r
651 for(Int_t iSec = 0; iSec < kNS; iSec++) {
\r
652 AliTPCCalROC * noiseROC;
\r
654 noiseROC =new AliTPCCalROC(iSec);//noise=0
\r
657 noiseROC = noiseTPC->GetCalROC(iSec); // noise per given sector
\r
659 Int_t nRows = 0; //number of rows in sector
\r
660 Int_t nDDLs = 0; //number of DDLs
\r
661 Int_t indexDDL = 0; //DDL index
\r
663 nRows = fParam->GetNRowLow();
\r
665 indexDDL = iSec * 2;
\r
667 nRows = fParam->GetNRowUp();
\r
669 indexDDL = (iSec-kNIS) * 4 + kNIS * 2;
\r
673 // Load the raw data for corresponding DDLs
\r
675 rawReader->Reset();
\r
676 rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1);
\r
680 // Allocate memory for rows in sector (pads(depends on row) x timebins)
\r
681 for(Int_t iRow = 0; iRow < nRows; iRow++) {
\r
682 digarr->CreateRow(iSec,iRow);
\r
683 }//end loop over rows
\r
685 rawReader->Reset();
\r
686 rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1);
\r
689 // Begin loop over altro data
\r
691 while (input.Next()) {
\r
693 //check sector consistency
\r
694 if (input.GetSector() != iSec)
\r
695 AliFatal(Form("Sector index mismatch ! Expected (%d), but got (%d) !",iSec,input.GetSector()));
\r
697 Int_t iRow = input.GetRow();
\r
698 Int_t iPad = input.GetPad();
\r
699 Int_t iTimeBin = input.GetTime();
\r
702 if(fDebugLevel==72){
\r
703 fHistoRow->Fill(iRow);
\r
704 fHistoPad->Fill(iPad);
\r
705 fHistoTime->Fill(iTimeBin);
\r
706 fHistoRowPad->Fill(iPad,iRow);
\r
707 }else if(fDebugLevel>=0&&fDebugLevel<72){
\r
708 if(iSec==fDebugLevel){
\r
709 fHistoRow->Fill(iRow);
\r
710 fHistoPad->Fill(iPad);
\r
711 fHistoTime->Fill(iTimeBin);
\r
712 fHistoRowPad->Fill(iPad,iRow);
\r
714 }else if(fDebugLevel==73){
\r
716 fHistoRow->Fill(iRow);
\r
717 fHistoPad->Fill(iPad);
\r
718 fHistoTime->Fill(iTimeBin);
\r
719 fHistoRowPad->Fill(iPad,iRow);
\r
721 }else if(fDebugLevel==74){
\r
723 fHistoRow->Fill(iRow);
\r
724 fHistoPad->Fill(iPad);
\r
725 fHistoTime->Fill(iTimeBin);
\r
726 fHistoRowPad->Fill(iPad,iRow);
\r
730 //check row consistency
\r
731 if (iRow < 0 ) continue;
\r
732 if (iRow < 0 || iRow >= nRows){
\r
733 AliError(Form("Pad-row index (%d) outside the range (%d -> %d) !",
\r
734 iRow, 0, nRows -1));
\r
738 //check pad consistency
\r
739 if (iPad < 0 || iPad >= (Int_t)(roc->GetNPads(iSec,iRow))) {
\r
740 AliError(Form("Pad index (%d) outside the range (%d -> %d) !",
\r
741 iPad, 0, roc->GetNPads(iSec,iRow) ));
\r
745 //check time consistency
\r
746 if ( iTimeBin < fRecoParam->GetFirstBin() || iTimeBin >= fRecoParam->GetLastBin()){
\r
747 //cout<<iTimeBin<<endl;
\r
749 AliFatal(Form("Timebin index (%d) outside the range (%d -> %d) !",
\r
750 iTimeBin, 0, fRecoParam->GetLastBin() -1));
\r
754 Int_t signal = input.GetSignal();
\r
755 if (signal <= fZeroSup ||
\r
756 iTimeBin < fFirstBin ||
\r
757 iTimeBin > fLastBin
\r
759 digarr->GetRow(iSec,iRow)->SetDigitFast(0,iTimeBin,iPad);
\r
762 if (!noiseROC) continue;
\r
763 Double_t noiseOnPad = noiseROC->GetValue(iRow,iPad);//noise on given pad and row in sector
\r
764 if (noiseOnPad > fMaxNoiseAbs){
\r
765 digarr->GetRow(iSec,iRow)->SetDigitFast(0,iTimeBin,iPad);
\r
766 continue; // consider noisy pad as dead
\r
768 if(signal <= fMaxNoiseSigma * noiseOnPad){
\r
769 digarr->GetRow(iSec,iRow)->SetDigitFast(0,iTimeBin,iPad);
\r
772 digarr->GetRow(iSec,iRow)->SetDigitFast(signal,iTimeBin,iPad);
\r
773 }//end of loop over altro data
\r
774 }//end of loop over sectors
\r
777 if(isAltro) FindClusterKrIO();
\r
783 void AliTPCclustererKr::CleanSector(Int_t sector){
\r
785 // clean isolated digits
\r
787 const Int_t kNRows=fParam->GetNRow(sector);//number of rows in sector
\r
788 for(Int_t iRow=0; iRow<kNRows; ++iRow){
\r
789 AliSimDigits *digrow;
\r
791 digrow = (AliSimDigits*)fDigarr->GetRow(sector,iRow);//real data
\r
793 digrow = (AliSimDigits*)fDigarr->LoadRow(sector,iRow);//MC
\r
795 if(!digrow) continue;
\r
796 digrow->ExpandBuffer(); //decrunch
\r
797 const Int_t kNPads = digrow->GetNCols(); // number of pads
\r
798 const Int_t kNTime = digrow->GetNRows(); // number of timebins
\r
799 for(Int_t iPad=1;iPad<kNPads-1;iPad++){
\r
800 Short_t* val = digrow->GetDigitsColumn(iPad);
\r
802 for(Int_t iTimeBin=1;iTimeBin<kNTime-1;iTimeBin++){
\r
803 if (val[iTimeBin]<=0) continue;
\r
804 if (val[iTimeBin-1]+val[iTimeBin+1]<fIsolCut) {val[iTimeBin]=0; continue;}
\r
805 if (val[iTimeBin-kNTime]+val[iTimeBin+kNTime]<fIsolCut) {val[iTimeBin]=0; continue;}
\r
807 if (val[iTimeBin-1-kNTime]+val[iTimeBin+1+kNTime]<fIsolCut) {val[iTimeBin]=0; continue;}
\r
808 if (val[iTimeBin+1-kNTime]+val[iTimeBin-1+kNTime]<fIsolCut) {val[iTimeBin]=0; continue;}
\r
816 ////____________________________________________________________________________
\r
817 Int_t AliTPCclustererKr::FindClusterKrIO()
\r
821 //fParam and fDigarr must be set to run this method
\r
824 Int_t clusterCounter=0;
\r
825 const Int_t nTotalSector=fParam->GetNSector();//number of sectors
\r
826 for(Int_t iSec=0; iSec<nTotalSector; ++iSec){
\r
829 //vector of maxima for each sector
\r
830 //std::vector<AliPadMax*> maximaInSector;
\r
831 TObjArray *maximaInSector=new TObjArray();//to store AliPadMax*
\r
834 // looking for the maxima on the pad
\r
837 const Int_t kNRows=fParam->GetNRow(iSec);//number of rows in sector
\r
838 for(Int_t iRow=0; iRow<kNRows; ++iRow){
\r
839 AliSimDigits *digrow;
\r
841 digrow = (AliSimDigits*)fDigarr->GetRow(iSec,iRow);//real data
\r
843 digrow = (AliSimDigits*)fDigarr->LoadRow(iSec,iRow);//MC
\r
845 if(digrow){//if pointer exist
\r
846 digrow->ExpandBuffer(); //decrunch
\r
847 const Int_t kNPads = digrow->GetNCols(); // number of pads
\r
848 const Int_t kNTime = digrow->GetNRows(); // number of timebins
\r
849 for(Int_t iPad=0;iPad<kNPads;iPad++){
\r
851 Int_t timeBinMax=-1;//timebin of maximum
\r
852 Int_t valueMaximum=-1;//value of maximum in adc
\r
853 Int_t increaseBegin=-1;//timebin when increase starts
\r
854 Int_t sumAdc=0;//sum of adc on the pad in maximum surrounding
\r
855 bool ifIncreaseBegin=true;//flag - check if increasing started
\r
856 bool ifMaximum=false;//flag - check if it could be maximum
\r
857 Short_t* val = digrow->GetDigitsColumn(iPad);
\r
858 for(Int_t iTimeBin=1;iTimeBin<kNTime-1;iTimeBin++){
\r
860 if (val[iTimeBin]==-1) break; // 0 until the end
\r
861 for( ; iTimeBin<kNTime-2&&val[iTimeBin]<fMinAdc ;iTimeBin++) {}
\r
864 Short_t adc = val[iTimeBin];
\r
866 if(adc<fMinAdc){//standard was 3 for fMinAdc
\r
868 if(iTimeBin-increaseBegin<fMinTimeBins){//at least 2 time bins
\r
873 ifIncreaseBegin=true;
\r
877 //insert maximum, default values and set flags
\r
878 //Double_t xCord,yCord;
\r
879 //GetXY(iSec,iRow,iPad,xCord,yCord);
\r
880 Double_t x[]={iRow,iPad,iTimeBin};
\r
882 AliTPCTransform trafo;
\r
883 trafo.Transform(x,i,0,1);
\r
885 AliPadMax *oneMaximum = new AliPadMax(AliTPCvtpr(valueMaximum,
\r
891 x[2]/*timeBinMax*/),
\r
895 maximaInSector->AddLast(oneMaximum);
\r
901 ifIncreaseBegin=true;
\r
912 if(ifIncreaseBegin){
\r
913 ifIncreaseBegin=false;
\r
914 increaseBegin=iTimeBin;
\r
917 if(adc>valueMaximum){
\r
918 timeBinMax=iTimeBin;
\r
923 if(iTimeBin==kNTime-1 && ifMaximum && kNTime-increaseBegin>fMinTimeBins){//on the edge
\r
924 //at least 3 timebins
\r
925 //insert maximum, default values and set flags
\r
926 //Double_t xCord,yCord;
\r
927 //GetXY(iSec,iRow,iPad,xCord,yCord);
\r
928 Double_t x[]={iRow,iPad,iTimeBin};
\r
930 AliTPCTransform trafo;
\r
931 trafo.Transform(x,i,0,1);
\r
932 AliPadMax *oneMaximum = new AliPadMax(AliTPCvtpr(valueMaximum,
\r
938 x[2]/*timeBinMax*/),
\r
942 maximaInSector->AddLast(oneMaximum);
\r
948 ifIncreaseBegin=true;
\r
953 }//end loop over timebins
\r
954 }//end loop over pads
\r
956 // cout<<"Pointer does not exist!!"<<endl;
\r
957 }//end if poiner exists
\r
958 }//end loop over rows
\r
960 MakeClusters(maximaInSector,iSec,clusterCounter);
\r
962 maximaInSector->SetOwner(kTRUE);
\r
963 maximaInSector->Delete();
\r
964 delete maximaInSector;
\r
966 cout<<"Number of clusters in event: "<<clusterCounter<<endl;
\r
970 void AliTPCclustererKr::MakeClusters(TObjArray * maximaInSector, Int_t iSec, Int_t &clusterCounter){
\r
977 Int_t maxTimeBin=0;
\r
983 Int_t entriesArr = maximaInSector->GetEntriesFast();
\r
984 for(Int_t it1 = 0; it1 < entriesArr; ++it1 ) {
\r
986 AliPadMax *mp1=(AliPadMax *)maximaInSector->UncheckedAt(it1);
\r
987 if (!mp1) continue;
\r
988 AliTPCclusterKr clusterKr;
\r
991 Int_t clusterValue=0;
\r
992 clusterValue+=(mp1)->GetSum();
\r
993 list<Int_t> nUsedRows;
\r
994 nUsedRows.push_back((mp1)->GetRow());
\r
996 maxDig =(mp1)->GetAdc() ;
\r
997 maxSumAdc =(mp1)->GetSum() ;
\r
998 maxTimeBin =(mp1)->GetTime();
\r
999 maxPad =(mp1)->GetPad() ;
\r
1000 maxRow =(mp1)->GetRow() ;
\r
1001 maxX =(mp1)->GetX();
\r
1002 maxY =(mp1)->GetY();
\r
1003 maxT =(mp1)->GetT();
\r
1005 AliSimDigits *digrowTmp;
\r
1007 digrowTmp = (AliSimDigits*)fDigarr->GetRow(iSec,(mp1)->GetRow());
\r
1009 digrowTmp = (AliSimDigits*)fDigarr->LoadRow(iSec,(mp1)->GetRow());
\r
1012 digrowTmp->ExpandBuffer(); //decrunch
\r
1014 for(Int_t itb=(mp1)->GetBegin(); itb<((mp1)->GetEnd())+1; itb++){
\r
1015 Int_t adcTmp = digrowTmp->GetDigitUnchecked(itb,(mp1)->GetPad());
\r
1016 AliTPCvtpr *vtpr=new AliTPCvtpr(adcTmp,itb,(mp1)->GetPad(),(mp1)->GetRow(),(mp1)->GetX(),(mp1)->GetY(),(mp1)->GetT());
\r
1017 clusterKr.AddDigitToCluster(vtpr);
\r
1019 clusterKr.SetCenter();//set centr of the cluster
\r
1021 for(Int_t it2 = it1+1; it2 < entriesArr; ++it2 ) {
\r
1022 AliPadMax *mp2=(AliPadMax *)maximaInSector->UncheckedAt(it2);
\r
1023 if (!mp2) continue;
\r
1024 if (TMath::Abs(clusterKr.GetCenterX() - (mp2)->GetX()) > fMaxPadRangeCm) continue;
\r
1025 if (TMath::Abs(clusterKr.GetCenterY() - (mp2)->GetY()) > fMaxRowRangeCm) continue;
\r
1026 if (TMath::Abs(clusterKr.GetCenterT() - (mp2)->GetT()) > fMaxTimeRange) continue;
\r
1029 clusterValue+=(mp2)->GetSum();
\r
1032 nUsedRows.push_back((mp2)->GetRow());
\r
1034 AliSimDigits *digrowTmp1;
\r
1036 digrowTmp1 = (AliSimDigits*)fDigarr->GetRow(iSec,(mp2)->GetRow());
\r
1038 digrowTmp1 = (AliSimDigits*)fDigarr->LoadRow(iSec,(mp2)->GetRow());
\r
1040 digrowTmp1->ExpandBuffer(); //decrunch
\r
1042 for(Int_t itb=(mp2)->GetBegin(); itb<(mp2)->GetEnd()+1; itb++){
\r
1043 Int_t adcTmp = digrowTmp1->GetDigitUnchecked(itb,(mp2)->GetPad());
\r
1044 AliTPCvtpr *vtpr=new AliTPCvtpr(adcTmp,itb,(mp2)->GetPad(),(mp2)->GetRow(),(mp2)->GetX(),(mp2)->GetY(),(mp2)->GetT());
\r
1045 clusterKr.AddDigitToCluster(vtpr);
\r
1048 clusterKr.SetCenter();//set center of the cluster
\r
1050 //which one is bigger
\r
1051 if( (mp2)->GetAdc() > maxDig ){
\r
1052 maxDig =(mp2)->GetAdc() ;
\r
1053 maxSumAdc =(mp2)->GetSum() ;
\r
1054 maxTimeBin =(mp2)->GetTime();
\r
1055 maxPad =(mp2)->GetPad() ;
\r
1056 maxRow =(mp2)->GetRow() ;
\r
1057 maxX =(mp2)->GetX() ;
\r
1058 maxY =(mp2)->GetY() ;
\r
1059 maxT =(mp2)->GetT() ;
\r
1060 } else if ( (mp2)->GetAdc() == maxDig ){
\r
1061 if( (mp2)->GetSum() > maxSumAdc){
\r
1062 maxDig =(mp2)->GetAdc() ;
\r
1063 maxSumAdc =(mp2)->GetSum() ;
\r
1064 maxTimeBin =(mp2)->GetTime();
\r
1065 maxPad =(mp2)->GetPad() ;
\r
1066 maxRow =(mp2)->GetRow() ;
\r
1067 maxX =(mp2)->GetX() ;
\r
1068 maxY =(mp2)->GetY() ;
\r
1069 maxT =(mp2)->GetT() ;
\r
1072 delete maximaInSector->RemoveAt(it2);
\r
1075 delete maximaInSector->RemoveAt(it1);
\r
1076 clusterKr.SetSize();
\r
1077 //through out clusters on the edge and noise
\r
1078 //if(clusterValue/clusterKr.fCluster.size()<fValueToSize)continue;
\r
1079 if(clusterValue/(clusterKr.GetSize())<fValueToSize)continue;
\r
1081 clusterKr.SetADCcluster(clusterValue);
\r
1082 clusterKr.SetNPads(nUsedPads);
\r
1083 clusterKr.SetMax(AliTPCvtpr(maxDig,maxTimeBin,maxPad,maxRow,maxX,maxY,maxT));
\r
1084 clusterKr.SetSec(iSec);
\r
1085 clusterKr.SetSize();
\r
1088 nUsedRows.unique();
\r
1089 clusterKr.SetNRows(nUsedRows.size());
\r
1090 clusterKr.SetCenter();
\r
1092 clusterKr.SetRMS();//Set pad,row,timebin RMS
\r
1093 clusterKr.Set1D();//Set size in pads and timebins
\r
1095 clusterKr.SetTimeStamp(fTimeStamp);
\r
1100 //save each cluster into file
\r
1102 (*fOutput)<<"Kr"<<
\r
1103 "Cl.="<<&clusterKr<<
\r
1106 //end of save each cluster into file adc.root
\r
1112 ////____________________________________________________________________________
\r
1115 void AliTPCclustererKr::GetXY(Int_t sec,Int_t row,Int_t pad,Double_t& xGlob,Double_t& yGlob){
\r
1117 //gives global XY coordinate of the pad
\r
1120 Double_t yLocal = fParam->GetPadRowRadii(sec,row);//radius of row in sector in cm
\r
1122 Int_t padmax = fParam->GetNPads(sec,row);//number of pads in a given row
\r
1124 if(sec<fParam->GetNInnerSector())padXSize=0.4;
\r
1125 else padXSize=0.6;
\r
1126 Double_t xLocal=(pad+0.5-padmax/2.)*padXSize;//x-value of the center of pad
\r
1129 fParam->AdjustCosSin((Int_t)sec,cos,sin);//return sinus and cosinus of the sector
\r
1131 Double_t xGlob1 = xLocal * cos + yLocal * sin;
\r
1132 Double_t yGlob1 = -xLocal * sin + yLocal * cos;
\r
1136 rot=TMath::Pi()/2.;
\r
1138 xGlob = xGlob1 * TMath::Cos(rot) + yGlob1 * TMath::Sin(rot);
\r
1139 yGlob = -xGlob1 * TMath::Sin(rot) + yGlob1 * TMath::Cos(rot);
\r
1142 if(sec<18||(sec>=36&&sec<54)) xGlob =-1*xGlob;
\r