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