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