Small bug fix (A.Marin)
[u/mrichter/AliRoot.git] / PWG4 / CaloCalib / AliAnalysisTaskEMCALClusterizeFast.cxx
CommitLineData
2f7259cf 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
44ClassImp(AliAnalysisTaskEMCALClusterizeFast)
45
46//________________________________________________________________________
2431b20f 47AliAnalysisTaskEMCALClusterizeFast::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{
68 // Constructor
69}
70
71//________________________________________________________________________
2f7259cf 72AliAnalysisTaskEMCALClusterizeFast::AliAnalysisTaskEMCALClusterizeFast(const char *name)
73 : AliAnalysisTaskSE(name),
74 fRun(0),
75 fDigitsArr(0),
76 fClusterArr(0),
77 fRecParam(new AliEMCALRecParam),
78 fClusterizer(0),
79 fUnfolder(0),
80 fJustUnfold(kFALSE),
81 fGeomName("EMCAL_FIRSTYEARV1"),
82 fGeomMatrixSet(kFALSE),
83 fLoadGeomMatrices(kFALSE),
84 fOCDBpath(),
85 fCalibData(0),
86 fPedestalData(0),
87 fOutputAODBranch(0),
88 fOutputAODBrName(),
10d33986 89 fRecoUtils(0),
90 fLoadCalib(0),
91 fLoadPed(0)
2f7259cf 92{
93 // Constructor
94
10922c99 95 fBranchNames = "ESD:AliESDHeader.,EMCALCells. AOD:header,emcalCells";
2f7259cf 96 for(Int_t i = 0; i < 10; ++i)
97 fGeomMatrix[i] = 0;
98}
99
100//________________________________________________________________________
101AliAnalysisTaskEMCALClusterizeFast::~AliAnalysisTaskEMCALClusterizeFast()
102{
103 // Destructor.
104
105 delete fDigitsArr;
106 delete fClusterizer;
107 delete fUnfolder;
108 delete fRecoUtils;
109}
110
111//-------------------------------------------------------------------
112void AliAnalysisTaskEMCALClusterizeFast::UserCreateOutputObjects()
113{
114 // Create output objects.
115
116 if (!fOutputAODBrName.IsNull()) {
117 fOutputAODBranch = new TClonesArray("AliAODCaloCluster", 0);
118 fOutputAODBranch->SetName(fOutputAODBrName);
119 AddAODBranch("TClonesArray", &fOutputAODBranch);
120 AliInfo(Form("Created Branch: %s",fOutputAODBrName.Data()));
121 }
122}
123
124//________________________________________________________________________
125void AliAnalysisTaskEMCALClusterizeFast::UserExec(Option_t *)
126{
127 // Main loop, called for each event
128
129 // remove the contents of output list set in the previous event
130 if (fOutputAODBranch)
131 fOutputAODBranch->Clear("C");
132
133 AliESDEvent *esdevent = dynamic_cast<AliESDEvent*>(InputEvent());
134 AliAODEvent *aodevent = dynamic_cast<AliAODEvent*>(InputEvent());
135
136 if (!esdevent&&!aodevent) {
137 Error("UserExec","Event not available");
138 return;
139 }
140
141 LoadBranches();
142
143 Init();
144
145 if (fJustUnfold) {
146 AliWarning("Unfolding not implemented");
147 } else {
148 if (esdevent)
149 FillDigitsArray(esdevent);
150 else
151 FillDigitsArray(aodevent);
152 fClusterizer->Digits2Clusters("");
2431b20f 153 if (esdevent && fRecoUtils)
2f7259cf 154 fRecoUtils->FindMatches(esdevent,fClusterArr);
2431b20f 155 if (fOutputAODBranch) {
2f7259cf 156 RecPoints2Clusters();
2431b20f 157 }
158 if (esdevent) {
159 UpdateCells(esdevent);
160 } else {
161 UpdateCells(aodevent);
162 }
163 }
164}
165
166//________________________________________________________________________
167void AliAnalysisTaskEMCALClusterizeFast::UpdateCells(AliAODEvent *event)
168{
7099ec08 169 // Update cells in case re-calibration was done.
170
2431b20f 171 if (!fCalibData)
172 return;
173
174 AliAODCaloCells *cells = event->GetEMCALCells();
175 Int_t ncells = cells->GetNumberOfCells();
176 for (Int_t icell = 0, idigit = 0; icell < ncells; ++icell) {
177 Double_t cellAmplitude=0, cellTime=0;
178 Short_t cellNumber=0;
179 if (cells->GetCell(icell, cellNumber, cellAmplitude, cellTime) != kTRUE)
180 break;
181 AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(idigit));
182 cellAmplitude = digit->GetCalibAmp();
183 cells->SetCell(icell, cellNumber, cellAmplitude, cellTime);
184 idigit++;
185 }
186 cells->Sort();
187}
188
189//________________________________________________________________________
190void AliAnalysisTaskEMCALClusterizeFast::UpdateCells(AliESDEvent *event)
191{
7099ec08 192 // Update cells in case re-calibration was done.
193
2431b20f 194 if (!fCalibData)
195 return;
196
197 AliESDCaloCells *cells = event->GetEMCALCells();
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->At(idigit));
205 cellAmplitude = digit->GetCalibAmp();
206 cells->SetCell(icell, cellNumber, cellAmplitude, cellTime);
207 idigit++;
2f7259cf 208 }
2431b20f 209 cells->Sort();
2f7259cf 210}
211
212//________________________________________________________________________
213void AliAnalysisTaskEMCALClusterizeFast::FillDigitsArray(AliAODEvent *event)
214{
215 // Fill digits from cells.
216
217 fDigitsArr->Clear("C");
218 AliAODCaloCells *cells = event->GetEMCALCells();
219 Int_t ncells = cells->GetNumberOfCells();
220 if (ncells>fDigitsArr->GetSize())
221 fDigitsArr->Expand(2*ncells);
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 = (AliEMCALDigit*) fDigitsArr->New(idigit);
228 digit->SetId(cellNumber);
229 digit->SetAmplitude(cellAmplitude);
230 digit->SetTime(cellTime);
231 digit->SetTimeR(cellTime);
232 digit->SetIndexInList(idigit);
233 digit->SetType(AliEMCALDigit::kHG);
234 idigit++;
235 }
236}
237
238//________________________________________________________________________
239void AliAnalysisTaskEMCALClusterizeFast::FillDigitsArray(AliESDEvent *event)
240{
241 // Fill digits from cells.
242
243 fDigitsArr->Clear("C");
244 AliESDCaloCells *cells = event->GetEMCALCells();
245 Int_t ncells = cells->GetNumberOfCells();
246 if (ncells>fDigitsArr->GetSize())
247 fDigitsArr->Expand(2*ncells);
248 for (Int_t icell = 0, idigit = 0; icell < ncells; ++icell) {
249 Double_t cellAmplitude=0, cellTime=0;
250 Short_t cellNumber=0;
251 if (cells->GetCell(icell, cellNumber, cellAmplitude, cellTime) != kTRUE)
252 break;
253 AliEMCALDigit *digit = (AliEMCALDigit*) fDigitsArr->New(idigit);
254 digit->SetId(cellNumber);
255 digit->SetAmplitude(cellAmplitude);
256 digit->SetTime(cellTime);
257 digit->SetTimeR(cellTime);
258 digit->SetIndexInList(idigit);
259 digit->SetType(AliEMCALDigit::kHG);
260 idigit++;
261 }
262}
263
264//________________________________________________________________________________________
265void AliAnalysisTaskEMCALClusterizeFast::RecPoints2Clusters()
266{
267 // Cluster energy, global position, cells and their amplitude fractions are restored.
268
269 AliESDEvent *esdevent = dynamic_cast<AliESDEvent*>(InputEvent());
270
271 Int_t Ncls = fClusterArr->GetEntriesFast();
272 for(Int_t i=0, nout=0; i < Ncls; ++i) {
273 AliEMCALRecPoint *recpoint = static_cast<AliEMCALRecPoint*>(fClusterArr->At(i));
274 Int_t ncells_true = 0;
275 const Int_t ncells = recpoint->GetMultiplicity();
276 UShort_t absIds[ncells];
277 Double32_t ratios[ncells];
278 Int_t *dlist = recpoint->GetDigitsList();
279 Float_t *elist = recpoint->GetEnergiesList();
2f7259cf 280 for (Int_t c = 0; c < ncells; ++c) {
281 AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(dlist[c]));
282 absIds[ncells_true] = digit->GetId();
283 ratios[ncells_true] = elist[c]/digit->GetAmplitude();
284 if (ratios[ncells_true] > 0.001)
285 ++ncells_true;
286 }
287
288 if (ncells_true < 1) {
289 AliWarning("Skipping cluster with no cells");
290 continue;
291 }
292
293 // calculate new cluster position
294 TVector3 gpos;
295 Float_t g[3];
296
297 recpoint->EvalGlobalPosition(fRecParam->GetW0(), fDigitsArr);
298 recpoint->GetGlobalPosition(gpos);
299 gpos.GetXYZ(g);
300
301 AliAODCaloCluster *clus = static_cast<AliAODCaloCluster*>(fOutputAODBranch->New(nout++));
302 clus->SetType(AliVCluster::kEMCALClusterv1);
303 clus->SetE(recpoint->GetEnergy());
304 clus->SetPosition(g);
305 clus->SetNCells(ncells_true);
306 clus->SetCellsAbsId(absIds);
307 clus->SetCellsAmplitudeFraction(ratios);
308 clus->SetDispersion(recpoint->GetDispersion());
309 clus->SetChi2(-1); //not yet implemented
310 clus->SetTOF(recpoint->GetTime()) ; //time-of-flight
311 clus->SetNExMax(recpoint->GetNExMax()); //number of local maxima
312 Float_t elipAxis[2];
313 recpoint->GetElipsAxis(elipAxis);
314 clus->SetM02(elipAxis[0]*elipAxis[0]) ;
315 clus->SetM20(elipAxis[1]*elipAxis[1]) ;
316 clus->SetDistToBadChannel(recpoint->GetDistanceToBadTower());
317 if (esdevent && fRecoUtils) {
318 Int_t trackIndex = fRecoUtils->GetMatchedTrackIndex(i);
319 if(trackIndex >= 0) {
320 clus->AddTrackMatched(esdevent->GetTrack(trackIndex));
321 if(DebugLevel() > 1)
322 AliInfo(Form("Matched Track index %d to new cluster %d\n",trackIndex,i));
323 }
324 }
325 }
326}
327
328//________________________________________________________________________________________
329void AliAnalysisTaskEMCALClusterizeFast::Init()
330{
331 //Select clusterization/unfolding algorithm and set all the needed parameters
332
333 AliVEvent * event = InputEvent();
334 if (!event) {
335 AliWarning("Event not available!!!");
336 return;
337 }
338
339 if (event->GetRunNumber()==fRun)
340 return;
341 fRun = event->GetRunNumber();
342
343 if (fJustUnfold){
344 // init the unfolding afterburner
345 delete fUnfolder;
346 fUnfolder = new AliEMCALAfterBurnerUF(fRecParam->GetW0(),fRecParam->GetLocMaxCut());
347 return;
348 }
349
350 AliCDBManager *cdb = AliCDBManager::Instance();
351 if (!cdb->IsDefaultStorageSet() && !fOCDBpath.IsNull())
352 cdb->SetDefaultStorage(fOCDBpath);
353 if (fRun!=cdb->GetRun())
354 cdb->SetRun(fRun);
355
356 AliEMCALGeometry *geometry = AliEMCALGeometry::GetInstance(fGeomName);
357 if (!geometry) {
358 AliFatal("Geometry not available!!!");
359 return;
360 }
361
362 if (!fGeomMatrixSet) {
363 if (fLoadGeomMatrices) {
364 for(Int_t mod=0; mod < (geometry->GetEMCGeometry())->GetNumberOfSuperModules(); ++mod) {
365 if(fGeomMatrix[mod]){
366 if(DebugLevel() > 2)
367 fGeomMatrix[mod]->Print();
368 geometry->SetMisalMatrix(fGeomMatrix[mod],mod);
369 }
370 }
371 } else {
372 for(Int_t mod=0; mod < geometry->GetEMCGeometry()->GetNumberOfSuperModules(); ++mod) {
373 if(event->GetEMCALMatrix(mod)) {
374 if(DebugLevel() > 2)
375 event->GetEMCALMatrix(mod)->Print();
376 geometry->SetMisalMatrix(event->GetEMCALMatrix(mod),mod);
377 }
378 }
379 }
380 fGeomMatrixSet=kTRUE;
381 }
382
383 // setup digit array if needed
384 if (!fDigitsArr) {
385 fDigitsArr = new TClonesArray("AliEMCALDigit", 1000);
386 fDigitsArr->SetOwner(1);
387 }
388
389 // then setup clusterizer
390 delete fClusterizer;
391 if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv1)
392 fClusterizer = new AliEMCALClusterizerv1(geometry);
393 else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerNxN)
394 fClusterizer = new AliEMCALClusterizerNxN(geometry);
395 else if(fRecParam->GetClusterizerFlag() > AliEMCALRecParam::kClusterizerNxN) {
396 AliEMCALClusterizerNxN *clusterizer = new AliEMCALClusterizerNxN(geometry);
397 clusterizer->SetNRowDiff(2);
398 clusterizer->SetNColDiff(2);
399 fClusterizer = clusterizer;
400 } else{
401 AliFatal(Form("Clusterizer < %d > not available", fRecParam->GetClusterizerFlag()));
402 }
403 fClusterizer->InitParameters(fRecParam);
10d33986 404 if (!fCalibData&&fLoadCalib) {
2f7259cf 405 AliCDBEntry *entry = static_cast<AliCDBEntry*>(AliCDBManager::Instance()->Get("EMCAL/Calib/Data"));
406 if (entry)
407 fCalibData = static_cast<AliEMCALCalibData*>(entry->GetObject());
408 if (!fCalibData)
409 AliFatal("Calibration parameters not found in CDB!");
410 }
10d33986 411 if (!fPedestalData&&fLoadPed) {
2f7259cf 412 AliCDBEntry *entry = static_cast<AliCDBEntry*>(AliCDBManager::Instance()->Get("EMCAL/Calib/Pedestals"));
413 if (entry)
414 fPedestalData = static_cast<AliCaloCalibPedestal*>(entry->GetObject());
415 }
10d33986 416 if (fCalibData) {
417 fClusterizer->SetInputCalibrated(kFALSE);
418 fClusterizer->SetCalibrationParameters(fCalibData);
419 fClusterizer->SetCaloCalibPedestal(fPedestalData);
420 } else {
421 fClusterizer->SetInputCalibrated(kTRUE);
422 }
423 fClusterizer->SetJustClusters(kTRUE);
2f7259cf 424 fClusterizer->SetDigitsArr(fDigitsArr);
425 fClusterizer->SetOutput(0);
426 fClusterArr = const_cast<TObjArray *>(fClusterizer->GetRecPoints());
427}