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