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