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