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 ,fUpdateCell(kFALSE)
  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)
  82 ,fBasePath("")
  83 ,fReClusterize(kFALSE)
  84 ,fClusterizer(0)
  85 ,fGeomMatrixSet(kFALSE)
  86 ,fLoadGeomMatrices(kFALSE)
  87 ,fRecParam(0x0)
  88 ,fRecParamSet(kFALSE)
  89 ,fOCDBpath(" ")
  90 ,fUnfolder(0)
  91 ,fDigitsArr(0)
  92 ,fClusterArr(0)
  93 ,fMisalignSurvey(kdefault)
  94 {
  95   // Default constructor.
  96   for(Int_t i = 0; i < 10; i++) fEMCALMatrix[i] = 0 ;
  97 }
  98
  99 //_____________________________________________________
 100 AliEMCALTenderSupply::AliEMCALTenderSupply(const char *name, const AliTender *tender) :
 101 AliTenderSupply(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)
 111 ,fUpdateCell(kFALSE)
 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)
 126 ,fBasePath("")
 127 ,fReClusterize(kFALSE)
 128 ,fClusterizer(0)
 129 ,fGeomMatrixSet(kFALSE)
 130 ,fLoadGeomMatrices(kFALSE)
 131 ,fRecParam(0x0)
 132 ,fRecParamSet(kFALSE)
 133 ,fOCDBpath(" ")
 134 ,fUnfolder(0)
 135 ,fDigitsArr(0)
 136 ,fClusterArr(0)
 137 ,fMisalignSurvey(kdefault)
 138 {
 139   // Named constructor
 140
 141   for(Int_t i = 0; i < 10; i++) fEMCALMatrix[i] = 0 ;
 142
 143   fEMCALRecoUtils  = new AliEMCALRecoUtils();
 144   fRecParam        = new AliEMCALRecParam;
 145   fDigitsArr       = new TClonesArray("AliEMCALDigit",1000);
 146 }
 147
 148 //_____________________________________________________
 149 AliEMCALTenderSupply::~AliEMCALTenderSupply()
 150 {
 151   //Destructor
 152
 153   delete fEMCALRecoUtils;
 154   delete fClusterizer;
 155   delete fUnfolder;
 156   if (fDigitsArr){
 157     fDigitsArr->Clear("C");
 158     delete fDigitsArr;
 159   }
 160 }
 161
 162 //_____________________________________________________
 163 void AliEMCALTenderSupply::Init()
 164 {
 165   // Initialise EMCAL tender.
 166
 167   if (fDebugLevel>0)
 168     AliInfo("Init EMCAL Tender supply");
 169
 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   }
 198
 199   // Init goemetry
 200   fEMCALGeo = AliEMCALGeometry::GetInstance(fEMCALGeoName) ;
 201
 202   // Initialising Non linearity parameters
 203   fEMCALRecoUtils->SetNonLinearityThreshold(fNonLinearThreshold);
 204   fEMCALRecoUtils->SetNonLinearityFunction(fNonLinearFunc);
 205
 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   }
 217 }
 218
 219 //_____________________________________________________
 220 void AliEMCALTenderSupply::ProcessEvent()
 221 {
 222   // Event loop.
 223
 224   AliESDEvent *event = fTender->GetEvent();
 225   if (!event) {
 226     AliError("ESD event ptr = 0, returning");
 227     return;
 228   }
 229 /*
 230   if(!fRecParamSet)
 231    {
 232     if(!fRecParam)
 233     {
 234       InitRecParam();
 235       fRecParamSet = kTRUE;
 236     }
 237    }
 238  */
 239   // Initialising parameters once per run number
 240   if(fTender->RunChanged()){
 241     GetPass();
 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         }
 268     }
 269
 270     if (fRecalClusPos || fReClusterize || fUpdateCell) {
 271       if (!InitMisalignMatrix()) {
 272         AliError("InitMisalignmentMatrix returned false, returning");
 273         return;
 274       }
 275     }
 276
 277     if (fReClusterize || fUpdateCell) {
 278       if (!InitClusterization()) {
 279         AliError("InitClusterization returned false, returning");
 280         return;
 281       }
 282     }
 283
 284     if(fDebugLevel>1)
 285       fEMCALRecoUtils->Print("");
 286
 287   }
 288
 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   }
 295
 296   // Recalibrate cells
 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
 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();
 330
 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;
 337
 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     }
 342
 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     }
 349
 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();
 359
 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   }
 371
 372   fEMCALRecoUtils->FindMatches(event,0x0,fEMCALGeo);
 373
 374   SetClusterMatchedToTrack(event);
 375   SetTracksMatchedToCluster(event);
 376
 377 }
 378
 379 //_____________________________________________________
 380 void AliEMCALTenderSupply::SetClusterMatchedToTrack(AliESDEvent *event)
 381 {
 382   // Checks if tracks are matched to EMC clusters and set the matched EMCAL cluster index to ESD track.
 383
 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);
 395     else
 396       track->ResetStatus(AliESDtrack::kEMCALmatch);
 397   }
 398   if (fDebugLevel>2)
 399     AliInfo("Track matched to closest cluster");
 400 }
 401
 402 //_____________________________________________________
 403 void AliEMCALTenderSupply::SetTracksMatchedToCluster(AliESDEvent *event)
 404 {
 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.
 407
 408   for (Int_t iClus=0; iClus < event->GetNumberOfCaloClusters(); ++iClus) {
 409     AliESDCaloCluster *cluster = event->GetCaloCluster(iClus);
 410     if (!cluster->IsEMCAL())
 411       continue;
 412
 413     Int_t nTracks = event->GetNumberOfTracks();
 414     TArrayI arrayTrackMatched(nTracks);
 415
 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     }
 423
 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     }
 433
 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   }
 442
 443   if (fDebugLevel>2)
 444     AliInfo("Cluster matched to tracks");
 445 }
 446
 447 //_____________________________________________________
 448 Bool_t AliEMCALTenderSupply::InitMisalignMatrix()
 449 {
 450   // Initialising misalignment matrices
 451
 452   AliESDEvent *event = fTender->GetEvent();
 453   if (!event)
 454     return kFALSE;
 455
 456   if (fGeomMatrixSet) {
 457     AliInfo("Misalignment matrix already set");
 458     return kTRUE;
 459   }
 460
 461   if (fDebugLevel>0)
 462     AliInfo("Initialising misalignment matrix");
 463
 464   fEMCALRecoUtils->SetPositionAlgorithm(AliEMCALRecoUtils::kPosTowerGlobal);
 465
 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   }
 477
 478   Int_t runGM = event->GetRunNumber();
 479   TObjArray *mobj = 0;
 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
 488   if (runGM <= 140000) { //2010 data
 489     AliOADBContainer emcalgeoCont(Form("emcal2010"));
 490     emcalgeoCont.InitFromFile("$ALICE_ROOT/OADB/EMCAL/EMCALlocal2master.root",Form("AliEMCALgeo"));
 491     mobj=(TObjArray*)emcalgeoCont.GetObject(100,"survey10");
 492
 493   } else if (runGM>140000) { // 2011 LHC11a pass1 data
 494     AliOADBContainer emcalgeoCont(Form("emcal2011"));
 495     emcalgeoCont.InitFromFile("$ALICE_ROOT/OADB/EMCAL/EMCALlocal2master.root",Form("AliEMCALgeo"));
 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");
 510   }
 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  }
 518   return kTRUE;
 519 }
 520
 521 //_____________________________________________________
 522 Int_t AliEMCALTenderSupply::InitBadChannels()
 523 {
 524   // Initialising bad channel maps
 525   AliESDEvent *event = fTender->GetEvent();
 526   if (!event)
 527     return 0;
 528
 529   if (fDebugLevel>0)
 530     AliInfo("Initialising Bad channel map");
 531
 532   if (fFiducial){
 533     fEMCALRecoUtils->SetNumberOfCellsFromEMCALBorder(fNCellsFromEMCALBorder);
 534     fEMCALRecoUtils->SwitchOnNoFiducialBorderInEMCALEta0();
 535   }
 536
 537   fEMCALRecoUtils->SwitchOnBadChannelsRemoval();
 538   if (fRecalDistToBadChannels)
 539     fEMCALRecoUtils->SwitchOnDistToBadChannelRecalculation();
 540
 541   Int_t runBC = event->GetRunNumber();
 542
 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");
 563     }
 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;
 569   }
 570
 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);
 592   }
 593
 594   return 1;
 595 }
 596
 597
 598 //_____________________________________________________
 599 Int_t AliEMCALTenderSupply::InitRecalib()
 600 {
 601   // Initialising recalibration factors.
 602
 603   AliESDEvent *event = fTender->GetEvent();
 604   if (!event)
 605     return 0;
 606
 607   if (fDebugLevel>0)
 608     AliInfo("Initialising recalibration factors");
 609
 610   fEMCALRecoUtils->SwitchOnRecalibration();
 611   Int_t runRC = event->GetRunNumber();
 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");
 617     if (fRecalib->IsZombie()) {
 618       AliFatal(Form("EMCALRecalib.root not found in %s, aborting",fBasePath.Data()));
 619       return 0;
 620     }
 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;
 630       }
 631       if(fRecalib) delete fRecalib;
 632       contRF->InitFromFile("$ALICE_ROOT/OADB/EMCAL/EMCALRecalib.root","AliEMCALRecalib");
 633     }
 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;
 664     }
 665     h->SetDirectory(0);
 666     fEMCALRecoUtils->SetEMCALChannelRecalibrationFactors(i,h);
 667   }
 668   return 1;
 669 }
 670
 671 //_____________________________________________________
 672 void AliEMCALTenderSupply::RecalibrateCells()
 673 {
 674   // Recalibrate cells.
 675
 676   AliESDCaloCells *cells = fTender->GetEvent()->GetEMCALCells();
 677
 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   }
 686 }
 687
 688
 689 //_____________________________________________________
 690 void 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 //_____________________________________________________
 713 void AliEMCALTenderSupply::InitRecParam()
 714 {
 715   if (fDebugLevel>0)
 716     AliInfo("Initialize the recParam");
 717   fRecParam = new AliEMCALRecParam;
 718 }
 719
 720
 721 //_____________________________________________________
 722 Bool_t AliEMCALTenderSupply::InitClusterization()
 723 {
 724   // Initialing clusterization/unfolding algorithm and set all the needed parameters.
 725
 726   AliESDEvent *event=fTender->GetEvent();
 727   if (!event)
 728     return kFALSE;
 729
 730   if (fDebugLevel>0)
 731     AliInfo(Form("Initialising reclustering parameters: Clusterizer-%d",fRecParam->GetClusterizerFlag()));
 732
 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   }
 748
 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                               );
 760
 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   }
 772
 773   fClusterizer->SetDigitsArr(fDigitsArr);
 774   fClusterizer->SetOutput(0);
 775   fClusterArr = const_cast<TObjArray *>(fClusterizer->GetRecPoints());
 776   return kTRUE;
 777 }
 778
 779 //_____________________________________________________
 780 void AliEMCALTenderSupply::FillDigitsArray()
 781 {
 782   // Fill digits from cells to a TClonesArray.
 783
 784   AliESDEvent *event = fTender->GetEvent();
 785   if (!event)
 786     return;
 787
 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   }
 805 }
 806
 807 //_____________________________________________________
 808 void AliEMCALTenderSupply::Clusterize()
 809 {
 810   // Clusterize.
 811
 812   fClusterizer->Digits2Clusters("");
 813 }
 814
 815 //_____________________________________________________
 816 void AliEMCALTenderSupply::UpdateClusters()
 817 {
 818   // Update ESD cluster list.
 819
 820   AliESDEvent *event = fTender->GetEvent();
 821   if (!event)
 822     return;
 823
 824   TClonesArray *clus = dynamic_cast<TClonesArray*>(event->FindListObject("caloClusters"));
 825   if (!clus)
 826     clus = dynamic_cast<TClonesArray*>(event->FindListObject("CaloClusters"));
 827   if (!clus) {
 828     AliError(" Null pointer to calo clusters array, returning");
 829     return;
 830   }
 831
 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);
 842 }
 843
 844 //_____________________________________________________
 845 void AliEMCALTenderSupply::RecPoints2Clusters(TClonesArray *clus)
 846 {
 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;
 853
 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 }
 903
 904 //_____________________________________________________
 905 void AliEMCALTenderSupply::GetPass()
 906 {
 907   // Get passx from filename.
 908
 909   AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
 910   fInputTree = mgr->GetTree();
 911
 912   if (!fInputTree) {
 913     AliError("Pointer to tree = 0, returning");
 914     return;
 915   }
 916
 917   fInputFile = fInputTree->GetCurrentFile();
 918   if (!fInputFile) {
 919     AliError("Null pointer input file, returning");
 920     return;
 921   }
 922
 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   }
 931 }
 932