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 cluster class
\r
21 // Origin: Adam Matyja, INP PAN, adam.matyja@ifj.edu.pl
\r
22 //-----------------------------------------------------------------
\r
24 #include "AliTPCclustererKr.h"
\r
25 #include "AliTPCclusterKr.h"
\r
28 #include "TObject.h"
\r
29 #include "AliPadMax.h"
\r
30 #include "AliSimDigits.h"
\r
31 #include "AliTPCv4.h"
\r
32 #include "AliTPCParam.h"
\r
33 #include "AliTPCDigitsArray.h"
\r
34 #include "AliTPCvtpr.h"
\r
35 #include "AliTPCClustersRow.h"
\r
38 //used in raw data finder
\r
39 #include "AliTPCROC.h"
\r
40 #include "AliTPCCalPad.h"
\r
41 #include "AliTPCAltroMapping.h"
\r
42 #include "AliTPCcalibDB.h"
\r
43 #include "AliTPCRawStream.h"
\r
44 #include "AliTPCRecoParam.h"
\r
45 #include "AliTPCReconstructor.h"
\r
46 #include "AliRawReader.h"
\r
47 #include "AliTPCCalROC.h"
\r
49 ClassImp(AliTPCclustererKr)
\r
52 AliTPCclustererKr::AliTPCclustererKr()
\r
61 fIsOldRCUFormat(kFALSE),
\r
73 fMaxPadRangeCm(2.5),
\r
77 // default constructor
\r
81 AliTPCclustererKr::AliTPCclustererKr(const AliTPCclustererKr ¶m)
\r
90 fIsOldRCUFormat(kFALSE),
\r
102 fMaxPadRangeCm(2.5),
\r
103 fMaxRowRangeCm(3.5)
\r
106 // copy constructor
\r
108 fParam = param.fParam;
\r
109 fRecoParam = param.fRecoParam;
\r
110 fIsOldRCUFormat = param.fIsOldRCUFormat;
\r
111 fRawData = param.fRawData;
\r
112 fRowCl = param.fRowCl ;
\r
113 fInput = param.fInput ;
\r
114 fOutput = param.fOutput;
\r
115 fDigarr = param.fDigarr;
\r
116 fZeroSup = param.fZeroSup ;
\r
117 fFirstBin = param.fFirstBin ;
\r
118 fLastBin = param.fLastBin ;
\r
119 fMaxNoiseAbs = param.fMaxNoiseAbs ;
\r
120 fMaxNoiseSigma = param.fMaxNoiseSigma ;
\r
121 fMinAdc = param.fMinAdc;
\r
122 fMinTimeBins = param.fMinTimeBins;
\r
123 // fMaxPadRange = param.fMaxPadRange ;
\r
124 // fMaxRowRange = param.fMaxRowRange ;
\r
125 fMaxTimeRange = param.fMaxTimeRange;
\r
126 fValueToSize = param.fValueToSize;
\r
127 fMaxPadRangeCm = param.fMaxPadRangeCm;
\r
128 fMaxRowRangeCm = param.fMaxRowRangeCm;
\r
131 AliTPCclustererKr & AliTPCclustererKr::operator = (const AliTPCclustererKr & param)
\r
133 fParam = param.fParam;
\r
134 fRecoParam = param.fRecoParam;
\r
135 fIsOldRCUFormat = param.fIsOldRCUFormat;
\r
136 fRawData = param.fRawData;
\r
137 fRowCl = param.fRowCl ;
\r
138 fInput = param.fInput ;
\r
139 fOutput = param.fOutput;
\r
140 fDigarr = param.fDigarr;
\r
141 fZeroSup = param.fZeroSup ;
\r
142 fFirstBin = param.fFirstBin ;
\r
143 fLastBin = param.fLastBin ;
\r
144 fMaxNoiseAbs = param.fMaxNoiseAbs ;
\r
145 fMaxNoiseSigma = param.fMaxNoiseSigma ;
\r
146 fMinAdc = param.fMinAdc;
\r
147 fMinTimeBins = param.fMinTimeBins;
\r
148 // fMaxPadRange = param.fMaxPadRange ;
\r
149 // fMaxRowRange = param.fMaxRowRange ;
\r
150 fMaxTimeRange = param.fMaxTimeRange;
\r
151 fValueToSize = param.fValueToSize;
\r
152 fMaxPadRangeCm = param.fMaxPadRangeCm;
\r
153 fMaxRowRangeCm = param.fMaxRowRangeCm;
\r
157 AliTPCclustererKr::~AliTPCclustererKr()
\r
164 void AliTPCclustererKr::SetRecoParam(AliTPCRecoParam *recoParam)
\r
167 fRecoParam = recoParam;
\r
169 //set default parameters if not specified
\r
170 fRecoParam = AliTPCReconstructor::GetRecoParam();
\r
171 if (!fRecoParam) fRecoParam = AliTPCRecoParam::GetLowFluxParam();
\r
177 ////____________________________________________________________________________
\r
179 void AliTPCclustererKr::SetInput(TTree * tree)
\r
182 // set input tree with digits
\r
185 if (!fInput->GetBranch("Segment")){
\r
186 cerr<<"AliTPCclusterKr::FindClusterKr(): no proper input tree !\n";
\r
192 void AliTPCclustererKr::SetOutput(TTree * tree)
\r
198 AliTPCClustersRow clrow;
\r
199 AliTPCClustersRow *pclrow=&clrow;
\r
200 clrow.SetClass("AliTPCclusterKr");
\r
201 clrow.SetArray(1); // to make Clones array
\r
202 fOutput->Branch("Segment","AliTPCClustersRow",&pclrow,32000,200);
\r
205 ////____________________________________________________________________________
\r
207 Int_t AliTPCclustererKr::FinderIO()
\r
209 // Krypton cluster finder for simulated events from MC
\r
212 Error("Digits2Clusters", "input tree not initialised");
\r
217 Error("Digits2Clusters", "output tree not initialised");
\r
225 Int_t AliTPCclustererKr::FinderIO(AliRawReader* rawReader)
\r
227 // Krypton cluster finder for the TPC raw data
\r
229 // fParam must be defined before
\r
231 // consider noiceROC or not
\r
233 if(rawReader)fRawData=kTRUE; //set flag to data
\r
236 Error("Digits2Clusters", "output tree not initialised");
\r
240 fParam->SetMaxTBin(fRecoParam->GetLastBin());//set number of timebins from reco -> param
\r
241 // used later for memory allocation
\r
243 Bool_t isAltro=kFALSE;
\r
245 AliTPCROC * roc = AliTPCROC::Instance();
\r
246 AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
\r
247 AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping();
\r
249 AliTPCRawStream input(rawReader,(AliAltroMapping**)mapping);
\r
251 const Int_t kNIS = fParam->GetNInnerSector();//number of inner sectors
\r
252 const Int_t kNOS = fParam->GetNOuterSector();//number of outer sectors
\r
253 const Int_t kNS = kNIS + kNOS;//all sectors
\r
255 // //set cluster finder parameters
\r
256 // SetZeroSup(fParam->GetZeroSup());//zero suppression parameter
\r
257 // SetFirstBin(60);
\r
258 // SetLastBin(950);
\r
259 // SetMaxNoiseAbs(2);
\r
260 // SetMaxNoiseSigma(3);
\r
263 AliTPCDigitsArray *digarr=new AliTPCDigitsArray(kFALSE);//data not sim
\r
264 digarr->Setup(fParam);//as usually parameters
\r
267 // Loop over sectors
\r
269 for(Int_t iSec = 0; iSec < kNS; iSec++) {
\r
270 AliTPCCalROC * noiseROC;
\r
272 noiseROC =new AliTPCCalROC(iSec);//noise=0
\r
275 noiseROC = noiseTPC->GetCalROC(iSec); // noise per given sector
\r
277 Int_t nRows = 0; //number of rows in sector
\r
278 Int_t nDDLs = 0; //number of DDLs
\r
279 Int_t indexDDL = 0; //DDL index
\r
281 nRows = fParam->GetNRowLow();
\r
283 indexDDL = iSec * 2;
\r
285 nRows = fParam->GetNRowUp();
\r
287 indexDDL = (iSec-kNIS) * 4 + kNIS * 2;
\r
291 // Load the raw data for corresponding DDLs
\r
293 rawReader->Reset();
\r
294 input.SetOldRCUFormat(fIsOldRCUFormat);
\r
295 rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1);
\r
299 // Allocate memory for rows in sector (pads(depends on row) x timebins)
\r
300 for(Int_t iRow = 0; iRow < nRows; iRow++) {
\r
301 digarr->CreateRow(iSec,iRow);
\r
302 }//end loop over rows
\r
304 rawReader->Reset();
\r
305 input.SetOldRCUFormat(fIsOldRCUFormat);
\r
306 rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1);
\r
309 // Begin loop over altro data
\r
311 while (input.Next()) {
\r
313 //check sector consistency
\r
314 if (input.GetSector() != iSec)
\r
315 AliFatal(Form("Sector index mismatch ! Expected (%d), but got (%d) !",iSec,input.GetSector()));
\r
317 //check row consistency
\r
318 Short_t iRow = input.GetRow();
\r
319 if (iRow < 0 || iRow >= nRows){
\r
320 AliError(Form("Pad-row index (%d) outside the range (%d -> %d) !",
\r
321 iRow, 0, nRows -1));
\r
325 //check pad consistency
\r
326 Short_t iPad = input.GetPad();
\r
327 if (iPad < 0 || iPad >= (Short_t)(roc->GetNPads(iSec,iRow))) {
\r
328 AliError(Form("Pad index (%d) outside the range (%d -> %d) !",
\r
329 iPad, 0, roc->GetNPads(iSec,iRow) ));
\r
333 //check time consistency
\r
334 Short_t iTimeBin = input.GetTime();
\r
335 if ( iTimeBin < fRecoParam->GetFirstBin() || iTimeBin >= fRecoParam->GetLastBin()){
\r
336 //cout<<iTimeBin<<endl;
\r
338 AliFatal(Form("Timebin index (%d) outside the range (%d -> %d) !",
\r
339 iTimeBin, 0, fRecoParam->GetLastBin() -1));
\r
343 Int_t signal = input.GetSignal();
\r
344 if (signal <= fZeroSup ||
\r
345 iTimeBin < fFirstBin ||
\r
346 iTimeBin > fLastBin
\r
348 digarr->GetRow(iSec,iRow)->SetDigitFast(0,iTimeBin,iPad);
\r
351 if (!noiseROC) continue;
\r
352 Double_t noiseOnPad = noiseROC->GetValue(iRow,iPad);//noise on given pad and row in sector
\r
353 if (noiseOnPad > fMaxNoiseAbs) continue; // consider noisy pad as dead
\r
355 if(signal <= fMaxNoiseSigma * noiseOnPad){
\r
356 digarr->GetRow(iSec,iRow)->SetDigitFast(0,iTimeBin,iPad);
\r
360 digarr->GetRow(iSec,iRow)->SetDigitFast(signal,iTimeBin,iPad);
\r
362 }//end of loop over altro data
\r
363 }//end of loop over sectors
\r
366 if(isAltro) FindClusterKrIO();
\r
372 ////____________________________________________________________________________
\r
373 Int_t AliTPCclustererKr::FindClusterKrIO()
\r
375 //fParam and fDigarr must be set to run this method
\r
377 // //set parameters
\r
378 // SetMinAdc(3);//usually is 3
\r
379 // SetMinTimeBins(2);//should be 2 - the best result of shape in MC
\r
380 //// SetMaxPadRange(4);
\r
381 //// SetMaxRowRange(3);
\r
382 // SetMaxTimeRange(7);
\r
383 // SetValueToSize(3.5);//3.5
\r
384 // SetMaxPadRangeCm(2.5);
\r
385 // SetMaxRowRangeCm(3.5);
\r
387 Int_t clusterCounter=0;
\r
388 const Short_t nTotalSector=fParam->GetNSector();//number of sectors
\r
389 for(Short_t iSec=0; iSec<nTotalSector; ++iSec){
\r
391 //vector of maxima for each sector
\r
392 std::vector<AliPadMax*> maximaInSector;
\r
395 // looking for the maxima on the pad
\r
398 const Short_t nRows=fParam->GetNRow(iSec);//number of rows in sector
\r
399 for(Short_t iRow=0; iRow<nRows; ++iRow){
\r
400 AliSimDigits *digrow;
\r
402 digrow = (AliSimDigits*)fDigarr->GetRow(iSec,iRow);//real data
\r
404 digrow = (AliSimDigits*)fDigarr->LoadRow(iSec,iRow);//MC
\r
406 if(digrow){//if pointer exist
\r
407 digrow->ExpandBuffer(); //decrunch
\r
408 const Short_t nPads = digrow->GetNCols(); // number of pads
\r
409 const Short_t nTime = digrow->GetNRows(); // number of timebins
\r
410 for(Short_t iPad=0;iPad<nPads;iPad++){
\r
412 Short_t timeBinMax=-1;//timebin of maximum
\r
413 Short_t valueMaximum=-1;//value of maximum in adc
\r
414 Short_t increaseBegin=-1;//timebin when increase starts
\r
415 Short_t sumAdc=0;//sum of adc on the pad in maximum surrounding
\r
416 bool ifIncreaseBegin=true;//flag - check if increasing started
\r
417 bool ifMaximum=false;//flag - check if it could be maximum
\r
419 for(Short_t iTimeBin=0;iTimeBin<nTime;iTimeBin++){
\r
420 Short_t adc = digrow->GetDigitFast(iTimeBin,iPad);
\r
421 if(adc<fMinAdc){//standard was 3
\r
423 if(iTimeBin-increaseBegin<fMinTimeBins){//at least 2 time bins
\r
428 ifIncreaseBegin=true;
\r
432 //insert maximum, default values and set flags
\r
433 Double_t xCord,yCord;
\r
434 GetXY(iSec,iRow,iPad,xCord,yCord);
\r
435 AliPadMax *oneMaximum = new AliPadMax(AliTPCvtpr(valueMaximum,
\r
445 maximaInSector.push_back(oneMaximum);
\r
451 ifIncreaseBegin=true;
\r
457 if(ifIncreaseBegin){
\r
458 ifIncreaseBegin=false;
\r
459 increaseBegin=iTimeBin;
\r
462 if(adc>valueMaximum){
\r
463 timeBinMax=iTimeBin;
\r
468 if(iTimeBin==nTime-1 && ifMaximum && nTime-increaseBegin>fMinTimeBins){//on the edge
\r
469 //at least 3 timebins
\r
470 //insert maximum, default values and set flags
\r
471 Double_t xCord,yCord;
\r
472 GetXY(iSec,iRow,iPad,xCord,yCord);
\r
473 AliPadMax *oneMaximum = new AliPadMax(AliTPCvtpr(valueMaximum,
\r
483 maximaInSector.push_back(oneMaximum);
\r
489 ifIncreaseBegin=true;
\r
494 }//end loop over timebins
\r
495 }//end loop over pads
\r
497 // cout<<"Pointer does not exist!!"<<endl;
\r
498 }//end if poiner exists
\r
499 }//end loop over rows
\r
506 Short_t maxSumAdc=0;
\r
507 Short_t maxTimeBin=0;
\r
513 for( std::vector<AliPadMax*>::iterator mp1 = maximaInSector.begin();
\r
514 mp1 != maximaInSector.end(); ++mp1 ) {
\r
516 AliTPCclusterKr *tmp=new AliTPCclusterKr();
\r
518 Short_t nUsedPads=1;
\r
519 Int_t clusterValue=0;
\r
520 clusterValue+=(*mp1)->GetSum();
\r
521 list<Short_t> nUsedRows;
\r
522 nUsedRows.push_back((*mp1)->GetRow());
\r
524 maxDig =(*mp1)->GetAdc() ;
\r
525 maxSumAdc =(*mp1)->GetSum() ;
\r
526 maxTimeBin =(*mp1)->GetTime();
\r
527 maxPad =(*mp1)->GetPad() ;
\r
528 maxRow =(*mp1)->GetRow() ;
\r
529 maxX =(*mp1)->GetX();
\r
530 maxY =(*mp1)->GetY();
\r
532 AliSimDigits *digrowTmp;
\r
534 digrowTmp = (AliSimDigits*)fDigarr->GetRow(iSec,(*mp1)->GetRow());
\r
536 digrowTmp = (AliSimDigits*)fDigarr->LoadRow(iSec,(*mp1)->GetRow());
\r
539 digrowTmp->ExpandBuffer(); //decrunch
\r
541 for(Short_t itb=(*mp1)->GetBegin(); itb<((*mp1)->GetEnd())+1; itb++){
\r
542 Short_t adcTmp = digrowTmp->GetDigitFast(itb,(*mp1)->GetPad());
\r
543 AliTPCvtpr *vtpr=new AliTPCvtpr(adcTmp,itb,(*mp1)->GetPad(),(*mp1)->GetRow(),(*mp1)->GetX(),(*mp1)->GetY(),itb);
\r
544 tmp->fCluster.push_back(vtpr);
\r
547 tmp->SetCenter();//set centr of the cluster
\r
549 maximaInSector.erase(mp1);
\r
552 for( std::vector<AliPadMax*>::iterator mp2 = maximaInSector.begin();
\r
553 mp2 != maximaInSector.end(); ++mp2 ) {
\r
554 // if(abs(maxRow - (*mp2)->GetRow()) < fMaxRowRange && //3
\r
555 // abs(maxPad - (*mp2)->GetPad()) < fMaxPadRange && //4
\r
557 if(TMath::Abs(tmp->GetCenterX() - (*mp2)->GetX()) < fMaxPadRangeCm &&
\r
558 TMath::Abs(tmp->GetCenterY() - (*mp2)->GetY()) < fMaxRowRangeCm &&
\r
559 TMath::Abs(tmp->GetCenterT() - (*mp2)->GetT()) < fMaxTimeRange){//7
\r
561 clusterValue+=(*mp2)->GetSum();
\r
564 nUsedRows.push_back((*mp2)->GetRow());
\r
566 AliSimDigits *digrowTmp1;
\r
568 digrowTmp1 = (AliSimDigits*)fDigarr->GetRow(iSec,(*mp2)->GetRow());
\r
570 digrowTmp1 = (AliSimDigits*)fDigarr->LoadRow(iSec,(*mp2)->GetRow());
\r
572 digrowTmp1->ExpandBuffer(); //decrunch
\r
574 for(Short_t itb=(*mp2)->GetBegin(); itb<(*mp2)->GetEnd()+1; itb++){
\r
575 Short_t adcTmp = digrowTmp1->GetDigitFast(itb,(*mp2)->GetPad());
\r
576 AliTPCvtpr *vtpr=new AliTPCvtpr(adcTmp,itb,(*mp2)->GetPad(),(*mp2)->GetRow(),(*mp2)->GetX(),(*mp2)->GetY(),itb);
\r
577 tmp->fCluster.push_back(vtpr);
\r
580 tmp->SetCenter();//set centr of the cluster
\r
582 //which one is bigger
\r
583 if( (*mp2)->GetAdc() > maxDig ){
\r
584 maxDig =(*mp2)->GetAdc() ;
\r
585 maxSumAdc =(*mp2)->GetSum() ;
\r
586 maxTimeBin =(*mp2)->GetTime();
\r
587 maxPad =(*mp2)->GetPad() ;
\r
588 maxRow =(*mp2)->GetRow() ;
\r
589 maxX =(*mp2)->GetX() ;
\r
590 maxY =(*mp2)->GetY() ;
\r
591 } else if ( (*mp2)->GetAdc() == maxDig ){
\r
592 if( (*mp2)->GetSum() > maxSumAdc){
\r
593 maxDig =(*mp2)->GetAdc() ;
\r
594 maxSumAdc =(*mp2)->GetSum() ;
\r
595 maxTimeBin =(*mp2)->GetTime();
\r
596 maxPad =(*mp2)->GetPad() ;
\r
597 maxRow =(*mp2)->GetRow() ;
\r
598 maxX =(*mp2)->GetX() ;
\r
599 maxY =(*mp2)->GetY() ;
\r
602 maximaInSector.erase(mp2);
\r
607 //through out ADC=6,7 on 1 pad, 2 tb and ADC=12 on 2 pads,2 tb
\r
608 //if(nUsedPads==1 && clusterValue/tmp->fCluster.size()<3.6)continue;
\r
609 //if(nUsedPads==2 && clusterValue/tmp->fCluster.size()<3.1)continue;
\r
611 //through out clusters on the edge and noise
\r
612 if(clusterValue/tmp->fCluster.size()<fValueToSize)continue;
\r
614 tmp->SetADCcluster(clusterValue);
\r
615 tmp->SetNPads(nUsedPads);
\r
616 tmp->SetMax(AliTPCvtpr(maxDig,maxTimeBin,maxPad,maxRow,maxX,maxY,maxTimeBin));
\r
621 nUsedRows.unique();
\r
622 tmp->SetNRows(nUsedRows.size());
\r
627 //save each cluster into file
\r
629 AliTPCClustersRow *clrow= new AliTPCClustersRow();
\r
631 clrow->SetClass("AliTPCclusterKr");
\r
632 clrow->SetArray(1);
\r
633 fOutput->GetBranch("Segment")->SetAddress(&clrow);
\r
635 Int_t tmpCluster=0;
\r
636 TClonesArray * arr = fRowCl->GetArray();
\r
637 AliTPCclusterKr * cl = new ((*arr)[tmpCluster]) AliTPCclusterKr(*tmp);
\r
642 //end of save each cluster into file adc.root
\r
646 cout<<"Number of clusters in event: "<<clusterCounter<<endl;
\r
651 ////____________________________________________________________________________
\r
654 void AliTPCclustererKr::GetXY(Short_t sec,Short_t row,Short_t pad,Double_t& xGlob,Double_t& yGlob){
\r
656 Double_t yLocal = fParam->GetPadRowRadii(sec,row);//radius of row in sector in cm
\r
658 Short_t padmax = fParam->GetNPads(sec,row);//number of pads in a given row
\r
660 if(sec<fParam->GetNInnerSector())padXSize=0.4;
\r
662 Double_t xLocal=(pad+0.5-padmax)*padXSize;//x-value of the center of pad
\r
665 fParam->AdjustCosSin((Int_t)sec,cos,sin);//return sinus and cosinus of the sector
\r
667 xGlob = xLocal * cos + yLocal * sin;
\r
668 yGlob = -xLocal * sin + yLocal * cos;
\r