memory leak fixed
[u/mrichter/AliRoot.git] / PWG4 / CaloCalib / AliAnalysisTaskEMCALClusterizeFast.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 // --- Root ---
19 #include <TClonesArray.h>
20 #include <TGeoManager.h>
21 #include <TObjArray.h>
22 #include <TString.h>
23 #include <TTree.h>
24 #include "AliAODCaloCluster.h"
25 #include "AliAODEvent.h"
26 #include "AliAnalysisManager.h"
27 #include "AliCDBEntry.h"
28 #include "AliCDBManager.h"
29 #include "AliCaloCalibPedestal.h"
30 #include "AliEMCALAfterBurnerUF.h"
31 #include "AliEMCALCalibData.h"
32 #include "AliEMCALClusterizerNxN.h"
33 #include "AliEMCALClusterizerv1.h"
34 #include "AliEMCALClusterizerv2.h"
35 #include "AliEMCALDigit.h"
36 #include "AliEMCALGeometry.h"
37 #include "AliEMCALRecParam.h"
38 #include "AliEMCALRecPoint.h"
39 #include "AliEMCALRecoUtils.h"
40 #include "AliESDEvent.h"
41 #include "AliInputEventHandler.h"
42 #include "AliLog.h"
43
44 #include "AliAnalysisTaskEMCALClusterizeFast.h"
45
46 ClassImp(AliAnalysisTaskEMCALClusterizeFast)
47
48 //________________________________________________________________________
49 AliAnalysisTaskEMCALClusterizeFast::AliAnalysisTaskEMCALClusterizeFast() 
50   : AliAnalysisTaskSE(), 
51     fRun(-1),
52     fDigitsArr(0),       
53     fClusterArr(0),       
54     fRecParam(0),
55     fClusterizer(0),
56     fUnfolder(0),
57     fJustUnfold(kFALSE),
58     fGeomName(),
59     fGeomMatrixSet(kFALSE), 
60     fLoadGeomMatrices(kFALSE),
61     fOCDBpath(),
62     fCalibData(0),
63     fPedestalData(0),
64     fOutputAODBranch(0),
65     fOutputAODBrName(),
66     fRecoUtils(0),
67     fLoadCalib(0),
68     fLoadPed(0),
69     fAttachClusters(0),
70     fRecalibOnly(0),
71     fSubBackground(0)
72
73   // Constructor
74 }
75
76 //________________________________________________________________________
77 AliAnalysisTaskEMCALClusterizeFast::AliAnalysisTaskEMCALClusterizeFast(const char *name) 
78   : AliAnalysisTaskSE(name), 
79     fRun(-1),
80     fDigitsArr(0),       
81     fClusterArr(0),       
82     fRecParam(new AliEMCALRecParam),
83     fClusterizer(0),
84     fUnfolder(0),
85     fJustUnfold(kFALSE),
86     fGeomName("EMCAL_FIRSTYEARV1"),
87     fGeomMatrixSet(kFALSE), 
88     fLoadGeomMatrices(kFALSE),
89     fOCDBpath(),
90     fCalibData(0),
91     fPedestalData(0),
92     fOutputAODBranch(0),
93     fOutputAODBrName(),
94     fRecoUtils(0),
95     fLoadCalib(0),
96     fLoadPed(0),
97     fAttachClusters(0),
98     fRecalibOnly(0),
99     fSubBackground(0)
100
101   // Constructor
102
103   fBranchNames     = "ESD:AliESDHeader.,AliESDRun.,EMCALCells. AOD:header,emcalCells";
104   for(Int_t i = 0; i < 12; ++i) 
105     fGeomMatrix[i] = 0;
106 }
107
108 //________________________________________________________________________
109 AliAnalysisTaskEMCALClusterizeFast::~AliAnalysisTaskEMCALClusterizeFast()
110 {
111   // Destructor.
112
113   delete fDigitsArr; 
114   delete fClusterizer;
115   delete fUnfolder;   
116   delete fRecoUtils;
117 }
118
119 //-------------------------------------------------------------------
120 void AliAnalysisTaskEMCALClusterizeFast::UserCreateOutputObjects()
121 {
122   // Create output objects.
123
124   if (!fOutputAODBrName.IsNull()) {
125     fOutputAODBranch = new TClonesArray("AliAODCaloCluster", 0);
126     fOutputAODBranch->SetName(fOutputAODBrName);
127     AddAODBranch("TClonesArray", &fOutputAODBranch);
128     AliInfo(Form("Created Branch: %s",fOutputAODBrName.Data()));
129   }
130 }
131
132 //________________________________________________________________________
133 void AliAnalysisTaskEMCALClusterizeFast::UserExec(Option_t *) 
134 {
135   // Main loop, called for each event
136
137   // remove the contents of output list set in the previous event 
138   if (fOutputAODBranch)
139     fOutputAODBranch->Clear("C");
140
141   AliESDEvent *esdevent = dynamic_cast<AliESDEvent*>(InputEvent());
142   AliAODEvent *aodevent = dynamic_cast<AliAODEvent*>(InputEvent());
143
144   if (!esdevent&&!aodevent) {
145     Error("UserExec","Event not available");
146     return;
147   }
148
149   LoadBranches();
150
151   UInt_t offtrigger = 0;
152   if (esdevent) {
153     UInt_t mask1 = esdevent->GetESDRun()->GetDetectorsInDAQ();
154     UInt_t mask2 = esdevent->GetESDRun()->GetDetectorsInReco();
155     Bool_t desc1 = (mask1 >> 18) & 0x1;
156     Bool_t desc2 = (mask2 >> 18) & 0x1;
157     if (desc1==0 || desc2==0) { //AliDAQ::OfflineModuleName(180=="EMCAL"
158       AliError(Form("EMCAL not in DAQ/RECO: %u (%u)/%u (%u)", 
159                     mask1, esdevent->GetESDRun()->GetDetectorsInReco(),
160                     mask2, esdevent->GetESDRun()->GetDetectorsInDAQ()));
161       return;
162     }
163     AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
164     offtrigger = ((AliInputEventHandler*)(am->GetInputEventHandler()))->IsEventSelected();
165   } else if (aodevent) {
166     offtrigger =  aodevent->GetHeader()->GetOfflineTrigger();
167   }
168   if (offtrigger & AliVEvent::kFastOnly) {
169     AliWarning(Form("EMCAL not in fast only partition"));
170     return;
171   }
172   
173   Init();
174
175   if (fJustUnfold) {
176     AliWarning("Unfolding not implemented");
177     return;
178   }
179
180   FillDigitsArray();
181
182   if (fRecalibOnly) {
183     UpdateCells();
184     return; // not requested to run clusterizer
185   }
186
187   Clusterize();
188   UpdateCells();
189   UpdateClusters();
190
191   if (fOutputAODBranch)
192     RecPoints2Clusters(fOutputAODBranch);
193 }
194
195 //________________________________________________________________________
196 void AliAnalysisTaskEMCALClusterizeFast::Clusterize()
197 {
198   // Clusterize
199
200   if (fSubBackground) {
201     fClusterizer->SetInputCalibrated(kTRUE);   
202     fClusterizer->SetCalibrationParameters(0);
203   }
204
205   fClusterizer->Digits2Clusters("");
206   if (fSubBackground) {
207     if (fCalibData) {
208       fClusterizer->SetInputCalibrated(kFALSE);   
209       fClusterizer->SetCalibrationParameters(fCalibData);
210     }
211   }
212 }
213
214 //________________________________________________________________________
215 void AliAnalysisTaskEMCALClusterizeFast::FillDigitsArray()
216 {
217   // Fill digits from cells.
218
219   fDigitsArr->Clear("C");
220   AliVCaloCells *cells = InputEvent()->GetEMCALCells();
221   Double_t avgE = 0; // for background subtraction
222   Int_t ncells = cells->GetNumberOfCells();
223   for (Int_t icell = 0, idigit = 0; icell < ncells; ++icell) {
224     Double_t cellAmplitude=0, cellTime=0;
225     Short_t cellNumber=0;
226     if (cells->GetCell(icell, cellNumber, cellAmplitude, cellTime) != kTRUE)
227       break;
228     AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->New(idigit));
229     digit->SetId(cellNumber);
230     digit->SetTime(cellTime);
231     digit->SetTimeR(cellTime);
232     digit->SetIndexInList(idigit);
233     digit->SetType(AliEMCALDigit::kHG);
234     if (fRecalibOnly||fSubBackground) {
235       Float_t energy = cellAmplitude;
236       Float_t time    = cellTime;
237       fClusterizer->Calibrate(energy,time,cellNumber);
238       digit->SetAmplitude(energy);
239       avgE += energy;
240     } else {
241       digit->SetAmplitude(cellAmplitude);
242     }
243     idigit++;
244   }
245
246   if (fSubBackground) {
247     avgE /= AliEMCALGeometry::GetInstance(fGeomName)->GetNumberOfSuperModules()*48*24;
248     Int_t ndigis = fDigitsArr->GetEntries();
249     for (Int_t i = 0; i < ndigis; ++i) {
250       AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(i));
251       Double_t energy = digit->GetAmplitude() - avgE;
252       if (energy<=0.001) {
253         digit->SetAmplitude(0);
254       } else {
255         digit->SetAmplitude(energy);
256       }
257     }
258   }
259 }
260
261 //________________________________________________________________________________________
262 void AliAnalysisTaskEMCALClusterizeFast::RecPoints2Clusters(TClonesArray *clus)
263 {
264   // Cluster energy, global position, cells and their amplitude fractions are restored.
265
266   Bool_t esdobjects = 0;
267   if (strcmp(clus->GetClass()->GetName(),"AliESDCaloCluster")==0)
268     esdobjects = 1;
269
270   AliVCaloCells *cells = InputEvent()->GetEMCALCells();
271   AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance(fGeomName);
272   
273   Int_t Ncls = fClusterArr->GetEntriesFast();
274   for(Int_t i=0, nout=clus->GetEntries(); i < Ncls; ++i) {
275     AliEMCALRecPoint *recpoint = static_cast<AliEMCALRecPoint*>(fClusterArr->At(i));
276     Int_t ncells_true = 0;
277     const Int_t ncells = recpoint->GetMultiplicity();
278     UShort_t   absIds[ncells];  
279     Double32_t ratios[ncells];
280     Int_t *dlist = recpoint->GetDigitsList();
281     Float_t *elist = recpoint->GetEnergiesList();
282     for (Int_t c = 0; c < ncells; ++c) {
283       AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(dlist[c]));
284       absIds[ncells_true] = digit->GetId();
285       ratios[ncells_true] = elist[c]/digit->GetAmplitude();
286       ++ncells_true;
287     }
288     
289     if (ncells_true < 1) {
290       AliWarning("Skipping cluster with no cells");
291       continue;
292     }
293     
294     // calculate new cluster position
295     TVector3 gpos;
296     recpoint->GetGlobalPosition(gpos);
297     Float_t g[3];
298     gpos.GetXYZ(g);
299     
300     AliVCluster *c = static_cast<AliVCluster*>(clus->New(nout++));
301     c->SetType(AliVCluster::kEMCALClusterv1);
302     c->SetE(recpoint->GetEnergy());
303     c->SetPosition(g);
304     c->SetNCells(ncells_true);
305     c->SetDispersion(recpoint->GetDispersion());
306     c->SetEmcCpvDistance(-1);            //not yet implemented
307     c->SetChi2(-1);                      //not yet implemented
308     c->SetTOF(recpoint->GetTime()) ;     //time-of-flight
309     c->SetNExMax(recpoint->GetNExMax()); //number of local maxima
310     Float_t elipAxis[2];
311     recpoint->GetElipsAxis(elipAxis);
312     c->SetM02(elipAxis[0]*elipAxis[0]) ;
313     c->SetM20(elipAxis[1]*elipAxis[1]) ;
314     if (fRecoUtils && fRecoUtils->IsBadChannelsRemovalSwitchedOn()) {
315       fRecoUtils->RecalculateClusterDistanceToBadChannel(geom, cells, c);
316     } else {
317       if (fPedestalData) 
318         recpoint->EvalDistanceToBadChannels(fPedestalData);
319       c->SetDistanceToBadChannel(recpoint->GetDistanceToBadTower()); 
320     }
321
322     if (esdobjects) {
323       AliESDCaloCluster *cesd = static_cast<AliESDCaloCluster*>(c);
324       cesd->SetCellsAbsId(absIds);
325       cesd->SetCellsAmplitudeFraction(ratios);
326     } else {
327       AliAODCaloCluster *caod = static_cast<AliAODCaloCluster*>(c);
328       caod->SetCellsAbsId(absIds);
329       caod->SetCellsAmplitudeFraction(ratios);
330     }
331   }
332  
333   AliESDEvent *esdevent = dynamic_cast<AliESDEvent*>(InputEvent());
334   if (!esdevent)
335     return;
336   if (!fRecoUtils)
337     return;
338
339   AliAnalysisManager::GetAnalysisManager()->LoadBranch("Tracks");
340   fRecoUtils->FindMatches(esdevent,clus);
341   
342   if (!esdobjects) {
343     Int_t Nclus = clus->GetEntries();
344     for(Int_t i=0; i < Nclus; ++i) {
345       AliAODCaloCluster *c = static_cast<AliAODCaloCluster*>(clus->At(i));
346       Int_t trackIndex = fRecoUtils->GetMatchedTrackIndex(i);
347       if(trackIndex >= 0) {
348         Float_t dR, dZ;
349         fRecoUtils->GetMatchedResiduals(i,dR, dZ);
350         c->AddTrackMatched(0x0); //esdevent->GetTrack(trackIndex));
351         c->SetTrackDistance(dR,dZ); // not implemented
352         c->SetEmcCpvDistance(dR);
353         c->SetChi2(dZ);
354         if(DebugLevel() > 1) 
355           AliInfo(Form("Matched Track index %d to new cluster %d\n",trackIndex,i));
356       }
357     }
358   } else {
359     Int_t Nclus = clus->GetEntries();
360     for(Int_t i=0; i < Nclus; ++i) {
361       AliESDCaloCluster *c = static_cast<AliESDCaloCluster*>(clus->At(i));
362       Int_t trackIndex = fRecoUtils->GetMatchedTrackIndex(i);
363       if(trackIndex >= 0) {
364         Float_t dR, dZ;
365         fRecoUtils->GetMatchedResiduals(i,dR, dZ);
366         c->SetTrackDistance(dR,dZ);
367         c->SetEmcCpvDistance(dR); //to be consistent with AODs
368         c->SetChi2(dZ);           //to be consistent with AODs
369         TArrayI tm(1,&trackIndex);
370         c->AddTracksMatched(tm);
371         if(DebugLevel() > 1) 
372           AliInfo(Form("Matched Track index %d to new cluster %d\n",trackIndex,i));
373       }
374     }
375   }
376 }
377
378 //________________________________________________________________________
379 void AliAnalysisTaskEMCALClusterizeFast::UpdateCells()
380 {
381   // Update cells in case re-calibration was done.
382
383   if (!fCalibData&&!fSubBackground)
384     return;
385
386   AliVCaloCells *cells = InputEvent()->GetEMCALCells();
387   Int_t ncells = cells->GetNumberOfCells();
388   Int_t ndigis = fDigitsArr->GetEntries();
389
390   if (ncells!=ndigis) {
391     cells->DeleteContainer();
392     cells->CreateContainer(ndigis);
393   }
394   for (Int_t idigit = 0; idigit < ndigis; ++idigit) {
395     AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(idigit));
396     Double_t cellAmplitude = digit->GetCalibAmp();
397     Short_t cellNumber = digit->GetId();
398     Double_t cellTime = digit->GetTime();
399     cells->SetCell(idigit, cellNumber, cellAmplitude, cellTime);
400   }
401 }
402
403 //________________________________________________________________________
404 void AliAnalysisTaskEMCALClusterizeFast::UpdateClusters()
405 {
406   // Update cells in case re-calibration was done.
407
408   if (!fAttachClusters)
409     return;
410
411   TClonesArray *clus = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("caloClusters"));
412   if (!clus)
413     clus = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("CaloClusters"));
414   if(!clus)
415     return;
416
417   Int_t nents = clus->GetEntries();
418   for (Int_t i=0;i<nents;++i) {
419     AliVCluster *c = static_cast<AliVCluster*>(clus->At(i));
420     if (!c)
421       continue;
422     if (c->IsEMCAL())
423       delete clus->RemoveAt(i);
424   }
425   clus->Compress();
426   RecPoints2Clusters(clus);
427 }
428
429 //________________________________________________________________________________________
430 void AliAnalysisTaskEMCALClusterizeFast::Init()
431 {
432   //Select clusterization/unfolding algorithm and set all the needed parameters
433
434   AliVEvent * event = InputEvent();
435   if (!event) {
436     AliWarning("Event not available!!!");
437     return;
438   }
439
440   if (event->GetRunNumber()==fRun)
441     return;
442   fRun = event->GetRunNumber();
443
444   if (fJustUnfold){
445     // init the unfolding afterburner 
446     delete fUnfolder;
447     fUnfolder = new AliEMCALAfterBurnerUF(fRecParam->GetW0(),fRecParam->GetLocMaxCut());
448     return;
449   }
450
451   AliEMCALGeometry *geometry = AliEMCALGeometry::GetInstance(fGeomName);
452   if (!geometry) {
453     AliFatal("Geometry not available!!!");
454     return;
455   }
456
457   if (!fGeomMatrixSet) {
458     if (fLoadGeomMatrices) {
459       for(Int_t mod=0; mod < geometry->GetNumberOfSuperModules(); ++mod) {
460         if(fGeomMatrix[mod]){
461           if(DebugLevel() > 2) 
462             fGeomMatrix[mod]->Print();
463           geometry->SetMisalMatrix(fGeomMatrix[mod],mod);  
464         }
465       }
466     } else { // get matrix from file (work around bug in aliroot)
467       for(Int_t mod=0; mod < geometry->GetEMCGeometry()->GetNumberOfSuperModules(); ++mod) {
468         const TGeoHMatrix *gm = 0;
469         AliESDEvent *esdevent = dynamic_cast<AliESDEvent*>(event);
470         if (esdevent) {
471           gm = esdevent->GetEMCALMatrix(mod);
472         } else {
473           AliAODHeader *aodheader = dynamic_cast<AliAODHeader*>(event->GetHeader());
474           if (aodheader) {
475             gm = aodheader->GetEMCALMatrix(mod);
476           }
477         }
478         if (gm) {
479           if(DebugLevel() > 2) 
480             gm->Print();
481           geometry->SetMisalMatrix(gm,mod);
482         }
483       }
484     }
485     fGeomMatrixSet=kTRUE;
486   }
487   
488   // setup digit array if needed
489   if (!fDigitsArr) {
490     fDigitsArr = new TClonesArray("AliEMCALDigit", 1000);
491     fDigitsArr->SetOwner(1);
492   }
493
494   // then setup clusterizer
495   delete fClusterizer;
496   if     (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv1)
497     fClusterizer = new AliEMCALClusterizerv1(geometry);
498   else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerNxN) 
499     fClusterizer = new AliEMCALClusterizerNxN(geometry);
500   else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv2) 
501     fClusterizer = new AliEMCALClusterizerv2(geometry);
502   else if(fRecParam->GetClusterizerFlag() > AliEMCALRecParam::kClusterizerNxN) {
503    AliEMCALClusterizerNxN *clusterizer = new AliEMCALClusterizerNxN(geometry);
504    clusterizer->SetNRowDiff(2);
505    clusterizer->SetNColDiff(2);
506    fClusterizer = clusterizer;
507   } else{
508     AliFatal(Form("Clusterizer < %d > not available", fRecParam->GetClusterizerFlag()));
509   }
510   fClusterizer->InitParameters(fRecParam);
511
512   if ((!fCalibData&&fLoadCalib) || (!fPedestalData&&fLoadPed)) {
513     AliCDBManager *cdb = AliCDBManager::Instance();
514     if (!cdb->IsDefaultStorageSet() && !fOCDBpath.IsNull())
515       cdb->SetDefaultStorage(fOCDBpath);
516     if (fRun!=cdb->GetRun())
517       cdb->SetRun(fRun);
518   }
519   if (!fCalibData&&fLoadCalib&&fRun>0) {
520     AliCDBEntry *entry = static_cast<AliCDBEntry*>(AliCDBManager::Instance()->Get("EMCAL/Calib/Data"));
521     if (entry) 
522       fCalibData =  static_cast<AliEMCALCalibData*>(entry->GetObject());
523     if (!fCalibData)
524       AliFatal("Calibration parameters not found in CDB!");
525   }
526   if (!fPedestalData&&fLoadPed&&fRun>0) {
527     AliCDBEntry *entry = static_cast<AliCDBEntry*>(AliCDBManager::Instance()->Get("EMCAL/Calib/Pedestals"));
528     if (entry) 
529       fPedestalData =  static_cast<AliCaloCalibPedestal*>(entry->GetObject());
530   }
531   if (fCalibData) {
532     fClusterizer->SetInputCalibrated(kFALSE);   
533     fClusterizer->SetCalibrationParameters(fCalibData);
534   } else {
535     fClusterizer->SetInputCalibrated(kTRUE);   
536   }
537   fClusterizer->SetCaloCalibPedestal(fPedestalData);
538   fClusterizer->SetJustClusters(kTRUE);
539   fClusterizer->SetDigitsArr(fDigitsArr);
540   fClusterizer->SetOutput(0);
541   fClusterArr = const_cast<TObjArray *>(fClusterizer->GetRecPoints());
542 }