PWG/CaloTrackCorrBase/AliCalorimeterUtils.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 // Class utility for Calorimeter specific selection methods                ///
  18 //
  19 //
  20 //
  21 //-- Author: Gustavo Conesa (LPSC-Grenoble)
  22 //////////////////////////////////////////////////////////////////////////////
  23
  24
  25 // --- ROOT system ---
  26 #include <TGeoManager.h>
  27 #include <TH2F.h>
  28 #include <TCanvas.h>
  29 #include <TStyle.h>
  30 #include <TPad.h>
  31 #include <TFile.h>
  32 #include <TParticle.h>
  33
  34 // --- ANALYSIS system ---
  35 #include "AliCalorimeterUtils.h"
  36 #include "AliESDEvent.h"
  37 #include "AliMCEvent.h"
  38 #include "AliStack.h"
  39 #include "AliAODPWG4Particle.h"
  40 #include "AliVCluster.h"
  41 #include "AliVCaloCells.h"
  42 #include "AliMixedEvent.h"
  43 #include "AliAODCaloCluster.h"
  44 #include "AliOADBContainer.h"
  45 #include "AliAnalysisManager.h"
  46 #include "AliAODMCParticle.h"
  47
  48 // --- Detector ---
  49 #include "AliEMCALGeometry.h"
  50 #include "AliPHOSGeoUtils.h"
  51
  52 ClassImp(AliCalorimeterUtils)
  53
  54
  55 //____________________________________________
  56   AliCalorimeterUtils::AliCalorimeterUtils() :
  57     TObject(), fDebug(0),
  58     fEMCALGeoName(""),
  59     fPHOSGeoName (""),
  60     fEMCALGeo(0x0),                   fPHOSGeo(0x0),
  61     fEMCALGeoMatrixSet(kFALSE),       fPHOSGeoMatrixSet(kFALSE),
  62     fLoadEMCALMatrices(kFALSE),       fLoadPHOSMatrices(kFALSE),
  63     fRemoveBadChannels(kFALSE),       fPHOSBadChannelMap(0x0),
  64     fNCellsFromPHOSBorder(0),
  65     fNMaskCellColumns(0),             fMaskCellColumns(0x0),
  66     fRecalibration(kFALSE),           fRunDependentCorrection(kFALSE),
  67     fPHOSRecalibrationFactors(),      fEMCALRecoUtils(new AliEMCALRecoUtils),
  68     fRecalculatePosition(kFALSE),     fCorrectELinearity(kFALSE),
  69     fRecalculateMatching(kFALSE),
  70     fCutR(20),                        fCutZ(20),
  71     fCutEta(20),                      fCutPhi(20),
  72     fLocMaxCutE(0),                   fLocMaxCutEDiff(0),
  73     fPlotCluster(0),                  fOADBSet(kFALSE),
  74     fOADBForEMCAL(kFALSE),            fOADBForPHOS(kFALSE),
  75     fOADBFilePathEMCAL(""),           fOADBFilePathPHOS(""),
  76     fImportGeometryFromFile(0),       fImportGeometryFilePath(""),
  77     fNSuperModulesUsed(0),
  78     fMCECellClusFracCorrOn(0),        fMCECellClusFracCorrParam()
  79 {
  80   //Ctor
  81
  82   //Initialize parameters
  83   InitParameters();
  84   for(Int_t i = 0; i < 12; i++) fEMCALMatrix[i] = 0 ;
  85   for(Int_t i = 0; i < 5 ; i++) fPHOSMatrix [i] = 0 ;
  86
  87 }
  88
  89 //_________________________________________
  90 AliCalorimeterUtils::~AliCalorimeterUtils()
  91 {
  92   //Dtor
  93
  94   //if(fPHOSGeo)  delete fPHOSGeo  ;
  95   if(fEMCALGeo) delete fEMCALGeo ;
  96
  97   if(fPHOSBadChannelMap) {
  98     fPHOSBadChannelMap->Clear();
  99     delete  fPHOSBadChannelMap;
 100   }
 101
 102   if(fPHOSRecalibrationFactors) {
 103     fPHOSRecalibrationFactors->Clear();
 104     delete  fPHOSRecalibrationFactors;
 105   }
 106
 107   if(fEMCALRecoUtils)   delete fEMCALRecoUtils ;
 108   if(fNMaskCellColumns) delete [] fMaskCellColumns;
 109
 110 }
 111
 112 //____________________________________________________
 113 void AliCalorimeterUtils::AccessOADB(AliVEvent* event)
 114 {
 115   // Set the AODB calibration, bad channels etc. parameters at least once
 116   // alignment matrices from OADB done in SetGeometryMatrices
 117
 118   //Set it only once
 119   if(fOADBSet) return ;
 120
 121   Int_t   runnumber = event->GetRunNumber() ;
 122   TString pass      = GetPass();
 123
 124   // EMCAL
 125   if(fOADBForEMCAL)
 126   {
 127     printf("AliCalorimeterUtils::SetOADBParameters() - Get AODB parameters from EMCAL in %s for run %d, and <%s> \n",fOADBFilePathEMCAL.Data(),runnumber,pass.Data());
 128
 129     Int_t nSM = fEMCALGeo->GetNumberOfSuperModules();
 130
 131     // Bad map
 132     if(fRemoveBadChannels)
 133     {
 134       AliOADBContainer *contBC=new AliOADBContainer("");
 135       contBC->InitFromFile(Form("%s/EMCALBadChannels.root",fOADBFilePathEMCAL.Data()),"AliEMCALBadChannels");
 136
 137       TObjArray *arrayBC=(TObjArray*)contBC->GetObject(runnumber);
 138
 139       if(arrayBC)
 140       {
 141         SwitchOnDistToBadChannelRecalculation();
 142         printf("AliCalorimeterUtils::SetOADBParameters() - Remove EMCAL bad cells \n");
 143
 144         for (Int_t i=0; i<nSM; ++i)
 145         {
 146           TH2I *hbm = GetEMCALChannelStatusMap(i);
 147
 148           if (hbm)
 149             delete hbm;
 150
 151           hbm=(TH2I*)arrayBC->FindObject(Form("EMCALBadChannelMap_Mod%d",i));
 152
 153           if (!hbm)
 154           {
 155             AliError(Form("Can not get EMCALBadChannelMap_Mod%d",i));
 156             continue;
 157           }
 158
 159           hbm->SetDirectory(0);
 160           SetEMCALChannelStatusMap(i,hbm);
 161
 162         } // loop
 163       } else printf("AliCalorimeterUtils::SetOADBParameters() - Do NOT remove EMCAL bad channels\n"); // run array
 164     }  // Remove bad
 165
 166     // Energy Recalibration
 167     if(fRecalibration)
 168     {
 169       AliOADBContainer *contRF=new AliOADBContainer("");
 170
 171       contRF->InitFromFile(Form("%s/EMCALRecalib.root",fOADBFilePathEMCAL.Data()),"AliEMCALRecalib");
 172
 173       TObjArray *recal=(TObjArray*)contRF->GetObject(runnumber);
 174
 175       if(recal)
 176       {
 177         TObjArray *recalpass=(TObjArray*)recal->FindObject(pass);
 178
 179         if(recalpass)
 180         {
 181           TObjArray *recalib=(TObjArray*)recalpass->FindObject("Recalib");
 182
 183           if(recalib)
 184           {
 185             printf("AliCalorimeterUtils::SetOADBParameters() - Recalibrate EMCAL \n");
 186             for (Int_t i=0; i<nSM; ++i)
 187             {
 188               TH2F *h = GetEMCALChannelRecalibrationFactors(i);
 189
 190               if (h)
 191                 delete h;
 192
 193               h = (TH2F*)recalib->FindObject(Form("EMCALRecalFactors_SM%d",i));
 194
 195               if (!h)
 196               {
 197                 AliError(Form("Could not load EMCALRecalFactors_SM%d",i));
 198                 continue;
 199               }
 200
 201               h->SetDirectory(0);
 202
 203               SetEMCALChannelRecalibrationFactors(i,h);
 204             } // SM loop
 205           }else printf("AliCalorimeterUtils::SetOADBParameters() - Do NOT recalibrate EMCAL, no params object array\n"); // array ok
 206         }else printf("AliCalorimeterUtils::SetOADBParameters() - Do NOT recalibrate EMCAL, no params for pass\n"); // array pass ok
 207       }else printf("AliCalorimeterUtils::SetOADBParameters() - Do NOT recalibrate EMCAL, no params for run\n");  // run number array ok
 208
 209       // once set, apply run dependent corrections if requested
 210       //fEMCALRecoUtils->SetRunDependentCorrections(runnumber);
 211
 212     } // Recalibration on
 213
 214     // Energy Recalibration, apply on top of previous calibration factors
 215     if(fRunDependentCorrection)
 216     {
 217       AliOADBContainer *contRFTD=new AliOADBContainer("");
 218
 219       contRFTD->InitFromFile(Form("%s/EMCALTemperatureCorrCalib.root",fOADBFilePathEMCAL.Data()),"AliEMCALRunDepTempCalibCorrections");
 220
 221       TH1S *htd=(TH1S*)contRFTD->GetObject(runnumber);
 222
 223       //If it did not exist for this run, get closes one
 224       if (!htd)
 225       {
 226         AliWarning(Form("No TemperatureCorrCalib Objects for run: %d",runnumber));
 227         // let's get the closest runnumber instead then..
 228         Int_t lower = 0;
 229         Int_t ic = 0;
 230         Int_t maxEntry = contRFTD->GetNumberOfEntries();
 231
 232         while ( (ic < maxEntry) && (contRFTD->UpperLimit(ic) < runnumber) ) {
 233           lower = ic;
 234           ic++;
 235         }
 236
 237         Int_t closest = lower;
 238         if ( (ic<maxEntry) &&
 239             (contRFTD->LowerLimit(ic)-runnumber) < (runnumber - contRFTD->UpperLimit(lower)) ) {
 240           closest = ic;
 241         }
 242
 243         AliWarning(Form("TemperatureCorrCalib Objects found closest id %d from run: %d", closest, contRFTD->LowerLimit(closest)));
 244         htd = (TH1S*) contRFTD->GetObjectByIndex(closest);
 245       }
 246
 247       if(htd)
 248       {
 249         printf("AliCalorimeterUtils::SetOADBParameters() - Recalibrate (Temperature) EMCAL \n");
 250
 251         for (Int_t ism=0; ism<nSM; ++ism)
 252         {
 253           for (Int_t icol=0; icol<48; ++icol)
 254           {
 255             for (Int_t irow=0; irow<24; ++irow)
 256             {
 257               Float_t factor = GetEMCALChannelRecalibrationFactor(ism,icol,irow);
 258
 259               Int_t absID = fEMCALGeo->GetAbsCellIdFromCellIndexes(ism, irow, icol); // original calibration factor
 260               factor *= htd->GetBinContent(absID) / 10000. ; // correction dependent on T
 261               //printf("\t ism %d, icol %d, irow %d,absID %d, corrA %2.3f, corrB %2.3f, corrAB %2.3f\n",ism, icol, irow, absID,
 262               //      GetEMCALChannelRecalibrationFactor(ism,icol,irow) , htd->GetBinContent(absID) / 10000., factor);
 263               SetEMCALChannelRecalibrationFactor(ism,icol,irow,factor);
 264             } // columns
 265           } // rows
 266         } // SM loop
 267       }else printf("AliCalorimeterUtils::SetOADBParameters() - Do NOT recalibrate EMCAL with T variations, no params TH1 \n");
 268     } // Run by Run T calibration
 269
 270     // Time Recalibration
 271     if(fEMCALRecoUtils->IsTimeRecalibrationOn())
 272     {
 273       AliOADBContainer *contTRF=new AliOADBContainer("");
 274
 275       contTRF->InitFromFile(Form("%s/EMCALTimeCalib.root",fOADBFilePathEMCAL.Data()),"AliEMCALTimeCalib");
 276
 277       TObjArray *trecal=(TObjArray*)contTRF->GetObject(runnumber);
 278
 279       if(trecal)
 280       {
 281         TString passM = pass;
 282         if(pass=="spc_calo") passM = "pass1";
 283         TObjArray *trecalpass=(TObjArray*)trecal->FindObject(passM);
 284
 285         if(trecalpass)
 286         {
 287           printf("AliCalorimeterUtils::SetOADBParameters() - Time Recalibrate EMCAL \n");
 288           for (Int_t ibc = 0; ibc < 4; ++ibc)
 289           {
 290             TH1F *h = GetEMCALChannelTimeRecalibrationFactors(ibc);
 291
 292             if (h)
 293               delete h;
 294
 295             h = (TH1F*)trecalpass->FindObject(Form("hAllTimeAvBC%d",ibc));
 296
 297             if (!h)
 298             {
 299               AliError(Form("Could not load hAllTimeAvBC%d",ibc));
 300               continue;
 301             }
 302
 303             h->SetDirectory(0);
 304
 305             SetEMCALChannelTimeRecalibrationFactors(ibc,h);
 306           } // bunch crossing loop
 307         }else printf("AliCalorimeterUtils::SetOADBParameters() - Do NOT recalibrate time EMCAL, no params for pass\n"); // array pass ok
 308       }else printf("AliCalorimeterUtils::SetOADBParameters() - Do NOT recalibrate time EMCAL, no params for run\n");  // run number array ok
 309
 310     } // Recalibration on
 311
 312   }// EMCAL
 313
 314   // PHOS
 315   if(fOADBForPHOS)
 316   {
 317     printf("AliCalorimeterUtils::SetOADBParameters() - Get AODB parameters from PHOS in %s for run %d, and <%s> \n",fOADBFilePathPHOS.Data(),runnumber,pass.Data());
 318
 319     // Bad map
 320     if(fRemoveBadChannels)
 321     {
 322       AliOADBContainer badmapContainer(Form("phosBadMap"));
 323       TString fileName="$ALICE_ROOT/OADB/PHOS/PHOSBadMaps.root";
 324       badmapContainer.InitFromFile(Form("%s/PHOSBadMaps.root",fOADBFilePathPHOS.Data()),"phosBadMap");
 325
 326       //Use a fixed run number from year 2010, this year not available yet.
 327       TObjArray *maps = (TObjArray*)badmapContainer.GetObject(139000,"phosBadMap");
 328       if(!maps)
 329       {
 330         printf("AliCalorimeterUtils::SetOADBParameters() - Can not read PHOS bad map for run %d.\n",runnumber) ;
 331       }
 332       else
 333       {
 334         printf("AliCalorimeterUtils::SetOADBParameters() - Setting PHOS bad map with name %s \n",maps->GetName()) ;
 335         for(Int_t mod=1; mod<5;mod++)
 336         {
 337           TH2I *hbmPH = GetPHOSChannelStatusMap(mod);
 338
 339           if(hbmPH)
 340             delete hbmPH ;
 341
 342           hbmPH = (TH2I*)maps->At(mod);
 343
 344           if(hbmPH) hbmPH->SetDirectory(0);
 345
 346           SetPHOSChannelStatusMap(mod-1,hbmPH);
 347
 348         } // modules loop
 349       } // maps exist
 350     } // Remove bad channels
 351   } // PHOS
 352
 353   // Parameters already set once, so do not it again
 354   fOADBSet = kTRUE;
 355
 356 }
 357
 358 //_____________________________________________________________
 359 void AliCalorimeterUtils::AccessGeometry(AliVEvent* inputEvent)
 360 {
 361   //Set the calorimeters transformation matrices and init geometry
 362
 363   // First init the geometry, a priory not done before
 364   Int_t runnumber = inputEvent->GetRunNumber() ;
 365   InitPHOSGeometry (runnumber);
 366   InitEMCALGeometry(runnumber);
 367
 368   //Get the EMCAL transformation geometry matrices from ESD
 369   if(!fEMCALGeoMatrixSet && fEMCALGeo)
 370   {
 371     if(fLoadEMCALMatrices)
 372     {
 373       printf("AliCalorimeterUtils::AccessGeometry() - Load user defined EMCAL geometry matrices\n");
 374
 375       // OADB if available
 376       AliOADBContainer emcGeoMat("AliEMCALgeo");
 377       emcGeoMat.InitFromFile(Form("%s/EMCALlocal2master.root",fOADBFilePathEMCAL.Data()),"AliEMCALgeo");
 378       TObjArray *matEMCAL=(TObjArray*)emcGeoMat.GetObject(runnumber,"EmcalMatrices");
 379
 380       for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++)
 381       {
 382         if (!fEMCALMatrix[mod]) // Get it from OADB
 383         {
 384           if(fDebug > 1 )
 385             printf("AliCalorimeterUtils::AccessGeometry() - EMCAL matrices SM %d, %p\n",
 386                    mod,((TGeoHMatrix*) matEMCAL->At(mod)));
 387           //((TGeoHMatrix*) matEMCAL->At(mod))->Print();
 388
 389           fEMCALMatrix[mod] = (TGeoHMatrix*) matEMCAL->At(mod) ;
 390         }
 391
 392         if(fEMCALMatrix[mod])
 393         {
 394           if(fDebug > 1)
 395             fEMCALMatrix[mod]->Print();
 396
 397           fEMCALGeo->SetMisalMatrix(fEMCALMatrix[mod],mod) ;
 398         }
 399       }//SM loop
 400
 401       fEMCALGeoMatrixSet = kTRUE;//At least one, so good
 402
 403     }//Load matrices
 404     else if (!gGeoManager)
 405     {
 406       if(fDebug > 1)
 407         printf(" AliCalorimeterUtils::AccessGeometry() - Load EMCAL misalignment matrices. \n");
 408       if(!strcmp(inputEvent->GetName(),"AliESDEvent"))
 409       {
 410         for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++)
 411         {
 412           //printf("Load ESD matrix %d, %p\n",mod,((AliESDEvent*)inputEvent)->GetEMCALMatrix(mod));
 413           if(((AliESDEvent*)inputEvent)->GetEMCALMatrix(mod))
 414           {
 415             fEMCALGeo->SetMisalMatrix(((AliESDEvent*)inputEvent)->GetEMCALMatrix(mod),mod) ;
 416           }
 417         }// loop over super modules
 418
 419         fEMCALGeoMatrixSet = kTRUE;//At least one, so good
 420
 421       }//ESD as input
 422       else
 423       {
 424         if(fDebug > 1)
 425           printf("AliCalorimeterUtils::SetGeometryTransformationMatrices() - Setting of EMCAL transformation matrixes for AODs not implemented yet. \n Import geometry.root file\n");
 426       }//AOD as input
 427     }//Get matrix from data
 428     else if(gGeoManager)
 429     {
 430       fEMCALGeoMatrixSet = kTRUE;
 431     }
 432   }//EMCAL geo && no geoManager
 433
 434         //Get the PHOS transformation geometry matrices from ESD
 435   if(!fPHOSGeoMatrixSet && fPHOSGeo)
 436   {
 437     if(fLoadPHOSMatrices)
 438     {
 439       printf("AliCalorimeterUtils::SetGeometryTransformationMatrices() - Load user defined PHOS geometry matrices\n");
 440
 441       // OADB if available
 442       AliOADBContainer geomContainer("phosGeo");
 443       geomContainer.InitFromFile(Form("%s/PHOSGeometry.root",fOADBFilePathPHOS.Data()),"PHOSRotationMatrixes");
 444       TObjArray *matPHOS = (TObjArray*)geomContainer.GetObject(139000,"PHOSRotationMatrixes");
 445
 446       for(Int_t mod = 0 ; mod < 5 ; mod++)
 447       {
 448         if (!fPHOSMatrix[mod]) // Get it from OADB
 449         {
 450           if(fDebug > 1 )
 451             printf("AliCalorimeterUtils::SetGeometryTransformationMatrices() - PHOS matrices module %d, %p\n",
 452                    mod,((TGeoHMatrix*) matPHOS->At(mod)));
 453           //((TGeoHMatrix*) matPHOS->At(mod))->Print();
 454
 455           fPHOSMatrix[mod] = (TGeoHMatrix*) matPHOS->At(mod) ;
 456         }
 457
 458         // Set it, if it exists
 459         if(fPHOSMatrix[mod])
 460         {
 461           if(fDebug > 1 )
 462             fPHOSMatrix[mod]->Print();
 463
 464           fPHOSGeo->SetMisalMatrix(fPHOSMatrix[mod],mod) ;
 465         }
 466       }// SM loop
 467
 468       fPHOSGeoMatrixSet = kTRUE;//At least one, so good
 469
 470     }//Load matrices
 471     else if (!gGeoManager)
 472     {
 473       if(fDebug > 1)
 474         printf(" AliCalorimeterUtils::SetGeometryTransformationMatrices() - Load PHOS misalignment matrices. \n");
 475                         if(!strcmp(inputEvent->GetName(),"AliESDEvent"))
 476       {
 477                                 for(Int_t mod = 0; mod < 5; mod++)
 478         {
 479                                         if( ((AliESDEvent*)inputEvent)->GetPHOSMatrix(mod))
 480           {
 481                                                 //printf("PHOS: mod %d, matrix %p\n",mod, ((AliESDEvent*)inputEvent)->GetPHOSMatrix(mod));
 482                                                 fPHOSGeo->SetMisalMatrix( ((AliESDEvent*)inputEvent)->GetPHOSMatrix(mod),mod) ;
 483                                         }
 484                                 }// loop over modules
 485         fPHOSGeoMatrixSet  = kTRUE; //At least one so good
 486                         }//ESD as input
 487                         else
 488       {
 489                                 if(fDebug > 1)
 490                                         printf("AliCalorimeterUtils::SetGeometryTransformationMatrices() - Setting of EMCAL transformation matrixes for AODs not implemented yet. \n Import geometry.root file\n");
 491       }//AOD as input
 492     }// get matrix from data
 493     else if(gGeoManager)
 494     {
 495       fPHOSGeoMatrixSet = kTRUE;
 496     }
 497         }//PHOS geo     and  geoManager was not set
 498
 499 }
 500
 501 //________________________________________________________________________________________
 502 Bool_t AliCalorimeterUtils::AreNeighbours(TString calo, Int_t absId1, Int_t absId2 ) const
 503 {
 504   // Tells if (true) or not (false) two cells are neighbours
 505   // A neighbour is defined as being two cells which share a side or corner
 506
 507   Bool_t areNeighbours = kFALSE ;
 508
 509   Int_t iRCU1 = -1, irow1 = -1, icol1 = -1;
 510   Int_t iRCU2 = -1, irow2 = -1, icol2 = -1;
 511
 512   Int_t rowdiff =  0, coldiff =  0;
 513
 514   Int_t nSupMod1 = GetModuleNumberCellIndexes(absId1, calo, icol1, irow1, iRCU1);
 515   Int_t nSupMod2 = GetModuleNumberCellIndexes(absId2, calo, icol2, irow2, iRCU2);
 516
 517   if(calo=="EMCAL" && nSupMod1!=nSupMod2)
 518   {
 519     // In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2-1
 520     // C Side impair SM, nSupMod%2=1; A side pair SM nSupMod%2=0
 521     if(nSupMod1%2) icol1+=AliEMCALGeoParams::fgkEMCALCols;
 522     else           icol2+=AliEMCALGeoParams::fgkEMCALCols;
 523         }
 524
 525   rowdiff = TMath::Abs( irow1 - irow2 ) ;
 526   coldiff = TMath::Abs( icol1 - icol2 ) ;
 527
 528   //if (( coldiff <= 1 )  && ( rowdiff <= 1 ) && (coldiff + rowdiff > 0))
 529   if ((coldiff + rowdiff == 1 ))
 530     areNeighbours = kTRUE ;
 531
 532   return areNeighbours;
 533
 534 }
 535 //_____________________________________________________________________________________
 536 Bool_t AliCalorimeterUtils::IsClusterSharedByTwoSuperModules(const AliEMCALGeometry * geom,
 537                                                              AliVCluster* cluster)
 538 {
 539   //Method that checks if any of the cells in the cluster belongs to a different SM
 540
 541   Int_t    iSupMod = -1;
 542   Int_t    iSM0    = -1;
 543   Int_t    iTower  = -1;
 544   Int_t    iIphi   = -1;
 545   Int_t    iIeta   = -1;
 546   Int_t    iphi    = -1;
 547   Int_t    ieta    = -1;
 548
 549   for(Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++)
 550   {
 551     //Get from the absid the supermodule, tower and eta/phi numbers
 552     geom->GetCellIndex(cluster->GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta);
 553     geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
 554
 555     //Check if there are cells of different SM
 556     if     (iDigit == 0   ) iSM0 = iSupMod;
 557     else if(iSupMod!= iSM0) return kTRUE;
 558   }
 559
 560   return kFALSE;
 561
 562 }
 563
 564 //_____________________________________________________________________________________
 565 Bool_t AliCalorimeterUtils::CheckCellFiducialRegion(AliVCluster* cluster,
 566                                                     AliVCaloCells* cells,
 567                                                     AliVEvent * event, Int_t iev) const
 568 {
 569
 570         // Given the list of AbsId of the cluster, get the maximum cell and
 571         // check if there are fNCellsFromBorder from the calorimeter border
 572
 573   //If the distance to the border is 0 or negative just exit accept all clusters
 574         if(cells->GetType()==AliVCaloCells::kEMCALCell && fEMCALRecoUtils->GetNumberOfCellsFromEMCALBorder() <= 0 ) return kTRUE;
 575         if(cells->GetType()==AliVCaloCells::kPHOSCell  && fNCellsFromPHOSBorder  <= 0 ) return kTRUE;
 576
 577   Int_t absIdMax        = -1;
 578         Float_t ampMax  = -1;
 579
 580   AliMixedEvent * mixEvent = dynamic_cast<AliMixedEvent*> (event);
 581   Int_t nMixedEvents = 0 ;
 582   Int_t * cellsCumul = NULL ;
 583   Int_t numberOfCells = 0 ;
 584   if (mixEvent){
 585     nMixedEvents = mixEvent->GetNumberOfEvents() ;
 586     if (cells->GetType()==AliVCaloCells::kEMCALCell) {
 587       cellsCumul =  mixEvent->GetEMCALCellsCumul() ;
 588       numberOfCells = mixEvent->GetNumberOfEMCALCells() ;
 589     }
 590
 591     else if (cells->GetType()==AliVCaloCells::kPHOSCell) {
 592       cellsCumul =  mixEvent->GetPHOSCellsCumul() ;
 593       numberOfCells = mixEvent->GetNumberOfPHOSCells() ;
 594     }
 595
 596     if(cellsCumul){
 597
 598       Int_t startCell = cellsCumul[iev] ;
 599       Int_t endCell   = (iev+1 < nMixedEvents)?cellsCumul[iev+1]:numberOfCells;
 600       //Find cells with maximum amplitude
 601       for(Int_t i = 0; i < cluster->GetNCells() ; i++){
 602         Int_t absId = cluster->GetCellAbsId(i) ;
 603         for (Int_t j = startCell; j < endCell ;  j++) {
 604           Short_t cellNumber;
 605           Int_t mclabel;
 606           Double_t amp, time, efrac;
 607           cells->GetCell(j, cellNumber, amp, time,mclabel,efrac) ;
 608           if (absId == cellNumber) {
 609             if(amp > ampMax){
 610               ampMax   = amp;
 611               absIdMax = absId;
 612             }
 613           }
 614         }
 615       }//loop on cluster cells
 616     }// cells cumul available
 617     else {
 618       printf("AliCalorimeterUtils::CheckCellFiducialRegion() - CellsCumul is NULL!!!\n");
 619       abort();
 620     }
 621   } else {//Normal SE Events
 622     for(Int_t i = 0; i < cluster->GetNCells() ; i++){
 623       Int_t absId = cluster->GetCellAbsId(i) ;
 624       Float_t amp       = cells->GetCellAmplitude(absId);
 625       if(amp > ampMax){
 626         ampMax   = amp;
 627         absIdMax = absId;
 628       }
 629     }
 630   }
 631
 632         if(fDebug > 1)
 633                 printf("AliCalorimeterUtils::CheckCellFiducialRegion() - Cluster Max AbsId %d, Cell Energy %2.2f, Cluster Energy %2.2f\n",
 634            absIdMax, ampMax, cluster->E());
 635
 636         if(absIdMax==-1) return kFALSE;
 637
 638         //Check if the cell is close to the borders:
 639         Bool_t okrow = kFALSE;
 640         Bool_t okcol = kFALSE;
 641
 642         if(cells->GetType()==AliVCaloCells::kEMCALCell){
 643
 644                 Int_t iTower = -1, iIphi = -1, iIeta = -1, iphi = -1, ieta = -1, iSM = -1;
 645                 fEMCALGeo->GetCellIndex(absIdMax,iSM,iTower,iIphi,iIeta);
 646                 fEMCALGeo->GetCellPhiEtaIndexInSModule(iSM,iTower,iIphi, iIeta,iphi,ieta);
 647                 if(iSM < 0 || iphi < 0 || ieta < 0 ) {
 648       Fatal("CheckCellFidutialRegion","Negative value for super module: %d, or cell ieta: %d, or cell iphi: %d, check EMCAL geometry name\n",iSM,ieta,iphi);
 649     }
 650
 651                 //Check rows/phi
 652     Int_t nborder = fEMCALRecoUtils->GetNumberOfCellsFromEMCALBorder();
 653                 if(iSM < 10)
 654     {
 655                         if(iphi >= nborder && iphi < 24-nborder) okrow =kTRUE;
 656     }
 657                 else
 658     {
 659       if(fEMCALGeoName.Contains("12SM")) // 1/3 SM
 660       {
 661         if(iphi >= nborder && iphi < 8-nborder) okrow =kTRUE;
 662       }
 663       else // 1/2 SM
 664       {
 665         if(iphi >= nborder && iphi <12-nborder) okrow =kTRUE;
 666       }
 667                 }
 668
 669                 //Check columns/eta
 670                 if(!fEMCALRecoUtils->IsEMCALNoBorderAtEta0())
 671     {
 672                         if(ieta  > nborder && ieta < 48-nborder) okcol =kTRUE;
 673                 }
 674                 else
 675     {
 676                         if(iSM%2==0)
 677       {
 678                                 if(ieta >= nborder)     okcol = kTRUE;
 679                         }
 680                         else
 681       {
 682                                 if(ieta <  48-nborder)  okcol = kTRUE;
 683                         }
 684                 }//eta 0 not checked
 685
 686                 if(fDebug > 1)
 687                 {
 688                         printf("AliCalorimeterUtils::CheckCellFiducialRegion() - EMCAL Cluster in %d cells fiducial volume: ieta %d, iphi %d, SM %d ?",
 689              nborder, ieta, iphi, iSM);
 690                         if (okcol && okrow ) printf(" YES \n");
 691                         else  printf(" NO: column ok? %d, row ok? %d \n",okcol,okrow);
 692                 }
 693         }//EMCAL
 694         else if(cells->GetType()==AliVCaloCells::kPHOSCell){
 695                 Int_t relId[4];
 696                 Int_t irow = -1, icol = -1;
 697                 fPHOSGeo->AbsToRelNumbering(absIdMax,relId);
 698                 irow = relId[2];
 699                 icol = relId[3];
 700                 //imod = relId[0]-1;
 701                 if(irow >= fNCellsFromPHOSBorder && irow < 64-fNCellsFromPHOSBorder) okrow =kTRUE;
 702                 if(icol >= fNCellsFromPHOSBorder && icol < 56-fNCellsFromPHOSBorder) okcol =kTRUE;
 703                 if(fDebug > 1)
 704                 {
 705                         printf("AliCalorimeterUtils::CheckCellFiducialRegion() - PHOS Cluster in %d cells fiducial volume: icol %d, irow %d, Module %d?",
 706              fNCellsFromPHOSBorder, icol, irow, relId[0]-1);
 707                         if (okcol && okrow ) printf(" YES \n");
 708                         else  printf(" NO: column ok? %d, row ok? %d \n",okcol,okrow);
 709                 }
 710         }//PHOS
 711
 712         if (okcol && okrow) return kTRUE;
 713         else                return kFALSE;
 714
 715 }
 716
 717 //__________________________________________________________________________________________________________
 718 Bool_t AliCalorimeterUtils::ClusterContainsBadChannel(TString calorimeter, UShort_t* cellList, Int_t nCells)
 719 {
 720         // Check that in the cluster cells, there is no bad channel of those stored
 721         // in fEMCALBadChannelMap or fPHOSBadChannelMap
 722
 723         if (!fRemoveBadChannels) return kFALSE;
 724         //printf("fEMCALBadChannelMap %p, fPHOSBadChannelMap %p \n",fEMCALBadChannelMap,fPHOSBadChannelMap);
 725         if(calorimeter == "EMCAL" && !fEMCALRecoUtils->GetEMCALChannelStatusMap(0)) return kFALSE;
 726         if(calorimeter == "PHOS"  && !fPHOSBadChannelMap)  return kFALSE;
 727
 728         Int_t icol = -1;
 729         Int_t irow = -1;
 730         Int_t imod = -1;
 731         for(Int_t iCell = 0; iCell<nCells; iCell++){
 732
 733                 //Get the column and row
 734                 if(calorimeter == "EMCAL"){
 735       return fEMCALRecoUtils->ClusterContainsBadChannel((AliEMCALGeometry*)fEMCALGeo,cellList,nCells);
 736                 }
 737                 else if(calorimeter=="PHOS"){
 738                         Int_t    relId[4];
 739                         fPHOSGeo->AbsToRelNumbering(cellList[iCell],relId);
 740                         irow = relId[2];
 741                         icol = relId[3];
 742                         imod = relId[0]-1;
 743                         if(fPHOSBadChannelMap->GetEntries() <= imod)continue;
 744       //printf("PHOS bad channels imod %d, icol %d, irow %d\n",imod, irow, icol);
 745                         if(GetPHOSChannelStatus(imod, irow, icol)) return kTRUE;
 746                 }
 747                 else return kFALSE;
 748
 749         }// cell cluster loop
 750
 751         return kFALSE;
 752
 753 }
 754
 755 //_______________________________________________________________
 756 void AliCalorimeterUtils::CorrectClusterEnergy(AliVCluster *clus)
 757 {
 758   // Correct cluster energy non linearity
 759
 760   clus->SetE(fEMCALRecoUtils->CorrectClusterEnergyLinearity(clus));
 761
 762 }
 763
 764 //________________________________________________________________________________________
 765 Int_t  AliCalorimeterUtils::GetMaxEnergyCell(AliVCaloCells* cells, const AliVCluster* clu,
 766                                              Float_t & clusterFraction) const
 767 {
 768
 769   //For a given CaloCluster gets the absId of the cell
 770   //with maximum energy deposit.
 771
 772   if( !clu || !cells ){
 773     AliInfo("Cluster or cells pointer is null!");
 774     return -1;
 775   }
 776
 777   Double_t eMax        =-1.;
 778   Double_t eTot        = 0.;
 779   Double_t eCell       =-1.;
 780   Float_t  fraction    = 1.;
 781   Float_t  recalFactor = 1.;
 782   Int_t    cellAbsId   =-1 , absId =-1 ;
 783   Int_t    iSupMod     =-1 , ieta  =-1 , iphi = -1, iRCU = -1;
 784
 785   TString           calo = "EMCAL";
 786   if(clu->IsPHOS()) calo = "PHOS";
 787
 788   for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
 789
 790     cellAbsId = clu->GetCellAbsId(iDig);
 791
 792     fraction  = clu->GetCellAmplitudeFraction(iDig);
 793     if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
 794
 795     iSupMod = GetModuleNumberCellIndexes(cellAbsId, calo, ieta, iphi, iRCU);
 796
 797     if(IsRecalibrationOn()) {
 798       if(calo=="EMCAL") recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
 799       else              recalFactor = GetPHOSChannelRecalibrationFactor (iSupMod,iphi,ieta);
 800     }
 801
 802     eCell  = cells->GetCellAmplitude(cellAbsId)*fraction*recalFactor;
 803
 804     if(eCell > eMax)  {
 805       eMax  = eCell;
 806       absId = cellAbsId;
 807     }
 808
 809     eTot+=eCell;
 810
 811   }// cell loop
 812
 813   if(eTot > 0.1)
 814     clusterFraction = (eTot-eMax)/eTot; //Do not use cluster energy in case it was corrected for non linearity.
 815   else
 816     clusterFraction =1.;
 817
 818   return absId;
 819
 820 }
 821
 822 //___________________________________________________________________________________
 823 AliVTrack * AliCalorimeterUtils::GetMatchedTrack(AliVCluster* cluster,
 824                                                  AliVEvent* event, Int_t index) const
 825 {
 826   // Get the matched track given its index, usually just the first match
 827   // Since it is different for ESDs and AODs here it is a wrap method to do it
 828
 829   AliVTrack *track = 0;
 830
 831   // EMCAL case only when matching is recalculated
 832   if(cluster->IsEMCAL() && IsRecalculationOfClusterTrackMatchingOn())
 833   {
 834     Int_t trackIndex = fEMCALRecoUtils->GetMatchedTrackIndex(cluster->GetID());
 835     //printf("track index %d, cluster ID %d \n ",trackIndex,cluster->GetID());
 836
 837     if(trackIndex < 0 )
 838     {
 839       printf("AliCalorimeterUtils::GetMatchedTrack() - Wrong track index %d, from recalculation\n", trackIndex);
 840     }
 841     else
 842     {
 843       track = dynamic_cast<AliVTrack*> (event->GetTrack(trackIndex));
 844     }
 845
 846     return track ;
 847
 848   }
 849
 850   // Normal case, get info from ESD or AOD
 851   // ESDs
 852   if(!strcmp("AliESDCaloCluster",Form("%s",cluster->ClassName())))
 853   {
 854     Int_t iESDtrack = cluster->GetTrackMatchedIndex();
 855
 856     if(iESDtrack < 0 )
 857     {
 858       printf("AliCalorimeterUtils::GetMatchedTrack() - Wrong track index %d\n", index);
 859       return 0x0;
 860     }
 861
 862     track = dynamic_cast<AliVTrack*> (event->GetTrack(iESDtrack));
 863
 864   }
 865   else // AODs
 866   {
 867     if(cluster->GetNTracksMatched() > 0 )
 868       track = dynamic_cast<AliVTrack*>(cluster->GetTrackMatched(index));
 869   }
 870
 871   return track ;
 872
 873 }
 874 //______________________________________________________________________________________________
 875 Float_t AliCalorimeterUtils::GetMCECellClusFracCorrection(Float_t eCell, Float_t eCluster) const
 876 {
 877   // Correction factor for cell energy in cluster to temptatively match Data and MC
 878   if( eCluster <= 0 || eCluster < eCell )
 879   {
 880     printf("AliCalorimeterUtils::GetMCECellClusFracCorrection() - Bad values eCell=%f, eCluster %f\n",eCell,eCluster);
 881     return 1;
 882   }
 883
 884   Float_t frac       = eCell / eCluster;
 885
 886   Float_t correction = fMCECellClusFracCorrParam[0] +
 887                        TMath::Exp( frac*fMCECellClusFracCorrParam[2]+fMCECellClusFracCorrParam[1] ) +
 888                        fMCECellClusFracCorrParam[3]/TMath::Sqrt(frac);
 889
 890 //  printf("AliCalorimeterUtils::GetMCECellClusFracCorrection(eCell=%f, eCluster %f, frac %f) = %f\n",eCell, eCluster, frac, correction);
 891 //  printf("\t %2.2f + TMath::Exp( %2.3f*%2.2f + %2.2f ) + %2.2f/TMath::Sqrt(%2.3f)) = %f\n",
 892 //         fMCECellClusFracCorrParam[0],frac,fMCECellClusFracCorrParam[2],fMCECellClusFracCorrParam[1],fMCECellClusFracCorrParam[3], frac, correction);
 893
 894   return correction;
 895 }
 896
 897 //_____________________________________________________________________________________________________
 898 Int_t AliCalorimeterUtils::GetModuleNumber(AliAODPWG4Particle * particle, AliVEvent * inputEvent) const
 899 {
 900         //Get the EMCAL/PHOS module number that corresponds to this particle
 901
 902         Int_t absId = -1;
 903         if(particle->GetDetector()=="EMCAL")
 904   {
 905                 fEMCALGeo->GetAbsCellIdFromEtaPhi(particle->Eta(),particle->Phi(), absId);
 906
 907                 if(fDebug > 2)
 908                   printf("AliCalorimeterUtils::GetModuleNumber(PWG4AOD) - EMCAL: cluster eta %f, phi %f, absid %d, SuperModule %d\n",
 909              particle->Eta(), particle->Phi()*TMath::RadToDeg(),absId, fEMCALGeo->GetSuperModuleNumber(absId));
 910
 911                 return fEMCALGeo->GetSuperModuleNumber(absId) ;
 912         }//EMCAL
 913         else if(particle->GetDetector()=="PHOS")
 914   {
 915     // In case we use the MC reader, the input are TParticles,
 916     // in this case use the corresponing method in PHOS Geometry to get the particle.
 917     if(strcmp(inputEvent->ClassName(), "AliMCEvent") == 0 )
 918     {
 919       Int_t mod =-1;
 920       Double_t z = 0., x=0.;
 921       TParticle* primary = 0x0;
 922       AliStack * stack = ((AliMCEvent*)inputEvent)->Stack();
 923       if(stack)
 924       {
 925         primary = stack->Particle(particle->GetCaloLabel(0));
 926       }
 927       else
 928       {
 929         Fatal("GetModuleNumber(PWG4AOD)", "Stack not available, stop!");
 930       }
 931
 932       if(primary)
 933       {
 934         fPHOSGeo->ImpactOnEmc(primary,mod,z,x) ;
 935       }
 936       else
 937       {
 938         Fatal("GetModuleNumber(PWG4AOD)", "Primary not available, stop!");
 939       }
 940       return mod;
 941     }
 942     // Input are ESDs or AODs, get the PHOS module number like this.
 943     else
 944     {
 945       //FIXME
 946       //AliVCluster *cluster = inputEvent->GetCaloCluster(particle->GetCaloLabel(0));
 947       //return GetModuleNumber(cluster);
 948       //MEFIX
 949       return -1;
 950     }
 951         }//PHOS
 952
 953         return -1;
 954 }
 955
 956 //_____________________________________________________________________
 957 Int_t AliCalorimeterUtils::GetModuleNumber(AliVCluster * cluster) const
 958 {
 959         //Get the EMCAL/PHOS module number that corresponds to this cluster
 960
 961   if(!cluster)
 962   {
 963     if(fDebug > 1) printf("AliCalorimeterUtils::GetModuleNumber() - NUL Cluster, please check!!!");
 964
 965     return -1;
 966   }
 967
 968   if ( cluster->GetNCells() <= 0 ) return -1;
 969
 970         Int_t absId = cluster->GetCellAbsId(0);
 971
 972   if ( absId < 0 ) return -1;
 973
 974         if( cluster->IsEMCAL() )
 975   {
 976                 if(fDebug > 2) printf("AliCalorimeterUtils::GetModuleNumber() - EMCAL absid %d, SuperModule %d\n",absId, fEMCALGeo->GetSuperModuleNumber(absId));
 977
 978                 return fEMCALGeo->GetSuperModuleNumber(absId) ;
 979         }//EMCAL
 980         else if ( cluster->IsPHOS() )
 981   {
 982                 Int_t    relId[4];
 983     fPHOSGeo->AbsToRelNumbering(absId,relId);
 984
 985     if(fDebug > 2) printf("AliCalorimeterUtils::GetModuleNumber() - PHOS absid %d Module %d\n",absId, relId[0]-1);
 986
 987     return relId[0]-1;
 988   }//PHOS
 989
 990         return -1;
 991 }
 992
 993 //___________________________________________________________________________________________________
 994 Int_t AliCalorimeterUtils::GetModuleNumberCellIndexes(Int_t absId, TString calo,
 995                                                       Int_t & icol, Int_t & irow, Int_t & iRCU) const
 996 {
 997         //Get the EMCAL/PHOS module, columns, row and RCU number that corresponds to this absId
 998
 999         Int_t imod = -1;
1000
1001         if ( absId < 0) return -1 ;
1002
1003   if ( calo == "EMCAL" )
1004   {
1005     Int_t iTower = -1, iIphi = -1, iIeta = -1;
1006     fEMCALGeo->GetCellIndex(absId,imod,iTower,iIphi,iIeta);
1007     fEMCALGeo->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
1008
1009     if(imod < 0 || irow < 0 || icol < 0 )
1010     {
1011       Fatal("GetModuleNumberCellIndexes()","Negative value for super module: %d, or cell icol: %d, or cell irow: %d, check EMCAL geometry name\n",imod,icol,irow);
1012     }
1013
1014     //RCU0
1015     if(imod < 10 )
1016     {
1017       if      (0<=irow&&irow<8)                       iRCU=0; // first cable row
1018       else if (8<=irow&&irow<16 &&  0<=icol&&icol<24) iRCU=0; // first half;
1019       //second cable row
1020       //RCU1
1021       else if (8<=irow&&irow<16 && 24<=icol&&icol<48) iRCU=1; // second half;
1022       //second cable row
1023       else if (16<=irow&&irow<24)                     iRCU=1; // third cable row
1024
1025       if (imod%2==1) iRCU = 1 - iRCU; // swap for odd=C side, to allow us to cable both sides the same
1026     }
1027     else
1028     {
1029       // Last 2 SM have one single SRU, just assign RCU 0
1030       iRCU = 0 ;
1031     }
1032
1033     if (iRCU<0)
1034     {
1035       Fatal("GetModuleNumberCellIndexes()","Wrong EMCAL RCU number = %d\n", iRCU);
1036     }
1037
1038     return imod ;
1039
1040   }//EMCAL
1041   else //PHOS
1042   {
1043     Int_t    relId[4];
1044     fPHOSGeo->AbsToRelNumbering(absId,relId);
1045     irow = relId[2];
1046     icol = relId[3];
1047     imod = relId[0]-1;
1048     iRCU= (Int_t)(relId[2]-1)/16 ;
1049     //Int_t iBranch= (Int_t)(relid[3]-1)/28 ; //0 to 1
1050     if (iRCU >= 4)
1051     {
1052       Fatal("GetModuleNumberCellIndexes()","Wrong PHOS RCU number = %d\n", iRCU);
1053     }
1054     return imod;
1055   }//PHOS
1056
1057         return -1;
1058 }
1059
1060 //___________________________________________________________________________________________
1061 Int_t AliCalorimeterUtils::GetNumberOfLocalMaxima(AliVCluster* cluster, AliVCaloCells* cells)
1062 {
1063   // Find local maxima in cluster
1064
1065   const Int_t   nc = cluster->GetNCells();
1066   Int_t   absIdList[nc];
1067   Float_t maxEList[nc];
1068
1069   Int_t nMax = GetNumberOfLocalMaxima(cluster, cells, absIdList, maxEList);
1070
1071   return nMax;
1072
1073 }
1074
1075 //___________________________________________________________________________________________
1076 Int_t AliCalorimeterUtils::GetNumberOfLocalMaxima(AliVCluster* cluster, AliVCaloCells* cells,
1077                                                   Int_t *absIdList,     Float_t *maxEList)
1078 {
1079   // Find local maxima in cluster
1080
1081   Int_t iDigitN = 0 ;
1082   Int_t iDigit  = 0 ;
1083   Int_t absId1 = -1 ;
1084   Int_t absId2 = -1 ;
1085   const Int_t nCells = cluster->GetNCells();
1086
1087   Float_t eCluster =  RecalibrateClusterEnergy(cluster, cells);// recalculate cluster energy, avoid non lin correction.
1088
1089   Float_t simuTotWeight = 0;
1090   if(fMCECellClusFracCorrOn)
1091   {
1092     simuTotWeight =  RecalibrateClusterEnergyWeightCell(cluster, cells,eCluster);// same but apply a weight
1093     simuTotWeight/= eCluster;
1094   }
1095
1096   TString calorimeter = "EMCAL";
1097   if(!cluster->IsEMCAL()) calorimeter = "PHOS";
1098
1099   //printf("cluster : ncells %d \n",nCells);
1100
1101   Float_t emax  = 0;
1102   Int_t   idmax =-1;
1103   for(iDigit = 0; iDigit < nCells ; iDigit++)
1104   {
1105     absIdList[iDigit] = cluster->GetCellsAbsId()[iDigit]  ;
1106     Float_t en = cells->GetCellAmplitude(absIdList[iDigit]);
1107     RecalibrateCellAmplitude(en,calorimeter,absIdList[iDigit]);
1108
1109     if(fMCECellClusFracCorrOn)
1110       en*=GetMCECellClusFracCorrection(en,eCluster)/simuTotWeight;
1111
1112     if( en > emax )
1113     {
1114       emax  = en ;
1115       idmax = absIdList[iDigit] ;
1116     }
1117     //Int_t icol = -1, irow = -1, iRCU = -1;
1118     //Int_t sm = GetModuleNumberCellIndexes(absIdList[iDigit], calorimeter, icol, irow, iRCU) ;
1119     //printf("\t cell %d, id %d, sm %d, col %d, row %d, e %f\n", iDigit, absIdList[iDigit], sm, icol, irow, en );
1120   }
1121
1122   for(iDigit = 0 ; iDigit < nCells; iDigit++)
1123   {
1124     if(absIdList[iDigit]>=0)
1125     {
1126       absId1 = cluster->GetCellsAbsId()[iDigit];
1127
1128       Float_t en1 = cells->GetCellAmplitude(absId1);
1129       RecalibrateCellAmplitude(en1,calorimeter,absId1);
1130
1131       if(fMCECellClusFracCorrOn)
1132         en1*=GetMCECellClusFracCorrection(en1,eCluster)/simuTotWeight;
1133
1134       //printf("%d : absIDi %d, E %f\n",iDigit, absId1,en1);
1135
1136       for(iDigitN = 0; iDigitN < nCells; iDigitN++)
1137       {
1138         absId2 = cluster->GetCellsAbsId()[iDigitN] ;
1139
1140         if(absId2==-1 || absId2==absId1) continue;
1141
1142         //printf("\t %d : absIDj %d\n",iDigitN, absId2);
1143
1144         Float_t en2 = cells->GetCellAmplitude(absId2);
1145         RecalibrateCellAmplitude(en2,calorimeter,absId2);
1146
1147         if(fMCECellClusFracCorrOn)
1148           en2*=GetMCECellClusFracCorrection(en2,eCluster)/simuTotWeight;
1149
1150         //printf("\t %d : absIDj %d, E %f\n",iDigitN, absId2,en2);
1151
1152         if ( AreNeighbours(calorimeter, absId1, absId2) )
1153         {
1154           // printf("\t \t Neighbours \n");
1155           if (en1 > en2 )
1156           {
1157             absIdList[iDigitN] = -1 ;
1158             //printf("\t \t indexN %d not local max\n",iDigitN);
1159             // but may be digit too is not local max ?
1160             if(en1 < en2 + fLocMaxCutEDiff) {
1161               //printf("\t \t index %d not local max cause locMaxCutEDiff\n",iDigit);
1162               absIdList[iDigit] = -1 ;
1163             }
1164           }
1165           else
1166           {
1167             absIdList[iDigit] = -1 ;
1168             //printf("\t \t index %d not local max\n",iDigitN);
1169             // but may be digitN too is not local max ?
1170             if(en1 > en2 - fLocMaxCutEDiff)
1171             {
1172               absIdList[iDigitN] = -1 ;
1173               //printf("\t \t indexN %d not local max cause locMaxCutEDiff\n",iDigit);
1174             }
1175           }
1176         } // if Are neighbours
1177         //else printf("\t \t NOT Neighbours \n");
1178       } // while digitN
1179     } // slot not empty
1180   } // while digit
1181
1182   iDigitN = 0 ;
1183   for(iDigit = 0; iDigit < nCells; iDigit++)
1184   {
1185     if(absIdList[iDigit]>=0 )
1186     {
1187       absIdList[iDigitN] = absIdList[iDigit] ;
1188
1189       Float_t en = cells->GetCellAmplitude(absIdList[iDigit]);
1190       RecalibrateCellAmplitude(en,calorimeter,absIdList[iDigit]);
1191
1192       if(fMCECellClusFracCorrOn)
1193         en*=GetMCECellClusFracCorrection(en,eCluster)/simuTotWeight;
1194
1195       if(en < fLocMaxCutE) continue; // Maxima only with seed energy at least
1196
1197       maxEList[iDigitN] = en ;
1198
1199       //printf("Local max %d, id %d, en %f\n", iDigit,absIdList[iDigitN],en);
1200       iDigitN++ ;
1201     }
1202   }
1203
1204   if(iDigitN == 0)
1205   {
1206     if(fDebug > 0)
1207       printf("AliCalorimeterUtils::GetNumberOfLocalMaxima() - No local maxima found, assign highest energy cell as maxima, id %d, en cell %2.2f, en cluster %2.2f\n",
1208              idmax,emax,cluster->E());
1209     iDigitN      = 1     ;
1210     maxEList[0]  = emax  ;
1211     absIdList[0] = idmax ;
1212   }
1213
1214   if(fDebug > 0)
1215   {
1216     printf("AliCalorimeterUtils::GetNumberOfLocalMaxima() - In cluster E %2.2f (wth non lin. %2.2f), M02 %2.2f, M20 %2.2f, N maxima %d \n",
1217            cluster->E(),eCluster, cluster->GetM02(),cluster->GetM20(), iDigitN);
1218
1219     if(fDebug > 1) for(Int_t imax = 0; imax < iDigitN; imax++)
1220     {
1221       printf(" \t i %d, absId %d, Ecell %f\n",imax,absIdList[imax],maxEList[imax]);
1222     }
1223   }
1224
1225   return iDigitN ;
1226
1227 }
1228
1229 //____________________________________
1230 TString AliCalorimeterUtils::GetPass()
1231 {
1232   // Get passx from filename.
1233
1234   if (!AliAnalysisManager::GetAnalysisManager()->GetTree())
1235   {
1236     AliError("AliCalorimeterUtils::GetPass() - Pointer to tree = 0, returning null\n");
1237     return TString("");
1238   }
1239
1240   if (!AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile())
1241   {
1242     AliError("AliCalorimeterUtils::GetPass() - Null pointer input file, returning null\n");
1243     return TString("");
1244   }
1245
1246   TString pass(AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile()->GetName());
1247   if      (pass.Contains("ass1")) return TString("pass1");
1248   else if (pass.Contains("ass2")) return TString("pass2");
1249   else if (pass.Contains("ass3")) return TString("pass3");
1250   else if (pass.Contains("ass4")) return TString("pass4");
1251   else if (pass.Contains("ass5")) return TString("pass5");
1252   else if (pass.Contains("LHC11c") && pass.Contains("spc_calo") ) return TString("spc_calo");
1253   else if (pass.Contains("calo") || pass.Contains("high_lumi"))
1254   {
1255     printf("AliCalorimeterUtils::GetPass() - Path contains <calo> or <high-lumi>, set as <pass1>\n");
1256     return TString("pass1");
1257   }
1258
1259   // No condition fullfilled, give a default value
1260   printf("AliCalorimeterUtils::GetPass() - Pass number string not found \n");
1261   return TString("");
1262
1263 }
1264
1265 //________________________________________
1266 void AliCalorimeterUtils::InitParameters()
1267 {
1268   //Initialize the parameters of the analysis.
1269
1270   fEMCALGeoName         = "";
1271   fPHOSGeoName          = "";
1272
1273   fEMCALGeoMatrixSet    = kFALSE;
1274   fPHOSGeoMatrixSet     = kFALSE;
1275
1276   fRemoveBadChannels    = kFALSE;
1277
1278   fNCellsFromPHOSBorder = 0;
1279
1280   fLocMaxCutE           = 0.1 ;
1281   fLocMaxCutEDiff       = 0.0 ;
1282
1283   //  fMaskCellColumns = new Int_t[fNMaskCellColumns];
1284   //  fMaskCellColumns[0] = 6 ;  fMaskCellColumns[1] = 7 ;  fMaskCellColumns[2] = 8 ;
1285   //  fMaskCellColumns[3] = 35;  fMaskCellColumns[4] = 36;  fMaskCellColumns[5] = 37;
1286   //  fMaskCellColumns[6] = 12+AliEMCALGeoParams::fgkEMCALCols; fMaskCellColumns[7] = 13+AliEMCALGeoParams::fgkEMCALCols;
1287   //  fMaskCellColumns[8] = 40+AliEMCALGeoParams::fgkEMCALCols; fMaskCellColumns[9] = 41+AliEMCALGeoParams::fgkEMCALCols;
1288   //  fMaskCellColumns[10]= 42+AliEMCALGeoParams::fgkEMCALCols;
1289
1290   fOADBSet      = kFALSE;
1291   fOADBForEMCAL = kTRUE ;
1292   fOADBForPHOS  = kFALSE;
1293
1294   fOADBFilePathEMCAL = "$ALICE_ROOT/OADB/EMCAL" ;
1295   fOADBFilePathPHOS  = "$ALICE_ROOT/OADB/PHOS"  ;
1296
1297   fImportGeometryFromFile = kTRUE;
1298   fImportGeometryFilePath = "";
1299
1300   fNSuperModulesUsed = 22;
1301
1302   fMCECellClusFracCorrParam[0] = 0.78;
1303   fMCECellClusFracCorrParam[1] =-1.8;
1304   fMCECellClusFracCorrParam[2] =-6.3;
1305   fMCECellClusFracCorrParam[3] = 0.014;
1306
1307 }
1308
1309
1310 //_____________________________________________________
1311 void AliCalorimeterUtils::InitPHOSBadChannelStatusMap()
1312 {
1313   //Init PHOS bad channels map
1314   if(fDebug > 0 )printf("AliCalorimeterUtils::InitPHOSBadChannelStatusMap()\n");
1315   //In order to avoid rewriting the same histograms
1316   Bool_t oldStatus = TH1::AddDirectoryStatus();
1317   TH1::AddDirectory(kFALSE);
1318
1319   fPHOSBadChannelMap = new TObjArray(5);
1320   for (int i = 0; i < 5; i++)fPHOSBadChannelMap->Add(new TH2I(Form("PHOS_BadMap_mod%d",i),Form("PHOS_BadMap_mod%d",i), 64, 0, 64, 56, 0, 56));
1321
1322   fPHOSBadChannelMap->SetOwner(kTRUE);
1323   fPHOSBadChannelMap->Compress();
1324
1325   //In order to avoid rewriting the same histograms
1326   TH1::AddDirectory(oldStatus);
1327 }
1328
1329 //______________________________________________________
1330 void AliCalorimeterUtils::InitPHOSRecalibrationFactors()
1331 {
1332         //Init EMCAL recalibration factors
1333         if(fDebug > 0 )printf("AliCalorimeterUtils::InitPHOSRecalibrationFactors()\n");
1334         //In order to avoid rewriting the same histograms
1335         Bool_t oldStatus = TH1::AddDirectoryStatus();
1336         TH1::AddDirectory(kFALSE);
1337
1338         fPHOSRecalibrationFactors = new TObjArray(5);
1339         for (int i = 0; i < 5; i++)fPHOSRecalibrationFactors->Add(new TH2F(Form("PHOSRecalFactors_Mod%d",i),Form("PHOSRecalFactors_Mod%d",i), 64, 0, 64, 56, 0, 56));
1340         //Init the histograms with 1
1341         for (Int_t m = 0; m < 5; m++) {
1342                 for (Int_t i = 0; i < 56; i++) {
1343                         for (Int_t j = 0; j < 64; j++) {
1344                                 SetPHOSChannelRecalibrationFactor(m,j,i,1.);
1345                         }
1346                 }
1347         }
1348         fPHOSRecalibrationFactors->SetOwner(kTRUE);
1349         fPHOSRecalibrationFactors->Compress();
1350
1351         //In order to avoid rewriting the same histograms
1352         TH1::AddDirectory(oldStatus);
1353 }
1354
1355
1356 //__________________________________________________________
1357 void AliCalorimeterUtils::InitEMCALGeometry(Int_t runnumber)
1358 {
1359         //Initialize EMCAL geometry if it did not exist previously
1360
1361         if (!fEMCALGeo)
1362   {
1363     if(fEMCALGeoName=="")
1364     {
1365       if     (runnumber <  140000 &&
1366               runnumber >= 100000)   fEMCALGeoName = "EMCAL_FIRSTYEARV1";
1367       else if(runnumber >= 140000 &&
1368               runnumber <  171000)   fEMCALGeoName = "EMCAL_COMPLETEV1";
1369       else                           fEMCALGeoName = "EMCAL_COMPLETE12SMV1";
1370       printf("AliCalorimeterUtils::InitEMCALGeometry() - Set EMCAL geometry name to <%s> for run %d\n",fEMCALGeoName.Data(),runnumber);
1371     }
1372
1373                 fEMCALGeo = AliEMCALGeometry::GetInstance(fEMCALGeoName);
1374
1375     // Init geometry, I do not like much to do it like this ...
1376     if(fImportGeometryFromFile && !gGeoManager)
1377     {
1378       if(fImportGeometryFilePath=="") // If not specified, set location depending on run number
1379       {
1380         // "$ALICE_ROOT/EVE/alice-data/default_geo.root"
1381         if     (runnumber <  140000 &&
1382                 runnumber >= 100000) fImportGeometryFilePath = "$ALICE_ROOT/OADB/EMCAL/geometry_2010.root";
1383         if     (runnumber >= 140000 &&
1384                 runnumber <  171000) fImportGeometryFilePath = "$ALICE_ROOT/OADB/EMCAL/geometry_2011.root";
1385         else                         fImportGeometryFilePath = "$ALICE_ROOT/OADB/EMCAL/geometry_2012.root"; // 2012-2013
1386
1387       }
1388       printf("AliCalorimeterUtils::InitEMCALGeometry() - Import %s\n",fImportGeometryFilePath.Data());
1389       TGeoManager::Import(fImportGeometryFilePath) ; // default need file "geometry.root" in local dir!!!!
1390     }
1391
1392
1393                 if(fDebug > 0)
1394     {
1395                         printf("AliCalorimeterUtils::InitEMCALGeometry(run=%d)",runnumber);
1396                         if (!gGeoManager) printf(" - Careful!, gGeoManager not loaded, load misalign matrices");
1397                         printf("\n");
1398                 }
1399         }
1400 }
1401
1402 //_________________________________________________________
1403 void AliCalorimeterUtils::InitPHOSGeometry(Int_t runnumber)
1404 {
1405         //Initialize PHOS geometry if it did not exist previously
1406
1407         if (!fPHOSGeo)
1408   {
1409     if(fPHOSGeoName=="") fPHOSGeoName = "PHOSgeo";
1410
1411                 fPHOSGeo = new AliPHOSGeoUtils(fPHOSGeoName);
1412
1413                 if(fDebug > 0)
1414     {
1415                         printf("AliCalorimeterUtils::InitPHOSGeometry(run=%d)",runnumber);
1416                         if (!gGeoManager) printf(" - Careful!, gGeoManager not loaded, load misalign matrices");
1417                         printf("\n");
1418                 }
1419         }
1420 }
1421
1422 //_______________________________________________________________________________________________
1423 Bool_t AliCalorimeterUtils::IsMCParticleInCalorimeterAcceptance(TString calo, TParticle* particle)
1424 {
1425   // Check that a MC ESD is in the calorimeter acceptance
1426
1427   if( (!IsPHOSGeoMatrixSet() && calo == "PHOS" ) ||
1428      (!IsPHOSGeoMatrixSet() && calo == "EMCAL")   )
1429   {
1430     AliFatal(Form("Careful Geo Matrix for %s is not set, use AliFidutialCut instead \n",calo.Data()));
1431     return kFALSE ;
1432   }
1433
1434   if(calo == "PHOS" )
1435   {
1436     Int_t mod = 0 ;
1437     Double_t x = 0, z = 0 ;
1438
1439     return GetPHOSGeometry()->ImpactOnEmc( particle, mod, z, x);
1440   }
1441   else if(calo == "EMCAL")
1442   {
1443     Int_t absID = 0 ;
1444
1445     GetEMCALGeometry()->GetAbsCellIdFromEtaPhi( particle->Eta(), particle->Phi(), absID);
1446
1447     if( absID >= 0) return kTRUE  ;
1448     else            return kFALSE ;
1449   }
1450
1451   return kFALSE ;
1452 }
1453
1454 //______________________________________________________________________________________________________
1455 Bool_t AliCalorimeterUtils::IsMCParticleInCalorimeterAcceptance(TString calo, AliAODMCParticle* particle)
1456 {
1457   // Check that a MC AOD is in the calorimeter acceptance
1458
1459   if( (!IsPHOSGeoMatrixSet() && calo == "PHOS" ) ||
1460       (!IsPHOSGeoMatrixSet() && calo == "EMCAL")   )
1461   {
1462     AliFatal(Form("Careful Geo Matrix for %s is not set, use AliFidutialCut instead \n",calo.Data()));
1463     return kFALSE ;
1464   }
1465
1466   if(calo == "PHOS" )
1467   {
1468     Int_t mod = 0 ;
1469     Double_t x = 0, z = 0 ;
1470
1471     Double_t vtx[]={ particle->Xv(), particle->Yv(), particle->Zv() } ;
1472
1473     return GetPHOSGeometry()->ImpactOnEmc(vtx, particle->Theta(), particle->Phi(), mod, z, x) ;
1474   }
1475   else if(calo == "EMCAL")
1476   {
1477     Int_t absID = 0 ;
1478
1479     GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(particle->Eta(),particle->Phi(),absID);
1480
1481     if( absID >= 0) return kTRUE  ;
1482     else            return kFALSE ;
1483   }
1484
1485   return kFALSE ;
1486 }
1487
1488 //_______________________________________________________________________
1489 Bool_t AliCalorimeterUtils::MaskFrameCluster(Int_t iSM, Int_t ieta) const
1490 {
1491   //Check if cell is in one of the regions where we have significant amount
1492   //of material in front. Only EMCAL
1493
1494   Int_t icol = ieta;
1495   if(iSM%2) icol+=48; // Impair SM, shift index [0-47] to [48-96]
1496
1497   if (fNMaskCellColumns && fMaskCellColumns)
1498   {
1499     for (Int_t imask = 0; imask < fNMaskCellColumns; imask++)
1500     {
1501       if(icol==fMaskCellColumns[imask]) return kTRUE;
1502     }
1503   }
1504
1505   return kFALSE;
1506
1507 }
1508
1509 //_________________________________________________________
1510 void AliCalorimeterUtils::Print(const Option_t * opt) const
1511 {
1512
1513   //Print some relevant parameters set for the analysis
1514   if(! opt)
1515     return;
1516
1517   printf("***** Print: %s %s ******\n", GetName(), GetTitle() ) ;
1518   printf("Remove Clusters with bad channels? %d\n",fRemoveBadChannels);
1519   printf("Remove Clusters with max cell at less than %d cells from EMCAL border and %d cells from PHOS border\n",
1520          fEMCALRecoUtils->GetNumberOfCellsFromEMCALBorder(), fNCellsFromPHOSBorder);
1521   if(fEMCALRecoUtils->IsEMCALNoBorderAtEta0()) printf("Do not remove EMCAL clusters at Eta = 0\n");
1522   printf("Recalibrate Clusters? %d, run by run  %d\n",fRecalibration,fRunDependentCorrection);
1523   printf("Recalculate Clusters Position? %d\n",fRecalculatePosition);
1524   printf("Recalculate Clusters Energy? %d\n",fCorrectELinearity);
1525   printf("Matching criteria: dR < %2.2f[cm], dZ < %2.2f[cm]\n",fCutR,fCutZ);
1526
1527   printf("Loc. Max. E > %2.2f\n",       fLocMaxCutE);
1528   printf("Loc. Max. E Diff > %2.2f\n",  fLocMaxCutEDiff);
1529
1530   printf("    \n") ;
1531 }
1532
1533 //_____________________________________________________________________________________________
1534 void AliCalorimeterUtils::RecalibrateCellAmplitude(Float_t & amp, TString calo, Int_t id) const
1535 {
1536   //Recaculate cell energy if recalibration factor
1537
1538   Int_t icol     = -1; Int_t irow     = -1; Int_t iRCU     = -1;
1539   Int_t nModule  = GetModuleNumberCellIndexes(id,calo, icol, irow, iRCU);
1540
1541   if (IsRecalibrationOn())
1542   {
1543     if(calo == "PHOS")
1544     {
1545       amp *= GetPHOSChannelRecalibrationFactor(nModule,icol,irow);
1546     }
1547     else
1548     {
1549       amp *= GetEMCALChannelRecalibrationFactor(nModule,icol,irow);
1550     }
1551   }
1552 }
1553
1554 //____________________________________________________________________________________________________
1555 void AliCalorimeterUtils::RecalibrateCellTime(Double_t & time, TString calo, Int_t id, Int_t bc) const
1556 {
1557   // Recalculate time if time recalibration available for EMCAL
1558   // not ready for PHOS
1559
1560   if(calo == "EMCAL" && GetEMCALRecoUtils()->IsTimeRecalibrationOn())
1561   {
1562     GetEMCALRecoUtils()->RecalibrateCellTime(id,bc,time);
1563   }
1564
1565 }
1566
1567 //__________________________________________________________________________
1568 Float_t AliCalorimeterUtils::RecalibrateClusterEnergy(AliVCluster * cluster,
1569                                                       AliVCaloCells * cells)
1570 {
1571         // Recalibrate the cluster energy, considering the recalibration map and the energy of the cells that compose the cluster.
1572
1573   //Initialize some used variables
1574         Float_t frac  = 0., energy = 0.;
1575
1576         if(cells)
1577   {
1578     //Get the cluster number of cells and list of absId, check what kind of cluster do we have.
1579
1580     UShort_t * index    = cluster->GetCellsAbsId() ;
1581     Double_t * fraction = cluster->GetCellsAmplitudeFraction() ;
1582
1583     Int_t ncells     = cluster->GetNCells();
1584
1585     TString calo     = "EMCAL";
1586     if(cluster->IsPHOS())
1587       calo = "PHOS";
1588
1589     //Loop on the cells, get the cell amplitude and recalibration factor, multiply and and to the new energy
1590     for(Int_t icell = 0; icell < ncells; icell++){
1591
1592       Int_t absId = index[icell];
1593
1594       frac =  fraction[icell];
1595       if(frac < 1e-3) frac = 1; //in case of EMCAL, this is set as 0, not used.
1596
1597       Float_t amp = cells->GetCellAmplitude(absId);
1598       RecalibrateCellAmplitude(amp,calo, absId);
1599
1600       if(fDebug>2)
1601         printf("AliCalorimeterUtils::RecalibrateClusterEnergy() - recalibrate cell: %s, cell fraction %f, cell energy %f\n",
1602                calo.Data(),frac,cells->GetCellAmplitude(absId));
1603
1604       energy += amp*frac;
1605     }
1606
1607     if(fDebug>1)
1608       printf("AliCalorimeterUtils::RecalibrateClusterEnergy() - Energy before %f, after %f\n",cluster->E(),energy);
1609
1610         }// cells available
1611   else
1612   {
1613     Fatal("RecalibrateClusterEnergy()","Cells pointer does not exist!");
1614   }
1615
1616         return energy;
1617 }
1618
1619 //__________________________________________________________________________
1620 Float_t AliCalorimeterUtils::RecalibrateClusterEnergyWeightCell(AliVCluster * cluster,
1621                                                                 AliVCaloCells * cells, Float_t energyOrg)
1622 {
1623         // Recalibrate the cluster energy, considering the recalibration map and the energy of the cells that compose the cluster.
1624   // Also consider reweighting of cells energy
1625
1626   //Initialize some used variables
1627         Float_t frac  = 0., energy = 0.;
1628
1629         if(cells)
1630   {
1631     //Get the cluster number of cells and list of absId, check what kind of cluster do we have.
1632
1633     UShort_t * index    = cluster->GetCellsAbsId() ;
1634     Double_t * fraction = cluster->GetCellsAmplitudeFraction() ;
1635
1636     Int_t ncells     = cluster->GetNCells();
1637
1638     TString calo     = "EMCAL";
1639     if(cluster->IsPHOS())
1640       calo = "PHOS";
1641
1642     //Loop on the cells, get the cell amplitude and recalibration factor, multiply and and to the new energy
1643     for(Int_t icell = 0; icell < ncells; icell++){
1644
1645       Int_t absId = index[icell];
1646
1647       frac =  fraction[icell];
1648       if(frac < 1e-3) frac = 1; //in case of EMCAL, this is set as 0, not used.
1649
1650       Float_t amp = cells->GetCellAmplitude(absId);
1651       RecalibrateCellAmplitude(amp,calo, absId);
1652
1653       amp*=GetMCECellClusFracCorrection(amp,energyOrg);
1654
1655       if(fDebug>2)
1656         printf("AliCalorimeterUtils::RecalibrateClusterEnergy() - recalibrate cell: %s, cell fraction %f, cell energy %f\n",
1657                calo.Data(),frac,cells->GetCellAmplitude(absId));
1658
1659       energy += amp*frac;
1660     }
1661
1662     if(fDebug>1)
1663       printf("AliCalorimeterUtils::RecalibrateClusterEnergy() - Energy before %f, after %f\n",cluster->E(),energy);
1664
1665         }// cells available
1666   else
1667   {
1668     Fatal("RecalibrateClusterEnergy()","Cells pointer does not exist!");
1669   }
1670
1671         return energy;
1672 }
1673
1674
1675 //__________________________________________________________________________________________
1676 void AliCalorimeterUtils::RecalculateClusterPosition(AliVCaloCells* cells, AliVCluster* clu)
1677 {
1678
1679   //Recalculate EMCAL cluster position
1680
1681   fEMCALRecoUtils->RecalculateClusterPosition((AliEMCALGeometry*)fEMCALGeo, cells,clu);
1682
1683 }
1684
1685 //________________________________________________________________________________
1686 void AliCalorimeterUtils::RecalculateClusterTrackMatching(AliVEvent * event,
1687                                                           TObjArray* clusterArray)
1688 {
1689   //Recalculate track matching
1690
1691   if (fRecalculateMatching)
1692   {
1693     fEMCALRecoUtils->FindMatches(event,clusterArray,fEMCALGeo)   ;
1694     //AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (event);
1695     //if(esdevent){
1696     //  fEMCALRecoUtils->SetClusterMatchedToTrack(esdevent)  ;
1697     //  fEMCALRecoUtils->SetTracksMatchedToCluster(esdevent) ;
1698     //}
1699   }
1700 }
1701
1702 //___________________________________________________________________________
1703 void AliCalorimeterUtils::SplitEnergy(Int_t absId1, Int_t absId2,
1704                                       AliVCluster* cluster,
1705                                       AliVCaloCells* cells,
1706                                       //Float_t & e1, Float_t & e2,
1707                                       AliAODCaloCluster* cluster1,
1708                                       AliAODCaloCluster* cluster2,
1709                                       Int_t nMax, Int_t eventNumber)
1710 {
1711
1712   // Split energy of cluster between the 2 local maxima, sum energy on 3x3, and if the 2
1713   // maxima are too close and have common cells, split the energy between the 2
1714
1715   TH2F* hClusterMap    = 0 ;
1716   TH2F* hClusterLocMax = 0 ;
1717   TH2F* hCluster1      = 0 ;
1718   TH2F* hCluster2      = 0 ;
1719
1720   if(fPlotCluster)
1721   {
1722     hClusterMap    = new TH2F("hClusterMap","Cluster Map",48,0,48,24,0,24);
1723     hClusterLocMax = new TH2F("hClusterLocMax","Cluster 2 highest local maxima",48,0,48,24,0,24);
1724     hCluster1      = new TH2F("hCluster1","Cluster 1",48,0,48,24,0,24);
1725     hCluster2      = new TH2F("hCluster2","Cluster 2",48,0,48,24,0,24);
1726     hClusterMap    ->SetXTitle("column");
1727     hClusterMap    ->SetYTitle("row");
1728     hClusterLocMax ->SetXTitle("column");
1729     hClusterLocMax ->SetYTitle("row");
1730     hCluster1      ->SetXTitle("column");
1731     hCluster1      ->SetYTitle("row");
1732     hCluster2      ->SetXTitle("column");
1733     hCluster2      ->SetYTitle("row");
1734   }
1735
1736   TString calorimeter = "EMCAL";
1737   if(cluster->IsPHOS())
1738   {
1739     calorimeter="PHOS";
1740     printf("AliCalorimeterUtils::SplitEnerg() Not supported for PHOS yet \n");
1741     return;
1742   }
1743
1744   const Int_t ncells  = cluster->GetNCells();
1745   Int_t absIdList[ncells];
1746
1747   Float_t e1 = 0,  e2   = 0 ;
1748   Int_t icol = -1, irow = -1, iRCU = -1, sm = -1;
1749   Float_t eCluster = 0;
1750   Float_t minCol = 100, minRow = 100, maxCol = -1, maxRow = -1;
1751   for(Int_t iDigit  = 0; iDigit < ncells; iDigit++ )
1752   {
1753     absIdList[iDigit] = cluster->GetCellsAbsId()[iDigit];
1754
1755     Float_t ec = cells->GetCellAmplitude(absIdList[iDigit]);
1756     RecalibrateCellAmplitude(ec,calorimeter, absIdList[iDigit]);
1757     eCluster+=ec;
1758
1759     if(fPlotCluster)
1760     {
1761       //printf("iDigit %d, absId %d, Ecell %f\n",iDigit,absIdList[iDigit], cells->GetCellAmplitude(absIdList[iDigit]));
1762       sm = GetModuleNumberCellIndexes(absIdList[iDigit], calorimeter, icol, irow, iRCU) ;
1763       if(sm > -1 && sm  < 12) // just to avoid compilation warning
1764       {
1765         if(icol > maxCol) maxCol = icol;
1766         if(icol < minCol) minCol = icol;
1767         if(irow > maxRow) maxRow = irow;
1768         if(irow < minRow) minRow = irow;
1769         hClusterMap->Fill(icol,irow,ec);
1770       }
1771     }
1772
1773   }
1774
1775   // Init counters and variables
1776   Int_t ncells1 = 1 ;
1777   UShort_t absIdList1[9] ;
1778   Double_t fracList1 [9] ;
1779   absIdList1[0] = absId1 ;
1780   fracList1 [0] = 1. ;
1781
1782   Float_t ecell1 = cells->GetCellAmplitude(absId1);
1783   RecalibrateCellAmplitude(ecell1, calorimeter, absId1);
1784   e1 =  ecell1;
1785
1786   Int_t ncells2 = 1 ;
1787   UShort_t absIdList2[9] ;
1788   Double_t fracList2 [9] ;
1789   absIdList2[0] = absId2 ;
1790   fracList2 [0] = 1. ;
1791
1792   Float_t ecell2 = cells->GetCellAmplitude(absId2);
1793   RecalibrateCellAmplitude(ecell2, calorimeter, absId2);
1794   e2 =  ecell2;
1795
1796   if(fPlotCluster)
1797   {
1798     Int_t icol1 = -1, irow1 = -1, icol2 = -1, irow2 = -1;
1799     sm = GetModuleNumberCellIndexes(absId1, calorimeter, icol1, irow1, iRCU) ;
1800     hClusterLocMax->Fill(icol1,irow1,ecell1);
1801     sm = GetModuleNumberCellIndexes(absId2, calorimeter, icol2, irow2, iRCU) ;
1802     hClusterLocMax->Fill(icol2,irow2,ecell2);
1803   }
1804
1805   // Very rough way to share the cluster energy
1806   Float_t eRemain = (eCluster-ecell1-ecell2)/2;
1807   Float_t shareFraction1 = ecell1/eCluster+eRemain/eCluster;
1808   Float_t shareFraction2 = ecell2/eCluster+eRemain/eCluster;
1809
1810   for(Int_t iDigit = 0; iDigit < ncells; iDigit++)
1811   {
1812     Int_t absId = absIdList[iDigit];
1813
1814     if(absId==absId1 || absId==absId2 || absId < 0) continue;
1815
1816     Float_t ecell = cells->GetCellAmplitude(absId);
1817     RecalibrateCellAmplitude(ecell, calorimeter, absId);
1818
1819     if(AreNeighbours(calorimeter, absId1,absId ))
1820     {
1821       absIdList1[ncells1]= absId;
1822
1823       if(AreNeighbours(calorimeter, absId2,absId ))
1824       {
1825         fracList1[ncells1] = shareFraction1;
1826         e1 += ecell*shareFraction1;
1827       }
1828       else
1829       {
1830         fracList1[ncells1] = 1.;
1831         e1 += ecell;
1832       }
1833
1834       ncells1++;
1835
1836     } // neigbour to cell1
1837
1838     if(AreNeighbours(calorimeter, absId2,absId ))
1839     {
1840       absIdList2[ncells2]= absId;
1841
1842       if(AreNeighbours(calorimeter, absId1,absId ))
1843       {
1844         fracList2[ncells2] = shareFraction2;
1845         e2 += ecell*shareFraction2;
1846       }
1847       else
1848       {
1849         fracList2[ncells2] = 1.;
1850         e2 += ecell;
1851       }
1852
1853       ncells2++;
1854
1855     } // neigbour to cell2
1856
1857   }
1858
1859   if(GetDebug() > 1) printf("AliCalorimeterUtils::SplitEnergy() - n Local Max %d, Cluster energy  = %f, Ecell1 = %f, Ecell2 = %f, Enew1 = %f, Enew2 = %f, Remain %f, \n ncells %d, ncells1 %d, ncells2 %d, f1 %f, f2  %f, sum f12 = %f \n",
1860                             nMax, eCluster,ecell1,ecell2,e1,e2,eCluster-e1-e2,ncells,ncells1,ncells2,shareFraction1,shareFraction2,shareFraction1+shareFraction2);
1861
1862   cluster1->SetE(e1);
1863   cluster2->SetE(e2);
1864
1865   cluster1->SetNCells(ncells1);
1866   cluster2->SetNCells(ncells2);
1867
1868   cluster1->SetCellsAbsId(absIdList1);
1869   cluster2->SetCellsAbsId(absIdList2);
1870
1871   cluster1->SetCellsAmplitudeFraction(fracList1);
1872   cluster2->SetCellsAmplitudeFraction(fracList2);
1873
1874   //Correct linearity
1875   CorrectClusterEnergy(cluster1) ;
1876   CorrectClusterEnergy(cluster2) ;
1877
1878   if(calorimeter=="EMCAL")
1879   {
1880     GetEMCALRecoUtils()->RecalculateClusterPosition(GetEMCALGeometry(), cells, cluster1);
1881     GetEMCALRecoUtils()->RecalculateClusterPosition(GetEMCALGeometry(), cells, cluster2);
1882   }
1883
1884   if(fPlotCluster)
1885   {
1886     //printf("Cells of cluster1: ");
1887     for(Int_t iDigit  = 0; iDigit < ncells1; iDigit++ )
1888     {
1889       //printf(" %d ",absIdList1[iDigit]);
1890
1891       sm = GetModuleNumberCellIndexes(absIdList1[iDigit], calorimeter, icol, irow, iRCU) ;
1892
1893       Float_t ecell = cells->GetCellAmplitude(absIdList1[iDigit]);
1894       RecalibrateCellAmplitude(ecell, calorimeter, absIdList1[iDigit]);
1895
1896       if( AreNeighbours(calorimeter, absId2,absIdList1[iDigit]) && absId1!=absIdList1[iDigit])
1897         hCluster1->Fill(icol,irow,ecell*shareFraction1);
1898       else
1899         hCluster1->Fill(icol,irow,ecell);
1900     }
1901
1902     //printf(" \n ");
1903     //printf("Cells of cluster2: ");
1904
1905     for(Int_t iDigit  = 0; iDigit < ncells2; iDigit++ )
1906     {
1907       //printf(" %d ",absIdList2[iDigit]);
1908
1909       sm = GetModuleNumberCellIndexes(absIdList2[iDigit], calorimeter, icol, irow, iRCU) ;
1910
1911       Float_t ecell = cells->GetCellAmplitude(absIdList2[iDigit]);
1912       RecalibrateCellAmplitude(ecell, calorimeter, absIdList2[iDigit]);
1913
1914       if( AreNeighbours(calorimeter, absId1,absIdList2[iDigit])  && absId2!=absIdList2[iDigit])
1915         hCluster2->Fill(icol,irow,ecell*shareFraction2);
1916       else
1917         hCluster2->Fill(icol,irow,ecell);
1918
1919     }
1920     //printf(" \n ");
1921
1922     gStyle->SetPadRightMargin(0.1);
1923     gStyle->SetPadLeftMargin(0.1);
1924     gStyle->SetOptStat(0);
1925     gStyle->SetOptFit(000000);
1926
1927     if(maxCol-minCol > maxRow-minRow)
1928     {
1929       maxRow+= maxCol-minCol;
1930     }
1931     else
1932     {
1933       maxCol+= maxRow-minRow;
1934     }
1935
1936     TCanvas  * c= new TCanvas("canvas", "canvas", 4000, 4000) ;
1937     c->Divide(2,2);
1938     c->cd(1);
1939     gPad->SetGridy();
1940     gPad->SetGridx();
1941     gPad->SetLogz();
1942     hClusterMap    ->SetAxisRange(minCol, maxCol,"X");
1943     hClusterMap    ->SetAxisRange(minRow, maxRow,"Y");
1944     hClusterMap    ->Draw("colz TEXT");
1945     c->cd(2);
1946     gPad->SetGridy();
1947     gPad->SetGridx();
1948     gPad->SetLogz();
1949     hClusterLocMax ->SetAxisRange(minCol, maxCol,"X");
1950     hClusterLocMax ->SetAxisRange(minRow, maxRow,"Y");
1951     hClusterLocMax ->Draw("colz TEXT");
1952     c->cd(3);
1953     gPad->SetGridy();
1954     gPad->SetGridx();
1955     gPad->SetLogz();
1956     hCluster1      ->SetAxisRange(minCol, maxCol,"X");
1957     hCluster1      ->SetAxisRange(minRow, maxRow,"Y");
1958     hCluster1      ->Draw("colz TEXT");
1959     c->cd(4);
1960     gPad->SetGridy();
1961     gPad->SetGridx();
1962     gPad->SetLogz();
1963     hCluster2      ->SetAxisRange(minCol, maxCol,"X");
1964     hCluster2      ->SetAxisRange(minRow, maxRow,"Y");
1965     hCluster2      ->Draw("colz TEXT");
1966
1967     if(eCluster > 6 )c->Print(Form("clusterFigures/Event%d_E%1.0f_nMax%d_NCell1_%d_NCell2_%d.eps",
1968                                    eventNumber,cluster->E(),nMax,ncells1,ncells2));
1969
1970     delete c;
1971     delete hClusterMap;
1972     delete hClusterLocMax;
1973     delete hCluster1;
1974     delete hCluster2;
1975   }
1976 }
1977