Classes for handling Kr cluster finder
[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 cluster class
20 //
21 // Origin: Adam Matyja, INP PAN, adam.matyja@ifj.edu.pl
22 //-----------------------------------------------------------------
23
24 #include "AliTPCclustererKr.h"
25 #include "AliTPCclusterKr.h"
26 #include <vector>
27 #include "TObject.h"
28 #include "AliPadMax.h"
29 #include "AliSimDigits.h"
30 #include "AliTPCv4.h"
31 #include "AliTPCParam.h"
32 #include "AliTPCDigitsArray.h"
33 #include "AliTPCvtpr.h"
34 #include "AliTPCClustersRow.h"
35 #include "TTree.h"
36
37 //used in raw data finder
38 #include "AliTPCROC.h"
39 #include "AliTPCCalPad.h"
40 #include "AliTPCAltroMapping.h"
41 #include "AliTPCcalibDB.h"
42 #include "AliTPCRawStream.h"
43 #include "AliTPCRecoParam.h"
44 #include "AliTPCReconstructor.h"
45 #include "AliRawReader.h"
46 #include "AliTPCCalROC.h"
47
48 ClassImp(AliTPCclustererKr)
49
50
51 AliTPCclustererKr::AliTPCclustererKr()
52   :TObject(),
53   fRawData(kFALSE),
54   fRowCl(0),
55   fInput(0),
56   fOutput(0),
57   fParam(0),
58   fDigarr(0),
59   fRecoParam(0),
60   fIsOldRCUFormat(kFALSE)
61 {
62 //
63 // default constructor
64 //
65 }
66
67 AliTPCclustererKr::AliTPCclustererKr(const AliTPCclustererKr &param)
68   :TObject(),
69   fRawData(kFALSE),
70   fRowCl(0),
71   fInput(0),
72   fOutput(0),
73   fParam(0),
74   fDigarr(0),
75   fRecoParam(0),
76   fIsOldRCUFormat(kFALSE)
77 {
78 //
79 // copy constructor
80 //
81   fParam = param.fParam;
82   fRecoParam=param.fRecoParam;
83   fIsOldRCUFormat=param.fIsOldRCUFormat;
84   fRawData=param.fRawData;
85   fRowCl =param.fRowCl ;
86   fInput =param.fInput ;
87   fOutput=param.fOutput;
88   fDigarr=param.fDigarr;
89
90
91 AliTPCclustererKr & AliTPCclustererKr::operator = (const AliTPCclustererKr & param)
92 {
93   fParam = param.fParam;
94   fRecoParam=param.fRecoParam;
95   fIsOldRCUFormat=param.fIsOldRCUFormat;
96   fRawData=param.fRawData;
97   fRowCl =param.fRowCl ;
98   fInput =param.fInput ;
99   fOutput=param.fOutput;
100   fDigarr=param.fDigarr;
101   return (*this);
102 }
103
104 AliTPCclustererKr::~AliTPCclustererKr()
105 {
106   //
107   // destructor
108   //
109 }
110
111 void AliTPCclustererKr::SetRecoParam(AliTPCRecoParam *recoParam)
112 {
113   if (recoParam) {
114     fRecoParam = recoParam;
115   }else{
116     //set default parameters if not specified
117     fRecoParam = AliTPCReconstructor::GetRecoParam();
118     if (!fRecoParam)  fRecoParam = AliTPCRecoParam::GetLowFluxParam();
119   }
120   return;
121 }
122
123
124 ////____________________________________________________________________________
125 ////       I/O
126 void AliTPCclustererKr::SetInput(TTree * tree)
127 {
128   //
129   // set input tree with digits
130   //
131   fInput = tree;  
132   if  (!fInput->GetBranch("Segment")){
133     cerr<<"AliTPCclusterKr::findClusterKr(): no proper input tree !\n";
134     fInput=0;
135     return;
136   }
137 }
138
139 void AliTPCclustererKr::SetOutput(TTree * tree) 
140 {
141   //
142   // set output
143   //
144   fOutput= tree;  
145   AliTPCClustersRow clrow;
146   AliTPCClustersRow *pclrow=&clrow;  
147   clrow.SetClass("AliTPCclusterKr");
148   clrow.SetArray(1); // to make Clones array
149   fOutput->Branch("Segment","AliTPCClustersRow",&pclrow,32000,200);    
150 }
151
152 ////____________________________________________________________________________
153 //// with new I/O
154 Int_t AliTPCclustererKr::finderIO()
155 {
156   // Krypton cluster finder for simulated events from MC
157
158   if (!fInput) { 
159     Error("Digits2Clusters", "input tree not initialised");
160     return 10;
161   }
162   
163   if (!fOutput) {
164     Error("Digits2Clusters", "output tree not initialised");
165     return 11;
166   }
167
168   findClusterKrIO();
169   return 0;
170 }
171
172 Int_t AliTPCclustererKr::finderIO(AliRawReader* rawReader)
173 {
174   // Krypton cluster finder for the TPC raw data
175   //
176   // fParam must be defined before
177
178   // consider noiceROC or not
179
180   if(rawReader)fRawData=kTRUE; //set flag to data
181
182   if (!fOutput) {
183     Error("Digits2Clusters", "output tree not initialised");
184     return 11;
185   }
186
187   fParam->SetMaxTBin(fRecoParam->GetLastBin());//set number of timebins from reco -> param
188   //   used later for memory allocation
189
190   Bool_t isAltro=kFALSE;
191
192   AliTPCROC * roc = AliTPCROC::Instance();
193   //AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
194   AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping();
195   //
196   AliTPCRawStream input(rawReader,(AliAltroMapping**)mapping);
197
198   const Int_t kNIS = fParam->GetNInnerSector();//number of inner sectors
199   const Int_t kNOS = fParam->GetNOuterSector();//number of outer sectors
200   const Int_t kNS = kNIS + kNOS;//all sectors
201   Int_t zeroSup = fParam->GetZeroSup();//zero suppression parameter
202
203   //crate TPC view
204   AliTPCDigitsArray *digarr=new AliTPCDigitsArray(kFALSE);//data not sim
205   digarr->Setup(fParam);//as usually parameters
206
207   //
208   // Loop over sectors
209   //
210   for(Int_t sec = 0; sec < kNS; sec++) {
211     //AliTPCCalROC * noiseROC   = noiseTPC->GetCalROC(sec);  // noise per given sector
212  
213     Int_t nRows = 0; //number of rows in sector
214     Int_t nDDLs = 0; //number of DDLs
215     Int_t indexDDL = 0; //DDL index
216     if (sec < kNIS) {
217       nRows = fParam->GetNRowLow();
218       nDDLs = 2;
219       indexDDL = sec * 2;
220     }else {
221       nRows = fParam->GetNRowUp();
222       nDDLs = 4;
223       indexDDL = (sec-kNIS) * 4 + kNIS * 2;
224     }
225
226     //
227     // Load the raw data for corresponding DDLs
228     //
229     rawReader->Reset();
230     input.SetOldRCUFormat(fIsOldRCUFormat);
231     rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1);
232
233     if(input.Next()) {
234       isAltro=kTRUE;
235       // Allocate memory for rows in sector (pads(depends on row) x timebins)
236       for(Int_t row = 0; row < nRows; row++) {
237         digarr->CreateRow(sec,row);
238       }//end loop over rows
239     }
240     rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1);
241
242     //
243     // Begin loop over altro data
244     //
245     while (input.Next()) {
246
247       //check sector consistency
248       if (input.GetSector() != sec)
249         AliFatal(Form("Sector index mismatch ! Expected (%d), but got (%d) !",sec,input.GetSector()));
250       
251       //check row consistency
252       Short_t iRow = input.GetRow();
253       if (iRow < 0 || iRow >= nRows){
254         AliError(Form("Pad-row index (%d) outside the range (%d -> %d) !",
255                       iRow, 0, nRows -1));
256         continue;
257       }
258
259       //check pad consistency
260       Short_t iPad = input.GetPad();
261       if (iPad < 0 || iPad >= (Short_t)(roc->GetNPads(sec,iRow))) {
262         AliError(Form("Pad index (%d) outside the range (%d -> %d) !",
263                       iPad, 0, roc->GetNPads(sec,iRow) ));
264         continue;
265       }
266
267       //check time consistency
268       Short_t iTimeBin = input.GetTime();
269       if ( iTimeBin < fRecoParam->GetFirstBin() || iTimeBin >= fRecoParam->GetLastBin()){
270         //cout<<iTimeBin<<endl;
271         continue;
272         AliFatal(Form("Timebin index (%d) outside the range (%d -> %d) !",
273                       iTimeBin, 0, fRecoParam->GetLastBin() -1));
274       }
275
276       //signal
277       Int_t signal = input.GetSignal();
278       if (signal <= zeroSup) {
279         //continue;
280         digarr->GetRow(sec,iRow)->SetDigitFast(0,iTimeBin,iPad);
281       }
282       (//(AliSimDigits*)
283        digarr->GetRow(sec,iRow))->SetDigitFast(signal,iTimeBin,iPad);
284     }//end of loop over altro data
285   }//end of loop over sectors
286   
287   SetDigArr(digarr);
288   if(isAltro) findClusterKrIO();
289   delete digarr;
290
291   return 0;
292 }
293
294 ////____________________________________________________________________________
295 Int_t AliTPCclustererKr::findClusterKrIO()
296 {
297   //fParam and  fDigarr must be set to run this method
298
299   Int_t cluster_counter=0;
300   const Short_t NTotalSector=fParam->GetNSector();//number of sectors
301   for(Short_t sec=0; sec<NTotalSector; sec++){
302     
303     //vector of maxima for each sector
304     std::vector<AliPadMax*> maxima_in_sector;
305     
306     const Short_t Nrows=fParam->GetNRow(sec);//number of rows in sector
307     for(Short_t row=0; row<Nrows; row++){
308       AliSimDigits *digrow;
309       if(fRawData){
310         digrow = (AliSimDigits*)fDigarr->GetRow(sec,row);//real data
311       }else{
312         digrow = (AliSimDigits*)fDigarr->LoadRow(sec,row);//MC
313       }
314       if(digrow){//if pointer exist
315         digrow->ExpandBuffer(); //decrunch
316         const Short_t npads = digrow->GetNCols();  // number of pads
317         const Short_t ntime = digrow->GetNRows(); // number of timebins
318         for(Short_t np=0;np<npads;np++){
319           
320           Short_t tb_max=-1;//timebin of maximum 
321           Short_t value_maximum=-1;//value of maximum in adc
322           Short_t increase_begin=-1;//timebin when increase starts
323           Short_t sum_adc=0;//sum of adc on the pad in maximum surrounding
324           bool if_increase_begin=true;//flag - check if increasing start
325           bool if_maximum=false;//flag - check if it could be maximum
326
327           for(Short_t nt=0;nt<ntime;nt++){
328             Short_t adc = digrow->GetDigitFast(nt,np);
329             if(adc<3){
330               if(if_maximum){
331                 if(nt-1-increase_begin<1){//at least 2 time bins
332                   tb_max=-1;
333                   value_maximum=-1;
334                   increase_begin=-1;
335                   sum_adc=0;
336                   if_increase_begin=true;
337                   if_maximum=false;
338                   continue;
339                 }
340                 //insert maximum, default values and set flags
341                 AliPadMax *one_maximum = new AliPadMax(AliTPCvtpr(value_maximum,
342                                                                   tb_max,
343                                                                   np,
344                                                                   row),
345                                                        increase_begin,
346                                                        nt-1,
347                                                        sum_adc);
348                 maxima_in_sector.push_back(one_maximum);
349                 
350                 tb_max=-1;
351                 value_maximum=-1;
352                 increase_begin=-1;
353                 sum_adc=0;
354                 if_increase_begin=true;
355                 if_maximum=false;
356               }
357               continue;
358             }
359             
360             if(if_increase_begin){
361               if_increase_begin=false;
362               increase_begin=nt;
363             }
364             
365             if(adc>value_maximum){
366               tb_max=nt;
367               value_maximum=adc;
368               if_maximum=true;
369             }
370             sum_adc+=adc;
371             if(nt==ntime-1 && if_maximum && ntime-1-increase_begin>2){//on the edge
372               //insert maximum, default values and set flags
373               AliPadMax *one_maximum = new AliPadMax(AliTPCvtpr(value_maximum,
374                                                                 tb_max,
375                                                                 np,
376                                                                 row),
377                                                      increase_begin,
378                                                      nt-1,
379                                                      sum_adc);
380               maxima_in_sector.push_back(one_maximum);
381               
382               tb_max=-1;
383               value_maximum=-1;
384               increase_begin=-1;
385               sum_adc=0;
386               if_increase_begin=true;
387               if_maximum=false;
388               continue;
389             }
390             
391           }//end loop over timebins
392         }//end loop over pads
393 //      }else{
394 //      cout<<"Pointer does not exist!!"<<endl;
395       }//end if poiner exists
396     }//end loop over rows
397     
398
399     Short_t max_dig=0;
400     Short_t max_sum_adc=0;
401     Short_t max_nt=0;
402     Short_t max_np=0;
403     Short_t max_row=0;
404     
405     for( std::vector<AliPadMax*>::iterator mp1  = maxima_in_sector.begin();
406          mp1 != maxima_in_sector.end(); ++mp1 ) {
407       
408       AliTPCclusterKr *tmp=new AliTPCclusterKr();
409       
410       Short_t n_used_pads=1;
411       Short_t cluster_value=0;
412       cluster_value+=(*mp1)->GetSum();
413       
414       max_dig      =(*mp1)->GetAdc() ;
415       max_sum_adc  =(*mp1)->GetSum() ;
416       max_nt       =(*mp1)->GetTime();
417       max_np       =(*mp1)->GetPad() ;
418       max_row      =(*mp1)->GetRow() ;
419       
420       AliSimDigits *digrow_tmp;
421       if(fRawData){
422         digrow_tmp = (AliSimDigits*)fDigarr->GetRow(sec,(*mp1)->GetRow());
423       }else{
424         digrow_tmp = (AliSimDigits*)fDigarr->LoadRow(sec,(*mp1)->GetRow());
425       }
426
427       digrow_tmp->ExpandBuffer(); //decrunch
428       for(Short_t tb=(*mp1)->GetBegin(); tb<((*mp1)->GetEnd())+1; tb++){
429         Short_t adc_tmp = digrow_tmp->GetDigitFast(tb,(*mp1)->GetPad());
430         AliTPCvtpr *vtpr=new AliTPCvtpr(adc_tmp,tb,(*mp1)->GetPad(),(*mp1)->GetRow());
431         tmp->fCluster.push_back(vtpr);
432       }
433       
434       maxima_in_sector.erase(mp1);
435       mp1--;
436       
437       for( std::vector<AliPadMax*>::iterator mp2  = maxima_in_sector.begin();
438            mp2 != maxima_in_sector.end(); ++mp2 ) {
439         
440         Short_t two_pad_max=4;
441         Short_t two_row_max=3;
442         
443         if(abs(max_row - (*mp2)->GetRow())<two_row_max && 
444            abs(max_np - (*mp2)->GetPad())<two_pad_max  &&
445            abs(max_nt - (*mp2)->GetTime())<7){
446           
447           cluster_value+=(*mp2)->GetSum();
448           n_used_pads++;
449           
450           AliSimDigits *digrow_tmp1;
451           if(fRawData){
452             digrow_tmp1 = (AliSimDigits*)fDigarr->GetRow(sec,(*mp2)->GetRow());
453           }else{
454             digrow_tmp1 = (AliSimDigits*)fDigarr->LoadRow(sec,(*mp2)->GetRow());
455           }
456           digrow_tmp1->ExpandBuffer(); //decrunch
457           
458           for(Short_t tb=(*mp2)->GetBegin(); tb<(*mp2)->GetEnd()+1; tb++){
459             Short_t adc_tmp = digrow_tmp1->GetDigitFast(tb,(*mp2)->GetPad());
460             AliTPCvtpr *vtpr=new AliTPCvtpr(adc_tmp,tb,(*mp2)->GetPad(),(*mp2)->GetRow());
461               tmp->fCluster.push_back(vtpr);
462           }
463           
464           //which one is bigger
465           if( (*mp2)->GetAdc() > max_dig ){
466             max_dig      =(*mp2)->GetAdc() ;
467             max_sum_adc  =(*mp2)->GetSum() ;
468             max_nt       =(*mp2)->GetTime();
469             max_np       =(*mp2)->GetPad() ;
470             max_row      =(*mp2)->GetRow() ;
471           } else if ( (*mp2)->GetAdc() == max_dig ){
472             if( (*mp2)->GetSum() > max_sum_adc){
473               max_dig      =(*mp2)->GetAdc() ;
474               max_sum_adc  =(*mp2)->GetSum() ;
475               max_nt       =(*mp2)->GetTime();
476               max_np       =(*mp2)->GetPad() ;
477               max_row      =(*mp2)->GetRow() ;
478             }
479           }
480           maxima_in_sector.erase(mp2);
481           mp2--;
482         }
483       }//inside loop
484       
485       //through out ADC=6,7 on 1 pad, 2 tb and ADC=12 on 2 pads,2 tb
486       if(n_used_pads==1 && cluster_value/tmp->fCluster.size()<3.6)continue;
487       if(n_used_pads==2 && cluster_value/tmp->fCluster.size()<3.1)continue;
488       
489       tmp->SetADCcluster(cluster_value);
490       tmp->SetNpads(n_used_pads);
491       tmp->SetMax(AliTPCvtpr(max_dig,max_nt,max_np,max_row));
492       tmp->SetSec(sec);
493       tmp->SetSize();
494       
495       cluster_counter++;
496       
497       //save each cluster into file
498
499       AliTPCClustersRow *clrow= new AliTPCClustersRow();
500       fRowCl=clrow;
501       clrow->SetClass("AliTPCclusterKr");
502       clrow->SetArray(1);
503       fOutput->GetBranch("Segment")->SetAddress(&clrow);
504
505       Int_t tmp_cluster=0;
506       TClonesArray * arr = fRowCl->GetArray();
507       AliTPCclusterKr * cl = new ((*arr)[tmp_cluster]) AliTPCclusterKr(*tmp);
508       
509       fOutput->Fill(); 
510       delete clrow;
511
512       //end of save each cluster into file adc.root
513     }//outer loop
514
515   }//end sector for
516   cout<<"Number of clusters in event: "<<cluster_counter<<endl;
517
518   return 0;
519 }
520
521 ////____________________________________________________________________________