]> git.uio.no Git - u/mrichter/AliRoot.git/blame - ANALYSIS/TenderSupplies/AliEMCALTenderSupply.cxx
Recover recalibration and bad channels map from OADB when needed (Marcel, Deepa,...
[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)
acf53135 67,fUpdateCell(kFALSE)
6b7df55d 68,fRecalClusPos(kFALSE)
69,fFiducial(kFALSE)
70,fNCellsFromEMCALBorder(1)
71,fRecalDistToBadChannels(kFALSE)
72,fInputTree(0)
73,fInputFile(0)
74,fFilepass(0)
75,fMass(0.139)
76,fStep(50)
77,fCutEtaPhiSum(kTRUE)
78,fCutEtaPhiSeparate(kFALSE)
79,fRcut(0.05)
80,fEtacut(0.025)
81,fPhicut(0.05)
acf53135 82,fBasePath("")
6b7df55d 83,fReClusterize(kFALSE)
84,fClusterizer(0)
85,fGeomMatrixSet(kFALSE)
86,fLoadGeomMatrices(kFALSE)
acf53135 87,fRecParam(0x0)
88,fRecParamSet(kFALSE)
6b7df55d 89,fOCDBpath(" ")
90,fUnfolder(0)
91,fDigitsArr(0)
92,fClusterArr(0)
acf53135 93,fMisalignSurvey(kdefault)
6d67b5a7 94{
766cc9de 95 // Default constructor.
6b7df55d 96 for(Int_t i = 0; i < 10; i++) fEMCALMatrix[i] = 0 ;
6d67b5a7 97}
98
99//_____________________________________________________
100AliEMCALTenderSupply::AliEMCALTenderSupply(const char *name, const AliTender *tender) :
6b7df55d 101AliTenderSupply(name,tender)
102,fEMCALGeo(0x0)
103,fEMCALGeoName("EMCAL_COMPLETEV1")
104,fEMCALRecoUtils(0)
105,fConfigName("")
106,fDebugLevel(0)
107,fNonLinearFunc(AliEMCALRecoUtils::kNoCorrection)
108,fNonLinearThreshold(30)
109,fReCalibCluster(kFALSE)
110,fReCalibCell(kFALSE)
acf53135 111,fUpdateCell(kFALSE)
6b7df55d 112,fRecalClusPos(kFALSE)
113,fFiducial(kFALSE)
114,fNCellsFromEMCALBorder(1)
115,fRecalDistToBadChannels(kFALSE)
116,fInputTree(0)
117,fInputFile(0)
118,fFilepass(0)
119,fMass(0.139)
120,fStep(50)
121,fCutEtaPhiSum(kTRUE)
122,fCutEtaPhiSeparate(kFALSE)
123,fRcut(0.05)
124,fEtacut(0.025)
125,fPhicut(0.05)
acf53135 126,fBasePath("")
6b7df55d 127,fReClusterize(kFALSE)
128,fClusterizer(0)
129,fGeomMatrixSet(kFALSE)
130,fLoadGeomMatrices(kFALSE)
acf53135 131,fRecParam(0x0)
132,fRecParamSet(kFALSE)
6b7df55d 133,fOCDBpath(" ")
134,fUnfolder(0)
135,fDigitsArr(0)
136,fClusterArr(0)
acf53135 137,fMisalignSurvey(kdefault)
6d67b5a7 138{
766cc9de 139 // Named constructor
6b7df55d 140
766cc9de 141 for(Int_t i = 0; i < 10; i++) fEMCALMatrix[i] = 0 ;
6b7df55d 142
766cc9de 143 fEMCALRecoUtils = new AliEMCALRecoUtils();
acf53135 144 fRecParam = new AliEMCALRecParam;
766cc9de 145 fDigitsArr = new TClonesArray("AliEMCALDigit",1000);
6d67b5a7 146}
147
148//_____________________________________________________
149AliEMCALTenderSupply::~AliEMCALTenderSupply()
150{
766cc9de 151 //Destructor
6b7df55d 152
766cc9de 153 delete fEMCALRecoUtils;
154 delete fClusterizer;
155 delete fUnfolder;
156 if (fDigitsArr){
157 fDigitsArr->Clear("C");
158 delete fDigitsArr;
159 }
6d67b5a7 160}
161
162//_____________________________________________________
163void AliEMCALTenderSupply::Init()
164{
766cc9de 165 // Initialise EMCAL tender.
acf53135 166
766cc9de 167 if (fDebugLevel>0)
168 AliInfo("Init EMCAL Tender supply");
6b7df55d 169
766cc9de 170 if(fConfigName.Length()>0 && gROOT->LoadMacro(fConfigName) >=0) {
171 AliDebug(1, Form("Loading settings from macro %s", fConfigName.Data()));
172 AliEMCALTenderSupply *tender = (AliEMCALTenderSupply*)gInterpreter->ProcessLine("ConfigEMCALTenderSupply()");
173 fDebugLevel = tender->fDebugLevel;
174 fEMCALGeoName = tender->fEMCALGeoName;
175 delete fEMCALRecoUtils;
176 fEMCALRecoUtils = tender->fEMCALRecoUtils;
177 fConfigName = tender->fConfigName;
178 fNonLinearFunc = tender->fNonLinearFunc;
179 fNonLinearThreshold = tender->fNonLinearThreshold;
180 fReCalibCluster = tender->fReCalibCluster;
181 fReCalibCell = tender->fReCalibCell;
182 fRecalClusPos = tender->fRecalClusPos;
183 fFiducial = tender->fFiducial;
184 fNCellsFromEMCALBorder = tender->fNCellsFromEMCALBorder;
185 fRecalDistToBadChannels = tender->fRecalDistToBadChannels;
186 fMass = tender->fMass;
187 fStep = tender->fStep;
188 fRcut = tender->fRcut;
189 fEtacut = tender->fEtacut;
190 fPhicut = tender->fPhicut;
191 fReClusterize = tender->fReClusterize;
192 fLoadGeomMatrices = tender->fLoadGeomMatrices;
193 fRecParam = tender->fRecParam;
194 fOCDBpath = tender->fOCDBpath;
195 for(Int_t i = 0; i < 10; i++)
196 fEMCALMatrix[i] = tender->fEMCALMatrix[i] ;
197 }
6b7df55d 198
766cc9de 199 // Init goemetry
200 fEMCALGeo = AliEMCALGeometry::GetInstance(fEMCALGeoName) ;
6b7df55d 201
766cc9de 202 // Initialising Non linearity parameters
203 fEMCALRecoUtils->SetNonLinearityThreshold(fNonLinearThreshold);
204 fEMCALRecoUtils->SetNonLinearityFunction(fNonLinearFunc);
6b7df55d 205
766cc9de 206 // Setting mass, step size and residual cut
207 fEMCALRecoUtils->SetMass(fMass);
208 fEMCALRecoUtils->SetStep(fStep);
209 if(fCutEtaPhiSum){
210 fEMCALRecoUtils->SwitchOnCutEtaPhiSum();
211 fEMCALRecoUtils->SetCutR(fRcut);
212 } else if (fCutEtaPhiSeparate) {
213 fEMCALRecoUtils->SwitchOnCutEtaPhiSeparate();
214 fEMCALRecoUtils->SetCutEta(fEtacut);
215 fEMCALRecoUtils->SetCutPhi(fPhicut);
216 }
6d67b5a7 217}
218
219//_____________________________________________________
220void AliEMCALTenderSupply::ProcessEvent()
221{
766cc9de 222 // Event loop.
6b7df55d 223
766cc9de 224 AliESDEvent *event = fTender->GetEvent();
225 if (!event) {
226 AliError("ESD event ptr = 0, returning");
227 return;
228 }
acf53135 229/*
230 if(!fRecParamSet)
231 {
232 if(!fRecParam)
233 {
234 InitRecParam();
235 fRecParamSet = kTRUE;
236 }
237 }
238 */
766cc9de 239 // Initialising parameters once per run number
240 if(fTender->RunChanged()){
241 GetPass();
acf53135 242 // Init bad channels
243 Int_t fInitBC=InitBadChannels();
244
245 if (fInitBC==0)
246 {
247 AliError("InitBadChannels returned false, returning");
248 return;
249 }
250 if(fInitBC>1)
251 {
252 AliInfo(Form("No external hot channel set: %d - %s", event->GetRunNumber(), fFilepass.Data()));
253 }
254
255 if (fReCalibCluster || fReCalibCell || fUpdateCell) {
256 Int_t fInitRecalib=InitRecalib();
257 if (fInitRecalib==0)
258 {
259 AliError("InitRecalib returned false, returning");
260 return;
261 }
262 if(fInitRecalib >1)
263 {
264 AliInfo(Form("No recalibration available: %d - %s", event->GetRunNumber(), fFilepass.Data()));
265 fReCalibCell=kFALSE;
266 fReCalibCluster=kFALSE;
267 }
766cc9de 268 }
acf53135 269
270 if (fRecalClusPos || fReClusterize || fUpdateCell) {
766cc9de 271 if (!InitMisalignMatrix()) {
272 AliError("InitMisalignmentMatrix returned false, returning");
273 return;
274 }
275 }
acf53135 276
277 if (fReClusterize || fUpdateCell) {
766cc9de 278 if (!InitClusterization()) {
279 AliError("InitClusterization returned false, returning");
280 return;
281 }
282 }
6b7df55d 283
766cc9de 284 if(fDebugLevel>1)
285 fEMCALRecoUtils->Print("");
acf53135 286
766cc9de 287 }
acf53135 288
766cc9de 289 // Test if cells present
290 AliESDCaloCells *cells= event->GetEMCALCells();
291 if (cells->GetNumberOfCells()<=0) {
292 AliWarning(Form("Number of EMCAL cells = %d, returning", cells->GetNumberOfCells()));
293 return;
294 }
6b7df55d 295
766cc9de 296 // Recalibrate cells
acf53135 297 if (fReCalibCell || fUpdateCell)
298 {
299 RecalibrateCells();
300 fReCalibCluster=kFALSE;
301 }
302 if(fDebugLevel>2)
303 AliInfo(Form("Re-calibrate cluster %d\n",fReCalibCluster));
304
305 if(fUpdateCell)
306 {
307 printf("Update cells\n");
308 UpdateCells();
309 fReClusterize=kTRUE;
310 }
311
312 // Reclusterize
313 if(fReClusterize)
314 {
315
316 FillDigitsArray();
317 Clusterize();
318 UpdateClusters();
319 }
320
766cc9de 321 // Store good clusters
322 TClonesArray *clusArr = dynamic_cast<TClonesArray*>(event->FindListObject("caloClusters"));
323 if (!clusArr)
324 clusArr = dynamic_cast<TClonesArray*>(event->FindListObject("CaloClusters"));
325 if (!clusArr) {
326 AliWarning(Form("No cluster array, number of cells in event = %d, returning", cells->GetNumberOfCells()));
327 return;
328 }
329 Int_t nclusters = clusArr->GetEntriesFast();
acf53135 330
766cc9de 331 for (Int_t icluster=0; icluster < nclusters; ++icluster) {
332 AliVCluster *clust = static_cast<AliVCluster*>(clusArr->At(icluster));
333 if (!clust)
334 continue;
335 if (!clust->IsEMCAL())
336 continue;
6b7df55d 337
766cc9de 338 if (fEMCALRecoUtils->ClusterContainsBadChannel(fEMCALGeo, clust->GetCellsAbsId(), clust->GetNCells())) {
339 delete clusArr->RemoveAt(icluster);
340 continue; //todo is it really needed to remove it? Or should we flag it?
341 }
6b7df55d 342
766cc9de 343 if (fFiducial){
344 if (!fEMCALRecoUtils->CheckCellFiducialRegion(fEMCALGeo, clust, cells)){
345 delete clusArr->RemoveAt(icluster);
346 continue; // todo it would be nice to store the distance
347 }
348 }
6b7df55d 349
766cc9de 350 fEMCALRecoUtils->CorrectClusterEnergyLinearity(clust);
351 if(fRecalDistToBadChannels)
352 fEMCALRecoUtils->RecalculateClusterDistanceToBadChannel(fEMCALGeo, cells, clust);
353 if(fReCalibCluster)
354 fEMCALRecoUtils->RecalibrateClusterEnergy(fEMCALGeo, clust, cells);
355 if(fRecalClusPos)
356 fEMCALRecoUtils->RecalculateClusterPosition(fEMCALGeo, cells, clust);
357 }
358 clusArr->Compress();
6b7df55d 359
766cc9de 360 // Track matching
361 if (!TGeoGlobalMagField::Instance()->GetField()) {
362 const AliESDRun *erun = event->GetESDRun();
363 AliMagF *field = AliMagF::CreateFieldMap(erun->GetCurrentL3(),
364 erun->GetCurrentDip(),
365 AliMagF::kConvLHC,
366 kFALSE,
367 erun->GetBeamEnergy(),
368 erun->GetBeamType());
369 TGeoGlobalMagField::Instance()->SetField(field);
370 }
a6a717c7 371
766cc9de 372 fEMCALRecoUtils->FindMatches(event,0x0,fEMCALGeo);
a6a717c7 373
374 SetClusterMatchedToTrack(event);
375 SetTracksMatchedToCluster(event);
376
6d67b5a7 377}
378
379//_____________________________________________________
380void AliEMCALTenderSupply::SetClusterMatchedToTrack(AliESDEvent *event)
381{
766cc9de 382 // Checks if tracks are matched to EMC clusters and set the matched EMCAL cluster index to ESD track.
6b7df55d 383
766cc9de 384 Int_t nTracks = event->GetNumberOfTracks();
385 for (Int_t iTrack = 0; iTrack < nTracks; ++iTrack) {
386 AliESDtrack* track = event->GetTrack(iTrack);
387 if (!track) {
388 AliWarning(Form("Could not receive track %d", iTrack));
389 continue;
390 }
391 Int_t matchClusIndex = fEMCALRecoUtils->GetMatchedClusterIndex(iTrack);
392 track->SetEMCALcluster(matchClusIndex); //sets -1 if track not matched within residual
393 if(matchClusIndex != -1)
394 track->SetStatus(AliESDtrack::kEMCALmatch);
acf53135 395 else
396 track->ResetStatus(AliESDtrack::kEMCALmatch);
766cc9de 397 }
398 if (fDebugLevel>2)
399 AliInfo("Track matched to closest cluster");
6d67b5a7 400}
401
402//_____________________________________________________
403void AliEMCALTenderSupply::SetTracksMatchedToCluster(AliESDEvent *event)
404{
766cc9de 405 // Checks if EMC clusters are matched to ESD track.
406 // Adds track indexes of all the tracks matched to a cluster withing residuals in ESDCalocluster.
6b7df55d 407
766cc9de 408 for (Int_t iClus=0; iClus < event->GetNumberOfCaloClusters(); ++iClus) {
409 AliESDCaloCluster *cluster = event->GetCaloCluster(iClus);
410 if (!cluster->IsEMCAL())
411 continue;
6b7df55d 412
766cc9de 413 Int_t nTracks = event->GetNumberOfTracks();
414 TArrayI arrayTrackMatched(nTracks);
6b7df55d 415
766cc9de 416 // Get the closest track matched to the cluster
417 Int_t nMatched = 0;
418 Int_t matchTrackIndex = fEMCALRecoUtils->GetMatchedTrackIndex(iClus);
419 if (matchTrackIndex != -1) {
420 arrayTrackMatched[nMatched] = matchTrackIndex;
421 nMatched++;
422 }
6b7df55d 423
766cc9de 424 // Get all other tracks matched to the cluster
425 for(Int_t iTrk=0; iTrk<nTracks; ++iTrk) {
426 AliESDtrack* track = event->GetTrack(iTrk);
427 if(iTrk == matchTrackIndex) continue;
428 if(track->GetEMCALcluster() == iClus){
429 arrayTrackMatched[nMatched] = iTrk;
430 ++nMatched;
431 }
432 }
6b7df55d 433
766cc9de 434 arrayTrackMatched.Set(nMatched);
435 cluster->AddTracksMatched(arrayTrackMatched);
436
437 Float_t eta= -999, phi = -999;
438 if (matchTrackIndex != -1)
439 fEMCALRecoUtils->GetMatchedResiduals(iClus, eta, phi);
440 cluster->SetTrackDistance(phi, eta);
441 }
6b7df55d 442
766cc9de 443 if (fDebugLevel>2)
444 AliInfo("Cluster matched to tracks");
6d67b5a7 445}
446
447//_____________________________________________________
448Bool_t AliEMCALTenderSupply::InitMisalignMatrix()
449{
766cc9de 450 // Initialising misalignment matrices
6b7df55d 451
766cc9de 452 AliESDEvent *event = fTender->GetEvent();
453 if (!event)
454 return kFALSE;
6b7df55d 455
766cc9de 456 if (fGeomMatrixSet) {
457 AliInfo("Misalignment matrix already set");
458 return kTRUE;
459 }
6b7df55d 460
766cc9de 461 if (fDebugLevel>0)
462 AliInfo("Initialising misalignment matrix");
6b7df55d 463
766cc9de 464 fEMCALRecoUtils->SetPositionAlgorithm(AliEMCALRecoUtils::kPosTowerGlobal);
6b7df55d 465
766cc9de 466 if (fLoadGeomMatrices) {
467 for(Int_t mod=0; mod < fEMCALGeo->GetNumberOfSuperModules(); ++mod) {
468 if (fEMCALMatrix[mod]){
469 if(fDebugLevel > 2)
470 fEMCALMatrix[mod]->Print();
471 fEMCALGeo->SetMisalMatrix(fEMCALMatrix[mod],mod);
472 }
473 }
474 fGeomMatrixSet = kTRUE;
475 return kTRUE;
476 }
6b7df55d 477
766cc9de 478 Int_t runGM = event->GetRunNumber();
6b7df55d 479 TObjArray *mobj = 0;
acf53135 480
481 if(fMisalignSurvey == kdefault){ //take default alignment corresponding to run no
482 AliOADBContainer emcalgeoCont(Form("emcal"));
483 emcalgeoCont.InitFromFile("$ALICE_ROOT/OADB/EMCAL/EMCALlocal2master.root",Form("AliEMCALgeo"));
484 mobj=(TObjArray*)emcalgeoCont.GetObject(runGM,"EmcalMatrices");
485 }
486
487 if(fMisalignSurvey == kSurveybyS){ //take alignment at sector level
766cc9de 488 if (runGM <= 140000) { //2010 data
6b7df55d 489 AliOADBContainer emcalgeoCont(Form("emcal2010"));
490 emcalgeoCont.InitFromFile("$ALICE_ROOT/OADB/EMCAL/EMCALlocal2master.root",Form("AliEMCALgeo"));
491 mobj=(TObjArray*)emcalgeoCont.GetObject(100,"survey10");
492
acf53135 493 } else if (runGM>140000) { // 2011 LHC11a pass1 data
766cc9de 494 AliOADBContainer emcalgeoCont(Form("emcal2011"));
6b7df55d 495 emcalgeoCont.InitFromFile("$ALICE_ROOT/OADB/EMCAL/EMCALlocal2master.root",Form("AliEMCALgeo"));
acf53135 496 mobj=(TObjArray*)emcalgeoCont.GetObject(100,"survey11byS");
497 }
498 }
499
500 if(fMisalignSurvey == kSurveybyS){ //take alignment at module level
501 if (runGM <= 140000) { //2010 data
502 AliOADBContainer emcalgeoCont(Form("emcal2010"));
503 emcalgeoCont.InitFromFile("$ALICE_ROOT/OADB/EMCAL/EMCALlocal2master.root",Form("AliEMCALgeo"));
504 mobj=(TObjArray*)emcalgeoCont.GetObject(100,"survey10");
505
506 } else if (runGM>140000) { // 2011 LHC11a pass1 data
507 AliOADBContainer emcalgeoCont(Form("emcal2011"));
508 emcalgeoCont.InitFromFile("$ALICE_ROOT/OADB/EMCAL/EMCALlocal2master.root",Form("AliEMCALgeo"));
509 mobj=(TObjArray*)emcalgeoCont.GetObject(100,"survey11byM");
766cc9de 510 }
acf53135 511 }
512
513 for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
514 fEMCALMatrix[mod] = (TGeoHMatrix*) mobj->At(mod);
515 fEMCALGeo->SetMisalMatrix(fEMCALMatrix[mod],mod);
516 fEMCALMatrix[mod]->Print();
517 }
766cc9de 518 return kTRUE;
6d67b5a7 519}
520
521//_____________________________________________________
acf53135 522Int_t AliEMCALTenderSupply::InitBadChannels()
6d67b5a7 523{
766cc9de 524 // Initialising bad channel maps
766cc9de 525 AliESDEvent *event = fTender->GetEvent();
526 if (!event)
acf53135 527 return 0;
6b7df55d 528
766cc9de 529 if (fDebugLevel>0)
530 AliInfo("Initialising Bad channel map");
6b7df55d 531
766cc9de 532 if (fFiducial){
533 fEMCALRecoUtils->SetNumberOfCellsFromEMCALBorder(fNCellsFromEMCALBorder);
534 fEMCALRecoUtils->SwitchOnNoFiducialBorderInEMCALEta0();
535 }
6b7df55d 536
766cc9de 537 fEMCALRecoUtils->SwitchOnBadChannelsRemoval();
538 if (fRecalDistToBadChannels)
539 fEMCALRecoUtils->SwitchOnDistToBadChannelRecalculation();
6b7df55d 540
766cc9de 541 Int_t runBC = event->GetRunNumber();
6b7df55d 542
acf53135 543 AliOADBContainer *contBC=new AliOADBContainer("");
544 if(fBasePath!=""){ //if fBasePath specified in the ->SetBasePath()
545 if (fDebugLevel>0) AliInfo(Form("Loading Bad Channels OADB from given path %s",fBasePath.Data()));
546 TFile *fbad=new TFile(Form("%s/EMCALBadChannels.root",fBasePath.Data()),"read");
547 if(fbad->IsZombie()){
548 AliFatal(Form("EMCALBadChannels.root was not found in the path provided: %s, aborting",fBasePath.Data()));
549 return 0;
550 }
551 if(fbad) delete fbad;
552 contBC->InitFromFile(Form("%s/EMCALBadChannels.root",fBasePath.Data()),"AliEMCALBadChannels");
553 }
554 else { // Else choose the one in the $ALICE_ROOT directory
555 if (fDebugLevel>0) AliInfo("Loading Bad Channels OADB from $ALICE_ROOT/OADB/EMCAL");
556 TFile *fbad=new TFile("$ALICE_ROOT/OADB/EMCAL/EMCALBadChannels.root","read");
557 if(fbad->IsZombie()){
558 AliFatal("$ALICE_ROOT/OADB/EMCAL/EMCALBadChannels.root was not found, aborting");
559 return 0;
560 }
561 if(fbad) delete fbad;
562 contBC->InitFromFile("$ALICE_ROOT/OADB/EMCAL/EMCALBadChannels.root","AliEMCALBadChannels");
766cc9de 563 }
acf53135 564
565 TObjArray *arrayBC=(TObjArray*)contBC->GetObject(runBC);
566 if(!arrayBC){
567 AliError(Form("No external hot channel set for run number: %d", runBC));
568 return 2;
766cc9de 569 }
6b7df55d 570
acf53135 571 TObjArray *arrayBCpass=(TObjArray*)arrayBC->FindObject(fFilepass); // Here, it looks for a specific pass
572 if(!arrayBCpass){
573 AliError(Form("No external hot channel set for: %d -%s", runBC,fFilepass.Data()));
574 return 2;
575 }
576
577 if(fDebugLevel>0) arrayBCpass->Print();
578
579 Int_t sms = fEMCALGeo->GetEMCGeometry()->GetNumberOfSuperModules();
580 for (Int_t i=0; i<sms; ++i) {
581 TH2I *h = fEMCALRecoUtils->GetEMCALChannelStatusMap(i);
582 if (h)
583 delete h;
584 h=(TH2I*)arrayBCpass->FindObject(Form("EMCALBadChannelMap_Mod%d",i));
585
586 if (!h) {
587 AliError(Form("Can not get EMCALBadChannelMap_Mod%d",i));
588 continue;
589 }
590 h->SetDirectory(0);
591 fEMCALRecoUtils->SetEMCALChannelStatusMap(i,h);
766cc9de 592 }
6b7df55d 593
acf53135 594 return 1;
6d67b5a7 595}
596
acf53135 597
6d67b5a7 598//_____________________________________________________
acf53135 599Int_t AliEMCALTenderSupply::InitRecalib()
6d67b5a7 600{
766cc9de 601 // Initialising recalibration factors.
6b7df55d 602
766cc9de 603 AliESDEvent *event = fTender->GetEvent();
604 if (!event)
acf53135 605 return 0;
6b7df55d 606
766cc9de 607 if (fDebugLevel>0)
608 AliInfo("Initialising recalibration factors");
6b7df55d 609
766cc9de 610 fEMCALRecoUtils->SwitchOnRecalibration();
611 Int_t runRC = event->GetRunNumber();
acf53135 612
613 AliOADBContainer *contRF=new AliOADBContainer("");
614 if(fBasePath!="") { //if fBasePath specified in the ->SetBasePath()
615 if (fDebugLevel>0) AliInfo(Form("Loading Recalib OADB from given path %s",fBasePath.Data()));
616 TFile *fRecalib= new TFile(Form("%s/EMCALRecalib.root",fBasePath.Data()),"read");
766cc9de 617 if (fRecalib->IsZombie()) {
acf53135 618 AliFatal(Form("EMCALRecalib.root not found in %s, aborting",fBasePath.Data()));
619 return 0;
766cc9de 620 }
acf53135 621 if(fRecalib) delete fRecalib;
622 contRF->InitFromFile(Form("%s/EMCALRecalib.root",fBasePath.Data()),"AliEMCALRecalib");
623 }
624 else{ // Else choose the one in the $ALICE_ROOT directory
625 if (fDebugLevel>0) AliInfo("Loading Recalib OADB from $ALICE_ROOT/OADB/EMCAL");
626 TFile *fRecalib= new TFile("$ALICE_ROOT/OADB/EMCAL/EMCALRecalib.root","read");
627 if (fRecalib->IsZombie()) {
628 AliFatal("$ALICE_ROOT/OADB/EMCAL/EMCALRecalib.root was not found, aborting");
629 return 0;
766cc9de 630 }
acf53135 631 if(fRecalib) delete fRecalib;
632 contRF->InitFromFile("$ALICE_ROOT/OADB/EMCAL/EMCALRecalib.root","AliEMCALRecalib");
766cc9de 633 }
acf53135 634
635 TObjArray *recal=(TObjArray*)contRF->GetObject(runRC); //GetObject(int runnumber)
636 if(!recal){
637 AliError(Form("No Objects for run: %d",runRC));
638 return 2;
639 }
640
641 TObjArray *recalpass=(TObjArray*)recal->FindObject(fFilepass);
642 if(!recalpass){
643 AliError(Form("No Objects for run: %d - %s",runRC,fFilepass.Data()));
644 return 2;
645 }
646
647 TObjArray *recalib=(TObjArray*)recalpass->FindObject("Recalib");
648 if(!recalib){
649 AliError(Form("No Recalib histos found for %d - %s",runRC,fFilepass.Data()));
650 return 2;
651 }
652
653 if(fDebugLevel>0) recalib->Print();
654
655 Int_t sms = fEMCALGeo->GetEMCGeometry()->GetNumberOfSuperModules();
656 for (Int_t i=0; i<sms; ++i) {
657 TH2F *h = fEMCALRecoUtils->GetEMCALChannelRecalibrationFactors(i);
658 if (h)
659 delete h;
660 h = (TH2F*)recalib->FindObject(Form("EMCALRecalFactors_SM%d",i));
661 if (!h) {
662 AliError(Form("Could not load EMCALRecalFactors_SM%d",i));
663 continue;
766cc9de 664 }
acf53135 665 h->SetDirectory(0);
666 fEMCALRecoUtils->SetEMCALChannelRecalibrationFactors(i,h);
766cc9de 667 }
acf53135 668 return 1;
c958a2f7 669}
670
671//_____________________________________________________
672void AliEMCALTenderSupply::RecalibrateCells()
673{
766cc9de 674 // Recalibrate cells.
6b7df55d 675
766cc9de 676 AliESDCaloCells *cells = fTender->GetEvent()->GetEMCALCells();
6b7df55d 677
766cc9de 678 Int_t nEMCcell = cells->GetNumberOfCells();
679 for(Int_t icell=0; icell<nEMCcell; ++icell) {
680 Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
681 fEMCALGeo->GetCellIndex(cells->GetCellNumber(icell),imod,iTower,iIphi,iIeta);
682 fEMCALGeo->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
683 Double_t calibFactor = fEMCALRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi);
684 cells->SetCell(icell,cells->GetCellNumber(icell),cells->GetAmplitude(icell)*calibFactor,cells->GetTime(icell));
685 }
b20bc239 686}
687
acf53135 688
689//_____________________________________________________
690void AliEMCALTenderSupply::UpdateCells()
691{
692 //Remove bad cells from the cell list
693 //This is required for later reclusterization
694
695 AliESDCaloCells *cells = fTender->GetEvent()->GetEMCALCells();
696
697 Int_t nEMCcell = cells->GetNumberOfCells();
698 for(Int_t icell=0; icell<nEMCcell; ++icell)
699 {
700 Int_t cellId = cells->GetCellNumber(icell);
701 Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
702 fEMCALGeo->GetCellIndex(cellId,imod,iTower,iIphi,iIeta);
703 fEMCALGeo->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
704 if(fEMCALRecoUtils->GetEMCALChannelStatus(imod,ieta,iphi))
705 {
706 cells->SetCell(icell,cellId,(-1)*cells->GetAmplitude(icell),cells->GetCellTime(cellId));
707 }
708 }
709}
710
711
712//_____________________________________________________
713void AliEMCALTenderSupply::InitRecParam()
714{
715 if (fDebugLevel>0)
716 AliInfo("Initialize the recParam");
717 fRecParam = new AliEMCALRecParam;
718}
719
720
b20bc239 721//_____________________________________________________
766cc9de 722Bool_t AliEMCALTenderSupply::InitClusterization()
b20bc239 723{
766cc9de 724 // Initialing clusterization/unfolding algorithm and set all the needed parameters.
6b7df55d 725
766cc9de 726 AliESDEvent *event=fTender->GetEvent();
727 if (!event)
728 return kFALSE;
6b7df55d 729
766cc9de 730 if (fDebugLevel>0)
731 AliInfo(Form("Initialising reclustering parameters: Clusterizer-%d",fRecParam->GetClusterizerFlag()));
6b7df55d 732
766cc9de 733 //---setup clusterizer
734 delete fClusterizer;
735 if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv1)
736 fClusterizer = new AliEMCALClusterizerv1 (fEMCALGeo);
737 else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv2)
738 fClusterizer = new AliEMCALClusterizerv2(fEMCALGeo);
739 else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerNxN) {
740 AliEMCALClusterizerNxN *clusterizer = new AliEMCALClusterizerNxN(fEMCALGeo);
741 clusterizer->SetNRowDiff(fRecParam->GetNRowDiff());
742 clusterizer->SetNColDiff(fRecParam->GetNColDiff());
743 fClusterizer = clusterizer;
744 } else {
745 AliFatal(Form("Clusterizer < %d > not available", fRecParam->GetClusterizerFlag()));
746 return kFALSE;
747 }
6b7df55d 748
766cc9de 749 // Set the clustering parameters
750 fClusterizer->SetECAClusteringThreshold( fRecParam->GetClusteringThreshold() );
751 fClusterizer->SetECALogWeight ( fRecParam->GetW0() );
752 fClusterizer->SetMinECut ( fRecParam->GetMinECut() );
753 fClusterizer->SetUnfolding ( fRecParam->GetUnfold() );
754 fClusterizer->SetECALocalMaxCut ( fRecParam->GetLocMaxCut() );
755 fClusterizer->SetTimeCut ( fRecParam->GetTimeCut() );
756 fClusterizer->SetTimeMin ( fRecParam->GetTimeMin() );
757 fClusterizer->SetTimeMax ( fRecParam->GetTimeMax() );
758 fClusterizer->SetInputCalibrated ( kTRUE );
759 fClusterizer->SetJustClusters ( kTRUE );
6b7df55d 760
766cc9de 761 // In case of unfolding after clusterization is requested, set the corresponding parameters
762 if (fRecParam->GetUnfold()) {
763 for (Int_t i = 0; i < 8; ++i) {
764 fClusterizer->SetSSPars(i, fRecParam->GetSSPars(i));
765 }
766 for (Int_t i = 0; i < 3; ++i) {
767 fClusterizer->SetPar5 (i, fRecParam->GetPar5(i));
768 fClusterizer->SetPar6 (i, fRecParam->GetPar6(i));
769 }
770 fClusterizer->InitClusterUnfolding();
771 }
6b7df55d 772
766cc9de 773 fClusterizer->SetDigitsArr(fDigitsArr);
774 fClusterizer->SetOutput(0);
775 fClusterArr = const_cast<TObjArray *>(fClusterizer->GetRecPoints());
776 return kTRUE;
b20bc239 777}
766cc9de 778
b20bc239 779//_____________________________________________________
780void AliEMCALTenderSupply::FillDigitsArray()
781{
766cc9de 782 // Fill digits from cells to a TClonesArray.
6b7df55d 783
766cc9de 784 AliESDEvent *event = fTender->GetEvent();
785 if (!event)
786 return;
6b7df55d 787
766cc9de 788 fDigitsArr->Clear("C");
789 AliESDCaloCells *cells = event->GetEMCALCells();
790 Int_t ncells = cells->GetNumberOfCells();
791 for (Int_t icell = 0, idigit = 0; icell < ncells; ++icell) {
792 Double_t cellAmplitude=0, cellTime=0;
793 Short_t cellNumber=0;
794 if (cells->GetCell(icell, cellNumber, cellAmplitude, cellTime) != kTRUE)
795 break;
796 AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->New(idigit));
797 digit->SetId(cellNumber);
798 digit->SetTime(cellTime);
799 digit->SetTimeR(cellTime);
800 digit->SetIndexInList(idigit);
801 digit->SetType(AliEMCALDigit::kHG);
802 digit->SetAmplitude(cellAmplitude);
803 idigit++;
804 }
b20bc239 805}
806
807//_____________________________________________________
808void AliEMCALTenderSupply::Clusterize()
809{
766cc9de 810 // Clusterize.
6b7df55d 811
766cc9de 812 fClusterizer->Digits2Clusters("");
b20bc239 813}
814
815//_____________________________________________________
816void AliEMCALTenderSupply::UpdateClusters()
817{
766cc9de 818 // Update ESD cluster list.
6b7df55d 819
766cc9de 820 AliESDEvent *event = fTender->GetEvent();
821 if (!event)
822 return;
6b7df55d 823
766cc9de 824 TClonesArray *clus = dynamic_cast<TClonesArray*>(event->FindListObject("caloClusters"));
825 if (!clus)
826 clus = dynamic_cast<TClonesArray*>(event->FindListObject("CaloClusters"));
827 if (!clus) {
6b7df55d 828 AliError(" Null pointer to calo clusters array, returning");
766cc9de 829 return;
830 }
6b7df55d 831
766cc9de 832 Int_t nents = clus->GetEntriesFast();
833 for (Int_t i=0; i < nents; ++i) {
834 AliESDCaloCluster *c = static_cast<AliESDCaloCluster*>(clus->At(i));
835 if (!c)
836 continue;
837 if (c->IsEMCAL())
838 delete clus->RemoveAt(i);
839 }
840 clus->Compress();
841 RecPoints2Clusters(clus);
b20bc239 842}
843
844//_____________________________________________________
845void AliEMCALTenderSupply::RecPoints2Clusters(TClonesArray *clus)
846{
766cc9de 847 // Convert AliEMCALRecoPoints to AliESDCaloClusters.
848 // Cluster energy, global position, cells and their amplitude fractions are restored.
849
850 AliESDEvent *event = fTender->GetEvent();
851 if (!event)
852 return;
6b7df55d 853
766cc9de 854 Int_t ncls = fClusterArr->GetEntriesFast();
855 for(Int_t i=0, nout=clus->GetEntriesFast(); i < ncls; ++i) {
856 AliEMCALRecPoint *recpoint = static_cast<AliEMCALRecPoint*>(fClusterArr->At(i));
857
858 Int_t ncells_true = 0;
859 const Int_t ncells = recpoint->GetMultiplicity();
860 UShort_t absIds[ncells];
861 Double32_t ratios[ncells];
862 Int_t *dlist = recpoint->GetDigitsList();
863 Float_t *elist = recpoint->GetEnergiesList();
864 for (Int_t c = 0; c < ncells; ++c) {
865 AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(dlist[c]));
866 absIds[ncells_true] = digit->GetId();
867 ratios[ncells_true] = elist[c]/digit->GetAmplitude();
868 if (ratios[ncells_true] < 0.001)
869 continue;
870 ++ncells_true;
871 }
872
873 if (ncells_true < 1) {
874 AliWarning("Skipping cluster with no cells");
875 continue;
876 }
877
878 // calculate new cluster position
879 TVector3 gpos;
880 recpoint->GetGlobalPosition(gpos);
881 Float_t g[3];
882 gpos.GetXYZ(g);
883
884 AliESDCaloCluster *c = static_cast<AliESDCaloCluster*>(clus->New(nout++));
885 c->SetType(AliVCluster::kEMCALClusterv1);
886 c->SetE(recpoint->GetEnergy());
887 c->SetPosition(g);
888 c->SetNCells(ncells_true);
889 c->SetDispersion(recpoint->GetDispersion());
890 c->SetEmcCpvDistance(-1); //not yet implemented
891 c->SetChi2(-1); //not yet implemented
892 c->SetTOF(recpoint->GetTime()) ; //time-of-flight
893 c->SetNExMax(recpoint->GetNExMax()); //number of local maxima
894 Float_t elipAxis[2];
895 recpoint->GetElipsAxis(elipAxis);
896 c->SetM02(elipAxis[0]*elipAxis[0]) ;
897 c->SetM20(elipAxis[1]*elipAxis[1]) ;
898 AliESDCaloCluster *cesd = static_cast<AliESDCaloCluster*>(c);
899 cesd->SetCellsAbsId(absIds);
900 cesd->SetCellsAmplitudeFraction(ratios);
901 }
902}
b20bc239 903
766cc9de 904//_____________________________________________________
905void AliEMCALTenderSupply::GetPass()
906{
907 // Get passx from filename.
6b7df55d 908
766cc9de 909 AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
910 fInputTree = mgr->GetTree();
6b7df55d 911
766cc9de 912 if (!fInputTree) {
6b7df55d 913 AliError("Pointer to tree = 0, returning");
766cc9de 914 return;
915 }
6b7df55d 916
766cc9de 917 fInputFile = fInputTree->GetCurrentFile();
918 if (!fInputFile) {
6b7df55d 919 AliError("Null pointer input file, returning");
766cc9de 920 return;
921 }
6b7df55d 922
766cc9de 923 TString fname(fInputFile->GetName());
924 if (fname.Contains("pass1")) fFilepass = TString("pass1");
925 else if(fname.Contains("pass2")) fFilepass = TString("pass2");
926 else if(fname.Contains("pass3")) fFilepass = TString("pass3");
927 else {
928 AliError(Form("Pass number string not found: %s", fname.Data()));
929 return;
930 }
b20bc239 931}
932