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