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