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