ANALYSIS/TenderSupplies/AliEMCALTenderSupply.cxx

   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                        //
  20 //  and do track matching.                                                   //
  21 //  Author: Deepa Thomas (Utrecht University)                                //
  22 //                                                                           //
  23 ///////////////////////////////////////////////////////////////////////////////
  24
  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>
  32 #include <AliLog.h>
  33 #include <AliESDEvent.h>
  34 #include <AliAnalysisManager.h>
  35 #include <AliTender.h>
  36 #include "AliOADBContainer.h"
  37 #include "AliCDBManager.h"
  38 #include "AliCDBStorage.h"
  39 #include "AliCDBEntry.h"
  40 #include "AliMagF.h"
  41 #include "AliESDCaloCluster.h"
  42 #include "AliEMCALTenderSupply.h"
  43 #include "AliEMCALGeometry.h"
  44 #include "AliEMCALRecoUtils.h"
  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"
  53
  54 ClassImp(AliEMCALTenderSupply)
  55
  56 AliEMCALTenderSupply::AliEMCALTenderSupply() :
  57 AliTenderSupply()
  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)
  91 {
  92   // Default constructor.
  93   for(Int_t i = 0; i < 10; i++) fEMCALMatrix[i] = 0 ;
  94 }
  95
  96 //_____________________________________________________
  97 AliEMCALTenderSupply::AliEMCALTenderSupply(const char *name, const AliTender *tender) :
  98 AliTenderSupply(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)
 132 {
 133   // Named constructor
 134
 135   for(Int_t i = 0; i < 10; i++) fEMCALMatrix[i] = 0 ;
 136
 137   fRecParam        = new AliEMCALRecParam;
 138   fEMCALRecoUtils  = new AliEMCALRecoUtils();
 139   fDigitsArr       = new TClonesArray("AliEMCALDigit",1000);
 140 }
 141
 142 //_____________________________________________________
 143 AliEMCALTenderSupply::~AliEMCALTenderSupply()
 144 {
 145   //Destructor
 146
 147   delete fEMCALRecoUtils;
 148   delete fClusterizer;
 149   delete fUnfolder;
 150   if (fDigitsArr){
 151     fDigitsArr->Clear("C");
 152     delete fDigitsArr;
 153   }
 154 }
 155
 156 //_____________________________________________________
 157 void AliEMCALTenderSupply::Init()
 158 {
 159   // Initialise EMCAL tender.
 160
 161   if (fDebugLevel>0)
 162     AliInfo("Init EMCAL Tender supply");
 163
 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   }
 192
 193   // Init goemetry
 194   fEMCALGeo = AliEMCALGeometry::GetInstance(fEMCALGeoName) ;
 195
 196   // Initialising Non linearity parameters
 197   fEMCALRecoUtils->SetNonLinearityThreshold(fNonLinearThreshold);
 198   fEMCALRecoUtils->SetNonLinearityFunction(fNonLinearFunc);
 199
 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   }
 211 }
 212
 213 //_____________________________________________________
 214 void AliEMCALTenderSupply::ProcessEvent()
 215 {
 216   // Event loop.
 217
 218   AliESDEvent *event = fTender->GetEvent();
 219   if (!event) {
 220     AliError("ESD event ptr = 0, returning");
 221     return;
 222   }
 223
 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     }
 249
 250     if(fDebugLevel>1)
 251       fEMCALRecoUtils->Print("");
 252   }
 253
 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   }
 260
 261   // Recalibrate cells
 262   if (fReCalibCell)
 263     RecalibrateCells();
 264
 265   // Reclusterize
 266   if(fReClusterize){
 267     FillDigitsArray();
 268     Clusterize();
 269     UpdateClusters();
 270   }
 271
 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;
 287
 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     }
 292
 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     }
 299
 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();
 309
 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   }
 327 }
 328
 329 //_____________________________________________________
 330 void AliEMCALTenderSupply::SetClusterMatchedToTrack(AliESDEvent *event)
 331 {
 332   // Checks if tracks are matched to EMC clusters and set the matched EMCAL cluster index to ESD track.
 333
 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");
 348 }
 349
 350 //_____________________________________________________
 351 void AliEMCALTenderSupply::SetTracksMatchedToCluster(AliESDEvent *event)
 352 {
 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.
 355
 356   for (Int_t iClus=0; iClus < event->GetNumberOfCaloClusters(); ++iClus) {
 357     AliESDCaloCluster *cluster = event->GetCaloCluster(iClus);
 358     if (!cluster->IsEMCAL())
 359       continue;
 360
 361     Int_t nTracks = event->GetNumberOfTracks();
 362     TArrayI arrayTrackMatched(nTracks);
 363
 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     }
 371
 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     }
 381
 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   }
 390
 391   if (fDebugLevel>2)
 392     AliInfo("Cluster matched to tracks");
 393 }
 394
 395 //_____________________________________________________
 396 Bool_t AliEMCALTenderSupply::InitMisalignMatrix()
 397 {
 398   // Initialising misalignment matrices
 399
 400   AliESDEvent *event = fTender->GetEvent();
 401   if (!event)
 402     return kFALSE;
 403
 404   if (fGeomMatrixSet) {
 405     AliInfo("Misalignment matrix already set");
 406     return kTRUE;
 407   }
 408
 409   if (fDebugLevel>0)
 410     AliInfo("Initialising misalignment matrix");
 411
 412   fEMCALRecoUtils->SetPositionAlgorithm(AliEMCALRecoUtils::kPosTowerGlobal);
 413
 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   }
 425
 426   Int_t runGM = event->GetRunNumber();
 427   TObjArray *mobj = 0;
 428
 429   if (runGM <= 140000) { //2010 data
 430     AliOADBContainer emcalgeoCont(Form("emcal2010"));
 431     emcalgeoCont.InitFromFile("$ALICE_ROOT/OADB/EMCAL/EMCALlocal2master.root",Form("AliEMCALgeo"));
 432     mobj=(TObjArray*)emcalgeoCont.GetObject(100,"survey10");
 433
 434   } else if (runGM>140000 && runGM <148531 && (fFilepass == "pass1")) { // 2011 LHC11a pass1 data
 435     AliOADBContainer emcalgeoCont(Form("emcal2011"));
 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"))){
 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;
 457 }
 458
 459 //_____________________________________________________
 460 Bool_t AliEMCALTenderSupply::InitBadChannels()
 461 {
 462   // Initialising bad channel maps
 463
 464   AliESDEvent *event = fTender->GetEvent();
 465   if (!event)
 466     return kFALSE;
 467
 468   if (fDebugLevel>0)
 469     AliInfo("Initialising Bad channel map");
 470
 471   if (fFiducial){
 472     fEMCALRecoUtils->SetNumberOfCellsFromEMCALBorder(fNCellsFromEMCALBorder);
 473     fEMCALRecoUtils->SwitchOnNoFiducialBorderInEMCALEta0();
 474   }
 475
 476   fEMCALRecoUtils->SwitchOnBadChannelsRemoval();
 477   if (fRecalDistToBadChannels)
 478     fEMCALRecoUtils->SwitchOnDistToBadChannelRecalculation();
 479
 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   }
 507
 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
 511
 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   }
 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};
 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   }
 542
 543   AliInfo(Form("No external hot channel set: %d - %s", runBC, fFilepass.Data()));
 544   return kTRUE;
 545 }
 546
 547 //_____________________________________________________
 548 Bool_t AliEMCALTenderSupply::InitRecalib()
 549 {
 550   // Initialising recalibration factors.
 551
 552   AliESDEvent *event = fTender->GetEvent();
 553   if (!event)
 554     return kFALSE;
 555
 556   if (fDebugLevel>0)
 557     AliInfo("Initialising recalibration factors");
 558
 559   fEMCALRecoUtils->SwitchOnRecalibration();
 560   Int_t runRC = event->GetRunNumber();
 561
 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");
 570
 571         if(     (runRC >= 114737 && runRC <= 117223 && (fFilepass = "pass2")) || //LHC10b pass2
 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
 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     }
 610
 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;
 627 }
 628
 629 //_____________________________________________________
 630 void AliEMCALTenderSupply::RecalibrateCells()
 631 {
 632   // Recalibrate cells.
 633
 634   AliESDCaloCells *cells = fTender->GetEvent()->GetEMCALCells();
 635
 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   }
 644 }
 645
 646 //_____________________________________________________
 647 Bool_t AliEMCALTenderSupply::InitClusterization()
 648 {
 649   // Initialing clusterization/unfolding algorithm and set all the needed parameters.
 650
 651   AliESDEvent *event=fTender->GetEvent();
 652   if (!event)
 653     return kFALSE;
 654
 655   if (fDebugLevel>0)
 656     AliInfo(Form("Initialising reclustering parameters: Clusterizer-%d",fRecParam->GetClusterizerFlag()));
 657
 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   }
 673
 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                               );
 685
 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   }
 697
 698   fClusterizer->SetDigitsArr(fDigitsArr);
 699   fClusterizer->SetOutput(0);
 700   fClusterArr = const_cast<TObjArray *>(fClusterizer->GetRecPoints());
 701   return kTRUE;
 702 }
 703
 704 //_____________________________________________________
 705 void AliEMCALTenderSupply::FillDigitsArray()
 706 {
 707   // Fill digits from cells to a TClonesArray.
 708
 709   AliESDEvent *event = fTender->GetEvent();
 710   if (!event)
 711     return;
 712
 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   }
 730 }
 731
 732 //_____________________________________________________
 733 void AliEMCALTenderSupply::Clusterize()
 734 {
 735   // Clusterize.
 736
 737   fClusterizer->Digits2Clusters("");
 738 }
 739
 740 //_____________________________________________________
 741 void AliEMCALTenderSupply::UpdateClusters()
 742 {
 743   // Update ESD cluster list.
 744
 745   AliESDEvent *event = fTender->GetEvent();
 746   if (!event)
 747     return;
 748
 749   TClonesArray *clus = dynamic_cast<TClonesArray*>(event->FindListObject("caloClusters"));
 750   if (!clus)
 751     clus = dynamic_cast<TClonesArray*>(event->FindListObject("CaloClusters"));
 752   if (!clus) {
 753     AliError(" Null pointer to calo clusters array, returning");
 754     return;
 755   }
 756
 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);
 767 }
 768
 769 //_____________________________________________________
 770 void AliEMCALTenderSupply::RecPoints2Clusters(TClonesArray *clus)
 771 {
 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;
 778
 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 }
 828
 829 //_____________________________________________________
 830 void AliEMCALTenderSupply::GetPass()
 831 {
 832   // Get passx from filename.
 833
 834   AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
 835   fInputTree = mgr->GetTree();
 836
 837   if (!fInputTree) {
 838     AliError("Pointer to tree = 0, returning");
 839     return;
 840   }
 841
 842   fInputFile = fInputTree->GetCurrentFile();
 843   if (!fInputFile) {
 844     AliError("Null pointer input file, returning");
 845     return;
 846   }
 847
 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   }
 856 }
 857