Adding include path to allow compilation of CleanGeom task
[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.),fRejectBelowThreshold(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.),fRejectBelowThreshold(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.),fRejectBelowThreshold(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   //Return energy with default parameters if calibration is not available
202   if (!fCalibData && !fCaloPed) {
203     if (fIsInputCalibrated == kTRUE)
204     {
205       AliDebug(10, Form("Input already calibrated!"));
206       return ;
207     }
208     else{
209       AliFatal("OCDB calibration and bad map parameters are not available");
210       return;
211     }   
212   }
213   
214   if (fGeom==0)
215     AliFatal("Did not get geometry from EMCALLoader") ;
216     
217   Int_t iSupMod = -1;
218   Int_t nModule = -1;
219   Int_t nIphi   = -1;
220   Int_t nIeta   = -1;
221   Int_t iphi    = -1;
222   Int_t ieta    = -1;
223     
224   Bool_t bCell = fGeom->GetCellIndex(absId, iSupMod, nModule, nIphi, nIeta) ;
225   if(!bCell) {
226     fGeom->PrintGeometry();
227     AliError(Form("Wrong cell id number : %i", absId));
228     //assert(0); // GCB: This aborts reconstruction of raw simulations 
229     //where simulation had more SM than default geometry, 
230     //change to return 0, to avoid aborting good generations.
231     amp  = 0;
232     time = 0;
233     return ;
234   }
235     
236   fGeom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
237           
238   // Check if channel is bad (dead or hot), in this case return 0.      
239   // Gustavo: 15-12-09 In case of RAW data this selection is already done, but not in simulation.
240   // for the moment keep it here but remember to do the selection at the sdigitizer level 
241   // and remove it from here
242   if (fCaloPed) {
243     Int_t channelStatus = (Int_t)(fCaloPed->GetDeadMap(iSupMod))->GetBinContent(ieta,iphi);
244     if(channelStatus == AliCaloCalibPedestal::kHot || channelStatus == AliCaloCalibPedestal::kDead) {
245       AliDebug(2,Form("Tower from SM %d, ieta %d, iphi %d is BAD : status %d !!!",iSupMod,ieta,iphi, channelStatus));
246       amp  = 0 ;
247       time = 0 ;
248       return ;
249     }
250   }
251     
252   if (fIsInputCalibrated || !fCalibData)
253   {
254     AliDebug(10, Form("Input already calibrated!"));
255     return ;
256   }
257           
258   Int_t bc = 0; // Get somehow the bunch crossing number
259
260   fADCchannelECA  = fCalibData->GetADCchannel (iSupMod,ieta,iphi);
261   fADCpedestalECA = fCalibData->GetADCpedestal(iSupMod,ieta,iphi);
262   fTimeECA        = fCalibData->GetTimeChannel(iSupMod,ieta,iphi, bc);
263   
264   time -= fTimeECA ;
265   amp   = amp * fADCchannelECA - fADCpedestalECA ;  
266   
267 }
268
269 //____________________________________________________________________________
270 void AliEMCALClusterizer::GetCalibrationParameters() 
271 {
272   // Set calibration parameters:
273   // If calibration database exists, they are read from database,
274   // otherwise, they are taken from digitizer.
275   // It is a user responsilibity to open CDB before reconstruction, 
276   // for example: 
277   // AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB");
278
279   if (fIsInputCalibrated)
280     return;
281   
282   //Check if calibration is stored in data base
283   if(!fCalibData)
284   {
285     AliCDBEntry *entry = (AliCDBEntry*) 
286     AliCDBManager::Instance()->Get("EMCAL/Calib/Data");
287     if (entry) fCalibData =  (AliEMCALCalibData*) entry->GetObject();
288   }
289   
290   if(!fCalibData)
291     AliFatal("Calibration parameters not found in CDB!");
292 }
293
294 //____________________________________________________________________________
295 void AliEMCALClusterizer::GetCaloCalibPedestal() 
296 {
297   // Set calibration parameters:
298   // if calibration database exists, they are read from database,
299   // otherwise, they are taken from digitizer.
300   //
301   // It is a user responsilibity to open CDB before reconstruction, 
302   // for example: 
303   // AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB");
304
305   if (fIsInputCalibrated)
306     return;
307   
308   // Check if calibration is stored in data base
309   if(!fCaloPed)
310     {
311       AliCDBEntry *entry = (AliCDBEntry*) 
312         AliCDBManager::Instance()->Get("EMCAL/Calib/Pedestals");
313       if (entry) fCaloPed =  (AliCaloCalibPedestal*) entry->GetObject();
314     }
315   
316   if(!fCaloPed)
317     AliFatal("Pedestal info not found in CDB!");
318 }
319
320 //____________________________________________________________________________
321 void AliEMCALClusterizer::Init()
322 {
323   // Make all memory allocations which can not be done in default constructor.
324   // Attach the Clusterizer task to the list of EMCAL tasks
325   
326   AliRunLoader *rl = AliRunLoader::Instance();
327   if (rl->GetAliRun()){
328     AliEMCAL* emcal = dynamic_cast<AliEMCAL*>(rl->GetAliRun()->GetDetector("EMCAL"));
329     if(emcal)fGeom = emcal->GetGeometry();
330   }
331   
332   if(!fGeom){ 
333     fGeom = AliEMCALGeometry::GetInstance(AliEMCALGeometry::GetDefaultGeometryName());
334   }
335   
336   AliDebug(1,Form("geom %p",fGeom));
337   
338   if(!gMinuit) 
339     gMinuit = new TMinuit(100) ;
340   
341   Int_t i=0;
342   for (i = 0; i < 8; i++)
343     fSSPars[i] = 0.;
344   for (i = 0; i < 3; i++) {
345     fPar5[i] = 0.;
346     fPar6[i] = 0.;
347   }
348 }
349
350 //____________________________________________________________________________
351 void AliEMCALClusterizer::InitParameters()
352
353   // Initializes the parameters for the Clusterizer from AliEMCALReconstructor::GetRecParam().
354
355   return InitParameters(AliEMCALReconstructor::GetRecParam());
356 }
357
358 //____________________________________________________________________________
359 void AliEMCALClusterizer::InitParameters(const AliEMCALRecParam* recParam)
360
361   // Initializes the parameters for the Clusterizer
362
363   fNumberOfECAClusters = 0 ;
364   fCalibData           = 0 ;
365   fCaloPed             = 0 ;
366         
367   if(!recParam) {
368     AliFatal("Reconstruction parameters for EMCAL not set!");
369   } 
370
371   fECAClusteringThreshold = recParam->GetClusteringThreshold();
372   fECAW0                  = recParam->GetW0();
373   fMinECut                = recParam->GetMinECut();    
374   fRejectBelowThreshold   = recParam->GetRejectBelowThreshold();
375   fToUnfold               = recParam->GetUnfold();
376   fECALocMaxCut           = recParam->GetLocMaxCut();
377   fTimeCut                = recParam->GetTimeCut();
378   fTimeMin                = recParam->GetTimeMin();
379   fTimeMax                = recParam->GetTimeMax();
380   
381   //For NxN
382   SetNRowDiff(recParam->GetNRowDiff());
383   SetNColDiff(recParam->GetNColDiff());
384   
385   AliDebug(1,Form("Reconstruction parameters: fECAClusteringThreshold=%.3f GeV, fECAW=%.3f, fMinECut=%.3f GeV, "
386                   "fToUnfold=%d, fECALocMaxCut=%.3f GeV, fTimeCut=%e s,fTimeMin=%e s,fTimeMax=%e s,fRejectBelowThreshold=%d",
387                   fECAClusteringThreshold,fECAW0,fMinECut,fToUnfold,fECALocMaxCut,fTimeCut, fTimeMin, fTimeMax, 
388                   fRejectBelowThreshold));
389
390   if (fToUnfold) {
391     Int_t i=0;
392     for (i = 0; i < 8; i++) {
393       fSSPars[i] = recParam->GetSSPars(i);
394     } //end of loop over parameters
395     for (i = 0; i < 3; i++) {
396       fPar5[i] = recParam->GetPar5(i);
397       fPar6[i] = recParam->GetPar6(i);
398     } //end of loop over parameters
399       
400     InitClusterUnfolding();
401       
402     for (i = 0; i < 8; i++) {
403       AliDebug(1,Form("unfolding shower shape parameters: fSSPars=%f \n",fSSPars[i]));
404     }
405     for (i = 0; i < 3; i++) {
406       AliDebug(1,Form("unfolding parameter 5: fPar5=%f \n",fPar5[i]));
407       AliDebug(1,Form("unfolding parameter 6: fPar6=%f \n",fPar6[i]));
408     }
409   } // to unfold
410 }
411
412 //____________________________________________________________________________
413 void AliEMCALClusterizer::Print(Option_t * /*option*/)const
414 {
415   // Print clusterizer parameters
416   
417   TString message("\n") ; 
418   
419   if (strcmp(GetName(),"") == 0) {
420     printf("AliEMCALClusterizer not initialized\n");
421     return;
422   }
423     
424   // Print parameters
425   TString taskName(Version()) ;
426     
427   printf("--------------- "); 
428   printf("%s",taskName.Data()) ; 
429   printf(" "); 
430   printf("Clusterizing digits: "); 
431   printf("\n                       ECA Local Maximum cut    = %f", fECALocMaxCut); 
432   printf("\n                       ECA Logarithmic weight   = %f", fECAW0); 
433   if (fToUnfold) {
434     printf("\nUnfolding on\n");
435     printf("Unfolding parameters: fSSpars: \n");
436     Int_t i=0;
437     for (i = 0; i < 8; i++) {
438       printf("fSSPars[%d] = %f \n", i, fSSPars[i]);
439     }
440     printf("Unfolding parameter 5 and 6: fPar5 and fPar6: \n");
441     for (i = 0; i < 3; i++) {
442       printf("fPar5[%d] = %f \n", i, fPar5[i]);
443       printf("fPar6[%d] = %f \n", i, fPar6[i]);
444     }
445   }
446   else
447     printf("\nUnfolding off\n");
448     
449   printf("------------------------------------------------------------------"); 
450 }
451
452 //____________________________________________________________________________
453 void AliEMCALClusterizer::PrintRecPoints(Option_t * option)
454 {
455   // Prints list of RecPoints produced at the current pass of AliEMCALClusterizer
456
457   if (strstr(option,"deb")) {
458     printf("PrintRecPoints: Clusterization result:") ; 
459     printf("           Found %d ECA Rec Points\n ", 
460            fRecPoints->GetEntriesFast()) ; 
461   }
462   
463   if (strstr(option,"all")) {
464     if (strstr(option,"deb")) {
465       printf("\n-----------------------------------------------------------------------\n") ;
466       printf("Clusters in ECAL section\n") ;
467       printf("Index    Ene(GeV) Multi Module     GX    GY   GZ  lX    lY   lZ   Dispersion Lambda 1   Lambda 2  # of prim  Primaries list\n") ;
468     }
469     Int_t index; 
470     for (index =  0 ; index < fRecPoints->GetEntries() ; index++) {
471       AliEMCALRecPoint * rp = dynamic_cast<AliEMCALRecPoint * >(fRecPoints->At(index)) ; 
472       if (!rp) 
473         continue;
474
475       TVector3  globalpos;  
476       //rp->GetGlobalPosition(globalpos);
477       TVector3  localpos;  
478       rp->GetLocalPosition(localpos);
479       Float_t lambda[2]; 
480       rp->GetElipsAxis(lambda);
481         
482       Int_t nprimaries=0;
483       Int_t * primaries = rp->GetPrimaries(nprimaries);
484       if(strstr(option,"deb")) 
485         printf("\n%6d  %8.4f  %3d     %4.1f    %4.1f %4.1f  %4.1f %4.1f %4.1f    %4.1f   %4f  %4f    %2d     : ", 
486                rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(),
487                globalpos.X(), globalpos.Y(), globalpos.Z(), localpos.X(), localpos.Y(), localpos.Z(), 
488                rp->GetDispersion(), lambda[0], lambda[1], nprimaries) ; 
489       if(strstr(option,"deb")){ 
490         for (Int_t iprimary=0; iprimary<nprimaries; iprimary++) {
491           printf("%d ", primaries[iprimary] ) ; 
492         }
493       }
494     }
495     
496     if(strstr(option,"deb"))
497       printf("\n-----------------------------------------------------------------------\n");
498   }
499 }
500
501 //___________________________________________________________________
502 void  AliEMCALClusterizer::PrintRecoInfo()
503 {
504   // Print reco version
505
506   printf(" AliEMCALClusterizer::PrintRecoInfo() : version %s \n", Version() );
507 }
508
509 //____________________________________________________________________________
510 void AliEMCALClusterizer::SetInput(TTree *digitsTree)
511 {
512   // Read the digits from the input tree
513
514   TBranch *branch = digitsTree->GetBranch("EMCAL");
515   if (!branch)
516   {
517     AliError("can't get the branch with the EMCAL digits !");
518     return;
519   }
520   
521   if (!fDigitsArr) fDigitsArr = new TClonesArray("AliEMCALDigit",100);
522   else             fDigitsArr->Clear("C"); // avoid leak
523     
524   branch->SetAddress(&fDigitsArr);
525   branch->GetEntry(0);
526 }
527
528 //____________________________________________________________________________
529 void AliEMCALClusterizer::SetOutput(TTree *clustersTree)
530 {
531   // Read the digits from the input tree
532
533   AliDebug(9, "Making array for EMCAL clusters");
534   fRecPoints = new TObjArray(1000);
535   if (clustersTree) {
536     fTreeR = clustersTree;
537     Int_t split = 0;
538     Int_t bufsize = 32000;
539     fTreeR->Branch("EMCALECARP", "TObjArray", &fRecPoints, bufsize, split);
540   }
541 }
542
543 //___________________________________________________________________
544 void AliEMCALClusterizer::SetInputCalibrated(Bool_t val)
545 {
546   // Flag to indicate that input is calibrated - the case when we run already on ESD
547
548   fIsInputCalibrated = val;
549 }
550
551 //___________________________________________________________________
552 void AliEMCALClusterizer::SetJustClusters(Bool_t val)
553 {
554   // Flag to indicate that we are running on ESDs, when calling 
555   // rp->EvalAll(fECAW0,fDigitsArr,fJustClusters); in derived classes
556
557   fJustClusters = val;
558 }