Modifications to allow reclusterization during analysis
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALClusterizer.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$ */
17
18 //_________________________________________________________________________
19 //  Base class for the clusterization algorithm (pure abstract)
20 //*--
21 //*-- Author: Yves Schutz  SUBATECH 
22 // 
23 //   Clusterization mother class. Contains common methods/data members of different 
24 //   clusterizers. GCB 2010
25 //////////////////////////////////////////////////////////////////////////////
26
27 // --- ROOT system ---
28 #include "TClonesArray.h"
29 #include "TTree.h"
30 #include <TFile.h> 
31 class TFolder;
32 #include <TMath.h> 
33 #include <TMinuit.h>
34 #include <TTree.h> 
35 class TSystem; 
36 #include <TBenchmark.h>
37 #include <TBrowser.h>
38 #include <TROOT.h>
39
40 // --- Standard library ---
41 #include <cassert>
42
43 // --- AliRoot header files ---
44 #include "AliEMCALClusterizer.h"
45 #include "AliEMCALReconstructor.h"
46 #include "AliRunLoader.h"
47 #include "AliRun.h"
48 #include "AliLog.h"
49 #include "AliEMCAL.h"
50 #include "AliEMCALRecPoint.h"
51 #include "AliEMCALRecParam.h"
52 #include "AliEMCALGeometry.h"
53 #include "AliEMCALRecParam.h"
54 #include "AliCDBManager.h"
55 #include "AliCaloCalibPedestal.h"
56 #include "AliEMCALCalibData.h"
57 class AliCDBStorage;
58 #include "AliCDBEntry.h"
59
60 ClassImp(AliEMCALClusterizer)
61
62 Bool_t AliEMCALClusterizer::fgkIsInputCalibrated = kFALSE;
63
64
65 //____________________________________________________________________________
66 AliEMCALClusterizer::AliEMCALClusterizer():
67   fDigitsArr(NULL),
68   fTreeR(NULL),
69   fRecPoints(NULL),
70   fGeom(NULL),
71   fCalibData(NULL), 
72   fCaloPed(NULL),
73   fADCchannelECA(0.),fADCpedestalECA(0.),
74   fTimeMin(-1.),fTimeMax(1.),fTimeCut(1.),
75   fDefaultInit(kFALSE),fToUnfold(kFALSE),
76   fNumberOfECAClusters(0), fECAClusteringThreshold(0.),
77   fECALocMaxCut(0.),fECAW0(0.),fMinECut(0.),
78   fClusterUnfolding(NULL)
79 {
80   // ctor
81   
82   Init();
83   
84 }
85
86 //____________________________________________________________________________
87 AliEMCALClusterizer::AliEMCALClusterizer(AliEMCALGeometry* geometry): 
88   fDigitsArr(NULL),
89   fTreeR(NULL),
90   fRecPoints(NULL),
91   fGeom(geometry),
92   fCalibData(NULL), 
93   fCaloPed(NULL),
94   fADCchannelECA(0.),fADCpedestalECA(0.),
95   fTimeMin(-1.),fTimeMax(1.),fTimeCut(1.),
96   fDefaultInit(kFALSE),fToUnfold(kFALSE),
97   fNumberOfECAClusters(0), fECAClusteringThreshold(0.),
98   fECALocMaxCut(0.),fECAW0(0.),fMinECut(0.),
99   fClusterUnfolding(NULL)
100 {
101   // ctor with the indication of the file where header Tree and digits Tree are stored
102   // use this contructor to avoid usage of Init() which uses runloader
103   // change needed by HLT - MP
104   
105   // Note for the future: the use on runloader should be avoided or optional at least
106   // another way is to make Init virtual and protected at least such that the deriving classes can overload
107   // Init() ;
108   //
109   
110   if (!fGeom)
111   {
112     AliFatal("Geometry not initialized.");
113   }
114   Int_t i=0;
115   for (i = 0; i < 8; i++)
116     fSSPars[i] = 0.;
117   for (i = 0; i < 3; i++) {
118     fPar5[i] = 0.;
119     fPar6[i] = 0.;
120   }
121 }
122
123 //____________________________________________________________________________
124 AliEMCALClusterizer::AliEMCALClusterizer(AliEMCALGeometry* geometry, AliEMCALCalibData * calib, AliCaloCalibPedestal * caloped): 
125   fDigitsArr(NULL),
126   fTreeR(NULL),
127   fRecPoints(NULL),
128   fGeom(geometry),
129   fCalibData(calib),
130   fCaloPed(caloped),
131   fADCchannelECA(0.),fADCpedestalECA(0.),
132   fTimeMin(-1.),fTimeMax(1.),fTimeCut(1.),
133   fDefaultInit(kFALSE),fToUnfold(kFALSE),
134   fNumberOfECAClusters(0), fECAClusteringThreshold(0.),
135   fECALocMaxCut(0.),fECAW0(0.),fMinECut(0.),
136   fClusterUnfolding(NULL)
137 {
138   // ctor, geometry and calibration are initialized elsewhere.
139   
140   if (!fGeom)
141     AliFatal("Geometry not initialized.");
142   
143   Int_t i=0;
144   for (i = 0; i < 8; i++)
145     fSSPars[i] = 0.;
146   for (i = 0; i < 3; i++) {
147     fPar5[i] = 0.;
148     fPar6[i] = 0.;
149   }
150   
151 }
152
153
154 //____________________________________________________________________________
155 AliEMCALClusterizer::~AliEMCALClusterizer()
156 {
157   // dtor
158   //Already deleted in AliEMCALReconstructor.
159
160 //   if(fGeom)      delete fGeom;
161 //   if(fCalibData) delete fCalibData;
162 //   if(fCaloPed)   delete fCaloPed;
163
164 //   if (fDigitsArr) {
165 //     fDigitsArr->Clear("C");
166 //     delete fDigitsArr;
167 //   }
168 //   if (fRecPoints) {
169 //     fRecPoints->Delete();
170 //     delete fRecPoints;
171 //  }
172   
173   if(fClusterUnfolding) delete fClusterUnfolding;
174 }
175
176 //____________________________________________________________________________
177 Float_t  AliEMCALClusterizer::Calibrate(const Float_t amp, const Float_t time, const Int_t absId) 
178 {
179   
180   // Convert digitized amplitude into energy.
181   // Calibration parameters are taken from calibration data base for raw data,
182   // or from digitizer parameters for simulated data.
183   if(fCalibData){
184     
185     if (fGeom==0)
186       AliFatal("Did not get geometry from EMCALLoader") ;
187     
188     Int_t iSupMod = -1;
189     Int_t nModule = -1;
190     Int_t nIphi   = -1;
191     Int_t nIeta   = -1;
192     Int_t iphi    = -1;
193     Int_t ieta    = -1;
194     
195     Bool_t bCell = fGeom->GetCellIndex(absId, iSupMod, nModule, nIphi, nIeta) ;
196     if(!bCell) {
197       fGeom->PrintGeometry();
198       Error("Calibrate()"," Wrong cell id number : %i", absId);
199       assert(0);
200     }
201     
202     fGeom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
203           
204     // Check if channel is bad (dead or hot), in this case return 0.    
205     // Gustavo: 15-12-09 In case of RAW data this selection is already done, but not in simulation.
206     // for the moment keep it here but remember to do the selection at the sdigitizer level 
207     // and remove it from here
208     Int_t channelStatus = (Int_t)(fCaloPed->GetDeadMap(iSupMod))->GetBinContent(ieta,iphi);
209     if(channelStatus == AliCaloCalibPedestal::kHot || channelStatus == AliCaloCalibPedestal::kDead) {
210                   AliDebug(2,Form("Tower from SM %d, ieta %d, iphi %d is BAD : status %d !!!",iSupMod,ieta,iphi, channelStatus));
211                   return 0;
212     }
213     //Check if time is too large or too small, indication of a noisy channel, remove in this case
214     if(time > fTimeMax || time < fTimeMin) return 0;
215     
216     if (AliEMCALClusterizer::fgkIsInputCalibrated == kTRUE)
217     {
218       AliDebug(10, Form("Input already calibrated!"));
219       return amp;
220     }
221           
222     fADCchannelECA  = fCalibData->GetADCchannel (iSupMod,ieta,iphi);
223     fADCpedestalECA = fCalibData->GetADCpedestal(iSupMod,ieta,iphi);
224     
225     
226     return -fADCpedestalECA + amp * fADCchannelECA ;        
227     
228   }
229   else{ //Return energy with default parameters if calibration is not available
230     
231     if (AliEMCALClusterizer::fgkIsInputCalibrated == kTRUE)
232     {
233       AliDebug(10, Form("Input already calibrated!"));
234       return amp;
235     }    
236     
237     return -fADCpedestalECA + amp * fADCchannelECA ; 
238   }
239   
240 }
241
242 //____________________________________________________________________________
243 void AliEMCALClusterizer::GetCalibrationParameters() 
244 {
245   // Set calibration parameters:
246   // if calibration database exists, they are read from database,
247   // otherwise, they are taken from digitizer.
248   //
249   // It is a user responsilibity to open CDB before reconstruction, 
250   // for example: 
251   // AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB");
252   
253   //Check if calibration is stored in data base
254   
255   if(!fCalibData)
256   {
257     AliCDBEntry *entry = (AliCDBEntry*) 
258     AliCDBManager::Instance()->Get("EMCAL/Calib/Data");
259     if (entry) fCalibData =  (AliEMCALCalibData*) entry->GetObject();
260   }
261   
262   if(!fCalibData)
263     AliFatal("Calibration parameters not found in CDB!");
264   
265 }
266
267 //____________________________________________________________________________
268 void AliEMCALClusterizer::GetCaloCalibPedestal() 
269 {
270   // Set calibration parameters:
271   // if calibration database exists, they are read from database,
272   // otherwise, they are taken from digitizer.
273   //
274   // It is a user responsilibity to open CDB before reconstruction, 
275   // for example: 
276   // AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB");
277   
278   //Check if calibration is stored in data base
279   
280   if(!fCaloPed)
281     {
282       AliCDBEntry *entry = (AliCDBEntry*) 
283         AliCDBManager::Instance()->Get("EMCAL/Calib/Pedestals");
284       if (entry) fCaloPed =  (AliCaloCalibPedestal*) entry->GetObject();
285     }
286   
287   if(!fCaloPed)
288     AliFatal("Pedestal info not found in CDB!");
289   
290 }
291
292 //____________________________________________________________________________
293 void AliEMCALClusterizer::Init()
294 {
295   // Make all memory allocations which can not be done in default constructor.
296   // Attach the Clusterizer task to the list of EMCAL tasks
297   
298   AliRunLoader *rl = AliRunLoader::Instance();
299   if (rl->GetAliRun()){
300     AliEMCAL* emcal = dynamic_cast<AliEMCAL*>(rl->GetAliRun()->GetDetector("EMCAL"));
301     if(emcal)fGeom = emcal->GetGeometry();
302   }
303   
304   if(!fGeom){ 
305     fGeom =  AliEMCALGeometry::GetInstance(AliEMCALGeometry::GetDefaultGeometryName());
306   }
307   
308   AliDebug(1,Form("geom %p",fGeom));
309   
310   if(!gMinuit) 
311     gMinuit = new TMinuit(100) ;
312   
313   Int_t i=0;
314   for (i = 0; i < 8; i++)
315     fSSPars[i] = 0.;
316   for (i = 0; i < 3; i++) {
317     fPar5[i] = 0.;
318     fPar6[i] = 0.;
319   }
320   
321 }
322
323 //____________________________________________________________________________
324 void AliEMCALClusterizer::InitParameters()
325
326   // Initializes the parameters for the Clusterizer
327   fNumberOfECAClusters = 0 ;
328   fCalibData           = 0 ;
329   fCaloPed             = 0 ;
330         
331   const AliEMCALRecParam* recParam = AliEMCALReconstructor::GetRecParam();
332   if(!recParam) {
333     AliFatal("Reconstruction parameters for EMCAL not set!");
334   } else {
335     fECAClusteringThreshold = recParam->GetClusteringThreshold();
336     fECAW0                  = recParam->GetW0();
337     fMinECut                = recParam->GetMinECut();    
338     fToUnfold               = recParam->GetUnfold();
339     fECALocMaxCut           = recParam->GetLocMaxCut();
340     fTimeCut                = recParam->GetTimeCut();
341     fTimeMin                = recParam->GetTimeMin();
342     fTimeMax                = recParam->GetTimeMax();
343     
344     AliDebug(1,Form("Reconstruction parameters: fECAClusteringThreshold=%.3f GeV, fECAW=%.3f, fMinECut=%.3f GeV, fToUnfold=%d, fECALocMaxCut=%.3f GeV, fTimeCut=%e s,fTimeMin=%e s,fTimeMax=%e s",
345                     fECAClusteringThreshold,fECAW0,fMinECut,fToUnfold,fECALocMaxCut,fTimeCut, fTimeMin, fTimeMax));
346     
347     if(fToUnfold){
348       
349       Int_t i=0;
350       for (i = 0; i < 8; i++) {
351         fSSPars[i] = recParam->GetSSPars(i);
352       }//end of loop over parameters
353       for (i = 0; i < 3; i++) {
354         fPar5[i] = recParam->GetPar5(i);
355         fPar6[i] = recParam->GetPar6(i);
356       }//end of loop over parameters
357       
358       InitClusterUnfolding();
359       
360       for (i = 0; i < 8; i++) {
361         AliDebug(1,Form("unfolding shower shape parameters: fSSPars=%f \n",fSSPars[i]));
362       }
363       for (i = 0; i < 3; i++) {
364         AliDebug(1,Form("unfolding parameter 5: fPar5=%f \n",fPar5[i]));
365         AliDebug(1,Form("unfolding parameter 6: fPar6=%f \n",fPar6[i]));
366       }
367       
368     }// to unfold
369   }// recparam not null
370   
371 }
372
373
374 //____________________________________________________________________________
375 void AliEMCALClusterizer::Print(Option_t * /*option*/)const
376 {
377   // Print clusterizer parameters
378   
379   TString message("\n") ; 
380   
381   if( strcmp(GetName(), "") !=0 ){
382     
383     // Print parameters
384     
385     TString taskName(Version()) ;
386     
387     printf("--------------- "); 
388     printf("%s",taskName.Data()) ; 
389     printf(" "); 
390     printf("Clusterizing digits: "); 
391     printf("\n                       ECA Local Maximum cut    = %f", fECALocMaxCut); 
392     printf("\n                       ECA Logarithmic weight   = %f", fECAW0); 
393     if(fToUnfold){
394       printf("\nUnfolding on\n");
395       printf("Unfolding parameters: fSSpars: \n");
396       Int_t i=0;
397       for (i = 0; i < 8; i++) {
398         printf("fSSPars[%d] = %f \n", i, fSSPars[i]);
399       }
400       printf("Unfolding parameter 5 and 6: fPar5 and fPar6: \n");
401       for (i = 0; i < 3; i++) {
402         printf("fPar5[%d] = %f \n", i, fPar5[i]);
403         printf("fPar6[%d] = %f \n", i, fPar6[i]);
404       }
405     }
406     else
407       printf("\nUnfolding off\n");
408     
409     printf("------------------------------------------------------------------"); 
410   }
411   else
412     printf("AliEMCALClusterizer not initialized ") ;
413 }
414
415 //____________________________________________________________________________
416 void AliEMCALClusterizer::PrintRecPoints(Option_t * option)
417 {
418   // Prints list of RecPoints produced at the current pass of AliEMCALClusterizer
419   if(strstr(option,"deb")) {
420     printf("PrintRecPoints: Clusterization result:") ; 
421     
422     printf("           Found %d ECA Rec Points\n ", 
423            fRecPoints->GetEntriesFast()) ; 
424   }
425   
426   if(strstr(option,"all")) {
427     if(strstr(option,"deb")) {
428       printf("\n-----------------------------------------------------------------------\n") ;
429       printf("Clusters in ECAL section\n") ;
430       printf("Index    Ene(GeV) Multi Module     GX    GY   GZ  lX    lY   lZ   Dispersion Lambda 1   Lambda 2  # of prim  Primaries list\n") ;
431     }
432     Int_t index; 
433     for (index =  0 ; index < fRecPoints->GetEntries() ; index++) {
434       AliEMCALRecPoint * rp = dynamic_cast<AliEMCALRecPoint * >(fRecPoints->At(index)) ; 
435       if(rp){
436         TVector3  globalpos;  
437         //rp->GetGlobalPosition(globalpos);
438         TVector3  localpos;  
439         rp->GetLocalPosition(localpos);
440         Float_t lambda[2]; 
441         rp->GetElipsAxis(lambda);
442         
443         Int_t nprimaries=0;
444         Int_t * primaries = rp->GetPrimaries(nprimaries);
445         
446         if(strstr(option,"deb")) 
447           printf("\n%6d  %8.4f  %3d     %4.1f    %4.1f %4.1f  %4.1f %4.1f %4.1f    %4.1f   %4f  %4f    %2d     : ", 
448                  rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(),
449                  globalpos.X(), globalpos.Y(), globalpos.Z(), localpos.X(), localpos.Y(), localpos.Z(), 
450                  rp->GetDispersion(), lambda[0], lambda[1], nprimaries) ; 
451         if(strstr(option,"deb")){ 
452           for (Int_t iprimary=0; iprimary<nprimaries; iprimary++) {
453             printf("%d ", primaries[iprimary] ) ; 
454           }
455         }
456       }
457     }
458     
459     if(strstr(option,"deb"))
460       printf("\n-----------------------------------------------------------------------\n");
461   }
462 }
463
464 //___________________________________________________________________
465 void  AliEMCALClusterizer::PrintRecoInfo()
466 {
467   //print reco version
468   printf(" AliEMCALClusterizer::PrintRecoInfo() : version %s \n", Version() );
469   
470 }
471
472 //____________________________________________________________________________
473 void AliEMCALClusterizer::SetInput(TTree *digitsTree)
474 {
475   // Read the digits from the input tree
476   TBranch *branch = digitsTree->GetBranch("EMCAL");
477   if (!branch) { 
478     AliError("can't get the branch with the EMCAL digits !");
479     return;
480   }
481   if (!fDigitsArr)
482     fDigitsArr = new TClonesArray("AliEMCALDigit",100);
483   branch->SetAddress(&fDigitsArr);
484   branch->GetEntry(0);
485 }
486
487 //____________________________________________________________________________
488 void AliEMCALClusterizer::SetOutput(TTree *clustersTree)
489 {
490   // Read the digits from the input tree
491   fTreeR = clustersTree;
492   
493   AliDebug(9, "Making array for EMCAL clusters");
494   fRecPoints = new TObjArray(100) ;
495   Int_t split = 0;
496   Int_t bufsize = 32000;
497   fTreeR->Branch("EMCALECARP", "TObjArray", &fRecPoints, bufsize, split);
498 }
499
500 //___________________________________________________________________
501 void   AliEMCALClusterizer::SetInputCalibrated(Bool_t val)
502 {
503   //
504   // input is calibrated - the case when we run already on ESD
505   //
506   AliEMCALClusterizer::fgkIsInputCalibrated = val;
507 }
508
509
510