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
447 if(rawReader)fRawData=kTRUE; //set flag to data
\r
450 Error("Digits2Clusters", "output tree not initialised");
\r
454 fParam->SetMaxTBin(fRecoParam->GetLastBin());//set number of timebins from reco -> param
\r
455 // used later for memory allocation
\r
457 AliRawEventHeaderBase* eventHeader = (AliRawEventHeaderBase*)rawReader->GetEventHeader();
\r
459 fTimeStamp = eventHeader->Get("Timestamp");
\r
460 fRun = rawReader->GetRunNumber();
\r
464 Bool_t isAltro=kFALSE;
\r
466 AliTPCROC * roc = AliTPCROC::Instance();
\r
467 AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
\r
468 AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping();
\r
470 AliTPCRawStreamV3 input(rawReader,(AliAltroMapping**)mapping);
\r
472 const Int_t kNIS = fParam->GetNInnerSector();//number of inner sectors
\r
473 const Int_t kNOS = fParam->GetNOuterSector();//number of outer sectors
\r
474 const Int_t kNS = kNIS + kNOS;//all sectors
\r
478 AliTPCDigitsArray *digarr=new AliTPCDigitsArray(kFALSE);//data not sim
\r
479 digarr->Setup(fParam);//as usually parameters
\r
481 for(Int_t iSec = 0; iSec < kNS; iSec++) {
\r
482 AliTPCCalROC * noiseROC;
\r
483 AliTPCCalROC noiseDummy(iSec);
\r
485 noiseROC = &noiseDummy;//noise=0
\r
487 noiseROC = noiseTPC->GetCalROC(iSec); // noise per given sector
\r
489 Int_t nRows = 0; //number of rows in sector
\r
490 Int_t nDDLs = 0; //number of DDLs
\r
491 Int_t indexDDL = 0; //DDL index
\r
493 nRows = fParam->GetNRowLow();
\r
495 indexDDL = iSec * 2;
\r
497 nRows = fParam->GetNRowUp();
\r
499 indexDDL = (iSec-kNIS) * 4 + kNIS * 2;
\r
503 // Load the raw data for corresponding DDLs
\r
505 rawReader->Reset();
\r
506 rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1);
\r
509 while (input.NextDDL()){
\r
510 // Allocate memory for rows in sector (pads(depends on row) x timebins)
\r
511 if (!digarr->GetRow(iSec,0)){
\r
512 for(Int_t iRow = 0; iRow < nRows; iRow++) {
\r
513 digarr->CreateRow(iSec,iRow);
\r
514 }//end loop over rows
\r
517 while ( input.NextChannel() ) {
\r
518 Int_t iRow = input.GetRow();
\r
519 Int_t iPad = input.GetPad();
\r
520 //check row consistency
\r
521 if (iRow < 0 ) continue;
\r
522 if (iRow < 0 || iRow >= nRows){
\r
523 AliError(Form("Pad-row index (%d) outside the range (%d -> %d) !",
\r
524 iRow, 0, nRows -1));
\r
528 //check pad consistency
\r
529 if (iPad < 0 || iPad >= (Int_t)(roc->GetNPads(iSec,iRow))) {
\r
530 AliError(Form("Pad index (%d) outside the range (%d -> %d) !",
\r
531 iPad, 0, roc->GetNPads(iSec,iRow) ));
\r
535 //loop over bunches
\r
536 while ( input.NextBunch() ){
\r
537 Int_t startTbin = (Int_t)input.GetStartTimeBin();
\r
538 Int_t bunchlength = (Int_t)input.GetBunchLength();
\r
539 const UShort_t *sig = input.GetSignals();
\r
541 for (Int_t iTime = 0; iTime<bunchlength; iTime++){
\r
542 Int_t iTimeBin=startTbin-iTime;
\r
544 if(fDebugLevel==72){
\r
545 fHistoRow->Fill(iRow);
\r
546 fHistoPad->Fill(iPad);
\r
547 fHistoTime->Fill(iTimeBin);
\r
548 fHistoRowPad->Fill(iPad,iRow);
\r
549 }else if(fDebugLevel>=0&&fDebugLevel<72){
\r
550 if(iSec==fDebugLevel){
\r
551 fHistoRow->Fill(iRow);
\r
552 fHistoPad->Fill(iPad);
\r
553 fHistoTime->Fill(iTimeBin);
\r
554 fHistoRowPad->Fill(iPad,iRow);
\r
556 }else if(fDebugLevel==73){
\r
558 fHistoRow->Fill(iRow);
\r
559 fHistoPad->Fill(iPad);
\r
560 fHistoTime->Fill(iTimeBin);
\r
561 fHistoRowPad->Fill(iPad,iRow);
\r
563 }else if(fDebugLevel==74){
\r
565 fHistoRow->Fill(iRow);
\r
566 fHistoPad->Fill(iPad);
\r
567 fHistoTime->Fill(iTimeBin);
\r
568 fHistoRowPad->Fill(iPad,iRow);
\r
572 //check time consistency
\r
573 if ( iTimeBin < fRecoParam->GetFirstBin() || iTimeBin >= fRecoParam->GetLastBin()){
\r
574 //cout<<iTimeBin<<endl;
\r
576 AliFatal(Form("Timebin index (%d) outside the range (%d -> %d) !",
\r
577 iTimeBin, 0, fRecoParam->GetLastBin() -1));
\r
580 Float_t signal=(Float_t)sig[iTime];
\r
581 if (signal <= fZeroSup ||
\r
582 iTimeBin < fFirstBin ||
\r
583 iTimeBin > fLastBin
\r
585 digarr->GetRow(iSec,iRow)->SetDigitFast(0,iTimeBin,iPad);
\r
588 if (!noiseROC) continue;
\r
589 Double_t noiseOnPad = noiseROC->GetValue(iRow,iPad);//noise on given pad and row in sector
\r
590 if (noiseOnPad > fMaxNoiseAbs){
\r
591 digarr->GetRow(iSec,iRow)->SetDigitFast(0,iTimeBin,iPad);
\r
592 continue; // consider noisy pad as dead
\r
594 if(signal <= fMaxNoiseSigma * noiseOnPad){
\r
595 digarr->GetRow(iSec,iRow)->SetDigitFast(0,iTimeBin,iPad);
\r
598 digarr->GetRow(iSec,iRow)->SetDigitFast(TMath::Nint(signal),iTimeBin,iPad);
\r
599 }// end loop signals in bunch
\r
600 }// end loop bunches
\r
603 }// end sector loop
\r
605 if(isAltro) FindClusterKrIO();
\r
615 Int_t AliTPCclustererKr::FinderIOold(AliRawReader* rawReader)
\r
617 // Krypton cluster finder for the TPC raw data
\r
619 // fParam must be defined before
\r
621 if(rawReader)fRawData=kTRUE; //set flag to data
\r
624 Error("Digits2Clusters", "output tree not initialised");
\r
628 fParam->SetMaxTBin(fRecoParam->GetLastBin());//set number of timebins from reco -> param
\r
629 // used later for memory allocation
\r
631 AliRawEventHeaderBase* eventHeader = (AliRawEventHeaderBase*)rawReader->GetEventHeader();
\r
633 fTimeStamp = eventHeader->Get("Timestamp");
\r
634 fRun = rawReader->GetRunNumber();
\r
637 Bool_t isAltro=kFALSE;
\r
639 AliTPCROC * roc = AliTPCROC::Instance();
\r
640 AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
\r
641 AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping();
\r
643 AliTPCRawStream input(rawReader,(AliAltroMapping**)mapping);
\r
645 const Int_t kNIS = fParam->GetNInnerSector();//number of inner sectors
\r
646 const Int_t kNOS = fParam->GetNOuterSector();//number of outer sectors
\r
647 const Int_t kNS = kNIS + kNOS;//all sectors
\r
651 AliTPCDigitsArray *digarr=new AliTPCDigitsArray(kFALSE);//data not sim
\r
652 digarr->Setup(fParam);//as usually parameters
\r
655 // Loop over sectors
\r
657 for(Int_t iSec = 0; iSec < kNS; iSec++) {
\r
658 AliTPCCalROC * noiseROC;
\r
660 noiseROC =new AliTPCCalROC(iSec);//noise=0
\r
663 noiseROC = noiseTPC->GetCalROC(iSec); // noise per given sector
\r
665 Int_t nRows = 0; //number of rows in sector
\r
666 Int_t nDDLs = 0; //number of DDLs
\r
667 Int_t indexDDL = 0; //DDL index
\r
669 nRows = fParam->GetNRowLow();
\r
671 indexDDL = iSec * 2;
\r
673 nRows = fParam->GetNRowUp();
\r
675 indexDDL = (iSec-kNIS) * 4 + kNIS * 2;
\r
679 // Load the raw data for corresponding DDLs
\r
681 rawReader->Reset();
\r
682 rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1);
\r
686 // Allocate memory for rows in sector (pads(depends on row) x timebins)
\r
687 for(Int_t iRow = 0; iRow < nRows; iRow++) {
\r
688 digarr->CreateRow(iSec,iRow);
\r
689 }//end loop over rows
\r
691 rawReader->Reset();
\r
692 rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1);
\r
695 // Begin loop over altro data
\r
697 while (input.Next()) {
\r
699 //check sector consistency
\r
700 if (input.GetSector() != iSec)
\r
701 AliFatal(Form("Sector index mismatch ! Expected (%d), but got (%d) !",iSec,input.GetSector()));
\r
703 Int_t iRow = input.GetRow();
\r
704 Int_t iPad = input.GetPad();
\r
705 Int_t iTimeBin = input.GetTime();
\r
708 if(fDebugLevel==72){
\r
709 fHistoRow->Fill(iRow);
\r
710 fHistoPad->Fill(iPad);
\r
711 fHistoTime->Fill(iTimeBin);
\r
712 fHistoRowPad->Fill(iPad,iRow);
\r
713 }else if(fDebugLevel>=0&&fDebugLevel<72){
\r
714 if(iSec==fDebugLevel){
\r
715 fHistoRow->Fill(iRow);
\r
716 fHistoPad->Fill(iPad);
\r
717 fHistoTime->Fill(iTimeBin);
\r
718 fHistoRowPad->Fill(iPad,iRow);
\r
720 }else if(fDebugLevel==73){
\r
722 fHistoRow->Fill(iRow);
\r
723 fHistoPad->Fill(iPad);
\r
724 fHistoTime->Fill(iTimeBin);
\r
725 fHistoRowPad->Fill(iPad,iRow);
\r
727 }else if(fDebugLevel==74){
\r
729 fHistoRow->Fill(iRow);
\r
730 fHistoPad->Fill(iPad);
\r
731 fHistoTime->Fill(iTimeBin);
\r
732 fHistoRowPad->Fill(iPad,iRow);
\r
736 //check row consistency
\r
737 if (iRow < 0 ) continue;
\r
738 if (iRow < 0 || iRow >= nRows){
\r
739 AliError(Form("Pad-row index (%d) outside the range (%d -> %d) !",
\r
740 iRow, 0, nRows -1));
\r
744 //check pad consistency
\r
745 if (iPad < 0 || iPad >= (Int_t)(roc->GetNPads(iSec,iRow))) {
\r
746 AliError(Form("Pad index (%d) outside the range (%d -> %d) !",
\r
747 iPad, 0, roc->GetNPads(iSec,iRow) ));
\r
751 //check time consistency
\r
752 if ( iTimeBin < fRecoParam->GetFirstBin() || iTimeBin >= fRecoParam->GetLastBin()){
\r
753 //cout<<iTimeBin<<endl;
\r
755 AliFatal(Form("Timebin index (%d) outside the range (%d -> %d) !",
\r
756 iTimeBin, 0, fRecoParam->GetLastBin() -1));
\r
760 Int_t signal = input.GetSignal();
\r
761 if (signal <= fZeroSup ||
\r
762 iTimeBin < fFirstBin ||
\r
763 iTimeBin > fLastBin
\r
765 digarr->GetRow(iSec,iRow)->SetDigitFast(0,iTimeBin,iPad);
\r
768 if (!noiseROC) continue;
\r
769 Double_t noiseOnPad = noiseROC->GetValue(iRow,iPad);//noise on given pad and row in sector
\r
770 if (noiseOnPad > fMaxNoiseAbs){
\r
771 digarr->GetRow(iSec,iRow)->SetDigitFast(0,iTimeBin,iPad);
\r
772 continue; // consider noisy pad as dead
\r
774 if(signal <= fMaxNoiseSigma * noiseOnPad){
\r
775 digarr->GetRow(iSec,iRow)->SetDigitFast(0,iTimeBin,iPad);
\r
778 digarr->GetRow(iSec,iRow)->SetDigitFast(signal,iTimeBin,iPad);
\r
779 }//end of loop over altro data
\r
780 }//end of loop over sectors
\r
783 if(isAltro) FindClusterKrIO();
\r
789 void AliTPCclustererKr::CleanSector(Int_t sector){
\r
791 // clean isolated digits
\r
793 const Int_t kNRows=fParam->GetNRow(sector);//number of rows in sector
\r
794 for(Int_t iRow=0; iRow<kNRows; ++iRow){
\r
795 AliSimDigits *digrow;
\r
797 digrow = (AliSimDigits*)fDigarr->GetRow(sector,iRow);//real data
\r
799 digrow = (AliSimDigits*)fDigarr->LoadRow(sector,iRow);//MC
\r
801 if(!digrow) continue;
\r
802 digrow->ExpandBuffer(); //decrunch
\r
803 const Int_t kNPads = digrow->GetNCols(); // number of pads
\r
804 const Int_t kNTime = digrow->GetNRows(); // number of timebins
\r
805 for(Int_t iPad=1;iPad<kNPads-1;iPad++){
\r
806 Short_t* val = digrow->GetDigitsColumn(iPad);
\r
808 for(Int_t iTimeBin=1;iTimeBin<kNTime-1;iTimeBin++){
\r
809 if (val[iTimeBin]<=0) continue;
\r
810 if (val[iTimeBin-1]+val[iTimeBin+1]<fIsolCut) {val[iTimeBin]=0; continue;}
\r
811 if (val[iTimeBin-kNTime]+val[iTimeBin+kNTime]<fIsolCut) {val[iTimeBin]=0; continue;}
\r
813 if (val[iTimeBin-1-kNTime]+val[iTimeBin+1+kNTime]<fIsolCut) {val[iTimeBin]=0; continue;}
\r
814 if (val[iTimeBin+1-kNTime]+val[iTimeBin-1+kNTime]<fIsolCut) {val[iTimeBin]=0; continue;}
\r
822 ////____________________________________________________________________________
\r
823 Int_t AliTPCclustererKr::FindClusterKrIO()
\r
827 //fParam and fDigarr must be set to run this method
\r
830 Int_t clusterCounter=0;
\r
831 const Int_t nTotalSector=fParam->GetNSector();//number of sectors
\r
832 for(Int_t iSec=0; iSec<nTotalSector; ++iSec){
\r
835 //vector of maxima for each sector
\r
836 //std::vector<AliPadMax*> maximaInSector;
\r
837 TObjArray *maximaInSector=new TObjArray();//to store AliPadMax*
\r
840 // looking for the maxima on the pad
\r
843 const Int_t kNRows=fParam->GetNRow(iSec);//number of rows in sector
\r
844 for(Int_t iRow=0; iRow<kNRows; ++iRow){
\r
845 AliSimDigits *digrow;
\r
847 digrow = (AliSimDigits*)fDigarr->GetRow(iSec,iRow);//real data
\r
849 digrow = (AliSimDigits*)fDigarr->LoadRow(iSec,iRow);//MC
\r
851 if(digrow){//if pointer exist
\r
852 digrow->ExpandBuffer(); //decrunch
\r
853 const Int_t kNPads = digrow->GetNCols(); // number of pads
\r
854 const Int_t kNTime = digrow->GetNRows(); // number of timebins
\r
855 for(Int_t iPad=0;iPad<kNPads;iPad++){
\r
857 Int_t timeBinMax=-1;//timebin of maximum
\r
858 Int_t valueMaximum=-1;//value of maximum in adc
\r
859 Int_t increaseBegin=-1;//timebin when increase starts
\r
860 Int_t sumAdc=0;//sum of adc on the pad in maximum surrounding
\r
861 bool ifIncreaseBegin=true;//flag - check if increasing started
\r
862 bool ifMaximum=false;//flag - check if it could be maximum
\r
863 Short_t* val = digrow->GetDigitsColumn(iPad);
\r
864 for(Int_t iTimeBin=1;iTimeBin<kNTime-1;iTimeBin++){
\r
866 if (val[iTimeBin]==-1) break; // 0 until the end
\r
867 for( ; iTimeBin<kNTime-2&&val[iTimeBin]<fMinAdc ;iTimeBin++) {}
\r
870 Short_t adc = val[iTimeBin];
\r
872 if(adc<fMinAdc){//standard was 3 for fMinAdc
\r
874 if(iTimeBin-increaseBegin<fMinTimeBins){//at least 2 time bins
\r
879 ifIncreaseBegin=true;
\r
883 //insert maximum, default values and set flags
\r
884 //Double_t xCord,yCord;
\r
885 //GetXY(iSec,iRow,iPad,xCord,yCord);
\r
886 Double_t x[]={iRow,iPad,iTimeBin};
\r
888 AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
\r
890 transform->Transform(x,i,0,1);
\r
892 AliPadMax *oneMaximum = new AliPadMax(AliTPCvtpr(valueMaximum,
\r
898 x[2]/*timeBinMax*/),
\r
902 maximaInSector->AddLast(oneMaximum);
\r
908 ifIncreaseBegin=true;
\r
919 if(ifIncreaseBegin){
\r
920 ifIncreaseBegin=false;
\r
921 increaseBegin=iTimeBin;
\r
924 if(adc>valueMaximum){
\r
925 timeBinMax=iTimeBin;
\r
930 if(iTimeBin==kNTime-1 && ifMaximum && kNTime-increaseBegin>fMinTimeBins){//on the edge
\r
931 //at least 3 timebins
\r
932 //insert maximum, default values and set flags
\r
933 //Double_t xCord,yCord;
\r
934 //GetXY(iSec,iRow,iPad,xCord,yCord);
\r
935 Double_t x[]={iRow,iPad,iTimeBin};
\r
937 //AliTPCTransform trafo;
\r
938 //trafo.Transform(x,i,0,1);
\r
940 AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
\r
942 transform->Transform(x,i,0,1);
\r
944 AliPadMax *oneMaximum = new AliPadMax(AliTPCvtpr(valueMaximum,
\r
950 x[2]/*timeBinMax*/),
\r
954 maximaInSector->AddLast(oneMaximum);
\r
960 ifIncreaseBegin=true;
\r
965 }//end loop over timebins
\r
966 }//end loop over pads
\r
968 // cout<<"Pointer does not exist!!"<<endl;
\r
969 }//end if poiner exists
\r
970 }//end loop over rows
\r
972 MakeClusters(maximaInSector,iSec,clusterCounter);
\r
974 maximaInSector->SetOwner(kTRUE);
\r
975 maximaInSector->Delete();
\r
976 delete maximaInSector;
\r
978 cout<<"Number of clusters in event: "<<clusterCounter<<endl;
\r
982 void AliTPCclustererKr::MakeClusters(TObjArray * maximaInSector, Int_t iSec, Int_t &clusterCounter){
\r
989 Int_t maxTimeBin=0;
\r
995 Int_t entriesArr = maximaInSector->GetEntriesFast();
\r
996 for(Int_t it1 = 0; it1 < entriesArr; ++it1 ) {
\r
998 AliPadMax *mp1=(AliPadMax *)maximaInSector->UncheckedAt(it1);
\r
999 if (!mp1) continue;
\r
1000 AliTPCclusterKr clusterKr;
\r
1002 Int_t nUsedPads=1;
\r
1003 Int_t clusterValue=0;
\r
1004 clusterValue+=(mp1)->GetSum();
\r
1005 list<Int_t> nUsedRows;
\r
1006 nUsedRows.push_back((mp1)->GetRow());
\r
1008 maxDig =(mp1)->GetAdc() ;
\r
1009 maxSumAdc =(mp1)->GetSum() ;
\r
1010 maxTimeBin =(mp1)->GetTime();
\r
1011 maxPad =(mp1)->GetPad() ;
\r
1012 maxRow =(mp1)->GetRow() ;
\r
1013 maxX =(mp1)->GetX();
\r
1014 maxY =(mp1)->GetY();
\r
1015 maxT =(mp1)->GetT();
\r
1017 AliSimDigits *digrowTmp;
\r
1019 digrowTmp = (AliSimDigits*)fDigarr->GetRow(iSec,(mp1)->GetRow());
\r
1021 digrowTmp = (AliSimDigits*)fDigarr->LoadRow(iSec,(mp1)->GetRow());
\r
1024 digrowTmp->ExpandBuffer(); //decrunch
\r
1026 for(Int_t itb=(mp1)->GetBegin(); itb<((mp1)->GetEnd())+1; itb++){
\r
1027 Int_t adcTmp = digrowTmp->GetDigitUnchecked(itb,(mp1)->GetPad());
\r
1028 AliTPCvtpr *vtpr=new AliTPCvtpr(adcTmp,itb,(mp1)->GetPad(),(mp1)->GetRow(),(mp1)->GetX(),(mp1)->GetY(),(mp1)->GetT());
\r
1029 clusterKr.AddDigitToCluster(vtpr);
\r
1031 clusterKr.SetCenter();//set centr of the cluster
\r
1033 for(Int_t it2 = it1+1; it2 < entriesArr; ++it2 ) {
\r
1034 AliPadMax *mp2=(AliPadMax *)maximaInSector->UncheckedAt(it2);
\r
1035 if (!mp2) continue;
\r
1036 if (TMath::Abs(clusterKr.GetCenterX() - (mp2)->GetX()) > fMaxPadRangeCm) continue;
\r
1037 if (TMath::Abs(clusterKr.GetCenterY() - (mp2)->GetY()) > fMaxRowRangeCm) continue;
\r
1038 if (TMath::Abs(clusterKr.GetCenterT() - (mp2)->GetT()) > fMaxTimeRange) continue;
\r
1041 clusterValue+=(mp2)->GetSum();
\r
1044 nUsedRows.push_back((mp2)->GetRow());
\r
1046 AliSimDigits *digrowTmp1;
\r
1048 digrowTmp1 = (AliSimDigits*)fDigarr->GetRow(iSec,(mp2)->GetRow());
\r
1050 digrowTmp1 = (AliSimDigits*)fDigarr->LoadRow(iSec,(mp2)->GetRow());
\r
1052 digrowTmp1->ExpandBuffer(); //decrunch
\r
1054 for(Int_t itb=(mp2)->GetBegin(); itb<(mp2)->GetEnd()+1; itb++){
\r
1055 Int_t adcTmp = digrowTmp1->GetDigitUnchecked(itb,(mp2)->GetPad());
\r
1056 AliTPCvtpr *vtpr=new AliTPCvtpr(adcTmp,itb,(mp2)->GetPad(),(mp2)->GetRow(),(mp2)->GetX(),(mp2)->GetY(),(mp2)->GetT());
\r
1057 clusterKr.AddDigitToCluster(vtpr);
\r
1060 clusterKr.SetCenter();//set center of the cluster
\r
1062 //which one is bigger
\r
1063 if( (mp2)->GetAdc() > maxDig ){
\r
1064 maxDig =(mp2)->GetAdc() ;
\r
1065 maxSumAdc =(mp2)->GetSum() ;
\r
1066 maxTimeBin =(mp2)->GetTime();
\r
1067 maxPad =(mp2)->GetPad() ;
\r
1068 maxRow =(mp2)->GetRow() ;
\r
1069 maxX =(mp2)->GetX() ;
\r
1070 maxY =(mp2)->GetY() ;
\r
1071 maxT =(mp2)->GetT() ;
\r
1072 } else if ( (mp2)->GetAdc() == maxDig ){
\r
1073 if( (mp2)->GetSum() > maxSumAdc){
\r
1074 maxDig =(mp2)->GetAdc() ;
\r
1075 maxSumAdc =(mp2)->GetSum() ;
\r
1076 maxTimeBin =(mp2)->GetTime();
\r
1077 maxPad =(mp2)->GetPad() ;
\r
1078 maxRow =(mp2)->GetRow() ;
\r
1079 maxX =(mp2)->GetX() ;
\r
1080 maxY =(mp2)->GetY() ;
\r
1081 maxT =(mp2)->GetT() ;
\r
1084 delete maximaInSector->RemoveAt(it2);
\r
1087 delete maximaInSector->RemoveAt(it1);
\r
1088 clusterKr.SetSize();
\r
1089 //through out clusters on the edge and noise
\r
1090 //if(clusterValue/clusterKr.fCluster.size()<fValueToSize)continue;
\r
1091 if(clusterValue/(clusterKr.GetSize())<fValueToSize)continue;
\r
1093 clusterKr.SetADCcluster(clusterValue);
\r
1094 clusterKr.SetNPads(nUsedPads);
\r
1095 clusterKr.SetMax(AliTPCvtpr(maxDig,maxTimeBin,maxPad,maxRow,maxX,maxY,maxT));
\r
1096 clusterKr.SetSec(iSec);
\r
1097 clusterKr.SetSize();
\r
1100 nUsedRows.unique();
\r
1101 clusterKr.SetNRows(nUsedRows.size());
\r
1102 clusterKr.SetCenter();
\r
1104 clusterKr.SetRMS();//Set pad,row,timebin RMS
\r
1105 clusterKr.Set1D();//Set size in pads and timebins
\r
1107 clusterKr.SetTimeStamp(fTimeStamp);
\r
1108 clusterKr.SetRun(fRun);
\r
1113 //save each cluster into file
\r
1115 (*fOutput)<<"Kr"<<
\r
1116 "Cl.="<<&clusterKr<<
\r
1119 //end of save each cluster into file adc.root
\r
1125 ////____________________________________________________________________________
\r
1128 void AliTPCclustererKr::GetXY(Int_t sec,Int_t row,Int_t pad,Double_t& xGlob,Double_t& yGlob){
\r
1130 //gives global XY coordinate of the pad
\r
1133 Double_t yLocal = fParam->GetPadRowRadii(sec,row);//radius of row in sector in cm
\r
1135 Int_t padmax = fParam->GetNPads(sec,row);//number of pads in a given row
\r
1137 if(sec<fParam->GetNInnerSector())padXSize=0.4;
\r
1138 else padXSize=0.6;
\r
1139 Double_t xLocal=(pad+0.5-padmax/2.)*padXSize;//x-value of the center of pad
\r
1142 fParam->AdjustCosSin((Int_t)sec,cos,sin);//return sinus and cosinus of the sector
\r
1144 Double_t xGlob1 = xLocal * cos + yLocal * sin;
\r
1145 Double_t yGlob1 = -xLocal * sin + yLocal * cos;
\r
1149 rot=TMath::Pi()/2.;
\r
1151 xGlob = xGlob1 * TMath::Cos(rot) + yGlob1 * TMath::Sin(rot);
\r
1152 yGlob = -xGlob1 * TMath::Sin(rot) + yGlob1 * TMath::Cos(rot);
\r
1155 if(sec<18||(sec>=36&&sec<54)) xGlob =-1*xGlob;
\r