revert previous commit on this class
[u/mrichter/AliRoot.git] / PWGGA / EMCALTasks / 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>
7a7d5a64 24
25// --- AliRoot ---
2f7259cf 26#include "AliAODCaloCluster.h"
27#include "AliAODEvent.h"
28#include "AliAnalysisManager.h"
29#include "AliCDBEntry.h"
30#include "AliCDBManager.h"
31#include "AliCaloCalibPedestal.h"
32#include "AliEMCALAfterBurnerUF.h"
33#include "AliEMCALCalibData.h"
34#include "AliEMCALClusterizerNxN.h"
35#include "AliEMCALClusterizerv1.h"
3786256f 36#include "AliEMCALClusterizerv2.h"
7a7d5a64 37#include "AliEMCALClusterizerFixedWindow.h"
2f7259cf 38#include "AliEMCALDigit.h"
39#include "AliEMCALGeometry.h"
40#include "AliEMCALRecParam.h"
41#include "AliEMCALRecPoint.h"
42#include "AliEMCALRecoUtils.h"
43#include "AliESDEvent.h"
6b7ed7bd 44#include "AliInputEventHandler.h"
2f7259cf 45#include "AliLog.h"
46
47#include "AliAnalysisTaskEMCALClusterizeFast.h"
48
49ClassImp(AliAnalysisTaskEMCALClusterizeFast)
50
51//________________________________________________________________________
2431b20f 52AliAnalysisTaskEMCALClusterizeFast::AliAnalysisTaskEMCALClusterizeFast()
53 : AliAnalysisTaskSE(),
349bf69e 54 fRun(-1),
2431b20f 55 fDigitsArr(0),
56 fClusterArr(0),
57 fRecParam(0),
58 fClusterizer(0),
59 fUnfolder(0),
60 fJustUnfold(kFALSE),
61 fGeomName(),
62 fGeomMatrixSet(kFALSE),
63 fLoadGeomMatrices(kFALSE),
64 fOCDBpath(),
65 fCalibData(0),
66 fPedestalData(0),
67 fOutputAODBranch(0),
68 fOutputAODBrName(),
69 fRecoUtils(0),
70 fLoadCalib(0),
4cf9204b 71 fLoadPed(0),
e50a444f 72 fAttachClusters(1),
95f69c66 73 fRecalibOnly(0),
7a7d5a64 74 fSubBackground(0),
75 fCreatePattern(0),
e50a444f 76 fOverwrite(0),
7a7d5a64 77 fNewClusterArrayName("newCaloClusters"),
78 fNPhi(4),
79 fNEta(4),
3c56da8f 80 fShiftPhi(2),
81 fShiftEta(2),
82 fTRUShift(0),
0408f54f 83 fClusterizeFastORs(0),
be268574 84 fTrackName(),
85 fCutL0Times(kTRUE)
2431b20f 86{
87 // Constructor
c6d60de7 88
89 for(Int_t i = 0; i < 12; ++i)
90 fGeomMatrix[i] = 0;
2431b20f 91}
92
93//________________________________________________________________________
2f7259cf 94AliAnalysisTaskEMCALClusterizeFast::AliAnalysisTaskEMCALClusterizeFast(const char *name)
95 : AliAnalysisTaskSE(name),
349bf69e 96 fRun(-1),
2f7259cf 97 fDigitsArr(0),
98 fClusterArr(0),
99 fRecParam(new AliEMCALRecParam),
100 fClusterizer(0),
101 fUnfolder(0),
102 fJustUnfold(kFALSE),
103 fGeomName("EMCAL_FIRSTYEARV1"),
104 fGeomMatrixSet(kFALSE),
105 fLoadGeomMatrices(kFALSE),
106 fOCDBpath(),
107 fCalibData(0),
108 fPedestalData(0),
109 fOutputAODBranch(0),
110 fOutputAODBrName(),
10d33986 111 fRecoUtils(0),
112 fLoadCalib(0),
4cf9204b 113 fLoadPed(0),
e50a444f 114 fAttachClusters(1),
95f69c66 115 fRecalibOnly(0),
7a7d5a64 116 fSubBackground(0),
117 fCreatePattern(0),
e50a444f 118 fOverwrite(0),
7a7d5a64 119 fNewClusterArrayName("newCaloClusters"),
120 fNPhi(4),
121 fNEta(4),
3c56da8f 122 fShiftPhi(2),
123 fShiftEta(2),
124 fTRUShift(0),
0408f54f 125 fClusterizeFastORs(0),
be268574 126 fTrackName(),
127 fCutL0Times(kTRUE)
2f7259cf 128{
129 // Constructor
130
0a7eb283 131 fBranchNames = "ESD:AliESDHeader.,AliESDRun.,EMCALCells.,EMCALTrigger. AOD:header,emcalCells";
95f69c66 132 for(Int_t i = 0; i < 12; ++i)
2f7259cf 133 fGeomMatrix[i] = 0;
134}
135
136//________________________________________________________________________
137AliAnalysisTaskEMCALClusterizeFast::~AliAnalysisTaskEMCALClusterizeFast()
138{
139 // Destructor.
140
2f7259cf 141 delete fClusterizer;
142 delete fUnfolder;
143 delete fRecoUtils;
144}
145
146//-------------------------------------------------------------------
147void AliAnalysisTaskEMCALClusterizeFast::UserCreateOutputObjects()
148{
149 // Create output objects.
150
151 if (!fOutputAODBrName.IsNull()) {
152 fOutputAODBranch = new TClonesArray("AliAODCaloCluster", 0);
153 fOutputAODBranch->SetName(fOutputAODBrName);
154 AddAODBranch("TClonesArray", &fOutputAODBranch);
155 AliInfo(Form("Created Branch: %s",fOutputAODBrName.Data()));
156 }
157}
158
159//________________________________________________________________________
160void AliAnalysisTaskEMCALClusterizeFast::UserExec(Option_t *)
161{
162 // Main loop, called for each event
163
164 // remove the contents of output list set in the previous event
165 if (fOutputAODBranch)
166 fOutputAODBranch->Clear("C");
167
168 AliESDEvent *esdevent = dynamic_cast<AliESDEvent*>(InputEvent());
169 AliAODEvent *aodevent = dynamic_cast<AliAODEvent*>(InputEvent());
170
171 if (!esdevent&&!aodevent) {
172 Error("UserExec","Event not available");
173 return;
174 }
175
176 LoadBranches();
6b7ed7bd 177
178 UInt_t offtrigger = 0;
179 if (esdevent) {
9809ed8c 180 UInt_t mask1 = esdevent->GetESDRun()->GetDetectorsInDAQ();
181 UInt_t mask2 = esdevent->GetESDRun()->GetDetectorsInReco();
182 Bool_t desc1 = (mask1 >> 18) & 0x1;
183 Bool_t desc2 = (mask2 >> 18) & 0x1;
184 if (desc1==0 || desc2==0) { //AliDAQ::OfflineModuleName(180=="EMCAL"
185 AliError(Form("EMCAL not in DAQ/RECO: %u (%u)/%u (%u)",
186 mask1, esdevent->GetESDRun()->GetDetectorsInReco(),
187 mask2, esdevent->GetESDRun()->GetDetectorsInDAQ()));
6b7ed7bd 188 return;
189 }
190 AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
191 offtrigger = ((AliInputEventHandler*)(am->GetInputEventHandler()))->IsEventSelected();
192 } else if (aodevent) {
193 offtrigger = aodevent->GetHeader()->GetOfflineTrigger();
194 }
179c3e32 195
2ecad9b5 196 if (esdevent) {
179c3e32 197 AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
198 Bool_t mcmode = 0;
199 if (am->GetMCtruthEventHandler())
200 mcmode = 1;
201 if (!mcmode) {
179c3e32 202 if (offtrigger & AliVEvent::kFastOnly) {
203 AliWarning(Form("EMCAL not in fast only partition"));
204 return;
205 }
206 }
6b7ed7bd 207 }
2f7259cf 208
209 Init();
210
211 if (fJustUnfold) {
212 AliWarning("Unfolding not implemented");
84ccad86 213 return;
214 }
215
95f69c66 216 FillDigitsArray();
84ccad86 217
95f69c66 218 if (fRecalibOnly) {
219 UpdateCells();
220 return; // not requested to run clusterizer
84ccad86 221 }
7099ec08 222
95f69c66 223 Clusterize();
224 UpdateCells();
225 UpdateClusters();
2431b20f 226
95f69c66 227 if (fOutputAODBranch)
228 RecPoints2Clusters(fOutputAODBranch);
2431b20f 229}
230
231//________________________________________________________________________
95f69c66 232void AliAnalysisTaskEMCALClusterizeFast::Clusterize()
2431b20f 233{
95f69c66 234 // Clusterize
7099ec08 235
95f69c66 236 if (fSubBackground) {
237 fClusterizer->SetInputCalibrated(kTRUE);
238 fClusterizer->SetCalibrationParameters(0);
2f7259cf 239 }
d3b5a5ec 240
95f69c66 241 fClusterizer->Digits2Clusters("");
242 if (fSubBackground) {
4889ed8b 243 if (fCalibData) {
244 fClusterizer->SetInputCalibrated(kFALSE);
245 fClusterizer->SetCalibrationParameters(fCalibData);
246 }
4cf9204b 247 }
4cf9204b 248}
249
250//________________________________________________________________________
95f69c66 251void AliAnalysisTaskEMCALClusterizeFast::FillDigitsArray()
2f7259cf 252{
6265e793 253 // Fill digits array
55a8e3e4 254
6265e793 255 fDigitsArr->Clear("C");
256
5133b922 257 if (fCreatePattern) { // Fill digits from a pattern
7a7d5a64 258 AliEMCALGeometry *fGeom = AliEMCALGeometry::GetInstance(fGeomName);
7a7d5a64 259 Int_t maxd = fGeom->GetNCells() / 4;
6265e793 260 for (Int_t idigit = 0; idigit < maxd; idigit++){
7a7d5a64 261 if (idigit % 24 == 12) idigit += 12;
262 AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->New(idigit));
263 digit->SetId(idigit * 4);
264 digit->SetTime(600);
265 digit->SetTimeR(600);
266 digit->SetIndexInList(idigit);
267 digit->SetType(AliEMCALDigit::kHG);
268 digit->SetAmplitude(0.1);
b0e354fe 269 }
06ca320f 270
5133b922 271 } else if (fClusterizeFastORs) { // Fill digits from FastORs
6265e793 272
273 AliEMCALGeometry *fGeom = AliEMCALGeometry::GetInstance(fGeomName);
06ca320f 274
6265e793 275 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
276 if (!esd){
277 AliError("Cannot get the ESD event");
278 return;
279 }
280
281 AliESDCaloTrigger *triggers = esd->GetCaloTrigger("EMCAL");
282
283 if (!triggers || !(triggers->GetEntries() > 0))
284 return;
285
286 Int_t idigit = 0;
287 triggers->Reset();
288
5133b922 289 while ((triggers->Next())) {
6265e793 290 Float_t triggerAmplitude = 0;
291 triggers->GetAmplitude(triggerAmplitude);
292 if (triggerAmplitude <= 0)
293 continue;
294
295 Int_t triggerTime = 0;
296 Int_t ntimes = 0;
297 triggers->GetNL0Times(ntimes);
be268574 298 if (!(ntimes > 0) && fCutL0Times)
299 continue;
300
301 Int_t trgtimes[25];
302 triggers->GetL0Times(trgtimes);
303 triggerTime = trgtimes[0];
304
6265e793 305 Int_t triggerCol = 0, triggerRow = 0;
306 triggers->GetPosition(triggerCol, triggerRow);
307
308 Int_t find = -1;
309 fGeom->GetAbsFastORIndexFromPositionInEMCAL(triggerCol, triggerRow, find);
310
311 if (find<0)
312 continue;
313
314 Int_t cidx[4] = {-1};
315 Bool_t ret = fGeom->GetCellIndexFromFastORIndex(find, cidx);
316
317 if (!ret)
318 continue;
319
5133b922 320 for (Int_t idxpos = 0; idxpos < 4; idxpos++) {
6265e793 321 Int_t triggerNumber = cidx[idxpos];
322 AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->New(idigit));
323 digit->SetId(triggerNumber);
324 digit->SetTime(triggerTime);
325 digit->SetTimeR(triggerTime);
326 digit->SetIndexInList(idigit);
327 digit->SetType(AliEMCALDigit::kHG);
328 digit->SetAmplitude(triggerAmplitude);
329 idigit++;
330 }
331 }
06ca320f 332
333 } else { // Fill digits from cells.
7a7d5a64 334
7a7d5a64 335 AliVCaloCells *cells = InputEvent()->GetEMCALCells();
55a8e3e4 336 Double_t avgE = 0; // for background subtraction
337 const Int_t ncells = cells->GetNumberOfCells();
7a7d5a64 338 for (Int_t icell = 0, idigit = 0; icell < ncells; ++icell) {
339 Double_t cellAmplitude=0, cellTime=0;
340 Short_t cellNumber=0;
341 if (cells->GetCell(icell, cellNumber, cellAmplitude, cellTime) != kTRUE)
342 break;
343 AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->New(idigit));
344 digit->SetId(cellNumber);
345 digit->SetTime(cellTime);
346 digit->SetTimeR(cellTime);
347 digit->SetIndexInList(idigit);
348 digit->SetType(AliEMCALDigit::kHG);
349 if (fRecalibOnly||fSubBackground) {
350 Float_t energy = cellAmplitude;
55a8e3e4 351 Float_t time = cellTime;
7a7d5a64 352 fClusterizer->Calibrate(energy,time,cellNumber);
95f69c66 353 digit->SetAmplitude(energy);
7a7d5a64 354 avgE += energy;
355 } else {
356 digit->SetAmplitude(cellAmplitude);
95f69c66 357 }
7a7d5a64 358 idigit++;
b0e354fe 359 }
7a7d5a64 360
55a8e3e4 361 fDigitsArr->Sort();
362
7a7d5a64 363 if (fSubBackground) {
364 avgE /= AliEMCALGeometry::GetInstance(fGeomName)->GetNumberOfSuperModules()*48*24;
365 Int_t ndigis = fDigitsArr->GetEntries();
366 for (Int_t i = 0; i < ndigis; ++i) {
367 AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(i));
368 Double_t energy = digit->GetAmplitude() - avgE;
369 if (energy<=0.001) {
370 digit->SetAmplitude(0);
371 } else {
372 digit->SetAmplitude(energy);
373 }
374 }
375 }
06ca320f 376 }
2f7259cf 377}
378
379//________________________________________________________________________________________
95f69c66 380void AliAnalysisTaskEMCALClusterizeFast::RecPoints2Clusters(TClonesArray *clus)
4cf9204b 381{
382 // Cluster energy, global position, cells and their amplitude fractions are restored.
383
95f69c66 384 Bool_t esdobjects = 0;
385 if (strcmp(clus->GetClass()->GetName(),"AliESDCaloCluster")==0)
386 esdobjects = 1;
387
b98487bc 388 AliVCaloCells *cells = InputEvent()->GetEMCALCells();
389 AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance(fGeomName);
0408f54f 390 AliESDEvent *esdevent = dynamic_cast<AliESDEvent*>(InputEvent());
31f5b395 391 if (!esdevent) {
392 AliError("Cannot get the ESD event");
393 return;
394 }
0408f54f 395
396 // tracks array for track/cluster matching
397 TClonesArray *tarr = 0;
398 if (!fTrackName.IsNull()) {
399 tarr = dynamic_cast<TClonesArray*>(esdevent->FindListObject(fTrackName));
400 if (!tarr) {
401 AliError(Form("Cannot get tracks named %s", fTrackName.Data()));
402 }
403 }
e93ec1f7 404
eb33b7f8 405 const Int_t Ncls = fClusterArr->GetEntries();
0408f54f 406 AliDebug(1, Form("total no of clusters %d", Ncls));
4cf9204b 407 for(Int_t i=0, nout=clus->GetEntries(); i < Ncls; ++i) {
408 AliEMCALRecPoint *recpoint = static_cast<AliEMCALRecPoint*>(fClusterArr->At(i));
409 Int_t ncells_true = 0;
410 const Int_t ncells = recpoint->GetMultiplicity();
411 UShort_t absIds[ncells];
412 Double32_t ratios[ncells];
413 Int_t *dlist = recpoint->GetDigitsList();
414 Float_t *elist = recpoint->GetEnergiesList();
415 for (Int_t c = 0; c < ncells; ++c) {
416 AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(dlist[c]));
417 absIds[ncells_true] = digit->GetId();
418 ratios[ncells_true] = elist[c]/digit->GetAmplitude();
84ccad86 419 ++ncells_true;
4cf9204b 420 }
421
422 if (ncells_true < 1) {
423 AliWarning("Skipping cluster with no cells");
424 continue;
425 }
426
427 // calculate new cluster position
428 TVector3 gpos;
4cf9204b 429 recpoint->GetGlobalPosition(gpos);
84ccad86 430 Float_t g[3];
4cf9204b 431 gpos.GetXYZ(g);
432
e93ec1f7 433 AliDebug(1, Form("energy %f", recpoint->GetEnergy()));
434
95f69c66 435 AliVCluster *c = static_cast<AliVCluster*>(clus->New(nout++));
4cf9204b 436 c->SetType(AliVCluster::kEMCALClusterv1);
437 c->SetE(recpoint->GetEnergy());
438 c->SetPosition(g);
439 c->SetNCells(ncells_true);
eb33b7f8 440 if (esdobjects) {
441 AliESDCaloCluster *cesd = static_cast<AliESDCaloCluster*>(c);
442 cesd->SetCellsAbsId(absIds);
443 cesd->SetCellsAmplitudeFraction(ratios);
444 cesd->SetID(recpoint->GetUniqueID());
445 } else {
446 AliAODCaloCluster *caod = static_cast<AliAODCaloCluster*>(c);
447 caod->SetCellsAbsId(absIds);
448 caod->SetCellsAmplitudeFraction(ratios);
449 }
4cf9204b 450 c->SetDispersion(recpoint->GetDispersion());
0408f54f 451 c->SetEmcCpvDistance(-1);
452 c->SetChi2(-1);
4cf9204b 453 c->SetTOF(recpoint->GetTime()) ; //time-of-flight
454 c->SetNExMax(recpoint->GetNExMax()); //number of local maxima
455 Float_t elipAxis[2];
456 recpoint->GetElipsAxis(elipAxis);
e93ec1f7 457 c->SetM02(elipAxis[0]*elipAxis[0]);
458 c->SetM20(elipAxis[1]*elipAxis[1]);
459 if (fPedestalData) {
b98487bc 460 c->SetDistanceToBadChannel(recpoint->GetDistanceToBadTower());
eb33b7f8 461 } else {
e93ec1f7 462 if (fRecoUtils && fRecoUtils->IsBadChannelsRemovalSwitchedOn()) {
463 fRecoUtils->RecalculateClusterDistanceToBadChannel(geom, cells, c);
464 }
b98487bc 465 }
84ccad86 466
0408f54f 467 if (tarr) {
468 Double_t dEtaMin = 1e9;
469 Double_t dPhiMin = 1e9;
470 Int_t imin = -1;
471 Double_t ceta = gpos.Eta();
472 Double_t cphi = gpos.Phi();
473 const Int_t ntrks = tarr->GetEntries();
474 for(Int_t t = 0; t<ntrks; ++t) {
475 AliVTrack *track = static_cast<AliVTrack*>(tarr->At(t));
476 if (!track)
477 continue;
478 const AliExternalTrackParam *outp = track->GetOuterParam();
479 if (!outp)
480 continue;
481 Double_t trkPos[3] = {0.,0.,0.};
482 if (!outp->GetXYZ(trkPos))
483 continue;
484 TVector3 vec(trkPos[0],trkPos[1],trkPos[2]);
485 Double_t veta = vec.Eta();
486 Double_t vphi = vec.Phi();
487 if(vphi<0)
488 vphi += 2*TMath::Pi();
489 if (TMath::Abs(veta)>0.75 || (vphi<70*TMath::DegToRad()) || (vphi>190*TMath::DegToRad()))
490 continue;
491 Double_t dR = vec.DeltaR(gpos);
492 if(dR > 25)
493 continue;
494 Float_t tmpEta=0, tmpPhi=0;
495 if (0) {
496 AliExternalTrackParam trkParTemp(*outp); // retrieve the starting point every time before the extrapolation
497 Bool_t ret = fRecoUtils->ExtrapolateTrackToCluster(&trkParTemp, c, fRecoUtils->GetMass(), fRecoUtils->GetStep(), tmpEta, tmpPhi);
498 if (!ret)
499 continue;
500 } else {
501 tmpEta = ceta - veta;
502 tmpPhi = cphi - vphi;
503 }
504 if (TMath::Abs(tmpEta)<TMath::Abs(dEtaMin) && TMath::Abs(tmpPhi)<TMath::Abs(dPhiMin)) {
505 dEtaMin = tmpEta;
506 dPhiMin = tmpPhi;
507 imin = t;
508 }
95f69c66 509 }
eb33b7f8 510 c->SetEmcCpvDistance(imin);
511 c->SetTrackDistance(dPhiMin, dEtaMin);
2f7259cf 512 }
513 }
514}
515
95f69c66 516//________________________________________________________________________
517void AliAnalysisTaskEMCALClusterizeFast::UpdateCells()
518{
519 // Update cells in case re-calibration was done.
520
521 if (!fCalibData&&!fSubBackground)
522 return;
523
524 AliVCaloCells *cells = InputEvent()->GetEMCALCells();
55a8e3e4 525 const Int_t ncells = cells->GetNumberOfCells();
526 const Int_t ndigis = fDigitsArr->GetEntries();
95f69c66 527 if (ncells!=ndigis) {
528 cells->DeleteContainer();
529 cells->CreateContainer(ndigis);
530 }
531 for (Int_t idigit = 0; idigit < ndigis; ++idigit) {
532 AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(idigit));
533 Double_t cellAmplitude = digit->GetCalibAmp();
55a8e3e4 534 Short_t cellNumber = digit->GetId();
535 Double_t cellTime = digit->GetTime();
95f69c66 536 cells->SetCell(idigit, cellNumber, cellAmplitude, cellTime);
537 }
538}
539
540//________________________________________________________________________
541void AliAnalysisTaskEMCALClusterizeFast::UpdateClusters()
542{
543 // Update cells in case re-calibration was done.
7a7d5a64 544
95f69c66 545 if (!fAttachClusters)
546 return;
7a7d5a64 547
06ca320f 548 TClonesArray *clus = 0;
7a7d5a64 549
e50a444f 550 if (fOverwrite) {
7a7d5a64 551 clus = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("caloClusters"));
552 if (!clus)
553 clus = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("CaloClusters"));
06ca320f 554 if (!clus)
7a7d5a64 555 return;
556
55a8e3e4 557 const Int_t nents = clus->GetEntries();
7a7d5a64 558 for (Int_t i=0;i<nents;++i) {
559 AliVCluster *c = static_cast<AliVCluster*>(clus->At(i));
560 if (!c)
561 continue;
562 if (c->IsEMCAL())
563 delete clus->RemoveAt(i);
564 }
06ca320f 565 } else {
7a7d5a64 566 clus = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fNewClusterArrayName));
06ca320f 567 if (!clus) {
7a7d5a64 568 clus = new TClonesArray("AliESDCaloCluster");
569 clus->SetName(fNewClusterArrayName);
570 InputEvent()->AddObject(clus);
06ca320f 571 } else {
0a7eb283 572 clus->Delete();
7a7d5a64 573 }
574 }
575
95f69c66 576 RecPoints2Clusters(clus);
577}
578
2f7259cf 579//________________________________________________________________________________________
580void AliAnalysisTaskEMCALClusterizeFast::Init()
581{
582 //Select clusterization/unfolding algorithm and set all the needed parameters
01bc1ce8 583
2f7259cf 584 AliVEvent * event = InputEvent();
585 if (!event) {
586 AliWarning("Event not available!!!");
587 return;
588 }
589
590 if (event->GetRunNumber()==fRun)
591 return;
592 fRun = event->GetRunNumber();
593
594 if (fJustUnfold){
595 // init the unfolding afterburner
596 delete fUnfolder;
1759f743 597 fUnfolder = new AliEMCALAfterBurnerUF(fRecParam->GetW0(),fRecParam->GetLocMaxCut(),fRecParam->GetMinECut());
2f7259cf 598 return;
599 }
600
2f7259cf 601 AliEMCALGeometry *geometry = AliEMCALGeometry::GetInstance(fGeomName);
602 if (!geometry) {
603 AliFatal("Geometry not available!!!");
604 return;
605 }
606
607 if (!fGeomMatrixSet) {
608 if (fLoadGeomMatrices) {
e51bc3f5 609 for(Int_t mod=0; mod < geometry->GetNumberOfSuperModules(); ++mod) {
06ca320f 610 if (fGeomMatrix[mod]){
611 if (DebugLevel() > 2)
2f7259cf 612 fGeomMatrix[mod]->Print();
613 geometry->SetMisalMatrix(fGeomMatrix[mod],mod);
614 }
615 }
29b496d5 616 } else { // get matrix from file (work around bug in aliroot)
2f7259cf 617 for(Int_t mod=0; mod < geometry->GetEMCGeometry()->GetNumberOfSuperModules(); ++mod) {
29b496d5 618 const TGeoHMatrix *gm = 0;
441c4d3c 619 AliESDEvent *esdevent = dynamic_cast<AliESDEvent*>(event);
29b496d5 620 if (esdevent) {
621 gm = esdevent->GetEMCALMatrix(mod);
622 } else {
623 AliAODHeader *aodheader = dynamic_cast<AliAODHeader*>(event->GetHeader());
624 if (aodheader) {
625 gm = aodheader->GetEMCALMatrix(mod);
626 }
627 }
628 if (gm) {
06ca320f 629 if (DebugLevel() > 2)
29b496d5 630 gm->Print();
631 geometry->SetMisalMatrix(gm,mod);
2f7259cf 632 }
633 }
634 }
635 fGeomMatrixSet=kTRUE;
636 }
637
638 // setup digit array if needed
639 if (!fDigitsArr) {
640 fDigitsArr = new TClonesArray("AliEMCALDigit", 1000);
641 fDigitsArr->SetOwner(1);
642 }
643
644 // then setup clusterizer
645 delete fClusterizer;
646 if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv1)
647 fClusterizer = new AliEMCALClusterizerv1(geometry);
06ca320f 648 else if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerNxN) {
2f7259cf 649 AliEMCALClusterizerNxN *clusterizer = new AliEMCALClusterizerNxN(geometry);
3c56da8f 650 clusterizer->SetNRowDiff(fRecParam->GetNRowDiff());
651 clusterizer->SetNColDiff(fRecParam->GetNColDiff());
652 fClusterizer = clusterizer;
06ca320f 653 } else if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv2)
3c56da8f 654 fClusterizer = new AliEMCALClusterizerv2(geometry);
06ca320f 655 else if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerFW){
7a7d5a64 656 AliEMCALClusterizerFixedWindow *clusterizer = new AliEMCALClusterizerFixedWindow(geometry);
0fda97fe 657 clusterizer->SetNphi(fNPhi);
658 clusterizer->SetNeta(fNEta);
659 clusterizer->SetShiftPhi(fShiftPhi);
660 clusterizer->SetShiftEta(fShiftEta);
3c56da8f 661 clusterizer->SetTRUshift(fTRUShift);
7a7d5a64 662 fClusterizer = clusterizer;
663 }
664 else{
2f7259cf 665 AliFatal(Form("Clusterizer < %d > not available", fRecParam->GetClusterizerFlag()));
666 }
667 fClusterizer->InitParameters(fRecParam);
01bc1ce8 668
669 if ((!fCalibData&&fLoadCalib) || (!fPedestalData&&fLoadPed)) {
670 AliCDBManager *cdb = AliCDBManager::Instance();
671 if (!cdb->IsDefaultStorageSet() && !fOCDBpath.IsNull())
672 cdb->SetDefaultStorage(fOCDBpath);
673 if (fRun!=cdb->GetRun())
674 cdb->SetRun(fRun);
675 }
349bf69e 676 if (!fCalibData&&fLoadCalib&&fRun>0) {
2f7259cf 677 AliCDBEntry *entry = static_cast<AliCDBEntry*>(AliCDBManager::Instance()->Get("EMCAL/Calib/Data"));
678 if (entry)
679 fCalibData = static_cast<AliEMCALCalibData*>(entry->GetObject());
680 if (!fCalibData)
681 AliFatal("Calibration parameters not found in CDB!");
682 }
349bf69e 683 if (!fPedestalData&&fLoadPed&&fRun>0) {
2f7259cf 684 AliCDBEntry *entry = static_cast<AliCDBEntry*>(AliCDBManager::Instance()->Get("EMCAL/Calib/Pedestals"));
685 if (entry)
686 fPedestalData = static_cast<AliCaloCalibPedestal*>(entry->GetObject());
687 }
10d33986 688 if (fCalibData) {
689 fClusterizer->SetInputCalibrated(kFALSE);
690 fClusterizer->SetCalibrationParameters(fCalibData);
10d33986 691 } else {
692 fClusterizer->SetInputCalibrated(kTRUE);
693 }
e93ec1f7 694 fClusterizer->SetCaloCalibPedestal(fPedestalData);
10d33986 695 fClusterizer->SetJustClusters(kTRUE);
2f7259cf 696 fClusterizer->SetDigitsArr(fDigitsArr);
697 fClusterizer->SetOutput(0);
698 fClusterArr = const_cast<TObjArray *>(fClusterizer->GetRecPoints());
699}