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