]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG/EMCAL/AliAnalysisTaskEMCALClusterizeFast.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
New tracking efficiency task that uses the EMCal jet framework
[u/mrichter/AliRoot.git] / PWG / EMCAL / 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
cd231d42 16//_________________________________________________________________________
17// This analysis provides a new list of clusters to be used in other analysis
18//
19// Author: Constantin Loizides, Salvatore Aiola
20// Adapted from analysis class from Deepa Thomas
21//
22// $Id$
23//_________________________________________________________________________
2f7259cf 24
25// --- Root ---
26#include <TClonesArray.h>
27#include <TGeoManager.h>
28#include <TObjArray.h>
29#include <TString.h>
30#include <TTree.h>
507f74bc 31#include <TArrayI.h>
7a7d5a64 32
33// --- AliRoot ---
2f7259cf 34#include "AliAODCaloCluster.h"
35#include "AliAODEvent.h"
36#include "AliAnalysisManager.h"
37#include "AliCDBEntry.h"
38#include "AliCDBManager.h"
39#include "AliCaloCalibPedestal.h"
40#include "AliEMCALAfterBurnerUF.h"
41#include "AliEMCALCalibData.h"
42#include "AliEMCALClusterizerNxN.h"
43#include "AliEMCALClusterizerv1.h"
3786256f 44#include "AliEMCALClusterizerv2.h"
7a7d5a64 45#include "AliEMCALClusterizerFixedWindow.h"
2f7259cf 46#include "AliEMCALDigit.h"
47#include "AliEMCALGeometry.h"
48#include "AliEMCALRecParam.h"
49#include "AliEMCALRecPoint.h"
50#include "AliEMCALRecoUtils.h"
51#include "AliESDEvent.h"
6b7ed7bd 52#include "AliInputEventHandler.h"
2f7259cf 53#include "AliLog.h"
54
55#include "AliAnalysisTaskEMCALClusterizeFast.h"
56
57ClassImp(AliAnalysisTaskEMCALClusterizeFast)
58
2431b20f 59//________________________________________________________________________
507f74bc 60AliAnalysisTaskEMCALClusterizeFast::AliAnalysisTaskEMCALClusterizeFast() :
61 AliAnalysisTaskSE(),
62 fRun(-1),
63 fDigitsArr(0),
64 fClusterArr(0),
65 fRecParam(new AliEMCALRecParam),
66 fClusterizer(0),
67 fUnfolder(0),
68 fJustUnfold(kFALSE),
69 fGeomName(),
70 fGeomMatrixSet(kFALSE),
71 fLoadGeomMatrices(kFALSE),
72 fOCDBpath(),
73 fCalibData(0),
74 fPedestalData(0),
75 fOutputAODBranch(0),
76 fOutputAODBrName(),
77 fRecoUtils(0),
78 fLoadCalib(kFALSE),
79 fLoadPed(kFALSE),
80 fAttachClusters(kTRUE),
81 fSubBackground(kFALSE),
82 fNPhi(4),
83 fNEta(4),
84 fShiftPhi(2),
85 fShiftEta(2),
86 fTRUShift(0),
87 fInputCellType(kFEEData),
88 fTrackName(),
89 fCaloCellsName(),
90 fCaloClustersName("newCaloClusters"),
788b1866 91 fDoUpdateCells(kFALSE),
e4d5da1e 92 fDoClusterize(kTRUE),
788b1866 93 fClusterBadChannelCheck(kFALSE),
507f74bc 94 fRejectExoticClusters(kFALSE),
eff32c61 95 fRejectExoticCells(kFALSE),
507f74bc 96 fFiducial(kFALSE),
97 fDoNonLinearity(kFALSE),
788b1866 98 fRecalDistToBadChannels(kFALSE),
507f74bc 99 fCaloCells(0),
100 fCaloClusters(0),
101 fEsd(0),
102 fAod(0),
103 fGeom(0)
2431b20f 104{
105 // Constructor
c6d60de7 106
107 for(Int_t i = 0; i < 12; ++i)
108 fGeomMatrix[i] = 0;
2431b20f 109}
110
2f7259cf 111//________________________________________________________________________
507f74bc 112AliAnalysisTaskEMCALClusterizeFast::AliAnalysisTaskEMCALClusterizeFast(const char *name) :
113 AliAnalysisTaskSE(name),
114 fRun(-1),
115 fDigitsArr(0),
116 fClusterArr(0),
117 fRecParam(new AliEMCALRecParam),
118 fClusterizer(0),
119 fUnfolder(0),
120 fJustUnfold(kFALSE),
121 fGeomName(),
122 fGeomMatrixSet(kFALSE),
123 fLoadGeomMatrices(kFALSE),
124 fOCDBpath(),
125 fCalibData(0),
126 fPedestalData(0),
127 fOutputAODBranch(0),
128 fOutputAODBrName(),
129 fRecoUtils(0),
130 fLoadCalib(kFALSE),
131 fLoadPed(kFALSE),
132 fAttachClusters(kTRUE),
133 fSubBackground(kFALSE),
134 fNPhi(4),
135 fNEta(4),
136 fShiftPhi(2),
137 fShiftEta(2),
138 fTRUShift(0),
139 fInputCellType(kFEEData),
140 fTrackName(),
141 fCaloCellsName(),
142 fCaloClustersName("newCaloClusters"),
788b1866 143 fDoUpdateCells(kFALSE),
e4d5da1e 144 fDoClusterize(kTRUE),
788b1866 145 fClusterBadChannelCheck(kFALSE),
507f74bc 146 fRejectExoticClusters(kFALSE),
eff32c61 147 fRejectExoticCells(kFALSE),
507f74bc 148 fFiducial(kFALSE),
149 fDoNonLinearity(kFALSE),
788b1866 150 fRecalDistToBadChannels(kFALSE),
507f74bc 151 fCaloCells(0),
152 fCaloClusters(0),
153 fEsd(0),
154 fAod(0),
155 fGeom(0)
2f7259cf 156{
157 // Constructor
158
0a7eb283 159 fBranchNames = "ESD:AliESDHeader.,AliESDRun.,EMCALCells.,EMCALTrigger. AOD:header,emcalCells";
95f69c66 160 for(Int_t i = 0; i < 12; ++i)
2f7259cf 161 fGeomMatrix[i] = 0;
162}
163
164//________________________________________________________________________
165AliAnalysisTaskEMCALClusterizeFast::~AliAnalysisTaskEMCALClusterizeFast()
166{
167 // Destructor.
168
2f7259cf 169 delete fClusterizer;
170 delete fUnfolder;
171 delete fRecoUtils;
507f74bc 172 delete fRecParam;
2f7259cf 173}
174
507f74bc 175//________________________________________________________________________
2f7259cf 176void AliAnalysisTaskEMCALClusterizeFast::UserCreateOutputObjects()
177{
178 // Create output objects.
179
180 if (!fOutputAODBrName.IsNull()) {
181 fOutputAODBranch = new TClonesArray("AliAODCaloCluster", 0);
182 fOutputAODBranch->SetName(fOutputAODBrName);
183 AddAODBranch("TClonesArray", &fOutputAODBranch);
184 AliInfo(Form("Created Branch: %s",fOutputAODBrName.Data()));
185 }
186}
187
188//________________________________________________________________________
189void AliAnalysisTaskEMCALClusterizeFast::UserExec(Option_t *)
190{
191 // Main loop, called for each event
192
30820b25 193 // if no operation is requested, return
194 if (!fDoClusterize && !fDoUpdateCells)
195 return;
196
2f7259cf 197 // remove the contents of output list set in the previous event
198 if (fOutputAODBranch)
199 fOutputAODBranch->Clear("C");
200
507f74bc 201 fEsd = dynamic_cast<AliESDEvent*>(InputEvent());
202 fAod = dynamic_cast<AliAODEvent*>(InputEvent());
2f7259cf 203
507f74bc 204 if (!fEsd&&!fAod) {
2f7259cf 205 Error("UserExec","Event not available");
206 return;
207 }
208
209 LoadBranches();
6b7ed7bd 210
211 UInt_t offtrigger = 0;
507f74bc 212 if (fEsd) {
213 UInt_t mask1 = fEsd->GetESDRun()->GetDetectorsInDAQ();
214 UInt_t mask2 = fEsd->GetESDRun()->GetDetectorsInReco();
9809ed8c 215 Bool_t desc1 = (mask1 >> 18) & 0x1;
216 Bool_t desc2 = (mask2 >> 18) & 0x1;
217 if (desc1==0 || desc2==0) { //AliDAQ::OfflineModuleName(180=="EMCAL"
218 AliError(Form("EMCAL not in DAQ/RECO: %u (%u)/%u (%u)",
507f74bc 219 mask1, fEsd->GetESDRun()->GetDetectorsInReco(),
220 mask2, fEsd->GetESDRun()->GetDetectorsInDAQ()));
6b7ed7bd 221 return;
222 }
223 AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
224 offtrigger = ((AliInputEventHandler*)(am->GetInputEventHandler()))->IsEventSelected();
507f74bc 225 } else {
226 offtrigger = fAod->GetHeader()->GetOfflineTrigger();
6b7ed7bd 227 }
179c3e32 228
507f74bc 229 if (!MCEvent()) {
230 if (offtrigger & AliVEvent::kFastOnly) {
231 AliError(Form("EMCAL not in fast only partition"));
232 return;
179c3e32 233 }
6b7ed7bd 234 }
2f7259cf 235
236 Init();
237
238 if (fJustUnfold) {
239 AliWarning("Unfolding not implemented");
84ccad86 240 return;
241 }
242
95f69c66 243 FillDigitsArray();
84ccad86 244
507f74bc 245 if (fDoClusterize)
246 Clusterize();
247
248 if (fDoUpdateCells)
95f69c66 249 UpdateCells();
7099ec08 250
7eff06f2 251 if (!fDoClusterize || (!fAttachClusters && !fOutputAODBranch) || !fCaloClusters)
507f74bc 252 return;
253
95f69c66 254 UpdateClusters();
507f74bc 255 CalibrateClusters();
2431b20f 256
7eff06f2 257 if (fOutputAODBranch && fCaloClusters != fOutputAODBranch)
507f74bc 258 CopyClusters(fCaloClusters, fOutputAODBranch);
259}
260
261//________________________________________________________________________
262void AliAnalysisTaskEMCALClusterizeFast::CopyClusters(TClonesArray *orig, TClonesArray *dest)
263{
264 const Int_t Ncls = orig->GetEntries();
265
266 for(Int_t i=0; i < Ncls; ++i) {
267 AliVCluster *oc = static_cast<AliVCluster*>(orig->At(i));
268
269 if (!oc)
270 continue;
271
272 if (!oc->IsEMCAL())
273 continue;
274
275 AliVCluster *dc = static_cast<AliVCluster*>(dest->New(i));
276 dc->SetType(AliVCluster::kEMCALClusterv1);
277 dc->SetE(oc->E());
278 Float_t pos[3] = {0};
279 oc->GetPosition(pos);
280 dc->SetPosition(pos);
281 dc->SetNCells(oc->GetNCells());
282 dc->SetCellsAbsId(oc->GetCellsAbsId());
283 dc->SetCellsAmplitudeFraction(oc->GetCellsAmplitudeFraction());
284 dc->SetID(oc->GetID());
285 dc->SetDispersion(oc->GetDispersion());
286 dc->SetEmcCpvDistance(-1);
287 dc->SetChi2(-1);
288 dc->SetTOF(oc->GetTOF()); //time-of-flight
289 dc->SetNExMax(oc->GetNExMax()); //number of local maxima
290 dc->SetM02(oc->GetM02());
291 dc->SetM20(oc->GetM20());
292 dc->SetDistanceToBadChannel(oc->GetDistanceToBadChannel());
7eff06f2 293 dc->SetMCEnergyFraction(oc->GetMCEnergyFraction());
507f74bc 294
295 //MC
296 UInt_t nlabels = oc->GetNLabels();
297 Int_t *labels = oc->GetLabels();
298
299 if (nlabels == 0 || !labels)
300 continue;
301
302 AliESDCaloCluster *esdClus = dynamic_cast<AliESDCaloCluster*>(dc);
303 if (esdClus) {
304 TArrayI parents(nlabels, labels);
305 esdClus->AddLabels(parents);
306 }
307 else {
308 AliAODCaloCluster *aodClus = dynamic_cast<AliAODCaloCluster*>(dc);
309 if (aodClus)
310 aodClus->SetLabel(labels, nlabels);
311 }
312 }
2431b20f 313}
314
315//________________________________________________________________________
95f69c66 316void AliAnalysisTaskEMCALClusterizeFast::Clusterize()
2431b20f 317{
95f69c66 318 // Clusterize
7099ec08 319
95f69c66 320 if (fSubBackground) {
321 fClusterizer->SetInputCalibrated(kTRUE);
322 fClusterizer->SetCalibrationParameters(0);
2f7259cf 323 }
d3b5a5ec 324
95f69c66 325 fClusterizer->Digits2Clusters("");
fde82e42 326
95f69c66 327 if (fSubBackground) {
4889ed8b 328 if (fCalibData) {
329 fClusterizer->SetInputCalibrated(kFALSE);
330 fClusterizer->SetCalibrationParameters(fCalibData);
331 }
4cf9204b 332 }
4cf9204b 333}
334
2f7259cf 335//________________________________________________________________________
95f69c66 336void AliAnalysisTaskEMCALClusterizeFast::FillDigitsArray()
2f7259cf 337{
6265e793 338 // Fill digits array
55a8e3e4 339
3119d479 340 fDigitsArr->Clear("C");
3119d479 341 switch (fInputCellType) {
342
343 case kFEEData :
507f74bc 344 case kFEEDataMCOnly :
345 case kFEEDataExcludeMC :
3119d479 346 {
3119d479 347 Double_t avgE = 0; // for background subtraction
507f74bc 348 const Int_t ncells = fCaloCells->GetNumberOfCells();
3119d479 349 for (Int_t icell = 0, idigit = 0; icell < ncells; ++icell) {
350 Double_t cellAmplitude=0, cellTime=0, cellEFrac = 0;
60d77596 351 Short_t cellNumber=0;
352 Int_t cellMCLabel=-1;
507f74bc 353 if (fCaloCells->GetCell(icell, cellNumber, cellAmplitude, cellTime, cellMCLabel, cellEFrac) != kTRUE)
354 break;
5e1e18d1 355
356 if (cellMCLabel > 0 && cellEFrac < 1e-6)
357 cellEFrac = 1;
358
359 if (cellAmplitude < 1e-6 || cellNumber < 0)
360 continue;
361
fde82e42 362 if (fInputCellType == kFEEDataMCOnly) {
363 if (cellMCLabel <= 0)
364 continue;
365 else {
366 cellAmplitude *= cellEFrac;
367 cellEFrac = 1;
368 }
369 }
370 else if (fInputCellType == kFEEDataExcludeMC) {
371 if (cellMCLabel > 0)
372 continue;
5e1e18d1 373 else {
fde82e42 374 cellAmplitude *= 1 - cellEFrac;
5e1e18d1 375 cellEFrac = 0;
376 }
fde82e42 377 }
eff32c61 378
379 if(!AcceptCell(cellNumber)) continue;
380
507f74bc 381 AliEMCALDigit *digit = new((*fDigitsArr)[idigit]) AliEMCALDigit(cellMCLabel, cellMCLabel, cellNumber,
382 (Float_t)cellAmplitude, (Float_t)cellTime,
383 AliEMCALDigit::kHG,idigit, 0, 0, cellEFrac*cellAmplitude);
384
385 if (!fDoClusterize||fSubBackground) {
3119d479 386 Float_t energy = cellAmplitude;
387 Float_t time = cellTime;
388 fClusterizer->Calibrate(energy,time,cellNumber);
389 digit->SetAmplitude(energy);
507f74bc 390 avgE += energy;
fde82e42 391 }
3119d479 392 idigit++;
95f69c66 393 }
fde82e42 394
3119d479 395 if (fSubBackground) {
507f74bc 396 avgE /= fGeom->GetNumberOfSuperModules()*48*24;
3119d479 397 Int_t ndigis = fDigitsArr->GetEntries();
398 for (Int_t i = 0; i < ndigis; ++i) {
507f74bc 399 AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(i));
400 Double_t energy = digit->GetAmplitude() - avgE;
401 if (energy<=0.001) {
402 digit->SetAmplitude(0);
403 } else {
404 digit->SetAmplitude(energy);
405 }
3119d479 406 }
407 }
408 }
409 break;
410
411 case kPattern :
412 {
413 // Fill digits from a pattern
507f74bc 414 Int_t maxd = fGeom->GetNCells() / 4;
3119d479 415 for (Int_t idigit = 0; idigit < maxd; idigit++){
416 if (idigit % 24 == 12) idigit += 12;
417 AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->New(idigit));
418 digit->SetId(idigit * 4);
419 digit->SetTime(600);
420 digit->SetTimeR(600);
421 digit->SetIndexInList(idigit);
422 digit->SetType(AliEMCALDigit::kHG);
423 digit->SetAmplitude(0.1);
424 }
425 }
426 break;
427
428 case kL0FastORs :
429 case kL0FastORsTC :
430 case kL1FastORs :
431 {
432 // Fill digits from FastORs
433
434 AliVCaloTrigger *triggers = InputEvent()->GetCaloTrigger("EMCAL");
435
436 if (!triggers || !(triggers->GetEntries() > 0))
437 return;
438
439 Int_t idigit = 0;
440 triggers->Reset();
441
442 while ((triggers->Next())) {
443 Float_t L0Amplitude = 0;
444 triggers->GetAmplitude(L0Amplitude);
445
446 if (L0Amplitude <= 0 && fInputCellType != kL1FastORs)
447 continue;
448
449 Int_t L1Amplitude = 0;
450 triggers->GetL1TimeSum(L1Amplitude);
451
452 if (L1Amplitude <= 0 && fInputCellType == kL1FastORs)
453 continue;
454
455 Int_t triggerTime = 0;
456 Int_t ntimes = 0;
457 triggers->GetNL0Times(ntimes);
458
459 if (ntimes < 1 && fInputCellType == kL0FastORsTC)
460 continue;
461
462 if (ntimes > 0) {
463 Int_t trgtimes[25];
464 triggers->GetL0Times(trgtimes);
465 triggerTime = trgtimes[0];
466 }
467
468 Int_t triggerCol = 0, triggerRow = 0;
469 triggers->GetPosition(triggerCol, triggerRow);
470
471 Int_t find = -1;
507f74bc 472 fGeom->GetAbsFastORIndexFromPositionInEMCAL(triggerCol, triggerRow, find);
3119d479 473
474 if (find < 0)
475 continue;
476
477 Int_t cidx[4] = {-1};
507f74bc 478 Bool_t ret = fGeom->GetCellIndexFromFastORIndex(find, cidx);
3119d479 479
480 if (!ret)
481 continue;
482
483 Float_t triggerAmplitude = 0;
484
485 if (fInputCellType == kL1FastORs) {
486 triggerAmplitude = 0.25 * L1Amplitude; // it will add 4 cells for 1 amplitude
487 }
488 else {
489 triggerAmplitude = L0Amplitude; // 10 bit truncated, so it is already divided by 4
490 }
491
492 for (Int_t idxpos = 0; idxpos < 4; idxpos++) {
493 Int_t triggerNumber = cidx[idxpos];
494 AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->New(idigit));
495 digit->SetId(triggerNumber);
496 digit->SetTime(triggerTime);
497 digit->SetTimeR(triggerTime);
498 digit->SetIndexInList(idigit);
499 digit->SetType(AliEMCALDigit::kHG);
500 digit->SetAmplitude(triggerAmplitude);
501 idigit++;
502 }
7a7d5a64 503 }
504 }
3119d479 505 break;
06ca320f 506 }
2f7259cf 507}
508
eff32c61 509//________________________________________________________________________________________
510Bool_t AliAnalysisTaskEMCALClusterizeFast::AcceptCell(Int_t cellNumber) {
511
512 Bool_t accept = kTRUE;
513 if(fRejectExoticCells) {
514 //Remove exotic cells before making digits
a3e564aa 515 Bool_t exRemoval = fRecoUtils->IsRejectExoticCell();
eff32c61 516 fRecoUtils->SwitchOnRejectExoticCell();//switch on and off
517 Int_t bunchCrossNo = InputEvent()->GetBunchCrossNumber();
518 Bool_t isEx = fRecoUtils->IsExoticCell(cellNumber, fCaloCells, bunchCrossNo);
e4d5da1e 519 if(isEx) accept = kFALSE;
a3e564aa 520 if(!exRemoval) fRecoUtils->SwitchOffRejectExoticCell();//switch on and off
eff32c61 521 }
522 return accept;
523}
524
2f7259cf 525//________________________________________________________________________________________
507f74bc 526void AliAnalysisTaskEMCALClusterizeFast::CalibrateClusters()
4cf9204b 527{
a9edb9a3 528 // Go through clusters one by one and process separate correction
507f74bc 529 // as those were defined or not
a9edb9a3 530
507f74bc 531 Int_t nclusters = fCaloClusters->GetEntriesFast();
532 for (Int_t icluster=0; icluster < nclusters; ++icluster) {
533 AliVCluster *clust = static_cast<AliVCluster*>(fCaloClusters->At(icluster));
534 if (!clust)
535 continue;
536 if (!clust->IsEMCAL())
537 continue;
4cf9204b 538
507f74bc 539 // REMOVE CLUSTERS WITH BAD CELLS -----------------------------
540 if (fClusterBadChannelCheck) {
541 // careful, the the ClusterContainsBadChannel is dependent on
542 // SwitchOnBadChannelsRemoval, switching it ON automatically
543 // and returning to original value after processing
544 Bool_t badRemoval = fRecoUtils->IsBadChannelsRemovalSwitchedOn();
545 fRecoUtils->SwitchOnBadChannelsRemoval();
546
547 Bool_t badResult = fRecoUtils->ClusterContainsBadChannel(fGeom, clust->GetCellsAbsId(), clust->GetNCells());
548
549 // switch the bad channels removal back
550 if (!badRemoval)
551 fRecoUtils->SwitchOffBadChannelsRemoval();
552
553 if (badResult) {
554 delete fCaloClusters->RemoveAt(icluster);
555 continue; //TODO is it really needed to remove it? Or should we flag it?
556 }
557 }
558
559 // REMOVE EXOTIC CLUSTERS -------------------------------------
560 // does process local cell recalibration energy and time without replacing
561 // the global cell values, in case of no cell recalib done yet
562 if (fRejectExoticClusters) {
563 // careful, the IsExoticCluster is dependent on
564 // SwitchOnRejectExoticCell, switching it ON automatically
565 // and returning to original value after processing
566 Bool_t exRemoval = fRecoUtils->IsRejectExoticCell();
567 fRecoUtils->SwitchOnRejectExoticCell();
568
569 // get bunch crossing
570 Int_t bunchCrossNo = InputEvent()->GetBunchCrossNumber();
571
572 Bool_t exResult = fRecoUtils->IsExoticCluster(clust, fCaloCells, bunchCrossNo);
573
574 // switch the exotic channels removal back
575 if (!exRemoval)
576 fRecoUtils->SwitchOffRejectExoticCell();
577
578 if (exResult) {
579 delete fCaloClusters->RemoveAt(icluster);
580 continue; //TODO is it really needed to remove it? Or should we flag it?
581 }
582 }
583
584 // FIDUCIAL CUT -----------------------------------------------
585 if (fFiducial) {
586 // depends on SetNumberOfCellsFromEMCALBorder
587 // SwitchOnNoFiducialBorderInEMCALEta0
588 if (!fRecoUtils->CheckCellFiducialRegion(fGeom, clust, fCaloCells)){
589 delete fCaloClusters->RemoveAt(icluster);
590 continue; //TODO it would be nice to store the distance
591 }
592 }
593
594 // NONLINEARITY -----------------------------------------------
595 if (fDoNonLinearity) {
596 Float_t correctedEnergy = fRecoUtils->CorrectClusterEnergyLinearity(clust);
597 clust->SetE(correctedEnergy);
598 }
0408f54f 599
507f74bc 600 // DISTANCE TO BAD CHANNELS -----------------------------------
601 if (fRecalDistToBadChannels)
602 fRecoUtils->RecalculateClusterDistanceToBadChannel(fGeom, fCaloCells, clust);
603 }
604
605 fCaloClusters->Compress();
606}
607
608//________________________________________________________________________________________
609void AliAnalysisTaskEMCALClusterizeFast::TrackClusterMatching(AliVCluster *c, TClonesArray *tarr)
610{
611 Float_t g[3]={0};
612 c->GetPosition(g);
613 TVector3 gpos(g);
614
615 Double_t dEtaMin = 1e9;
616 Double_t dPhiMin = 1e9;
617 Int_t imin = -1;
618 Double_t ceta = gpos.Eta();
619 Double_t cphi = gpos.Phi();
620 const Int_t ntrks = tarr->GetEntries();
621 for(Int_t t = 0; t<ntrks; ++t) {
622 AliVTrack *track = static_cast<AliVTrack*>(tarr->At(t));
623 if (!track)
624 continue;
625 const AliExternalTrackParam *outp = track->GetOuterParam();
626 if (!outp)
627 continue;
628 Double_t trkPos[3] = {0.,0.,0.};
629 if (!outp->GetXYZ(trkPos))
630 continue;
631 TVector3 vec(trkPos[0],trkPos[1],trkPos[2]);
632 Double_t veta = vec.Eta();
633 Double_t vphi = vec.Phi();
634 if(vphi<0)
635 vphi += 2*TMath::Pi();
636 if (TMath::Abs(veta)>0.75 || (vphi<70*TMath::DegToRad()) || (vphi>190*TMath::DegToRad()))
637 continue;
638 Double_t dR = vec.DeltaR(gpos);
639 if(dR > 25)
640 continue;
641 Float_t tmpEta=0, tmpPhi=0;
642 if (0) {
643 AliExternalTrackParam trkParTemp(*outp); // retrieve the starting point every time before the extrapolation
644 Bool_t ret = fRecoUtils->ExtrapolateTrackToCluster(&trkParTemp, c, fRecoUtils->GetMass(), fRecoUtils->GetStep(), tmpEta, tmpPhi);
645 if (!ret)
646 continue;
647 } else {
648 tmpEta = ceta - veta;
649 tmpPhi = cphi - vphi;
650 }
651 if (TMath::Abs(tmpEta)<TMath::Abs(dEtaMin) && TMath::Abs(tmpPhi)<TMath::Abs(dPhiMin)) {
652 dEtaMin = tmpEta;
653 dPhiMin = tmpPhi;
654 imin = t;
655 }
656 }
657 c->SetEmcCpvDistance(imin);
658 c->SetTrackDistance(dPhiMin, dEtaMin);
659}
660
661//________________________________________________________________________________________
662void AliAnalysisTaskEMCALClusterizeFast::RecPoints2Clusters(TClonesArray *clus)
663{
664 // Cluster energy, global position, cells and their amplitude fractions are restored.
665
0408f54f 666 // tracks array for track/cluster matching
667 TClonesArray *tarr = 0;
668 if (!fTrackName.IsNull()) {
507f74bc 669 tarr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTrackName));
0408f54f 670 if (!tarr) {
671 AliError(Form("Cannot get tracks named %s", fTrackName.Data()));
672 }
673 }
fde82e42 674
eb33b7f8 675 const Int_t Ncls = fClusterArr->GetEntries();
0408f54f 676 AliDebug(1, Form("total no of clusters %d", Ncls));
4cf9204b 677 for(Int_t i=0, nout=clus->GetEntries(); i < Ncls; ++i) {
678 AliEMCALRecPoint *recpoint = static_cast<AliEMCALRecPoint*>(fClusterArr->At(i));
679 Int_t ncells_true = 0;
680 const Int_t ncells = recpoint->GetMultiplicity();
681 UShort_t absIds[ncells];
682 Double32_t ratios[ncells];
683 Int_t *dlist = recpoint->GetDigitsList();
684 Float_t *elist = recpoint->GetEnergiesList();
7eff06f2 685 Double_t mcEnergy = 0;
4cf9204b 686 for (Int_t c = 0; c < ncells; ++c) {
687 AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(dlist[c]));
688 absIds[ncells_true] = digit->GetId();
5ded2286 689 ratios[ncells_true] = elist[c]/recpoint->GetEnergy();
7eff06f2 690 if (digit->GetIparent(1) > 0)
691 mcEnergy += digit->GetDEParent(1)/recpoint->GetEnergy();
84ccad86 692 ++ncells_true;
4cf9204b 693 }
694
695 if (ncells_true < 1) {
696 AliWarning("Skipping cluster with no cells");
697 continue;
698 }
699
700 // calculate new cluster position
701 TVector3 gpos;
4cf9204b 702 recpoint->GetGlobalPosition(gpos);
84ccad86 703 Float_t g[3];
4cf9204b 704 gpos.GetXYZ(g);
705
e93ec1f7 706 AliDebug(1, Form("energy %f", recpoint->GetEnergy()));
707
95f69c66 708 AliVCluster *c = static_cast<AliVCluster*>(clus->New(nout++));
4cf9204b 709 c->SetType(AliVCluster::kEMCALClusterv1);
710 c->SetE(recpoint->GetEnergy());
711 c->SetPosition(g);
712 c->SetNCells(ncells_true);
a55e4f1d 713 c->SetCellsAbsId(absIds);
714 c->SetCellsAmplitudeFraction(ratios);
715 c->SetID(recpoint->GetUniqueID());
4cf9204b 716 c->SetDispersion(recpoint->GetDispersion());
0408f54f 717 c->SetEmcCpvDistance(-1);
718 c->SetChi2(-1);
4cf9204b 719 c->SetTOF(recpoint->GetTime()) ; //time-of-flight
720 c->SetNExMax(recpoint->GetNExMax()); //number of local maxima
721 Float_t elipAxis[2];
722 recpoint->GetElipsAxis(elipAxis);
e93ec1f7 723 c->SetM02(elipAxis[0]*elipAxis[0]);
724 c->SetM20(elipAxis[1]*elipAxis[1]);
7eff06f2 725 c->SetMCEnergyFraction(mcEnergy);
507f74bc 726
727 //MC
728 AliESDCaloCluster *esdClus = dynamic_cast<AliESDCaloCluster*>(c);
729 if (esdClus) {
730 Int_t parentMult = 0;
731 Int_t *parentList = recpoint->GetParents(parentMult);
732 if (parentMult > 0) {
733 TArrayI parents(parentMult, parentList);
734 esdClus->AddLabels(parents);
735 }
736 }
737 else {
738 AliAODCaloCluster *aodClus = dynamic_cast<AliAODCaloCluster*>(c);
739 if (aodClus) {
740 Int_t parentMult = 0;
741 Int_t *parentList = recpoint->GetParents(parentMult);
742 aodClus->SetLabel(parentList, parentMult);
95f69c66 743 }
2f7259cf 744 }
507f74bc 745
746 if (tarr)
747 TrackClusterMatching(c, tarr);
2f7259cf 748 }
749}
750
95f69c66 751//________________________________________________________________________
752void AliAnalysisTaskEMCALClusterizeFast::UpdateCells()
753{
754 // Update cells in case re-calibration was done.
95f69c66 755 if (!fCalibData&&!fSubBackground)
756 return;
757
507f74bc 758 const Int_t ncells = fCaloCells->GetNumberOfCells();
55a8e3e4 759 const Int_t ndigis = fDigitsArr->GetEntries();
95f69c66 760 if (ncells!=ndigis) {
507f74bc 761 fCaloCells->DeleteContainer();
762 fCaloCells->CreateContainer(ndigis);
95f69c66 763 }
764 for (Int_t idigit = 0; idigit < ndigis; ++idigit) {
765 AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(idigit));
766 Double_t cellAmplitude = digit->GetCalibAmp();
55a8e3e4 767 Short_t cellNumber = digit->GetId();
768 Double_t cellTime = digit->GetTime();
507f74bc 769 fCaloCells->SetCell(idigit, cellNumber, cellAmplitude, cellTime);
95f69c66 770 }
771}
772
773//________________________________________________________________________
774void AliAnalysisTaskEMCALClusterizeFast::UpdateClusters()
775{
776 // Update cells in case re-calibration was done.
7a7d5a64 777
507f74bc 778 const Int_t nents = fCaloClusters->GetEntries();
779 for (Int_t i=0;i<nents;++i) {
780 AliVCluster *c = static_cast<AliVCluster*>(fCaloClusters->At(i));
781 if (!c)
782 continue;
783 if (c->IsEMCAL())
784 delete fCaloClusters->RemoveAt(i);
7a7d5a64 785 }
507f74bc 786
787 fCaloClusters->Compress();
7a7d5a64 788
507f74bc 789 RecPoints2Clusters(fCaloClusters);
95f69c66 790}
791
2f7259cf 792//________________________________________________________________________________________
793void AliAnalysisTaskEMCALClusterizeFast::Init()
794{
4f66d2ba 795 // Select clusterization/unfolding algorithm and set all the needed parameters.
01bc1ce8 796
507f74bc 797 if (InputEvent()->GetRunNumber()==fRun)
2f7259cf 798 return;
507f74bc 799 fRun = InputEvent()->GetRunNumber();
2f7259cf 800
801 if (fJustUnfold){
802 // init the unfolding afterburner
803 delete fUnfolder;
1759f743 804 fUnfolder = new AliEMCALAfterBurnerUF(fRecParam->GetW0(),fRecParam->GetLocMaxCut(),fRecParam->GetMinECut());
2f7259cf 805 return;
806 }
807
4f66d2ba 808 if (fGeomName.Length()>0)
507f74bc 809 fGeom = AliEMCALGeometry::GetInstance(fGeomName);
4f66d2ba 810 else
507f74bc 811 fGeom = AliEMCALGeometry::GetInstance();
812 if (!fGeom) {
2f7259cf 813 AliFatal("Geometry not available!!!");
814 return;
815 }
816
817 if (!fGeomMatrixSet) {
818 if (fLoadGeomMatrices) {
507f74bc 819 for(Int_t mod=0; mod < fGeom->GetNumberOfSuperModules(); ++mod) {
06ca320f 820 if (fGeomMatrix[mod]){
821 if (DebugLevel() > 2)
2f7259cf 822 fGeomMatrix[mod]->Print();
507f74bc 823 fGeom->SetMisalMatrix(fGeomMatrix[mod],mod);
2f7259cf 824 }
825 }
29b496d5 826 } else { // get matrix from file (work around bug in aliroot)
507f74bc 827 for(Int_t mod=0; mod < fGeom->GetEMCGeometry()->GetNumberOfSuperModules(); ++mod) {
29b496d5 828 const TGeoHMatrix *gm = 0;
507f74bc 829 if (fEsd) {
830 gm = fEsd->GetEMCALMatrix(mod);
29b496d5 831 } else {
507f74bc 832 AliAODHeader *aodheader = fAod->GetHeader();
29b496d5 833 if (aodheader) {
834 gm = aodheader->GetEMCALMatrix(mod);
835 }
836 }
837 if (gm) {
06ca320f 838 if (DebugLevel() > 2)
29b496d5 839 gm->Print();
507f74bc 840 fGeom->SetMisalMatrix(gm,mod);
2f7259cf 841 }
842 }
843 }
844 fGeomMatrixSet=kTRUE;
845 }
846
847 // setup digit array if needed
848 if (!fDigitsArr) {
849 fDigitsArr = new TClonesArray("AliEMCALDigit", 1000);
850 fDigitsArr->SetOwner(1);
851 }
852
853 // then setup clusterizer
cac0e10b 854 if (fClusterizer) {
855 // avoid to delete digits array
856 fClusterizer->SetDigitsArr(0);
857 delete fClusterizer;
858 }
507f74bc 859 if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv1)
860 fClusterizer = new AliEMCALClusterizerv1(fGeom);
06ca320f 861 else if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerNxN) {
80dbd4f7 862 AliEMCALClusterizerNxN *clusterizer = new AliEMCALClusterizerNxN(fGeom);
863 clusterizer->SetNRowDiff(fRecParam->GetNRowDiff()); //MV: already done in AliEMCALClusterizer::InitParameters
864 clusterizer->SetNColDiff(fRecParam->GetNColDiff()); //MV: already done in AliEMCALClusterizer::InitParameters
3c56da8f 865 fClusterizer = clusterizer;
507f74bc 866 }
867 else if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv2)
868 fClusterizer = new AliEMCALClusterizerv2(fGeom);
869 else if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerFW) {
870 AliEMCALClusterizerFixedWindow *clusterizer = new AliEMCALClusterizerFixedWindow(fGeom);
0fda97fe 871 clusterizer->SetNphi(fNPhi);
872 clusterizer->SetNeta(fNEta);
873 clusterizer->SetShiftPhi(fShiftPhi);
874 clusterizer->SetShiftEta(fShiftEta);
3c56da8f 875 clusterizer->SetTRUshift(fTRUShift);
7a7d5a64 876 fClusterizer = clusterizer;
877 }
507f74bc 878 else {
2f7259cf 879 AliFatal(Form("Clusterizer < %d > not available", fRecParam->GetClusterizerFlag()));
880 }
881 fClusterizer->InitParameters(fRecParam);
01bc1ce8 882
883 if ((!fCalibData&&fLoadCalib) || (!fPedestalData&&fLoadPed)) {
884 AliCDBManager *cdb = AliCDBManager::Instance();
885 if (!cdb->IsDefaultStorageSet() && !fOCDBpath.IsNull())
886 cdb->SetDefaultStorage(fOCDBpath);
887 if (fRun!=cdb->GetRun())
888 cdb->SetRun(fRun);
889 }
349bf69e 890 if (!fCalibData&&fLoadCalib&&fRun>0) {
2f7259cf 891 AliCDBEntry *entry = static_cast<AliCDBEntry*>(AliCDBManager::Instance()->Get("EMCAL/Calib/Data"));
892 if (entry)
893 fCalibData = static_cast<AliEMCALCalibData*>(entry->GetObject());
894 if (!fCalibData)
895 AliFatal("Calibration parameters not found in CDB!");
896 }
349bf69e 897 if (!fPedestalData&&fLoadPed&&fRun>0) {
2f7259cf 898 AliCDBEntry *entry = static_cast<AliCDBEntry*>(AliCDBManager::Instance()->Get("EMCAL/Calib/Pedestals"));
899 if (entry)
900 fPedestalData = static_cast<AliCaloCalibPedestal*>(entry->GetObject());
901 }
10d33986 902 if (fCalibData) {
903 fClusterizer->SetInputCalibrated(kFALSE);
904 fClusterizer->SetCalibrationParameters(fCalibData);
10d33986 905 } else {
906 fClusterizer->SetInputCalibrated(kTRUE);
907 }
e93ec1f7 908 fClusterizer->SetCaloCalibPedestal(fPedestalData);
10d33986 909 fClusterizer->SetJustClusters(kTRUE);
2f7259cf 910 fClusterizer->SetDigitsArr(fDigitsArr);
911 fClusterizer->SetOutput(0);
912 fClusterArr = const_cast<TObjArray *>(fClusterizer->GetRecPoints());
507f74bc 913
914 // Get the emcal cells
fde82e42 915 if ((fInputCellType == kFEEData || fInputCellType == kFEEDataMCOnly || fInputCellType == kFEEDataExcludeMC) && !fCaloCells) {
507f74bc 916 if (fCaloCellsName.IsNull()) {
917 fCaloCells = InputEvent()->GetEMCALCells();
918 }
919 else {
920 fCaloCells = dynamic_cast<AliVCaloCells*>(InputEvent()->FindListObject(fCaloCellsName));
921 if (!fCaloCells)
922 AliError(Form("%s: Could not retrieve cells %s!", GetName(), fCaloCellsName.Data()));
923 }
924 if (!fCaloCells)
925 AliFatal("Could not get EMCal cells!");
926 }
927
928 // Set output clusters collection
929 if (!fAttachClusters) {
930 fCaloClusters = fOutputAODBranch;
931 return;
932 }
933
934 if (!fCaloClusters) {
935 if (fCaloClustersName.IsNull()) { //overwrite mode
936 if (fEsd)
937 fCaloClusters = static_cast<TClonesArray*>(InputEvent()->FindListObject("CaloClusters"));
938 else if (fAod)
939 fCaloClusters = static_cast<TClonesArray*>(InputEvent()->FindListObject("caloClusters"));
940 }
941 else {
942 fCaloClusters = static_cast<TClonesArray*>(InputEvent()->FindListObject(fCaloClustersName));
943
944 if (!fCaloClusters) {
945 if (fEsd)
946 fCaloClusters = new TClonesArray("AliESDCaloCluster");
947 else if (fAod)
948 fCaloClusters = new TClonesArray("AliAODCaloCluster");
949
950 fCaloClusters->SetName(fCaloClustersName);
951 InputEvent()->AddObject(fCaloClusters);
952 }
953 }
954
955 if (!fCaloClusters)
956 AliFatal("Could not get cluster collection!");
957
958 TClass *cl = fCaloClusters->GetClass();
959 if (!cl->GetBaseClass("AliVCluster")) {
960 AliFatal(Form("%s: Collection %s does not contain AliVCluster objects!", GetName(), fCaloClusters->GetName()));
961 fCaloClusters = 0;
962 return;
963 }
964 }
2f7259cf 965}