649954fb3cad679766b15225ce2ce249625bb49a
[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.,AliESDRun.,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     if (fCalibData) {
183       fClusterizer->SetInputCalibrated(kFALSE);   
184       fClusterizer->SetCalibrationParameters(fCalibData);
185     }
186   }
187 }
188
189 //________________________________________________________________________
190 void AliAnalysisTaskEMCALClusterizeFast::FillDigitsArray()
191 {
192   // Fill digits from cells.
193
194   fDigitsArr->Clear("C");
195   AliVCaloCells *cells = InputEvent()->GetEMCALCells();
196   Double_t avgE = 0; // for background subtraction
197   Int_t ncells = cells->GetNumberOfCells();
198   for (Int_t icell = 0, idigit = 0; icell < ncells; ++icell) {
199     Double_t cellAmplitude=0, cellTime=0;
200     Short_t cellNumber=0;
201     if (cells->GetCell(icell, cellNumber, cellAmplitude, cellTime) != kTRUE)
202       break;
203     AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->New(idigit));
204     digit->SetId(cellNumber);
205     digit->SetTime(cellTime);
206     digit->SetTimeR(cellTime);
207     digit->SetIndexInList(idigit);
208     digit->SetType(AliEMCALDigit::kHG);
209     if (fRecalibOnly||fSubBackground) {
210       Double_t energy = fClusterizer->Calibrate(cellAmplitude,cellTime,cellNumber);
211       digit->SetAmplitude(energy);
212       avgE += energy;
213     } else {
214       digit->SetAmplitude(cellAmplitude);
215     }
216     idigit++;
217   }
218
219   if (fSubBackground) {
220     avgE /= AliEMCALGeometry::GetInstance(fGeomName)->GetNumberOfSuperModules()*48*24;
221     Int_t ndigis = fDigitsArr->GetEntries();
222     for (Int_t i = 0; i < ndigis; ++i) {
223       AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(i));
224       Double_t energy = digit->GetAmplitude() - avgE;
225       if (energy<=0.001) {
226         //fDigitsArr->RemoveAt(i);
227         digit->SetAmplitude(0);
228       } else {
229         digit->SetAmplitude(energy);
230       }
231     }
232   }
233   fDigitsArr->Compress();
234   fDigitsArr->Sort();
235 }
236
237 //________________________________________________________________________________________
238 void AliAnalysisTaskEMCALClusterizeFast::RecPoints2Clusters(TClonesArray *clus)
239 {
240   // Cluster energy, global position, cells and their amplitude fractions are restored.
241
242   Bool_t esdobjects = 0;
243   if (strcmp(clus->GetClass()->GetName(),"AliESDCaloCluster")==0)
244     esdobjects = 1;
245
246   Int_t Ncls = fClusterArr->GetEntriesFast();
247   for(Int_t i=0, nout=clus->GetEntries(); i < Ncls; ++i) {
248     AliEMCALRecPoint *recpoint = static_cast<AliEMCALRecPoint*>(fClusterArr->At(i));
249     Int_t ncells_true = 0;
250     const Int_t ncells = recpoint->GetMultiplicity();
251     UShort_t   absIds[ncells];  
252     Double32_t ratios[ncells];
253     Int_t *dlist = recpoint->GetDigitsList();
254     Float_t *elist = recpoint->GetEnergiesList();
255     for (Int_t c = 0; c < ncells; ++c) {
256       AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(dlist[c]));
257       absIds[ncells_true] = digit->GetId();
258       ratios[ncells_true] = elist[c]/digit->GetAmplitude();
259       if (ratios[ncells_true] < 0.001) 
260         continue;
261       ++ncells_true;
262     }
263     
264     if (ncells_true < 1) {
265       AliWarning("Skipping cluster with no cells");
266       continue;
267     }
268     
269     // calculate new cluster position
270     TVector3 gpos;
271     recpoint->GetGlobalPosition(gpos);
272     Float_t g[3];
273     gpos.GetXYZ(g);
274     
275     AliVCluster *c = static_cast<AliVCluster*>(clus->New(nout++));
276     c->SetType(AliVCluster::kEMCALClusterv1);
277     c->SetE(recpoint->GetEnergy());
278     c->SetPosition(g);
279     c->SetNCells(ncells_true);
280     c->SetDispersion(recpoint->GetDispersion());
281     c->SetEmcCpvDistance(-1);            //not yet implemented
282     c->SetChi2(-1);                      //not yet implemented
283     c->SetTOF(recpoint->GetTime()) ;     //time-of-flight
284     c->SetNExMax(recpoint->GetNExMax()); //number of local maxima
285     Float_t elipAxis[2];
286     recpoint->GetElipsAxis(elipAxis);
287     c->SetM02(elipAxis[0]*elipAxis[0]) ;
288     c->SetM20(elipAxis[1]*elipAxis[1]) ;
289     c->SetDistanceToBadChannel(recpoint->GetDistanceToBadTower()); 
290     if (esdobjects) {
291       AliESDCaloCluster *cesd = static_cast<AliESDCaloCluster*>(c);
292       cesd->SetCellsAbsId(absIds);
293       cesd->SetCellsAmplitudeFraction(ratios);
294     } else {
295       AliAODCaloCluster *caod = static_cast<AliAODCaloCluster*>(c);
296       caod->SetCellsAbsId(absIds);
297       caod->SetCellsAmplitudeFraction(ratios);
298     }
299   }
300  
301   AliESDEvent *esdevent = dynamic_cast<AliESDEvent*>(InputEvent());
302   if (!esdevent)
303     return;
304   if (!fRecoUtils)
305     return;
306
307   AliAnalysisManager::GetAnalysisManager()->LoadBranch("Tracks");
308   fRecoUtils->FindMatches(esdevent,clus);
309   
310   if (!esdobjects) {
311     Int_t Nclus = clus->GetEntries();
312     for(Int_t i=0; i < Nclus; ++i) {
313       AliAODCaloCluster *c = static_cast<AliAODCaloCluster*>(clus->At(i));
314       Int_t trackIndex = fRecoUtils->GetMatchedTrackIndex(i);
315       if(trackIndex >= 0) {
316         Float_t dR, dZ;
317         fRecoUtils->GetMatchedResiduals(i,dR, dZ);
318         c->AddTrackMatched(0x0); //esdevent->GetTrack(trackIndex));
319         c->SetTrackDistance(dR,dZ); // not implemented
320         c->SetEmcCpvDistance(dR);
321         c->SetChi2(dZ);
322         if(DebugLevel() > 1) 
323           AliInfo(Form("Matched Track index %d to new cluster %d\n",trackIndex,i));
324       }
325     }
326   } else {
327     Int_t Nclus = clus->GetEntries();
328     for(Int_t i=0; i < Nclus; ++i) {
329       AliESDCaloCluster *c = static_cast<AliESDCaloCluster*>(clus->At(i));
330       Int_t trackIndex = fRecoUtils->GetMatchedTrackIndex(i);
331       if(trackIndex >= 0) {
332         Float_t dR, dZ;
333         fRecoUtils->GetMatchedResiduals(i,dR, dZ);
334         c->SetTrackDistance(dR,dZ);
335         c->SetEmcCpvDistance(dR); //to be consistent with AODs
336         c->SetChi2(dZ);           //to be consistent with AODs
337         TArrayI tm(1,&trackIndex);
338         c->AddTracksMatched(tm);
339         if(DebugLevel() > 1) 
340           AliInfo(Form("Matched Track index %d to new cluster %d\n",trackIndex,i));
341       }
342     }
343   }
344 }
345
346 //________________________________________________________________________
347 void AliAnalysisTaskEMCALClusterizeFast::UpdateCells()
348 {
349   // Update cells in case re-calibration was done.
350
351   if (!fCalibData&&!fSubBackground)
352     return;
353
354   AliVCaloCells *cells = InputEvent()->GetEMCALCells();
355   Int_t ncells = cells->GetNumberOfCells();
356   Int_t ndigis = fDigitsArr->GetEntries();
357
358   if (ncells!=ndigis) {
359     cells->DeleteContainer();
360     cells->CreateContainer(ndigis);
361   }
362   for (Int_t idigit = 0; idigit < ndigis; ++idigit) {
363     AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(idigit));
364     Double_t cellAmplitude = digit->GetCalibAmp();
365     Short_t cellNumber = digit->GetId();
366     Double_t cellTime = digit->GetTime();
367     cells->SetCell(idigit, cellNumber, cellAmplitude, cellTime);
368   }
369 }
370
371 //________________________________________________________________________
372 void AliAnalysisTaskEMCALClusterizeFast::UpdateClusters()
373 {
374   // Update cells in case re-calibration was done.
375
376   if (!fAttachClusters)
377     return;
378
379   TClonesArray *clus = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("caloClusters"));
380   if (!clus)
381     clus = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("CaloClusters"));
382   if(!clus)
383     return;
384
385   Int_t nents = clus->GetEntries();
386   for (Int_t i=0;i<nents;++i) {
387     AliVCluster *c = static_cast<AliVCluster*>(clus->At(i));
388     if (!c)
389       continue;
390     if (c->IsEMCAL())
391       clus->RemoveAt(i);
392   }
393   clus->Compress();
394   RecPoints2Clusters(clus);
395 }
396
397 //________________________________________________________________________________________
398 void AliAnalysisTaskEMCALClusterizeFast::Init()
399 {
400   //Select clusterization/unfolding algorithm and set all the needed parameters
401
402   AliVEvent * event = InputEvent();
403   if (!event) {
404     AliWarning("Event not available!!!");
405     return;
406   }
407
408   if (event->GetRunNumber()==fRun)
409     return;
410   fRun = event->GetRunNumber();
411
412   if (fJustUnfold){
413     // init the unfolding afterburner 
414     delete fUnfolder;
415     fUnfolder = new AliEMCALAfterBurnerUF(fRecParam->GetW0(),fRecParam->GetLocMaxCut());
416     return;
417   }
418
419   AliEMCALGeometry *geometry = AliEMCALGeometry::GetInstance(fGeomName);
420   if (!geometry) {
421     AliFatal("Geometry not available!!!");
422     return;
423   }
424
425   if (!fGeomMatrixSet) {
426     if (fLoadGeomMatrices) {
427       for(Int_t mod=0; mod < geometry->GetNumberOfSuperModules(); ++mod) {
428         if(fGeomMatrix[mod]){
429           if(DebugLevel() > 2) 
430             fGeomMatrix[mod]->Print();
431           geometry->SetMisalMatrix(fGeomMatrix[mod],mod);  
432         }
433       }
434     } else { // get matrix from file (work around bug in aliroot)
435       for(Int_t mod=0; mod < geometry->GetEMCGeometry()->GetNumberOfSuperModules(); ++mod) {
436         const TGeoHMatrix *gm = 0;
437         AliESDEvent *esdevent = dynamic_cast<AliESDEvent*>(event);
438         if (esdevent) {
439           gm = esdevent->GetEMCALMatrix(mod);
440         } else {
441           AliAODHeader *aodheader = dynamic_cast<AliAODHeader*>(event->GetHeader());
442           if (aodheader) {
443             gm = aodheader->GetEMCALMatrix(mod);
444           }
445         }
446         if (gm) {
447           if(DebugLevel() > 2) 
448             gm->Print();
449           geometry->SetMisalMatrix(gm,mod);
450         }
451       }
452     }
453     fGeomMatrixSet=kTRUE;
454   }
455   
456   // setup digit array if needed
457   if (!fDigitsArr) {
458     fDigitsArr = new TClonesArray("AliEMCALDigit", 1000);
459     fDigitsArr->SetOwner(1);
460   }
461
462   // then setup clusterizer
463   delete fClusterizer;
464   if     (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv1)
465     fClusterizer = new AliEMCALClusterizerv1(geometry);
466   else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerNxN) 
467     fClusterizer = new AliEMCALClusterizerNxN(geometry);
468   else if(fRecParam->GetClusterizerFlag() > AliEMCALRecParam::kClusterizerNxN) {
469    AliEMCALClusterizerNxN *clusterizer = new AliEMCALClusterizerNxN(geometry);
470    clusterizer->SetNRowDiff(2);
471    clusterizer->SetNColDiff(2);
472    fClusterizer = clusterizer;
473   } else{
474     AliFatal(Form("Clusterizer < %d > not available", fRecParam->GetClusterizerFlag()));
475   }
476   fClusterizer->InitParameters(fRecParam);
477
478   if ((!fCalibData&&fLoadCalib) || (!fPedestalData&&fLoadPed)) {
479     AliCDBManager *cdb = AliCDBManager::Instance();
480     if (!cdb->IsDefaultStorageSet() && !fOCDBpath.IsNull())
481       cdb->SetDefaultStorage(fOCDBpath);
482     if (fRun!=cdb->GetRun())
483       cdb->SetRun(fRun);
484   }
485   if (!fCalibData&&fLoadCalib) {
486     AliCDBEntry *entry = static_cast<AliCDBEntry*>(AliCDBManager::Instance()->Get("EMCAL/Calib/Data"));
487     if (entry) 
488       fCalibData =  static_cast<AliEMCALCalibData*>(entry->GetObject());
489     if (!fCalibData)
490       AliFatal("Calibration parameters not found in CDB!");
491   }
492   if (!fPedestalData&&fLoadPed) {
493     AliCDBEntry *entry = static_cast<AliCDBEntry*>(AliCDBManager::Instance()->Get("EMCAL/Calib/Pedestals"));
494     if (entry) 
495       fPedestalData =  static_cast<AliCaloCalibPedestal*>(entry->GetObject());
496   }
497   if (fCalibData) {
498     fClusterizer->SetInputCalibrated(kFALSE);   
499     fClusterizer->SetCalibrationParameters(fCalibData);
500   } else {
501     fClusterizer->SetInputCalibrated(kTRUE);   
502   }
503   fClusterizer->SetCaloCalibPedestal(fPedestalData);
504   fClusterizer->SetJustClusters(kTRUE);
505   fClusterizer->SetDigitsArr(fDigitsArr);
506   fClusterizer->SetOutput(0);
507   fClusterArr = const_cast<TObjArray *>(fClusterizer->GetRecPoints());
508 }