]> 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
Added SPD outlier trigger bit
[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 {
0a918d8d 226 offtrigger = ((AliVAODHeader*)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
46b40b53 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);
46b40b53 259
507f74bc 260}
261
262//________________________________________________________________________
263void AliAnalysisTaskEMCALClusterizeFast::CopyClusters(TClonesArray *orig, TClonesArray *dest)
264{
265 const Int_t Ncls = orig->GetEntries();
266
267 for(Int_t i=0; i < Ncls; ++i) {
268 AliVCluster *oc = static_cast<AliVCluster*>(orig->At(i));
269
270 if (!oc)
271 continue;
272
273 if (!oc->IsEMCAL())
274 continue;
275
276 AliVCluster *dc = static_cast<AliVCluster*>(dest->New(i));
277 dc->SetType(AliVCluster::kEMCALClusterv1);
278 dc->SetE(oc->E());
279 Float_t pos[3] = {0};
280 oc->GetPosition(pos);
281 dc->SetPosition(pos);
282 dc->SetNCells(oc->GetNCells());
283 dc->SetCellsAbsId(oc->GetCellsAbsId());
284 dc->SetCellsAmplitudeFraction(oc->GetCellsAmplitudeFraction());
285 dc->SetID(oc->GetID());
286 dc->SetDispersion(oc->GetDispersion());
287 dc->SetEmcCpvDistance(-1);
288 dc->SetChi2(-1);
289 dc->SetTOF(oc->GetTOF()); //time-of-flight
290 dc->SetNExMax(oc->GetNExMax()); //number of local maxima
291 dc->SetM02(oc->GetM02());
292 dc->SetM20(oc->GetM20());
293 dc->SetDistanceToBadChannel(oc->GetDistanceToBadChannel());
7eff06f2 294 dc->SetMCEnergyFraction(oc->GetMCEnergyFraction());
507f74bc 295
296 //MC
297 UInt_t nlabels = oc->GetNLabels();
298 Int_t *labels = oc->GetLabels();
299
300 if (nlabels == 0 || !labels)
301 continue;
302
303 AliESDCaloCluster *esdClus = dynamic_cast<AliESDCaloCluster*>(dc);
304 if (esdClus) {
305 TArrayI parents(nlabels, labels);
306 esdClus->AddLabels(parents);
307 }
308 else {
309 AliAODCaloCluster *aodClus = dynamic_cast<AliAODCaloCluster*>(dc);
310 if (aodClus)
311 aodClus->SetLabel(labels, nlabels);
312 }
313 }
2431b20f 314}
315
316//________________________________________________________________________
95f69c66 317void AliAnalysisTaskEMCALClusterizeFast::Clusterize()
2431b20f 318{
95f69c66 319 // Clusterize
7099ec08 320
95f69c66 321 if (fSubBackground) {
322 fClusterizer->SetInputCalibrated(kTRUE);
323 fClusterizer->SetCalibrationParameters(0);
2f7259cf 324 }
d3b5a5ec 325
95f69c66 326 fClusterizer->Digits2Clusters("");
fde82e42 327
95f69c66 328 if (fSubBackground) {
4889ed8b 329 if (fCalibData) {
330 fClusterizer->SetInputCalibrated(kFALSE);
331 fClusterizer->SetCalibrationParameters(fCalibData);
332 }
4cf9204b 333 }
4cf9204b 334}
335
2f7259cf 336//________________________________________________________________________
95f69c66 337void AliAnalysisTaskEMCALClusterizeFast::FillDigitsArray()
2f7259cf 338{
6265e793 339 // Fill digits array
55a8e3e4 340
3119d479 341 fDigitsArr->Clear("C");
3119d479 342 switch (fInputCellType) {
343
344 case kFEEData :
507f74bc 345 case kFEEDataMCOnly :
346 case kFEEDataExcludeMC :
3119d479 347 {
3119d479 348 Double_t avgE = 0; // for background subtraction
507f74bc 349 const Int_t ncells = fCaloCells->GetNumberOfCells();
3119d479 350 for (Int_t icell = 0, idigit = 0; icell < ncells; ++icell) {
351 Double_t cellAmplitude=0, cellTime=0, cellEFrac = 0;
60d77596 352 Short_t cellNumber=0;
353 Int_t cellMCLabel=-1;
507f74bc 354 if (fCaloCells->GetCell(icell, cellNumber, cellAmplitude, cellTime, cellMCLabel, cellEFrac) != kTRUE)
355 break;
5e1e18d1 356
357 if (cellMCLabel > 0 && cellEFrac < 1e-6)
358 cellEFrac = 1;
359
360 if (cellAmplitude < 1e-6 || cellNumber < 0)
361 continue;
362
fde82e42 363 if (fInputCellType == kFEEDataMCOnly) {
364 if (cellMCLabel <= 0)
365 continue;
366 else {
367 cellAmplitude *= cellEFrac;
368 cellEFrac = 1;
369 }
370 }
371 else if (fInputCellType == kFEEDataExcludeMC) {
372 if (cellMCLabel > 0)
373 continue;
5e1e18d1 374 else {
fde82e42 375 cellAmplitude *= 1 - cellEFrac;
5e1e18d1 376 cellEFrac = 0;
377 }
fde82e42 378 }
eff32c61 379
380 if(!AcceptCell(cellNumber)) continue;
381
507f74bc 382 AliEMCALDigit *digit = new((*fDigitsArr)[idigit]) AliEMCALDigit(cellMCLabel, cellMCLabel, cellNumber,
383 (Float_t)cellAmplitude, (Float_t)cellTime,
384 AliEMCALDigit::kHG,idigit, 0, 0, cellEFrac*cellAmplitude);
385
386 if (!fDoClusterize||fSubBackground) {
3119d479 387 Float_t energy = cellAmplitude;
388 Float_t time = cellTime;
389 fClusterizer->Calibrate(energy,time,cellNumber);
390 digit->SetAmplitude(energy);
507f74bc 391 avgE += energy;
fde82e42 392 }
3119d479 393 idigit++;
95f69c66 394 }
fde82e42 395
3119d479 396 if (fSubBackground) {
507f74bc 397 avgE /= fGeom->GetNumberOfSuperModules()*48*24;
3119d479 398 Int_t ndigis = fDigitsArr->GetEntries();
399 for (Int_t i = 0; i < ndigis; ++i) {
507f74bc 400 AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(i));
401 Double_t energy = digit->GetAmplitude() - avgE;
402 if (energy<=0.001) {
403 digit->SetAmplitude(0);
404 } else {
405 digit->SetAmplitude(energy);
406 }
3119d479 407 }
408 }
409 }
410 break;
411
412 case kPattern :
413 {
414 // Fill digits from a pattern
507f74bc 415 Int_t maxd = fGeom->GetNCells() / 4;
3119d479 416 for (Int_t idigit = 0; idigit < maxd; idigit++){
417 if (idigit % 24 == 12) idigit += 12;
418 AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->New(idigit));
419 digit->SetId(idigit * 4);
420 digit->SetTime(600);
421 digit->SetTimeR(600);
422 digit->SetIndexInList(idigit);
423 digit->SetType(AliEMCALDigit::kHG);
424 digit->SetAmplitude(0.1);
425 }
426 }
427 break;
428
429 case kL0FastORs :
430 case kL0FastORsTC :
431 case kL1FastORs :
432 {
433 // Fill digits from FastORs
434
435 AliVCaloTrigger *triggers = InputEvent()->GetCaloTrigger("EMCAL");
436
437 if (!triggers || !(triggers->GetEntries() > 0))
438 return;
439
440 Int_t idigit = 0;
441 triggers->Reset();
442
443 while ((triggers->Next())) {
444 Float_t L0Amplitude = 0;
445 triggers->GetAmplitude(L0Amplitude);
446
447 if (L0Amplitude <= 0 && fInputCellType != kL1FastORs)
448 continue;
449
450 Int_t L1Amplitude = 0;
451 triggers->GetL1TimeSum(L1Amplitude);
452
453 if (L1Amplitude <= 0 && fInputCellType == kL1FastORs)
454 continue;
455
456 Int_t triggerTime = 0;
457 Int_t ntimes = 0;
458 triggers->GetNL0Times(ntimes);
459
460 if (ntimes < 1 && fInputCellType == kL0FastORsTC)
461 continue;
462
463 if (ntimes > 0) {
464 Int_t trgtimes[25];
465 triggers->GetL0Times(trgtimes);
466 triggerTime = trgtimes[0];
467 }
468
469 Int_t triggerCol = 0, triggerRow = 0;
470 triggers->GetPosition(triggerCol, triggerRow);
471
472 Int_t find = -1;
507f74bc 473 fGeom->GetAbsFastORIndexFromPositionInEMCAL(triggerCol, triggerRow, find);
3119d479 474
475 if (find < 0)
476 continue;
477
478 Int_t cidx[4] = {-1};
507f74bc 479 Bool_t ret = fGeom->GetCellIndexFromFastORIndex(find, cidx);
3119d479 480
481 if (!ret)
482 continue;
483
484 Float_t triggerAmplitude = 0;
485
486 if (fInputCellType == kL1FastORs) {
487 triggerAmplitude = 0.25 * L1Amplitude; // it will add 4 cells for 1 amplitude
488 }
489 else {
490 triggerAmplitude = L0Amplitude; // 10 bit truncated, so it is already divided by 4
491 }
492
493 for (Int_t idxpos = 0; idxpos < 4; idxpos++) {
494 Int_t triggerNumber = cidx[idxpos];
495 AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->New(idigit));
496 digit->SetId(triggerNumber);
497 digit->SetTime(triggerTime);
498 digit->SetTimeR(triggerTime);
499 digit->SetIndexInList(idigit);
500 digit->SetType(AliEMCALDigit::kHG);
501 digit->SetAmplitude(triggerAmplitude);
502 idigit++;
503 }
7a7d5a64 504 }
505 }
3119d479 506 break;
06ca320f 507 }
2f7259cf 508}
509
eff32c61 510//________________________________________________________________________________________
511Bool_t AliAnalysisTaskEMCALClusterizeFast::AcceptCell(Int_t cellNumber) {
512
513 Bool_t accept = kTRUE;
514 if(fRejectExoticCells) {
515 //Remove exotic cells before making digits
a3e564aa 516 Bool_t exRemoval = fRecoUtils->IsRejectExoticCell();
eff32c61 517 fRecoUtils->SwitchOnRejectExoticCell();//switch on and off
518 Int_t bunchCrossNo = InputEvent()->GetBunchCrossNumber();
519 Bool_t isEx = fRecoUtils->IsExoticCell(cellNumber, fCaloCells, bunchCrossNo);
e4d5da1e 520 if(isEx) accept = kFALSE;
a3e564aa 521 if(!exRemoval) fRecoUtils->SwitchOffRejectExoticCell();//switch on and off
eff32c61 522 }
523 return accept;
524}
525
2f7259cf 526//________________________________________________________________________________________
507f74bc 527void AliAnalysisTaskEMCALClusterizeFast::CalibrateClusters()
4cf9204b 528{
a9edb9a3 529 // Go through clusters one by one and process separate correction
507f74bc 530 // as those were defined or not
a9edb9a3 531
507f74bc 532 Int_t nclusters = fCaloClusters->GetEntriesFast();
533 for (Int_t icluster=0; icluster < nclusters; ++icluster) {
534 AliVCluster *clust = static_cast<AliVCluster*>(fCaloClusters->At(icluster));
535 if (!clust)
536 continue;
537 if (!clust->IsEMCAL())
538 continue;
4cf9204b 539
507f74bc 540 // REMOVE CLUSTERS WITH BAD CELLS -----------------------------
541 if (fClusterBadChannelCheck) {
542 // careful, the the ClusterContainsBadChannel is dependent on
543 // SwitchOnBadChannelsRemoval, switching it ON automatically
544 // and returning to original value after processing
545 Bool_t badRemoval = fRecoUtils->IsBadChannelsRemovalSwitchedOn();
546 fRecoUtils->SwitchOnBadChannelsRemoval();
547
548 Bool_t badResult = fRecoUtils->ClusterContainsBadChannel(fGeom, clust->GetCellsAbsId(), clust->GetNCells());
549
550 // switch the bad channels removal back
551 if (!badRemoval)
552 fRecoUtils->SwitchOffBadChannelsRemoval();
553
554 if (badResult) {
555 delete fCaloClusters->RemoveAt(icluster);
556 continue; //TODO is it really needed to remove it? Or should we flag it?
557 }
558 }
559
560 // REMOVE EXOTIC CLUSTERS -------------------------------------
561 // does process local cell recalibration energy and time without replacing
562 // the global cell values, in case of no cell recalib done yet
563 if (fRejectExoticClusters) {
564 // careful, the IsExoticCluster is dependent on
565 // SwitchOnRejectExoticCell, switching it ON automatically
566 // and returning to original value after processing
567 Bool_t exRemoval = fRecoUtils->IsRejectExoticCell();
568 fRecoUtils->SwitchOnRejectExoticCell();
569
570 // get bunch crossing
571 Int_t bunchCrossNo = InputEvent()->GetBunchCrossNumber();
572
573 Bool_t exResult = fRecoUtils->IsExoticCluster(clust, fCaloCells, bunchCrossNo);
574
575 // switch the exotic channels removal back
576 if (!exRemoval)
577 fRecoUtils->SwitchOffRejectExoticCell();
578
579 if (exResult) {
580 delete fCaloClusters->RemoveAt(icluster);
581 continue; //TODO is it really needed to remove it? Or should we flag it?
582 }
583 }
584
585 // FIDUCIAL CUT -----------------------------------------------
586 if (fFiducial) {
587 // depends on SetNumberOfCellsFromEMCALBorder
588 // SwitchOnNoFiducialBorderInEMCALEta0
589 if (!fRecoUtils->CheckCellFiducialRegion(fGeom, clust, fCaloCells)){
590 delete fCaloClusters->RemoveAt(icluster);
591 continue; //TODO it would be nice to store the distance
592 }
593 }
594
595 // NONLINEARITY -----------------------------------------------
596 if (fDoNonLinearity) {
597 Float_t correctedEnergy = fRecoUtils->CorrectClusterEnergyLinearity(clust);
598 clust->SetE(correctedEnergy);
599 }
0408f54f 600
507f74bc 601 // DISTANCE TO BAD CHANNELS -----------------------------------
602 if (fRecalDistToBadChannels)
603 fRecoUtils->RecalculateClusterDistanceToBadChannel(fGeom, fCaloCells, clust);
604 }
605
606 fCaloClusters->Compress();
607}
608
609//________________________________________________________________________________________
610void AliAnalysisTaskEMCALClusterizeFast::TrackClusterMatching(AliVCluster *c, TClonesArray *tarr)
611{
612 Float_t g[3]={0};
613 c->GetPosition(g);
614 TVector3 gpos(g);
615
616 Double_t dEtaMin = 1e9;
617 Double_t dPhiMin = 1e9;
618 Int_t imin = -1;
619 Double_t ceta = gpos.Eta();
620 Double_t cphi = gpos.Phi();
621 const Int_t ntrks = tarr->GetEntries();
622 for(Int_t t = 0; t<ntrks; ++t) {
623 AliVTrack *track = static_cast<AliVTrack*>(tarr->At(t));
624 if (!track)
625 continue;
626 const AliExternalTrackParam *outp = track->GetOuterParam();
627 if (!outp)
628 continue;
629 Double_t trkPos[3] = {0.,0.,0.};
630 if (!outp->GetXYZ(trkPos))
631 continue;
632 TVector3 vec(trkPos[0],trkPos[1],trkPos[2]);
633 Double_t veta = vec.Eta();
634 Double_t vphi = vec.Phi();
635 if(vphi<0)
636 vphi += 2*TMath::Pi();
637 if (TMath::Abs(veta)>0.75 || (vphi<70*TMath::DegToRad()) || (vphi>190*TMath::DegToRad()))
638 continue;
639 Double_t dR = vec.DeltaR(gpos);
640 if(dR > 25)
641 continue;
642 Float_t tmpEta=0, tmpPhi=0;
643 if (0) {
644 AliExternalTrackParam trkParTemp(*outp); // retrieve the starting point every time before the extrapolation
645 Bool_t ret = fRecoUtils->ExtrapolateTrackToCluster(&trkParTemp, c, fRecoUtils->GetMass(), fRecoUtils->GetStep(), tmpEta, tmpPhi);
646 if (!ret)
647 continue;
648 } else {
649 tmpEta = ceta - veta;
650 tmpPhi = cphi - vphi;
651 }
652 if (TMath::Abs(tmpEta)<TMath::Abs(dEtaMin) && TMath::Abs(tmpPhi)<TMath::Abs(dPhiMin)) {
653 dEtaMin = tmpEta;
654 dPhiMin = tmpPhi;
655 imin = t;
656 }
657 }
658 c->SetEmcCpvDistance(imin);
659 c->SetTrackDistance(dPhiMin, dEtaMin);
660}
661
662//________________________________________________________________________________________
663void AliAnalysisTaskEMCALClusterizeFast::RecPoints2Clusters(TClonesArray *clus)
664{
665 // Cluster energy, global position, cells and their amplitude fractions are restored.
666
0408f54f 667 // tracks array for track/cluster matching
668 TClonesArray *tarr = 0;
669 if (!fTrackName.IsNull()) {
507f74bc 670 tarr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTrackName));
0408f54f 671 if (!tarr) {
672 AliError(Form("Cannot get tracks named %s", fTrackName.Data()));
673 }
674 }
fde82e42 675
eb33b7f8 676 const Int_t Ncls = fClusterArr->GetEntries();
0408f54f 677 AliDebug(1, Form("total no of clusters %d", Ncls));
4cf9204b 678 for(Int_t i=0, nout=clus->GetEntries(); i < Ncls; ++i) {
679 AliEMCALRecPoint *recpoint = static_cast<AliEMCALRecPoint*>(fClusterArr->At(i));
680 Int_t ncells_true = 0;
681 const Int_t ncells = recpoint->GetMultiplicity();
682 UShort_t absIds[ncells];
683 Double32_t ratios[ncells];
684 Int_t *dlist = recpoint->GetDigitsList();
685 Float_t *elist = recpoint->GetEnergiesList();
7eff06f2 686 Double_t mcEnergy = 0;
4cf9204b 687 for (Int_t c = 0; c < ncells; ++c) {
688 AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(dlist[c]));
689 absIds[ncells_true] = digit->GetId();
5ded2286 690 ratios[ncells_true] = elist[c]/recpoint->GetEnergy();
7eff06f2 691 if (digit->GetIparent(1) > 0)
692 mcEnergy += digit->GetDEParent(1)/recpoint->GetEnergy();
84ccad86 693 ++ncells_true;
4cf9204b 694 }
695
696 if (ncells_true < 1) {
697 AliWarning("Skipping cluster with no cells");
698 continue;
699 }
700
701 // calculate new cluster position
702 TVector3 gpos;
4cf9204b 703 recpoint->GetGlobalPosition(gpos);
84ccad86 704 Float_t g[3];
4cf9204b 705 gpos.GetXYZ(g);
706
e93ec1f7 707 AliDebug(1, Form("energy %f", recpoint->GetEnergy()));
708
95f69c66 709 AliVCluster *c = static_cast<AliVCluster*>(clus->New(nout++));
4cf9204b 710 c->SetType(AliVCluster::kEMCALClusterv1);
711 c->SetE(recpoint->GetEnergy());
712 c->SetPosition(g);
713 c->SetNCells(ncells_true);
a55e4f1d 714 c->SetCellsAbsId(absIds);
715 c->SetCellsAmplitudeFraction(ratios);
716 c->SetID(recpoint->GetUniqueID());
4cf9204b 717 c->SetDispersion(recpoint->GetDispersion());
0408f54f 718 c->SetEmcCpvDistance(-1);
719 c->SetChi2(-1);
4cf9204b 720 c->SetTOF(recpoint->GetTime()) ; //time-of-flight
721 c->SetNExMax(recpoint->GetNExMax()); //number of local maxima
722 Float_t elipAxis[2];
723 recpoint->GetElipsAxis(elipAxis);
e93ec1f7 724 c->SetM02(elipAxis[0]*elipAxis[0]);
725 c->SetM20(elipAxis[1]*elipAxis[1]);
7eff06f2 726 c->SetMCEnergyFraction(mcEnergy);
507f74bc 727
728 //MC
729 AliESDCaloCluster *esdClus = dynamic_cast<AliESDCaloCluster*>(c);
730 if (esdClus) {
731 Int_t parentMult = 0;
732 Int_t *parentList = recpoint->GetParents(parentMult);
733 if (parentMult > 0) {
734 TArrayI parents(parentMult, parentList);
735 esdClus->AddLabels(parents);
736 }
737 }
738 else {
739 AliAODCaloCluster *aodClus = dynamic_cast<AliAODCaloCluster*>(c);
740 if (aodClus) {
741 Int_t parentMult = 0;
742 Int_t *parentList = recpoint->GetParents(parentMult);
743 aodClus->SetLabel(parentList, parentMult);
95f69c66 744 }
2f7259cf 745 }
507f74bc 746
747 if (tarr)
748 TrackClusterMatching(c, tarr);
2f7259cf 749 }
750}
751
95f69c66 752//________________________________________________________________________
753void AliAnalysisTaskEMCALClusterizeFast::UpdateCells()
754{
755 // Update cells in case re-calibration was done.
95f69c66 756 if (!fCalibData&&!fSubBackground)
757 return;
758
507f74bc 759 const Int_t ncells = fCaloCells->GetNumberOfCells();
55a8e3e4 760 const Int_t ndigis = fDigitsArr->GetEntries();
95f69c66 761 if (ncells!=ndigis) {
507f74bc 762 fCaloCells->DeleteContainer();
763 fCaloCells->CreateContainer(ndigis);
95f69c66 764 }
765 for (Int_t idigit = 0; idigit < ndigis; ++idigit) {
766 AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(idigit));
767 Double_t cellAmplitude = digit->GetCalibAmp();
55a8e3e4 768 Short_t cellNumber = digit->GetId();
769 Double_t cellTime = digit->GetTime();
507f74bc 770 fCaloCells->SetCell(idigit, cellNumber, cellAmplitude, cellTime);
95f69c66 771 }
772}
773
774//________________________________________________________________________
775void AliAnalysisTaskEMCALClusterizeFast::UpdateClusters()
776{
777 // Update cells in case re-calibration was done.
7a7d5a64 778
507f74bc 779 const Int_t nents = fCaloClusters->GetEntries();
780 for (Int_t i=0;i<nents;++i) {
781 AliVCluster *c = static_cast<AliVCluster*>(fCaloClusters->At(i));
782 if (!c)
783 continue;
784 if (c->IsEMCAL())
785 delete fCaloClusters->RemoveAt(i);
7a7d5a64 786 }
507f74bc 787
788 fCaloClusters->Compress();
7a7d5a64 789
507f74bc 790 RecPoints2Clusters(fCaloClusters);
95f69c66 791}
792
2f7259cf 793//________________________________________________________________________________________
794void AliAnalysisTaskEMCALClusterizeFast::Init()
795{
4f66d2ba 796 // Select clusterization/unfolding algorithm and set all the needed parameters.
01bc1ce8 797
507f74bc 798 if (InputEvent()->GetRunNumber()==fRun)
2f7259cf 799 return;
507f74bc 800 fRun = InputEvent()->GetRunNumber();
2f7259cf 801
802 if (fJustUnfold){
803 // init the unfolding afterburner
804 delete fUnfolder;
1759f743 805 fUnfolder = new AliEMCALAfterBurnerUF(fRecParam->GetW0(),fRecParam->GetLocMaxCut(),fRecParam->GetMinECut());
2f7259cf 806 return;
807 }
808
4f66d2ba 809 if (fGeomName.Length()>0)
507f74bc 810 fGeom = AliEMCALGeometry::GetInstance(fGeomName);
4f66d2ba 811 else
507f74bc 812 fGeom = AliEMCALGeometry::GetInstance();
813 if (!fGeom) {
2f7259cf 814 AliFatal("Geometry not available!!!");
815 return;
816 }
817
818 if (!fGeomMatrixSet) {
819 if (fLoadGeomMatrices) {
507f74bc 820 for(Int_t mod=0; mod < fGeom->GetNumberOfSuperModules(); ++mod) {
06ca320f 821 if (fGeomMatrix[mod]){
822 if (DebugLevel() > 2)
2f7259cf 823 fGeomMatrix[mod]->Print();
507f74bc 824 fGeom->SetMisalMatrix(fGeomMatrix[mod],mod);
2f7259cf 825 }
826 }
29b496d5 827 } else { // get matrix from file (work around bug in aliroot)
507f74bc 828 for(Int_t mod=0; mod < fGeom->GetEMCGeometry()->GetNumberOfSuperModules(); ++mod) {
29b496d5 829 const TGeoHMatrix *gm = 0;
507f74bc 830 if (fEsd) {
831 gm = fEsd->GetEMCALMatrix(mod);
29b496d5 832 } else {
0a918d8d 833 AliAODHeader *aodheader = dynamic_cast<AliAODHeader*>(fAod->GetHeader());
834 if(!aodheader) AliFatal("Not a standard AOD");
29b496d5 835 if (aodheader) {
836 gm = aodheader->GetEMCALMatrix(mod);
837 }
838 }
839 if (gm) {
06ca320f 840 if (DebugLevel() > 2)
29b496d5 841 gm->Print();
507f74bc 842 fGeom->SetMisalMatrix(gm,mod);
2f7259cf 843 }
844 }
845 }
846 fGeomMatrixSet=kTRUE;
847 }
848
849 // setup digit array if needed
850 if (!fDigitsArr) {
851 fDigitsArr = new TClonesArray("AliEMCALDigit", 1000);
852 fDigitsArr->SetOwner(1);
853 }
854
855 // then setup clusterizer
cac0e10b 856 if (fClusterizer) {
857 // avoid to delete digits array
858 fClusterizer->SetDigitsArr(0);
859 delete fClusterizer;
860 }
507f74bc 861 if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv1)
862 fClusterizer = new AliEMCALClusterizerv1(fGeom);
06ca320f 863 else if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerNxN) {
80dbd4f7 864 AliEMCALClusterizerNxN *clusterizer = new AliEMCALClusterizerNxN(fGeom);
865 clusterizer->SetNRowDiff(fRecParam->GetNRowDiff()); //MV: already done in AliEMCALClusterizer::InitParameters
866 clusterizer->SetNColDiff(fRecParam->GetNColDiff()); //MV: already done in AliEMCALClusterizer::InitParameters
3c56da8f 867 fClusterizer = clusterizer;
507f74bc 868 }
869 else if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv2)
870 fClusterizer = new AliEMCALClusterizerv2(fGeom);
871 else if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerFW) {
872 AliEMCALClusterizerFixedWindow *clusterizer = new AliEMCALClusterizerFixedWindow(fGeom);
0fda97fe 873 clusterizer->SetNphi(fNPhi);
874 clusterizer->SetNeta(fNEta);
875 clusterizer->SetShiftPhi(fShiftPhi);
876 clusterizer->SetShiftEta(fShiftEta);
3c56da8f 877 clusterizer->SetTRUshift(fTRUShift);
7a7d5a64 878 fClusterizer = clusterizer;
879 }
507f74bc 880 else {
2f7259cf 881 AliFatal(Form("Clusterizer < %d > not available", fRecParam->GetClusterizerFlag()));
882 }
883 fClusterizer->InitParameters(fRecParam);
01bc1ce8 884
885 if ((!fCalibData&&fLoadCalib) || (!fPedestalData&&fLoadPed)) {
886 AliCDBManager *cdb = AliCDBManager::Instance();
887 if (!cdb->IsDefaultStorageSet() && !fOCDBpath.IsNull())
888 cdb->SetDefaultStorage(fOCDBpath);
889 if (fRun!=cdb->GetRun())
890 cdb->SetRun(fRun);
891 }
349bf69e 892 if (!fCalibData&&fLoadCalib&&fRun>0) {
2f7259cf 893 AliCDBEntry *entry = static_cast<AliCDBEntry*>(AliCDBManager::Instance()->Get("EMCAL/Calib/Data"));
894 if (entry)
895 fCalibData = static_cast<AliEMCALCalibData*>(entry->GetObject());
896 if (!fCalibData)
897 AliFatal("Calibration parameters not found in CDB!");
898 }
349bf69e 899 if (!fPedestalData&&fLoadPed&&fRun>0) {
2f7259cf 900 AliCDBEntry *entry = static_cast<AliCDBEntry*>(AliCDBManager::Instance()->Get("EMCAL/Calib/Pedestals"));
901 if (entry)
902 fPedestalData = static_cast<AliCaloCalibPedestal*>(entry->GetObject());
903 }
10d33986 904 if (fCalibData) {
905 fClusterizer->SetInputCalibrated(kFALSE);
906 fClusterizer->SetCalibrationParameters(fCalibData);
10d33986 907 } else {
908 fClusterizer->SetInputCalibrated(kTRUE);
909 }
e93ec1f7 910 fClusterizer->SetCaloCalibPedestal(fPedestalData);
10d33986 911 fClusterizer->SetJustClusters(kTRUE);
2f7259cf 912 fClusterizer->SetDigitsArr(fDigitsArr);
913 fClusterizer->SetOutput(0);
914 fClusterArr = const_cast<TObjArray *>(fClusterizer->GetRecPoints());
507f74bc 915
916 // Get the emcal cells
fde82e42 917 if ((fInputCellType == kFEEData || fInputCellType == kFEEDataMCOnly || fInputCellType == kFEEDataExcludeMC) && !fCaloCells) {
507f74bc 918 if (fCaloCellsName.IsNull()) {
919 fCaloCells = InputEvent()->GetEMCALCells();
920 }
921 else {
922 fCaloCells = dynamic_cast<AliVCaloCells*>(InputEvent()->FindListObject(fCaloCellsName));
923 if (!fCaloCells)
924 AliError(Form("%s: Could not retrieve cells %s!", GetName(), fCaloCellsName.Data()));
925 }
926 if (!fCaloCells)
927 AliFatal("Could not get EMCal cells!");
928 }
929
930 // Set output clusters collection
931 if (!fAttachClusters) {
932 fCaloClusters = fOutputAODBranch;
933 return;
934 }
935
936 if (!fCaloClusters) {
937 if (fCaloClustersName.IsNull()) { //overwrite mode
938 if (fEsd)
939 fCaloClusters = static_cast<TClonesArray*>(InputEvent()->FindListObject("CaloClusters"));
940 else if (fAod)
941 fCaloClusters = static_cast<TClonesArray*>(InputEvent()->FindListObject("caloClusters"));
942 }
943 else {
944 fCaloClusters = static_cast<TClonesArray*>(InputEvent()->FindListObject(fCaloClustersName));
945
946 if (!fCaloClusters) {
947 if (fEsd)
948 fCaloClusters = new TClonesArray("AliESDCaloCluster");
949 else if (fAod)
950 fCaloClusters = new TClonesArray("AliAODCaloCluster");
951
952 fCaloClusters->SetName(fCaloClustersName);
953 InputEvent()->AddObject(fCaloClusters);
954 }
955 }
956
957 if (!fCaloClusters)
958 AliFatal("Could not get cluster collection!");
959
960 TClass *cl = fCaloClusters->GetClass();
961 if (!cl->GetBaseClass("AliVCluster")) {
962 AliFatal(Form("%s: Collection %s does not contain AliVCluster objects!", GetName(), fCaloClusters->GetName()));
963 fCaloClusters = 0;
964 return;
965 }
966 }
2f7259cf 967}