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