update cells in case calibration was done
[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{
169 if (!fCalibData)
170 return;
171
172 AliAODCaloCells *cells = event->GetEMCALCells();
173 Int_t ncells = cells->GetNumberOfCells();
174 for (Int_t icell = 0, idigit = 0; icell < ncells; ++icell) {
175 Double_t cellAmplitude=0, cellTime=0;
176 Short_t cellNumber=0;
177 if (cells->GetCell(icell, cellNumber, cellAmplitude, cellTime) != kTRUE)
178 break;
179 AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(idigit));
180 cellAmplitude = digit->GetCalibAmp();
181 cells->SetCell(icell, cellNumber, cellAmplitude, cellTime);
182 idigit++;
183 }
184 cells->Sort();
185}
186
187//________________________________________________________________________
188void AliAnalysisTaskEMCALClusterizeFast::UpdateCells(AliESDEvent *event)
189{
190 if (!fCalibData)
191 return;
192
193 AliESDCaloCells *cells = event->GetEMCALCells();
194 Int_t ncells = cells->GetNumberOfCells();
195 for (Int_t icell = 0, idigit = 0; icell < ncells; ++icell) {
196 Double_t cellAmplitude=0, cellTime=0;
197 Short_t cellNumber=0;
198 if (cells->GetCell(icell, cellNumber, cellAmplitude, cellTime) != kTRUE)
199 break;
200 AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(idigit));
201 cellAmplitude = digit->GetCalibAmp();
202 cells->SetCell(icell, cellNumber, cellAmplitude, cellTime);
203 idigit++;
2f7259cf 204 }
2431b20f 205 cells->Sort();
2f7259cf 206}
207
208//________________________________________________________________________
209void AliAnalysisTaskEMCALClusterizeFast::FillDigitsArray(AliAODEvent *event)
210{
211 // Fill digits from cells.
212
213 fDigitsArr->Clear("C");
214 AliAODCaloCells *cells = event->GetEMCALCells();
215 Int_t ncells = cells->GetNumberOfCells();
216 if (ncells>fDigitsArr->GetSize())
217 fDigitsArr->Expand(2*ncells);
218 for (Int_t icell = 0, idigit = 0; icell < ncells; ++icell) {
219 Double_t cellAmplitude=0, cellTime=0;
220 Short_t cellNumber=0;
221 if (cells->GetCell(icell, cellNumber, cellAmplitude, cellTime) != kTRUE)
222 break;
223 AliEMCALDigit *digit = (AliEMCALDigit*) fDigitsArr->New(idigit);
224 digit->SetId(cellNumber);
225 digit->SetAmplitude(cellAmplitude);
226 digit->SetTime(cellTime);
227 digit->SetTimeR(cellTime);
228 digit->SetIndexInList(idigit);
229 digit->SetType(AliEMCALDigit::kHG);
230 idigit++;
231 }
232}
233
234//________________________________________________________________________
235void AliAnalysisTaskEMCALClusterizeFast::FillDigitsArray(AliESDEvent *event)
236{
237 // Fill digits from cells.
238
239 fDigitsArr->Clear("C");
240 AliESDCaloCells *cells = event->GetEMCALCells();
241 Int_t ncells = cells->GetNumberOfCells();
242 if (ncells>fDigitsArr->GetSize())
243 fDigitsArr->Expand(2*ncells);
244 for (Int_t icell = 0, idigit = 0; icell < ncells; ++icell) {
245 Double_t cellAmplitude=0, cellTime=0;
246 Short_t cellNumber=0;
247 if (cells->GetCell(icell, cellNumber, cellAmplitude, cellTime) != kTRUE)
248 break;
249 AliEMCALDigit *digit = (AliEMCALDigit*) fDigitsArr->New(idigit);
250 digit->SetId(cellNumber);
251 digit->SetAmplitude(cellAmplitude);
252 digit->SetTime(cellTime);
253 digit->SetTimeR(cellTime);
254 digit->SetIndexInList(idigit);
255 digit->SetType(AliEMCALDigit::kHG);
256 idigit++;
257 }
258}
259
260//________________________________________________________________________________________
261void AliAnalysisTaskEMCALClusterizeFast::RecPoints2Clusters()
262{
263 // Cluster energy, global position, cells and their amplitude fractions are restored.
264
265 AliESDEvent *esdevent = dynamic_cast<AliESDEvent*>(InputEvent());
266
267 Int_t Ncls = fClusterArr->GetEntriesFast();
268 for(Int_t i=0, nout=0; i < Ncls; ++i) {
269 AliEMCALRecPoint *recpoint = static_cast<AliEMCALRecPoint*>(fClusterArr->At(i));
270 Int_t ncells_true = 0;
271 const Int_t ncells = recpoint->GetMultiplicity();
272 UShort_t absIds[ncells];
273 Double32_t ratios[ncells];
274 Int_t *dlist = recpoint->GetDigitsList();
275 Float_t *elist = recpoint->GetEnergiesList();
2f7259cf 276 for (Int_t c = 0; c < ncells; ++c) {
277 AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(dlist[c]));
278 absIds[ncells_true] = digit->GetId();
279 ratios[ncells_true] = elist[c]/digit->GetAmplitude();
280 if (ratios[ncells_true] > 0.001)
281 ++ncells_true;
282 }
283
284 if (ncells_true < 1) {
285 AliWarning("Skipping cluster with no cells");
286 continue;
287 }
288
289 // calculate new cluster position
290 TVector3 gpos;
291 Float_t g[3];
292
293 recpoint->EvalGlobalPosition(fRecParam->GetW0(), fDigitsArr);
294 recpoint->GetGlobalPosition(gpos);
295 gpos.GetXYZ(g);
296
297 AliAODCaloCluster *clus = static_cast<AliAODCaloCluster*>(fOutputAODBranch->New(nout++));
298 clus->SetType(AliVCluster::kEMCALClusterv1);
299 clus->SetE(recpoint->GetEnergy());
300 clus->SetPosition(g);
301 clus->SetNCells(ncells_true);
302 clus->SetCellsAbsId(absIds);
303 clus->SetCellsAmplitudeFraction(ratios);
304 clus->SetDispersion(recpoint->GetDispersion());
305 clus->SetChi2(-1); //not yet implemented
306 clus->SetTOF(recpoint->GetTime()) ; //time-of-flight
307 clus->SetNExMax(recpoint->GetNExMax()); //number of local maxima
308 Float_t elipAxis[2];
309 recpoint->GetElipsAxis(elipAxis);
310 clus->SetM02(elipAxis[0]*elipAxis[0]) ;
311 clus->SetM20(elipAxis[1]*elipAxis[1]) ;
312 clus->SetDistToBadChannel(recpoint->GetDistanceToBadTower());
313 if (esdevent && fRecoUtils) {
314 Int_t trackIndex = fRecoUtils->GetMatchedTrackIndex(i);
315 if(trackIndex >= 0) {
316 clus->AddTrackMatched(esdevent->GetTrack(trackIndex));
317 if(DebugLevel() > 1)
318 AliInfo(Form("Matched Track index %d to new cluster %d\n",trackIndex,i));
319 }
320 }
321 }
322}
323
324//________________________________________________________________________________________
325void AliAnalysisTaskEMCALClusterizeFast::Init()
326{
327 //Select clusterization/unfolding algorithm and set all the needed parameters
328
329 AliVEvent * event = InputEvent();
330 if (!event) {
331 AliWarning("Event not available!!!");
332 return;
333 }
334
335 if (event->GetRunNumber()==fRun)
336 return;
337 fRun = event->GetRunNumber();
338
339 if (fJustUnfold){
340 // init the unfolding afterburner
341 delete fUnfolder;
342 fUnfolder = new AliEMCALAfterBurnerUF(fRecParam->GetW0(),fRecParam->GetLocMaxCut());
343 return;
344 }
345
346 AliCDBManager *cdb = AliCDBManager::Instance();
347 if (!cdb->IsDefaultStorageSet() && !fOCDBpath.IsNull())
348 cdb->SetDefaultStorage(fOCDBpath);
349 if (fRun!=cdb->GetRun())
350 cdb->SetRun(fRun);
351
352 AliEMCALGeometry *geometry = AliEMCALGeometry::GetInstance(fGeomName);
353 if (!geometry) {
354 AliFatal("Geometry not available!!!");
355 return;
356 }
357
358 if (!fGeomMatrixSet) {
359 if (fLoadGeomMatrices) {
360 for(Int_t mod=0; mod < (geometry->GetEMCGeometry())->GetNumberOfSuperModules(); ++mod) {
361 if(fGeomMatrix[mod]){
362 if(DebugLevel() > 2)
363 fGeomMatrix[mod]->Print();
364 geometry->SetMisalMatrix(fGeomMatrix[mod],mod);
365 }
366 }
367 } else {
368 for(Int_t mod=0; mod < geometry->GetEMCGeometry()->GetNumberOfSuperModules(); ++mod) {
369 if(event->GetEMCALMatrix(mod)) {
370 if(DebugLevel() > 2)
371 event->GetEMCALMatrix(mod)->Print();
372 geometry->SetMisalMatrix(event->GetEMCALMatrix(mod),mod);
373 }
374 }
375 }
376 fGeomMatrixSet=kTRUE;
377 }
378
379 // setup digit array if needed
380 if (!fDigitsArr) {
381 fDigitsArr = new TClonesArray("AliEMCALDigit", 1000);
382 fDigitsArr->SetOwner(1);
383 }
384
385 // then setup clusterizer
386 delete fClusterizer;
387 if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv1)
388 fClusterizer = new AliEMCALClusterizerv1(geometry);
389 else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerNxN)
390 fClusterizer = new AliEMCALClusterizerNxN(geometry);
391 else if(fRecParam->GetClusterizerFlag() > AliEMCALRecParam::kClusterizerNxN) {
392 AliEMCALClusterizerNxN *clusterizer = new AliEMCALClusterizerNxN(geometry);
393 clusterizer->SetNRowDiff(2);
394 clusterizer->SetNColDiff(2);
395 fClusterizer = clusterizer;
396 } else{
397 AliFatal(Form("Clusterizer < %d > not available", fRecParam->GetClusterizerFlag()));
398 }
399 fClusterizer->InitParameters(fRecParam);
10d33986 400 if (!fCalibData&&fLoadCalib) {
2f7259cf 401 AliCDBEntry *entry = static_cast<AliCDBEntry*>(AliCDBManager::Instance()->Get("EMCAL/Calib/Data"));
402 if (entry)
403 fCalibData = static_cast<AliEMCALCalibData*>(entry->GetObject());
404 if (!fCalibData)
405 AliFatal("Calibration parameters not found in CDB!");
406 }
10d33986 407 if (!fPedestalData&&fLoadPed) {
2f7259cf 408 AliCDBEntry *entry = static_cast<AliCDBEntry*>(AliCDBManager::Instance()->Get("EMCAL/Calib/Pedestals"));
409 if (entry)
410 fPedestalData = static_cast<AliCaloCalibPedestal*>(entry->GetObject());
411 }
10d33986 412 if (fCalibData) {
413 fClusterizer->SetInputCalibrated(kFALSE);
414 fClusterizer->SetCalibrationParameters(fCalibData);
415 fClusterizer->SetCaloCalibPedestal(fPedestalData);
416 } else {
417 fClusterizer->SetInputCalibrated(kTRUE);
418 }
419 fClusterizer->SetJustClusters(kTRUE);
2f7259cf 420 fClusterizer->SetDigitsArr(fDigitsArr);
421 fClusterizer->SetOutput(0);
422 fClusterArr = const_cast<TObjArray *>(fClusterizer->GetRecPoints());
423}