Correct leak in AliCaloCalibPedestal
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALReconstructor.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 //--
20 //-- Yves Schutz (SUBATECH) 
21 // Reconstruction class. Redesigned from the old AliReconstructionner class and 
22 // derived from STEER/AliReconstructor. 
23 // 
24 //-- Aleksei Pavlinov : added staf for EMCAL jet trigger 9Apr 25, 2008)
25 //                    : fgDigitsArr should read just once at event
26
27 // --- ROOT system ---
28 #include <TList.h>
29 #include <TClonesArray.h>
30 #include <TH2.h>
31 #include "TGeoManager.h"
32 #include "TGeoMatrix.h"
33
34 // --- Standard library ---
35
36 // --- AliRoot header files ---
37 #include "AliEMCALReconstructor.h"
38
39 #include "AliCodeTimer.h"
40 #include "AliCaloCalibPedestal.h"
41 #include "AliEMCALCalibData.h"
42 #include "AliESDEvent.h"
43 #include "AliESDCaloCluster.h"
44 #include "AliESDCaloCells.h"
45 #include "AliESDtrack.h"
46 #include "AliEMCALLoader.h"
47 #include "AliEMCALRawUtils.h"
48 #include "AliEMCALDigit.h"
49 #include "AliEMCALClusterizerv1.h"
50 #include "AliEMCALClusterizerNxN.h"
51 #include "AliEMCALRecPoint.h"
52 #include "AliEMCALPID.h"
53 #include "AliEMCALTrigger.h"
54 #include "AliRawReader.h"
55 #include "AliCDBEntry.h"
56 #include "AliCDBManager.h"
57 #include "AliEMCALGeometry.h"
58 #include "AliEMCAL.h"
59 #include "AliESDVZERO.h"
60 #include "AliCDBManager.h"
61 #include "AliRunLoader.h"
62 #include "AliRun.h"
63 #include "AliEMCALTriggerData.h"
64 #include "AliEMCALTriggerElectronics.h"
65 #include "AliEMCALTriggerDCSConfigDB.h"
66 #include "AliEMCALTriggerDCSConfig.h"
67
68 ClassImp(AliEMCALReconstructor) 
69
70 const AliEMCALRecParam*     AliEMCALReconstructor::fgkRecParam        = 0;   // EMCAL rec. parameters
71 AliEMCALRawUtils*           AliEMCALReconstructor::fgRawUtils         = 0;   // EMCAL raw utilities class
72 AliEMCALClusterizer*        AliEMCALReconstructor::fgClusterizer      = 0;   // EMCAL clusterizer class
73 TClonesArray*               AliEMCALReconstructor::fgDigitsArr        = 0;   // list of digits, to be used multiple times
74 TObjArray*                  AliEMCALReconstructor::fgClustersArr      = 0;   // list of clusters, to be used multiple times
75 AliEMCALTriggerElectronics* AliEMCALReconstructor::fgTriggerProcessor = 0x0;
76 //____________________________________________________________________________
77 AliEMCALReconstructor::AliEMCALReconstructor() 
78   : fDebug(kFALSE), fList(0), fGeom(0),fCalibData(0),fPedestalData(0) 
79 {
80   // ctor
81
82   fgRawUtils = new AliEMCALRawUtils;
83
84   //To make sure we match with the geometry in a simulation file,
85   //let's try to get it first.  If not, take the default geometry
86   AliRunLoader *rl = AliRunLoader::Instance();
87   if (rl->GetAliRun() && rl->GetAliRun()->GetDetector("EMCAL")) {
88     fGeom = dynamic_cast<AliEMCAL*>(rl->GetAliRun()->GetDetector("EMCAL"))->GetGeometry();
89   } else {
90     AliInfo(Form("Using default geometry in reconstruction"));
91     fGeom =  AliEMCALGeometry::GetInstance(AliEMCALGeometry::GetDefaultGeometryName());
92   }
93
94   //Get calibration parameters  
95   if(!fCalibData)
96     {
97                 AliCDBEntry *entry = (AliCDBEntry*) 
98                 AliCDBManager::Instance()->Get("EMCAL/Calib/Data");
99                 if (entry) fCalibData =  (AliEMCALCalibData*) entry->GetObject();
100     }
101         
102   if(!fCalibData)
103                 AliFatal("Calibration parameters not found in CDB!");
104         
105   //Get calibration parameters  
106   if(!fPedestalData)
107     {
108                 AliCDBEntry *entry = (AliCDBEntry*) 
109                 AliCDBManager::Instance()->Get("EMCAL/Calib/Pedestals");
110                 if (entry) fPedestalData =  (AliCaloCalibPedestal*) entry->GetObject();
111     }
112         
113   if(!fPedestalData)
114     AliFatal("Dead map not found in CDB!");
115         
116   InitClusterizer();
117                 
118   if(!fGeom) AliFatal(Form("Could not get geometry!"));
119
120   AliEMCALTriggerDCSConfigDB* dcsConfigDB = AliEMCALTriggerDCSConfigDB::Instance();
121
122   const AliEMCALTriggerDCSConfig* dcsConfig = dcsConfigDB->GetTriggerDCSConfig();
123
124   if (!dcsConfig) AliFatal("No Trigger DCS Configuration from OCDB!");
125   fgTriggerProcessor = new AliEMCALTriggerElectronics( dcsConfig );
126   
127   //Init temporary list of digits
128   fgDigitsArr   = new TClonesArray("AliEMCALDigit",1000);
129   fgClustersArr = new TObjArray(1000);
130
131
132
133 //____________________________________________________________________________
134 AliEMCALReconstructor::~AliEMCALReconstructor()
135 {
136   // dtor
137
138   if(fGeom)              delete fGeom;
139   if(fCalibData)         delete fCalibData;
140   if(fPedestalData)      delete fPedestalData;
141
142   if(fgDigitsArr){
143     fgDigitsArr->Clear("C");
144     delete fgDigitsArr; 
145   }
146   
147   if(fgClustersArr){
148     fgClustersArr->Clear();
149     delete fgClustersArr; 
150   }
151   
152   if(fgRawUtils)         delete fgRawUtils;
153   if(fgClusterizer)      delete fgClusterizer;
154   if(fgTriggerProcessor) delete fgTriggerProcessor;
155   
156   AliCodeTimer::Instance()->Print();
157
158
159 // //____________________________________________________________________________
160 // void AliEMCALReconstructor::Init()
161 // {
162 //   // Trigger hists - Oct 24, 2007
163 //    fList = AliEMCALHistoUtilities::GetTriggersListOfHists(kTRUE);
164 // }
165
166 //____________________________________________________________________________
167 void AliEMCALReconstructor::InitClusterizer() 
168 {
169   //Init the clusterizer with geometry and calibration pointers, avoid doing it twice.
170   
171   AliEMCALRecParam *recParam = NULL;
172   AliCDBEntry *entry = (AliCDBEntry*) 
173   AliCDBManager::Instance()->Get("EMCAL/Calib/RecoParam");
174   //Get The reco param for the default event specie
175   if (entry) 
176    recParam = (AliEMCALRecParam*)((TObjArray *) entry->GetObject())->At(0);
177     
178   if(!recParam)  
179     AliFatal("RecoParam not found in CDB!");
180     
181   if (recParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv1)
182   {
183     fgClusterizer = new AliEMCALClusterizerv1(fGeom, fCalibData,fPedestalData); 
184   }
185   else
186   {
187     fgClusterizer = new AliEMCALClusterizerNxN(fGeom, fCalibData,fPedestalData); 
188   }
189   
190 }
191
192 //____________________________________________________________________________
193 void AliEMCALReconstructor::Reconstruct(TTree* digitsTree, TTree* clustersTree) const
194 {
195   // method called by AliReconstruction; 
196   // Only the clusterization is performed,; the rest of the reconstruction is done in FillESD because the track
197   // segment maker needs access to the AliESD object to retrieve the tracks reconstructed by 
198   // the global tracking.
199   // Works on the current event.
200
201   AliCodeTimerAuto("",0)
202
203   ReadDigitsArrayFromTree(digitsTree);
204   
205   fgClusterizer->InitParameters();
206   fgClusterizer->SetOutput(clustersTree);
207  
208   AliEMCALTriggerData* trgData = new AliEMCALTriggerData();
209         
210   Int_t bufferSize = 32000;
211
212   if (TBranch* triggerBranch = clustersTree->GetBranch("EMTRG"))
213                 triggerBranch->SetAddress(&trgData);
214         else
215                 clustersTree->Branch("EMTRG","AliEMCALTriggerData",&trgData,bufferSize);
216
217   TClonesArray *trgDigits = new TClonesArray("AliEMCALRawDigit",1000);
218   TBranch *branchdig = digitsTree->GetBranch("EMTRG");
219   if (!branchdig) 
220   { 
221           AliError("Can't get the branch with the EMCAL trigger digits !");
222           return;
223   }
224
225   branchdig->SetAddress(&trgDigits);
226   branchdig->GetEntry(0);
227         
228   //Skip clusterization of LED events
229   if (GetRecParam()->GetEventSpecie()!=AliRecoParam::kCalib){
230
231           Int_t v0M[2] = {0,0};
232           fgTriggerProcessor->Digits2Trigger(trgDigits, v0M, trgData);
233         
234  
235           if(fgDigitsArr && fgDigitsArr->GetEntries()) {
236
237                   fgClusterizer->SetInput(digitsTree);
238     
239                   if(Debug())
240                           fgClusterizer->Digits2Clusters("deb all") ;
241                   else
242                           fgClusterizer->Digits2Clusters("");
243     
244                   fgClusterizer->Clear();
245
246           }//digits array exists and has somethind
247   }//not a LED event
248         
249   clustersTree->Fill(); 
250   trgDigits->Delete();
251   delete trgDigits; trgDigits = 0x0;
252   delete trgData;   trgData   = 0x0;
253
254 }
255
256 //____________________________________________________________________________
257 void AliEMCALReconstructor::ConvertDigits(AliRawReader* rawReader, TTree* digitsTree) const
258
259 {
260   // Conversion from raw data to
261   // EMCAL digits.
262   // Works on a single-event basis
263
264   rawReader->Reset() ; 
265   if(fgDigitsArr) fgDigitsArr->Clear("C");
266   
267   TClonesArray *digitsTrg = new TClonesArray("AliEMCALRawDigit", 200);
268
269   Int_t bufsize = 32000;
270   digitsTree->Branch("EMCAL", &fgDigitsArr, bufsize);
271   digitsTree->Branch("EMTRG", &digitsTrg, bufsize);
272         
273   //Skip calibration events do the rest
274   Bool_t doFit = kTRUE;
275   if ( !(GetRecParam()->FitLEDEvents()) && GetRecParam()->GetEventSpecie()==AliRecoParam::kCalib) doFit = kFALSE;
276   if (doFit){
277           //must be done here because, in constructor, option is not yet known
278           fgRawUtils->SetOption(GetOption());
279
280           fgRawUtils->SetRawFormatHighLowGainFactor(GetRecParam()->GetHighLowGainFactor());
281           fgRawUtils->SetRawFormatOrder(GetRecParam()->GetOrderParameter());
282           fgRawUtils->SetRawFormatTau(GetRecParam()->GetTau());
283           fgRawUtils->SetNoiseThreshold(GetRecParam()->GetNoiseThreshold());
284           fgRawUtils->SetNPedSamples(GetRecParam()->GetNPedSamples());
285           fgRawUtils->SetRemoveBadChannels(GetRecParam()->GetRemoveBadChannels());
286           fgRawUtils->SetFittingAlgorithm(GetRecParam()->GetFittingAlgorithm());
287           fgRawUtils->SetFALTROUsage(GetRecParam()->UseFALTRO());
288           fgRawUtils->SetTimeMin(GetRecParam()->GetTimeMin());
289           fgRawUtils->SetTimeMax(GetRecParam()->GetTimeMax());
290         
291           fgRawUtils->Raw2Digits(rawReader,fgDigitsArr,fPedestalData,digitsTrg);
292   }//skip calibration event
293   else{
294         AliDebug(1," Calibration Event, skip!");
295   }
296         
297   digitsTree->Fill();
298   digitsTrg->Delete();
299   delete digitsTrg;
300
301 }
302
303
304 //____________________________________________________________________________
305 void AliEMCALReconstructor::FillESD(TTree* digitsTree, TTree* clustersTree, 
306                                     AliESDEvent* esd) const
307 {
308   // Called by AliReconstruct after Reconstruct() and global tracking and vertexing 
309   // and V0 
310   // Works on the current event
311   //  printf(" ## AliEMCALReconstructor::FillESD() is started ### \n ");
312   //return;
313
314   //########################################
315   //##############Fill CaloCells###############
316   //########################################
317   ReadDigitsArrayFromTree(digitsTree);
318
319 //  TClonesArray *digits = new TClonesArray("AliEMCALDigit",1000);
320 //  TBranch *branchdig = digitsTree->GetBranch("EMCAL");
321 //  if (!branchdig) { 
322 //    AliError("can't get the branch with the EMCAL digits !");
323 //    return;
324 //  }
325 //  branchdig->SetAddress(&digits);
326 //  digitsTree->GetEvent(0);
327   Int_t nDigits = fgDigitsArr->GetEntries(), idignew = 0 ;
328   AliDebug(1,Form("%d digits",nDigits));
329
330   AliESDCaloCells &emcCells = *(esd->GetEMCALCells());
331   emcCells.CreateContainer(nDigits);
332   emcCells.SetType(AliVCaloCells::kEMCALCell);
333   Float_t energy = 0;
334   for (Int_t idig = 0 ; idig < nDigits ; idig++) {
335     const AliEMCALDigit * dig = (const AliEMCALDigit*)fgDigitsArr->At(idig);
336     if(dig->GetAmplitude() > 0 ){
337           energy = fgClusterizer->Calibrate(dig->GetAmplitude(),dig->GetTime(),dig->GetId()); //TimeR or Time?
338           if(energy > 0){ //Digits tagged as bad (dead, hot, not alive) are set to 0 in calibrate, remove them  
339                   emcCells.SetCell(idignew,dig->GetId(),energy, dig->GetTime());   
340                   idignew++;
341           }
342     }
343   }
344   emcCells.SetNumberOfCells(idignew);
345   emcCells.Sort();
346
347   //------------------------------------------------------------
348   //-----------------CLUSTERS-----------------------------
349   //------------------------------------------------------------
350   clustersTree->SetBranchStatus("*",0); //disable all branches
351   clustersTree->SetBranchStatus("EMCALECARP",1); //Enable only the branch we need
352   if(fgClustersArr) fgClustersArr->Clear();
353   TBranch *branch = clustersTree->GetBranch("EMCALECARP");
354   branch->SetAddress(&fgClustersArr);
355   branch->GetEntry(0);
356   //clustersTree->GetEvent(0);
357
358   Int_t nClusters = fgClustersArr->GetEntries(),  nClustersNew=0;
359   AliDebug(1,Form("%d clusters",nClusters));
360
361   //######################################################
362   //#######################TRACK MATCHING###############
363   //######################################################
364   //Fill list of integers, each one is index of track to which the cluster belongs.
365
366   // step 1 - initialize array of matched track indexes
367   Int_t *matchedTrack = new Int_t[nClusters];
368   for (Int_t iclus = 0; iclus < nClusters; iclus++)
369     matchedTrack[iclus] = -1;  // neg. index --> no matched track
370   
371   // step 2, change the flag for all matched clusters found in tracks
372   Int_t iemcalMatch = -1;
373   Int_t endtpc = esd->GetNumberOfTracks();
374   for (Int_t itrack = 0; itrack < endtpc; itrack++) {
375     AliESDtrack * track = esd->GetTrack(itrack) ; // retrieve track
376     iemcalMatch = track->GetEMCALcluster();
377     if(iemcalMatch >= 0) matchedTrack[iemcalMatch] = itrack;
378   } 
379   
380   //########################################
381   //##############Fill CaloClusters#############
382   //########################################
383   for (Int_t iClust = 0 ; iClust < nClusters ; iClust++) {
384     const AliEMCALRecPoint * clust = (const AliEMCALRecPoint*)fgClustersArr->At(iClust);
385     //if(clust->GetClusterType()== AliVCluster::kEMCALClusterv1) nRP++; else nPC++;
386     if (Debug()) clust->Print();
387     // Get information from EMCAL reconstruction points
388     Float_t xyz[3];
389     TVector3 gpos;
390     clust->GetGlobalPosition(gpos);
391     for (Int_t ixyz=0; ixyz<3; ixyz++)
392       xyz[ixyz] = gpos[ixyz];
393     Float_t elipAxis[2];
394     clust->GetElipsAxis(elipAxis);
395        //Create digits lists
396     Int_t cellMult = clust->GetMultiplicity();
397     //TArrayS digiList(digitMult);
398     Float_t *amplFloat = clust->GetEnergiesList();
399     Int_t   *digitInts = clust->GetAbsId();
400     TArrayS absIdList(cellMult);
401     TArrayD fracList(cellMult);
402
403     Int_t newCellMult = 0;
404     for (Int_t iCell=0; iCell<cellMult; iCell++) {
405       if (amplFloat[iCell] > 0) {
406       absIdList[newCellMult] = (UShort_t)(digitInts[iCell]);
407       //Uncomment when unfolding is done
408       //if(emcCells.GetCellAmplitude(digitInts[iCell])>0)
409       //fracList[newCellMult] = amplFloat[iCell]/(emcCells.GetCellAmplitude(digitInts[iCell])*calibration);//get cell calibration value 
410       //else
411       fracList[newCellMult] = 0; 
412       newCellMult++;
413       }
414     }
415
416     absIdList.Set(newCellMult);
417     fracList.Set(newCellMult);
418     
419     if(newCellMult > 0) { // accept cluster if it has some digit
420       nClustersNew++;
421       //Primaries
422       Int_t  parentMult  = 0;
423       Int_t *parentList =  clust->GetParents(parentMult);
424       // fills the ESDCaloCluster
425       AliESDCaloCluster * ec = new AliESDCaloCluster() ;
426       ec->SetType(AliVCluster::kEMCALClusterv1);
427       ec->SetPosition(xyz);
428       ec->SetE(clust->GetEnergy());
429                 
430       //Distance to the nearest bad crystal
431       ec->SetDistanceToBadChannel(clust->GetDistanceToBadTower()); 
432
433       ec->SetNCells(newCellMult);
434       //Change type of list from short to ushort
435       UShort_t *newAbsIdList  = new UShort_t[newCellMult];
436       Double_t *newFracList   = new Double_t[newCellMult];
437       for(Int_t i = 0; i < newCellMult ; i++) {
438         newAbsIdList[i]=absIdList[i];
439         newFracList[i]=fracList[i];
440       }
441       ec->SetCellsAbsId(newAbsIdList);
442       ec->SetCellsAmplitudeFraction(newFracList);
443       ec->SetDispersion(clust->GetDispersion());
444       ec->SetChi2(-1); //not yet implemented
445       ec->SetM02(elipAxis[0]*elipAxis[0]) ;
446       ec->SetM20(elipAxis[1]*elipAxis[1]) ;
447       ec->SetTOF(clust->GetTime()) ; //time-of-fligh
448       ec->SetNExMax(clust->GetNExMax());          //number of local maxima
449       TArrayI arrayTrackMatched(1);// Only one track, temporal solution.
450       arrayTrackMatched[0]= matchedTrack[iClust];
451       ec->AddTracksMatched(arrayTrackMatched);
452
453       TArrayI arrayParents(parentMult,parentList);
454       ec->AddLabels(arrayParents);
455
456       // add the cluster to the esd object
457       esd->AddCaloCluster(ec);
458       delete ec;
459       delete [] newAbsIdList ;
460       delete [] newFracList ;
461    }
462  } // cycle on clusters
463
464  delete [] matchedTrack;
465
466  //Fill ESDCaloCluster with PID weights
467  AliEMCALPID *pid = new AliEMCALPID;
468  //pid->SetPrintInfo(kTRUE);
469  pid->SetReconstructor(kTRUE);
470  pid->RunPID(esd);
471  delete pid;
472   
473  //delete digits;
474  //delete clusters;
475   
476  //Store EMCAL misalignment matrixes
477  FillMisalMatrixes(esd) ;
478
479 }
480
481 //==================================================================================
482 void AliEMCALReconstructor::FillMisalMatrixes(AliESDEvent* esd)const{
483         //Store EMCAL matrixes in ESD Header
484         
485         //Check, if matrixes was already stored
486         for(Int_t sm = 0 ; sm < fGeom->GetNumberOfSuperModules(); sm++){
487                 if(esd->GetEMCALMatrix(sm)!=0)
488                         return ;
489         }
490         
491         //Create and store matrixes
492         if(!gGeoManager){
493                 AliError("Can not store misal. matrixes: no gGeoManager! \n") ;
494                 return ;
495         }
496         //Note, that owner of copied marixes will be header
497         char path[255] ;
498         TGeoHMatrix * m = 0x0;
499         for(Int_t sm = 0; sm < fGeom->GetNumberOfSuperModules(); sm++){
500                 sprintf(path,"/ALIC_1/XEN1_1/SMOD_%d",sm+1) ; //In Geometry modules numbered 1,2,.,5
501                 if(sm >= 10) sprintf(path,"/ALIC_1/XEN1_1/SM10_%d",sm-10+1) ;
502                 
503                 if (gGeoManager->CheckPath(path)){
504                         gGeoManager->cd(path);
505                         m = gGeoManager->GetCurrentMatrix() ;
506 //                      printf("================================================= \n");
507 //                      printf("AliEMCALReconstructor::FixMisalMatrixes(), sm %d, \n",sm);
508 //                      m->Print("");
509                         esd->SetEMCALMatrix(new TGeoHMatrix(*m),sm) ;
510 //                      printf("================================================= \n");
511                 }
512                 else{
513                         esd->SetEMCALMatrix(NULL,sm) ;
514                 }
515         }
516 }
517
518
519
520 //__________________________________________________________________________
521 void AliEMCALReconstructor::ReadDigitsArrayFromTree(TTree *digitsTree) const
522 {
523   // See AliEMCALClusterizer::SetInput(TTree *digitsTree);
524   if(fgDigitsArr) {
525     // Clear previous digits 
526     fgDigitsArr->Clear("C");
527     //delete fgDigitsArr;
528   }
529   // Read the digits from the input tree
530   TBranch *branch = digitsTree->GetBranch("EMCAL");
531   if (!branch) { 
532     AliError("can't get the branch with the EMCAL digits !");
533     return;
534   }
535   fgDigitsArr = new TClonesArray("AliEMCALDigit",100);
536   branch->SetAddress(&fgDigitsArr);
537   branch->GetEntry(0);
538 }
539
540