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 if (this == ¶m) return (*this);
\r
344 fParam = param.fParam;
\r
345 fRecoParam = param.fRecoParam;
\r
346 fRawData = param.fRawData;
\r
347 fInput = param.fInput ;
\r
348 fOutput = param.fOutput;
\r
349 fDigarr = param.fDigarr;
\r
350 fZeroSup = param.fZeroSup ;
\r
351 fFirstBin = param.fFirstBin ;
\r
352 fLastBin = param.fLastBin ;
\r
353 fMaxNoiseAbs = param.fMaxNoiseAbs ;
\r
354 fMaxNoiseSigma = param.fMaxNoiseSigma ;
\r
355 fMinAdc = param.fMinAdc;
\r
356 fMinTimeBins = param.fMinTimeBins;
\r
357 // fMaxPadRange = param.fMaxPadRange ;
\r
358 // fMaxRowRange = param.fMaxRowRange ;
\r
359 fMaxTimeRange = param.fMaxTimeRange;
\r
360 fValueToSize = param.fValueToSize;
\r
361 fMaxPadRangeCm = param.fMaxPadRangeCm;
\r
362 fMaxRowRangeCm = param.fMaxRowRangeCm;
\r
363 fIsolCut = param.fIsolCut;
\r
364 fDebugLevel = param.fDebugLevel;
\r
365 fHistoRow = param.fHistoRow ;
\r
366 fHistoPad = param.fHistoPad ;
\r
367 fHistoTime = param.fHistoTime;
\r
368 fHistoRowPad = param.fHistoRowPad;
\r
369 fTimeStamp = param.fTimeStamp;
\r
374 AliTPCclustererKr::~AliTPCclustererKr()
\r
382 void AliTPCclustererKr::SetRecoParam(AliTPCRecoParam *recoParam)
\r
385 // set reconstruction parameters
\r
388 fRecoParam = recoParam;
\r
390 //set default parameters if not specified
\r
391 fRecoParam = AliTPCReconstructor::GetRecoParam();
\r
392 if (!fRecoParam) fRecoParam = AliTPCRecoParam::GetLowFluxParam();
\r
398 ////____________________________________________________________________________
\r
400 void AliTPCclustererKr::SetInput(TTree * tree)
\r
403 // set input tree with digits
\r
406 if (!fInput->GetBranch("Segment")){
\r
407 cerr<<"AliTPCclusterKr::FindClusterKr(): no proper input tree !\n";
\r
413 void AliTPCclustererKr::SetOutput(TTree * /*tree*/)
\r
418 fOutput = new TTreeSRedirector("Krypton.root");
\r
421 ////____________________________________________________________________________
\r
423 Int_t AliTPCclustererKr::FinderIO()
\r
425 // Krypton cluster finder for simulated events from MC
\r
428 Error("Digits2Clusters", "input tree not initialised");
\r
433 Error("Digits2Clusters", "output tree not initialised");
\r
443 Int_t AliTPCclustererKr::FinderIO(AliRawReader* rawReader)
\r
445 // Krypton cluster finder for the TPC raw data
\r
446 // this method is unsing AliAltroRawStreamV3
\r
447 // fParam must be defined before
\r
448 if (!rawReader) return 1;
\r
450 fRawData=kTRUE; //set flag to data
\r
453 Error("Digits2Clusters", "output tree not initialised");
\r
457 fParam->SetMaxTBin(fRecoParam->GetLastBin());//set number of timebins from reco -> param
\r
458 // used later for memory allocation
\r
460 AliRawEventHeaderBase* eventHeader = (AliRawEventHeaderBase*)rawReader->GetEventHeader();
\r
462 fTimeStamp = eventHeader->Get("Timestamp");
\r
463 fRun = rawReader->GetRunNumber();
\r
467 Bool_t isAltro=kFALSE;
\r
469 AliTPCROC * roc = AliTPCROC::Instance();
\r
470 AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
\r
471 AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping();
\r
473 AliTPCRawStreamV3 input(rawReader,(AliAltroMapping**)mapping);
\r
475 const Int_t kNIS = fParam->GetNInnerSector();//number of inner sectors
\r
476 const Int_t kNOS = fParam->GetNOuterSector();//number of outer sectors
\r
477 const Int_t kNS = kNIS + kNOS;//all sectors
\r
481 AliTPCDigitsArray *digarr=new AliTPCDigitsArray(kFALSE);//data not sim
\r
482 digarr->Setup(fParam);//as usually parameters
\r
484 for(Int_t iSec = 0; iSec < kNS; iSec++) {
\r
485 AliTPCCalROC * noiseROC;
\r
486 AliTPCCalROC noiseDummy(iSec);
\r
488 noiseROC = &noiseDummy;//noise=0
\r
490 noiseROC = noiseTPC->GetCalROC(iSec); // noise per given sector
\r
492 Int_t nRows = 0; //number of rows in sector
\r
493 Int_t nDDLs = 0; //number of DDLs
\r
494 Int_t indexDDL = 0; //DDL index
\r
496 nRows = fParam->GetNRowLow();
\r
498 indexDDL = iSec * 2;
\r
500 nRows = fParam->GetNRowUp();
\r
502 indexDDL = (iSec-kNIS) * 4 + kNIS * 2;
\r
506 // Load the raw data for corresponding DDLs
\r
508 rawReader->Reset();
\r
509 rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1);
\r
512 while (input.NextDDL()){
\r
513 // Allocate memory for rows in sector (pads(depends on row) x timebins)
\r
514 if (!digarr->GetRow(iSec,0)){
\r
515 for(Int_t iRow = 0; iRow < nRows; iRow++) {
\r
516 digarr->CreateRow(iSec,iRow);
\r
517 }//end loop over rows
\r
520 while ( input.NextChannel() ) {
\r
521 Int_t iRow = input.GetRow();
\r
522 Int_t iPad = input.GetPad();
\r
523 //check row consistency
\r
524 if (iRow < 0 ) continue;
\r
525 if (iRow < 0 || iRow >= nRows){
\r
526 AliError(Form("Pad-row index (%d) outside the range (%d -> %d) !",
\r
527 iRow, 0, nRows -1));
\r
531 //check pad consistency
\r
532 if (iPad < 0 || iPad >= (Int_t)(roc->GetNPads(iSec,iRow))) {
\r
533 AliError(Form("Pad index (%d) outside the range (%d -> %d) !",
\r
534 iPad, 0, roc->GetNPads(iSec,iRow) ));
\r
538 //loop over bunches
\r
539 while ( input.NextBunch() ){
\r
540 Int_t startTbin = (Int_t)input.GetStartTimeBin();
\r
541 Int_t bunchlength = (Int_t)input.GetBunchLength();
\r
542 const UShort_t *sig = input.GetSignals();
\r
544 for (Int_t iTime = 0; iTime<bunchlength; iTime++){
\r
545 Int_t iTimeBin=startTbin-iTime;
\r
547 if(fDebugLevel==72){
\r
548 fHistoRow->Fill(iRow);
\r
549 fHistoPad->Fill(iPad);
\r
550 fHistoTime->Fill(iTimeBin);
\r
551 fHistoRowPad->Fill(iPad,iRow);
\r
552 }else if(fDebugLevel>=0&&fDebugLevel<72){
\r
553 if(iSec==fDebugLevel){
\r
554 fHistoRow->Fill(iRow);
\r
555 fHistoPad->Fill(iPad);
\r
556 fHistoTime->Fill(iTimeBin);
\r
557 fHistoRowPad->Fill(iPad,iRow);
\r
559 }else if(fDebugLevel==73){
\r
561 fHistoRow->Fill(iRow);
\r
562 fHistoPad->Fill(iPad);
\r
563 fHistoTime->Fill(iTimeBin);
\r
564 fHistoRowPad->Fill(iPad,iRow);
\r
566 }else if(fDebugLevel==74){
\r
568 fHistoRow->Fill(iRow);
\r
569 fHistoPad->Fill(iPad);
\r
570 fHistoTime->Fill(iTimeBin);
\r
571 fHistoRowPad->Fill(iPad,iRow);
\r
575 //check time consistency
\r
576 if ( iTimeBin < fRecoParam->GetFirstBin() || iTimeBin >= fRecoParam->GetLastBin()){
\r
577 //cout<<iTimeBin<<endl;
\r
579 AliFatal(Form("Timebin index (%d) outside the range (%d -> %d) !",
\r
580 iTimeBin, 0, fRecoParam->GetLastBin() -1));
\r
583 Float_t signal=(Float_t)sig[iTime];
\r
584 if (signal <= fZeroSup ||
\r
585 iTimeBin < fFirstBin ||
\r
586 iTimeBin > fLastBin
\r
588 digarr->GetRow(iSec,iRow)->SetDigitFast(0,iTimeBin,iPad);
\r
591 if (!noiseROC) continue;
\r
592 Double_t noiseOnPad = noiseROC->GetValue(iRow,iPad);//noise on given pad and row in sector
\r
593 if (noiseOnPad > fMaxNoiseAbs){
\r
594 digarr->GetRow(iSec,iRow)->SetDigitFast(0,iTimeBin,iPad);
\r
595 continue; // consider noisy pad as dead
\r
597 if(signal <= fMaxNoiseSigma * noiseOnPad){
\r
598 digarr->GetRow(iSec,iRow)->SetDigitFast(0,iTimeBin,iPad);
\r
601 digarr->GetRow(iSec,iRow)->SetDigitFast(TMath::Nint(signal),iTimeBin,iPad);
\r
602 }// end loop signals in bunch
\r
603 }// end loop bunches
\r
606 }// end sector loop
\r
608 if(isAltro) FindClusterKrIO();
\r
618 Int_t AliTPCclustererKr::FinderIOold(AliRawReader* rawReader)
\r
620 // Krypton cluster finder for the TPC raw data
\r
622 // fParam must be defined before
\r
623 if (!rawReader) return 1;
\r
625 if(rawReader)fRawData=kTRUE; //set flag to data
\r
628 Error("Digits2Clusters", "output tree not initialised");
\r
632 fParam->SetMaxTBin(fRecoParam->GetLastBin());//set number of timebins from reco -> param
\r
633 // used later for memory allocation
\r
635 AliRawEventHeaderBase* eventHeader = (AliRawEventHeaderBase*)rawReader->GetEventHeader();
\r
637 fTimeStamp = eventHeader->Get("Timestamp");
\r
638 fRun = rawReader->GetRunNumber();
\r
641 Bool_t isAltro=kFALSE;
\r
643 AliTPCROC * roc = AliTPCROC::Instance();
\r
644 AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
\r
645 AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping();
\r
647 AliTPCRawStream input(rawReader,(AliAltroMapping**)mapping);
\r
649 const Int_t kNIS = fParam->GetNInnerSector();//number of inner sectors
\r
650 const Int_t kNOS = fParam->GetNOuterSector();//number of outer sectors
\r
651 const Int_t kNS = kNIS + kNOS;//all sectors
\r
655 AliTPCDigitsArray *digarr=new AliTPCDigitsArray(kFALSE);//data not sim
\r
656 digarr->Setup(fParam);//as usually parameters
\r
659 // Loop over sectors
\r
661 for(Int_t iSec = 0; iSec < kNS; iSec++) {
\r
662 AliTPCCalROC * noiseROC;
\r
664 noiseROC =new AliTPCCalROC(iSec);//noise=0
\r
667 noiseROC = noiseTPC->GetCalROC(iSec); // noise per given sector
\r
669 Int_t nRows = 0; //number of rows in sector
\r
670 Int_t nDDLs = 0; //number of DDLs
\r
671 Int_t indexDDL = 0; //DDL index
\r
673 nRows = fParam->GetNRowLow();
\r
675 indexDDL = iSec * 2;
\r
677 nRows = fParam->GetNRowUp();
\r
679 indexDDL = (iSec-kNIS) * 4 + kNIS * 2;
\r
683 // Load the raw data for corresponding DDLs
\r
685 rawReader->Reset();
\r
686 rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1);
\r
690 // Allocate memory for rows in sector (pads(depends on row) x timebins)
\r
691 for(Int_t iRow = 0; iRow < nRows; iRow++) {
\r
692 digarr->CreateRow(iSec,iRow);
\r
693 }//end loop over rows
\r
695 rawReader->Reset();
\r
696 rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1);
\r
699 // Begin loop over altro data
\r
701 while (input.Next()) {
\r
703 //check sector consistency
\r
704 if (input.GetSector() != iSec)
\r
705 AliFatal(Form("Sector index mismatch ! Expected (%d), but got (%d) !",iSec,input.GetSector()));
\r
707 Int_t iRow = input.GetRow();
\r
708 Int_t iPad = input.GetPad();
\r
709 Int_t iTimeBin = input.GetTime();
\r
712 if(fDebugLevel==72){
\r
713 fHistoRow->Fill(iRow);
\r
714 fHistoPad->Fill(iPad);
\r
715 fHistoTime->Fill(iTimeBin);
\r
716 fHistoRowPad->Fill(iPad,iRow);
\r
717 }else if(fDebugLevel>=0&&fDebugLevel<72){
\r
718 if(iSec==fDebugLevel){
\r
719 fHistoRow->Fill(iRow);
\r
720 fHistoPad->Fill(iPad);
\r
721 fHistoTime->Fill(iTimeBin);
\r
722 fHistoRowPad->Fill(iPad,iRow);
\r
724 }else if(fDebugLevel==73){
\r
726 fHistoRow->Fill(iRow);
\r
727 fHistoPad->Fill(iPad);
\r
728 fHistoTime->Fill(iTimeBin);
\r
729 fHistoRowPad->Fill(iPad,iRow);
\r
731 }else if(fDebugLevel==74){
\r
733 fHistoRow->Fill(iRow);
\r
734 fHistoPad->Fill(iPad);
\r
735 fHistoTime->Fill(iTimeBin);
\r
736 fHistoRowPad->Fill(iPad,iRow);
\r
740 //check row consistency
\r
741 if (iRow < 0 ) continue;
\r
742 if (iRow < 0 || iRow >= nRows){
\r
743 AliError(Form("Pad-row index (%d) outside the range (%d -> %d) !",
\r
744 iRow, 0, nRows -1));
\r
748 //check pad consistency
\r
749 if (iPad < 0 || iPad >= (Int_t)(roc->GetNPads(iSec,iRow))) {
\r
750 AliError(Form("Pad index (%d) outside the range (%d -> %d) !",
\r
751 iPad, 0, roc->GetNPads(iSec,iRow) ));
\r
755 //check time consistency
\r
756 if ( iTimeBin < fRecoParam->GetFirstBin() || iTimeBin >= fRecoParam->GetLastBin()){
\r
757 //cout<<iTimeBin<<endl;
\r
759 AliFatal(Form("Timebin index (%d) outside the range (%d -> %d) !",
\r
760 iTimeBin, 0, fRecoParam->GetLastBin() -1));
\r
764 Int_t signal = input.GetSignal();
\r
765 if (signal <= fZeroSup ||
\r
766 iTimeBin < fFirstBin ||
\r
767 iTimeBin > fLastBin
\r
769 digarr->GetRow(iSec,iRow)->SetDigitFast(0,iTimeBin,iPad);
\r
772 if (!noiseROC) continue;
\r
773 Double_t noiseOnPad = noiseROC->GetValue(iRow,iPad);//noise on given pad and row in sector
\r
774 if (noiseOnPad > fMaxNoiseAbs){
\r
775 digarr->GetRow(iSec,iRow)->SetDigitFast(0,iTimeBin,iPad);
\r
776 continue; // consider noisy pad as dead
\r
778 if(signal <= fMaxNoiseSigma * noiseOnPad){
\r
779 digarr->GetRow(iSec,iRow)->SetDigitFast(0,iTimeBin,iPad);
\r
782 digarr->GetRow(iSec,iRow)->SetDigitFast(signal,iTimeBin,iPad);
\r
783 }//end of loop over altro data
\r
784 }//end of loop over sectors
\r
787 if(isAltro) FindClusterKrIO();
\r
793 void AliTPCclustererKr::CleanSector(Int_t sector){
\r
795 // clean isolated digits
\r
797 const Int_t kNRows=fParam->GetNRow(sector);//number of rows in sector
\r
798 for(Int_t iRow=0; iRow<kNRows; ++iRow){
\r
799 AliSimDigits *digrow;
\r
801 digrow = (AliSimDigits*)fDigarr->GetRow(sector,iRow);//real data
\r
803 digrow = (AliSimDigits*)fDigarr->LoadRow(sector,iRow);//MC
\r
805 if(!digrow) continue;
\r
806 digrow->ExpandBuffer(); //decrunch
\r
807 const Int_t kNPads = digrow->GetNCols(); // number of pads
\r
808 const Int_t kNTime = digrow->GetNRows(); // number of timebins
\r
809 for(Int_t iPad=1;iPad<kNPads-1;iPad++){
\r
810 Short_t* val = digrow->GetDigitsColumn(iPad);
\r
812 for(Int_t iTimeBin=1;iTimeBin<kNTime-1;iTimeBin++){
\r
813 if (val[iTimeBin]<=0) continue;
\r
814 if (val[iTimeBin-1]+val[iTimeBin+1]<fIsolCut) {val[iTimeBin]=0; continue;}
\r
815 if (val[iTimeBin-kNTime]+val[iTimeBin+kNTime]<fIsolCut) {val[iTimeBin]=0; continue;}
\r
817 if (val[iTimeBin-1-kNTime]+val[iTimeBin+1+kNTime]<fIsolCut) {val[iTimeBin]=0; continue;}
\r
818 if (val[iTimeBin+1-kNTime]+val[iTimeBin-1+kNTime]<fIsolCut) {val[iTimeBin]=0; continue;}
\r
826 ////____________________________________________________________________________
\r
827 Int_t AliTPCclustererKr::FindClusterKrIO()
\r
831 //fParam and fDigarr must be set to run this method
\r
834 Int_t clusterCounter=0;
\r
835 const Int_t nTotalSector=fParam->GetNSector();//number of sectors
\r
836 for(Int_t iSec=0; iSec<nTotalSector; ++iSec){
\r
839 //vector of maxima for each sector
\r
840 //std::vector<AliPadMax*> maximaInSector;
\r
841 TObjArray *maximaInSector=new TObjArray();//to store AliPadMax*
\r
844 // looking for the maxima on the pad
\r
847 const Int_t kNRows=fParam->GetNRow(iSec);//number of rows in sector
\r
848 for(Int_t iRow=0; iRow<kNRows; ++iRow){
\r
849 AliSimDigits *digrow;
\r
851 digrow = (AliSimDigits*)fDigarr->GetRow(iSec,iRow);//real data
\r
853 digrow = (AliSimDigits*)fDigarr->LoadRow(iSec,iRow);//MC
\r
855 if(digrow){//if pointer exist
\r
856 digrow->ExpandBuffer(); //decrunch
\r
857 const Int_t kNPads = digrow->GetNCols(); // number of pads
\r
858 const Int_t kNTime = digrow->GetNRows(); // number of timebins
\r
859 for(Int_t iPad=0;iPad<kNPads;iPad++){
\r
861 Int_t timeBinMax=-1;//timebin of maximum
\r
862 Int_t valueMaximum=-1;//value of maximum in adc
\r
863 Int_t increaseBegin=-1;//timebin when increase starts
\r
864 Int_t sumAdc=0;//sum of adc on the pad in maximum surrounding
\r
865 bool ifIncreaseBegin=true;//flag - check if increasing started
\r
866 bool ifMaximum=false;//flag - check if it could be maximum
\r
867 Short_t* val = digrow->GetDigitsColumn(iPad);
\r
868 for(Int_t iTimeBin=1;iTimeBin<kNTime-1;iTimeBin++){
\r
870 if (val[iTimeBin]==-1) break; // 0 until the end
\r
871 for( ; iTimeBin<kNTime-2&&val[iTimeBin]<fMinAdc ;iTimeBin++) {}
\r
874 Short_t adc = val[iTimeBin];
\r
876 if(adc<fMinAdc){//standard was 3 for fMinAdc
\r
878 if(iTimeBin-increaseBegin<fMinTimeBins){//at least 2 time bins
\r
883 ifIncreaseBegin=true;
\r
887 //insert maximum, default values and set flags
\r
888 //Double_t xCord,yCord;
\r
889 //GetXY(iSec,iRow,iPad,xCord,yCord);
\r
890 Double_t x[]={iRow,iPad,iTimeBin};
\r
892 AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
\r
894 transform->Transform(x,i,0,1);
\r
896 AliPadMax *oneMaximum = new AliPadMax(AliTPCvtpr(valueMaximum,
\r
902 x[2]/*timeBinMax*/),
\r
906 maximaInSector->AddLast(oneMaximum);
\r
912 ifIncreaseBegin=true;
\r
923 if(ifIncreaseBegin){
\r
924 ifIncreaseBegin=false;
\r
925 increaseBegin=iTimeBin;
\r
928 if(adc>valueMaximum){
\r
929 timeBinMax=iTimeBin;
\r
934 if(iTimeBin==kNTime-1 && ifMaximum && kNTime-increaseBegin>fMinTimeBins){//on the edge
\r
935 //at least 3 timebins
\r
936 //insert maximum, default values and set flags
\r
937 //Double_t xCord,yCord;
\r
938 //GetXY(iSec,iRow,iPad,xCord,yCord);
\r
939 Double_t x[]={iRow,iPad,iTimeBin};
\r
941 //AliTPCTransform trafo;
\r
942 //trafo.Transform(x,i,0,1);
\r
944 AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
\r
946 transform->Transform(x,i,0,1);
\r
948 AliPadMax *oneMaximum = new AliPadMax(AliTPCvtpr(valueMaximum,
\r
954 x[2]/*timeBinMax*/),
\r
958 maximaInSector->AddLast(oneMaximum);
\r
964 ifIncreaseBegin=true;
\r
969 }//end loop over timebins
\r
970 }//end loop over pads
\r
972 // cout<<"Pointer does not exist!!"<<endl;
\r
973 }//end if poiner exists
\r
974 }//end loop over rows
\r
976 MakeClusters(maximaInSector,iSec,clusterCounter);
\r
978 maximaInSector->SetOwner(kTRUE);
\r
979 maximaInSector->Delete();
\r
980 delete maximaInSector;
\r
982 cout<<"Number of clusters in event: "<<clusterCounter<<endl;
\r
986 void AliTPCclustererKr::MakeClusters(TObjArray * maximaInSector, Int_t iSec, Int_t &clusterCounter){
\r
993 Int_t maxTimeBin=0;
\r
999 Int_t entriesArr = maximaInSector->GetEntriesFast();
\r
1000 for(Int_t it1 = 0; it1 < entriesArr; ++it1 ) {
\r
1002 AliPadMax *mp1=(AliPadMax *)maximaInSector->UncheckedAt(it1);
\r
1003 if (!mp1) continue;
\r
1004 AliTPCclusterKr clusterKr;
\r
1006 Int_t nUsedPads=1;
\r
1007 Int_t clusterValue=0;
\r
1008 clusterValue+=(mp1)->GetSum();
\r
1009 list<Int_t> nUsedRows;
\r
1010 nUsedRows.push_back((mp1)->GetRow());
\r
1012 maxDig =(mp1)->GetAdc() ;
\r
1013 maxSumAdc =(mp1)->GetSum() ;
\r
1014 maxTimeBin =(mp1)->GetTime();
\r
1015 maxPad =(mp1)->GetPad() ;
\r
1016 maxRow =(mp1)->GetRow() ;
\r
1017 maxX =(mp1)->GetX();
\r
1018 maxY =(mp1)->GetY();
\r
1019 maxT =(mp1)->GetT();
\r
1021 AliSimDigits *digrowTmp;
\r
1023 digrowTmp = (AliSimDigits*)fDigarr->GetRow(iSec,(mp1)->GetRow());
\r
1025 digrowTmp = (AliSimDigits*)fDigarr->LoadRow(iSec,(mp1)->GetRow());
\r
1028 digrowTmp->ExpandBuffer(); //decrunch
\r
1030 for(Int_t itb=(mp1)->GetBegin(); itb<((mp1)->GetEnd())+1; itb++){
\r
1031 Int_t adcTmp = digrowTmp->GetDigitUnchecked(itb,(mp1)->GetPad());
\r
1032 AliTPCvtpr *vtpr=new AliTPCvtpr(adcTmp,itb,(mp1)->GetPad(),(mp1)->GetRow(),(mp1)->GetX(),(mp1)->GetY(),(mp1)->GetT());
\r
1033 clusterKr.AddDigitToCluster(vtpr);
\r
1035 clusterKr.SetCenter();//set centr of the cluster
\r
1037 for(Int_t it2 = it1+1; it2 < entriesArr; ++it2 ) {
\r
1038 AliPadMax *mp2=(AliPadMax *)maximaInSector->UncheckedAt(it2);
\r
1039 if (!mp2) continue;
\r
1040 if (TMath::Abs(clusterKr.GetCenterX() - (mp2)->GetX()) > fMaxPadRangeCm) continue;
\r
1041 if (TMath::Abs(clusterKr.GetCenterY() - (mp2)->GetY()) > fMaxRowRangeCm) continue;
\r
1042 if (TMath::Abs(clusterKr.GetCenterT() - (mp2)->GetT()) > fMaxTimeRange) continue;
\r
1045 clusterValue+=(mp2)->GetSum();
\r
1048 nUsedRows.push_back((mp2)->GetRow());
\r
1050 AliSimDigits *digrowTmp1;
\r
1052 digrowTmp1 = (AliSimDigits*)fDigarr->GetRow(iSec,(mp2)->GetRow());
\r
1054 digrowTmp1 = (AliSimDigits*)fDigarr->LoadRow(iSec,(mp2)->GetRow());
\r
1056 digrowTmp1->ExpandBuffer(); //decrunch
\r
1058 for(Int_t itb=(mp2)->GetBegin(); itb<(mp2)->GetEnd()+1; itb++){
\r
1059 Int_t adcTmp = digrowTmp1->GetDigitUnchecked(itb,(mp2)->GetPad());
\r
1060 AliTPCvtpr *vtpr=new AliTPCvtpr(adcTmp,itb,(mp2)->GetPad(),(mp2)->GetRow(),(mp2)->GetX(),(mp2)->GetY(),(mp2)->GetT());
\r
1061 clusterKr.AddDigitToCluster(vtpr);
\r
1064 clusterKr.SetCenter();//set center of the cluster
\r
1066 //which one is bigger
\r
1067 if( (mp2)->GetAdc() > maxDig ){
\r
1068 maxDig =(mp2)->GetAdc() ;
\r
1069 maxSumAdc =(mp2)->GetSum() ;
\r
1070 maxTimeBin =(mp2)->GetTime();
\r
1071 maxPad =(mp2)->GetPad() ;
\r
1072 maxRow =(mp2)->GetRow() ;
\r
1073 maxX =(mp2)->GetX() ;
\r
1074 maxY =(mp2)->GetY() ;
\r
1075 maxT =(mp2)->GetT() ;
\r
1076 } else if ( (mp2)->GetAdc() == maxDig ){
\r
1077 if( (mp2)->GetSum() > maxSumAdc){
\r
1078 maxDig =(mp2)->GetAdc() ;
\r
1079 maxSumAdc =(mp2)->GetSum() ;
\r
1080 maxTimeBin =(mp2)->GetTime();
\r
1081 maxPad =(mp2)->GetPad() ;
\r
1082 maxRow =(mp2)->GetRow() ;
\r
1083 maxX =(mp2)->GetX() ;
\r
1084 maxY =(mp2)->GetY() ;
\r
1085 maxT =(mp2)->GetT() ;
\r
1088 delete maximaInSector->RemoveAt(it2);
\r
1091 delete maximaInSector->RemoveAt(it1);
\r
1092 clusterKr.SetSize();
\r
1093 //through out clusters on the edge and noise
\r
1094 //if(clusterValue/clusterKr.fCluster.size()<fValueToSize)continue;
\r
1095 if(clusterValue/(clusterKr.GetSize())<fValueToSize)continue;
\r
1097 clusterKr.SetADCcluster(clusterValue);
\r
1098 clusterKr.SetNPads(nUsedPads);
\r
1099 clusterKr.SetMax(AliTPCvtpr(maxDig,maxTimeBin,maxPad,maxRow,maxX,maxY,maxT));
\r
1100 clusterKr.SetSec(iSec);
\r
1101 clusterKr.SetSize();
\r
1104 nUsedRows.unique();
\r
1105 clusterKr.SetNRows(nUsedRows.size());
\r
1106 clusterKr.SetCenter();
\r
1108 clusterKr.SetRMS();//Set pad,row,timebin RMS
\r
1109 clusterKr.Set1D();//Set size in pads and timebins
\r
1111 clusterKr.SetTimeStamp(fTimeStamp);
\r
1112 clusterKr.SetRun(fRun);
\r
1117 //save each cluster into file
\r
1119 (*fOutput)<<"Kr"<<
\r
1120 "Cl.="<<&clusterKr<<
\r
1123 //end of save each cluster into file adc.root
\r
1129 ////____________________________________________________________________________
\r
1132 void AliTPCclustererKr::GetXY(Int_t sec,Int_t row,Int_t pad,Double_t& xGlob,Double_t& yGlob){
\r
1134 //gives global XY coordinate of the pad
\r
1137 Double_t yLocal = fParam->GetPadRowRadii(sec,row);//radius of row in sector in cm
\r
1139 Int_t padmax = fParam->GetNPads(sec,row);//number of pads in a given row
\r
1141 if(sec<fParam->GetNInnerSector())padXSize=0.4;
\r
1142 else padXSize=0.6;
\r
1143 Double_t xLocal=(pad+0.5-padmax/2.)*padXSize;//x-value of the center of pad
\r
1146 fParam->AdjustCosSin((Int_t)sec,cos,sin);//return sinus and cosinus of the sector
\r
1148 Double_t xGlob1 = xLocal * cos + yLocal * sin;
\r
1149 Double_t yGlob1 = -xLocal * sin + yLocal * cos;
\r
1153 rot=TMath::Pi()/2.;
\r
1155 xGlob = xGlob1 * TMath::Cos(rot) + yGlob1 * TMath::Sin(rot);
\r
1156 yGlob = -xGlob1 * TMath::Sin(rot) + yGlob1 * TMath::Cos(rot);
\r
1159 if(sec<18||(sec>=36&&sec<54)) xGlob =-1*xGlob;
\r