PWG/CaloTrackCorrBase/AliCalorimeterUtils.cxx

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