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