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