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