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