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