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