]> 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
add TPC-only track cuts support
[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"),
0a9eac09 91 fDoUpdateCells(kFALSE),
648ade84 92 fDoClusterize(kTRUE),
0a9eac09 93 fClusterBadChannelCheck(kFALSE),
507f74bc 94 fRejectExoticClusters(kFALSE),
eff32c61 95 fRejectExoticCells(kFALSE),
507f74bc 96 fFiducial(kFALSE),
97 fDoNonLinearity(kFALSE),
0a9eac09 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"),
0a9eac09 143 fDoUpdateCells(kFALSE),
648ade84 144 fDoClusterize(kTRUE),
0a9eac09 145 fClusterBadChannelCheck(kFALSE),
507f74bc 146 fRejectExoticClusters(kFALSE),
eff32c61 147 fRejectExoticCells(kFALSE),
507f74bc 148 fFiducial(kFALSE),
149 fDoNonLinearity(kFALSE),
0a9eac09 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
515 fRecoUtils->SwitchOnRejectExoticCell();//switch on and off
516 Int_t bunchCrossNo = InputEvent()->GetBunchCrossNumber();
517 Bool_t isEx = fRecoUtils->IsExoticCell(cellNumber, fCaloCells, bunchCrossNo);
648ade84 518 if(isEx) accept = kFALSE;
eff32c61 519 fRecoUtils->SwitchOffRejectExoticCell();//switch on and off
520 }
521 return accept;
522}
523
2f7259cf 524//________________________________________________________________________________________
507f74bc 525void AliAnalysisTaskEMCALClusterizeFast::CalibrateClusters()
4cf9204b 526{
a9edb9a3 527 // Go through clusters one by one and process separate correction
507f74bc 528 // as those were defined or not
a9edb9a3 529
507f74bc 530 Int_t nclusters = fCaloClusters->GetEntriesFast();
531 for (Int_t icluster=0; icluster < nclusters; ++icluster) {
532 AliVCluster *clust = static_cast<AliVCluster*>(fCaloClusters->At(icluster));
533 if (!clust)
534 continue;
535 if (!clust->IsEMCAL())
536 continue;
4cf9204b 537
507f74bc 538 // REMOVE CLUSTERS WITH BAD CELLS -----------------------------
539 if (fClusterBadChannelCheck) {
540 // careful, the the ClusterContainsBadChannel is dependent on
541 // SwitchOnBadChannelsRemoval, switching it ON automatically
542 // and returning to original value after processing
543 Bool_t badRemoval = fRecoUtils->IsBadChannelsRemovalSwitchedOn();
544 fRecoUtils->SwitchOnBadChannelsRemoval();
545
546 Bool_t badResult = fRecoUtils->ClusterContainsBadChannel(fGeom, clust->GetCellsAbsId(), clust->GetNCells());
547
548 // switch the bad channels removal back
549 if (!badRemoval)
550 fRecoUtils->SwitchOffBadChannelsRemoval();
551
552 if (badResult) {
553 delete fCaloClusters->RemoveAt(icluster);
554 continue; //TODO is it really needed to remove it? Or should we flag it?
555 }
556 }
557
558 // REMOVE EXOTIC CLUSTERS -------------------------------------
559 // does process local cell recalibration energy and time without replacing
560 // the global cell values, in case of no cell recalib done yet
561 if (fRejectExoticClusters) {
562 // careful, the IsExoticCluster is dependent on
563 // SwitchOnRejectExoticCell, switching it ON automatically
564 // and returning to original value after processing
565 Bool_t exRemoval = fRecoUtils->IsRejectExoticCell();
566 fRecoUtils->SwitchOnRejectExoticCell();
567
568 // get bunch crossing
569 Int_t bunchCrossNo = InputEvent()->GetBunchCrossNumber();
570
571 Bool_t exResult = fRecoUtils->IsExoticCluster(clust, fCaloCells, bunchCrossNo);
572
573 // switch the exotic channels removal back
574 if (!exRemoval)
575 fRecoUtils->SwitchOffRejectExoticCell();
576
577 if (exResult) {
578 delete fCaloClusters->RemoveAt(icluster);
579 continue; //TODO is it really needed to remove it? Or should we flag it?
580 }
581 }
582
583 // FIDUCIAL CUT -----------------------------------------------
584 if (fFiducial) {
585 // depends on SetNumberOfCellsFromEMCALBorder
586 // SwitchOnNoFiducialBorderInEMCALEta0
587 if (!fRecoUtils->CheckCellFiducialRegion(fGeom, clust, fCaloCells)){
588 delete fCaloClusters->RemoveAt(icluster);
589 continue; //TODO it would be nice to store the distance
590 }
591 }
592
593 // NONLINEARITY -----------------------------------------------
594 if (fDoNonLinearity) {
595 Float_t correctedEnergy = fRecoUtils->CorrectClusterEnergyLinearity(clust);
596 clust->SetE(correctedEnergy);
597 }
0408f54f 598
507f74bc 599 // DISTANCE TO BAD CHANNELS -----------------------------------
600 if (fRecalDistToBadChannels)
601 fRecoUtils->RecalculateClusterDistanceToBadChannel(fGeom, fCaloCells, clust);
602 }
603
604 fCaloClusters->Compress();
605}
606
607//________________________________________________________________________________________
608void AliAnalysisTaskEMCALClusterizeFast::TrackClusterMatching(AliVCluster *c, TClonesArray *tarr)
609{
610 Float_t g[3]={0};
611 c->GetPosition(g);
612 TVector3 gpos(g);
613
614 Double_t dEtaMin = 1e9;
615 Double_t dPhiMin = 1e9;
616 Int_t imin = -1;
617 Double_t ceta = gpos.Eta();
618 Double_t cphi = gpos.Phi();
619 const Int_t ntrks = tarr->GetEntries();
620 for(Int_t t = 0; t<ntrks; ++t) {
621 AliVTrack *track = static_cast<AliVTrack*>(tarr->At(t));
622 if (!track)
623 continue;
624 const AliExternalTrackParam *outp = track->GetOuterParam();
625 if (!outp)
626 continue;
627 Double_t trkPos[3] = {0.,0.,0.};
628 if (!outp->GetXYZ(trkPos))
629 continue;
630 TVector3 vec(trkPos[0],trkPos[1],trkPos[2]);
631 Double_t veta = vec.Eta();
632 Double_t vphi = vec.Phi();
633 if(vphi<0)
634 vphi += 2*TMath::Pi();
635 if (TMath::Abs(veta)>0.75 || (vphi<70*TMath::DegToRad()) || (vphi>190*TMath::DegToRad()))
636 continue;
637 Double_t dR = vec.DeltaR(gpos);
638 if(dR > 25)
639 continue;
640 Float_t tmpEta=0, tmpPhi=0;
641 if (0) {
642 AliExternalTrackParam trkParTemp(*outp); // retrieve the starting point every time before the extrapolation
643 Bool_t ret = fRecoUtils->ExtrapolateTrackToCluster(&trkParTemp, c, fRecoUtils->GetMass(), fRecoUtils->GetStep(), tmpEta, tmpPhi);
644 if (!ret)
645 continue;
646 } else {
647 tmpEta = ceta - veta;
648 tmpPhi = cphi - vphi;
649 }
650 if (TMath::Abs(tmpEta)<TMath::Abs(dEtaMin) && TMath::Abs(tmpPhi)<TMath::Abs(dPhiMin)) {
651 dEtaMin = tmpEta;
652 dPhiMin = tmpPhi;
653 imin = t;
654 }
655 }
656 c->SetEmcCpvDistance(imin);
657 c->SetTrackDistance(dPhiMin, dEtaMin);
658}
659
660//________________________________________________________________________________________
661void AliAnalysisTaskEMCALClusterizeFast::RecPoints2Clusters(TClonesArray *clus)
662{
663 // Cluster energy, global position, cells and their amplitude fractions are restored.
664
0408f54f 665 // tracks array for track/cluster matching
666 TClonesArray *tarr = 0;
667 if (!fTrackName.IsNull()) {
507f74bc 668 tarr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTrackName));
0408f54f 669 if (!tarr) {
670 AliError(Form("Cannot get tracks named %s", fTrackName.Data()));
671 }
672 }
fde82e42 673
eb33b7f8 674 const Int_t Ncls = fClusterArr->GetEntries();
0408f54f 675 AliDebug(1, Form("total no of clusters %d", Ncls));
4cf9204b 676 for(Int_t i=0, nout=clus->GetEntries(); i < Ncls; ++i) {
677 AliEMCALRecPoint *recpoint = static_cast<AliEMCALRecPoint*>(fClusterArr->At(i));
678 Int_t ncells_true = 0;
679 const Int_t ncells = recpoint->GetMultiplicity();
680 UShort_t absIds[ncells];
681 Double32_t ratios[ncells];
682 Int_t *dlist = recpoint->GetDigitsList();
683 Float_t *elist = recpoint->GetEnergiesList();
7eff06f2 684 Double_t mcEnergy = 0;
4cf9204b 685 for (Int_t c = 0; c < ncells; ++c) {
686 AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(dlist[c]));
687 absIds[ncells_true] = digit->GetId();
5ded2286 688 ratios[ncells_true] = elist[c]/recpoint->GetEnergy();
7eff06f2 689 if (digit->GetIparent(1) > 0)
690 mcEnergy += digit->GetDEParent(1)/recpoint->GetEnergy();
84ccad86 691 ++ncells_true;
4cf9204b 692 }
693
694 if (ncells_true < 1) {
695 AliWarning("Skipping cluster with no cells");
696 continue;
697 }
698
699 // calculate new cluster position
700 TVector3 gpos;
4cf9204b 701 recpoint->GetGlobalPosition(gpos);
84ccad86 702 Float_t g[3];
4cf9204b 703 gpos.GetXYZ(g);
704
e93ec1f7 705 AliDebug(1, Form("energy %f", recpoint->GetEnergy()));
706
95f69c66 707 AliVCluster *c = static_cast<AliVCluster*>(clus->New(nout++));
4cf9204b 708 c->SetType(AliVCluster::kEMCALClusterv1);
709 c->SetE(recpoint->GetEnergy());
710 c->SetPosition(g);
711 c->SetNCells(ncells_true);
a55e4f1d 712 c->SetCellsAbsId(absIds);
713 c->SetCellsAmplitudeFraction(ratios);
714 c->SetID(recpoint->GetUniqueID());
4cf9204b 715 c->SetDispersion(recpoint->GetDispersion());
0408f54f 716 c->SetEmcCpvDistance(-1);
717 c->SetChi2(-1);
4cf9204b 718 c->SetTOF(recpoint->GetTime()) ; //time-of-flight
719 c->SetNExMax(recpoint->GetNExMax()); //number of local maxima
720 Float_t elipAxis[2];
721 recpoint->GetElipsAxis(elipAxis);
e93ec1f7 722 c->SetM02(elipAxis[0]*elipAxis[0]);
723 c->SetM20(elipAxis[1]*elipAxis[1]);
7eff06f2 724 c->SetMCEnergyFraction(mcEnergy);
507f74bc 725
726 //MC
727 AliESDCaloCluster *esdClus = dynamic_cast<AliESDCaloCluster*>(c);
728 if (esdClus) {
729 Int_t parentMult = 0;
730 Int_t *parentList = recpoint->GetParents(parentMult);
731 if (parentMult > 0) {
732 TArrayI parents(parentMult, parentList);
733 esdClus->AddLabels(parents);
734 }
735 }
736 else {
737 AliAODCaloCluster *aodClus = dynamic_cast<AliAODCaloCluster*>(c);
738 if (aodClus) {
739 Int_t parentMult = 0;
740 Int_t *parentList = recpoint->GetParents(parentMult);
741 aodClus->SetLabel(parentList, parentMult);
95f69c66 742 }
2f7259cf 743 }
507f74bc 744
745 if (tarr)
746 TrackClusterMatching(c, tarr);
2f7259cf 747 }
748}
749
95f69c66 750//________________________________________________________________________
751void AliAnalysisTaskEMCALClusterizeFast::UpdateCells()
752{
753 // Update cells in case re-calibration was done.
95f69c66 754 if (!fCalibData&&!fSubBackground)
755 return;
756
507f74bc 757 const Int_t ncells = fCaloCells->GetNumberOfCells();
55a8e3e4 758 const Int_t ndigis = fDigitsArr->GetEntries();
95f69c66 759 if (ncells!=ndigis) {
507f74bc 760 fCaloCells->DeleteContainer();
761 fCaloCells->CreateContainer(ndigis);
95f69c66 762 }
763 for (Int_t idigit = 0; idigit < ndigis; ++idigit) {
764 AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(idigit));
765 Double_t cellAmplitude = digit->GetCalibAmp();
55a8e3e4 766 Short_t cellNumber = digit->GetId();
767 Double_t cellTime = digit->GetTime();
507f74bc 768 fCaloCells->SetCell(idigit, cellNumber, cellAmplitude, cellTime);
95f69c66 769 }
770}
771
772//________________________________________________________________________
773void AliAnalysisTaskEMCALClusterizeFast::UpdateClusters()
774{
775 // Update cells in case re-calibration was done.
7a7d5a64 776
507f74bc 777 const Int_t nents = fCaloClusters->GetEntries();
778 for (Int_t i=0;i<nents;++i) {
779 AliVCluster *c = static_cast<AliVCluster*>(fCaloClusters->At(i));
780 if (!c)
781 continue;
782 if (c->IsEMCAL())
783 delete fCaloClusters->RemoveAt(i);
7a7d5a64 784 }
507f74bc 785
786 fCaloClusters->Compress();
7a7d5a64 787
507f74bc 788 RecPoints2Clusters(fCaloClusters);
95f69c66 789}
790
2f7259cf 791//________________________________________________________________________________________
792void AliAnalysisTaskEMCALClusterizeFast::Init()
793{
4f66d2ba 794 // Select clusterization/unfolding algorithm and set all the needed parameters.
01bc1ce8 795
507f74bc 796 if (InputEvent()->GetRunNumber()==fRun)
2f7259cf 797 return;
507f74bc 798 fRun = InputEvent()->GetRunNumber();
2f7259cf 799
800 if (fJustUnfold){
801 // init the unfolding afterburner
802 delete fUnfolder;
1759f743 803 fUnfolder = new AliEMCALAfterBurnerUF(fRecParam->GetW0(),fRecParam->GetLocMaxCut(),fRecParam->GetMinECut());
2f7259cf 804 return;
805 }
806
4f66d2ba 807 if (fGeomName.Length()>0)
507f74bc 808 fGeom = AliEMCALGeometry::GetInstance(fGeomName);
4f66d2ba 809 else
507f74bc 810 fGeom = AliEMCALGeometry::GetInstance();
811 if (!fGeom) {
2f7259cf 812 AliFatal("Geometry not available!!!");
813 return;
814 }
815
816 if (!fGeomMatrixSet) {
817 if (fLoadGeomMatrices) {
507f74bc 818 for(Int_t mod=0; mod < fGeom->GetNumberOfSuperModules(); ++mod) {
06ca320f 819 if (fGeomMatrix[mod]){
820 if (DebugLevel() > 2)
2f7259cf 821 fGeomMatrix[mod]->Print();
507f74bc 822 fGeom->SetMisalMatrix(fGeomMatrix[mod],mod);
2f7259cf 823 }
824 }
29b496d5 825 } else { // get matrix from file (work around bug in aliroot)
507f74bc 826 for(Int_t mod=0; mod < fGeom->GetEMCGeometry()->GetNumberOfSuperModules(); ++mod) {
29b496d5 827 const TGeoHMatrix *gm = 0;
507f74bc 828 if (fEsd) {
829 gm = fEsd->GetEMCALMatrix(mod);
29b496d5 830 } else {
507f74bc 831 AliAODHeader *aodheader = fAod->GetHeader();
29b496d5 832 if (aodheader) {
833 gm = aodheader->GetEMCALMatrix(mod);
834 }
835 }
836 if (gm) {
06ca320f 837 if (DebugLevel() > 2)
29b496d5 838 gm->Print();
507f74bc 839 fGeom->SetMisalMatrix(gm,mod);
2f7259cf 840 }
841 }
842 }
843 fGeomMatrixSet=kTRUE;
844 }
845
846 // setup digit array if needed
847 if (!fDigitsArr) {
848 fDigitsArr = new TClonesArray("AliEMCALDigit", 1000);
849 fDigitsArr->SetOwner(1);
850 }
851
852 // then setup clusterizer
cac0e10b 853 if (fClusterizer) {
854 // avoid to delete digits array
855 fClusterizer->SetDigitsArr(0);
856 delete fClusterizer;
857 }
507f74bc 858 if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv1)
859 fClusterizer = new AliEMCALClusterizerv1(fGeom);
06ca320f 860 else if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerNxN) {
98dee1ae 861 AliEMCALClusterizerNxN *clusterizer = new AliEMCALClusterizerNxN(fGeom);
862 clusterizer->SetNRowDiff(fRecParam->GetNRowDiff()); //MV: already done in AliEMCALClusterizer::InitParameters
863 clusterizer->SetNColDiff(fRecParam->GetNColDiff()); //MV: already done in AliEMCALClusterizer::InitParameters
3c56da8f 864 fClusterizer = clusterizer;
507f74bc 865 }
866 else if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv2)
867 fClusterizer = new AliEMCALClusterizerv2(fGeom);
868 else if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerFW) {
869 AliEMCALClusterizerFixedWindow *clusterizer = new AliEMCALClusterizerFixedWindow(fGeom);
0fda97fe 870 clusterizer->SetNphi(fNPhi);
871 clusterizer->SetNeta(fNEta);
872 clusterizer->SetShiftPhi(fShiftPhi);
873 clusterizer->SetShiftEta(fShiftEta);
3c56da8f 874 clusterizer->SetTRUshift(fTRUShift);
7a7d5a64 875 fClusterizer = clusterizer;
876 }
507f74bc 877 else {
2f7259cf 878 AliFatal(Form("Clusterizer < %d > not available", fRecParam->GetClusterizerFlag()));
879 }
880 fClusterizer->InitParameters(fRecParam);
01bc1ce8 881
882 if ((!fCalibData&&fLoadCalib) || (!fPedestalData&&fLoadPed)) {
883 AliCDBManager *cdb = AliCDBManager::Instance();
884 if (!cdb->IsDefaultStorageSet() && !fOCDBpath.IsNull())
885 cdb->SetDefaultStorage(fOCDBpath);
886 if (fRun!=cdb->GetRun())
887 cdb->SetRun(fRun);
888 }
349bf69e 889 if (!fCalibData&&fLoadCalib&&fRun>0) {
2f7259cf 890 AliCDBEntry *entry = static_cast<AliCDBEntry*>(AliCDBManager::Instance()->Get("EMCAL/Calib/Data"));
891 if (entry)
892 fCalibData = static_cast<AliEMCALCalibData*>(entry->GetObject());
893 if (!fCalibData)
894 AliFatal("Calibration parameters not found in CDB!");
895 }
349bf69e 896 if (!fPedestalData&&fLoadPed&&fRun>0) {
2f7259cf 897 AliCDBEntry *entry = static_cast<AliCDBEntry*>(AliCDBManager::Instance()->Get("EMCAL/Calib/Pedestals"));
898 if (entry)
899 fPedestalData = static_cast<AliCaloCalibPedestal*>(entry->GetObject());
900 }
10d33986 901 if (fCalibData) {
902 fClusterizer->SetInputCalibrated(kFALSE);
903 fClusterizer->SetCalibrationParameters(fCalibData);
10d33986 904 } else {
905 fClusterizer->SetInputCalibrated(kTRUE);
906 }
e93ec1f7 907 fClusterizer->SetCaloCalibPedestal(fPedestalData);
10d33986 908 fClusterizer->SetJustClusters(kTRUE);
2f7259cf 909 fClusterizer->SetDigitsArr(fDigitsArr);
910 fClusterizer->SetOutput(0);
911 fClusterArr = const_cast<TObjArray *>(fClusterizer->GetRecPoints());
507f74bc 912
913 // Get the emcal cells
fde82e42 914 if ((fInputCellType == kFEEData || fInputCellType == kFEEDataMCOnly || fInputCellType == kFEEDataExcludeMC) && !fCaloCells) {
507f74bc 915 if (fCaloCellsName.IsNull()) {
916 fCaloCells = InputEvent()->GetEMCALCells();
917 }
918 else {
919 fCaloCells = dynamic_cast<AliVCaloCells*>(InputEvent()->FindListObject(fCaloCellsName));
920 if (!fCaloCells)
921 AliError(Form("%s: Could not retrieve cells %s!", GetName(), fCaloCellsName.Data()));
922 }
923 if (!fCaloCells)
924 AliFatal("Could not get EMCal cells!");
925 }
926
927 // Set output clusters collection
928 if (!fAttachClusters) {
929 fCaloClusters = fOutputAODBranch;
930 return;
931 }
932
933 if (!fCaloClusters) {
934 if (fCaloClustersName.IsNull()) { //overwrite mode
935 if (fEsd)
936 fCaloClusters = static_cast<TClonesArray*>(InputEvent()->FindListObject("CaloClusters"));
937 else if (fAod)
938 fCaloClusters = static_cast<TClonesArray*>(InputEvent()->FindListObject("caloClusters"));
939 }
940 else {
941 fCaloClusters = static_cast<TClonesArray*>(InputEvent()->FindListObject(fCaloClustersName));
942
943 if (!fCaloClusters) {
944 if (fEsd)
945 fCaloClusters = new TClonesArray("AliESDCaloCluster");
946 else if (fAod)
947 fCaloClusters = new TClonesArray("AliAODCaloCluster");
948
949 fCaloClusters->SetName(fCaloClustersName);
950 InputEvent()->AddObject(fCaloClusters);
951 }
952 }
953
954 if (!fCaloClusters)
955 AliFatal("Could not get cluster collection!");
956
957 TClass *cl = fCaloClusters->GetClass();
958 if (!cl->GetBaseClass("AliVCluster")) {
959 AliFatal(Form("%s: Collection %s does not contain AliVCluster objects!", GetName(), fCaloClusters->GetName()));
960 fCaloClusters = 0;
961 return;
962 }
963 }
2f7259cf 964}