]> git.uio.no Git - u/mrichter/AliRoot.git/blame - ANALYSIS/TenderSupplies/AliEMCALTenderSupply.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Coverity: Uninitialized data member in default ctor, fEMCALMatrix
[u/mrichter/AliRoot.git] / ANALYSIS / TenderSupplies / AliEMCALTenderSupply.cxx
CommitLineData
6d67b5a7 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16
17///////////////////////////////////////////////////////////////////////////////
18// //
19// EMCAL tender, apply corrections to EMCAl clusters //
766cc9de 20// and do track matching. //
21// Author: Deepa Thomas (Utrecht University) //
6d67b5a7 22// //
23///////////////////////////////////////////////////////////////////////////////
24
766cc9de 25#include <TROOT.h>
26#include <TFile.h>
27#include <TTree.h>
28#include <TInterpreter.h>
29#include <TObjArray.h>
30#include <TClonesArray.h>
31#include <TGeoGlobalMagField.h>
6d67b5a7 32#include <AliLog.h>
33#include <AliESDEvent.h>
34#include <AliAnalysisManager.h>
35#include <AliTender.h>
36#include "AliOADBContainer.h"
b20bc239 37#include "AliCDBManager.h"
38#include "AliCDBStorage.h"
39#include "AliCDBEntry.h"
6d67b5a7 40#include "AliMagF.h"
6d67b5a7 41#include "AliESDCaloCluster.h"
42#include "AliEMCALTenderSupply.h"
6d67b5a7 43#include "AliEMCALGeometry.h"
44#include "AliEMCALRecoUtils.h"
b20bc239 45#include "AliEMCALClusterizer.h"
46#include "AliEMCALRecParam.h"
47#include "AliEMCALRecPoint.h"
48#include "AliEMCALAfterBurnerUF.h"
49#include "AliEMCALClusterizerNxN.h"
50#include "AliEMCALClusterizerv1.h"
51#include "AliEMCALClusterizerv2.h"
52#include "AliEMCALDigit.h"
6d67b5a7 53
3b7d451e 54ClassImp(AliEMCALTenderSupply)
55
6d67b5a7 56AliEMCALTenderSupply::AliEMCALTenderSupply() :
6b7df55d 57AliTenderSupply()
58,fEMCALGeo(0x0)
59,fEMCALGeoName("EMCAL_COMPLETEV1")
60,fEMCALRecoUtils(0)
61,fConfigName("")
62,fDebugLevel(0)
63,fNonLinearFunc(AliEMCALRecoUtils::kNoCorrection)
64,fNonLinearThreshold(30)
65,fReCalibCluster(kFALSE)
66,fReCalibCell(kFALSE)
67,fRecalClusPos(kFALSE)
68,fFiducial(kFALSE)
69,fNCellsFromEMCALBorder(1)
70,fRecalDistToBadChannels(kFALSE)
71,fInputTree(0)
72,fInputFile(0)
73,fFilepass(0)
74,fMass(0.139)
75,fStep(50)
76,fCutEtaPhiSum(kTRUE)
77,fCutEtaPhiSeparate(kFALSE)
78,fRcut(0.05)
79,fEtacut(0.025)
80,fPhicut(0.05)
81,fBasePath(".")
82,fReClusterize(kFALSE)
83,fClusterizer(0)
84,fGeomMatrixSet(kFALSE)
85,fLoadGeomMatrices(kFALSE)
86,fRecParam(0)
87,fOCDBpath(" ")
88,fUnfolder(0)
89,fDigitsArr(0)
90,fClusterArr(0)
6d67b5a7 91{
766cc9de 92 // Default constructor.
6b7df55d 93 for(Int_t i = 0; i < 10; i++) fEMCALMatrix[i] = 0 ;
6d67b5a7 94}
95
96//_____________________________________________________
97AliEMCALTenderSupply::AliEMCALTenderSupply(const char *name, const AliTender *tender) :
6b7df55d 98AliTenderSupply(name,tender)
99,fEMCALGeo(0x0)
100,fEMCALGeoName("EMCAL_COMPLETEV1")
101,fEMCALRecoUtils(0)
102,fConfigName("")
103,fDebugLevel(0)
104,fNonLinearFunc(AliEMCALRecoUtils::kNoCorrection)
105,fNonLinearThreshold(30)
106,fReCalibCluster(kFALSE)
107,fReCalibCell(kFALSE)
108,fRecalClusPos(kFALSE)
109,fFiducial(kFALSE)
110,fNCellsFromEMCALBorder(1)
111,fRecalDistToBadChannels(kFALSE)
112,fInputTree(0)
113,fInputFile(0)
114,fFilepass(0)
115,fMass(0.139)
116,fStep(50)
117,fCutEtaPhiSum(kTRUE)
118,fCutEtaPhiSeparate(kFALSE)
119,fRcut(0.05)
120,fEtacut(0.025)
121,fPhicut(0.05)
122,fBasePath(".")
123,fReClusterize(kFALSE)
124,fClusterizer(0)
125,fGeomMatrixSet(kFALSE)
126,fLoadGeomMatrices(kFALSE)
127,fRecParam(0)
128,fOCDBpath(" ")
129,fUnfolder(0)
130,fDigitsArr(0)
131,fClusterArr(0)
6d67b5a7 132{
766cc9de 133 // Named constructor
6b7df55d 134
766cc9de 135 for(Int_t i = 0; i < 10; i++) fEMCALMatrix[i] = 0 ;
6b7df55d 136
766cc9de 137 fRecParam = new AliEMCALRecParam;
138 fEMCALRecoUtils = new AliEMCALRecoUtils();
139 fDigitsArr = new TClonesArray("AliEMCALDigit",1000);
6d67b5a7 140}
141
142//_____________________________________________________
143AliEMCALTenderSupply::~AliEMCALTenderSupply()
144{
766cc9de 145 //Destructor
6b7df55d 146
766cc9de 147 delete fEMCALRecoUtils;
148 delete fClusterizer;
149 delete fUnfolder;
150 if (fDigitsArr){
151 fDigitsArr->Clear("C");
152 delete fDigitsArr;
153 }
6d67b5a7 154}
155
156//_____________________________________________________
157void AliEMCALTenderSupply::Init()
158{
766cc9de 159 // Initialise EMCAL tender.
6b7df55d 160
766cc9de 161 if (fDebugLevel>0)
162 AliInfo("Init EMCAL Tender supply");
6b7df55d 163
766cc9de 164 if(fConfigName.Length()>0 && gROOT->LoadMacro(fConfigName) >=0) {
165 AliDebug(1, Form("Loading settings from macro %s", fConfigName.Data()));
166 AliEMCALTenderSupply *tender = (AliEMCALTenderSupply*)gInterpreter->ProcessLine("ConfigEMCALTenderSupply()");
167 fDebugLevel = tender->fDebugLevel;
168 fEMCALGeoName = tender->fEMCALGeoName;
169 delete fEMCALRecoUtils;
170 fEMCALRecoUtils = tender->fEMCALRecoUtils;
171 fConfigName = tender->fConfigName;
172 fNonLinearFunc = tender->fNonLinearFunc;
173 fNonLinearThreshold = tender->fNonLinearThreshold;
174 fReCalibCluster = tender->fReCalibCluster;
175 fReCalibCell = tender->fReCalibCell;
176 fRecalClusPos = tender->fRecalClusPos;
177 fFiducial = tender->fFiducial;
178 fNCellsFromEMCALBorder = tender->fNCellsFromEMCALBorder;
179 fRecalDistToBadChannels = tender->fRecalDistToBadChannels;
180 fMass = tender->fMass;
181 fStep = tender->fStep;
182 fRcut = tender->fRcut;
183 fEtacut = tender->fEtacut;
184 fPhicut = tender->fPhicut;
185 fReClusterize = tender->fReClusterize;
186 fLoadGeomMatrices = tender->fLoadGeomMatrices;
187 fRecParam = tender->fRecParam;
188 fOCDBpath = tender->fOCDBpath;
189 for(Int_t i = 0; i < 10; i++)
190 fEMCALMatrix[i] = tender->fEMCALMatrix[i] ;
191 }
6b7df55d 192
766cc9de 193 // Init goemetry
194 fEMCALGeo = AliEMCALGeometry::GetInstance(fEMCALGeoName) ;
6b7df55d 195
766cc9de 196 // Initialising Non linearity parameters
197 fEMCALRecoUtils->SetNonLinearityThreshold(fNonLinearThreshold);
198 fEMCALRecoUtils->SetNonLinearityFunction(fNonLinearFunc);
6b7df55d 199
766cc9de 200 // Setting mass, step size and residual cut
201 fEMCALRecoUtils->SetMass(fMass);
202 fEMCALRecoUtils->SetStep(fStep);
203 if(fCutEtaPhiSum){
204 fEMCALRecoUtils->SwitchOnCutEtaPhiSum();
205 fEMCALRecoUtils->SetCutR(fRcut);
206 } else if (fCutEtaPhiSeparate) {
207 fEMCALRecoUtils->SwitchOnCutEtaPhiSeparate();
208 fEMCALRecoUtils->SetCutEta(fEtacut);
209 fEMCALRecoUtils->SetCutPhi(fPhicut);
210 }
6d67b5a7 211}
212
213//_____________________________________________________
214void AliEMCALTenderSupply::ProcessEvent()
215{
766cc9de 216 // Event loop.
6b7df55d 217
766cc9de 218 AliESDEvent *event = fTender->GetEvent();
219 if (!event) {
220 AliError("ESD event ptr = 0, returning");
221 return;
222 }
6b7df55d 223
766cc9de 224 // Initialising parameters once per run number
225 if(fTender->RunChanged()){
226 GetPass();
227 if (!InitBadChannels()) {
228 AliError("InitBadChannels returned false, returning");
229 return;
230 }
231 if (fRecalClusPos || fReClusterize) {
232 if (!InitMisalignMatrix()) {
233 AliError("InitMisalignmentMatrix returned false, returning");
234 return;
235 }
236 }
237 if (fReCalibCluster || fReCalibCell) {
238 if (!InitRecalib()) {
239 AliError("InitRecalib returned false, returning");
240 return;
241 }
242 }
243 if (fReClusterize) {
244 if (!InitClusterization()) {
245 AliError("InitClusterization returned false, returning");
246 return;
247 }
248 }
6b7df55d 249
766cc9de 250 if(fDebugLevel>1)
251 fEMCALRecoUtils->Print("");
252 }
6b7df55d 253
766cc9de 254 // Test if cells present
255 AliESDCaloCells *cells= event->GetEMCALCells();
256 if (cells->GetNumberOfCells()<=0) {
257 AliWarning(Form("Number of EMCAL cells = %d, returning", cells->GetNumberOfCells()));
258 return;
259 }
6b7df55d 260
766cc9de 261 // Recalibrate cells
262 if (fReCalibCell)
263 RecalibrateCells();
6b7df55d 264
766cc9de 265 // Reclusterize
266 if(fReClusterize){
267 FillDigitsArray();
268 Clusterize();
269 UpdateClusters();
270 }
6b7df55d 271
766cc9de 272 // Store good clusters
273 TClonesArray *clusArr = dynamic_cast<TClonesArray*>(event->FindListObject("caloClusters"));
274 if (!clusArr)
275 clusArr = dynamic_cast<TClonesArray*>(event->FindListObject("CaloClusters"));
276 if (!clusArr) {
277 AliWarning(Form("No cluster array, number of cells in event = %d, returning", cells->GetNumberOfCells()));
278 return;
279 }
280 Int_t nclusters = clusArr->GetEntriesFast();
281 for (Int_t icluster=0; icluster < nclusters; ++icluster) {
282 AliVCluster *clust = static_cast<AliVCluster*>(clusArr->At(icluster));
283 if (!clust)
284 continue;
285 if (!clust->IsEMCAL())
286 continue;
6b7df55d 287
766cc9de 288 if (fEMCALRecoUtils->ClusterContainsBadChannel(fEMCALGeo, clust->GetCellsAbsId(), clust->GetNCells())) {
289 delete clusArr->RemoveAt(icluster);
290 continue; //todo is it really needed to remove it? Or should we flag it?
291 }
6b7df55d 292
766cc9de 293 if (fFiducial){
294 if (!fEMCALRecoUtils->CheckCellFiducialRegion(fEMCALGeo, clust, cells)){
295 delete clusArr->RemoveAt(icluster);
296 continue; // todo it would be nice to store the distance
297 }
298 }
6b7df55d 299
766cc9de 300 fEMCALRecoUtils->CorrectClusterEnergyLinearity(clust);
301 if(fRecalDistToBadChannels)
302 fEMCALRecoUtils->RecalculateClusterDistanceToBadChannel(fEMCALGeo, cells, clust);
303 if(fReCalibCluster)
304 fEMCALRecoUtils->RecalibrateClusterEnergy(fEMCALGeo, clust, cells);
305 if(fRecalClusPos)
306 fEMCALRecoUtils->RecalculateClusterPosition(fEMCALGeo, cells, clust);
307 }
308 clusArr->Compress();
6b7df55d 309
766cc9de 310 // Track matching
311 if (!TGeoGlobalMagField::Instance()->GetField()) {
312 const AliESDRun *erun = event->GetESDRun();
313 AliMagF *field = AliMagF::CreateFieldMap(erun->GetCurrentL3(),
314 erun->GetCurrentDip(),
315 AliMagF::kConvLHC,
316 kFALSE,
317 erun->GetBeamEnergy(),
318 erun->GetBeamType());
319 TGeoGlobalMagField::Instance()->SetField(field);
320 }
321 fEMCALRecoUtils->FindMatches(event,0x0,fEMCALGeo);
322 Int_t nTracks = event->GetNumberOfTracks();
323 if (nTracks>0) {
324 SetClusterMatchedToTrack(event);
325 SetTracksMatchedToCluster(event);
326 }
6d67b5a7 327}
328
329//_____________________________________________________
330void AliEMCALTenderSupply::SetClusterMatchedToTrack(AliESDEvent *event)
331{
766cc9de 332 // Checks if tracks are matched to EMC clusters and set the matched EMCAL cluster index to ESD track.
6b7df55d 333
766cc9de 334 Int_t nTracks = event->GetNumberOfTracks();
335 for (Int_t iTrack = 0; iTrack < nTracks; ++iTrack) {
336 AliESDtrack* track = event->GetTrack(iTrack);
337 if (!track) {
338 AliWarning(Form("Could not receive track %d", iTrack));
339 continue;
340 }
341 Int_t matchClusIndex = fEMCALRecoUtils->GetMatchedClusterIndex(iTrack);
342 track->SetEMCALcluster(matchClusIndex); //sets -1 if track not matched within residual
343 if(matchClusIndex != -1)
344 track->SetStatus(AliESDtrack::kEMCALmatch);
345 }
346 if (fDebugLevel>2)
347 AliInfo("Track matched to closest cluster");
6d67b5a7 348}
349
350//_____________________________________________________
351void AliEMCALTenderSupply::SetTracksMatchedToCluster(AliESDEvent *event)
352{
766cc9de 353 // Checks if EMC clusters are matched to ESD track.
354 // Adds track indexes of all the tracks matched to a cluster withing residuals in ESDCalocluster.
6b7df55d 355
766cc9de 356 for (Int_t iClus=0; iClus < event->GetNumberOfCaloClusters(); ++iClus) {
357 AliESDCaloCluster *cluster = event->GetCaloCluster(iClus);
358 if (!cluster->IsEMCAL())
359 continue;
6b7df55d 360
766cc9de 361 Int_t nTracks = event->GetNumberOfTracks();
362 TArrayI arrayTrackMatched(nTracks);
6b7df55d 363
766cc9de 364 // Get the closest track matched to the cluster
365 Int_t nMatched = 0;
366 Int_t matchTrackIndex = fEMCALRecoUtils->GetMatchedTrackIndex(iClus);
367 if (matchTrackIndex != -1) {
368 arrayTrackMatched[nMatched] = matchTrackIndex;
369 nMatched++;
370 }
6b7df55d 371
766cc9de 372 // Get all other tracks matched to the cluster
373 for(Int_t iTrk=0; iTrk<nTracks; ++iTrk) {
374 AliESDtrack* track = event->GetTrack(iTrk);
375 if(iTrk == matchTrackIndex) continue;
376 if(track->GetEMCALcluster() == iClus){
377 arrayTrackMatched[nMatched] = iTrk;
378 ++nMatched;
379 }
380 }
6b7df55d 381
766cc9de 382 arrayTrackMatched.Set(nMatched);
383 cluster->AddTracksMatched(arrayTrackMatched);
384
385 Float_t eta= -999, phi = -999;
386 if (matchTrackIndex != -1)
387 fEMCALRecoUtils->GetMatchedResiduals(iClus, eta, phi);
388 cluster->SetTrackDistance(phi, eta);
389 }
6b7df55d 390
766cc9de 391 if (fDebugLevel>2)
392 AliInfo("Cluster matched to tracks");
6d67b5a7 393}
394
395//_____________________________________________________
396Bool_t AliEMCALTenderSupply::InitMisalignMatrix()
397{
766cc9de 398 // Initialising misalignment matrices
6b7df55d 399
766cc9de 400 AliESDEvent *event = fTender->GetEvent();
401 if (!event)
402 return kFALSE;
6b7df55d 403
766cc9de 404 if (fGeomMatrixSet) {
405 AliInfo("Misalignment matrix already set");
406 return kTRUE;
407 }
6b7df55d 408
766cc9de 409 if (fDebugLevel>0)
410 AliInfo("Initialising misalignment matrix");
6b7df55d 411
766cc9de 412 fEMCALRecoUtils->SetPositionAlgorithm(AliEMCALRecoUtils::kPosTowerGlobal);
6b7df55d 413
766cc9de 414 if (fLoadGeomMatrices) {
415 for(Int_t mod=0; mod < fEMCALGeo->GetNumberOfSuperModules(); ++mod) {
416 if (fEMCALMatrix[mod]){
417 if(fDebugLevel > 2)
418 fEMCALMatrix[mod]->Print();
419 fEMCALGeo->SetMisalMatrix(fEMCALMatrix[mod],mod);
420 }
421 }
422 fGeomMatrixSet = kTRUE;
423 return kTRUE;
424 }
6b7df55d 425
766cc9de 426 Int_t runGM = event->GetRunNumber();
6b7df55d 427 TObjArray *mobj = 0;
428
766cc9de 429 if (runGM <= 140000) { //2010 data
6b7df55d 430 AliOADBContainer emcalgeoCont(Form("emcal2010"));
431 emcalgeoCont.InitFromFile("$ALICE_ROOT/OADB/EMCAL/EMCALlocal2master.root",Form("AliEMCALgeo"));
432 mobj=(TObjArray*)emcalgeoCont.GetObject(100,"survey10");
433
766cc9de 434 } else if (runGM>140000 && runGM <148531 && (fFilepass == "pass1")) { // 2011 LHC11a pass1 data
435 AliOADBContainer emcalgeoCont(Form("emcal2011"));
6b7df55d 436 // emcalgeoCont.InitFromFile(Form("%s/BetaGood.root",fBasePath.Data()),Form("AliEMCALgeo"));
437 emcalgeoCont.InitFromFile("$ALICE_ROOT/OADB/EMCAL/EMCALlocal2master.root",Form("AliEMCALgeo"));
438 mobj=(TObjArray*)emcalgeoCont.GetObject(100,"survey11byS"); }
439
440 if((runGM<=140000)||(runGM>140000 && runGM <148531 && (fFilepass == "pass1"))){
766cc9de 441 for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
442 fEMCALMatrix[mod] = (TGeoHMatrix*) mobj->At(mod);
443 fEMCALGeo->SetMisalMatrix(fEMCALMatrix[mod],mod);
444 fEMCALMatrix[mod]->Print();
445 }
446 } else {
447 AliInfo(Form("No external misalignment matrix set: %d - %s, loading from ESD event", runGM, fFilepass.Data()));
448 for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
449 if(fDebugLevel > 2)
450 event->GetEMCALMatrix(mod)->Print();
451 if(event->GetEMCALMatrix(mod))
452 fEMCALGeo->SetMisalMatrix(event->GetEMCALMatrix(mod),mod) ;
453 }
454 fGeomMatrixSet=kTRUE;
455 }
456 return kTRUE;
6d67b5a7 457}
458
459//_____________________________________________________
460Bool_t AliEMCALTenderSupply::InitBadChannels()
461{
766cc9de 462 // Initialising bad channel maps
6b7df55d 463
766cc9de 464 AliESDEvent *event = fTender->GetEvent();
465 if (!event)
466 return kFALSE;
6b7df55d 467
766cc9de 468 if (fDebugLevel>0)
469 AliInfo("Initialising Bad channel map");
6b7df55d 470
766cc9de 471 if (fFiducial){
472 fEMCALRecoUtils->SetNumberOfCellsFromEMCALBorder(fNCellsFromEMCALBorder);
473 fEMCALRecoUtils->SwitchOnNoFiducialBorderInEMCALEta0();
474 }
6b7df55d 475
766cc9de 476 fEMCALRecoUtils->SwitchOnBadChannelsRemoval();
477 if (fRecalDistToBadChannels)
478 fEMCALRecoUtils->SwitchOnDistToBadChannelRecalculation();
6b7df55d 479
766cc9de 480 Int_t runBC = event->GetRunNumber();
481 if (runBC <=140000) { //2010
482 TFile *fbad = new TFile(Form("%s/BadChannels.root",fBasePath.Data()),"read");
483 if (fbad->IsZombie()) {
484 TString fPath = TString(fBasePath.Data());
485 if (fPath.Contains("alien")) {
486 if (fDebugLevel>1)
487 AliInfo("Connecting to alien to get BadChannels.root");
488 delete fbad;
489 fbad = TFile::Open(Form("%s/BadChannelsDB/BadChannels.root",fBasePath.Data()));
490 }
491 }
492 for (Int_t i=0; i<4; ++i) {
493 TH2I *h = fEMCALRecoUtils->GetEMCALChannelStatusMap(i);
494 if (h)
495 delete h;
496 h = (TH2I*)fbad->Get(Form("EMCALBadChannelMap_Mod%d",i));
497 if (!h) {
498 AliError(Form("Can not get EMCALBadChannelMap_Mod%d",i));
499 continue;
500 }
501 h->SetDirectory(0);
502 fEMCALRecoUtils->SetEMCALChannelStatusMap(i,h);
503 }
504 delete fbad;
505 return kTRUE;
506 }
6b7df55d 507
766cc9de 508 //2011
509 Int_t nSupMod=-1, nModule=-1, nIphi=-1, nIeta=-1, iphi=-1, ieta=-1;
510 if (runBC>=144871 && runBC<=146860){ // LHC11a 2.76 TeV pp
6b7df55d 511
766cc9de 512 if (fFilepass == "pass1") { // pass1
513 const Int_t nTowers=89;
514 Int_t hotChannels[nTowers]={74, 103, 152, 320, 321, 322, 323, 324, 325, 326, 327, 328, 329, 330, 331, 332, 333, 334, 335, 368, 369, 370, 371, 372, 373, 374,375, 376, 377, 378, 379, 380, 381, 382, 383, 917, 1275, 1288, 1519, 1595, 1860, 1967, 2022, 2026, 2047, 2117, 2298, 2540, 2776, 3135, 3764, 6095, 6111, 6481, 6592, 6800, 6801, 6802, 6803, 6804, 6805, 6806, 6807, 6808, 6809, 6810, 6811, 6812, 6813, 6814, 6815, 7371, 7425, 7430, 7457, 7491, 7709, 8352, 8353, 8356, 8357, 8808, 8810, 8812, 8814, 9056, 9769, 9815, 9837};
515 for (Int_t i=0; i<nTowers; ++i) {
516 fEMCALGeo->GetCellIndex(hotChannels[i],nSupMod,nModule,nIphi,nIeta);
517 fEMCALGeo->GetCellPhiEtaIndexInSModule(nSupMod,nModule,nIphi,nIeta,iphi,ieta);
518 fEMCALRecoUtils->SetEMCALChannelStatus(nSupMod, ieta, iphi);
519 }
520 } else if (fFilepass == "pass2") { // pass2
521 const Int_t nTowers=24;
522 Int_t hotChannels[nTowers]= {74, 103, 152, 917, 1059, 1175, 1276, 1288, 1376, 1382, 1595, 2022, 2026, 2210, 2540, 2778, 2793, 3135, 3764, 5767, 6481, 7371, 7878, 9769};
523 for(Int_t i=0; i<nTowers; ++i) {
524 fEMCALGeo->GetCellIndex(hotChannels[i],nSupMod,nModule,nIphi,nIeta);
525 fEMCALGeo->GetCellPhiEtaIndexInSModule(nSupMod,nModule,nIphi,nIeta,iphi,ieta);
526 fEMCALRecoUtils->SetEMCALChannelStatus(nSupMod, ieta, iphi);
527 }
528 }
529 return kTRUE;
530 }
6b7df55d 531
532 if (runBC>=151636 && runBC<=155384) { //LHC11c
533 const Int_t nTowers=231;
534 Int_t hotChannels[nTowers]={74, 103, 152, 320, 321, 322, 323, 324, 325, 326, 327, 328, 329, 330, 331, 332, 333, 334, 335, 368, 369, 370, 371, 372, 373, 374,375, 376, 377, 378, 379, 380, 381, 382, 383, 917, 1059, 1160, 1263, 1275, 1276, 1288, 1376, 1382, 1384, 1414,1519, 1595, 1860, 1720, 1912, 1961, 1967, 2016, 2017, 2019, 2022, 2023, 2026, 2027, 2047, 2065, 2071, 2074, 2079, 2112, 2114, 2115, 2116, 2117, 2120, 2123, 2145, 2160, 2190, 2298, 2506, 2540, 2776, 2778, 3135, 3544, 3567, 3664, 3665, 3666, 3667, 3668, 3669, 3670, 3671, 3672, 3673, 3674, 3675, 3676, 3677, 3678, 3679, 3712, 3713, 3714, 3715, 3716, 3717, 3718, 3719, 3720, 3721, 3722, 3723, 3724, 3725, 3726, 3727, 3764, 4026, 4121, 4157, 4208, 4209, 4568, 5560, 5767, 5969, 6065, 6076, 6084, 6087, 6095, 6111, 6128, 6129, 6130, 6131, 6133, 6136, 6137, 6138, 6139, 6140, 6141, 6142, 6143, 6340, 6369, 6425, 6481, 6561, 6592, 6678, 6907, 6925, 6800, 6801, 6802, 6803, 6804, 6805, 6806, 6807, 6808, 6809, 6810, 6811, 6812, 6813, 6814, 6815, 7089, 7371, 7375, 7425, 7457, 7430, 7491, 7709, 7874, 8273, 8352, 8353, 8354, 8356, 8357, 8362, 8808, 8769, 8810, 8811, 8812, 8814, 9056, 9217, 9302, 9365, 9389, 9815, 9769, 9833, 9837, 9850, 9895, 9997, 10082, 10086, 10087, 10091, 10095, 10112, 10113, 10114, 10115, 10116, 10117, 10118, 10119, 10120, 10121, 10122, 10123, 10576, 10718, 10723, 10918, 10919, 10922, 10925, 10926, 10927, 11276, 1136};
766cc9de 535 for(Int_t i=0; i<nTowers; i++) {
536 fEMCALGeo->GetCellIndex(hotChannels[i],nSupMod,nModule,nIphi,nIeta);
537 fEMCALGeo->GetCellPhiEtaIndexInSModule(nSupMod,nModule,nIphi,nIeta,iphi,ieta);
538 fEMCALRecoUtils->SetEMCALChannelStatus(nSupMod, ieta, iphi);
539 }
540 return kTRUE;
541 }
6b7df55d 542
766cc9de 543 AliInfo(Form("No external hot channel set: %d - %s", runBC, fFilepass.Data()));
544 return kTRUE;
6d67b5a7 545}
546
547//_____________________________________________________
ab73393f 548Bool_t AliEMCALTenderSupply::InitRecalib()
6d67b5a7 549{
766cc9de 550 // Initialising recalibration factors.
6b7df55d 551
766cc9de 552 AliESDEvent *event = fTender->GetEvent();
553 if (!event)
554 return kFALSE;
6b7df55d 555
766cc9de 556 if (fDebugLevel>0)
557 AliInfo("Initialising recalibration factors");
6b7df55d 558
766cc9de 559 fEMCALRecoUtils->SwitchOnRecalibration();
560 Int_t runRC = event->GetRunNumber();
6b7df55d 561
766cc9de 562 TFile *fRecalib = 0;
563 if (runRC <= 140000) { // 2010
564 fRecalib = new TFile(Form("%s/RecalibrationFactors.root",fBasePath.Data()),"read");
565 if (fRecalib->IsZombie()) {
566 TString fPath = TString(fBasePath.Data());
567 if(fPath.Contains("alien")) {
568 if (fDebugLevel>1)
569 AliInfo("Connecting to alien to get RecalibrationFactors.root");
6b7df55d 570
766cc9de 571 if( (runRC >= 114737 && runRC <= 117223 && (fFilepass = "pass2")) || //LHC10b pass2
6b7df55d 572 (runRC >= 118503 && runRC <= 121040 && (fFilepass = "pass2")) || //LHC10c pass2
573 (runRC >= 118503 && runRC <= 121040 && (fFilepass = "pass3")) || //LHC10c pass3
574 (runRC >= 122195 && runRC <= 126437 && (fFilepass = "pass1")) || //LHC10d pass1
575 (runRC >= 127712 && runRC <= 130850 && (fFilepass = "pass1")) || //LHC10e pass1
576 (runRC >= 133004 && runRC < 134657 && (fFilepass = "pass1"))) //LHC10f pass1<134657
766cc9de 577 {
578 fRecalib = TFile::Open(Form("%s/RecalDB/summer_december_2010/RecalibrationFactors.root",fBasePath.Data()));
579 }
580 else if((runRC >= 122195 && runRC <= 126437 && (fFilepass = "pass2")) || //LHC10d pass2
581 (runRC >= 134657 && runRC <= 135029 && (fFilepass = "pass1")) || //LHC10f pass1 >= 134657
582 (runRC >= 135654 && runRC <= 136377 && (fFilepass = "pass1")) || //LHC10g pass1
583 (runRC >= 136851 && runRC < 137231 && (fFilepass = "pass1"))) //LHC10h pass1 untill christmas
584 {
585 fRecalib = TFile::Open(Form("%s/RecalDB/december2010/RecalibrationFactors.root",fBasePath.Data()));
586 }
587 else if(runRC >= 137231 && runRC <= 139517 && (fFilepass = "pass1")) //LHC10h pass1 from christmas
588 {
589 fRecalib = TFile::Open(Form("%s/RecalDB/summer2010/RecalibrationFactors.root",fBasePath.Data()));
590 }
591 else {
592 AliError(Form("Run number or pass number not found: %d - %s; RECALIBRATION NOT APPLIED", runRC, fFilepass.Data()));
593 return kTRUE;
594 }
595 }
596 }
597 } else if (runRC > 140000) { // 2011
598 fRecalib = new TFile(Form("%s/RecalibrationFactors.root",fBasePath.Data()),"read");
599 if (fRecalib->IsZombie()){
600 TString fPath = TString(fBasePath.Data());
601 if(fPath.Contains("alien")) {
602 if (fDebugLevel>1)
603 AliInfo("Connecting to alien to get RecalibrationFactors.root");
604 fRecalib = TFile::Open(Form("%s/RecalDB/2011_v0/RecalibrationFactors.root",fBasePath.Data()));
605 if (!fRecalib)
606 AliError("Recalibration file not found");
607 return kFALSE;
608 }
609 }
6b7df55d 610
766cc9de 611 Int_t sms = fEMCALGeo->GetEMCGeometry()->GetNumberOfSuperModules();
612 for (Int_t i=0; i<sms; ++i) {
613 TH2F *h = fEMCALRecoUtils->GetEMCALChannelRecalibrationFactors(i);
614 if (h)
615 delete h;
616 h = (TH2F*)fRecalib->Get(Form("EMCALRecalFactors_SM%d",i));
617 if (!h) {
618 AliError(Form("Could not load EMCALRecalFactors_SM%d",i));
619 continue;
620 }
621 h->SetDirectory(0);
622 fEMCALRecoUtils->SetEMCALChannelRecalibrationFactors(i,h);
623 }
624 }
625 delete fRecalib;
626 return kTRUE;
c958a2f7 627}
628
629//_____________________________________________________
630void AliEMCALTenderSupply::RecalibrateCells()
631{
766cc9de 632 // Recalibrate cells.
6b7df55d 633
766cc9de 634 AliESDCaloCells *cells = fTender->GetEvent()->GetEMCALCells();
6b7df55d 635
766cc9de 636 Int_t nEMCcell = cells->GetNumberOfCells();
637 for(Int_t icell=0; icell<nEMCcell; ++icell) {
638 Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
639 fEMCALGeo->GetCellIndex(cells->GetCellNumber(icell),imod,iTower,iIphi,iIeta);
640 fEMCALGeo->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
641 Double_t calibFactor = fEMCALRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi);
642 cells->SetCell(icell,cells->GetCellNumber(icell),cells->GetAmplitude(icell)*calibFactor,cells->GetTime(icell));
643 }
b20bc239 644}
645
646//_____________________________________________________
766cc9de 647Bool_t AliEMCALTenderSupply::InitClusterization()
b20bc239 648{
766cc9de 649 // Initialing clusterization/unfolding algorithm and set all the needed parameters.
6b7df55d 650
766cc9de 651 AliESDEvent *event=fTender->GetEvent();
652 if (!event)
653 return kFALSE;
6b7df55d 654
766cc9de 655 if (fDebugLevel>0)
656 AliInfo(Form("Initialising reclustering parameters: Clusterizer-%d",fRecParam->GetClusterizerFlag()));
6b7df55d 657
766cc9de 658 //---setup clusterizer
659 delete fClusterizer;
660 if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv1)
661 fClusterizer = new AliEMCALClusterizerv1 (fEMCALGeo);
662 else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv2)
663 fClusterizer = new AliEMCALClusterizerv2(fEMCALGeo);
664 else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerNxN) {
665 AliEMCALClusterizerNxN *clusterizer = new AliEMCALClusterizerNxN(fEMCALGeo);
666 clusterizer->SetNRowDiff(fRecParam->GetNRowDiff());
667 clusterizer->SetNColDiff(fRecParam->GetNColDiff());
668 fClusterizer = clusterizer;
669 } else {
670 AliFatal(Form("Clusterizer < %d > not available", fRecParam->GetClusterizerFlag()));
671 return kFALSE;
672 }
6b7df55d 673
766cc9de 674 // Set the clustering parameters
675 fClusterizer->SetECAClusteringThreshold( fRecParam->GetClusteringThreshold() );
676 fClusterizer->SetECALogWeight ( fRecParam->GetW0() );
677 fClusterizer->SetMinECut ( fRecParam->GetMinECut() );
678 fClusterizer->SetUnfolding ( fRecParam->GetUnfold() );
679 fClusterizer->SetECALocalMaxCut ( fRecParam->GetLocMaxCut() );
680 fClusterizer->SetTimeCut ( fRecParam->GetTimeCut() );
681 fClusterizer->SetTimeMin ( fRecParam->GetTimeMin() );
682 fClusterizer->SetTimeMax ( fRecParam->GetTimeMax() );
683 fClusterizer->SetInputCalibrated ( kTRUE );
684 fClusterizer->SetJustClusters ( kTRUE );
6b7df55d 685
766cc9de 686 // In case of unfolding after clusterization is requested, set the corresponding parameters
687 if (fRecParam->GetUnfold()) {
688 for (Int_t i = 0; i < 8; ++i) {
689 fClusterizer->SetSSPars(i, fRecParam->GetSSPars(i));
690 }
691 for (Int_t i = 0; i < 3; ++i) {
692 fClusterizer->SetPar5 (i, fRecParam->GetPar5(i));
693 fClusterizer->SetPar6 (i, fRecParam->GetPar6(i));
694 }
695 fClusterizer->InitClusterUnfolding();
696 }
6b7df55d 697
766cc9de 698 fClusterizer->SetDigitsArr(fDigitsArr);
699 fClusterizer->SetOutput(0);
700 fClusterArr = const_cast<TObjArray *>(fClusterizer->GetRecPoints());
701 return kTRUE;
b20bc239 702}
766cc9de 703
b20bc239 704//_____________________________________________________
705void AliEMCALTenderSupply::FillDigitsArray()
706{
766cc9de 707 // Fill digits from cells to a TClonesArray.
6b7df55d 708
766cc9de 709 AliESDEvent *event = fTender->GetEvent();
710 if (!event)
711 return;
6b7df55d 712
766cc9de 713 fDigitsArr->Clear("C");
714 AliESDCaloCells *cells = event->GetEMCALCells();
715 Int_t ncells = cells->GetNumberOfCells();
716 for (Int_t icell = 0, idigit = 0; icell < ncells; ++icell) {
717 Double_t cellAmplitude=0, cellTime=0;
718 Short_t cellNumber=0;
719 if (cells->GetCell(icell, cellNumber, cellAmplitude, cellTime) != kTRUE)
720 break;
721 AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->New(idigit));
722 digit->SetId(cellNumber);
723 digit->SetTime(cellTime);
724 digit->SetTimeR(cellTime);
725 digit->SetIndexInList(idigit);
726 digit->SetType(AliEMCALDigit::kHG);
727 digit->SetAmplitude(cellAmplitude);
728 idigit++;
729 }
b20bc239 730}
731
732//_____________________________________________________
733void AliEMCALTenderSupply::Clusterize()
734{
766cc9de 735 // Clusterize.
6b7df55d 736
766cc9de 737 fClusterizer->Digits2Clusters("");
b20bc239 738}
739
740//_____________________________________________________
741void AliEMCALTenderSupply::UpdateClusters()
742{
766cc9de 743 // Update ESD cluster list.
6b7df55d 744
766cc9de 745 AliESDEvent *event = fTender->GetEvent();
746 if (!event)
747 return;
6b7df55d 748
766cc9de 749 TClonesArray *clus = dynamic_cast<TClonesArray*>(event->FindListObject("caloClusters"));
750 if (!clus)
751 clus = dynamic_cast<TClonesArray*>(event->FindListObject("CaloClusters"));
752 if (!clus) {
6b7df55d 753 AliError(" Null pointer to calo clusters array, returning");
766cc9de 754 return;
755 }
6b7df55d 756
766cc9de 757 Int_t nents = clus->GetEntriesFast();
758 for (Int_t i=0; i < nents; ++i) {
759 AliESDCaloCluster *c = static_cast<AliESDCaloCluster*>(clus->At(i));
760 if (!c)
761 continue;
762 if (c->IsEMCAL())
763 delete clus->RemoveAt(i);
764 }
765 clus->Compress();
766 RecPoints2Clusters(clus);
b20bc239 767}
768
769//_____________________________________________________
770void AliEMCALTenderSupply::RecPoints2Clusters(TClonesArray *clus)
771{
766cc9de 772 // Convert AliEMCALRecoPoints to AliESDCaloClusters.
773 // Cluster energy, global position, cells and their amplitude fractions are restored.
774
775 AliESDEvent *event = fTender->GetEvent();
776 if (!event)
777 return;
6b7df55d 778
766cc9de 779 Int_t ncls = fClusterArr->GetEntriesFast();
780 for(Int_t i=0, nout=clus->GetEntriesFast(); i < ncls; ++i) {
781 AliEMCALRecPoint *recpoint = static_cast<AliEMCALRecPoint*>(fClusterArr->At(i));
782
783 Int_t ncells_true = 0;
784 const Int_t ncells = recpoint->GetMultiplicity();
785 UShort_t absIds[ncells];
786 Double32_t ratios[ncells];
787 Int_t *dlist = recpoint->GetDigitsList();
788 Float_t *elist = recpoint->GetEnergiesList();
789 for (Int_t c = 0; c < ncells; ++c) {
790 AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(dlist[c]));
791 absIds[ncells_true] = digit->GetId();
792 ratios[ncells_true] = elist[c]/digit->GetAmplitude();
793 if (ratios[ncells_true] < 0.001)
794 continue;
795 ++ncells_true;
796 }
797
798 if (ncells_true < 1) {
799 AliWarning("Skipping cluster with no cells");
800 continue;
801 }
802
803 // calculate new cluster position
804 TVector3 gpos;
805 recpoint->GetGlobalPosition(gpos);
806 Float_t g[3];
807 gpos.GetXYZ(g);
808
809 AliESDCaloCluster *c = static_cast<AliESDCaloCluster*>(clus->New(nout++));
810 c->SetType(AliVCluster::kEMCALClusterv1);
811 c->SetE(recpoint->GetEnergy());
812 c->SetPosition(g);
813 c->SetNCells(ncells_true);
814 c->SetDispersion(recpoint->GetDispersion());
815 c->SetEmcCpvDistance(-1); //not yet implemented
816 c->SetChi2(-1); //not yet implemented
817 c->SetTOF(recpoint->GetTime()) ; //time-of-flight
818 c->SetNExMax(recpoint->GetNExMax()); //number of local maxima
819 Float_t elipAxis[2];
820 recpoint->GetElipsAxis(elipAxis);
821 c->SetM02(elipAxis[0]*elipAxis[0]) ;
822 c->SetM20(elipAxis[1]*elipAxis[1]) ;
823 AliESDCaloCluster *cesd = static_cast<AliESDCaloCluster*>(c);
824 cesd->SetCellsAbsId(absIds);
825 cesd->SetCellsAmplitudeFraction(ratios);
826 }
827}
b20bc239 828
766cc9de 829//_____________________________________________________
830void AliEMCALTenderSupply::GetPass()
831{
832 // Get passx from filename.
6b7df55d 833
766cc9de 834 AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
835 fInputTree = mgr->GetTree();
6b7df55d 836
766cc9de 837 if (!fInputTree) {
6b7df55d 838 AliError("Pointer to tree = 0, returning");
766cc9de 839 return;
840 }
6b7df55d 841
766cc9de 842 fInputFile = fInputTree->GetCurrentFile();
843 if (!fInputFile) {
6b7df55d 844 AliError("Null pointer input file, returning");
766cc9de 845 return;
846 }
6b7df55d 847
766cc9de 848 TString fname(fInputFile->GetName());
849 if (fname.Contains("pass1")) fFilepass = TString("pass1");
850 else if(fname.Contains("pass2")) fFilepass = TString("pass2");
851 else if(fname.Contains("pass3")) fFilepass = TString("pass3");
852 else {
853 AliError(Form("Pass number string not found: %s", fname.Data()));
854 return;
855 }
b20bc239 856}
857