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