PWG/CaloTrackCorrBase/AliCalorimeterUtils.cxx

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