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
37 To run clusterizaton for MC type:
\r
40 If you don't want to use the standard selection criteria then you
\r
41 have to do following:
\r
43 // load the standard setup
\r
44 AliRunLoader* rl = AliRunLoader::Open("galice.root");
\r
45 AliTPCLoader *tpcl = (AliTPCLoader*)rl->GetLoader("TPCLoader");
\r
48 gAlice=rl->GetAliRun();
\r
49 TDirectory *cwd = gDirectory;
\r
50 AliTPCv4 *tpc = (AliTPCv4*)gAlice->GetDetector("TPC");
\r
51 Int_t ver = tpc->IsVersion();
\r
53 AliTPCParam *param=(AliTPCParamSR *)gDirectory->Get("75x40_100x60_150x60");
\r
54 AliTPCDigitsArray *digarr=new AliTPCDigitsArray;
\r
55 digarr->Setup(param);
\r
59 Int_t nevmax=rl->GetNumberOfEvents();
\r
60 for(Int_t nev=0;nev<nevmax ;nev++){
\r
62 TTree* input_tree= tpcl->TreeD();//load tree with digits
\r
63 digarr->ConnectTree(input_tree);
\r
64 TTree *output_tree =tpcl->TreeR();//load output tree
\r
66 AliTPCclustererKr *clusters = new AliTPCclustererKr();
\r
67 clusters->SetParam(param);
\r
68 clusters->SetInput(input_tree);
\r
69 clusters->SetOutput(output_tree);
\r
70 clusters->SetDigArr(digarr);
\r
72 //If you want to change the cluster finder parameters for MC there are
\r
75 //1. signal threshold (everything below the given number is treated as 0)
\r
76 clusters->SetMinAdc(3);
\r
78 //2. number of neighbouring timebins to be considered
\r
79 clusters->SetMinTimeBins(2);
\r
81 //3. distance of the cluster center to the center of a pad in pad-padrow plane
\r
82 //(in cm). Remenber that this is still quantified by pad size.
\r
83 clusters->SetMaxPadRangeCm(2.5);
\r
85 //4. distance of the cluster center to the center of a padrow in pad-padrow
\r
86 //plane (in cm). Remenber that this is still quantified by pad size.
\r
87 clusters->SetMaxRowRangeCm(3.5);
\r
89 //5. distance of the cluster center to the max time bin on a pad (in tackts)
\r
90 //ie. fabs(centerT - time)<7
\r
91 clusters->SetMaxTimeRange(7);
\r
93 //6. cut reduce peak at 0. There are noises which appear mostly as two
\r
94 //timebins on one pad.
\r
95 clusters->SetValueToSize(3.5);
\r
98 clusters->finderIO();
\r
99 tpcl->WriteRecPoints("OVERWRITE");
\r
101 delete rl;//cleans everything
\r
104 ********* DATA *********
\r
107 To run clusterizaton for DATA for file named raw_data.root type:
\r
108 .x FindKrClustersRaw.C("raw_data.root")
\r
110 If you want to change some criteria do the following:
\r
113 // remove Altro warnings
\r
115 AliLog::SetClassDebugLevel("AliTPCRawStream",-5);
\r
116 AliLog::SetClassDebugLevel("AliRawReaderDate",-5);
\r
117 AliLog::SetClassDebugLevel("AliTPCAltroMapping",-5);
\r
118 AliLog::SetModuleDebugLevel("RAW",-5);
\r
121 // Get database with noises
\r
123 // char *ocdbpath = gSystem->Getenv("OCDB_PATH");
\r
124 char *ocdbpath ="local:///afs/cern.ch/alice/tpctest/OCDB";
\r
126 ocdbpath="alien://folder=/alice/data/2007/LHC07w/OCDB/";
\r
128 AliCDBManager * man = AliCDBManager::Instance();
\r
129 man->SetDefaultStorage(ocdbpath);
\r
131 AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
\r
132 AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping();
\r
135 TFile *hfile=new TFile("adc.root","RECREATE","ADC file");
\r
136 // Create a ROOT Tree
\r
137 TTree *mytree = new TTree("Kr","Krypton cluster tree");
\r
139 //define infput file
\r
140 const char *fileName="data.root";
\r
141 AliRawReader *reader = new AliRawReaderRoot(fileName);
\r
142 //AliRawReader *reader = new AliRawReaderDate(fileName);
\r
144 AliAltroRawStreamFast* stream = new AliAltroRawStreamFast(reader);
\r
145 stream->SelectRawData("TPC");
\r
147 //one general output
\r
148 AliTPCclustererKr *clusters = new AliTPCclustererKr();
\r
149 clusters->SetOutput(mytree);
\r
150 clusters->SetRecoParam(0);//standard reco parameters
\r
151 AliTPCParamSR *param=new AliTPCParamSR();
\r
152 clusters->SetParam(param);//TPC parameters(sectors, timebins, etc.)
\r
154 //set cluster finder parameters (from data):
\r
155 //1. zero suppression parameter
\r
156 clusters->SetZeroSup(param->GetZeroSup());
\r
159 clusters->SetFirstBin(60);
\r
162 clusters->SetLastBin(950);
\r
165 clusters->SetMaxNoiseAbs(2);
\r
167 //5. maximal amount of sigma of noise
\r
168 clusters->SetMaxNoiseSigma(3);
\r
170 //The remaining parameters are the same paramters as for MC (see MC section
\r
172 clusters->SetMinAdc(3);
\r
173 clusters->SetMinTimeBins(2);
\r
174 clusters->SetMaxPadRangeCm(2.5);
\r
175 clusters->SetMaxRowRangeCm(3.5);
\r
176 clusters->SetMaxTimeRange(7);
\r
177 clusters->SetValueToSize(3.5);
\r
179 while (reader->NextEvent()) {
\r
180 clusters->FinderIO(reader);
\r
190 #include "AliTPCclustererKr.h"
\r
191 #include "AliTPCclusterKr.h"
\r
192 //#include <vector>
\r
194 #include "TObject.h"
\r
195 #include "AliPadMax.h"
\r
196 #include "AliSimDigits.h"
\r
197 #include "AliTPCv4.h"
\r
198 #include "AliTPCParam.h"
\r
199 #include "AliTPCDigitsArray.h"
\r
200 #include "AliTPCvtpr.h"
\r
201 #include "AliTPCClustersRow.h"
\r
206 //used in raw data finder
\r
207 #include "AliTPCROC.h"
\r
208 #include "AliTPCCalPad.h"
\r
209 #include "AliTPCAltroMapping.h"
\r
210 #include "AliTPCcalibDB.h"
\r
211 #include "AliTPCRawStream.h"
\r
212 #include "AliTPCRecoParam.h"
\r
213 #include "AliTPCReconstructor.h"
\r
214 #include "AliRawReader.h"
\r
215 #include "AliTPCCalROC.h"
\r
217 ClassImp(AliTPCclustererKr)
\r
220 AliTPCclustererKr::AliTPCclustererKr()
\r
236 // fMaxPadRange(4),
\r
237 // fMaxRowRange(3),
\r
240 fMaxPadRangeCm(2.5),
\r
241 fMaxRowRangeCm(3.5),
\r
249 // default constructor
\r
253 AliTPCclustererKr::AliTPCclustererKr(const AliTPCclustererKr ¶m)
\r
269 // fMaxPadRange(4),
\r
270 // fMaxRowRange(3),
\r
273 fMaxPadRangeCm(2.5),
\r
274 fMaxRowRangeCm(3.5),
\r
282 // copy constructor
\r
284 fParam = param.fParam;
\r
285 fRecoParam = param.fRecoParam;
\r
286 fRawData = param.fRawData;
\r
287 fRowCl = param.fRowCl ;
\r
288 fInput = param.fInput ;
\r
289 fOutput = param.fOutput;
\r
290 fDigarr = param.fDigarr;
\r
291 fZeroSup = param.fZeroSup ;
\r
292 fFirstBin = param.fFirstBin ;
\r
293 fLastBin = param.fLastBin ;
\r
294 fMaxNoiseAbs = param.fMaxNoiseAbs ;
\r
295 fMaxNoiseSigma = param.fMaxNoiseSigma ;
\r
296 fMinAdc = param.fMinAdc;
\r
297 fMinTimeBins = param.fMinTimeBins;
\r
298 // fMaxPadRange = param.fMaxPadRange ;
\r
299 // fMaxRowRange = param.fMaxRowRange ;
\r
300 fMaxTimeRange = param.fMaxTimeRange;
\r
301 fValueToSize = param.fValueToSize;
\r
302 fMaxPadRangeCm = param.fMaxPadRangeCm;
\r
303 fMaxRowRangeCm = param.fMaxRowRangeCm;
\r
304 fDebugLevel = param.fDebugLevel;
\r
305 fHistoRow = param.fHistoRow ;
\r
306 fHistoPad = param.fHistoPad ;
\r
307 fHistoTime = param.fHistoTime;
\r
308 fHistoRowPad = param.fHistoRowPad;
\r
312 AliTPCclustererKr & AliTPCclustererKr::operator = (const AliTPCclustererKr & param)
\r
315 // assignment operator
\r
317 fParam = param.fParam;
\r
318 fRecoParam = param.fRecoParam;
\r
319 fRawData = param.fRawData;
\r
320 fRowCl = param.fRowCl ;
\r
321 fInput = param.fInput ;
\r
322 fOutput = param.fOutput;
\r
323 fDigarr = param.fDigarr;
\r
324 fZeroSup = param.fZeroSup ;
\r
325 fFirstBin = param.fFirstBin ;
\r
326 fLastBin = param.fLastBin ;
\r
327 fMaxNoiseAbs = param.fMaxNoiseAbs ;
\r
328 fMaxNoiseSigma = param.fMaxNoiseSigma ;
\r
329 fMinAdc = param.fMinAdc;
\r
330 fMinTimeBins = param.fMinTimeBins;
\r
331 // fMaxPadRange = param.fMaxPadRange ;
\r
332 // fMaxRowRange = param.fMaxRowRange ;
\r
333 fMaxTimeRange = param.fMaxTimeRange;
\r
334 fValueToSize = param.fValueToSize;
\r
335 fMaxPadRangeCm = param.fMaxPadRangeCm;
\r
336 fMaxRowRangeCm = param.fMaxRowRangeCm;
\r
337 fDebugLevel = param.fDebugLevel;
\r
338 fHistoRow = param.fHistoRow ;
\r
339 fHistoPad = param.fHistoPad ;
\r
340 fHistoTime = param.fHistoTime;
\r
341 fHistoRowPad = param.fHistoRowPad;
\r
345 AliTPCclustererKr::~AliTPCclustererKr()
\r
352 void AliTPCclustererKr::SetRecoParam(AliTPCRecoParam *recoParam)
\r
355 // set reconstruction parameters
\r
358 fRecoParam = recoParam;
\r
360 //set default parameters if not specified
\r
361 fRecoParam = AliTPCReconstructor::GetRecoParam();
\r
362 if (!fRecoParam) fRecoParam = AliTPCRecoParam::GetLowFluxParam();
\r
368 ////____________________________________________________________________________
\r
370 void AliTPCclustererKr::SetInput(TTree * tree)
\r
373 // set input tree with digits
\r
376 if (!fInput->GetBranch("Segment")){
\r
377 cerr<<"AliTPCclusterKr::FindClusterKr(): no proper input tree !\n";
\r
383 void AliTPCclustererKr::SetOutput(TTree * tree)
\r
389 AliTPCClustersRow clrow;
\r
390 AliTPCClustersRow *pclrow=&clrow;
\r
391 clrow.SetClass("AliTPCclusterKr");
\r
392 clrow.SetArray(1); // to make Clones array
\r
393 fOutput->Branch("Segment","AliTPCClustersRow",&pclrow,32000,200);
\r
396 ////____________________________________________________________________________
\r
398 Int_t AliTPCclustererKr::FinderIO()
\r
400 // Krypton cluster finder for simulated events from MC
\r
403 Error("Digits2Clusters", "input tree not initialised");
\r
408 Error("Digits2Clusters", "output tree not initialised");
\r
416 Int_t AliTPCclustererKr::FinderIO(AliRawReader* rawReader)
\r
418 // Krypton cluster finder for the TPC raw data
\r
420 // fParam must be defined before
\r
422 if(rawReader)fRawData=kTRUE; //set flag to data
\r
425 Error("Digits2Clusters", "output tree not initialised");
\r
429 fParam->SetMaxTBin(fRecoParam->GetLastBin());//set number of timebins from reco -> param
\r
430 // used later for memory allocation
\r
432 Bool_t isAltro=kFALSE;
\r
434 AliTPCROC * roc = AliTPCROC::Instance();
\r
435 AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
\r
436 AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping();
\r
438 AliTPCRawStream input(rawReader,(AliAltroMapping**)mapping);
\r
440 const Int_t kNIS = fParam->GetNInnerSector();//number of inner sectors
\r
441 const Int_t kNOS = fParam->GetNOuterSector();//number of outer sectors
\r
442 const Int_t kNS = kNIS + kNOS;//all sectors
\r
444 // //set cluster finder parameters
\r
445 // SetZeroSup(fParam->GetZeroSup());//zero suppression parameter
\r
446 // SetFirstBin(60);
\r
447 // SetLastBin(950);
\r
448 // SetMaxNoiseAbs(2);
\r
449 // SetMaxNoiseSigma(3);
\r
452 AliTPCDigitsArray *digarr=new AliTPCDigitsArray(kFALSE);//data not sim
\r
453 digarr->Setup(fParam);//as usually parameters
\r
456 // Loop over sectors
\r
458 for(Int_t iSec = 0; iSec < kNS; iSec++) {
\r
459 AliTPCCalROC * noiseROC;
\r
461 noiseROC =new AliTPCCalROC(iSec);//noise=0
\r
464 noiseROC = noiseTPC->GetCalROC(iSec); // noise per given sector
\r
466 Int_t nRows = 0; //number of rows in sector
\r
467 Int_t nDDLs = 0; //number of DDLs
\r
468 Int_t indexDDL = 0; //DDL index
\r
470 nRows = fParam->GetNRowLow();
\r
472 indexDDL = iSec * 2;
\r
474 nRows = fParam->GetNRowUp();
\r
476 indexDDL = (iSec-kNIS) * 4 + kNIS * 2;
\r
480 // Load the raw data for corresponding DDLs
\r
482 rawReader->Reset();
\r
483 rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1);
\r
487 // Allocate memory for rows in sector (pads(depends on row) x timebins)
\r
488 for(Int_t iRow = 0; iRow < nRows; iRow++) {
\r
489 digarr->CreateRow(iSec,iRow);
\r
490 }//end loop over rows
\r
492 rawReader->Reset();
\r
493 rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1);
\r
496 // Begin loop over altro data
\r
498 while (input.Next()) {
\r
500 //check sector consistency
\r
501 if (input.GetSector() != iSec)
\r
502 AliFatal(Form("Sector index mismatch ! Expected (%d), but got (%d) !",iSec,input.GetSector()));
\r
504 Short_t iRow = input.GetRow();
\r
505 Short_t iPad = input.GetPad();
\r
506 Short_t iTimeBin = input.GetTime();
\r
508 if(fDebugLevel==72){
\r
509 fHistoRow->Fill(iRow);
\r
510 fHistoPad->Fill(iPad);
\r
511 fHistoTime->Fill(iTimeBin);
\r
512 fHistoRowPad->Fill(iPad,iRow);
\r
513 }else if(fDebugLevel>=0&&fDebugLevel<72){
\r
514 if(iSec==fDebugLevel){
\r
515 fHistoRow->Fill(iRow);
\r
516 fHistoPad->Fill(iPad);
\r
517 fHistoTime->Fill(iTimeBin);
\r
518 fHistoRowPad->Fill(iPad,iRow);
\r
520 }else if(fDebugLevel==73){
\r
522 fHistoRow->Fill(iRow);
\r
523 fHistoPad->Fill(iPad);
\r
524 fHistoTime->Fill(iTimeBin);
\r
525 fHistoRowPad->Fill(iPad,iRow);
\r
527 }else if(fDebugLevel==74){
\r
529 fHistoRow->Fill(iRow);
\r
530 fHistoPad->Fill(iPad);
\r
531 fHistoTime->Fill(iTimeBin);
\r
532 fHistoRowPad->Fill(iPad,iRow);
\r
536 //check row consistency
\r
537 if (iRow < 0 || iRow >= nRows){
\r
538 AliError(Form("Pad-row index (%d) outside the range (%d -> %d) !",
\r
539 iRow, 0, nRows -1));
\r
543 //check pad consistency
\r
544 if (iPad < 0 || iPad >= (Short_t)(roc->GetNPads(iSec,iRow))) {
\r
545 AliError(Form("Pad index (%d) outside the range (%d -> %d) !",
\r
546 iPad, 0, roc->GetNPads(iSec,iRow) ));
\r
550 //check time consistency
\r
551 if ( iTimeBin < fRecoParam->GetFirstBin() || iTimeBin >= fRecoParam->GetLastBin()){
\r
552 //cout<<iTimeBin<<endl;
\r
554 AliFatal(Form("Timebin index (%d) outside the range (%d -> %d) !",
\r
555 iTimeBin, 0, fRecoParam->GetLastBin() -1));
\r
559 Int_t signal = input.GetSignal();
\r
560 if (signal <= fZeroSup ||
\r
561 iTimeBin < fFirstBin ||
\r
562 iTimeBin > fLastBin
\r
564 digarr->GetRow(iSec,iRow)->SetDigitFast(0,iTimeBin,iPad);
\r
567 if (!noiseROC) continue;
\r
568 Double_t noiseOnPad = noiseROC->GetValue(iRow,iPad);//noise on given pad and row in sector
\r
569 if (noiseOnPad > fMaxNoiseAbs){
\r
570 digarr->GetRow(iSec,iRow)->SetDigitFast(0,iTimeBin,iPad);
\r
571 continue; // consider noisy pad as dead
\r
573 if(signal <= fMaxNoiseSigma * noiseOnPad){
\r
574 digarr->GetRow(iSec,iRow)->SetDigitFast(0,iTimeBin,iPad);
\r
578 digarr->GetRow(iSec,iRow)->SetDigitFast(signal,iTimeBin,iPad);
\r
580 }//end of loop over altro data
\r
581 }//end of loop over sectors
\r
584 if(isAltro) FindClusterKrIO();
\r
590 ////____________________________________________________________________________
\r
591 Int_t AliTPCclustererKr::FindClusterKrIO()
\r
593 //fParam and fDigarr must be set to run this method
\r
595 // //set parameters
\r
596 // SetMinAdc(3);//usually is 3
\r
597 // SetMinTimeBins(2);//should be 2 - the best result of shape in MC
\r
598 //// SetMaxPadRange(4);
\r
599 //// SetMaxRowRange(3);
\r
600 // SetMaxTimeRange(7);
\r
601 // SetValueToSize(3.5);//3.5
\r
602 // SetMaxPadRangeCm(2.5);
\r
603 // SetMaxRowRangeCm(3.5);
\r
605 Int_t clusterCounter=0;
\r
606 const Short_t nTotalSector=fParam->GetNSector();//number of sectors
\r
607 for(Short_t iSec=0; iSec<nTotalSector; ++iSec){
\r
609 //vector of maxima for each sector
\r
610 //std::vector<AliPadMax*> maximaInSector;
\r
611 TObjArray *maximaInSector=new TObjArray();//to store AliPadMax*
\r
614 // looking for the maxima on the pad
\r
617 const Short_t kNRows=fParam->GetNRow(iSec);//number of rows in sector
\r
618 for(Short_t iRow=0; iRow<kNRows; ++iRow){
\r
619 AliSimDigits *digrow;
\r
621 digrow = (AliSimDigits*)fDigarr->GetRow(iSec,iRow);//real data
\r
623 digrow = (AliSimDigits*)fDigarr->LoadRow(iSec,iRow);//MC
\r
625 if(digrow){//if pointer exist
\r
626 digrow->ExpandBuffer(); //decrunch
\r
627 const Short_t kNPads = digrow->GetNCols(); // number of pads
\r
628 const Short_t kNTime = digrow->GetNRows(); // number of timebins
\r
629 for(Short_t iPad=0;iPad<kNPads;iPad++){
\r
631 Short_t timeBinMax=-1;//timebin of maximum
\r
632 Short_t valueMaximum=-1;//value of maximum in adc
\r
633 Short_t increaseBegin=-1;//timebin when increase starts
\r
634 Short_t sumAdc=0;//sum of adc on the pad in maximum surrounding
\r
635 bool ifIncreaseBegin=true;//flag - check if increasing started
\r
636 bool ifMaximum=false;//flag - check if it could be maximum
\r
638 for(Short_t iTimeBin=0;iTimeBin<kNTime;iTimeBin++){
\r
639 Short_t adc = digrow->GetDigitFast(iTimeBin,iPad);
\r
640 if(adc<fMinAdc){//standard was 3
\r
642 if(iTimeBin-increaseBegin<fMinTimeBins){//at least 2 time bins
\r
647 ifIncreaseBegin=true;
\r
651 //insert maximum, default values and set flags
\r
652 Double_t xCord,yCord;
\r
653 GetXY(iSec,iRow,iPad,xCord,yCord);
\r
654 AliPadMax *oneMaximum = new AliPadMax(AliTPCvtpr(valueMaximum,
\r
664 //maximaInSector.push_back(oneMaximum);
\r
665 maximaInSector->AddLast(oneMaximum);
\r
671 ifIncreaseBegin=true;
\r
677 if(ifIncreaseBegin){
\r
678 ifIncreaseBegin=false;
\r
679 increaseBegin=iTimeBin;
\r
682 if(adc>valueMaximum){
\r
683 timeBinMax=iTimeBin;
\r
688 if(iTimeBin==kNTime-1 && ifMaximum && kNTime-increaseBegin>fMinTimeBins){//on the edge
\r
689 //at least 3 timebins
\r
690 //insert maximum, default values and set flags
\r
691 Double_t xCord,yCord;
\r
692 GetXY(iSec,iRow,iPad,xCord,yCord);
\r
693 AliPadMax *oneMaximum = new AliPadMax(AliTPCvtpr(valueMaximum,
\r
703 //maximaInSector.push_back(oneMaximum);
\r
704 maximaInSector->AddLast(oneMaximum);
\r
710 ifIncreaseBegin=true;
\r
715 }//end loop over timebins
\r
716 }//end loop over pads
\r
718 // cout<<"Pointer does not exist!!"<<endl;
\r
719 }//end if poiner exists
\r
720 }//end loop over rows
\r
722 // cout<<"EF" <<maximaInSector->GetEntriesFast()<<" E "<<maximaInSector->GetEntries()<<" "<<maximaInSector->GetLast()<<endl;
\r
723 maximaInSector->Compress();
\r
724 // GetEntriesFast() - liczba wejsc w array of maxima
\r
725 //cout<<"EF" <<maximaInSector->GetEntriesFast()<<" E"<<maximaInSector->GetEntries()<<endl;
\r
732 Short_t maxSumAdc=0;
\r
733 Short_t maxTimeBin=0;
\r
739 // for( std::vector<AliPadMax*>::iterator mp1 = maximaInSector.begin();
\r
740 // mp1 != maximaInSector.end(); ++mp1 ) {
\r
741 for(Int_t it1 = 0; it1 < maximaInSector->GetEntriesFast(); ++it1 ) {
\r
743 AliPadMax *mp1=(AliPadMax *)maximaInSector->At(it1);
\r
744 AliTPCclusterKr *tmp=new AliTPCclusterKr();
\r
746 Short_t nUsedPads=1;
\r
747 Int_t clusterValue=0;
\r
748 clusterValue+=(mp1)->GetSum();
\r
749 list<Short_t> nUsedRows;
\r
750 nUsedRows.push_back((mp1)->GetRow());
\r
752 maxDig =(mp1)->GetAdc() ;
\r
753 maxSumAdc =(mp1)->GetSum() ;
\r
754 maxTimeBin =(mp1)->GetTime();
\r
755 maxPad =(mp1)->GetPad() ;
\r
756 maxRow =(mp1)->GetRow() ;
\r
757 maxX =(mp1)->GetX();
\r
758 maxY =(mp1)->GetY();
\r
760 AliSimDigits *digrowTmp;
\r
762 digrowTmp = (AliSimDigits*)fDigarr->GetRow(iSec,(mp1)->GetRow());
\r
764 digrowTmp = (AliSimDigits*)fDigarr->LoadRow(iSec,(mp1)->GetRow());
\r
767 digrowTmp->ExpandBuffer(); //decrunch
\r
769 for(Short_t itb=(mp1)->GetBegin(); itb<((mp1)->GetEnd())+1; itb++){
\r
770 Short_t adcTmp = digrowTmp->GetDigitFast(itb,(mp1)->GetPad());
\r
771 AliTPCvtpr *vtpr=new AliTPCvtpr(adcTmp,itb,(mp1)->GetPad(),(mp1)->GetRow(),(mp1)->GetX(),(mp1)->GetY(),itb);
\r
772 //tmp->fCluster.push_back(vtpr);
\r
773 tmp->AddDigitToCluster(vtpr);
\r
775 tmp->SetCenter();//set centr of the cluster
\r
777 //maximaInSector.erase(mp1);
\r
779 //cout<<maximaInSector->GetEntriesFast()<<" ";
\r
780 maximaInSector->RemoveAt(it1);
\r
781 maximaInSector->Compress();
\r
783 // cout<<maximaInSector->GetEntriesFast()<<" "<<endl;
\r
785 // for( std::vector<AliPadMax*>::iterator mp2 = maximaInSector.begin();
\r
786 // mp2 != maximaInSector.end(); ++mp2 ) {
\r
787 for(Int_t it2 = 0; it2 < maximaInSector->GetEntriesFast(); ++it2 ) {
\r
788 AliPadMax *mp2=(AliPadMax *)maximaInSector->At(it2);
\r
790 // if(abs(maxRow - (*mp2)->GetRow()) < fMaxRowRange && //3
\r
791 // abs(maxPad - (*mp2)->GetPad()) < fMaxPadRange && //4
\r
793 if(TMath::Abs(tmp->GetCenterX() - (mp2)->GetX()) < fMaxPadRangeCm &&
\r
794 TMath::Abs(tmp->GetCenterY() - (mp2)->GetY()) < fMaxRowRangeCm &&
\r
795 TMath::Abs(tmp->GetCenterT() - (mp2)->GetT()) < fMaxTimeRange){//7
\r
797 clusterValue+=(mp2)->GetSum();
\r
800 nUsedRows.push_back((mp2)->GetRow());
\r
802 AliSimDigits *digrowTmp1;
\r
804 digrowTmp1 = (AliSimDigits*)fDigarr->GetRow(iSec,(mp2)->GetRow());
\r
806 digrowTmp1 = (AliSimDigits*)fDigarr->LoadRow(iSec,(mp2)->GetRow());
\r
808 digrowTmp1->ExpandBuffer(); //decrunch
\r
810 for(Short_t itb=(mp2)->GetBegin(); itb<(mp2)->GetEnd()+1; itb++){
\r
811 Short_t adcTmp = digrowTmp1->GetDigitFast(itb,(mp2)->GetPad());
\r
812 AliTPCvtpr *vtpr=new AliTPCvtpr(adcTmp,itb,(mp2)->GetPad(),(mp2)->GetRow(),(mp2)->GetX(),(mp2)->GetY(),itb);
\r
813 //tmp->fCluster.push_back(vtpr);
\r
814 tmp->AddDigitToCluster(vtpr);
\r
817 tmp->SetCenter();//set center of the cluster
\r
819 //which one is bigger
\r
820 if( (mp2)->GetAdc() > maxDig ){
\r
821 maxDig =(mp2)->GetAdc() ;
\r
822 maxSumAdc =(mp2)->GetSum() ;
\r
823 maxTimeBin =(mp2)->GetTime();
\r
824 maxPad =(mp2)->GetPad() ;
\r
825 maxRow =(mp2)->GetRow() ;
\r
826 maxX =(mp2)->GetX() ;
\r
827 maxY =(mp2)->GetY() ;
\r
828 } else if ( (mp2)->GetAdc() == maxDig ){
\r
829 if( (mp2)->GetSum() > maxSumAdc){
\r
830 maxDig =(mp2)->GetAdc() ;
\r
831 maxSumAdc =(mp2)->GetSum() ;
\r
832 maxTimeBin =(mp2)->GetTime();
\r
833 maxPad =(mp2)->GetPad() ;
\r
834 maxRow =(mp2)->GetRow() ;
\r
835 maxX =(mp2)->GetX() ;
\r
836 maxY =(mp2)->GetY() ;
\r
839 maximaInSector->RemoveAt(it2);
\r
840 maximaInSector->Compress();
\r
843 //maximaInSector.erase(mp2);
\r
849 //through out ADC=6,7 on 1 pad, 2 tb and ADC=12 on 2 pads,2 tb
\r
850 //if(nUsedPads==1 && clusterValue/tmp->fCluster.size()<3.6)continue;
\r
851 //if(nUsedPads==2 && clusterValue/tmp->fCluster.size()<3.1)continue;
\r
853 //through out clusters on the edge and noise
\r
854 //if(clusterValue/tmp->fCluster.size()<fValueToSize)continue;
\r
855 if(clusterValue/(tmp->GetSize())<fValueToSize)continue;
\r
857 tmp->SetADCcluster(clusterValue);
\r
858 tmp->SetNPads(nUsedPads);
\r
859 tmp->SetMax(AliTPCvtpr(maxDig,maxTimeBin,maxPad,maxRow,maxX,maxY,maxTimeBin));
\r
864 nUsedRows.unique();
\r
865 tmp->SetNRows(nUsedRows.size());
\r
871 //save each cluster into file
\r
873 AliTPCClustersRow *clrow= new AliTPCClustersRow();
\r
875 clrow->SetClass("AliTPCclusterKr");
\r
876 clrow->SetArray(1);
\r
877 fOutput->GetBranch("Segment")->SetAddress(&clrow);
\r
879 Int_t tmpCluster=0;
\r
880 TClonesArray * arr = fRowCl->GetArray();
\r
881 AliTPCclusterKr * cl = new ((*arr)[tmpCluster]) AliTPCclusterKr(*tmp);
\r
885 //end of save each cluster into file adc.root
\r
889 cout<<"Number of clusters in event: "<<clusterCounter<<endl;
\r
894 ////____________________________________________________________________________
\r
897 void AliTPCclustererKr::GetXY(Short_t sec,Short_t row,Short_t pad,Double_t& xGlob,Double_t& yGlob){
\r
899 //gives global XY coordinate of the pad
\r
902 Double_t yLocal = fParam->GetPadRowRadii(sec,row);//radius of row in sector in cm
\r
904 Short_t padmax = fParam->GetNPads(sec,row);//number of pads in a given row
\r
906 if(sec<fParam->GetNInnerSector())padXSize=0.4;
\r
908 Double_t xLocal=(pad+0.5-padmax/2.)*padXSize;//x-value of the center of pad
\r
911 fParam->AdjustCosSin((Int_t)sec,cos,sin);//return sinus and cosinus of the sector
\r
913 Double_t xGlob1 = xLocal * cos + yLocal * sin;
\r
914 Double_t yGlob1 = -xLocal * sin + yLocal * cos;
\r
918 rot=TMath::Pi()/2.;
\r
920 xGlob = xGlob1 * TMath::Cos(rot) + yGlob1 * TMath::Sin(rot);
\r
921 yGlob = -xGlob1 * TMath::Sin(rot) + yGlob1 * TMath::Cos(rot);
\r
924 if(sec<18||(sec>=36&&sec<54)) xGlob =-1*xGlob;
\r