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 //  EMCAL tender, apply corrections to EMCAL clusters and do track matching. //
  19 //                                                                           //
  20 //  Author: Deepa Thomas (Utrecht University)                                //
  21 //  Later mods/rewrite: Jiri Kral (University of Jyvaskyla)                  //
  22 //  S. Aiola, C. Loizides : Make it work for AODs                            //
  23 //                                                                           //
  24 ///////////////////////////////////////////////////////////////////////////////
  25
  26 #include <TROOT.h>
  27 #include <TFile.h>
  28 #include <TTree.h>
  29 #include <TInterpreter.h>
  30 #include <TObjArray.h>
  31 #include <TClonesArray.h>
  32 #include <TGeoGlobalMagField.h>
  33 #include "AliEMCALAfterBurnerUF.h"
  34 #include "AliEMCALClusterizer.h"
  35 #include "AliEMCALClusterizerNxN.h"
  36 #include "AliEMCALClusterizerv1.h"
  37 #include "AliEMCALClusterizerv2.h"
  38 #include "AliEMCALDigit.h"
  39 #include "AliEMCALGeometry.h"
  40 #include "AliEMCALRecParam.h"
  41 #include "AliEMCALRecParam.h"
  42 #include "AliEMCALRecPoint.h"
  43 #include "AliEMCALRecoUtils.h"
  44 #include "AliEMCALTenderSupply.h"
  45 #include "AliESDCaloCluster.h"
  46 #include "AliMagF.h"
  47 #include "AliOADBContainer.h"
  48 #include "AliAODEvent.h"
  49 #include "AliAnalysisManager.h"
  50 #include "AliESDEvent.h"
  51 #include "AliLog.h"
  52 #include "AliTender.h"
  53
  54 ClassImp(AliEMCALTenderSupply)
  55
  56 AliEMCALTenderSupply::AliEMCALTenderSupply() :
  57 AliTenderSupply()
  58 ,fTask(0)
  59 ,fRun(0)
  60 ,fEMCALGeo(0x0)
  61 ,fEMCALGeoName("")
  62 ,fEMCALRecoUtils(0)
  63 ,fConfigName("")
  64 ,fDebugLevel(0)
  65 ,fNonLinearFunc(-1)
  66 ,fNonLinearThreshold(-1)
  67 ,fReCalibCluster(kFALSE)
  68 ,fUpdateCell(kFALSE)
  69 ,fCalibrateEnergy(kFALSE)
  70 ,fCalibrateTime(kFALSE)
  71 ,fDoNonLinearity(kFALSE)
  72 ,fBadCellRemove(kFALSE)
  73 ,fRejectExoticCells(kFALSE)
  74 ,fRejectExoticClusters(kFALSE)
  75 ,fClusterBadChannelCheck(kFALSE)
  76 ,fRecalClusPos(kFALSE)
  77 ,fFiducial(kFALSE)
  78 ,fNCellsFromEMCALBorder(-1)
  79 ,fRecalDistToBadChannels(kFALSE)
  80 ,fRecalShowerShape(kFALSE)
  81 ,fInputTree(0)
  82 ,fInputFile(0)
  83 ,fGetPassFromFileName(kTRUE)
  84 ,fFilepass(0)
  85 ,fMass(-1)
  86 ,fStep(-1)
  87 ,fCutEtaPhiSum(kTRUE)
  88 ,fCutEtaPhiSeparate(kFALSE)
  89 ,fRcut(-1)
  90 ,fEtacut(-1)
  91 ,fPhicut(-1)
  92 ,fBasePath("")
  93 ,fReClusterize(kFALSE)
  94 ,fClusterizer(0)
  95 ,fGeomMatrixSet(kFALSE)
  96 ,fLoadGeomMatrices(kFALSE)
  97 ,fRecParam(0x0)
  98 ,fDoTrackMatch(kFALSE)
  99 ,fDoUpdateOnly(kFALSE)
 100 ,fUnfolder(0)
 101 ,fDigitsArr(0)
 102 ,fClusterArr(0)
 103 ,fMisalignSurvey(kdefault)
 104 ,fExoticCellFraction(-1)
 105 ,fExoticCellDiffTime(-1)
 106 ,fExoticCellMinAmplitude(-1)
 107 ,fSetCellMCLabelFromCluster(kFALSE)
 108 {
 109   // Default constructor.
 110
 111   for(Int_t i = 0; i < 12; i++) fEMCALMatrix[i] = 0 ;
 112 }
 113
 114 //_____________________________________________________
 115 AliEMCALTenderSupply::AliEMCALTenderSupply(const char *name, const AliTender *tender) :
 116 AliTenderSupply(name,tender)
 117 ,fTask(0)
 118 ,fRun(0)
 119 ,fEMCALGeo(0x0)
 120 ,fEMCALGeoName("")
 121 ,fEMCALRecoUtils(0)
 122 ,fConfigName("")
 123 ,fDebugLevel(0)
 124 ,fNonLinearFunc(-1)
 125 ,fNonLinearThreshold(-1)
 126 ,fReCalibCluster(kFALSE)
 127 ,fUpdateCell(kFALSE)
 128 ,fCalibrateEnergy(kFALSE)
 129 ,fCalibrateTime(kFALSE)
 130 ,fDoNonLinearity(kFALSE)
 131 ,fBadCellRemove(kFALSE)
 132 ,fRejectExoticCells(kFALSE)
 133 ,fRejectExoticClusters(kFALSE)
 134 ,fClusterBadChannelCheck(kFALSE)
 135 ,fRecalClusPos(kFALSE)
 136 ,fFiducial(kFALSE)
 137 ,fNCellsFromEMCALBorder(-1)
 138 ,fRecalDistToBadChannels(kFALSE)
 139 ,fRecalShowerShape(kFALSE)
 140 ,fInputTree(0)
 141 ,fInputFile(0)
 142 ,fGetPassFromFileName(kTRUE)
 143 ,fFilepass("")
 144 ,fMass(-1)
 145 ,fStep(-1)
 146 ,fCutEtaPhiSum(kTRUE)
 147 ,fCutEtaPhiSeparate(kFALSE)
 148 ,fRcut(-1)
 149 ,fEtacut(-1)
 150 ,fPhicut(-1)
 151 ,fBasePath("")
 152 ,fReClusterize(kFALSE)
 153 ,fClusterizer(0)
 154 ,fGeomMatrixSet(kFALSE)
 155 ,fLoadGeomMatrices(kFALSE)
 156 ,fRecParam(0x0)
 157 ,fDoTrackMatch(kFALSE)
 158 ,fDoUpdateOnly(kFALSE)
 159 ,fUnfolder(0)
 160 ,fDigitsArr(0)
 161 ,fClusterArr(0)
 162 ,fMisalignSurvey(kdefault)
 163 ,fExoticCellFraction(-1)
 164 ,fExoticCellDiffTime(-1)
 165 ,fExoticCellMinAmplitude(-1)
 166 ,fSetCellMCLabelFromCluster(kFALSE)
 167 {
 168   // Named constructor
 169
 170   for(Int_t i = 0; i < 12; i++) fEMCALMatrix[i] = 0 ;
 171 }
 172
 173 //_____________________________________________________
 174 AliEMCALTenderSupply::AliEMCALTenderSupply(const char *name, AliAnalysisTaskSE *task) :
 175 AliTenderSupply(name)
 176 ,fTask(task)
 177 ,fRun(0)
 178 ,fEMCALGeo(0x0)
 179 ,fEMCALGeoName("")
 180 ,fEMCALRecoUtils(0)
 181 ,fConfigName("")
 182 ,fDebugLevel(0)
 183 ,fNonLinearFunc(-1)
 184 ,fNonLinearThreshold(-1)
 185 ,fReCalibCluster(kFALSE)
 186 ,fUpdateCell(kFALSE)
 187 ,fCalibrateEnergy(kFALSE)
 188 ,fCalibrateTime(kFALSE)
 189 ,fDoNonLinearity(kFALSE)
 190 ,fBadCellRemove(kFALSE)
 191 ,fRejectExoticCells(kFALSE)
 192 ,fRejectExoticClusters(kFALSE)
 193 ,fClusterBadChannelCheck(kFALSE)
 194 ,fRecalClusPos(kFALSE)
 195 ,fFiducial(kFALSE)
 196 ,fNCellsFromEMCALBorder(-1)
 197 ,fRecalDistToBadChannels(kFALSE)
 198 ,fRecalShowerShape(kFALSE)
 199 ,fInputTree(0)
 200 ,fInputFile(0)
 201 ,fGetPassFromFileName(kTRUE)
 202 ,fFilepass("")
 203 ,fMass(-1)
 204 ,fStep(-1)
 205 ,fCutEtaPhiSum(kTRUE)
 206 ,fCutEtaPhiSeparate(kFALSE)
 207 ,fRcut(-1)
 208 ,fEtacut(-1)
 209 ,fPhicut(-1)
 210 ,fBasePath("")
 211 ,fReClusterize(kFALSE)
 212 ,fClusterizer(0)
 213 ,fGeomMatrixSet(kFALSE)
 214 ,fLoadGeomMatrices(kFALSE)
 215 ,fRecParam(0x0)
 216 ,fDoTrackMatch(kFALSE)
 217 ,fDoUpdateOnly(kFALSE)
 218 ,fUnfolder(0)
 219 ,fDigitsArr(0)
 220 ,fClusterArr(0)
 221 ,fMisalignSurvey(kdefault)
 222 ,fExoticCellFraction(-1)
 223 ,fExoticCellDiffTime(-1)
 224 ,fExoticCellMinAmplitude(-1)
 225 ,fSetCellMCLabelFromCluster(kFALSE)
 226 {
 227   // Named constructor.
 228
 229   for(Int_t i = 0; i < 12; i++) fEMCALMatrix[i] = 0 ;
 230 }
 231
 232 //_____________________________________________________
 233 AliEMCALTenderSupply::~AliEMCALTenderSupply()
 234 {
 235   //Destructor
 236
 237   if (!AliAnalysisManager::GetAnalysisManager()->IsProofMode())
 238   {
 239     delete fEMCALRecoUtils;
 240     delete fRecParam;
 241     delete fUnfolder;
 242     if (!fClusterizer)
 243     {
 244       fDigitsArr->Clear("C");
 245       delete fDigitsArr;
 246     } else
 247     {
 248       delete fClusterizer;
 249       fDigitsArr = 0;
 250     }
 251   }
 252 }
 253
 254 //_____________________________________________________
 255 void AliEMCALTenderSupply::SetDefaults()
 256 {
 257   // Set default settings.
 258
 259   SwitchOnBadCellRemove();
 260   SwitchOnExoticCellRemove();
 261   SwitchOnCalibrateEnergy();
 262   SwitchOnCalibrateTime();
 263   SwitchOnUpdateCell();
 264   SwitchOnReclustering();
 265   SwitchOnClusterBadChannelCheck();
 266   SwitchOnClusterExoticChannelCheck();
 267   SwitchOnTrackMatch();
 268 }
 269
 270 //_____________________________________________________
 271 Bool_t AliEMCALTenderSupply::RunChanged() const
 272 {
 273   // Get run number.
 274
 275   return (fTender && fTender->RunChanged()) || (fTask && fRun != fTask->InputEvent()->GetRunNumber());
 276 }
 277
 278 //_____________________________________________________
 279 void AliEMCALTenderSupply::Init()
 280 {
 281   // Initialise EMCAL tender.
 282
 283   if (fDebugLevel>0)
 284     AliWarning("Init EMCAL Tender supply");
 285
 286   if (fConfigName.Length()>0 && gROOT->LoadMacro(fConfigName) >=0) {
 287     AliDebug(1, Form("Loading settings from macro %s", fConfigName.Data()));
 288     AliEMCALTenderSupply *tender = (AliEMCALTenderSupply*)gInterpreter->ProcessLine("ConfigEMCALTenderSupply()");
 289     fDebugLevel             = tender->fDebugLevel;
 290     fEMCALGeoName           = tender->fEMCALGeoName;
 291     fEMCALRecoUtils         = tender->fEMCALRecoUtils;
 292     fConfigName             = tender->fConfigName;
 293     fNonLinearFunc          = tender->fNonLinearFunc;
 294     fNonLinearThreshold     = tender->fNonLinearThreshold;
 295     fReCalibCluster         = tender->fReCalibCluster;
 296     fUpdateCell             = tender->fUpdateCell;
 297     fRecalClusPos           = tender->fRecalClusPos;
 298     fCalibrateEnergy        = tender->fCalibrateEnergy;
 299     fCalibrateTime          = tender->fCalibrateTime;
 300     fFiducial               = tender->fFiducial;
 301     fNCellsFromEMCALBorder  = tender->fNCellsFromEMCALBorder;
 302     fRecalDistToBadChannels = tender->fRecalDistToBadChannels;
 303     fRecalShowerShape       = tender->fRecalShowerShape;
 304     fClusterBadChannelCheck = tender->fClusterBadChannelCheck;
 305     fBadCellRemove          = tender->fBadCellRemove;
 306     fRejectExoticCells      = tender->fRejectExoticCells;
 307     fRejectExoticClusters   = tender->fRejectExoticClusters;
 308     fMass                   = tender->fMass;
 309     fStep                   = tender->fStep;
 310     fCutEtaPhiSum           = tender->fCutEtaPhiSum;
 311     fCutEtaPhiSeparate      = tender->fCutEtaPhiSeparate;
 312     fRcut                   = tender->fRcut;
 313     fEtacut                 = tender->fEtacut;
 314     fPhicut                 = tender->fPhicut;
 315     fReClusterize           = tender->fReClusterize;
 316     fLoadGeomMatrices       = tender->fLoadGeomMatrices;
 317     fRecParam               = tender->fRecParam;
 318     fDoNonLinearity         = tender->fDoNonLinearity;
 319     fDoTrackMatch           = tender->fDoTrackMatch;
 320     fDoUpdateOnly           = tender->fDoUpdateOnly;
 321     fMisalignSurvey         = tender->fMisalignSurvey;
 322     fExoticCellFraction     = tender->fExoticCellFraction;
 323     fExoticCellDiffTime     = tender->fExoticCellDiffTime;
 324     fExoticCellMinAmplitude = tender->fExoticCellMinAmplitude;
 325
 326     for(Int_t i = 0; i < 12; i++)
 327       fEMCALMatrix[i] = tender->fEMCALMatrix[i] ;
 328   }
 329
 330   if (fDebugLevel>0){
 331     AliInfo( "Emcal Tender settings: ======================================" );
 332     AliInfo( "------------ Switches --------------------------" );
 333     AliInfo( Form( "BadCellRemove : %d", fBadCellRemove ));
 334     AliInfo( Form( "ExoticCellRemove : %d", fRejectExoticCells ));
 335     AliInfo( Form( "CalibrateEnergy : %d", fCalibrateEnergy ));
 336     AliInfo( Form( "CalibrateTime : %d", fCalibrateTime ));
 337     AliInfo( Form( "UpdateCell : %d", fUpdateCell ));
 338     AliInfo( Form( "DoUpdateOnly : %d", fDoUpdateOnly ));
 339     AliInfo( Form( "Reclustering : %d", fReClusterize ));
 340     AliInfo( Form( "ClusterBadChannelCheck : %d", fClusterBadChannelCheck ));
 341     AliInfo( Form( "ClusterExoticChannelCheck : %d", fRejectExoticClusters ));
 342     AliInfo( Form( "CellFiducialRegion : %d", fFiducial ));
 343     AliInfo( Form( "ReCalibrateCluster : %d", fReCalibCluster ));
 344     AliInfo( Form( "RecalculateClusPos : %d", fRecalClusPos ));
 345     AliInfo( Form( "RecalShowerShape : %d", fRecalShowerShape ));
 346     AliInfo( Form( "NonLinearityCorrection : %d", fDoNonLinearity ));
 347     AliInfo( Form( "RecalDistBadChannel : %d", fRecalDistToBadChannels ));
 348     AliInfo( Form( "TrackMatch : %d", fDoTrackMatch ));
 349     AliInfo( "------------ Variables -------------------------" );
 350     AliInfo( Form( "DebugLevel : %d", fDebugLevel ));
 351     AliInfo( Form( "BasePath : %s", fBasePath.Data() ));
 352     AliInfo( Form( "ConfigFileName : %s", fConfigName.Data() ));
 353     AliInfo( Form( "EMCALGeometryName : %s", fEMCALGeoName.Data() ));
 354     AliInfo( Form( "NonLinearityFunction : %d", fNonLinearFunc ));
 355     AliInfo( Form( "NonLinearityThreshold : %d", fNonLinearThreshold ));
 356     AliInfo( Form( "MisalignmentMatrixSurvey : %d", fMisalignSurvey ));
 357     AliInfo( Form( "NumberOfCellsFromEMCALBorder : %d", fNCellsFromEMCALBorder ));
 358     AliInfo( Form( "RCut : %f", fRcut ));
 359     AliInfo( Form( "Mass : %f", fMass ));
 360     AliInfo( Form( "Step : %f", fStep ));
 361     AliInfo( Form( "EtaCut : %f", fEtacut ));
 362     AliInfo( Form( "PhiCut : %f", fPhicut ));
 363     AliInfo( Form( "ExoticCellFraction : %f", fExoticCellFraction ));
 364     AliInfo( Form( "ExoticCellDiffTime : %f", fExoticCellDiffTime ));
 365     AliInfo( Form( "ExoticCellMinAmplitude : %f", fExoticCellMinAmplitude ));
 366     AliInfo( "=============================================================" );
 367   }
 368
 369   // init reco utils
 370
 371   if(!fEMCALRecoUtils)
 372     fEMCALRecoUtils  = new AliEMCALRecoUtils;
 373
 374   // init geometry if requested
 375   if (fEMCALGeoName.Length()>0)
 376     fEMCALGeo = AliEMCALGeometry::GetInstance(fEMCALGeoName) ;
 377
 378   // digits array
 379   fDigitsArr       = new TClonesArray("AliEMCALDigit",1000);
 380
 381   // initialising non-linearity parameters
 382   if( fNonLinearThreshold != -1 )
 383     fEMCALRecoUtils->SetNonLinearityThreshold(fNonLinearThreshold);
 384   if( fNonLinearFunc != -1 )
 385     fEMCALRecoUtils->SetNonLinearityFunction(fNonLinearFunc);
 386
 387   // missalignment function
 388   fEMCALRecoUtils->SetPositionAlgorithm(AliEMCALRecoUtils::kPosTowerGlobal);
 389
 390   // fiducial cut
 391   // do not do the eta0 fiducial cut
 392   if( fNCellsFromEMCALBorder != -1 )
 393     fEMCALRecoUtils->SetNumberOfCellsFromEMCALBorder(fNCellsFromEMCALBorder);
 394   fEMCALRecoUtils->SwitchOnNoFiducialBorderInEMCALEta0();
 395
 396   // exotic cell rejection
 397   if( fExoticCellFraction != -1 )
 398     fEMCALRecoUtils->SetExoticCellFractionCut( fExoticCellFraction );
 399   if( fExoticCellDiffTime != -1 )
 400     fEMCALRecoUtils->SetExoticCellDiffTimeCut( fExoticCellDiffTime );
 401   if( fExoticCellMinAmplitude != -1 )
 402     fEMCALRecoUtils->SetExoticCellMinAmplitudeCut( fExoticCellMinAmplitude );
 403
 404   // setting track matching parameters ... mass, step size and residual cut
 405   if( fMass != -1 )
 406     fEMCALRecoUtils->SetMass(fMass);
 407   if( fStep != -1 )
 408     fEMCALRecoUtils->SetStep(fStep);
 409
 410   // spatial cut based on separate eta/phi or common processing
 411   if(fCutEtaPhiSum){
 412     fEMCALRecoUtils->SwitchOnCutEtaPhiSum();
 413     if( fRcut != -1 )
 414       fEMCALRecoUtils->SetCutR(fRcut);
 415   } else if (fCutEtaPhiSeparate) {
 416     fEMCALRecoUtils->SwitchOnCutEtaPhiSeparate();
 417     if( fEtacut != -1 )
 418       fEMCALRecoUtils->SetCutEta(fEtacut);
 419     if( fPhicut != -1 )
 420       fEMCALRecoUtils->SetCutPhi(fPhicut);
 421   }
 422 }
 423
 424 //_____________________________________________________
 425 AliVEvent* AliEMCALTenderSupply::GetEvent()
 426 {
 427   // Return the event pointer.
 428
 429   if (fTender) {
 430     return fTender->GetEvent();
 431   }
 432   else if (fTask) {
 433     return fTask->InputEvent();
 434   }
 435
 436   return 0;
 437 }
 438
 439 //_____________________________________________________
 440 void AliEMCALTenderSupply::ProcessEvent()
 441 {
 442   // Event loop.
 443
 444   AliVEvent *event = GetEvent();
 445
 446   if (!event) {
 447     AliError("Event ptr = 0, returning");
 448     return;
 449   }
 450
 451  // Initialising parameters once per run number
 452   if (RunChanged()) {
 453
 454     AliWarning( "Run changed, initializing parameters" );
 455     fRun = event->GetRunNumber();
 456
 457     // init geometry if not already done
 458     if (fEMCALGeoName.Length()==0) {
 459       fEMCALGeoName = "EMCAL_FIRSTYEARV1";
 460       if (fRun>139517) {
 461         fEMCALGeoName = "EMCAL_COMPLETEV1";
 462       }
 463       if (fRun>170593) {
 464         fEMCALGeoName = "EMCAL_COMPLETE12SMV1";
 465       }
 466       fEMCALGeo = AliEMCALGeometry::GetInstance(fEMCALGeoName);
 467       if (!fEMCALGeo) {
 468         AliFatal(Form("Can not create geometry: %s", fEMCALGeoName.Data()));
 469         return;
 470       }
 471     }
 472
 473     // get pass
 474     if (fGetPassFromFileName)
 475       GetPass();
 476
 477     // define what recalib parameters are needed for various switches
 478     // this is based on implementation in AliEMCALRecoUtils
 479     Bool_t needRecoParam   = fReClusterize;
 480     Bool_t needBadChannels = fBadCellRemove   | fClusterBadChannelCheck | fRecalDistToBadChannels | fReClusterize;
 481     Bool_t needRecalib     = fCalibrateEnergy | fReClusterize;
 482     Bool_t needTimecalib   = fCalibrateTime   | fReClusterize;
 483     Bool_t needMisalign    = fRecalClusPos    | fReClusterize;
 484     Bool_t needClusterizer = fReClusterize;
 485
 486     // initiate reco params with some defaults
 487     // will not overwrite, if those have been provided by user
 488     if( needRecoParam ){
 489       Int_t initRC = InitRecParam();
 490
 491       if( initRC == 0 )
 492         AliInfo("Defaults reco params loaded.");
 493       if( initRC > 1 )
 494         AliWarning("User defined reco params.");
 495     }
 496
 497     // init bad channels
 498     if( needBadChannels ){
 499       Int_t fInitBC = InitBadChannels();
 500       if (fInitBC==0)
 501         AliError("InitBadChannels returned false, returning");
 502       if (fInitBC==1)
 503         AliWarning("InitBadChannels OK");
 504       if (fInitBC>1)
 505         AliWarning(Form("No external hot channel set: %d - %s", event->GetRunNumber(), fFilepass.Data()));
 506     }
 507
 508     // init recalibration factors
 509     if( needRecalib )
 510     {
 511       Int_t fInitRecalib = InitRecalib();
 512       if (fInitRecalib==0)
 513         AliError("InitRecalib returned false, returning");
 514       if (fInitRecalib==1)
 515         AliWarning("InitRecalib OK");
 516       if (fInitRecalib>1)
 517         AliWarning(Form("No recalibration available: %d - %s", event->GetRunNumber(), fFilepass.Data()));
 518
 519       Int_t fInitRunDepRecalib = InitRunDepRecalib();
 520       if (fInitRunDepRecalib==0)
 521         AliError("InitrunDepRecalib returned false, returning");
 522       if (fInitRunDepRecalib==1)
 523         AliWarning("InitRecalib OK");
 524       if (fInitRunDepRecalib>1)
 525         AliWarning(Form("No Temperature recalibration available: %d - %s", event->GetRunNumber(), fFilepass.Data()));
 526
 527     }
 528
 529
 530
 531     // init time calibration
 532     if( needTimecalib ){
 533       Int_t initTC = InitTimeCalibration();
 534       if ( !initTC )
 535         AliError("InitTimeCalibration returned false, returning");
 536       if (initTC==1)
 537         AliWarning("InitTimeCalib OK");
 538       if( initTC > 1 )
 539         AliWarning(Form("No external time calibration set: %d - %s", event->GetRunNumber(), fFilepass.Data()));
 540     }
 541
 542     // init misalignment matrix
 543     if( needMisalign ) {
 544       if (!InitMisalignMatrix())
 545         AliError("InitMisalignmentMatrix returned false, returning");
 546       else
 547         AliWarning("InitMisalignMatrix OK");
 548     }
 549
 550     // init clusterizer
 551     if( needClusterizer ) {
 552       if (!InitClusterization())
 553         AliError("InitClusterization returned false, returning");
 554       else
 555         AliWarning("InitClusterization OK");
 556     }
 557
 558     if(fDebugLevel>1)
 559       fEMCALRecoUtils->Print("");
 560   }
 561
 562   // disable implied switches -------------------------------------------------
 563   // AliEMCALRecoUtils or clusterizer functions alredy include some
 564   // recalibration so based on those implied callibration te switches
 565   // here are disabled to avoid duplication
 566
 567   // clusterizer does cluster energy recalibration, position recomputation
 568   // and shower shape
 569   if( fReClusterize ){
 570     fReCalibCluster   = kFALSE;
 571     fRecalClusPos     = kFALSE;
 572     fRecalShowerShape = kFALSE;
 573   }
 574
 575   // CONFIGURE THE RECO UTILS -------------------------------------------------
 576   // configure the reco utils
 577   // this option does energy recalibration
 578   if( fCalibrateEnergy )
 579     fEMCALRecoUtils->SwitchOnRecalibration();
 580   else
 581     fEMCALRecoUtils->SwitchOffRecalibration();
 582
 583   // allows time calibration
 584   if( fCalibrateTime )
 585     fEMCALRecoUtils->SwitchOnTimeRecalibration();
 586   else
 587     fEMCALRecoUtils->SwitchOffTimeRecalibration();
 588
 589   // allows to zero bad cells
 590   if( fBadCellRemove )
 591     fEMCALRecoUtils->SwitchOnBadChannelsRemoval();
 592   else
 593     fEMCALRecoUtils->SwitchOffBadChannelsRemoval();
 594
 595   // distance to bad channel recalibration
 596   if( fRecalDistToBadChannels )
 597     fEMCALRecoUtils->SwitchOnDistToBadChannelRecalculation();
 598   else
 599     fEMCALRecoUtils->SwitchOffDistToBadChannelRecalculation();
 600
 601   // exclude exotic cells
 602   if( fRejectExoticCells )
 603     fEMCALRecoUtils->SwitchOnRejectExoticCell();
 604   else
 605     fEMCALRecoUtils->SwitchOffRejectExoticCell();
 606
 607   // exclude clusters with exotic cells
 608   if( fRejectExoticClusters )
 609     fEMCALRecoUtils->SwitchOnRejectExoticCluster();
 610   else
 611     fEMCALRecoUtils->SwitchOffRejectExoticCluster();
 612
 613   // TODO: not implemented switches
 614   // SwitchOnClusterEnergySmearing
 615   // SwitchOnRunDepCorrection
 616
 617   // START PROCESSING ---------------------------------------------------------
 618   // Test if cells present
 619   AliVCaloCells *cells= event->GetEMCALCells();
 620   if (cells->GetNumberOfCells()<=0)
 621   {
 622     if(fDebugLevel>1)
 623       AliWarning(Form("Number of EMCAL cells = %d, returning", cells->GetNumberOfCells()));
 624     return;
 625   }
 626
 627   if (fDebugLevel>2)
 628     AliInfo(Form("Re-calibrate cluster %d\n",fReCalibCluster));
 629
 630   // mark the cells not recalibrated in case of selected
 631   // time, energy recalibration or bad channel removal
 632   if( fCalibrateEnergy || fCalibrateTime || fBadCellRemove )
 633     fEMCALRecoUtils->ResetCellsCalibrated();
 634
 635  // CELL RECALIBRATION -------------------------------------------------------
 636   // cell objects will be updated
 637   // the cell calibrations are also processed locally any time those are needed
 638   // in case that the cell objects are not to be updated here for later use
 639   if( fUpdateCell )
 640   {
 641     // do the update
 642     // includes exotic cell check in the UpdateCells function - is not provided
 643     // by the reco utils
 644     UpdateCells();
 645
 646     // switch off recalibrations so those are not done multiple times
 647     // this is just for safety, the recalibrated flag of cell object
 648     // should not allow for farther processing anyways
 649     fEMCALRecoUtils->SwitchOffRecalibration();
 650     fEMCALRecoUtils->SwitchOffTimeRecalibration();
 651
 652     if (fDoUpdateOnly)
 653       return;
 654   }
 655
 656   // RECLUSTERIZATION ---------------------------------------------------------
 657   if (fReClusterize)
 658   {
 659     FillDigitsArray();
 660     Clusterize();
 661     UpdateClusters();
 662   }
 663
 664   // Store good clusters
 665   TClonesArray *clusArr = dynamic_cast<TClonesArray*>(event->FindListObject("caloClusters"));
 666   if (!clusArr)
 667     clusArr = dynamic_cast<TClonesArray*>(event->FindListObject("CaloClusters"));
 668   if (!clusArr) {
 669     AliWarning(Form("No cluster array, number of cells in event = %d, returning", cells->GetNumberOfCells()));
 670     return;
 671   }
 672
 673   // PROCESS SINGLE CLUSTER RECALIBRATION -------------------------------------
 674   // now go through clusters one by one and process separate correction
 675   // as those were defined or not
 676   Int_t nclusters = clusArr->GetEntriesFast();
 677   for (Int_t icluster=0; icluster < nclusters; ++icluster)
 678   {
 679     AliVCluster *clust = static_cast<AliVCluster*>(clusArr->At(icluster));
 680     if (!clust)
 681       continue;
 682     if  (!clust->IsEMCAL())
 683       continue;
 684
 685     // REMOVE CLUSTERS WITH BAD CELLS -----------------------------
 686     if( fClusterBadChannelCheck )
 687     {
 688       // careful, the the ClusterContainsBadChannel is dependent on
 689       // SwitchOnBadChannelsRemoval, switching it ON automatically
 690       // and returning to original value after processing
 691       Bool_t badRemoval = fEMCALRecoUtils->IsBadChannelsRemovalSwitchedOn();
 692       fEMCALRecoUtils->SwitchOnBadChannelsRemoval();
 693
 694       Bool_t badResult = fEMCALRecoUtils->ClusterContainsBadChannel(fEMCALGeo, clust->GetCellsAbsId(), clust->GetNCells());
 695
 696       // switch the bad channels removal back
 697       if( ! badRemoval )
 698         fEMCALRecoUtils->SwitchOffBadChannelsRemoval();
 699
 700       if( badResult )
 701       {
 702         delete clusArr->RemoveAt(icluster);
 703         continue; //TODO is it really needed to remove it? Or should we flag it?
 704       }
 705     }
 706
 707     // REMOVE EXOTIC CLUSTERS -------------------------------------
 708     // does process local cell recalibration energy and time without replacing
 709     // the global cell values, in case of no cell recalib done yet
 710     if( fRejectExoticClusters )
 711     {
 712       // careful, the IsExoticCluster is dependent on
 713       // SwitchOnRejectExoticCell, switching it ON automatically
 714       // and returning to original value after processing
 715       Bool_t exRemoval = fEMCALRecoUtils->IsRejectExoticCell();
 716       fEMCALRecoUtils->SwitchOnRejectExoticCell();
 717
 718       // get bunch crossing
 719       Int_t bunchCrossNo = event->GetBunchCrossNumber();
 720
 721       Bool_t exResult = fEMCALRecoUtils->IsExoticCluster(clust, cells, bunchCrossNo );
 722
 723       // switch the exotic channels removal back
 724       if( ! exRemoval )
 725         fEMCALRecoUtils->SwitchOffRejectExoticCell();
 726
 727       if( exResult )
 728       {
 729         delete clusArr->RemoveAt(icluster);
 730         continue; //TODO is it really needed to remove it? Or should we flag it?
 731       }
 732     }
 733
 734     // FIDUCIAL CUT -----------------------------------------------
 735     if (fFiducial)
 736     {
 737       // depends on SetNumberOfCellsFromEMCALBorder
 738       // SwitchOnNoFiducialBorderInEMCALEta0
 739       if (!fEMCALRecoUtils->CheckCellFiducialRegion(fEMCALGeo, clust, cells)){
 740         delete clusArr->RemoveAt(icluster);
 741         continue; //TODO it would be nice to store the distance
 742       }
 743     }
 744
 745     // CLUSTER ENERGY ---------------------------------------------
 746     // does process local cell recalibration energy and time without replacing
 747     // the global cell values, in case of no cell recalib done yet
 748     if( fReCalibCluster ) {
 749       fEMCALRecoUtils->RecalibrateClusterEnergy(fEMCALGeo, clust, cells);
 750       if (clust->E() < 1e-9) {
 751         delete clusArr->RemoveAt(icluster);
 752         continue;
 753       }
 754     }
 755
 756     // CLUSTER POSITION -------------------------------------------
 757     // does process local cell energy recalibration, if enabled and cells
 758     // not calibrated yet
 759     if( fRecalClusPos )
 760       fEMCALRecoUtils->RecalculateClusterPosition(fEMCALGeo, cells, clust);
 761
 762     // SHOWER SHAPE -----------------------------------------------
 763     if( fRecalShowerShape )
 764       fEMCALRecoUtils->RecalculateClusterShowerShapeParameters(fEMCALGeo, cells, clust);
 765
 766     // NONLINEARITY -----------------------------------------------
 767     if( fDoNonLinearity )
 768     {
 769       Float_t correctedEnergy = fEMCALRecoUtils->CorrectClusterEnergyLinearity(clust);
 770       clust->SetE(correctedEnergy);
 771     }
 772
 773     // DISTANCE TO BAD CHANNELS -----------------------------------
 774     if( fRecalDistToBadChannels )
 775       fEMCALRecoUtils->RecalculateClusterDistanceToBadChannel(fEMCALGeo, cells, clust);
 776   }
 777
 778   clusArr->Compress();
 779
 780   if (!fDoTrackMatch)
 781     return;
 782
 783   // TRACK MATCHING -----------------------------------------------------------
 784   if (!TGeoGlobalMagField::Instance()->GetField()) {
 785     AliESDEvent *esd = dynamic_cast<AliESDEvent*>(event);
 786     if (esd)
 787       esd->InitMagneticField();
 788     else {
 789       AliAODEvent *aod = dynamic_cast<AliAODEvent*>(event);
 790       Double_t curSol = 30000*aod->GetMagneticField()/5.00668;
 791       Double_t curDip = 6000 *aod->GetMuonMagFieldScale();
 792       AliMagF *field  = AliMagF::CreateFieldMap(curSol,curDip);
 793       TGeoGlobalMagField::Instance()->SetField(field);
 794     }
 795   }
 796
 797   fEMCALRecoUtils->FindMatches(event,0x0,fEMCALGeo);
 798   fEMCALRecoUtils->SetClusterMatchedToTrack(event);
 799   fEMCALRecoUtils->SetTracksMatchedToCluster(event);
 800 }
 801
 802 //_____________________________________________________
 803 Bool_t AliEMCALTenderSupply::InitMisalignMatrix()
 804 {
 805   // Initialising misalignment matrices
 806
 807   AliVEvent *event = GetEvent();
 808
 809   if (!event)
 810     return kFALSE;
 811
 812   if (fGeomMatrixSet)
 813   {
 814     AliInfo("Misalignment matrix already set");
 815     return kTRUE;
 816   }
 817
 818   if (fDebugLevel>0)
 819     AliInfo("Initialising misalignment matrix");
 820
 821   if (fLoadGeomMatrices) {
 822     for(Int_t mod=0; mod < fEMCALGeo->GetNumberOfSuperModules(); ++mod)
 823     {
 824       if (fEMCALMatrix[mod]){
 825         if(fDebugLevel > 2)
 826           fEMCALMatrix[mod]->Print();
 827         fEMCALGeo->SetMisalMatrix(fEMCALMatrix[mod],mod);
 828       }
 829     }
 830     fGeomMatrixSet = kTRUE;
 831     return kTRUE;
 832   }
 833
 834   Int_t runGM = event->GetRunNumber();
 835   TObjArray *mobj = 0;
 836
 837   if(fMisalignSurvey == kdefault)
 838   { //take default alignment corresponding to run no
 839     AliOADBContainer emcalgeoCont(Form("emcal"));
 840     emcalgeoCont.InitFromFile("$ALICE_ROOT/OADB/EMCAL/EMCALlocal2master.root",Form("AliEMCALgeo"));
 841     mobj=(TObjArray*)emcalgeoCont.GetObject(runGM,"EmcalMatrices");
 842   }
 843
 844   if(fMisalignSurvey == kSurveybyS)
 845   { //take alignment at sector level
 846     if (runGM <= 140000) { //2010 data
 847       AliOADBContainer emcalgeoCont(Form("emcal2010"));
 848       emcalgeoCont.InitFromFile("$ALICE_ROOT/OADB/EMCAL/EMCALlocal2master.root",Form("AliEMCALgeo"));
 849       mobj=(TObjArray*)emcalgeoCont.GetObject(100,"survey10");
 850     }
 851     else if (runGM>140000)
 852     { // 2011 LHC11a pass1 data
 853       AliOADBContainer emcalgeoCont(Form("emcal2011"));
 854       emcalgeoCont.InitFromFile("$ALICE_ROOT/OADB/EMCAL/EMCALlocal2master.root",Form("AliEMCALgeo"));
 855       mobj=(TObjArray*)emcalgeoCont.GetObject(100,"survey11byS");
 856     }
 857   }
 858
 859   if(fMisalignSurvey == kSurveybyM)
 860   { //take alignment at module level
 861     if (runGM <= 140000) { //2010 data
 862       AliOADBContainer emcalgeoCont(Form("emcal2010"));
 863       emcalgeoCont.InitFromFile("$ALICE_ROOT/OADB/EMCAL/EMCALlocal2master.root",Form("AliEMCALgeo"));
 864       mobj=(TObjArray*)emcalgeoCont.GetObject(100,"survey10");
 865     }
 866     else if (runGM>140000)
 867     { // 2011 LHC11a pass1 data
 868       AliOADBContainer emcalgeoCont(Form("emcal2011"));
 869       emcalgeoCont.InitFromFile("$ALICE_ROOT/OADB/EMCAL/EMCALlocal2master.root",Form("AliEMCALgeo"));
 870       mobj=(TObjArray*)emcalgeoCont.GetObject(100,"survey11byM");
 871     }
 872   }
 873
 874   if(!mobj)
 875   {
 876     AliFatal("Geometry matrix array not found");
 877     return kFALSE;
 878   }
 879
 880   for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++)
 881   {
 882     fEMCALMatrix[mod] = (TGeoHMatrix*) mobj->At(mod);
 883     fEMCALGeo->SetMisalMatrix(fEMCALMatrix[mod],mod);
 884     fEMCALMatrix[mod]->Print();
 885   }
 886
 887   return kTRUE;
 888 }
 889
 890 //_____________________________________________________
 891 Int_t AliEMCALTenderSupply::InitBadChannels()
 892 {
 893   // Initialising bad channel maps
 894
 895   AliVEvent *event = GetEvent();
 896
 897   if (!event)
 898     return 0;
 899
 900   if (fDebugLevel>0)
 901     AliInfo("Initialising Bad channel map");
 902
 903   // init default maps first
 904   if( !fEMCALRecoUtils->GetEMCALBadChannelStatusMapArray() )
 905     fEMCALRecoUtils->InitEMCALBadChannelStatusMap() ;
 906
 907   Int_t runBC = event->GetRunNumber();
 908
 909   AliOADBContainer *contBC = new AliOADBContainer("");
 910   if (fBasePath!="")
 911   { //if fBasePath specified in the ->SetBasePath()
 912     if (fDebugLevel>0) AliInfo(Form("Loading Bad Channels OADB from given path %s",fBasePath.Data()));
 913
 914     TFile *fbad=new TFile(Form("%s/EMCALBadChannels.root",fBasePath.Data()),"read");
 915     if (!fbad || fbad->IsZombie())
 916     {
 917       AliFatal(Form("EMCALBadChannels.root was not found in the path provided: %s",fBasePath.Data()));
 918       return 0;
 919     }
 920
 921     if (fbad) delete fbad;
 922
 923     contBC->InitFromFile(Form("%s/EMCALBadChannels.root",fBasePath.Data()),"AliEMCALBadChannels");
 924   }
 925   else
 926   { // Else choose the one in the $ALICE_ROOT directory
 927     if (fDebugLevel>0) AliInfo("Loading Bad Channels OADB from $ALICE_ROOT/OADB/EMCAL");
 928
 929     TFile *fbad=new TFile("$ALICE_ROOT/OADB/EMCAL/EMCALBadChannels.root","read");
 930     if (!fbad || fbad->IsZombie())
 931     {
 932       AliFatal("$ALICE_ROOT/OADB/EMCAL/EMCALBadChannels.root was not found");
 933       return 0;
 934     }
 935
 936     if (fbad) delete fbad;
 937
 938     contBC->InitFromFile("$ALICE_ROOT/OADB/EMCAL/EMCALBadChannels.root","AliEMCALBadChannels");
 939   }
 940
 941   TObjArray *arrayBC=(TObjArray*)contBC->GetObject(runBC);
 942   if (!arrayBC)
 943   {
 944     AliError(Form("No external hot channel set for run number: %d", runBC));
 945     return 2;
 946   }
 947
 948   Int_t sms = fEMCALGeo->GetEMCGeometry()->GetNumberOfSuperModules();
 949   for (Int_t i=0; i<sms; ++i)
 950   {
 951     TH2I *h = fEMCALRecoUtils->GetEMCALChannelStatusMap(i);
 952     if (h)
 953       delete h;
 954     h=(TH2I*)arrayBC->FindObject(Form("EMCALBadChannelMap_Mod%d",i));
 955
 956     if (!h)
 957     {
 958       AliError(Form("Can not get EMCALBadChannelMap_Mod%d",i));
 959       continue;
 960     }
 961     h->SetDirectory(0);
 962     fEMCALRecoUtils->SetEMCALChannelStatusMap(i,h);
 963   }
 964   return 1;
 965 }
 966
 967 //_____________________________________________________
 968 Int_t AliEMCALTenderSupply::InitRecalib()
 969 {
 970   // Initialising recalibration factors.
 971
 972   AliVEvent *event = GetEvent();
 973
 974   if (!event)
 975     return 0;
 976
 977   if (fDebugLevel>0)
 978     AliInfo("Initialising recalibration factors");
 979
 980   // init default maps first
 981   if( !fEMCALRecoUtils->GetEMCALRecalibrationFactorsArray() )
 982     fEMCALRecoUtils->InitEMCALRecalibrationFactors() ;
 983
 984   Int_t runRC = event->GetRunNumber();
 985
 986   AliOADBContainer *contRF=new AliOADBContainer("");
 987   if (fBasePath!="")
 988   { //if fBasePath specified in the ->SetBasePath()
 989     if (fDebugLevel>0)  AliInfo(Form("Loading Recalib OADB from given path %s",fBasePath.Data()));
 990
 991     TFile *fRecalib= new TFile(Form("%s/EMCALRecalib.root",fBasePath.Data()),"read");
 992     if (!fRecalib || fRecalib->IsZombie())
 993     {
 994       AliFatal(Form("EMCALRecalib.root not found in %s",fBasePath.Data()));
 995       return 0;
 996     }
 997
 998     if (fRecalib) delete fRecalib;
 999
1000     contRF->InitFromFile(Form("%s/EMCALRecalib.root",fBasePath.Data()),"AliEMCALRecalib");
1001   }
1002   else
1003   { // Else choose the one in the $ALICE_ROOT directory
1004     if (fDebugLevel>0)  AliInfo("Loading Recalib OADB from $ALICE_ROOT/OADB/EMCAL");
1005
1006     TFile *fRecalib= new TFile("$ALICE_ROOT/OADB/EMCAL/EMCALRecalib.root","read");
1007     if (!fRecalib || fRecalib->IsZombie())
1008     {
1009       AliFatal("$ALICE_ROOT/OADB/EMCAL/EMCALRecalib.root was not found");
1010       return 0;
1011     }
1012
1013     if (fRecalib) delete fRecalib;
1014
1015     contRF->InitFromFile("$ALICE_ROOT/OADB/EMCAL/EMCALRecalib.root","AliEMCALRecalib");
1016   }
1017
1018   TObjArray *recal=(TObjArray*)contRF->GetObject(runRC);
1019   if (!recal)
1020   {
1021     AliError(Form("No Objects for run: %d",runRC));
1022     return 2;
1023   }
1024
1025   TObjArray *recalpass=(TObjArray*)recal->FindObject(fFilepass);
1026   if (!recalpass)
1027   {
1028     AliError(Form("No Objects for run: %d - %s",runRC,fFilepass.Data()));
1029     return 2;
1030   }
1031
1032   TObjArray *recalib=(TObjArray*)recalpass->FindObject("Recalib");
1033   if (!recalib)
1034   {
1035     AliError(Form("No Recalib histos found for  %d - %s",runRC,fFilepass.Data()));
1036     return 2;
1037   }
1038
1039   if (fDebugLevel>0) recalib->Print();
1040
1041   Int_t sms = fEMCALGeo->GetEMCGeometry()->GetNumberOfSuperModules();
1042   for (Int_t i=0; i<sms; ++i)
1043   {
1044     TH2F *h = fEMCALRecoUtils->GetEMCALChannelRecalibrationFactors(i);
1045     if (h)
1046       delete h;
1047     h = (TH2F*)recalib->FindObject(Form("EMCALRecalFactors_SM%d",i));
1048     if (!h)
1049     {
1050       AliError(Form("Could not load EMCALRecalFactors_SM%d",i));
1051       continue;
1052     }
1053     h->SetDirectory(0);
1054     fEMCALRecoUtils->SetEMCALChannelRecalibrationFactors(i,h);
1055   }
1056   return 1;
1057 }
1058
1059 //_____________________________________________________
1060 Int_t AliEMCALTenderSupply::InitRunDepRecalib()
1061 {
1062   // Initialising recalibration factors.
1063
1064   AliVEvent *event = GetEvent();
1065
1066   if (!event)
1067     return 0;
1068
1069   if (fDebugLevel>0)
1070     AliInfo("Initialising recalibration factors");
1071
1072   // init default maps first
1073   if( !fEMCALRecoUtils->GetEMCALRecalibrationFactorsArray() )
1074     fEMCALRecoUtils->InitEMCALRecalibrationFactors() ;
1075
1076   Int_t runRC = event->GetRunNumber();
1077
1078   AliOADBContainer *contRF=new AliOADBContainer("");
1079   if (fBasePath!="")
1080   { //if fBasePath specified in the ->SetBasePath()
1081     if (fDebugLevel>0)  AliInfo(Form("Loading Recalib OADB from given path %s",fBasePath.Data()));
1082
1083     TFile *fRunDepRecalib= new TFile(Form("%s/EMCALTemperatureCorrCalib.root",fBasePath.Data()),"read");
1084     if (!fRunDepRecalib || fRunDepRecalib->IsZombie())
1085     {
1086       AliFatal(Form("EMCALTemperatureCorrCalib.root not found in %s",fBasePath.Data()));
1087       return 0;
1088     }
1089
1090     if (fRunDepRecalib) delete fRunDepRecalib;
1091
1092     contRF->InitFromFile(Form("%s/EMCALTemperatureCorrCalib.root",fBasePath.Data()),"AliEMCALRunDepTempCalibCorrections");
1093   }
1094   else
1095   { // Else choose the one in the $ALICE_ROOT directory
1096     if (fDebugLevel>0)  AliInfo("Loading Recalib OADB from $ALICE_ROOT/OADB/EMCAL");
1097
1098     TFile *fRunDepRecalib= new TFile("$ALICE_ROOT/OADB/EMCAL/EMCALTemperatureCorrCalib.root","read");
1099     if (!fRunDepRecalib || fRunDepRecalib->IsZombie())
1100     {
1101       AliFatal("$ALICE_ROOT/OADB/EMCAL/EMCALTemperatureCorrCalib.root was not found");
1102       return 0;
1103     }
1104
1105     if (fRunDepRecalib) delete fRunDepRecalib;
1106
1107     contRF->InitFromFile("$ALICE_ROOT/OADB/EMCAL/EMCALTemperatureCorrCalib.root","AliEMCALRunDepTempCalibCorrections");
1108   }
1109
1110   TH1S *rundeprecal=(TH1S*)contRF->GetObject(runRC);
1111   if (!rundeprecal)
1112   {
1113     AliWarning(Form("No TemperatureCorrCalib Objects for run: %d",runRC));
1114     // let's get the closest runnumber instead then..
1115     Int_t lower = 0;
1116     Int_t ic = 0;
1117     Int_t maxEntry = contRF->GetNumberOfEntries();
1118
1119     while ( (ic < maxEntry) && (contRF->UpperLimit(ic) < runRC) ) {
1120       lower = ic;
1121       ic++;
1122     }
1123
1124     Int_t closest = lower;
1125     if ( (ic<maxEntry) &&
1126          (contRF->LowerLimit(ic)-runRC) < (runRC - contRF->UpperLimit(lower)) ) {
1127          closest = ic;
1128     }
1129
1130     AliWarning(Form("TemperatureCorrCalib Objects found closest id %d from run: %d", closest, contRF->LowerLimit(closest)));
1131     rundeprecal = (TH1S*) contRF->GetObjectByIndex(closest);
1132   }
1133
1134   if (fDebugLevel>0) rundeprecal->Print();
1135
1136   Int_t nSM = fEMCALGeo->GetEMCGeometry()->GetNumberOfSuperModules();
1137
1138   for (Int_t ism=0; ism<nSM; ++ism)
1139   {
1140     for (Int_t icol=0; icol<48; ++icol)
1141     {
1142       for (Int_t irow=0; irow<24; ++irow)
1143       {
1144         Float_t factor = fEMCALRecoUtils->GetEMCALChannelRecalibrationFactor(ism,icol,irow);
1145
1146         Int_t absID = fEMCALGeo->GetAbsCellIdFromCellIndexes(ism, irow, icol); // original calibration factor
1147         factor *= rundeprecal->GetBinContent(absID) / 10000. ; // correction dependent on T
1148         //printf("\t ism %d, icol %d, irow %d,absID %d, corrA %2.3f, corrB %2.3f, corrAB %2.3f\n",ism, icol, irow, absID,
1149         //      GetEMCALChannelRecalibrationFactor(ism,icol,irow) , rundeprecal->GetBinContent(absID) / 10000., factor);
1150         fEMCALRecoUtils->SetEMCALChannelRecalibrationFactor(ism,icol,irow,factor);
1151       } // columns
1152     } // rows
1153   } // SM loop
1154
1155   return 1;
1156 }
1157
1158
1159 //_____________________________________________________
1160 Int_t AliEMCALTenderSupply::InitTimeCalibration()
1161 {
1162   // Initialising bad channel maps
1163   AliVEvent *event = GetEvent();
1164
1165   if (!event)
1166     return 0;
1167
1168   if (fDebugLevel>0)
1169     AliInfo("Initialising time calibration map");
1170
1171   // init default maps first
1172   if ( !fEMCALRecoUtils->GetEMCALTimeRecalibrationFactorsArray() )
1173     fEMCALRecoUtils->InitEMCALTimeRecalibrationFactors() ;
1174
1175   Int_t runBC = event->GetRunNumber();
1176
1177   AliOADBContainer *contBC = new AliOADBContainer("");
1178   if (fBasePath!="")
1179   { //if fBasePath specified in the ->SetBasePath()
1180     if (fDebugLevel>0) AliInfo(Form("Loading time calibration OADB from given path %s",fBasePath.Data()));
1181
1182     TFile *fbad=new TFile(Form("%s/EMCALTimeCalib.root",fBasePath.Data()),"read");
1183     if (!fbad || fbad->IsZombie())
1184     {
1185       AliFatal(Form("EMCALTimeCalib.root was not found in the path provided: %s",fBasePath.Data()));
1186       return 0;
1187     }
1188
1189     if (fbad) delete fbad;
1190
1191     contBC->InitFromFile(Form("%s/EMCALTimeCalib.root",fBasePath.Data()),"AliEMCALTimeCalib");
1192   }
1193   else
1194   { // Else choose the one in the $ALICE_ROOT directory
1195     if (fDebugLevel>0) AliInfo("Loading time calibration OADB from $ALICE_ROOT/OADB/EMCAL");
1196
1197     TFile *fbad=new TFile("$ALICE_ROOT/OADB/EMCAL/EMCALTimeCalib.root","read");
1198     if (!fbad || fbad->IsZombie())
1199     {
1200       AliFatal("$ALICE_ROOT/OADB/EMCAL/EMCALTimeCalib.root was not found");
1201       return 0;
1202     }
1203
1204     if (fbad) delete fbad;
1205
1206     contBC->InitFromFile("$ALICE_ROOT/OADB/EMCAL/EMCALTimeCalib.root","AliEMCALTimeCalib");
1207   }
1208
1209   TObjArray *arrayBC=(TObjArray*)contBC->GetObject(runBC);
1210   if (!arrayBC)
1211   {
1212     AliError(Form("No external time calibration set for run number: %d", runBC));
1213     return 2;
1214   }
1215
1216   // Here, it looks for a specific pass
1217   TObjArray *arrayBCpass=(TObjArray*)arrayBC->FindObject(fFilepass);
1218   if (!arrayBCpass)
1219   {
1220     AliError(Form("No external time calibration set for: %d -%s", runBC,fFilepass.Data()));
1221     return 2;
1222   }
1223
1224   if (fDebugLevel>0) arrayBCpass->Print();
1225
1226   for( Int_t i = 0; i < 4; i++ )
1227   {
1228     TH1F *h = fEMCALRecoUtils->GetEMCALChannelTimeRecalibrationFactors( i );
1229     if( h )
1230       delete h;
1231
1232     h = (TH1F*)arrayBCpass->FindObject(Form("hAllTimeAvBC%d",i));
1233
1234     if (!h)
1235     {
1236       AliError(Form("Can not get hAllTimeAvBC%d",i));
1237       continue;
1238     }
1239     h->SetDirectory(0);
1240     fEMCALRecoUtils->SetEMCALChannelTimeRecalibrationFactors(i,h);
1241   }
1242   return 1;
1243 }
1244
1245 //_____________________________________________________
1246 void AliEMCALTenderSupply::UpdateCells()
1247 {
1248   //Remove bad cells from the cell list
1249   //Recalibrate energy and time cells
1250   //This is required for later reclusterization
1251
1252   AliVEvent *event = GetEvent();
1253
1254   if(!event) return ;
1255
1256   AliVCaloCells *cells = event->GetEMCALCells();
1257   Int_t bunchCrossNo = event->GetBunchCrossNumber();
1258
1259   fEMCALRecoUtils->RecalibrateCells(cells, bunchCrossNo);
1260
1261   // remove exotic cells - loop through cells and zero the exotic ones
1262   // just like with bad cell rejection in reco utils (inside RecalibrateCells)
1263   if( fRejectExoticCells )
1264   {
1265     Short_t  absId  =-1;
1266     Double_t ecell = 0;
1267     Double_t tcell = 0;
1268     Double_t efrac = 0;
1269     Int_t  mclabel = -1;
1270     Bool_t   isExot = kFALSE;
1271
1272     // loop through cells
1273     Int_t nEMcell  = cells->GetNumberOfCells() ;
1274     for (Int_t iCell = 0; iCell < nEMcell; iCell++)
1275     {
1276       cells->GetCell( iCell, absId, ecell, tcell, mclabel, efrac );
1277
1278       isExot = fEMCALRecoUtils->IsExoticCell( absId, cells, bunchCrossNo );
1279       // zero if exotic
1280       if( isExot )
1281         cells->SetCell( iCell, absId, 0.0, -1.0, mclabel, efrac );
1282     } // cell loop
1283   } // reject exotic cells
1284
1285   cells->Sort();
1286 }
1287
1288 //_____________________________________________________
1289 TString AliEMCALTenderSupply::GetBeamType()
1290 {
1291   // Get beam type : pp-AA-pA
1292   // ESDs have it directly, AODs get it from hardcoded run number ranges
1293
1294   AliVEvent *event = GetEvent();
1295
1296   if (!event) {
1297     AliError("Couldn't retrieve event!");
1298     return "";
1299   }
1300
1301   TString beamType;
1302
1303   AliESDEvent *esd = dynamic_cast<AliESDEvent*>(event);
1304   if (esd) {
1305     const AliESDRun *run = esd->GetESDRun();
1306     beamType = run->GetBeamType();
1307   }
1308   else
1309   {
1310     Int_t runNumber = event->GetRunNumber();
1311     if ((runNumber >= 136851 && runNumber <= 139517)  // LHC10h
1312   || (runNumber >= 166529 && runNumber <= 170593))  // LHC11h
1313     {
1314       beamType = "A-A";
1315     }
1316     else
1317     {
1318       beamType = "p-p";
1319     }
1320   }
1321
1322   return beamType;
1323 }
1324
1325 //_____________________________________________________
1326 Int_t AliEMCALTenderSupply::InitRecParam()
1327 {
1328   // Init reco params if not yet exist (probably shipped by the user already)
1329
1330   if( fRecParam != 0 )
1331     return 2;
1332
1333   TString beamType = GetBeamType();
1334
1335   // set some default reco params
1336   fRecParam = new AliEMCALRecParam();
1337   fRecParam->SetClusteringThreshold(0.100);
1338   fRecParam->SetMinECut(0.050);
1339   fRecParam->SetTimeCut(250);
1340   fRecParam->SetTimeMin(425);
1341   fRecParam->SetTimeMax(825);
1342   if ( beamType == "A-A") {
1343     fRecParam->SetClusterizerFlag(AliEMCALRecParam::kClusterizerv2);
1344   }
1345   else {
1346     fRecParam->SetClusterizerFlag(AliEMCALRecParam::kClusterizerv1);
1347   }
1348
1349   return 0;
1350 }
1351
1352 //_____________________________________________________
1353 Bool_t AliEMCALTenderSupply::InitClusterization()
1354 {
1355   // Initialing clusterization/unfolding algorithm and set all the needed parameters.
1356
1357   AliVEvent *event = GetEvent();
1358
1359   if (!event)
1360     return kFALSE;
1361
1362   if (fDebugLevel>0)
1363     AliInfo(Form("Initialising reclustering parameters: Clusterizer type: %d",fRecParam->GetClusterizerFlag()));
1364
1365   //---setup clusterizer
1366   if (fClusterizer) {
1367     // avoid deleting fDigitsArr
1368     fClusterizer->SetDigitsArr(0);
1369     delete fClusterizer;
1370     fClusterizer = 0;
1371   }
1372   if     (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv1)
1373     fClusterizer = new AliEMCALClusterizerv1 (fEMCALGeo);
1374   else if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv2)
1375     fClusterizer = new AliEMCALClusterizerv2(fEMCALGeo);
1376   else if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerNxN)
1377   {
1378     AliEMCALClusterizerNxN *clusterizer = new AliEMCALClusterizerNxN(fEMCALGeo);
1379     clusterizer->SetNRowDiff(fRecParam->GetNRowDiff());
1380     clusterizer->SetNColDiff(fRecParam->GetNColDiff());
1381     fClusterizer = clusterizer;
1382   }
1383   else
1384   {
1385     AliFatal(Form("Clusterizer < %d > not available", fRecParam->GetClusterizerFlag()));
1386     return kFALSE;
1387   }
1388
1389   // Set the clustering parameters
1390   fClusterizer->SetECAClusteringThreshold( fRecParam->GetClusteringThreshold() );
1391   fClusterizer->SetECALogWeight          ( fRecParam->GetW0()                  );
1392   fClusterizer->SetMinECut               ( fRecParam->GetMinECut()             );
1393   fClusterizer->SetUnfolding             ( fRecParam->GetUnfold()              );
1394   fClusterizer->SetECALocalMaxCut        ( fRecParam->GetLocMaxCut()           );
1395   fClusterizer->SetTimeCut               ( fRecParam->GetTimeCut()             );
1396   fClusterizer->SetTimeMin               ( fRecParam->GetTimeMin()             );
1397   fClusterizer->SetTimeMax               ( fRecParam->GetTimeMax()             );
1398   fClusterizer->SetInputCalibrated       ( kTRUE                               );
1399   fClusterizer->SetJustClusters          ( kTRUE                               );
1400
1401   // In case of unfolding after clusterization is requested, set the corresponding parameters
1402   if (fRecParam->GetUnfold())
1403   {
1404     for (Int_t i = 0; i < 8; ++i)
1405     {
1406       fClusterizer->SetSSPars(i, fRecParam->GetSSPars(i));
1407     }
1408     for (Int_t i = 0; i < 3; ++i)
1409     {
1410       fClusterizer->SetPar5  (i, fRecParam->GetPar5(i));
1411       fClusterizer->SetPar6  (i, fRecParam->GetPar6(i));
1412     }
1413     fClusterizer->InitClusterUnfolding();
1414   }
1415
1416   fClusterizer->SetDigitsArr(fDigitsArr);
1417   fClusterizer->SetOutput(0);
1418   fClusterArr = const_cast<TObjArray *>(fClusterizer->GetRecPoints());
1419   return kTRUE;
1420 }
1421
1422 //_____________________________________________________
1423 void AliEMCALTenderSupply::FillDigitsArray()
1424 {
1425   // Fill digits from cells to a TClonesArray.
1426
1427   AliVEvent *event = GetEvent();
1428
1429  if (!event)
1430     return;
1431
1432   // In case of MC productions done before aliroot tag v5-02-Rev09
1433   // assing the cluster label to all the cells belonging to this cluster
1434   // very rough
1435   Int_t cellLabels[12672];
1436   if(fSetCellMCLabelFromCluster)
1437   {
1438     for (Int_t i = 0; i < 12672; i++) cellLabels[i] = 0;
1439
1440     Int_t nClusters = event->GetNumberOfCaloClusters();
1441     for (Int_t i = 0; i < nClusters; i++)
1442     {
1443       AliVCluster *clus =  event->GetCaloCluster(i);
1444
1445       if(!clus) continue;
1446
1447       if(!clus->IsEMCAL()) continue ;
1448
1449       Int_t      label = clus->GetLabel();
1450       UShort_t * index = clus->GetCellsAbsId() ;
1451
1452       for(Int_t icell=0; icell < clus->GetNCells(); icell++ )
1453         cellLabels[index[icell]] = label;
1454
1455     }// cluster loop
1456   }
1457
1458   fDigitsArr->Clear("C");
1459   AliVCaloCells *cells = event->GetEMCALCells();
1460   Int_t ncells = cells->GetNumberOfCells();
1461   for (Int_t icell = 0, idigit = 0; icell < ncells; ++icell)
1462   {
1463     Double_t cellAmplitude=0, cellTime=0, efrac = 0;
1464     Short_t  cellNumber=0;
1465     Int_t mcLabel=-1;
1466
1467     if (cells->GetCell(icell, cellNumber, cellAmplitude, cellTime, mcLabel, efrac) != kTRUE)
1468       break;
1469
1470     // Do not add if energy already too low (some cells set to 0 if bad channels)
1471     if (cellAmplitude < fRecParam->GetMinECut())
1472       continue;
1473
1474     // If requested, do not include exotic cells
1475    if (fEMCALRecoUtils->IsExoticCell(cellNumber,cells,event->GetBunchCrossNumber()))
1476       continue;
1477
1478     if(fSetCellMCLabelFromCluster) mcLabel = cellLabels[cellNumber];
1479
1480     if (mcLabel > 0 && efrac == 0) efrac = 1;
1481
1482     new((*fDigitsArr)[idigit]) AliEMCALDigit(mcLabel, mcLabel, cellNumber,
1483                                              (Float_t)cellAmplitude, (Float_t)cellTime,
1484                                              AliEMCALDigit::kHG,idigit, 0, 0, efrac*cellAmplitude);
1485     idigit++;
1486   }
1487 }
1488
1489 //_____________________________________________________
1490 void AliEMCALTenderSupply::Clusterize()
1491 {
1492   // Clusterize.
1493
1494   fClusterizer->Digits2Clusters("");
1495 }
1496
1497 //_____________________________________________________
1498 void AliEMCALTenderSupply::UpdateClusters()
1499 {
1500   // Update ESD cluster list.
1501
1502   AliVEvent *event = GetEvent();
1503
1504   if (!event)
1505     return;
1506
1507   TClonesArray *clus = dynamic_cast<TClonesArray*>(event->FindListObject("caloClusters"));
1508   if (!clus)
1509     clus = dynamic_cast<TClonesArray*>(event->FindListObject("CaloClusters"));
1510   if (!clus)
1511   {
1512     AliError(" Null pointer to calo clusters array, returning");
1513     return;
1514   }
1515
1516   Int_t nents = clus->GetEntriesFast();
1517   for (Int_t i=0; i < nents; ++i)
1518   {
1519     AliVCluster *c = dynamic_cast<AliVCluster*>(clus->At(i));
1520     if (!c)
1521       continue;
1522     if (c->IsEMCAL())
1523       delete clus->RemoveAt(i);
1524   }
1525
1526   clus->Compress();
1527
1528   RecPoints2Clusters(clus);
1529 }
1530
1531 //_____________________________________________________
1532 void AliEMCALTenderSupply::RecPoints2Clusters(TClonesArray *clus)
1533 {
1534   // Convert AliEMCALRecoPoints to AliESDCaloClusters/AliAODCaloClusters.
1535   // Cluster energy, global position, cells and their amplitude fractions are restored.
1536
1537   AliVEvent *event = GetEvent();
1538
1539   if (!event)
1540     return;
1541
1542   Int_t ncls = fClusterArr->GetEntriesFast();
1543   for(Int_t i=0, nout=clus->GetEntriesFast(); i < ncls; ++i)
1544   {
1545     AliEMCALRecPoint *recpoint = static_cast<AliEMCALRecPoint*>(fClusterArr->At(i));
1546
1547     Int_t ncellsTrue = 0;
1548     const Int_t ncells = recpoint->GetMultiplicity();
1549     UShort_t   absIds[ncells];
1550     Double32_t ratios[ncells];
1551     Int_t *dlist = recpoint->GetDigitsList();
1552     Float_t *elist = recpoint->GetEnergiesList();
1553     for (Int_t c = 0; c < ncells; ++c)
1554     {
1555       AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(dlist[c]));
1556       absIds[ncellsTrue] = digit->GetId();
1557       ratios[ncellsTrue] = elist[c]/digit->GetAmplitude();
1558       if (ratios[ncellsTrue] < 0.001)
1559         continue;
1560       ++ncellsTrue;
1561     }
1562
1563     if (ncellsTrue < 1)
1564     {
1565       AliWarning("Skipping cluster with no cells");
1566       continue;
1567     }
1568
1569     // calculate new cluster position
1570     TVector3 gpos;
1571     recpoint->GetGlobalPosition(gpos);
1572     Float_t g[3];
1573     gpos.GetXYZ(g);
1574
1575     AliVCluster *c = static_cast<AliVCluster*>(clus->New(nout++));
1576     c->SetID(nout-1);
1577     c->SetType(AliVCluster::kEMCALClusterv1);
1578     c->SetE(recpoint->GetEnergy());
1579     c->SetPosition(g);
1580     c->SetNCells(ncellsTrue);
1581     c->SetDispersion(recpoint->GetDispersion());
1582     c->SetEmcCpvDistance(-1);            //not yet implemented
1583     c->SetChi2(-1);                      //not yet implemented
1584     c->SetTOF(recpoint->GetTime()) ;     //time-of-flight
1585     c->SetNExMax(recpoint->GetNExMax()); //number of local maxima
1586     Float_t elipAxis[2];
1587     recpoint->GetElipsAxis(elipAxis);
1588     c->SetM02(elipAxis[0]*elipAxis[0]) ;
1589     c->SetM20(elipAxis[1]*elipAxis[1]) ;
1590     c->SetCellsAbsId(absIds);
1591     c->SetCellsAmplitudeFraction(ratios);
1592
1593     //MC
1594     AliESDCaloCluster *esdClus = dynamic_cast<AliESDCaloCluster*>(c);
1595     if (esdClus) {
1596       Int_t  parentMult = 0;
1597       Int_t *parentList = recpoint->GetParents(parentMult);
1598       if (parentMult > 0) {
1599         TArrayI parents(parentMult, parentList);
1600         esdClus->AddLabels(parents);
1601       }
1602     }
1603     else {
1604       AliAODCaloCluster *aodClus = dynamic_cast<AliAODCaloCluster*>(c);
1605       if (aodClus) {
1606         Int_t  parentMult = 0;
1607         Int_t *parentList = recpoint->GetParents(parentMult);
1608         aodClus->SetLabel(parentList, parentMult);
1609       }
1610     }
1611   }
1612 }
1613
1614 //_____________________________________________________
1615 void AliEMCALTenderSupply::GetPass()
1616 {
1617   // Get passx from filename
1618
1619   AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
1620   fInputTree = mgr->GetTree();
1621
1622   if (!fInputTree)
1623   {
1624     AliError("Pointer to tree = 0, returning");
1625     return;
1626   }
1627
1628   fInputFile = fInputTree->GetCurrentFile();
1629   if (!fInputFile) {
1630     AliError("Null pointer input file, returning");
1631     return;
1632   }
1633
1634   TString fname(fInputFile->GetName());
1635   if      (fname.Contains("pass1")) fFilepass = TString("pass1");
1636   else if (fname.Contains("pass2")) fFilepass = TString("pass2");
1637   else if (fname.Contains("pass3")) fFilepass = TString("pass3");
1638   else if (fname.Contains("pass4")) fFilepass = TString("pass4");
1639   else
1640   {
1641     AliError(Form("Pass number string not found: %s", fname.Data()));
1642     return;
1643   }
1644 }