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