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