1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 ///////////////////////////////////////////////////////////////////////////////
18 // EMCAL tender, apply corrections to EMCAL clusters and do track matching. //
20 // Author: Deepa Thomas (Utrecht University) //
21 // Later mods/rewrite: Jiri Kral (University of Jyvaskyla) //
22 // S. Aiola, C. Loizides : Make it work for AODs //
24 ///////////////////////////////////////////////////////////////////////////////
26 #include <TClonesArray.h>
28 #include <TGeoGlobalMagField.h>
29 #include <TInterpreter.h>
30 #include <TObjArray.h>
33 #include "AliAODEvent.h"
34 #include "AliAODMCParticle.h"
35 #include "AliAnalysisManager.h"
36 #include "AliEMCALAfterBurnerUF.h"
37 #include "AliEMCALClusterizer.h"
38 #include "AliEMCALClusterizerNxN.h"
39 #include "AliEMCALClusterizerv1.h"
40 #include "AliEMCALClusterizerv2.h"
41 #include "AliEMCALDigit.h"
42 #include "AliEMCALGeometry.h"
43 #include "AliEMCALRecParam.h"
44 #include "AliEMCALRecParam.h"
45 #include "AliEMCALRecPoint.h"
46 #include "AliEMCALRecoUtils.h"
47 #include "AliESDCaloCluster.h"
48 #include "AliESDEvent.h"
51 #include "AliOADBContainer.h"
52 #include "AliTender.h"
53 #include "AliEMCALTenderSupply.h"
55 ClassImp(AliEMCALTenderSupply)
57 AliEMCALTenderSupply::AliEMCALTenderSupply() :
67 ,fNonLinearThreshold(-1)
68 ,fReCalibCluster(kFALSE)
70 ,fCalibrateEnergy(kFALSE)
71 ,fCalibrateTime(kFALSE)
72 ,fCalibrateTimeParamAvailable(kFALSE)
73 ,fDoNonLinearity(kFALSE)
74 ,fBadCellRemove(kFALSE)
75 ,fRejectExoticCells(kFALSE)
76 ,fRejectExoticClusters(kFALSE)
77 ,fClusterBadChannelCheck(kFALSE)
78 ,fRecalClusPos(kFALSE)
80 ,fNCellsFromEMCALBorder(-1)
81 ,fRecalDistToBadChannels(kFALSE)
82 ,fRecalShowerShape(kFALSE)
85 ,fGetPassFromFileName(kTRUE)
90 ,fCutEtaPhiSeparate(kFALSE)
95 ,fReClusterize(kFALSE)
97 ,fGeomMatrixSet(kFALSE)
98 ,fLoadGeomMatrices(kFALSE)
100 ,fDoTrackMatch(kFALSE)
101 ,fDoUpdateOnly(kFALSE)
105 ,fMisalignSurvey(kdefault)
106 ,fExoticCellFraction(-1)
107 ,fExoticCellDiffTime(-1)
108 ,fExoticCellMinAmplitude(-1)
109 ,fSetCellMCLabelFromCluster(0)
111 ,fRemapMCLabelForAODs(0)
113 // Default constructor.
115 for(Int_t i = 0; i < 12; i++) fEMCALMatrix[i] = 0 ;
116 for(Int_t j = 0; j < 12672; j++) fOrgClusterCellId[j] = -1;
119 //_____________________________________________________
120 AliEMCALTenderSupply::AliEMCALTenderSupply(const char *name, const AliTender *tender) :
121 AliTenderSupply(name,tender)
130 ,fNonLinearThreshold(-1)
131 ,fReCalibCluster(kFALSE)
133 ,fCalibrateEnergy(kFALSE)
134 ,fCalibrateTime(kFALSE)
135 ,fCalibrateTimeParamAvailable(kFALSE)
136 ,fDoNonLinearity(kFALSE)
137 ,fBadCellRemove(kFALSE)
138 ,fRejectExoticCells(kFALSE)
139 ,fRejectExoticClusters(kFALSE)
140 ,fClusterBadChannelCheck(kFALSE)
141 ,fRecalClusPos(kFALSE)
143 ,fNCellsFromEMCALBorder(-1)
144 ,fRecalDistToBadChannels(kFALSE)
145 ,fRecalShowerShape(kFALSE)
148 ,fGetPassFromFileName(kTRUE)
152 ,fCutEtaPhiSum(kTRUE)
153 ,fCutEtaPhiSeparate(kFALSE)
158 ,fReClusterize(kFALSE)
160 ,fGeomMatrixSet(kFALSE)
161 ,fLoadGeomMatrices(kFALSE)
163 ,fDoTrackMatch(kFALSE)
164 ,fDoUpdateOnly(kFALSE)
168 ,fMisalignSurvey(kdefault)
169 ,fExoticCellFraction(-1)
170 ,fExoticCellDiffTime(-1)
171 ,fExoticCellMinAmplitude(-1)
172 ,fSetCellMCLabelFromCluster(0)
174 ,fRemapMCLabelForAODs(0)
178 for(Int_t i = 0; i < 12; i++) fEMCALMatrix[i] = 0 ;
179 for(Int_t j = 0; j < 12672; j++) fOrgClusterCellId[j] = -1;
182 //_____________________________________________________
183 AliEMCALTenderSupply::AliEMCALTenderSupply(const char *name, AliAnalysisTaskSE *task) :
184 AliTenderSupply(name)
193 ,fNonLinearThreshold(-1)
194 ,fReCalibCluster(kFALSE)
196 ,fCalibrateEnergy(kFALSE)
197 ,fCalibrateTime(kFALSE)
198 ,fCalibrateTimeParamAvailable(kFALSE)
199 ,fDoNonLinearity(kFALSE)
200 ,fBadCellRemove(kFALSE)
201 ,fRejectExoticCells(kFALSE)
202 ,fRejectExoticClusters(kFALSE)
203 ,fClusterBadChannelCheck(kFALSE)
204 ,fRecalClusPos(kFALSE)
206 ,fNCellsFromEMCALBorder(-1)
207 ,fRecalDistToBadChannels(kFALSE)
208 ,fRecalShowerShape(kFALSE)
211 ,fGetPassFromFileName(kTRUE)
215 ,fCutEtaPhiSum(kTRUE)
216 ,fCutEtaPhiSeparate(kFALSE)
221 ,fReClusterize(kFALSE)
223 ,fGeomMatrixSet(kFALSE)
224 ,fLoadGeomMatrices(kFALSE)
226 ,fDoTrackMatch(kFALSE)
227 ,fDoUpdateOnly(kFALSE)
231 ,fMisalignSurvey(kdefault)
232 ,fExoticCellFraction(-1)
233 ,fExoticCellDiffTime(-1)
234 ,fExoticCellMinAmplitude(-1)
235 ,fSetCellMCLabelFromCluster(0)
237 ,fRemapMCLabelForAODs(0)
239 // Named constructor.
241 for(Int_t i = 0; i < 12; i++) fEMCALMatrix[i] = 0 ;
242 for(Int_t j = 0; j < 12672; j++) fOrgClusterCellId[j] = -1;
245 //_____________________________________________________
246 AliEMCALTenderSupply::~AliEMCALTenderSupply()
250 if (!AliAnalysisManager::GetAnalysisManager()) return;
252 if (!AliAnalysisManager::GetAnalysisManager()->IsProofMode())
254 delete fEMCALRecoUtils;
262 fDigitsArr->Clear("C");
274 //_____________________________________________________
275 void AliEMCALTenderSupply::SetDefaults()
277 // Set default settings.
279 SwitchOnReclustering();
280 SwitchOnTrackMatch();
283 //_____________________________________________________
284 Bool_t AliEMCALTenderSupply::RunChanged() const
288 return (fTender && fTender->RunChanged()) || (fTask && fRun != fTask->InputEvent()->GetRunNumber());
291 //_____________________________________________________
292 void AliEMCALTenderSupply::Init()
294 // Initialise EMCAL tender.
297 AliWarning("Init EMCAL Tender supply");
299 if (fConfigName.Length()>0 && gROOT->LoadMacro(fConfigName) >=0) {
300 AliDebug(1, Form("Loading settings from macro %s", fConfigName.Data()));
301 AliEMCALTenderSupply *tender = (AliEMCALTenderSupply*)gInterpreter->ProcessLine("ConfigEMCALTenderSupply()");
302 fDebugLevel = tender->fDebugLevel;
303 fEMCALGeoName = tender->fEMCALGeoName;
304 fEMCALRecoUtils = tender->fEMCALRecoUtils;
305 fConfigName = tender->fConfigName;
306 fNonLinearFunc = tender->fNonLinearFunc;
307 fNonLinearThreshold = tender->fNonLinearThreshold;
308 fReCalibCluster = tender->fReCalibCluster;
309 fUpdateCell = tender->fUpdateCell;
310 fRecalClusPos = tender->fRecalClusPos;
311 fCalibrateEnergy = tender->fCalibrateEnergy;
312 fCalibrateTime = tender->fCalibrateTime;
313 fCalibrateTimeParamAvailable = tender->fCalibrateTimeParamAvailable;
314 fFiducial = tender->fFiducial;
315 fNCellsFromEMCALBorder = tender->fNCellsFromEMCALBorder;
316 fRecalDistToBadChannels = tender->fRecalDistToBadChannels;
317 fRecalShowerShape = tender->fRecalShowerShape;
318 fClusterBadChannelCheck = tender->fClusterBadChannelCheck;
319 fBadCellRemove = tender->fBadCellRemove;
320 fRejectExoticCells = tender->fRejectExoticCells;
321 fRejectExoticClusters = tender->fRejectExoticClusters;
322 fMass = tender->fMass;
323 fStep = tender->fStep;
324 fCutEtaPhiSum = tender->fCutEtaPhiSum;
325 fCutEtaPhiSeparate = tender->fCutEtaPhiSeparate;
326 fRcut = tender->fRcut;
327 fEtacut = tender->fEtacut;
328 fPhicut = tender->fPhicut;
329 fReClusterize = tender->fReClusterize;
330 fLoadGeomMatrices = tender->fLoadGeomMatrices;
331 fRecParam = tender->fRecParam;
332 fDoNonLinearity = tender->fDoNonLinearity;
333 fDoTrackMatch = tender->fDoTrackMatch;
334 fDoUpdateOnly = tender->fDoUpdateOnly;
335 fMisalignSurvey = tender->fMisalignSurvey;
336 fExoticCellFraction = tender->fExoticCellFraction;
337 fExoticCellDiffTime = tender->fExoticCellDiffTime;
338 fExoticCellMinAmplitude = tender->fExoticCellMinAmplitude;
340 for(Int_t i = 0; i < 12; i++)
341 fEMCALMatrix[i] = tender->fEMCALMatrix[i] ;
345 AliInfo("Emcal Tender settings: ======================================");
346 AliInfo("------------ Switches --------------------------");
347 AliInfo(Form("BadCellRemove : %d", fBadCellRemove));
348 AliInfo(Form("ExoticCellRemove : %d", fRejectExoticCells));
349 AliInfo(Form("CalibrateEnergy : %d", fCalibrateEnergy));
350 AliInfo(Form("CalibrateTime : %d", fCalibrateTime));
351 AliInfo(Form("UpdateCell : %d", fUpdateCell));
352 AliInfo(Form("DoUpdateOnly : %d", fDoUpdateOnly));
353 AliInfo(Form("Reclustering : %d", fReClusterize));
354 AliInfo(Form("ClusterBadChannelCheck : %d", fClusterBadChannelCheck));
355 AliInfo(Form("ClusterExoticChannelCheck : %d", fRejectExoticClusters));
356 AliInfo(Form("CellFiducialRegion : %d", fFiducial));
357 AliInfo(Form("ReCalibrateCluster : %d", fReCalibCluster));
358 AliInfo(Form("RecalculateClusPos : %d", fRecalClusPos));
359 AliInfo(Form("RecalShowerShape : %d", fRecalShowerShape));
360 AliInfo(Form("NonLinearityCorrection : %d", fDoNonLinearity));
361 AliInfo(Form("RecalDistBadChannel : %d", fRecalDistToBadChannels));
362 AliInfo(Form("TrackMatch : %d", fDoTrackMatch));
363 AliInfo("------------ Variables -------------------------");
364 AliInfo(Form("DebugLevel : %d", fDebugLevel));
365 AliInfo(Form("BasePath : %s", fBasePath.Data()));
366 AliInfo(Form("ConfigFileName : %s", fConfigName.Data()));
367 AliInfo(Form("EMCALGeometryName : %s", fEMCALGeoName.Data()));
368 AliInfo(Form("NonLinearityFunction : %d", fNonLinearFunc));
369 AliInfo(Form("NonLinearityThreshold : %d", fNonLinearThreshold));
370 AliInfo(Form("MisalignmentMatrixSurvey : %d", fMisalignSurvey));
371 AliInfo(Form("NumberOfCellsFromEMCALBorder : %d", fNCellsFromEMCALBorder));
372 AliInfo(Form("RCut : %f", fRcut));
373 AliInfo(Form("Mass : %f", fMass));
374 AliInfo(Form("Step : %f", fStep));
375 AliInfo(Form("EtaCut : %f", fEtacut));
376 AliInfo(Form("PhiCut : %f", fPhicut));
377 AliInfo(Form("ExoticCellFraction : %f", fExoticCellFraction));
378 AliInfo(Form("ExoticCellDiffTime : %f", fExoticCellDiffTime));
379 AliInfo(Form("ExoticCellMinAmplitude : %f", fExoticCellMinAmplitude));
380 AliInfo("=============================================================");
385 if (!fEMCALRecoUtils)
386 fEMCALRecoUtils = new AliEMCALRecoUtils;
388 // init geometry if requested
389 if (fEMCALGeoName.Length()>0)
390 fEMCALGeo = AliEMCALGeometry::GetInstance(fEMCALGeoName) ;
393 fDigitsArr = new TClonesArray("AliEMCALDigit",1000);
395 // initialising non-linearity parameters
396 if (fNonLinearThreshold != -1)
397 fEMCALRecoUtils->SetNonLinearityThreshold(fNonLinearThreshold);
398 if (fNonLinearFunc != -1)
399 fEMCALRecoUtils->SetNonLinearityFunction(fNonLinearFunc);
401 // missalignment function
402 fEMCALRecoUtils->SetPositionAlgorithm(AliEMCALRecoUtils::kPosTowerGlobal);
405 // do not do the eta0 fiducial cut
406 if (fNCellsFromEMCALBorder != -1)
407 fEMCALRecoUtils->SetNumberOfCellsFromEMCALBorder(fNCellsFromEMCALBorder);
408 fEMCALRecoUtils->SwitchOnNoFiducialBorderInEMCALEta0();
410 // exotic cell rejection
411 if (fExoticCellFraction != -1)
412 fEMCALRecoUtils->SetExoticCellFractionCut(fExoticCellFraction);
413 if (fExoticCellDiffTime != -1)
414 fEMCALRecoUtils->SetExoticCellDiffTimeCut(fExoticCellDiffTime);
415 if (fExoticCellMinAmplitude != -1)
416 fEMCALRecoUtils->SetExoticCellMinAmplitudeCut(fExoticCellMinAmplitude);
418 // setting track matching parameters ... mass, step size and residual cut
420 fEMCALRecoUtils->SetMass(fMass);
422 fEMCALRecoUtils->SetStep(fStep);
424 // spatial cut based on separate eta/phi or common processing
426 fEMCALRecoUtils->SwitchOnCutEtaPhiSum();
428 fEMCALRecoUtils->SetCutR(fRcut);
429 } else if (fCutEtaPhiSeparate) {
430 fEMCALRecoUtils->SwitchOnCutEtaPhiSeparate();
432 fEMCALRecoUtils->SetCutEta(fEtacut);
434 fEMCALRecoUtils->SetCutPhi(fPhicut);
438 //_____________________________________________________
439 AliVEvent* AliEMCALTenderSupply::GetEvent()
441 // Return the event pointer.
444 return fTender->GetEvent();
446 return fTask->InputEvent();
452 //_____________________________________________________
453 void AliEMCALTenderSupply::ProcessEvent()
457 AliVEvent *event = GetEvent();
460 AliError("Event ptr = 0, returning");
464 // Initialising parameters once per run number
467 fRun = event->GetRunNumber();
468 AliWarning(Form("Run changed, initializing parameters for %d", fRun));
469 if (dynamic_cast<AliAODEvent*>(event)) {
470 AliWarning("=============================================================");
471 AliWarning("=== Running on AOD is not equivalent to running on ESD! ===");
472 AliWarning("=============================================================");
475 // init geometry if not already done
476 if (fEMCALGeoName.Length()==0) {
477 fEMCALGeoName = "EMCAL_FIRSTYEARV1";
479 fEMCALGeoName = "EMCAL_COMPLETEV1";
482 fEMCALGeoName = "EMCAL_COMPLETE12SMV1";
484 fEMCALGeo = AliEMCALGeometry::GetInstance(fEMCALGeoName);
486 AliFatal(Form("Can not create geometry: %s", fEMCALGeoName.Data()));
492 if (fGetPassFromFileName)
495 // define what recalib parameters are needed for various switches
496 // this is based on implementation in AliEMCALRecoUtils
497 Bool_t needRecoParam = fReClusterize;
498 Bool_t needBadChannels = fBadCellRemove | fClusterBadChannelCheck | fRecalDistToBadChannels | fReClusterize;
499 Bool_t needRecalib = fCalibrateEnergy | fReClusterize;
500 Bool_t needTimecalib = fCalibrateTime | fReClusterize;
501 Bool_t needMisalign = fRecalClusPos | fReClusterize;
502 Bool_t needClusterizer = fReClusterize;
505 if (needBadChannels) {
506 Int_t fInitBC = InitBadChannels();
508 AliError("InitBadChannels returned false, returning");
510 AliWarning("InitBadChannels OK");
512 AliWarning(Form("No external hot channel set: %d - %s", event->GetRunNumber(), fFilepass.Data()));
515 // init recalibration factors
517 Int_t fInitRecalib = InitRecalib();
519 AliError("InitRecalib returned false, returning");
521 AliWarning("InitRecalib OK");
523 AliWarning(Form("No recalibration available: %d - %s", event->GetRunNumber(), fFilepass.Data()));
525 Int_t fInitRunDepRecalib = InitRunDepRecalib();
526 if (fInitRunDepRecalib==0)
527 AliError("InitrunDepRecalib returned false, returning");
528 if (fInitRunDepRecalib==1)
529 AliWarning("InitRecalib OK");
530 if (fInitRunDepRecalib>1)
531 AliWarning(Form("No Temperature recalibration available: %d - %s", event->GetRunNumber(), fFilepass.Data()));
534 // init time calibration
536 Int_t initTC = InitTimeCalibration();
538 AliError("InitTimeCalibration returned false, returning");
540 fCalibrateTimeParamAvailable = kTRUE;
541 AliWarning("InitTimeCalib OK");
544 AliWarning(Form("No external time calibration available: %d - %s", event->GetRunNumber(), fFilepass.Data()));
547 // init misalignment matrix
549 if (!InitMisalignMatrix())
550 AliError("InitMisalignmentMatrix returned false, returning");
552 AliWarning("InitMisalignMatrix OK");
555 // initiate reco params with some defaults
556 // will not overwrite, if those have been provided by user
558 Int_t initRC = InitRecParam();
561 AliInfo("Defaults reco params loaded.");
563 AliWarning("User defined reco params.");
567 if (needClusterizer) {
568 if (!InitClusterization())
569 AliError("InitClusterization returned false, returning");
571 AliWarning("InitClusterization OK");
575 fEMCALRecoUtils->Print("");
578 // disable implied switches -------------------------------------------------
579 // AliEMCALRecoUtils or clusterizer functions alredy include some
580 // recalibration so based on those implied callibration te switches
581 // here are disabled to avoid duplication
583 // clusterizer does cluster energy recalibration, position recomputation
586 fReCalibCluster = kFALSE;
587 fRecalClusPos = kFALSE;
588 fRecalShowerShape = kFALSE;
591 // CONFIGURE THE RECO UTILS -------------------------------------------------
592 // configure the reco utils
593 // this option does energy recalibration
594 if (fCalibrateEnergy)
595 fEMCALRecoUtils->SwitchOnRecalibration();
597 fEMCALRecoUtils->SwitchOffRecalibration();
599 // allows time calibration
601 fEMCALRecoUtils->SwitchOnTimeRecalibration();
603 fEMCALRecoUtils->SwitchOffTimeRecalibration();
605 // allows to zero bad cells
607 fEMCALRecoUtils->SwitchOnBadChannelsRemoval();
609 fEMCALRecoUtils->SwitchOffBadChannelsRemoval();
611 // distance to bad channel recalibration
612 if (fRecalDistToBadChannels)
613 fEMCALRecoUtils->SwitchOnDistToBadChannelRecalculation();
615 fEMCALRecoUtils->SwitchOffDistToBadChannelRecalculation();
617 // exclude exotic cells
618 if (fRejectExoticCells)
619 fEMCALRecoUtils->SwitchOnRejectExoticCell();
621 fEMCALRecoUtils->SwitchOffRejectExoticCell();
623 // exclude clusters with exotic cells
624 if (fRejectExoticClusters)
625 fEMCALRecoUtils->SwitchOnRejectExoticCluster();
627 fEMCALRecoUtils->SwitchOffRejectExoticCluster();
629 // TODO: not implemented switches
630 // SwitchOnClusterEnergySmearing
631 // SwitchOnRunDepCorrection
633 // START PROCESSING ---------------------------------------------------------
634 // Test if cells present
635 AliVCaloCells *cells= event->GetEMCALCells();
636 if (cells->GetNumberOfCells()<=0)
639 AliWarning(Form("Number of EMCAL cells = %d, returning", cells->GetNumberOfCells()));
644 AliInfo(Form("Re-calibrate cluster %d\n",fReCalibCluster));
646 // mark the cells not recalibrated in case of selected
647 // time, energy recalibration or bad channel removal
648 if (fCalibrateEnergy || fCalibrateTime || fBadCellRemove)
649 fEMCALRecoUtils->ResetCellsCalibrated();
651 // CELL RECALIBRATION -------------------------------------------------------
652 // cell objects will be updated
653 // the cell calibrations are also processed locally any time those are needed
654 // in case that the cell objects are not to be updated here for later use
657 // includes exotic cell check in the UpdateCells function - is not provided
661 // switch off recalibrations so those are not done multiple times
662 // this is just for safety, the recalibrated flag of cell object
663 // should not allow for farther processing anyways
664 fEMCALRecoUtils->SwitchOffRecalibration();
665 fEMCALRecoUtils->SwitchOffTimeRecalibration();
671 // RECLUSTERIZATION ---------------------------------------------------------
679 // Store good clusters
680 TClonesArray *clusArr = dynamic_cast<TClonesArray*>(event->FindListObject("caloClusters"));
682 clusArr = dynamic_cast<TClonesArray*>(event->FindListObject("CaloClusters"));
684 AliWarning(Form("No cluster array, number of cells in event = %d, returning", cells->GetNumberOfCells()));
688 // PROCESS SINGLE CLUSTER RECALIBRATION -------------------------------------
689 // now go through clusters one by one and process separate correction
690 // as those were defined or not
691 Int_t nclusters = clusArr->GetEntriesFast();
692 for (Int_t icluster=0; icluster < nclusters; ++icluster)
694 AliVCluster *clust = static_cast<AliVCluster*>(clusArr->At(icluster));
697 if (!clust->IsEMCAL())
700 // REMOVE CLUSTERS WITH BAD CELLS -----------------------------
701 if (fClusterBadChannelCheck) {
702 // careful, the the ClusterContainsBadChannel is dependent on
703 // SwitchOnBadChannelsRemoval, switching it ON automatically
704 // and returning to original value after processing
705 Bool_t badRemoval = fEMCALRecoUtils->IsBadChannelsRemovalSwitchedOn();
706 fEMCALRecoUtils->SwitchOnBadChannelsRemoval();
708 Bool_t badResult = fEMCALRecoUtils->ClusterContainsBadChannel(fEMCALGeo, clust->GetCellsAbsId(), clust->GetNCells());
710 // switch the bad channels removal back
712 fEMCALRecoUtils->SwitchOffBadChannelsRemoval();
716 delete clusArr->RemoveAt(icluster);
717 continue; //TODO is it really needed to remove it? Or should we flag it?
721 // REMOVE EXOTIC CLUSTERS -------------------------------------
722 // does process local cell recalibration energy and time without replacing
723 // the global cell values, in case of no cell recalib done yet
724 if (fRejectExoticClusters)
726 // careful, the IsExoticCluster is dependent on
727 // SwitchOnRejectExoticCell, switching it ON automatically
728 // and returning to original value after processing
729 Bool_t exRemoval = fEMCALRecoUtils->IsRejectExoticCell();
730 fEMCALRecoUtils->SwitchOnRejectExoticCell();
732 // get bunch crossing
733 Int_t bunchCrossNo = event->GetBunchCrossNumber();
735 Bool_t exResult = fEMCALRecoUtils->IsExoticCluster(clust, cells, bunchCrossNo);
737 // switch the exotic channels removal back
739 fEMCALRecoUtils->SwitchOffRejectExoticCell();
742 delete clusArr->RemoveAt(icluster);
743 continue; //TODO is it really needed to remove it? Or should we flag it?
747 // FIDUCIAL CUT -----------------------------------------------
749 // depends on SetNumberOfCellsFromEMCALBorder
750 // SwitchOnNoFiducialBorderInEMCALEta0
751 if (!fEMCALRecoUtils->CheckCellFiducialRegion(fEMCALGeo, clust, cells)){
752 delete clusArr->RemoveAt(icluster);
753 continue; //TODO it would be nice to store the distance
757 // CLUSTER ENERGY ---------------------------------------------
758 // does process local cell recalibration energy and time without replacing
759 // the global cell values, in case of no cell recalib done yet
760 if (fReCalibCluster) {
761 fEMCALRecoUtils->RecalibrateClusterEnergy(fEMCALGeo, clust, cells);
762 if (clust->E() < 1e-9) {
763 delete clusArr->RemoveAt(icluster);
768 // CLUSTER POSITION -------------------------------------------
769 // does process local cell energy recalibration, if enabled and cells
770 // not calibrated yet
772 fEMCALRecoUtils->RecalculateClusterPosition(fEMCALGeo, cells, clust);
774 // SHOWER SHAPE -----------------------------------------------
775 if (fRecalShowerShape)
776 fEMCALRecoUtils->RecalculateClusterShowerShapeParameters(fEMCALGeo, cells, clust);
778 // NONLINEARITY -----------------------------------------------
779 if (fDoNonLinearity) {
780 Float_t correctedEnergy = fEMCALRecoUtils->CorrectClusterEnergyLinearity(clust);
781 clust->SetE(correctedEnergy);
784 // DISTANCE TO BAD CHANNELS -----------------------------------
785 if (fRecalDistToBadChannels)
786 fEMCALRecoUtils->RecalculateClusterDistanceToBadChannel(fEMCALGeo, cells, clust);
794 // TRACK MATCHING -----------------------------------------------------------
795 if (!TGeoGlobalMagField::Instance()->GetField()) {
796 event->InitMagneticField();
799 fEMCALRecoUtils->FindMatches(event,0x0,fEMCALGeo);
800 fEMCALRecoUtils->SetClusterMatchedToTrack(event);
801 fEMCALRecoUtils->SetTracksMatchedToCluster(event);
804 //_____________________________________________________
805 Bool_t AliEMCALTenderSupply::InitMisalignMatrix()
807 // Initialising misalignment matrices
809 AliVEvent *event = GetEvent();
816 AliInfo("Misalignment matrix already set");
821 AliInfo("Initialising misalignment matrix");
823 if (fLoadGeomMatrices) {
824 for(Int_t mod=0; mod < fEMCALGeo->GetNumberOfSuperModules(); ++mod)
826 if (fEMCALMatrix[mod]){
828 fEMCALMatrix[mod]->Print();
829 fEMCALGeo->SetMisalMatrix(fEMCALMatrix[mod],mod);
832 fGeomMatrixSet = kTRUE;
836 Int_t runGM = event->GetRunNumber();
839 if (fMisalignSurvey == kdefault)
840 { //take default alignment corresponding to run no
841 AliOADBContainer emcalgeoCont(Form("emcal"));
842 emcalgeoCont.InitFromFile("$ALICE_ROOT/OADB/EMCAL/EMCALlocal2master.root",Form("AliEMCALgeo"));
843 mobj=(TObjArray*)emcalgeoCont.GetObject(runGM,"EmcalMatrices");
846 if (fMisalignSurvey == kSurveybyS)
847 { //take alignment at sector level
848 if (runGM <= 140000) { //2010 data
849 AliOADBContainer emcalgeoCont(Form("emcal2010"));
850 emcalgeoCont.InitFromFile("$ALICE_ROOT/OADB/EMCAL/EMCALlocal2master.root",Form("AliEMCALgeo"));
851 mobj=(TObjArray*)emcalgeoCont.GetObject(100,"survey10");
853 else if (runGM>140000)
854 { // 2011 LHC11a pass1 data
855 AliOADBContainer emcalgeoCont(Form("emcal2011"));
856 emcalgeoCont.InitFromFile("$ALICE_ROOT/OADB/EMCAL/EMCALlocal2master.root",Form("AliEMCALgeo"));
857 mobj=(TObjArray*)emcalgeoCont.GetObject(100,"survey11byS");
861 if (fMisalignSurvey == kSurveybyM)
862 { //take alignment at module level
863 if (runGM <= 140000) { //2010 data
864 AliOADBContainer emcalgeoCont(Form("emcal2010"));
865 emcalgeoCont.InitFromFile("$ALICE_ROOT/OADB/EMCAL/EMCALlocal2master.root",Form("AliEMCALgeo"));
866 mobj=(TObjArray*)emcalgeoCont.GetObject(100,"survey10");
868 else if (runGM>140000)
869 { // 2011 LHC11a pass1 data
870 AliOADBContainer emcalgeoCont(Form("emcal2011"));
871 emcalgeoCont.InitFromFile("$ALICE_ROOT/OADB/EMCAL/EMCALlocal2master.root",Form("AliEMCALgeo"));
872 mobj=(TObjArray*)emcalgeoCont.GetObject(100,"survey11byM");
877 AliFatal("Geometry matrix array not found");
881 for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++)
883 fEMCALMatrix[mod] = (TGeoHMatrix*) mobj->At(mod);
884 fEMCALGeo->SetMisalMatrix(fEMCALMatrix[mod],mod);
885 fEMCALMatrix[mod]->Print();
891 //_____________________________________________________
892 Int_t AliEMCALTenderSupply::InitBadChannels()
894 // Initialising bad channel maps
896 AliVEvent *event = GetEvent();
902 AliInfo("Initialising Bad channel map");
904 // init default maps first
905 if (!fEMCALRecoUtils->GetEMCALBadChannelStatusMapArray())
906 fEMCALRecoUtils->InitEMCALBadChannelStatusMap() ;
908 Int_t runBC = event->GetRunNumber();
910 AliOADBContainer *contBC = new AliOADBContainer("");
912 { //if fBasePath specified in the ->SetBasePath()
913 if (fDebugLevel>0) AliInfo(Form("Loading Bad Channels OADB from given path %s",fBasePath.Data()));
915 TFile *fbad=new TFile(Form("%s/EMCALBadChannels.root",fBasePath.Data()),"read");
916 if (!fbad || fbad->IsZombie())
918 AliFatal(Form("EMCALBadChannels.root was not found in the path provided: %s",fBasePath.Data()));
922 if (fbad) delete fbad;
924 contBC->InitFromFile(Form("%s/EMCALBadChannels.root",fBasePath.Data()),"AliEMCALBadChannels");
927 { // Else choose the one in the $ALICE_ROOT directory
928 if (fDebugLevel>0) AliInfo("Loading Bad Channels OADB from $ALICE_ROOT/OADB/EMCAL");
930 TFile *fbad=new TFile("$ALICE_ROOT/OADB/EMCAL/EMCALBadChannels.root","read");
931 if (!fbad || fbad->IsZombie())
933 AliFatal("$ALICE_ROOT/OADB/EMCAL/EMCALBadChannels.root was not found");
937 if (fbad) delete fbad;
939 contBC->InitFromFile("$ALICE_ROOT/OADB/EMCAL/EMCALBadChannels.root","AliEMCALBadChannels");
942 TObjArray *arrayBC=(TObjArray*)contBC->GetObject(runBC);
945 AliError(Form("No external hot channel set for run number: %d", runBC));
949 Int_t sms = fEMCALGeo->GetEMCGeometry()->GetNumberOfSuperModules();
950 for (Int_t i=0; i<sms; ++i)
952 TH2I *h = fEMCALRecoUtils->GetEMCALChannelStatusMap(i);
955 h=(TH2I*)arrayBC->FindObject(Form("EMCALBadChannelMap_Mod%d",i));
959 AliError(Form("Can not get EMCALBadChannelMap_Mod%d",i));
963 fEMCALRecoUtils->SetEMCALChannelStatusMap(i,h);
968 //_____________________________________________________
969 Int_t AliEMCALTenderSupply::InitRecalib()
971 // Initialising recalibration factors.
973 AliVEvent *event = GetEvent();
979 AliInfo("Initialising recalibration factors");
981 // init default maps first
982 if (!fEMCALRecoUtils->GetEMCALRecalibrationFactorsArray())
983 fEMCALRecoUtils->InitEMCALRecalibrationFactors() ;
985 Int_t runRC = event->GetRunNumber();
987 AliOADBContainer *contRF=new AliOADBContainer("");
989 { //if fBasePath specified in the ->SetBasePath()
990 if (fDebugLevel>0) AliInfo(Form("Loading Recalib OADB from given path %s",fBasePath.Data()));
992 TFile *fRecalib= new TFile(Form("%s/EMCALRecalib.root",fBasePath.Data()),"read");
993 if (!fRecalib || fRecalib->IsZombie())
995 AliFatal(Form("EMCALRecalib.root not found in %s",fBasePath.Data()));
999 if (fRecalib) delete fRecalib;
1001 contRF->InitFromFile(Form("%s/EMCALRecalib.root",fBasePath.Data()),"AliEMCALRecalib");
1004 { // Else choose the one in the $ALICE_ROOT directory
1005 if (fDebugLevel>0) AliInfo("Loading Recalib OADB from $ALICE_ROOT/OADB/EMCAL");
1007 TFile *fRecalib= new TFile("$ALICE_ROOT/OADB/EMCAL/EMCALRecalib.root","read");
1008 if (!fRecalib || fRecalib->IsZombie())
1010 AliFatal("$ALICE_ROOT/OADB/EMCAL/EMCALRecalib.root was not found");
1014 if (fRecalib) delete fRecalib;
1016 contRF->InitFromFile("$ALICE_ROOT/OADB/EMCAL/EMCALRecalib.root","AliEMCALRecalib");
1019 TObjArray *recal=(TObjArray*)contRF->GetObject(runRC);
1022 AliError(Form("No Objects for run: %d",runRC));
1026 TObjArray *recalpass=(TObjArray*)recal->FindObject(fFilepass);
1029 AliError(Form("No Objects for run: %d - %s",runRC,fFilepass.Data()));
1033 TObjArray *recalib=(TObjArray*)recalpass->FindObject("Recalib");
1036 AliError(Form("No Recalib histos found for %d - %s",runRC,fFilepass.Data()));
1040 if (fDebugLevel>0) recalib->Print();
1042 Int_t sms = fEMCALGeo->GetEMCGeometry()->GetNumberOfSuperModules();
1043 for (Int_t i=0; i<sms; ++i)
1045 TH2F *h = fEMCALRecoUtils->GetEMCALChannelRecalibrationFactors(i);
1048 h = (TH2F*)recalib->FindObject(Form("EMCALRecalFactors_SM%d",i));
1051 AliError(Form("Could not load EMCALRecalFactors_SM%d",i));
1055 fEMCALRecoUtils->SetEMCALChannelRecalibrationFactors(i,h);
1060 //_____________________________________________________
1061 Int_t AliEMCALTenderSupply::InitRunDepRecalib()
1063 // Initialising recalibration factors.
1065 AliVEvent *event = GetEvent();
1071 AliInfo("Initialising recalibration factors");
1073 // init default maps first
1074 if (!fEMCALRecoUtils->GetEMCALRecalibrationFactorsArray())
1075 fEMCALRecoUtils->InitEMCALRecalibrationFactors() ;
1077 Int_t runRC = event->GetRunNumber();
1079 AliOADBContainer *contRF=new AliOADBContainer("");
1081 { //if fBasePath specified in the ->SetBasePath()
1082 if (fDebugLevel>0) AliInfo(Form("Loading Recalib OADB from given path %s",fBasePath.Data()));
1084 TFile *fRunDepRecalib= new TFile(Form("%s/EMCALTemperatureCorrCalib.root",fBasePath.Data()),"read");
1085 if (!fRunDepRecalib || fRunDepRecalib->IsZombie())
1087 AliFatal(Form("EMCALTemperatureCorrCalib.root not found in %s",fBasePath.Data()));
1091 if (fRunDepRecalib) delete fRunDepRecalib;
1093 contRF->InitFromFile(Form("%s/EMCALTemperatureCorrCalib.root",fBasePath.Data()),"AliEMCALRunDepTempCalibCorrections");
1096 { // Else choose the one in the $ALICE_ROOT directory
1097 if (fDebugLevel>0) AliInfo("Loading Recalib OADB from $ALICE_ROOT/OADB/EMCAL");
1099 TFile *fRunDepRecalib= new TFile("$ALICE_ROOT/OADB/EMCAL/EMCALTemperatureCorrCalib.root","read");
1100 if (!fRunDepRecalib || fRunDepRecalib->IsZombie())
1102 AliFatal("$ALICE_ROOT/OADB/EMCAL/EMCALTemperatureCorrCalib.root was not found");
1106 if (fRunDepRecalib) delete fRunDepRecalib;
1108 contRF->InitFromFile("$ALICE_ROOT/OADB/EMCAL/EMCALTemperatureCorrCalib.root","AliEMCALRunDepTempCalibCorrections");
1111 TH1S *rundeprecal=(TH1S*)contRF->GetObject(runRC);
1114 AliWarning(Form("No TemperatureCorrCalib Objects for run: %d",runRC));
1115 // let's get the closest runnumber instead then..
1118 Int_t maxEntry = contRF->GetNumberOfEntries();
1120 while ((ic < maxEntry) && (contRF->UpperLimit(ic) < runRC)) {
1125 Int_t closest = lower;
1126 if ((ic<maxEntry) &&
1127 (contRF->LowerLimit(ic)-runRC) < (runRC - contRF->UpperLimit(lower))) {
1131 AliWarning(Form("TemperatureCorrCalib Objects found closest id %d from run: %d", closest, contRF->LowerLimit(closest)));
1132 rundeprecal = (TH1S*) contRF->GetObjectByIndex(closest);
1135 if (fDebugLevel>0) rundeprecal->Print();
1137 Int_t nSM = fEMCALGeo->GetEMCGeometry()->GetNumberOfSuperModules();
1139 for (Int_t ism=0; ism<nSM; ++ism)
1141 for (Int_t icol=0; icol<48; ++icol)
1143 for (Int_t irow=0; irow<24; ++irow)
1145 Float_t factor = fEMCALRecoUtils->GetEMCALChannelRecalibrationFactor(ism,icol,irow);
1147 Int_t absID = fEMCALGeo->GetAbsCellIdFromCellIndexes(ism, irow, icol); // original calibration factor
1148 factor *= rundeprecal->GetBinContent(absID) / 10000. ; // correction dependent on T
1150 fEMCALRecoUtils->SetEMCALChannelRecalibrationFactor(ism,icol,irow,factor);
1159 //_____________________________________________________
1160 Int_t AliEMCALTenderSupply::InitTimeCalibration()
1162 // Initialising bad channel maps
1163 AliVEvent *event = GetEvent();
1169 AliInfo("Initialising time calibration map");
1171 // init default maps first
1172 if (!fEMCALRecoUtils->GetEMCALTimeRecalibrationFactorsArray())
1173 fEMCALRecoUtils->InitEMCALTimeRecalibrationFactors() ;
1175 Int_t runBC = event->GetRunNumber();
1177 AliOADBContainer *contBC = new AliOADBContainer("");
1179 { //if fBasePath specified in the ->SetBasePath()
1180 if (fDebugLevel>0) AliInfo(Form("Loading time calibration OADB from given path %s",fBasePath.Data()));
1182 TFile *fbad=new TFile(Form("%s/EMCALTimeCalib.root",fBasePath.Data()),"read");
1183 if (!fbad || fbad->IsZombie())
1185 AliFatal(Form("EMCALTimeCalib.root was not found in the path provided: %s",fBasePath.Data()));
1189 if (fbad) delete fbad;
1191 contBC->InitFromFile(Form("%s/EMCALTimeCalib.root",fBasePath.Data()),"AliEMCALTimeCalib");
1194 { // Else choose the one in the $ALICE_ROOT directory
1195 if (fDebugLevel>0) AliInfo("Loading time calibration OADB from $ALICE_ROOT/OADB/EMCAL");
1197 TFile *fbad=new TFile("$ALICE_ROOT/OADB/EMCAL/EMCALTimeCalib.root","read");
1198 if (!fbad || fbad->IsZombie())
1200 AliFatal("$ALICE_ROOT/OADB/EMCAL/EMCALTimeCalib.root was not found");
1204 if (fbad) delete fbad;
1206 contBC->InitFromFile("$ALICE_ROOT/OADB/EMCAL/EMCALTimeCalib.root","AliEMCALTimeCalib");
1209 TObjArray *arrayBC=(TObjArray*)contBC->GetObject(runBC);
1212 AliError(Form("No external time calibration set for run number: %d", runBC));
1216 // Here, it looks for a specific pass
1217 TString pass = fFilepass;
1218 if (fFilepass=="calo_spc") pass ="pass1";
1219 TObjArray *arrayBCpass=(TObjArray*)arrayBC->FindObject(pass);
1222 AliError(Form("No external time calibration set for: %d -%s", runBC,pass.Data()));
1226 if (fDebugLevel>0) arrayBCpass->Print();
1228 for(Int_t i = 0; i < 4; i++)
1230 TH1F *h = fEMCALRecoUtils->GetEMCALChannelTimeRecalibrationFactors(i);
1234 h = (TH1F*)arrayBCpass->FindObject(Form("hAllTimeAvBC%d",i));
1238 AliError(Form("Can not get hAllTimeAvBC%d",i));
1242 fEMCALRecoUtils->SetEMCALChannelTimeRecalibrationFactors(i,h);
1247 //_____________________________________________________
1248 void AliEMCALTenderSupply::UpdateCells()
1250 //Remove bad cells from the cell list
1251 //Recalibrate energy and time cells
1252 //This is required for later reclusterization
1254 AliVEvent *event = GetEvent();
1256 if (!event) return ;
1258 AliVCaloCells *cells = event->GetEMCALCells();
1259 Int_t bunchCrossNo = event->GetBunchCrossNumber();
1261 fEMCALRecoUtils->RecalibrateCells(cells, bunchCrossNo);
1263 // remove exotic cells - loop through cells and zero the exotic ones
1264 // just like with bad cell rejection in reco utils (inside RecalibrateCells)
1265 if (fRejectExoticCells)
1272 Bool_t isExot = kFALSE;
1274 // loop through cells
1275 Int_t nEMcell = cells->GetNumberOfCells() ;
1276 for (Int_t iCell = 0; iCell < nEMcell; iCell++)
1278 cells->GetCell(iCell, absId, ecell, tcell, mclabel, efrac);
1280 isExot = fEMCALRecoUtils->IsExoticCell(absId, cells, bunchCrossNo);
1283 cells->SetCell(iCell, absId, 0.0, -1.0, mclabel, efrac);
1285 } // reject exotic cells
1290 //_____________________________________________________
1291 TString AliEMCALTenderSupply::GetBeamType()
1293 // Get beam type : pp-AA-pA
1294 // ESDs have it directly, AODs get it from hardcoded run number ranges
1296 AliVEvent *event = GetEvent();
1299 AliError("Couldn't retrieve event!");
1305 AliESDEvent *esd = dynamic_cast<AliESDEvent*>(event);
1307 const AliESDRun *run = esd->GetESDRun();
1308 beamType = run->GetBeamType();
1312 Int_t runNumber = event->GetRunNumber();
1313 if ((runNumber >= 136851 && runNumber <= 139517) // LHC10h
1314 || (runNumber >= 166529 && runNumber <= 170593)) // LHC11h
1327 //_____________________________________________________
1328 Int_t AliEMCALTenderSupply::InitRecParam()
1330 // Init reco params if not yet exist (probably shipped by the user already)
1335 TString beamType = GetBeamType();
1337 // set some default reco params
1338 fRecParam = new AliEMCALRecParam();
1339 fRecParam->SetClusteringThreshold(0.100);
1340 fRecParam->SetMinECut(0.050);
1342 if (!fCalibrateTimeParamAvailable) {
1343 fRecParam->SetTimeCut(250*1.e-9);
1344 fRecParam->SetTimeMin(425*1.e-9);
1345 fRecParam->SetTimeMax(825*1.e-9);
1347 fRecParam->SetTimeCut(100*1.e-9);
1348 fRecParam->SetTimeMin(-50*1.e-9);
1349 fRecParam->SetTimeMax(50*1.e-9);
1352 if (beamType == "A-A") {
1353 fRecParam->SetClusterizerFlag(AliEMCALRecParam::kClusterizerv2);
1355 fRecParam->SetClusterizerFlag(AliEMCALRecParam::kClusterizerv1);
1361 //_____________________________________________________
1362 Bool_t AliEMCALTenderSupply::InitClusterization()
1364 // Initialing clusterization/unfolding algorithm and set all the needed parameters.
1366 AliVEvent *event = GetEvent();
1372 AliInfo(Form("Initialising reclustering parameters: Clusterizer type: %d",fRecParam->GetClusterizerFlag()));
1374 //---setup clusterizer
1376 // avoid deleting fDigitsArr
1377 fClusterizer->SetDigitsArr(0);
1378 delete fClusterizer;
1381 if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv1)
1382 fClusterizer = new AliEMCALClusterizerv1 (fEMCALGeo);
1383 else if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv2)
1384 fClusterizer = new AliEMCALClusterizerv2(fEMCALGeo);
1385 else if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerNxN)
1387 AliEMCALClusterizerNxN *clusterizer = new AliEMCALClusterizerNxN(fEMCALGeo);
1388 clusterizer->SetNRowDiff(fRecParam->GetNRowDiff());
1389 clusterizer->SetNColDiff(fRecParam->GetNColDiff());
1390 fClusterizer = clusterizer;
1394 AliFatal(Form("Clusterizer < %d > not available", fRecParam->GetClusterizerFlag()));
1398 // Set the clustering parameters
1399 fClusterizer->SetECAClusteringThreshold(fRecParam->GetClusteringThreshold());
1400 fClusterizer->SetECALogWeight (fRecParam->GetW0() );
1401 fClusterizer->SetMinECut (fRecParam->GetMinECut() );
1402 fClusterizer->SetUnfolding (fRecParam->GetUnfold() );
1403 fClusterizer->SetECALocalMaxCut (fRecParam->GetLocMaxCut() );
1404 fClusterizer->SetTimeCut (fRecParam->GetTimeCut() );
1405 fClusterizer->SetTimeMin (fRecParam->GetTimeMin() );
1406 fClusterizer->SetTimeMax (fRecParam->GetTimeMax() );
1407 fClusterizer->SetInputCalibrated (kTRUE );
1408 fClusterizer->SetJustClusters (kTRUE );
1410 // In case of unfolding after clusterization is requested, set the corresponding parameters
1411 if (fRecParam->GetUnfold())
1413 for (Int_t i = 0; i < 8; ++i)
1415 fClusterizer->SetSSPars(i, fRecParam->GetSSPars(i));
1417 for (Int_t i = 0; i < 3; ++i)
1419 fClusterizer->SetPar5 (i, fRecParam->GetPar5(i));
1420 fClusterizer->SetPar6 (i, fRecParam->GetPar6(i));
1422 fClusterizer->InitClusterUnfolding();
1425 fClusterizer->SetDigitsArr(fDigitsArr);
1426 fClusterizer->SetOutput(0);
1427 fClusterArr = const_cast<TObjArray *>(fClusterizer->GetRecPoints());
1431 //_____________________________________________________
1432 void AliEMCALTenderSupply::FillDigitsArray()
1434 // Fill digits from cells to a TClonesArray.
1436 AliVEvent *event = GetEvent();
1441 // In case of MC productions done before aliroot tag v5-02-Rev09
1442 // assing the cluster label to all the cells belonging to this cluster
1444 Int_t cellLabels[12672];
1445 if (fSetCellMCLabelFromCluster)
1447 for (Int_t i = 0; i < 12672; i++)
1449 cellLabels [i] = 0 ;
1450 fOrgClusterCellId[i] =-1 ;
1453 Int_t nClusters = event->GetNumberOfCaloClusters();
1454 for (Int_t i = 0; i < nClusters; i++)
1456 AliVCluster *clus = event->GetCaloCluster(i);
1458 if (!clus) continue;
1460 if (!clus->IsEMCAL()) continue ;
1462 Int_t label = clus->GetLabel();
1463 UShort_t * index = clus->GetCellsAbsId() ;
1465 for(Int_t icell=0; icell < clus->GetNCells(); icell++)
1467 cellLabels [index[icell]] = label;
1468 fOrgClusterCellId[index[icell]] = i ; // index of the original cluster
1473 fDigitsArr->Clear("C");
1474 AliVCaloCells *cells = event->GetEMCALCells();
1475 Int_t ncells = cells->GetNumberOfCells();
1476 for (Int_t icell = 0, idigit = 0; icell < ncells; ++icell)
1478 Double_t cellAmplitude=0, cellTime=0, efrac = 0;
1479 Short_t cellNumber=0;
1482 if (cells->GetCell(icell, cellNumber, cellAmplitude, cellTime, mcLabel, efrac) != kTRUE)
1485 // Do not add if energy already too low (some cells set to 0 if bad channels)
1486 if (cellAmplitude < fRecParam->GetMinECut())
1489 // If requested, do not include exotic cells
1490 if (fEMCALRecoUtils->IsExoticCell(cellNumber,cells,event->GetBunchCrossNumber()))
1493 if (fSetCellMCLabelFromCluster) mcLabel = cellLabels[cellNumber];
1494 else if (fRemapMCLabelForAODs ) RemapMCLabelForAODs(mcLabel);
1496 if (mcLabel > 0 && efrac < 1e-6) efrac = 1;
1498 new((*fDigitsArr)[idigit]) AliEMCALDigit(mcLabel, mcLabel, cellNumber,
1499 (Float_t)cellAmplitude, (Float_t)cellTime,
1500 AliEMCALDigit::kHG,idigit, 0, 0, efrac*cellAmplitude);
1505 //_____________________________________________________
1506 void AliEMCALTenderSupply::Clusterize()
1510 fClusterizer->Digits2Clusters("");
1513 //_____________________________________________________
1514 void AliEMCALTenderSupply::UpdateClusters()
1516 // Update ESD cluster list.
1518 AliVEvent *event = GetEvent();
1523 TClonesArray *clus = dynamic_cast<TClonesArray*>(event->FindListObject("caloClusters"));
1525 clus = dynamic_cast<TClonesArray*>(event->FindListObject("CaloClusters"));
1528 AliError(" Null pointer to calo clusters array, returning");
1532 // Before destroying the orignal list, assign to the rec points the MC labels
1533 // of the original clusters, if requested
1534 if (fSetCellMCLabelFromCluster == 2)
1535 SetClustersMCLabelFromOriginalClusters() ;
1537 Int_t nents = clus->GetEntriesFast();
1538 for (Int_t i=0; i < nents; ++i)
1540 AliVCluster *c = dynamic_cast<AliVCluster*>(clus->At(i));
1545 delete clus->RemoveAt(i);
1551 RecPoints2Clusters(clus);
1554 //_____________________________________________________
1555 void AliEMCALTenderSupply::RecPoints2Clusters(TClonesArray *clus)
1557 // Convert AliEMCALRecoPoints to AliESDCaloClusters/AliAODCaloClusters.
1558 // Cluster energy, global position, cells and their amplitude fractions are restored.
1560 AliVEvent *event = GetEvent();
1565 Int_t ncls = fClusterArr->GetEntriesFast();
1566 for(Int_t i=0, nout=clus->GetEntriesFast(); i < ncls; ++i)
1568 AliEMCALRecPoint *recpoint = static_cast<AliEMCALRecPoint*>(fClusterArr->At(i));
1570 Int_t ncellsTrue = 0;
1571 const Int_t ncells = recpoint->GetMultiplicity();
1572 UShort_t absIds[ncells];
1573 Double32_t ratios[ncells];
1574 Int_t *dlist = recpoint->GetDigitsList();
1575 Float_t *elist = recpoint->GetEnergiesList();
1576 for (Int_t c = 0; c < ncells; ++c)
1578 AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(dlist[c]));
1579 absIds[ncellsTrue] = digit->GetId();
1580 ratios[ncellsTrue] = elist[c]/digit->GetAmplitude();
1581 if (ratios[ncellsTrue] < 0.001)
1588 AliWarning("Skipping cluster with no cells");
1592 // calculate new cluster position
1594 recpoint->GetGlobalPosition(gpos);
1598 AliVCluster *c = static_cast<AliVCluster*>(clus->New(nout++));
1600 c->SetType(AliVCluster::kEMCALClusterv1);
1601 c->SetE(recpoint->GetEnergy());
1603 c->SetNCells(ncellsTrue);
1604 c->SetDispersion(recpoint->GetDispersion());
1605 c->SetEmcCpvDistance(-1); //not yet implemented
1606 c->SetChi2(-1); //not yet implemented
1607 c->SetTOF(recpoint->GetTime()) ; //time-of-flight
1608 c->SetNExMax(recpoint->GetNExMax()); //number of local maxima
1609 Float_t elipAxis[2];
1610 recpoint->GetElipsAxis(elipAxis);
1611 c->SetM02(elipAxis[0]*elipAxis[0]) ;
1612 c->SetM20(elipAxis[1]*elipAxis[1]) ;
1613 c->SetCellsAbsId(absIds);
1614 c->SetCellsAmplitudeFraction(ratios);
1617 Int_t parentMult = 0;
1618 Int_t *parentList = recpoint->GetParents(parentMult);
1619 if (parentMult > 0) c->SetLabel(parentList, parentMult);
1624 //___________________________________________________________
1625 void AliEMCALTenderSupply::RemapMCLabelForAODs(Int_t & label)
1627 // MC label for Cells not remapped after ESD filtering, do it here.
1629 if (label < 0) return;
1631 AliAODEvent * evt = dynamic_cast<AliAODEvent*> (GetEvent()) ;
1634 TClonesArray * arr = dynamic_cast<TClonesArray*>(evt->FindListObject("mcparticles")) ;
1637 if (label < arr->GetEntriesFast())
1639 AliAODMCParticle * particle = dynamic_cast<AliAODMCParticle *>(arr->At(label));
1640 if (!particle) return ;
1642 if (label == particle->Label()) return ; // label already OK
1645 // loop on the particles list and check if there is one with the same label
1646 for (Int_t ind = 0; ind < arr->GetEntriesFast(); ind++)
1648 AliAODMCParticle * particle = dynamic_cast<AliAODMCParticle *>(arr->At(ind));
1649 if (!particle) continue ;
1651 if (label == particle->Label())
1661 //_____________________________________________________________________________________________
1662 void AliEMCALTenderSupply::SetClustersMCLabelFromOriginalClusters()
1664 // Get the original clusters that contribute to the new rec point cluster,
1665 // assign the labels of such clusters to the new rec point cluster.
1666 // Only approximatedly valid when output are V1 clusters, or input/output clusterizer
1667 // are the same handle with care
1668 // Copy from same method in AliAnalysisTaskEMCALClusterize, but here modify the recpoint and
1669 // not the output calocluster
1671 Int_t ncls = fClusterArr->GetEntriesFast();
1672 for(Int_t irp=0; irp < ncls; ++irp)
1674 TArrayI clArray(300) ; //Weird if more than a few clusters are in the origin ...
1677 Int_t nLabTotOrg = 0;
1681 AliEMCALRecPoint *clus = static_cast<AliEMCALRecPoint*>(fClusterArr->At(irp));
1683 //Find the clusters that originally had the cells
1684 const Int_t ncells = clus->GetMultiplicity();
1685 Int_t *digList = clus->GetDigitsList();
1687 for (Int_t iLoopCell = 0 ; iLoopCell < ncells ; iLoopCell++)
1689 AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(digList[iLoopCell]));
1690 Int_t idCell = digit->GetId();
1694 Int_t idCluster = fOrgClusterCellId[idCell];
1696 for (Int_t icl =0; icl < nClu; icl++)
1698 if (((Int_t)clArray.GetAt(icl))==-1) continue;
1699 if (idCluster == ((Int_t)clArray.GetAt(icl))) set = kFALSE;
1701 if (set && idCluster >= 0)
1703 clArray.SetAt(idCluster,nClu++);
1704 nLabTotOrg+=(GetEvent()->GetCaloCluster(idCluster))->GetNLabels();
1706 //Search highest E cluster
1707 AliVCluster * clOrg = GetEvent()->GetCaloCluster(idCluster);
1708 if (emax < clOrg->E())
1717 // Put the first in the list the cluster with highest energy
1718 if (idMax != ((Int_t)clArray.GetAt(0))) // Max not at first position
1720 Int_t maxIndex = -1;
1721 Int_t firstCluster = ((Int_t)clArray.GetAt(0));
1722 for (Int_t iLoopCluster = 0 ; iLoopCluster < nClu ; iLoopCluster++)
1724 if (idMax == ((Int_t)clArray.GetAt(iLoopCluster))) maxIndex = iLoopCluster;
1727 if (firstCluster >=0 && idMax >=0)
1729 clArray.SetAt(idMax,0);
1730 clArray.SetAt(firstCluster,maxIndex);
1734 // Get the labels list in the original clusters, assign all to the new cluster
1735 TArrayI clMCArray(nLabTotOrg) ;
1739 for (Int_t iLoopCluster = 0 ; iLoopCluster < nClu ; iLoopCluster++)
1741 Int_t idCluster = (Int_t) clArray.GetAt(iLoopCluster);
1742 AliVCluster * clOrg = GetEvent()->GetCaloCluster(idCluster);
1743 Int_t nLab = clOrg->GetNLabels();
1745 for (Int_t iLab = 0 ; iLab < nLab ; iLab++)
1747 Int_t lab = clOrg->GetLabelAt(iLab) ;
1751 for(Int_t iLabTot =0; iLabTot < nLabTot; iLabTot++)
1753 if (lab == ((Int_t)clMCArray.GetAt(iLabTot))) set = kFALSE;
1755 if (set) clMCArray.SetAt(lab,nLabTot++);
1760 // Set the final list of labels to rec point
1762 Int_t *labels = new Int_t[nLabTot];
1763 for(Int_t il = 0; il < nLabTot; il++) labels[il] = clMCArray.GetArray()[il];
1764 clus->SetParents(nLabTot,labels);
1766 } // rec point array
1769 //_____________________________________________________
1770 void AliEMCALTenderSupply::GetPass()
1772 // Get passx from filename
1774 AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
1775 fInputTree = mgr->GetTree();
1779 AliError("Pointer to tree = 0, returning");
1783 fInputFile = fInputTree->GetCurrentFile();
1785 AliError("Null pointer input file, returning");
1789 TString fname(fInputFile->GetName());
1790 if (fname.Contains("pass1")) fFilepass = TString("pass1");
1791 else if (fname.Contains("pass2")) fFilepass = TString("pass2");
1792 else if (fname.Contains("pass3")) fFilepass = TString("pass3");
1793 else if (fname.Contains("pass4")) fFilepass = TString("pass4");
1794 else if (fname.Contains("pass5")) fFilepass = TString("pass5");
1795 else if (fname.Contains("LHC11c") &&
1796 fname.Contains("spc_calo")) fFilepass = TString("spc_calo");
1797 else if (fname.Contains("calo") || fname.Contains("high_lumi"))
1800 //printf("AliEMCALTenderSupply::GetPass() - Path contains <calo> or <high-lumi>, set as <pass1>\n");
1801 fFilepass = TString("pass1");
1805 AliError(Form("Pass number string not found: %s", fname.Data()));