AliTPCclustererKr.h - Adding debugging histograms
[u/mrichter/AliRoot.git] / TPC / AliTPCclustererKr.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /* $Id: AliTPCclustererKr.cxx,v 1.7 2008/02/06 17:24:53 matyja Exp $ */
17
18 //-----------------------------------------------------------------
19 //           Implementation of the TPC Kr cluster class
20 //
21 // Origin: Adam Matyja, INP PAN, adam.matyja@ifj.edu.pl
22 //-----------------------------------------------------------------
23
24 /*
25 Instruction - how to use that:
26 There are two macros prepared. One is for preparing clusters from MC 
27 samples:  FindKrClusters.C. The output is kept in TPC.RecPoints.root.
28 The other macro is prepared for data analysis: FindKrClustersRaw.C. 
29 The output is created for each processed file in root file named adc.root. 
30 For each data subsample the same named file is created. So be careful 
31 do not overwrite them. 
32
33 *
34 **** MC ****
35 *
36
37 To run clusterizaton for MC type:
38 .x FindKrClusters.C
39
40 If you don't want to use the standard selection criteria then you 
41 have to do following:
42
43 // load the standard setup
44 AliRunLoader* rl = AliRunLoader::Open("galice.root");
45 AliTPCLoader *tpcl = (AliTPCLoader*)rl->GetLoader("TPCLoader");
46 tpcl->LoadDigits();
47 rl->LoadgAlice();
48 gAlice=rl->GetAliRun();
49 TDirectory *cwd = gDirectory;
50 AliTPCv4 *tpc = (AliTPCv4*)gAlice->GetDetector("TPC");
51 Int_t ver = tpc->IsVersion();
52 rl->CdGAFile();
53 AliTPCParam *param=(AliTPCParamSR *)gDirectory->Get("75x40_100x60_150x60");
54 AliTPCDigitsArray *digarr=new AliTPCDigitsArray;
55 digarr->Setup(param);
56 cwd->cd();
57
58 //loop over events
59 Int_t nevmax=rl->GetNumberOfEvents();
60 for(Int_t nev=0;nev<nevmax ;nev++){
61   rl->GetEvent(nev);
62   TTree* input_tree= tpcl->TreeD();//load tree with digits
63   digarr->ConnectTree(input_tree);
64   TTree *output_tree =tpcl->TreeR();//load output tree
65
66   AliTPCclustererKr *clusters = new AliTPCclustererKr();
67   clusters->SetParam(param);
68   clusters->SetInput(input_tree);
69   clusters->SetOutput(output_tree);
70   clusters->SetDigArr(digarr);
71   
72 //If you want to change the cluster finder parameters for MC there are 
73 //several of them:
74
75 //1. signal threshold (everything below the given number is treated as 0)
76   clusters->SetMinAdc(3);
77
78 //2. number of neighbouring timebins to be considered
79   clusters->SetMinTimeBins(2);
80
81 //3. distance of the cluster center to the center of a pad in pad-padrow plane 
82 //(in cm). Remenber that this is still quantified by pad size.
83   clusters->SetMaxPadRangeCm(2.5);
84
85 //4. distance of the cluster center to the center of a padrow in pad-padrow 
86 //plane (in cm). Remenber that this is still quantified by pad size.
87   clusters->SetMaxRowRangeCm(3.5);
88
89 //5. distance of the cluster center to the max time bin on a pad (in tackts)
90 //ie. fabs(centerT - time)<7
91   clusters->SetMaxTimeRange(7);
92
93 //6. cut reduce peak at 0. There are noises which appear mostly as two 
94 //timebins on one pad.
95   clusters->SetValueToSize(3.5);
96
97
98   clusters->finderIO();
99   tpcl->WriteRecPoints("OVERWRITE");
100 }
101 delete rl;//cleans everything
102
103 *
104 ********* DATA *********
105 *
106
107 To run clusterizaton for DATA for file named raw_data.root type:
108 .x FindKrClustersRaw.C("raw_data.root")
109
110 If you want to change some criteria do the following:
111
112 //
113 // remove Altro warnings
114 //
115 AliLog::SetClassDebugLevel("AliTPCRawStream",-5);
116 AliLog::SetClassDebugLevel("AliRawReaderDate",-5);
117 AliLog::SetClassDebugLevel("AliTPCAltroMapping",-5);
118 AliLog::SetModuleDebugLevel("RAW",-5);
119
120 //
121 // Get database with noises
122 //
123 //  char *ocdbpath = gSystem->Getenv("OCDB_PATH");
124 char *ocdbpath ="local:///afs/cern.ch/alice/tpctest/OCDB";
125 if (ocdbpath==0){
126 ocdbpath="alien://folder=/alice/data/2007/LHC07w/OCDB/";
127 }
128 AliCDBManager * man = AliCDBManager::Instance();
129 man->SetDefaultStorage(ocdbpath);
130 man->SetRun(0);
131 AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
132 AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping();
133
134 //define tree
135 TFile *hfile=new TFile("adc.root","RECREATE","ADC file");
136 // Create a ROOT Tree
137 TTree *mytree = new TTree("Kr","Krypton cluster tree");
138
139 //define infput file
140 const char *fileName="data.root";
141 AliRawReader *reader = new AliRawReaderRoot(fileName);
142 //AliRawReader *reader = new AliRawReaderDate(fileName);
143 reader->Reset();
144 AliAltroRawStreamFast* stream = new AliAltroRawStreamFast(reader);
145 stream->SelectRawData("TPC");
146
147 //one general output
148 AliTPCclustererKr *clusters = new AliTPCclustererKr();
149 clusters->SetOutput(mytree);
150 clusters->SetRecoParam(0);//standard reco parameters
151 AliTPCParamSR *param=new AliTPCParamSR();
152 clusters->SetParam(param);//TPC parameters(sectors, timebins, etc.)
153
154 //set cluster finder parameters (from data):
155 //1. zero suppression parameter
156   clusters->SetZeroSup(param->GetZeroSup());
157
158 //2. first bin
159   clusters->SetFirstBin(60);
160
161 //3. last bin
162   clusters->SetLastBin(950);
163
164 //4. maximal noise
165   clusters->SetMaxNoiseAbs(2);
166
167 //5. maximal amount of sigma of noise
168   clusters->SetMaxNoiseSigma(3);
169
170 //The remaining parameters are the same paramters as for MC (see MC section 
171 //points 1-6)
172   clusters->SetMinAdc(3);
173   clusters->SetMinTimeBins(2);
174   clusters->SetMaxPadRangeCm(2.5);
175   clusters->SetMaxRowRangeCm(3.5);
176   clusters->SetMaxTimeRange(7);
177   clusters->SetValueToSize(3.5);
178
179 while (reader->NextEvent()) {
180   clusters->FinderIO(reader);
181 }
182
183 hfile->Write();
184 hfile->Close();
185 delete stream;
186
187
188 */
189
190 #include "AliTPCclustererKr.h"
191 #include "AliTPCclusterKr.h"
192 //#include <vector>
193 #include <list>
194 #include "TObject.h"
195 #include "AliPadMax.h"
196 #include "AliSimDigits.h"
197 #include "AliTPCv4.h"
198 #include "AliTPCParam.h"
199 #include "AliTPCDigitsArray.h"
200 #include "AliTPCvtpr.h"
201 #include "AliTPCClustersRow.h"
202 #include "TTree.h"
203 #include "TH1F.h"
204 #include "TH2F.h"
205
206 //used in raw data finder
207 #include "AliTPCROC.h"
208 #include "AliTPCCalPad.h"
209 #include "AliTPCAltroMapping.h"
210 #include "AliTPCcalibDB.h"
211 #include "AliTPCRawStream.h"
212 #include "AliTPCRecoParam.h"
213 #include "AliTPCReconstructor.h"
214 #include "AliRawReader.h"
215 #include "AliTPCCalROC.h"
216
217 ClassImp(AliTPCclustererKr)
218
219
220 AliTPCclustererKr::AliTPCclustererKr()
221   :TObject(),
222   fRawData(kFALSE),
223   fRowCl(0),
224   fInput(0),
225   fOutput(0),
226   fParam(0),
227   fDigarr(0),
228   fRecoParam(0),
229   fZeroSup(2),
230   fFirstBin(60),
231   fLastBin(950),
232   fMaxNoiseAbs(2),
233   fMaxNoiseSigma(3),
234   fMinAdc(3),
235   fMinTimeBins(2),
236 //  fMaxPadRange(4),
237 //  fMaxRowRange(3),
238   fMaxTimeRange(7),
239   fValueToSize(3.5),
240   fMaxPadRangeCm(2.5),
241   fMaxRowRangeCm(3.5),
242   fDebugLevel(-1),
243   fHistoRow(0),
244   fHistoPad(0),
245   fHistoTime(0),
246   fHistoRowPad(0)
247 {
248 //
249 // default constructor
250 //
251 }
252
253 AliTPCclustererKr::AliTPCclustererKr(const AliTPCclustererKr &param)
254   :TObject(),
255   fRawData(kFALSE),
256   fRowCl(0),
257   fInput(0),
258   fOutput(0),
259   fParam(0),
260   fDigarr(0),
261   fRecoParam(0),
262   fZeroSup(2),
263   fFirstBin(60),
264   fLastBin(950),
265   fMaxNoiseAbs(2),
266   fMaxNoiseSigma(3),
267   fMinAdc(3),
268   fMinTimeBins(2),
269 //  fMaxPadRange(4),
270 //  fMaxRowRange(3),
271   fMaxTimeRange(7),
272   fValueToSize(3.5),
273   fMaxPadRangeCm(2.5),
274   fMaxRowRangeCm(3.5),
275   fDebugLevel(-1),
276   fHistoRow(0),
277   fHistoPad(0),
278   fHistoTime(0),
279   fHistoRowPad(0)
280 {
281 //
282 // copy constructor
283 //
284   fParam = param.fParam;
285   fRecoParam = param.fRecoParam;
286   fRawData = param.fRawData;
287   fRowCl  = param.fRowCl ;
288   fInput  = param.fInput ;
289   fOutput = param.fOutput;
290   fDigarr = param.fDigarr;
291   fZeroSup       = param.fZeroSup       ;
292   fFirstBin      = param.fFirstBin      ;
293   fLastBin       = param.fLastBin       ;
294   fMaxNoiseAbs   = param.fMaxNoiseAbs   ;
295   fMaxNoiseSigma = param.fMaxNoiseSigma ;
296   fMinAdc = param.fMinAdc;
297   fMinTimeBins = param.fMinTimeBins;
298 //  fMaxPadRange  = param.fMaxPadRange ;
299 //  fMaxRowRange  = param.fMaxRowRange ;
300   fMaxTimeRange = param.fMaxTimeRange;
301   fValueToSize  = param.fValueToSize;
302   fMaxPadRangeCm = param.fMaxPadRangeCm;
303   fMaxRowRangeCm = param.fMaxRowRangeCm;
304   fDebugLevel = param.fDebugLevel;
305   fHistoRow    = param.fHistoRow   ;
306   fHistoPad    = param.fHistoPad  ;
307   fHistoTime   = param.fHistoTime;
308   fHistoRowPad = param.fHistoRowPad;
309
310
311
312 AliTPCclustererKr & AliTPCclustererKr::operator = (const AliTPCclustererKr & param)
313 {
314   //
315   // assignment operator
316   //
317   fParam = param.fParam;
318   fRecoParam = param.fRecoParam;
319   fRawData = param.fRawData;
320   fRowCl  = param.fRowCl ;
321   fInput  = param.fInput ;
322   fOutput = param.fOutput;
323   fDigarr = param.fDigarr;
324   fZeroSup       = param.fZeroSup       ;
325   fFirstBin      = param.fFirstBin      ;
326   fLastBin       = param.fLastBin       ;
327   fMaxNoiseAbs   = param.fMaxNoiseAbs   ;
328   fMaxNoiseSigma = param.fMaxNoiseSigma ;
329   fMinAdc = param.fMinAdc;
330   fMinTimeBins = param.fMinTimeBins;
331 //  fMaxPadRange  = param.fMaxPadRange ;
332 //  fMaxRowRange  = param.fMaxRowRange ;
333   fMaxTimeRange = param.fMaxTimeRange;
334   fValueToSize  = param.fValueToSize;
335   fMaxPadRangeCm = param.fMaxPadRangeCm;
336   fMaxRowRangeCm = param.fMaxRowRangeCm;
337   fDebugLevel = param.fDebugLevel;
338   fHistoRow    = param.fHistoRow   ;
339   fHistoPad    = param.fHistoPad  ;
340   fHistoTime   = param.fHistoTime;
341   fHistoRowPad = param.fHistoRowPad;
342   return (*this);
343 }
344
345 AliTPCclustererKr::~AliTPCclustererKr()
346 {
347   //
348   // destructor
349   //
350 }
351
352 void AliTPCclustererKr::SetRecoParam(AliTPCRecoParam *recoParam)
353 {
354   //
355   // set reconstruction parameters
356   //
357   if (recoParam) {
358     fRecoParam = recoParam;
359   }else{
360     //set default parameters if not specified
361     fRecoParam = AliTPCReconstructor::GetRecoParam();
362     if (!fRecoParam)  fRecoParam = AliTPCRecoParam::GetLowFluxParam();
363   }
364   return;
365 }
366
367
368 ////____________________________________________________________________________
369 ////       I/O
370 void AliTPCclustererKr::SetInput(TTree * tree)
371 {
372   //
373   // set input tree with digits
374   //
375   fInput = tree;  
376   if  (!fInput->GetBranch("Segment")){
377     cerr<<"AliTPCclusterKr::FindClusterKr(): no proper input tree !\n";
378     fInput=0;
379     return;
380   }
381 }
382
383 void AliTPCclustererKr::SetOutput(TTree * tree) 
384 {
385   //
386   // set output
387   //
388   fOutput= tree;  
389   AliTPCClustersRow clrow;
390   AliTPCClustersRow *pclrow=&clrow;  
391   clrow.SetClass("AliTPCclusterKr");
392   clrow.SetArray(1); // to make Clones array
393   fOutput->Branch("Segment","AliTPCClustersRow",&pclrow,32000,200);    
394 }
395
396 ////____________________________________________________________________________
397 //// with new I/O
398 Int_t AliTPCclustererKr::FinderIO()
399 {
400   // Krypton cluster finder for simulated events from MC
401
402   if (!fInput) { 
403     Error("Digits2Clusters", "input tree not initialised");
404     return 10;
405   }
406   
407   if (!fOutput) {
408     Error("Digits2Clusters", "output tree not initialised");
409     return 11;
410   }
411
412   FindClusterKrIO();
413   return 0;
414 }
415
416 Int_t AliTPCclustererKr::FinderIO(AliRawReader* rawReader)
417 {
418   // Krypton cluster finder for the TPC raw data
419   //
420   // fParam must be defined before
421
422   if(rawReader)fRawData=kTRUE; //set flag to data
423
424   if (!fOutput) {
425     Error("Digits2Clusters", "output tree not initialised");
426     return 11;
427   }
428
429   fParam->SetMaxTBin(fRecoParam->GetLastBin());//set number of timebins from reco -> param
430   //   used later for memory allocation
431
432   Bool_t isAltro=kFALSE;
433
434   AliTPCROC * roc = AliTPCROC::Instance();
435   AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
436   AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping();
437   //
438   AliTPCRawStream input(rawReader,(AliAltroMapping**)mapping);
439
440   const Int_t kNIS = fParam->GetNInnerSector();//number of inner sectors
441   const Int_t kNOS = fParam->GetNOuterSector();//number of outer sectors
442   const Int_t kNS = kNIS + kNOS;//all sectors
443
444 //  //set cluster finder parameters
445 //  SetZeroSup(fParam->GetZeroSup());//zero suppression parameter
446 //  SetFirstBin(60);
447 //  SetLastBin(950);
448 //  SetMaxNoiseAbs(2);
449 //  SetMaxNoiseSigma(3);
450
451   //crate TPC view
452   AliTPCDigitsArray *digarr=new AliTPCDigitsArray(kFALSE);//data not sim
453   digarr->Setup(fParam);//as usually parameters
454
455   //
456   // Loop over sectors
457   //
458   for(Int_t iSec = 0; iSec < kNS; iSec++) {
459     AliTPCCalROC * noiseROC;
460     if(noiseTPC==0x0){
461       noiseROC =new AliTPCCalROC(iSec);//noise=0
462     }
463     else{
464       noiseROC = noiseTPC->GetCalROC(iSec);  // noise per given sector
465     }
466     Int_t nRows = 0; //number of rows in sector
467     Int_t nDDLs = 0; //number of DDLs
468     Int_t indexDDL = 0; //DDL index
469     if (iSec < kNIS) {
470       nRows = fParam->GetNRowLow();
471       nDDLs = 2;
472       indexDDL = iSec * 2;
473     }else {
474       nRows = fParam->GetNRowUp();
475       nDDLs = 4;
476       indexDDL = (iSec-kNIS) * 4 + kNIS * 2;
477     }
478
479     //
480     // Load the raw data for corresponding DDLs
481     //
482     rawReader->Reset();
483     rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1);
484
485     if(input.Next()) {
486       isAltro=kTRUE;
487       // Allocate memory for rows in sector (pads(depends on row) x timebins)
488       for(Int_t iRow = 0; iRow < nRows; iRow++) {
489         digarr->CreateRow(iSec,iRow);
490       }//end loop over rows
491     }
492     rawReader->Reset();
493     rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1);
494
495     //
496     // Begin loop over altro data
497     //
498     while (input.Next()) {
499
500       //check sector consistency
501       if (input.GetSector() != iSec)
502         AliFatal(Form("Sector index mismatch ! Expected (%d), but got (%d) !",iSec,input.GetSector()));
503       
504       Short_t iRow = input.GetRow();
505       Short_t iPad = input.GetPad();
506       Short_t iTimeBin = input.GetTime();
507
508       if(fDebugLevel==72){
509         fHistoRow->Fill(iRow);
510         fHistoPad->Fill(iPad);
511         fHistoTime->Fill(iTimeBin);
512         fHistoRowPad->Fill(iPad,iRow);
513       }else if(fDebugLevel>=0&&fDebugLevel<72){
514         if(iSec==fDebugLevel){
515           fHistoRow->Fill(iRow);
516           fHistoPad->Fill(iPad);
517           fHistoTime->Fill(iTimeBin);
518           fHistoRowPad->Fill(iPad,iRow);
519         }
520       }else if(fDebugLevel==73){
521         if(iSec<36){
522           fHistoRow->Fill(iRow);
523           fHistoPad->Fill(iPad);
524           fHistoTime->Fill(iTimeBin);
525           fHistoRowPad->Fill(iPad,iRow);
526         }
527       }else if(fDebugLevel==74){
528         if(iSec>=36){
529           fHistoRow->Fill(iRow);
530           fHistoPad->Fill(iPad);
531           fHistoTime->Fill(iTimeBin);
532           fHistoRowPad->Fill(iPad,iRow);
533         }
534       }
535
536       //check row consistency
537       if (iRow < 0 || iRow >= nRows){
538         AliError(Form("Pad-row index (%d) outside the range (%d -> %d) !",
539                       iRow, 0, nRows -1));
540         continue;
541       }
542
543       //check pad consistency
544       if (iPad < 0 || iPad >= (Short_t)(roc->GetNPads(iSec,iRow))) {
545         AliError(Form("Pad index (%d) outside the range (%d -> %d) !",
546                       iPad, 0, roc->GetNPads(iSec,iRow) ));
547         continue;
548       }
549
550       //check time consistency
551       if ( iTimeBin < fRecoParam->GetFirstBin() || iTimeBin >= fRecoParam->GetLastBin()){
552         //cout<<iTimeBin<<endl;
553         continue;
554         AliFatal(Form("Timebin index (%d) outside the range (%d -> %d) !",
555                       iTimeBin, 0, fRecoParam->GetLastBin() -1));
556       }
557
558       //signal
559       Int_t signal = input.GetSignal();
560       if (signal <= fZeroSup || 
561           iTimeBin < fFirstBin ||
562           iTimeBin > fLastBin
563           ) {
564         digarr->GetRow(iSec,iRow)->SetDigitFast(0,iTimeBin,iPad);
565         continue;
566       }
567       if (!noiseROC) continue;
568       Double_t noiseOnPad = noiseROC->GetValue(iRow,iPad);//noise on given pad and row in sector
569       if (noiseOnPad > fMaxNoiseAbs){
570         digarr->GetRow(iSec,iRow)->SetDigitFast(0,iTimeBin,iPad);
571         continue; // consider noisy pad as dead
572       }
573       if(signal <= fMaxNoiseSigma * noiseOnPad){
574         digarr->GetRow(iSec,iRow)->SetDigitFast(0,iTimeBin,iPad);
575         continue;
576       }
577
578       digarr->GetRow(iSec,iRow)->SetDigitFast(signal,iTimeBin,iPad);
579
580     }//end of loop over altro data
581   }//end of loop over sectors
582   
583   SetDigArr(digarr);
584   if(isAltro) FindClusterKrIO();
585   delete digarr;
586
587   return 0;
588 }
589
590 ////____________________________________________________________________________
591 Int_t AliTPCclustererKr::FindClusterKrIO()
592 {
593   //fParam and  fDigarr must be set to run this method
594
595 //  //set parameters 
596 //  SetMinAdc(3);//usually is 3 
597 //  SetMinTimeBins(2);//should be 2 - the best result of shape in MC
598 ////  SetMaxPadRange(4);
599 ////  SetMaxRowRange(3);
600 //  SetMaxTimeRange(7);
601 //  SetValueToSize(3.5);//3.5
602 //  SetMaxPadRangeCm(2.5);
603 //  SetMaxRowRangeCm(3.5);
604
605   Int_t clusterCounter=0;
606   const Short_t nTotalSector=fParam->GetNSector();//number of sectors
607   for(Short_t iSec=0; iSec<nTotalSector; ++iSec){
608     
609     //vector of maxima for each sector
610     //std::vector<AliPadMax*> maximaInSector;
611     TObjArray *maximaInSector=new TObjArray();//to store AliPadMax*
612
613     //
614     //  looking for the maxima on the pad
615     //
616
617     const Short_t kNRows=fParam->GetNRow(iSec);//number of rows in sector
618     for(Short_t iRow=0; iRow<kNRows; ++iRow){
619       AliSimDigits *digrow;
620       if(fRawData){
621         digrow = (AliSimDigits*)fDigarr->GetRow(iSec,iRow);//real data
622       }else{
623         digrow = (AliSimDigits*)fDigarr->LoadRow(iSec,iRow);//MC
624       }
625       if(digrow){//if pointer exist
626         digrow->ExpandBuffer(); //decrunch
627         const Short_t kNPads = digrow->GetNCols();  // number of pads
628         const Short_t kNTime = digrow->GetNRows(); // number of timebins
629         for(Short_t iPad=0;iPad<kNPads;iPad++){
630           
631           Short_t timeBinMax=-1;//timebin of maximum 
632           Short_t valueMaximum=-1;//value of maximum in adc
633           Short_t increaseBegin=-1;//timebin when increase starts
634           Short_t sumAdc=0;//sum of adc on the pad in maximum surrounding
635           bool ifIncreaseBegin=true;//flag - check if increasing started
636           bool ifMaximum=false;//flag - check if it could be maximum
637
638           for(Short_t iTimeBin=0;iTimeBin<kNTime;iTimeBin++){
639             Short_t adc = digrow->GetDigitFast(iTimeBin,iPad);
640             if(adc<fMinAdc){//standard was 3
641               if(ifMaximum){
642                 if(iTimeBin-increaseBegin<fMinTimeBins){//at least 2 time bins
643                   timeBinMax=-1;
644                   valueMaximum=-1;
645                   increaseBegin=-1;
646                   sumAdc=0;
647                   ifIncreaseBegin=true;
648                   ifMaximum=false;
649                   continue;
650                 }
651                 //insert maximum, default values and set flags
652                 Double_t xCord,yCord;
653                 GetXY(iSec,iRow,iPad,xCord,yCord);
654                 AliPadMax *oneMaximum = new AliPadMax(AliTPCvtpr(valueMaximum,
655                                                                  timeBinMax,
656                                                                  iPad,
657                                                                  iRow,
658                                                                  xCord,
659                                                                  yCord,
660                                                                  timeBinMax),
661                                                       increaseBegin,
662                                                       iTimeBin-1,
663                                                       sumAdc);
664                 //maximaInSector.push_back(oneMaximum);
665                 maximaInSector->AddLast(oneMaximum);
666                 
667                 timeBinMax=-1;
668                 valueMaximum=-1;
669                 increaseBegin=-1;
670                 sumAdc=0;
671                 ifIncreaseBegin=true;
672                 ifMaximum=false;
673               }
674               continue;
675             }
676             
677             if(ifIncreaseBegin){
678               ifIncreaseBegin=false;
679               increaseBegin=iTimeBin;
680             }
681             
682             if(adc>valueMaximum){
683               timeBinMax=iTimeBin;
684               valueMaximum=adc;
685               ifMaximum=true;
686             }
687             sumAdc+=adc;
688             if(iTimeBin==kNTime-1 && ifMaximum && kNTime-increaseBegin>fMinTimeBins){//on the edge
689               //at least 3 timebins
690               //insert maximum, default values and set flags
691               Double_t xCord,yCord;
692               GetXY(iSec,iRow,iPad,xCord,yCord);
693               AliPadMax *oneMaximum = new AliPadMax(AliTPCvtpr(valueMaximum,
694                                                                timeBinMax,
695                                                                iPad,
696                                                                iRow,
697                                                                xCord,
698                                                                yCord,
699                                                                timeBinMax),
700                                                     increaseBegin,
701                                                     iTimeBin-1,
702                                                     sumAdc);
703               //maximaInSector.push_back(oneMaximum);
704               maximaInSector->AddLast(oneMaximum);
705                 
706               timeBinMax=-1;
707               valueMaximum=-1;
708               increaseBegin=-1;
709               sumAdc=0;
710               ifIncreaseBegin=true;
711               ifMaximum=false;
712               continue;
713             }
714             
715           }//end loop over timebins
716         }//end loop over pads
717 //      }else{
718 //      cout<<"Pointer does not exist!!"<<endl;
719       }//end if poiner exists
720     }//end loop over rows
721
722     //    cout<<"EF" <<maximaInSector->GetEntriesFast()<<" E "<<maximaInSector->GetEntries()<<" "<<maximaInSector->GetLast()<<endl;
723     maximaInSector->Compress();
724     // GetEntriesFast() - liczba wejsc w array of maxima
725     //cout<<"EF" <<maximaInSector->GetEntriesFast()<<" E"<<maximaInSector->GetEntries()<<endl;
726
727     //
728     // Making clusters
729     //
730
731     Short_t maxDig=0;
732     Short_t maxSumAdc=0;
733     Short_t maxTimeBin=0;
734     Short_t maxPad=0;
735     Short_t maxRow=0;
736     Double_t maxX=0;
737     Double_t maxY=0;
738     
739 //    for( std::vector<AliPadMax*>::iterator mp1  = maximaInSector.begin();
740 //       mp1 != maximaInSector.end(); ++mp1 ) {
741     for(Int_t it1 = 0; it1 < maximaInSector->GetEntriesFast(); ++it1 ) {
742
743       AliPadMax *mp1=(AliPadMax *)maximaInSector->At(it1);
744       AliTPCclusterKr *tmp=new AliTPCclusterKr();
745       
746       Short_t nUsedPads=1;
747       Int_t clusterValue=0;
748       clusterValue+=(mp1)->GetSum();
749       list<Short_t> nUsedRows;
750       nUsedRows.push_back((mp1)->GetRow());
751
752       maxDig      =(mp1)->GetAdc() ;
753       maxSumAdc   =(mp1)->GetSum() ;
754       maxTimeBin  =(mp1)->GetTime();
755       maxPad      =(mp1)->GetPad() ;
756       maxRow      =(mp1)->GetRow() ;
757       maxX        =(mp1)->GetX();
758       maxY        =(mp1)->GetY();
759
760       AliSimDigits *digrowTmp;
761       if(fRawData){
762         digrowTmp = (AliSimDigits*)fDigarr->GetRow(iSec,(mp1)->GetRow());
763       }else{
764         digrowTmp = (AliSimDigits*)fDigarr->LoadRow(iSec,(mp1)->GetRow());
765       }
766
767       digrowTmp->ExpandBuffer(); //decrunch
768
769       for(Short_t itb=(mp1)->GetBegin(); itb<((mp1)->GetEnd())+1; itb++){
770         Short_t adcTmp = digrowTmp->GetDigitFast(itb,(mp1)->GetPad());
771         AliTPCvtpr *vtpr=new AliTPCvtpr(adcTmp,itb,(mp1)->GetPad(),(mp1)->GetRow(),(mp1)->GetX(),(mp1)->GetY(),itb);
772         //tmp->fCluster.push_back(vtpr);
773         tmp->AddDigitToCluster(vtpr);
774       }
775       tmp->SetCenter();//set centr of the cluster
776       
777       //maximaInSector.erase(mp1);
778       //mp1--;
779       //cout<<maximaInSector->GetEntriesFast()<<" ";
780       maximaInSector->RemoveAt(it1);
781       maximaInSector->Compress();
782       it1--;
783       //     cout<<maximaInSector->GetEntriesFast()<<" "<<endl;
784
785 //      for( std::vector<AliPadMax*>::iterator mp2  = maximaInSector.begin();
786 //         mp2 != maximaInSector.end(); ++mp2 ) {
787       for(Int_t it2 = 0; it2 < maximaInSector->GetEntriesFast(); ++it2 ) {
788         AliPadMax *mp2=(AliPadMax *)maximaInSector->At(it2);
789
790 //      if(abs(maxRow - (*mp2)->GetRow()) < fMaxRowRange && //3
791 //         abs(maxPad - (*mp2)->GetPad()) < fMaxPadRange && //4
792           
793         if(TMath::Abs(tmp->GetCenterX() - (mp2)->GetX()) < fMaxPadRangeCm &&
794            TMath::Abs(tmp->GetCenterY() - (mp2)->GetY()) < fMaxRowRangeCm &&
795            TMath::Abs(tmp->GetCenterT() - (mp2)->GetT()) < fMaxTimeRange){//7
796           
797           clusterValue+=(mp2)->GetSum();
798
799           nUsedPads++;
800           nUsedRows.push_back((mp2)->GetRow());
801
802           AliSimDigits *digrowTmp1;
803           if(fRawData){
804             digrowTmp1 = (AliSimDigits*)fDigarr->GetRow(iSec,(mp2)->GetRow());
805           }else{
806             digrowTmp1 = (AliSimDigits*)fDigarr->LoadRow(iSec,(mp2)->GetRow());
807           }
808           digrowTmp1->ExpandBuffer(); //decrunch
809           
810           for(Short_t itb=(mp2)->GetBegin(); itb<(mp2)->GetEnd()+1; itb++){
811             Short_t adcTmp = digrowTmp1->GetDigitFast(itb,(mp2)->GetPad());
812             AliTPCvtpr *vtpr=new AliTPCvtpr(adcTmp,itb,(mp2)->GetPad(),(mp2)->GetRow(),(mp2)->GetX(),(mp2)->GetY(),itb);
813             //tmp->fCluster.push_back(vtpr);
814             tmp->AddDigitToCluster(vtpr);
815           }
816           
817           tmp->SetCenter();//set center of the cluster
818
819           //which one is bigger
820           if( (mp2)->GetAdc() > maxDig ){
821             maxDig      =(mp2)->GetAdc() ;
822             maxSumAdc   =(mp2)->GetSum() ;
823             maxTimeBin  =(mp2)->GetTime();
824             maxPad      =(mp2)->GetPad() ;
825             maxRow      =(mp2)->GetRow() ;
826             maxX        =(mp2)->GetX() ;
827             maxY        =(mp2)->GetY() ;
828           } else if ( (mp2)->GetAdc() == maxDig ){
829             if( (mp2)->GetSum() > maxSumAdc){
830               maxDig      =(mp2)->GetAdc() ;
831               maxSumAdc   =(mp2)->GetSum() ;
832               maxTimeBin  =(mp2)->GetTime();
833               maxPad      =(mp2)->GetPad() ;
834               maxRow      =(mp2)->GetRow() ;
835               maxX        =(mp2)->GetX() ;
836               maxY        =(mp2)->GetY() ;
837             }
838           }
839           maximaInSector->RemoveAt(it2);
840           maximaInSector->Compress();
841           it2--;
842
843           //maximaInSector.erase(mp2);
844           //mp2--;
845         }
846       }//inside loop
847
848       tmp->SetSize();
849       //through out ADC=6,7 on 1 pad, 2 tb and ADC=12 on 2 pads,2 tb
850       //if(nUsedPads==1 && clusterValue/tmp->fCluster.size()<3.6)continue;
851       //if(nUsedPads==2 && clusterValue/tmp->fCluster.size()<3.1)continue;
852
853       //through out clusters on the edge and noise
854       //if(clusterValue/tmp->fCluster.size()<fValueToSize)continue;
855       if(clusterValue/(tmp->GetSize())<fValueToSize)continue;
856
857       tmp->SetADCcluster(clusterValue);
858       tmp->SetNPads(nUsedPads);
859       tmp->SetMax(AliTPCvtpr(maxDig,maxTimeBin,maxPad,maxRow,maxX,maxY,maxTimeBin));
860       tmp->SetSec(iSec);
861       tmp->SetSize();
862
863       nUsedRows.sort();
864       nUsedRows.unique();
865       tmp->SetNRows(nUsedRows.size());
866       tmp->SetCenter();
867
868       clusterCounter++;
869       
870
871       //save each cluster into file
872
873       AliTPCClustersRow *clrow= new AliTPCClustersRow();
874       fRowCl=clrow;
875       clrow->SetClass("AliTPCclusterKr");
876       clrow->SetArray(1);
877       fOutput->GetBranch("Segment")->SetAddress(&clrow);
878
879       Int_t tmpCluster=0;
880       TClonesArray * arr = fRowCl->GetArray();
881       AliTPCclusterKr * cl = new ((*arr)[tmpCluster]) AliTPCclusterKr(*tmp);
882       
883       fOutput->Fill(); 
884       delete clrow;
885       //end of save each cluster into file adc.root
886     }//outer loop
887
888   }//end sector for
889   cout<<"Number of clusters in event: "<<clusterCounter<<endl;
890
891   return 0;
892 }
893
894 ////____________________________________________________________________________
895
896
897 void AliTPCclustererKr::GetXY(Short_t sec,Short_t row,Short_t pad,Double_t& xGlob,Double_t& yGlob){
898   //
899   //gives global XY coordinate of the pad
900   //
901
902   Double_t yLocal = fParam->GetPadRowRadii(sec,row);//radius of row in sector in cm
903
904   Short_t padmax = fParam->GetNPads(sec,row);//number of pads in a given row
905   Float_t padXSize;
906   if(sec<fParam->GetNInnerSector())padXSize=0.4;
907   else padXSize=0.6;
908   Double_t xLocal=(pad+0.5-padmax/2.)*padXSize;//x-value of the center of pad
909
910   Float_t sin,cos;
911   fParam->AdjustCosSin((Int_t)sec,cos,sin);//return sinus and cosinus of the sector
912
913   Double_t xGlob1 =  xLocal * cos + yLocal * sin;
914   Double_t yGlob1 = -xLocal * sin + yLocal * cos;
915
916
917   Double_t rot=0;
918   rot=TMath::Pi()/2.;
919
920   xGlob =  xGlob1 * TMath::Cos(rot) + yGlob1 * TMath::Sin(rot);
921   yGlob = -xGlob1 * TMath::Sin(rot) + yGlob1 * TMath::Cos(rot);
922
923    yGlob=-1*yGlob;
924    if(sec<18||(sec>=36&&sec<54)) xGlob =-1*xGlob;
925
926
927   return;
928 }