PWGGA/CaloTrackCorrelations/AliAnaCalorimeterQA.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 to check results from simulations or reconstructed real data.
  18 // Fill few histograms and do some checking plots
  19 //
  20 //-- Author: Gustavo Conesa (INFN-LNF)
  21 //_________________________________________________________________________
  22
  23
  24 // --- ROOT system ---
  25 #include <TObjArray.h>
  26 #include <TParticle.h>
  27 #include <TDatabasePDG.h>
  28 #include <TH3F.h>
  29 #include <TObjString.h>
  30
  31 //---- AliRoot system ----
  32 #include "AliAnaCalorimeterQA.h"
  33 #include "AliCaloTrackReader.h"
  34 #include "AliStack.h"
  35 #include "AliVCaloCells.h"
  36 #include "AliFiducialCut.h"
  37 #include "AliVCluster.h"
  38 #include "AliVTrack.h"
  39 #include "AliVEvent.h"
  40 #include "AliVEventHandler.h"
  41 #include "AliAODMCParticle.h"
  42 #include "AliMCAnalysisUtils.h"
  43
  44 // --- Detectors ---
  45 #include "AliPHOSGeoUtils.h"
  46 #include "AliEMCALGeometry.h"
  47
  48 ClassImp(AliAnaCalorimeterQA)
  49
  50 //________________________________________
  51 AliAnaCalorimeterQA::AliAnaCalorimeterQA() :
  52 AliAnaCaloTrackCorrBaseClass(),        fCalorimeter(""),
  53
  54 //Switches
  55 fFillAllCellTimeHisto(kTRUE),
  56 fFillAllPosHisto(kFALSE),              fFillAllPosHisto2(kTRUE),
  57 fFillAllTH12(kFALSE),                  fFillAllTH3(kTRUE),
  58 fFillAllTMHisto(kTRUE),                fFillAllPi0Histo(kTRUE),
  59 fCorrelate(kTRUE),                     fStudyBadClusters(kFALSE),
  60 fStudyClustersAsymmetry(kFALSE),       fStudyExotic(kFALSE),
  61 fStudyWeight(kFALSE),
  62
  63 //Parameters and cuts
  64 fNModules(12),                         fNRCU(2),
  65 fNMaxCols(48),                         fNMaxRows(24),
  66 fTimeCutMin(-10000),                   fTimeCutMax(10000),
  67 fEMCALCellAmpMin(0),                   fPHOSCellAmpMin(0),
  68
  69 // Exotic
  70 fExoNECrossCuts(0),                    fExoECrossCuts(),
  71 fExoNDTimeCuts(0),                     fExoDTimeCuts(),
  72
  73 //Histograms
  74 fhE(0),                                fhPt(0),
  75 fhPhi(0),                              fhEta(0),                               fhEtaPhiE(0),
  76 fhECharged(0),                         fhPtCharged(0),
  77 fhPhiCharged(0),                       fhEtaCharged(0),                        fhEtaPhiECharged(0),
  78
  79 //Invariant mass
  80 fhIM(0 ),                              fhAsym(0),
  81
  82 fhNCellsPerCluster(0),                 fhNCellsPerClusterNoCut(0),             fhNClusters(0),
  83
  84 //Timing
  85 fhClusterTimeEnergy(0),                fhCellTimeSpreadRespectToCellMax(0),
  86 fhCellIdCellLargeTimeSpread(0),        fhClusterPairDiffTimeE(0),
  87 fhClusterMaxCellCloseCellRatio(0),     fhClusterMaxCellCloseCellDiff(0),
  88 fhClusterMaxCellDiff(0),               fhClusterMaxCellDiffNoCut(0),
  89 fhClusterMaxCellDiffAverageTime(0),    fhClusterMaxCellDiffWeightedTime(0),
  90 fhClusterMaxCellECross(0),
  91 fhLambda0(0),                          fhLambda1(0),                           fhDispersion(0),
  92
  93 //bad clusters
  94 fhBadClusterEnergy(0),                 fhBadClusterTimeEnergy(0),
  95 fhBadClusterPairDiffTimeE(0),          fhBadCellTimeSpreadRespectToCellMax(0),
  96 fhBadClusterMaxCellCloseCellRatio(0),  fhBadClusterMaxCellCloseCellDiff(0),    fhBadClusterMaxCellDiff(0),
  97 fhBadClusterMaxCellDiffAverageTime(0), fhBadClusterMaxCellDiffWeightedTime(0),
  98 fhBadClusterMaxCellECross(0),
  99 fhBadClusterDeltaIEtaDeltaIPhiE0(0),   fhBadClusterDeltaIEtaDeltaIPhiE2(0),
 100 fhBadClusterDeltaIEtaDeltaIPhiE6(0),   fhBadClusterDeltaIA(0),
 101
 102 //Position
 103 fhRNCells(0),                          fhXNCells(0),
 104 fhYNCells(0),                          fhZNCells(0),
 105 fhRE(0),                               fhXE(0),
 106 fhYE(0),                               fhZE(0),
 107 fhXYZ(0),
 108 fhRCellE(0),                           fhXCellE(0),
 109 fhYCellE(0),                           fhZCellE(0),
 110 fhXYZCell(0),
 111 fhDeltaCellClusterRNCells(0),          fhDeltaCellClusterXNCells(0),
 112 fhDeltaCellClusterYNCells(0),          fhDeltaCellClusterZNCells(0),
 113 fhDeltaCellClusterRE(0),               fhDeltaCellClusterXE(0),
 114 fhDeltaCellClusterYE(0),               fhDeltaCellClusterZE(0),
 115
 116 // Cells
 117 fhNCells(0),                           fhAmplitude(0),
 118 fhAmpId(0),                            fhEtaPhiAmp(0),
 119 fhTime(0),                             fhTimeVz(0),
 120 fhTimeId(0),                           fhTimeAmp(0),
 121 fhCellECross(0),
 122 fhCaloCorrNClusters(0),                fhCaloCorrEClusters(0),
 123 fhCaloCorrNCells(0),                   fhCaloCorrECells(0),
 124 fhCaloV0SCorrNClusters(0),             fhCaloV0SCorrEClusters(0),
 125 fhCaloV0SCorrNCells(0),                fhCaloV0SCorrECells(0),
 126 fhCaloV0MCorrNClusters(0),             fhCaloV0MCorrEClusters(0),
 127 fhCaloV0MCorrNCells(0),                fhCaloV0MCorrECells(0),
 128 fhCaloTrackMCorrNClusters(0),          fhCaloTrackMCorrEClusters(0),
 129 fhCaloTrackMCorrNCells(0),             fhCaloTrackMCorrECells(0),
 130 fhCaloCenNClusters(0),                 fhCaloCenEClusters(0),
 131 fhCaloCenNCells(0),                    fhCaloCenECells(0),
 132 fhCaloEvPNClusters(0),                 fhCaloEvPEClusters(0),
 133 fhCaloEvPNCells(0),                    fhCaloEvPECells(0),
 134 //Super-Module dependent histgrams
 135 fhEMod(0),                             fhAmpMod(0),                            fhTimeMod(0),
 136 fhNClustersMod(0),                     fhNCellsMod(0),
 137 fhNCellsPerClusterMod(0),              fhNCellsPerClusterModNoCut(0),
 138
 139 fhGridCells(0),                        fhGridCellsE(0),                        fhGridCellsTime(0),
 140 fhTimeAmpPerRCU(0),                    fhIMMod(0),
 141
 142 // Weight studies
 143 fhECellClusterRatio(0),                fhECellClusterLogRatio(0),
 144 fhEMaxCellClusterRatio(0),             fhEMaxCellClusterLogRatio(0),
 145
 146 fhExoL0ECross(0),                      fhExoL1ECross(0),
 147
 148 // MC and reco
 149 fhRecoMCE(),                           fhRecoMCPhi(),                          fhRecoMCEta(),
 150 fhRecoMCDeltaE(),                      fhRecoMCRatioE(),
 151 fhRecoMCDeltaPhi(),                    fhRecoMCDeltaEta(),
 152
 153 // MC only
 154 fhGenMCE(),                            fhGenMCEtaPhi(),
 155 fhGenMCAccE(),                         fhGenMCAccEtaPhi(),
 156
 157 //matched MC
 158 fhEMVxyz(0),                           fhEMR(0),
 159 fhHaVxyz(0),                           fhHaR(0),
 160 fh1EOverP(0),                          fh2dR(0),
 161 fh2EledEdx(0),                         fh2MatchdEdx(0),
 162 fhMCEle1EOverP(0),                     fhMCEle1dR(0),                          fhMCEle2MatchdEdx(0),
 163 fhMCChHad1EOverP(0),                   fhMCChHad1dR(0),                        fhMCChHad2MatchdEdx(0),
 164 fhMCNeutral1EOverP(0),                 fhMCNeutral1dR(0),                      fhMCNeutral2MatchdEdx(0), fh1EOverPR02(0),
 165 fhMCEle1EOverPR02(0),                  fhMCChHad1EOverPR02(0),                 fhMCNeutral1EOverPR02(0),
 166 fh1EleEOverP(0),                       fhMCEle1EleEOverP(0),
 167 fhMCChHad1EleEOverP(0),                fhMCNeutral1EleEOverP(0),
 168 fhTrackMatchedDEta(0),                 fhTrackMatchedDPhi(0),                  fhTrackMatchedDEtaDPhi(0),
 169 fhTrackMatchedDEtaPos(0),              fhTrackMatchedDPhiPos(0),               fhTrackMatchedDEtaDPhiPos(0)
 170 {
 171   //Default Ctor
 172
 173   //Weight studies
 174   for(Int_t i =0; i < 14; i++){
 175     fhLambda0ForW0[i] = 0;
 176     //fhLambda1ForW0[i] = 0;
 177
 178     for(Int_t j = 0; j < 5; j++){
 179       fhLambda0ForW0MC[i][j] = 0;
 180       //fhLambda1ForW0MC[i][j] = 0;
 181     }
 182
 183   }
 184
 185   //Cluster size
 186   fhDeltaIEtaDeltaIPhiE0[0] = 0 ;         fhDeltaIEtaDeltaIPhiE2[0] = 0;          fhDeltaIEtaDeltaIPhiE6[0] = 0;
 187   fhDeltaIEtaDeltaIPhiE0[1] = 0 ;         fhDeltaIEtaDeltaIPhiE2[1] = 0;          fhDeltaIEtaDeltaIPhiE6[1] = 0;
 188   fhDeltaIA[0]              = 0 ;         fhDeltaIAL0[0]            = 0;          fhDeltaIAL1[0]            = 0;
 189   fhDeltaIA[1]              = 0 ;         fhDeltaIAL0[1]            = 0;          fhDeltaIAL1[1]            = 0;
 190   fhDeltaIANCells[0]        = 0 ;         fhDeltaIANCells[1]        = 0;
 191   fhDeltaIAMC[0]            = 0 ;         fhDeltaIAMC[1]            = 0;
 192   fhDeltaIAMC[2]            = 0 ;         fhDeltaIAMC[3]            = 0;
 193
 194   // Exotic
 195   for (Int_t ie = 0; ie < 10 ; ie++)
 196   {
 197     fhExoDTime[ie] = 0;
 198     for (Int_t idt = 0; idt < 5 ; idt++)
 199     {
 200       fhExoNCell    [ie][idt] = 0;
 201       fhExoL0       [ie][idt] = 0;
 202       fhExoL1       [ie][idt] = 0;
 203       fhExoECross   [ie][idt] = 0;
 204       fhExoTime     [ie][idt] = 0;
 205       fhExoL0NCell  [ie][idt] = 0;
 206       fhExoL1NCell  [ie][idt] = 0;
 207     }
 208   }
 209
 210   // MC
 211
 212   for(Int_t i = 0; i < 6; i++){
 213
 214     fhRecoMCE[i][0]         = 0; fhRecoMCE[i][1]        = 0;
 215     fhRecoMCPhi[i][0]       = 0; fhRecoMCPhi[i][1]      = 0;
 216     fhRecoMCEta[i][0]       = 0; fhRecoMCEta[i][1]      = 0;
 217     fhRecoMCDeltaE[i][0]    = 0; fhRecoMCDeltaE[i][1]   = 0;
 218     fhRecoMCRatioE[i][0]    = 0; fhRecoMCRatioE[i][1]   = 0;
 219     fhRecoMCDeltaPhi[i][0]  = 0; fhRecoMCDeltaPhi[i][1] = 0;
 220     fhRecoMCDeltaEta[i][0]  = 0; fhRecoMCDeltaEta[i][1] = 0;
 221
 222   }
 223
 224   //Initialize parameters
 225   InitParameters();
 226 }
 227
 228 //______________________________________________________________________________________________________________________
 229 void AliAnaCalorimeterQA::BadClusterHistograms(AliVCluster* clus, const TObjArray *caloClusters, AliVCaloCells * cells,
 230                                                Int_t absIdMax, Double_t maxCellFraction, Float_t eCrossFrac,
 231                                                Double_t tmax,  Double_t timeAverages[2] )
 232 {
 233   //Bad cluster histograms
 234
 235   //  printf("AliAnaCalorimeterQA::BadClusterHistograms() - Event %d - Calorimeter %s \n \t  E %f, n cells %d, max cell absId %d, maxCellFrac %f\n",
 236   //         GetReader()->GetEventNumber(), fCalorimeter.Data(),
 237   //         clus->E(),clus->GetNCells(),absIdMax,maxCellFraction);
 238
 239   fhBadClusterEnergy     ->Fill(clus->E());
 240   Double_t tof = clus->GetTOF()*1.e9;
 241   fhBadClusterTimeEnergy   ->Fill(clus->E(),tof);
 242   fhBadClusterMaxCellDiff  ->Fill(clus->E(),maxCellFraction);
 243   fhBadClusterMaxCellECross->Fill(clus->E(),eCrossFrac);
 244
 245   if(fStudyClustersAsymmetry) ClusterAsymmetryHistograms(clus,absIdMax,kFALSE);
 246
 247   //Clusters in event time differencem bad minus good
 248
 249   for(Int_t iclus2 = 0; iclus2 < caloClusters->GetEntriesFast(); iclus2++ ){
 250
 251     AliVCluster* clus2 =  (AliVCluster*)caloClusters->At(iclus2);
 252
 253     if(clus->GetID()==clus2->GetID()) continue;
 254
 255     Float_t maxCellFraction2 = 0.;
 256     Int_t absIdMax2 = GetCaloUtils()->GetMaxEnergyCell(cells, clus2,maxCellFraction2);
 257     if(IsGoodCluster(absIdMax2,cells)){
 258       Double_t tof2   = clus2->GetTOF()*1.e9;
 259       fhBadClusterPairDiffTimeE  ->Fill(clus->E(), (tof-tof2));
 260     }
 261
 262   } // loop
 263
 264   // Max cell compared to other cells in cluster
 265   if(fFillAllCellTimeHisto)
 266   {
 267     fhBadClusterMaxCellDiffAverageTime      ->Fill(clus->E(),tmax-timeAverages[0]);
 268     fhBadClusterMaxCellDiffWeightedTime     ->Fill(clus->E(),tmax-timeAverages[1]);
 269   }
 270
 271   for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++)
 272   {
 273     Int_t absId  = clus->GetCellsAbsId()[ipos];
 274     if(absId!=absIdMax && cells->GetCellAmplitude(absIdMax) > 0.01)
 275     {
 276       Float_t frac = cells->GetCellAmplitude(absId)/cells->GetCellAmplitude(absIdMax);
 277
 278       fhBadClusterMaxCellCloseCellRatio->Fill(clus->E(),frac);
 279       fhBadClusterMaxCellCloseCellDiff ->Fill(clus->E(),cells->GetCellAmplitude(absIdMax)-cells->GetCellAmplitude(absId));
 280
 281       if(fFillAllCellTimeHisto)
 282       {
 283         Double_t time  = cells->GetCellTime(absId);
 284         GetCaloUtils()->RecalibrateCellTime(time, fCalorimeter, absId,GetReader()->GetInputEvent()->GetBunchCrossNumber());
 285
 286         Float_t diff = (tmax-time*1e9);
 287         fhBadCellTimeSpreadRespectToCellMax->Fill(clus->E(), diff);
 288
 289       }
 290     }// Not max
 291   }//loop
 292
 293 }
 294
 295 //______________________________________________________________________
 296 void AliAnaCalorimeterQA::CalculateAverageTime(AliVCluster *clus,
 297                                                AliVCaloCells* cells,
 298                                                Double_t timeAverages[2])
 299 {
 300   // Calculate time averages and weights
 301
 302   // First recalculate energy in case non linearity was applied
 303   Float_t  energy = 0;
 304   Float_t  ampMax = 0, amp = 0;
 305 //  Int_t    absIdMax =-1;
 306   for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++)
 307   {
 308     Int_t id       = clus->GetCellsAbsId()[ipos];
 309
 310     //Recalibrate cell energy if needed
 311     amp = cells->GetCellAmplitude(id);
 312     GetCaloUtils()->RecalibrateCellAmplitude(amp,fCalorimeter, id);
 313
 314     energy    += amp;
 315
 316     if(amp> ampMax)
 317     {
 318       ampMax   = amp;
 319 //      absIdMax = id;
 320     }
 321
 322   } // energy loop
 323
 324   // Calculate average time of cells in cluster and weighted average
 325   Double_t aTime  = 0;
 326   Double_t wTime  = 0;
 327   Float_t  wTot   = 0;
 328   Double_t time   = 0;
 329   Int_t    id     =-1;
 330   Double_t w      = 0;
 331   Int_t    ncells = clus->GetNCells();
 332   for (Int_t ipos = 0; ipos < ncells; ipos++)
 333   {
 334     id   = clus ->GetCellsAbsId()[ipos];
 335     amp  = cells->GetCellAmplitude(id);
 336     time = cells->GetCellTime(id);
 337
 338     //Recalibrate energy and time
 339     GetCaloUtils()->RecalibrateCellAmplitude(amp , fCalorimeter, id);
 340     GetCaloUtils()->RecalibrateCellTime     (time, fCalorimeter, id, GetReader()->GetInputEvent()->GetBunchCrossNumber());
 341
 342     w      = GetCaloUtils()->GetEMCALRecoUtils()->GetCellWeight(cells->GetCellAmplitude(id),energy);
 343     aTime += time*1e9;
 344     wTime += time*1e9 * w;
 345     wTot  += w;
 346
 347   }
 348
 349   if(ncells > 0) aTime /= ncells;
 350   else           aTime  = 0;
 351
 352   if(wTot   > 0) wTime /= wTot;
 353   else           wTime  = 0;
 354
 355   timeAverages[0] = aTime;
 356   timeAverages[1] = wTime;
 357
 358 }
 359
 360 //____________________________________________________________
 361 void AliAnaCalorimeterQA::CellHistograms(AliVCaloCells *cells)
 362 {
 363   // Plot histograms related to cells only
 364
 365   Int_t ncells = cells->GetNumberOfCells();
 366
 367   if(GetDebug() > 0)
 368     printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - %s cell entries %d\n", fCalorimeter.Data(), ncells );
 369
 370   //Init arrays and used variables
 371   Int_t *nCellsInModule = new Int_t[fNModules];
 372   for(Int_t imod = 0; imod < fNModules; imod++ ) nCellsInModule[imod] = 0;
 373
 374   Int_t    icol   = -1;
 375   Int_t    irow   = -1;
 376   Int_t    iRCU   = -1;
 377   Float_t  amp    = 0.;
 378   Double_t time   = 0.;
 379   Int_t    id     = -1;
 380   Float_t  recalF = 1.;
 381   Int_t    bc     = GetReader()->GetInputEvent()->GetBunchCrossNumber();
 382
 383   for (Int_t iCell = 0; iCell < cells->GetNumberOfCells(); iCell++) {
 384     if(GetDebug() > 2)
 385       printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - Cell : amp %f, absId %d \n", cells->GetAmplitude(iCell), cells->GetCellNumber(iCell));
 386     Int_t nModule = GetModuleNumberCellIndexes(cells->GetCellNumber(iCell),fCalorimeter, icol, irow, iRCU);
 387     if(GetDebug() > 2)
 388       printf("\t module %d, column %d, row %d \n", nModule,icol,irow);
 389
 390     if(nModule < fNModules)
 391     {
 392       //Check if the cell is a bad channel
 393       if(GetCaloUtils()->IsBadChannelsRemovalSwitchedOn()){
 394         if(fCalorimeter=="EMCAL")
 395         {
 396           if(GetCaloUtils()->GetEMCALChannelStatus(nModule,icol,irow)) continue;
 397         }
 398         else
 399         {
 400           if(GetCaloUtils()->GetPHOSChannelStatus(nModule,icol,irow) ) continue;
 401         }
 402       } // use bad channel map
 403
 404       amp     = cells->GetAmplitude(iCell)*recalF;
 405       time    = cells->GetTime(iCell);
 406       id      = cells->GetCellNumber(iCell);
 407
 408       // Amplitude recalibration if set
 409       GetCaloUtils()->RecalibrateCellAmplitude(amp,  fCalorimeter, id);
 410
 411       // Time recalibration if set
 412       GetCaloUtils()->RecalibrateCellTime     (time, fCalorimeter, id, GetReader()->GetInputEvent()->GetBunchCrossNumber());
 413
 414       //Transform time to ns
 415       time *= 1.0e9;
 416
 417       if(time < fTimeCutMin || time > fTimeCutMax)
 418       {
 419           if(GetDebug() > 0 )
 420             printf("AliAnaCalorimeterQA - Remove cell with Time %f\n",time);
 421           continue;
 422       }
 423
 424       // Remove exotic cells, defined only for EMCAL
 425       if(fCalorimeter=="EMCAL" &&
 426          GetCaloUtils()->GetEMCALRecoUtils()->IsExoticCell(id, cells, bc)) continue;
 427
 428
 429       fhAmplitude->Fill(amp);
 430       fhAmpId    ->Fill(amp,id);
 431       fhAmpMod   ->Fill(amp,nModule);
 432
 433       if ((fCalorimeter=="EMCAL" && amp > fEMCALCellAmpMin) ||
 434           (fCalorimeter=="PHOS"  && amp > fPHOSCellAmpMin )   )
 435       {
 436
 437         //E cross for exotic cells
 438         if(amp > 0.01) fhCellECross->Fill(amp,1-GetECross(id,cells)/amp);
 439
 440         nCellsInModule[nModule]++ ;
 441
 442         Int_t icols = icol;
 443         Int_t irows = irow;
 444
 445         if(fCalorimeter=="EMCAL")
 446         {
 447           icols = (nModule % 2) ? icol + fNMaxCols : icol;
 448           if(nModule < 10 )
 449             irows = irow + fNMaxRows       * Int_t(nModule / 2);
 450           else // 1/3 SM
 451             irows = irow + (fNMaxRows / 3) * Int_t(nModule / 2);
 452         }
 453         else
 454         {
 455           irows = irow + fNMaxRows * nModule;
 456         }
 457
 458         fhGridCells ->Fill(icols,irows);
 459         fhGridCellsE->Fill(icols,irows,amp);
 460
 461         if(fFillAllCellTimeHisto)
 462         {
 463           //printf("%s: time %g\n",fCalorimeter.Data(), time);
 464
 465           Double_t v[3] = {0,0,0}; //vertex ;
 466           GetReader()->GetVertex(v);
 467           if(amp > 0.5) fhTimeVz   ->Fill(TMath::Abs(v[2]),time);
 468
 469           fhTime     ->Fill(time);
 470           fhTimeId   ->Fill(time,id);
 471           fhTimeAmp  ->Fill(amp,time);
 472           fhGridCellsTime->Fill(icols,irows,time);
 473           fhTimeMod  ->Fill(time,nModule);
 474           fhTimeAmpPerRCU  [nModule*fNRCU+iRCU]->Fill(amp, time);
 475
 476         }
 477       }
 478
 479       //Get Eta-Phi position of Cell
 480       if(fFillAllPosHisto)
 481       {
 482         if(fCalorimeter=="EMCAL" && GetCaloUtils()->IsEMCALGeoMatrixSet()){
 483           Float_t celleta = 0.;
 484           Float_t cellphi = 0.;
 485           GetEMCALGeometry()->EtaPhiFromIndex(id, celleta, cellphi);
 486
 487           fhEtaPhiAmp->Fill(celleta,cellphi,amp);
 488           Double_t cellpos[] = {0, 0, 0};
 489           GetEMCALGeometry()->GetGlobal(id, cellpos);
 490           fhXCellE->Fill(cellpos[0],amp)  ;
 491           fhYCellE->Fill(cellpos[1],amp)  ;
 492           fhZCellE->Fill(cellpos[2],amp)  ;
 493           Float_t rcell = TMath::Sqrt(cellpos[0]*cellpos[0]+cellpos[1]*cellpos[1]);//+cellpos[2]*cellpos[2]);
 494           fhRCellE->Fill(rcell,amp)  ;
 495           fhXYZCell->Fill(cellpos[0],cellpos[1],cellpos[2])  ;
 496         }//EMCAL Cells
 497         else if(fCalorimeter=="PHOS" && GetCaloUtils()->IsPHOSGeoMatrixSet()){
 498           TVector3 xyz;
 499           Int_t relId[4], module;
 500           Float_t xCell, zCell;
 501
 502           GetPHOSGeometry()->AbsToRelNumbering(id,relId);
 503           module = relId[0];
 504           GetPHOSGeometry()->RelPosInModule(relId,xCell,zCell);
 505           GetPHOSGeometry()->Local2Global(module,xCell,zCell,xyz);
 506           Float_t rcell = TMath::Sqrt(xyz.X()*xyz.X()+xyz.Y()*xyz.Y());
 507           fhXCellE ->Fill(xyz.X(),amp)  ;
 508           fhYCellE ->Fill(xyz.Y(),amp)  ;
 509           fhZCellE ->Fill(xyz.Z(),amp)  ;
 510           fhRCellE ->Fill(rcell  ,amp)  ;
 511           fhXYZCell->Fill(xyz.X(),xyz.Y(),xyz.Z())  ;
 512         }//PHOS cells
 513       }//fill cell position histograms
 514
 515       if     (fCalorimeter=="EMCAL" && amp > fEMCALCellAmpMin) ncells ++ ;
 516       else if(fCalorimeter=="PHOS"  && amp > fPHOSCellAmpMin)  ncells ++ ;
 517       //else
 518       //  printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - no %s CELLS passed the analysis cut\n",fCalorimeter.Data());
 519     }//nmodules
 520   }//cell loop
 521
 522   if(ncells > 0 )fhNCells->Fill(ncells) ; //fill the cells after the cut
 523
 524   //Number of cells per module
 525   for(Int_t imod = 0; imod < fNModules; imod++ ) {
 526
 527     if(GetDebug() > 1)
 528       printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - module %d calo %s cells %d\n", imod, fCalorimeter.Data(), nCellsInModule[imod]);
 529
 530     fhNCellsMod->Fill(nCellsInModule[imod],imod) ;
 531
 532   }
 533
 534   delete [] nCellsInModule;
 535
 536 }
 537
 538 //__________________________________________________________________________
 539 void AliAnaCalorimeterQA::CellInClusterPositionHistograms(AliVCluster* clus)
 540 {
 541   // Fill histograms releated to cell position
 542
 543
 544   Int_t nCaloCellsPerCluster = clus->GetNCells();
 545   UShort_t * indexList = clus->GetCellsAbsId();
 546   Float_t pos[3];
 547   clus->GetPosition(pos);
 548   Float_t clEnergy = clus->E();
 549
 550   //Loop on cluster cells
 551   for (Int_t ipos = 0; ipos < nCaloCellsPerCluster; ipos++) {
 552
 553     //  printf("Index %d\n",ipos);
 554     Int_t absId  = indexList[ipos];
 555
 556     //Get position of cell compare to cluster
 557
 558     if(fCalorimeter=="EMCAL" && GetCaloUtils()->IsEMCALGeoMatrixSet()){
 559
 560       Double_t cellpos[] = {0, 0, 0};
 561       GetEMCALGeometry()->GetGlobal(absId, cellpos);
 562
 563       fhDeltaCellClusterXNCells->Fill(pos[0]-cellpos[0],nCaloCellsPerCluster) ;
 564       fhDeltaCellClusterYNCells->Fill(pos[1]-cellpos[1],nCaloCellsPerCluster) ;
 565       fhDeltaCellClusterZNCells->Fill(pos[2]-cellpos[2],nCaloCellsPerCluster) ;
 566
 567       fhDeltaCellClusterXE->Fill(pos[0]-cellpos[0],clEnergy)  ;
 568       fhDeltaCellClusterYE->Fill(pos[1]-cellpos[1],clEnergy)  ;
 569       fhDeltaCellClusterZE->Fill(pos[2]-cellpos[2],clEnergy)  ;
 570
 571       Float_t r     = TMath::Sqrt(pos[0]    *pos[0]     + pos[1]    * pos[1]    );
 572       Float_t rcell = TMath::Sqrt(cellpos[0]*cellpos[0] + cellpos[1]* cellpos[1]);
 573
 574       fhDeltaCellClusterRNCells->Fill(r-rcell, nCaloCellsPerCluster) ;
 575       fhDeltaCellClusterRE     ->Fill(r-rcell, clEnergy)  ;
 576
 577     }//EMCAL and its matrices are available
 578     else if(fCalorimeter=="PHOS" && GetCaloUtils()->IsPHOSGeoMatrixSet()){
 579       TVector3 xyz;
 580       Int_t relId[4], module;
 581       Float_t xCell, zCell;
 582
 583       GetPHOSGeometry()->AbsToRelNumbering(absId,relId);
 584       module = relId[0];
 585       GetPHOSGeometry()->RelPosInModule(relId,xCell,zCell);
 586       GetPHOSGeometry()->Local2Global(module,xCell,zCell,xyz);
 587
 588       fhDeltaCellClusterXNCells->Fill(pos[0]-xyz.X(),nCaloCellsPerCluster) ;
 589       fhDeltaCellClusterYNCells->Fill(pos[1]-xyz.Y(),nCaloCellsPerCluster) ;
 590       fhDeltaCellClusterZNCells->Fill(pos[2]-xyz.Z(),nCaloCellsPerCluster) ;
 591
 592       fhDeltaCellClusterXE->Fill(pos[0]-xyz.X(),clEnergy)  ;
 593       fhDeltaCellClusterYE->Fill(pos[1]-xyz.Y(),clEnergy)  ;
 594       fhDeltaCellClusterZE->Fill(pos[2]-xyz.Z(),clEnergy)  ;
 595
 596       Float_t r     = TMath::Sqrt(pos[0]  * pos[0]  + pos[1]  * pos[1] );
 597       Float_t rcell = TMath::Sqrt(xyz.X() * xyz.X() + xyz.Y() * xyz.Y());
 598
 599       fhDeltaCellClusterRNCells->Fill(r-rcell, nCaloCellsPerCluster) ;
 600       fhDeltaCellClusterRE     ->Fill(r-rcell, clEnergy)  ;
 601
 602     }//PHOS and its matrices are available
 603   }// cluster cell loop
 604 }
 605
 606 //___________________________________________________________________________________________
 607 void AliAnaCalorimeterQA::ClusterAsymmetryHistograms(AliVCluster* clus, Int_t absIdMax,
 608                                                      Bool_t goodCluster)
 609 {
 610   // Study the shape of the cluster in cell units terms
 611
 612   //No use to study clusters with less than 4 cells
 613   if(clus->GetNCells() <=3 ) return;
 614
 615   Int_t dIeta = 0;
 616   Int_t dIphi = 0;
 617
 618   Int_t ietaMax=-1; Int_t iphiMax = 0; Int_t rcuMax = 0;
 619   Int_t smMax = GetModuleNumberCellIndexes(absIdMax,fCalorimeter, ietaMax, iphiMax, rcuMax);
 620
 621   for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++) {
 622
 623     Int_t absId = clus->GetCellsAbsId()[ipos];
 624
 625     Int_t ieta=-1; Int_t iphi = 0; Int_t rcu = 0;
 626     Int_t sm = GetModuleNumberCellIndexes(absId,fCalorimeter, ieta, iphi, rcu);
 627
 628     if(dIphi < TMath::Abs(iphi-iphiMax)) dIphi = TMath::Abs(iphi-iphiMax);
 629
 630     if(smMax==sm){
 631       if(dIeta < TMath::Abs(ieta-ietaMax)) dIeta = TMath::Abs(ieta-ietaMax);
 632     }
 633     else {
 634       Int_t ietaShift    = ieta;
 635       Int_t ietaMaxShift = ietaMax;
 636       if (ieta > ietaMax)  ietaMaxShift+=48;
 637       else                 ietaShift   +=48;
 638       if(dIeta < TMath::Abs(ietaShift-ietaMaxShift)) dIeta = TMath::Abs(ietaShift-ietaMaxShift);
 639     }
 640
 641   }// fill cell-cluster histogram loop
 642
 643
 644   Float_t dIA = 1.*(dIphi-dIeta)/(dIeta+dIphi);
 645
 646   if(goodCluster)
 647   {
 648
 649     // Was cluster matched?
 650     Bool_t matched = GetCaloPID()->IsTrackMatched(clus,GetCaloUtils(),GetReader()->GetInputEvent());
 651
 652     if     (clus->E() < 2 ) fhDeltaIEtaDeltaIPhiE0[matched]->Fill(dIeta,dIphi);
 653     else if(clus->E() < 6 ) fhDeltaIEtaDeltaIPhiE2[matched]->Fill(dIeta,dIphi);
 654     else                    fhDeltaIEtaDeltaIPhiE6[matched]->Fill(dIeta,dIphi);
 655
 656     fhDeltaIA[matched]->Fill(clus->E(),dIA);
 657
 658     if(clus->E() > 0.5){
 659
 660       fhDeltaIAL0[matched]    ->Fill(clus->GetM02(),dIA);
 661       fhDeltaIAL1[matched]    ->Fill(clus->GetM20(),dIA);
 662       fhDeltaIANCells[matched]->Fill(clus->GetNCells(),dIA);
 663
 664     }
 665
 666     // Origin of  clusters
 667     Int_t  nLabel = clus->GetNLabels();
 668     Int_t* labels = clus->GetLabels();
 669     if(IsDataMC()){
 670       Int_t tag = GetMCAnalysisUtils()->CheckOrigin(labels,nLabel, GetReader());
 671       if(   GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPhoton) &&
 672          !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPi0)    &&
 673          !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCEta)    &&
 674          !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion)        ){
 675         fhDeltaIAMC[0]->Fill(clus->E(),dIA);//Pure Photon
 676       }
 677       else if ( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCElectron) &&
 678                !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion)  ){
 679         fhDeltaIAMC[1]->Fill(clus->E(),dIA);//Pure electron
 680       }
 681       else if ( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPhoton) &&
 682                GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion)  ){
 683         fhDeltaIAMC[2]->Fill(clus->E(),dIA);//Converted cluster
 684       }
 685       else if(!GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPhoton)){
 686         fhDeltaIAMC[3]->Fill(clus->E(),dIA);//Hadrons
 687       }
 688
 689     }  // MC
 690   } // good cluster
 691   else
 692   {
 693     if     (clus->E() < 2 ) fhBadClusterDeltaIEtaDeltaIPhiE0->Fill(dIeta,dIphi);
 694     else if(clus->E() < 6 ) fhBadClusterDeltaIEtaDeltaIPhiE2->Fill(dIeta,dIphi);
 695     else                    fhBadClusterDeltaIEtaDeltaIPhiE6->Fill(dIeta,dIphi);
 696
 697     fhBadClusterDeltaIA->Fill(clus->E(),dIA);
 698
 699   }
 700 }
 701
 702 //__________________________________________________________________________________________________________________
 703 void AliAnaCalorimeterQA::ClusterHistograms(AliVCluster* clus, const TObjArray *caloClusters, AliVCaloCells * cells,
 704                                             Int_t absIdMax, Double_t maxCellFraction, Float_t eCrossFrac,
 705                                             Double_t tmax,  Double_t timeAverages[2])
 706 {
 707   //Fill CaloCluster related histograms
 708
 709   Double_t tof = clus->GetTOF()*1.e9;
 710
 711   fhLambda0             ->Fill(clus->E(),clus->GetM02());
 712   fhLambda1             ->Fill(clus->E(),clus->GetM20());
 713   fhDispersion          ->Fill(clus->E(),clus->GetDispersion());
 714
 715   fhClusterMaxCellDiff  ->Fill(clus->E(),maxCellFraction);
 716   fhClusterMaxCellECross->Fill(clus->E(),eCrossFrac);
 717   fhClusterTimeEnergy   ->Fill(clus->E(),tof);
 718
 719   if(fStudyClustersAsymmetry) ClusterAsymmetryHistograms(clus,absIdMax,kTRUE);
 720
 721   //Clusters in event time difference
 722   for(Int_t iclus2 = 0; iclus2 < caloClusters->GetEntriesFast(); iclus2++ ){
 723
 724     AliVCluster* clus2 =  (AliVCluster*) caloClusters->At(iclus2);
 725
 726     if(clus->GetID()==clus2->GetID()) continue;
 727
 728     if(clus->GetM02() > 0.01 && clus2->GetM02() > 0.01)
 729     {
 730       Double_t tof2   = clus2->GetTOF()*1.e9;
 731       fhClusterPairDiffTimeE  ->Fill(clus->E(), tof-tof2);
 732     }
 733   }
 734
 735   Int_t    nModule = GetModuleNumber(clus);
 736   Int_t    nCaloCellsPerCluster = clus->GetNCells();
 737
 738   if(nCaloCellsPerCluster > 1){
 739
 740     // check time of cells respect to max energy cell
 741
 742     if(fFillAllCellTimeHisto)
 743     {
 744       fhClusterMaxCellDiffAverageTime      ->Fill(clus->E(),tmax-timeAverages[0]);
 745       fhClusterMaxCellDiffWeightedTime     ->Fill(clus->E(),tmax-timeAverages[1]);
 746     }
 747
 748     for (Int_t ipos = 0; ipos < nCaloCellsPerCluster; ipos++)
 749     {
 750       Int_t absId  = clus->GetCellsAbsId()[ipos];
 751       if(absId == absIdMax || cells->GetCellAmplitude(absIdMax) < 0.01) continue;
 752
 753       Float_t frac = cells->GetCellAmplitude(absId)/cells->GetCellAmplitude(absIdMax);
 754       fhClusterMaxCellCloseCellRatio->Fill(clus->E(),frac);
 755       fhClusterMaxCellCloseCellDiff ->Fill(clus->E(),cells->GetCellAmplitude(absIdMax)-cells->GetCellAmplitude(absId));
 756
 757       if(fFillAllCellTimeHisto)
 758       {
 759         Double_t time  = cells->GetCellTime(absId);
 760         GetCaloUtils()->RecalibrateCellTime(time, fCalorimeter, absId,GetReader()->GetInputEvent()->GetBunchCrossNumber());
 761
 762         Float_t diff = (tmax-time*1.0e9);
 763         fhCellTimeSpreadRespectToCellMax->Fill(clus->E(), diff);
 764         if(TMath::Abs(TMath::Abs(diff) > 100) && clus->E() > 1 ) fhCellIdCellLargeTimeSpread->Fill(absId);
 765       }
 766
 767     }// fill cell-cluster histogram loop
 768
 769   }//check time and energy of cells respect to max energy cell if cluster of more than 1 cell
 770
 771
 772   // Get vertex for photon momentum calculation and event selection
 773   Double_t v[3] = {0,0,0}; //vertex ;
 774   //GetReader()->GetVertex(v); //
 775
 776   TLorentzVector mom  ;
 777   clus->GetMomentum(mom,v);
 778
 779   Float_t e   = mom.E();
 780   Float_t pt  = mom.Pt();
 781   Float_t eta = mom.Eta();
 782   Float_t phi = mom.Phi();
 783   if(phi < 0) phi +=TMath::TwoPi();
 784
 785   if(GetDebug() > 0) {
 786     printf("AliAnaCalorimeterQA::ClusterHistograms() - cluster: E %2.3f, pT %2.3f, eta %2.3f, phi %2.3f \n",e,pt,eta,phi*TMath::RadToDeg());
 787   }
 788
 789   fhE     ->Fill(e);
 790   if(nModule >=0 && nModule < fNModules) fhEMod->Fill(e,nModule);
 791   if(fFillAllTH12){
 792     fhPt     ->Fill(pt);
 793     fhPhi    ->Fill(phi);
 794     fhEta    ->Fill(eta);
 795   }
 796
 797   if(fFillAllTH3)
 798     fhEtaPhiE->Fill(eta,phi,e);
 799
 800   //Cells per cluster
 801   fhNCellsPerCluster   ->Fill(e, nCaloCellsPerCluster);
 802
 803   //Position
 804   if(fFillAllPosHisto2){
 805
 806     Float_t pos[3] ;
 807     clus->GetPosition(pos);
 808
 809     fhXE     ->Fill(pos[0],e);
 810     fhYE     ->Fill(pos[1],e);
 811     fhZE     ->Fill(pos[2],e);
 812     if(fFillAllTH3)
 813       fhXYZ    ->Fill(pos[0], pos[1],pos[2]);
 814
 815     fhXNCells->Fill(pos[0],nCaloCellsPerCluster);
 816     fhYNCells->Fill(pos[1],nCaloCellsPerCluster);
 817     fhZNCells->Fill(pos[2],nCaloCellsPerCluster);
 818     Float_t rxyz = TMath::Sqrt(pos[0]*pos[0]+pos[1]*pos[1]);//+pos[2]*pos[2]);
 819     fhRE     ->Fill(rxyz,e);
 820     fhRNCells->Fill(rxyz  ,nCaloCellsPerCluster);
 821   }
 822
 823   if(nModule >=0 && nModule < fNModules) fhNCellsPerClusterMod[nModule]->Fill(e, nCaloCellsPerCluster);
 824
 825 }
 826
 827 //____________________________________________________________________________
 828 void AliAnaCalorimeterQA::ClusterLoopHistograms(const TObjArray *caloClusters,
 829                                                 AliVCaloCells* cells)
 830 {
 831   // Fill clusters related histograms
 832
 833   TLorentzVector mom  ;
 834   Int_t  nLabel                = 0  ;
 835   Int_t *labels                = 0x0;
 836   Int_t  nCaloClusters         = caloClusters->GetEntriesFast() ;
 837   Int_t  nCaloClustersAccepted = 0  ;
 838   Int_t  nCaloCellsPerCluster  = 0  ;
 839   Bool_t matched               = kFALSE;
 840   Int_t  nModule               =-1  ;
 841
 842   // Get vertex for photon momentum calculation and event selection
 843   Double_t v[3] = {0,0,0}; //vertex ;
 844   //GetReader()->GetVertex(v);
 845
 846   Int_t *nClustersInModule     = new Int_t[fNModules];
 847   for(Int_t imod = 0; imod < fNModules; imod++ ) nClustersInModule[imod] = 0;
 848
 849   if(GetDebug() > 0)
 850     printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - In %s there are %d clusters \n", fCalorimeter.Data(), nCaloClusters);
 851
 852   // Loop over CaloClusters
 853   for(Int_t iclus = 0; iclus < nCaloClusters; iclus++){
 854
 855     if(GetDebug() > 0)
 856       printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - cluster: %d/%d, data %d \n",
 857              iclus+1,nCaloClusters,GetReader()->GetDataType());
 858
 859     AliVCluster* clus =  (AliVCluster*)caloClusters->At(iclus);
 860
 861     // Get the fraction of the cluster energy that carries the cell with highest energy and its absId
 862     Float_t maxCellFraction = 0.;
 863     Int_t absIdMax = GetCaloUtils()->GetMaxEnergyCell(cells, clus,maxCellFraction);
 864
 865     //Cut on time of clusters
 866     Double_t tof = clus->GetTOF()*1.e9;
 867     if(tof < fTimeCutMin || tof > fTimeCutMax)
 868     {
 869       if(GetDebug() > 0 )printf("AliAnaCalorimeterQA - Remove cluster with TOF %f\n",tof);
 870       continue;
 871     }
 872
 873     // Get cluster kinematics
 874     clus->GetMomentum(mom,v);
 875
 876     // Check only certain regions
 877     Bool_t in = kTRUE;
 878     if(IsFiducialCutOn()) in =  GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter) ;
 879     if(!in) continue;
 880
 881     // MC labels
 882     nLabel = clus->GetNLabels();
 883     labels = clus->GetLabels();
 884
 885     // SuperModule number of cluster
 886     nModule = GetModuleNumber(clus);
 887
 888     // Cells per cluster
 889     nCaloCellsPerCluster = clus->GetNCells();
 890
 891     // Cluster mathed with track?
 892     matched = GetCaloPID()->IsTrackMatched(clus,GetCaloUtils(), GetReader()->GetInputEvent());
 893
 894     // Get some time averages
 895     Double_t averTime[4] = {0.,0.,0.,0.};
 896     CalculateAverageTime(clus, cells, averTime);
 897
 898     //Get time of max cell
 899     Double_t tmax  = cells->GetCellTime(absIdMax);
 900     GetCaloUtils()->RecalibrateCellTime(tmax, fCalorimeter, absIdMax,GetReader()->GetInputEvent()->GetBunchCrossNumber());
 901     tmax*=1.e9;
 902
 903     // Fill histograms related to single cluster
 904
 905
 906     // Fill some histograms before applying the exotic cell / bad map cut
 907     fhNCellsPerClusterNoCut  ->Fill(clus->E(), nCaloCellsPerCluster);
 908     if(nModule >=0 && nModule < fNModules) fhNCellsPerClusterModNoCut[nModule]->Fill(clus->E(), nCaloCellsPerCluster);
 909
 910     fhClusterMaxCellDiffNoCut->Fill(clus->E(),maxCellFraction);
 911
 912     Float_t ampMax = cells->GetCellAmplitude(absIdMax);
 913     GetCaloUtils()->RecalibrateCellAmplitude(ampMax,fCalorimeter, absIdMax);
 914
 915     if(fStudyExotic) ExoticHistograms(absIdMax, ampMax, clus, cells);
 916
 917     //Check bad clusters if requested and rejection was not on
 918     Bool_t goodCluster = IsGoodCluster(absIdMax, cells);
 919
 920     Float_t eCrossFrac = 0;
 921     if(ampMax > 0.01) eCrossFrac = 1-GetECross(absIdMax,cells)/ampMax;
 922
 923     if(!goodCluster)
 924     {
 925       BadClusterHistograms(clus, caloClusters, cells, absIdMax,
 926                            maxCellFraction, eCrossFrac, tmax, averTime);
 927       continue;
 928     }
 929
 930     ClusterHistograms(clus, caloClusters, cells, absIdMax,
 931                       maxCellFraction, eCrossFrac, tmax, averTime);
 932
 933     nCaloClustersAccepted++;
 934     nModule = GetModuleNumber(clus);
 935     if(nModule >=0 && nModule < fNModules)
 936     {
 937       if     (fCalorimeter=="EMCAL" && mom.E() > 2*fEMCALCellAmpMin)  nClustersInModule[nModule]++;
 938       else if(fCalorimeter=="PHOS"  && mom.E() > 2*fPHOSCellAmpMin )  nClustersInModule[nModule]++;
 939     }
 940
 941     // Cluster weights
 942     if(fStudyWeight) WeightHistograms(clus, cells);
 943
 944     // Cells in cluster position
 945     if(fFillAllPosHisto) CellInClusterPositionHistograms(clus);
 946
 947     // Fill histograms related to single cluster, mc vs data
 948     Int_t  mcOK = kFALSE;
 949     Int_t  pdg  = -1;
 950     if(IsDataMC() && nLabel > 0 && labels)
 951       mcOK = ClusterMCHistograms(mom, matched, labels, nLabel, pdg);
 952
 953     // Matched clusters with tracks, also do some MC comparison, needs input from ClusterMCHistograms
 954     if( matched &&  fFillAllTMHisto)
 955       ClusterMatchedWithTrackHistograms(clus,mom,mcOK,pdg);
 956
 957     // Invariant mass
 958     // Try to reduce background with a mild shower shape cut and no more than 1 maxima
 959     // in cluster and remove low energy clusters
 960     if(fFillAllPi0Histo && nCaloClusters > 1 && nCaloCellsPerCluster > 1 &&
 961        GetCaloUtils()->GetNumberOfLocalMaxima(clus,cells) == 1 &&
 962        clus->GetM02() < 0.5 && clus->E() > 0.3)
 963       InvariantMassHistograms(iclus, mom, nModule, caloClusters,cells);
 964
 965   }//cluster loop
 966
 967   // Number of clusters histograms
 968   if(nCaloClustersAccepted > 0) fhNClusters->Fill(nCaloClustersAccepted);
 969
 970   // Number of clusters per module
 971   for(Int_t imod = 0; imod < fNModules; imod++ )
 972   {
 973     if(GetDebug() > 1)
 974       printf("AliAnaCalorimeterQA::ClusterLoopHistograms() - module %d calo %s clusters %d\n", imod, fCalorimeter.Data(), nClustersInModule[imod]);
 975     fhNClustersMod->Fill(nClustersInModule[imod],imod);
 976   }
 977
 978   delete [] nClustersInModule;
 979
 980 }
 981
 982 //__________________________________________________________________________________________
 983 Bool_t AliAnaCalorimeterQA::ClusterMCHistograms(TLorentzVector mom, Bool_t matched,
 984                                                 const Int_t * labels, Int_t nLabels, Int_t & pdg )
 985 {
 986
 987   //Fill histograms only possible when simulation
 988
 989   if(!labels || nLabels<=0)
 990   {
 991     if(GetDebug() > 1) printf("AliAnaCalorimeterQA::ClusterMCHistograms() - Strange, labels array %p, n labels %d \n", labels,nLabels);
 992     return kFALSE;
 993   }
 994
 995   if(GetDebug() > 1)
 996   {
 997     printf("AliAnaCalorimeterQA::ClusterMCHistograms() - Primaries: nlabels %d\n",nLabels);
 998   }
 999
1000   Float_t e   = mom.E();
1001   Float_t eta = mom.Eta();
1002   Float_t phi = mom.Phi();
1003   if(phi < 0) phi +=TMath::TwoPi();
1004
1005   AliAODMCParticle * aodprimary  = 0x0;
1006   TParticle * primary = 0x0;
1007
1008   //Play with the MC stack if available
1009   Int_t label = labels[0];
1010
1011   if(label < 0)
1012   {
1013     if(GetDebug() >= 0) printf("AliAnaCalorimeterQA::ClusterMCHistograms() *** bad label ***:  label %d \n", label);
1014     return kFALSE;
1015   }
1016
1017   Int_t pdg0  =-1;Int_t status = -1; Int_t iMother = -1; Int_t iParent = -1;
1018   Float_t vxMC= 0; Float_t vyMC = 0;
1019   Float_t eMC = 0; //Float_t ptMC= 0;
1020   Float_t phiMC =0; Float_t etaMC = 0;
1021   Int_t charge = 0;
1022
1023   //Check the origin.
1024   Int_t tag = GetMCAnalysisUtils()->CheckOrigin(labels,nLabels, GetReader());
1025
1026   if     ( GetReader()->ReadStack() &&
1027           !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCUnknown))
1028   { //it MC stack and known tag
1029
1030     if( label >= GetMCStack()->GetNtrack())
1031     {
1032       if(GetDebug() >= 0) printf("AliAnaCalorimeterQA::ClusterMCHistograms() *** large label ***:  label %d, n tracks %d \n", label, GetMCStack()->GetNtrack());
1033       return kFALSE;
1034     }
1035
1036     primary = GetMCStack()->Particle(label);
1037     iMother = label;
1038     pdg0    = TMath::Abs(primary->GetPdgCode());
1039     pdg     = pdg0;
1040     status  = primary->GetStatusCode();
1041     vxMC    = primary->Vx();
1042     vyMC    = primary->Vy();
1043     iParent = primary->GetFirstMother();
1044
1045     if(GetDebug() > 1 )
1046     {
1047       printf("AliAnaCalorimeterQA::ClusterMCHistograms() - Cluster most contributing mother: \n");
1048       printf("\t Mother label %d, pdg %d, %s, status %d, parent %d \n",iMother, pdg0, primary->GetName(),status, iParent);
1049     }
1050
1051     //Get final particle, no conversion products
1052     if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion))
1053     {
1054       //Get the parent
1055       primary = GetMCStack()->Particle(iParent);
1056       pdg = TMath::Abs(primary->GetPdgCode());
1057
1058       if(GetDebug() > 1 ) printf("AliAnaCalorimeterQA::ClusterMCHistograms() - Converted cluster!. Find before conversion: \n");
1059
1060       while((pdg == 22 || pdg == 11) && status != 1)
1061       {
1062         Int_t iMotherOrg = iMother;
1063         iMother = iParent;
1064         primary = GetMCStack()->Particle(iMother);
1065         status  = primary->GetStatusCode();
1066         pdg     = TMath::Abs(primary->GetPdgCode());
1067         iParent = primary->GetFirstMother();
1068
1069         // If gone too back and non stable, assign the decay photon/electron
1070         // there are other possible decays, ignore them for the moment
1071         if(pdg==111 || pdg==221)
1072         {
1073           primary = GetMCStack()->Particle(iMotherOrg);
1074           break;
1075         }
1076
1077         if( iParent < 0 )
1078         {
1079           iParent = iMother;
1080           break;
1081         }
1082
1083         if(GetDebug() > 1 )printf("\t pdg %d, index %d, %s, status %d \n",pdg, iMother,  primary->GetName(),status);
1084       }
1085
1086       if(GetDebug() > 1 )
1087       {
1088         printf("AliAnaCalorimeterQA::ClusterHistograms() - Converted Cluster mother before conversion: \n");
1089         printf("\t Mother label %d, pdg %d, %s, status %d, parent %d \n",iMother, pdg, primary->GetName(), status, iParent);
1090       }
1091
1092     }
1093
1094     //Overlapped pi0 (or eta, there will be very few), get the meson
1095     if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPi0) ||
1096        GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCEta))
1097     {
1098       if(GetDebug() > 1 ) printf("AliAnaCalorimeterQA::ClusterHistograms() - Overlapped Meson decay!, Find it: \n");
1099
1100       while(pdg != 111 && pdg != 221)
1101       {
1102         //printf("iMother %d, pdg %d, iParent %d, pdg %d\n",iMother,pdg,iParent,GetMCStack()->Particle(iParent)->GetPdgCode());
1103         iMother = iParent;
1104         primary = GetMCStack()->Particle(iMother);
1105         status  = primary->GetStatusCode();
1106         pdg     = TMath::Abs(primary->GetPdgCode());
1107         iParent = primary->GetFirstMother();
1108
1109         if( iParent < 0 )break;
1110
1111         if(GetDebug() > 1 ) printf("\t pdg %d, %s, index %d\n",pdg,  primary->GetName(),iMother);
1112
1113         if(iMother==-1)
1114         {
1115           printf("AliAnaCalorimeterQA::ClusterHistograms() - Tagged as Overlapped photon but meson not found, why?\n");
1116           //break;
1117         }
1118       }
1119
1120       if(GetDebug() > 2 ) printf("AliAnaCalorimeterQA::ClusterHistograms() - Overlapped %s decay, label %d \n",
1121                                  primary->GetName(),iMother);
1122     }
1123
1124     eMC    = primary->Energy();
1125     //ptMC   = primary->Pt();
1126     phiMC  = primary->Phi();
1127     etaMC  = primary->Eta();
1128     pdg    = TMath::Abs(primary->GetPdgCode());
1129     charge = (Int_t) TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
1130
1131   }
1132   else if( GetReader()->ReadAODMCParticles() &&
1133           !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCUnknown))
1134   {//it MC AOD and known tag
1135     //Get the list of MC particles
1136     if(!GetReader()->GetAODMCParticles())
1137       AliFatal("MCParticles not available!");
1138
1139     aodprimary = (AliAODMCParticle*) (GetReader()->GetAODMCParticles())->At(label);
1140     iMother = label;
1141     pdg0    = TMath::Abs(aodprimary->GetPdgCode());
1142     pdg     = pdg0;
1143     status  = aodprimary->IsPrimary();
1144     vxMC    = aodprimary->Xv();
1145     vyMC    = aodprimary->Yv();
1146     iParent = aodprimary->GetMother();
1147
1148     if(GetDebug() > 1 )
1149     {
1150       printf("AliAnaCalorimeterQA::ClusterHistograms() - Cluster most contributing mother: \n");
1151       printf("\t Mother label %d, pdg %d, Primary? %d, Physical Primary? %d, parent %d \n",
1152              iMother, pdg0, aodprimary->IsPrimary(), aodprimary->IsPhysicalPrimary(), iParent);
1153     }
1154
1155     //Get final particle, no conversion products
1156     if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion))
1157     {
1158       if(GetDebug() > 1 )
1159         printf("AliAnaCalorimeterQA::ClusterHistograms() - Converted cluster!. Find before conversion: \n");
1160       //Get the parent
1161       aodprimary = (AliAODMCParticle*)(GetReader()->GetAODMCParticles())->At(iParent);
1162       pdg = TMath::Abs(aodprimary->GetPdgCode());
1163       while ((pdg == 22 || pdg == 11) && !aodprimary->IsPhysicalPrimary())
1164       {
1165         Int_t iMotherOrg = iMother;
1166         iMother    = iParent;
1167         aodprimary = (AliAODMCParticle*)(GetReader()->GetAODMCParticles())->At(iMother);
1168         status     = aodprimary->IsPrimary();
1169         iParent    = aodprimary->GetMother();
1170         pdg        = TMath::Abs(aodprimary->GetPdgCode());
1171
1172         // If gone too back and non stable, assign the decay photon/electron
1173         // there are other possible decays, ignore them for the moment
1174         if(pdg==111 || pdg==221)
1175         {
1176           aodprimary = (AliAODMCParticle*)(GetReader()->GetAODMCParticles())->At(iMotherOrg);
1177           break;
1178         }
1179
1180         if(iParent < 0 )
1181         {
1182           iParent = iMother;
1183           break;
1184         }
1185
1186         if(GetDebug() > 1 )
1187           printf("\t pdg %d, index %d, Primary? %d, Physical Primary? %d \n",
1188                  pdg, iMother, aodprimary->IsPrimary(), aodprimary->IsPhysicalPrimary());
1189       }
1190
1191       if(GetDebug() > 1 )
1192       {
1193         printf("AliAnaCalorimeterQA::ClusterHistograms() - Converted Cluster mother before conversion: \n");
1194         printf("\t Mother label %d, pdg %d, parent %d, Primary? %d, Physical Primary? %d \n",
1195                iMother, pdg, iParent, aodprimary->IsPrimary(), aodprimary->IsPhysicalPrimary());
1196       }
1197
1198     }
1199
1200     //Overlapped pi0 (or eta, there will be very few), get the meson
1201     if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPi0) ||
1202        GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCEta))
1203     {
1204       if(GetDebug() > 1 ) printf("AliAnaCalorimeterQA::ClusterHistograms() - Overlapped Meson decay!, Find it: PDG %d, mom %d \n",pdg, iMother);
1205       while(pdg != 111 && pdg != 221)
1206       {
1207         iMother    = iParent;
1208         aodprimary = (AliAODMCParticle*)(GetReader()->GetAODMCParticles())->At(iMother);
1209         status     = aodprimary->IsPrimary();
1210         iParent    = aodprimary->GetMother();
1211         pdg        = TMath::Abs(aodprimary->GetPdgCode());
1212
1213         if(iParent < 0 ) break;
1214
1215         if(GetDebug() > 1 ) printf("\t pdg %d, index %d\n",pdg, iMother);
1216
1217         if(iMother==-1)
1218         {
1219           printf("AliAnaCalorimeterQA::ClusterHistograms() - Tagged as Overlapped photon but meson not found, why?\n");
1220           //break;
1221         }
1222       }
1223
1224       if(GetDebug() > 2 ) printf("AliAnaCalorimeterQA::ClusterHistograms() - Overlapped %s decay, label %d \n",
1225                                  aodprimary->GetName(),iMother);
1226     }
1227
1228     status = aodprimary->IsPrimary();
1229     eMC    = aodprimary->E();
1230     //ptMC   = aodprimary->Pt();
1231     phiMC  = aodprimary->Phi();
1232     etaMC  = aodprimary->Eta();
1233     pdg    = TMath::Abs(aodprimary->GetPdgCode());
1234     charge = aodprimary->Charge();
1235
1236   }
1237
1238   //Float_t vz = primary->Vz();
1239   Float_t rVMC = TMath::Sqrt(vxMC*vxMC + vyMC*vyMC);
1240   if((pdg == 22 || TMath::Abs(pdg)==11) && status!=1)
1241   {
1242     fhEMVxyz   ->Fill(vxMC,vyMC);//,vz);
1243     fhEMR      ->Fill(e,rVMC);
1244   }
1245
1246   //printf("reco e %f, pt %f, phi %f, eta %f \n", e, pt, phi, eta);
1247   //printf("prim e %f, pt %f, phi %f, eta %f \n", eMC,ptMC,phiMC ,etaMC );
1248   //printf("vertex: vx %f, vy %f, vz %f, r %f \n", vxMC, vyMC, vz, r);
1249
1250   //Overlapped pi0 (or eta, there will be very few)
1251   if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPi0))
1252   {
1253     fhRecoMCE  [kmcPi0][matched]     ->Fill(e,eMC);
1254     if(e > 0.5 && eMC > 0.5) fhRecoMCEta[kmcPi0][(matched)]->Fill(eta,etaMC);
1255     if(e > 0.5 && eMC > 0.5) fhRecoMCPhi[kmcPi0][(matched)]->Fill(phi,phiMC);
1256     if(eMC > 0) fhRecoMCRatioE  [kmcPi0][(matched)]->Fill(e,e/eMC);
1257     fhRecoMCDeltaE  [kmcPi0][(matched)]->Fill(e,eMC-e);
1258     fhRecoMCDeltaPhi[kmcPi0][(matched)]->Fill(e,phiMC-phi);
1259     fhRecoMCDeltaEta[kmcPi0][(matched)]->Fill(e,etaMC-eta);
1260   }//Overlapped pizero decay
1261   else     if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCEta))
1262   {
1263     fhRecoMCE  [kmcEta][(matched)]     ->Fill(e,eMC);
1264     if(e > 0.5 && eMC > 0.5) fhRecoMCEta[kmcEta][(matched)]->Fill(eta,etaMC);
1265     if(e > 0.5 && eMC > 0.5) fhRecoMCPhi[kmcEta][(matched)]->Fill(phi,phiMC);
1266     if(eMC > 0) fhRecoMCRatioE  [kmcEta][(matched)]->Fill(e,e/eMC);
1267     fhRecoMCDeltaE  [kmcEta][(matched)]->Fill(e,eMC-e);
1268     fhRecoMCDeltaPhi[kmcEta][(matched)]->Fill(e,phiMC-phi);
1269     fhRecoMCDeltaEta[kmcEta][(matched)]->Fill(e,etaMC-eta);
1270   }//Overlapped eta decay
1271   else if( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPhoton) &&
1272           !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion))
1273   {
1274     fhRecoMCE  [kmcPhoton][(matched)]     ->Fill(e,eMC);
1275     if(e > 0.5 && eMC > 0.5) fhRecoMCEta[kmcPhoton][(matched)]->Fill(eta,etaMC);
1276     if(e > 0.5 && eMC > 0.5) fhRecoMCPhi[kmcPhoton][(matched)]->Fill(phi,phiMC);
1277     if(eMC > 0) fhRecoMCRatioE  [kmcPhoton][(matched)]->Fill(e,e/eMC);
1278
1279     fhRecoMCDeltaE  [kmcPhoton][(matched)]->Fill(e,eMC-e);
1280     fhRecoMCDeltaPhi[kmcPhoton][(matched)]->Fill(e,phiMC-phi);
1281     fhRecoMCDeltaEta[kmcPhoton][(matched)]->Fill(e,etaMC-eta);
1282   }//photon
1283   else if( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCElectron) &&
1284           !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion))
1285   {
1286     fhRecoMCE  [kmcElectron][(matched)]     ->Fill(e,eMC);
1287     if(e > 0.5 && eMC > 0.5) fhRecoMCEta[kmcElectron][(matched)]->Fill(eta,etaMC);
1288     if(e > 0.5 && eMC > 0.5) fhRecoMCPhi[kmcElectron][(matched)]->Fill(phi,phiMC);
1289     if(eMC > 0) fhRecoMCRatioE  [kmcElectron][(matched)]->Fill(e,e/eMC);
1290     fhRecoMCDeltaE  [kmcElectron][(matched)]->Fill(e,eMC-e);
1291     fhRecoMCDeltaPhi[kmcElectron][(matched)]->Fill(e,phiMC-phi);
1292     fhRecoMCDeltaEta[kmcElectron][(matched)]->Fill(e,etaMC-eta);
1293     fhEMVxyz   ->Fill(vxMC,vyMC);//,vz);
1294     fhEMR      ->Fill(e,rVMC);
1295   }
1296   else if(charge == 0)
1297   {
1298     fhRecoMCE  [kmcNeHadron][(matched)]     ->Fill(e,eMC);
1299     if(e > 0.5 && eMC > 0.5) fhRecoMCEta[kmcNeHadron][(matched)]->Fill(eta,etaMC);
1300     if(e > 0.5 && eMC > 0.5) fhRecoMCPhi[kmcNeHadron][(matched)]->Fill(phi,phiMC);
1301     if(eMC > 0) fhRecoMCRatioE  [kmcNeHadron][(matched)]->Fill(e,e/eMC);
1302     fhRecoMCDeltaE  [kmcNeHadron][(matched)]->Fill(e,eMC-e);
1303     fhRecoMCDeltaPhi[kmcNeHadron][(matched)]->Fill(e,phiMC-phi);
1304     fhRecoMCDeltaEta[kmcNeHadron][(matched)]->Fill(e,etaMC-eta);
1305     fhHaVxyz     ->Fill(vxMC,vyMC);//,vz);
1306     fhHaR        ->Fill(e,rVMC);
1307   }
1308   else if(charge!=0)
1309   {
1310     fhRecoMCE  [kmcChHadron][(matched)]     ->Fill(e,eMC);
1311     if(e > 0.5 && eMC > 0.5) fhRecoMCEta[kmcChHadron][(matched)]->Fill(eta,etaMC);
1312     if(e > 0.5 && eMC > 0.5) fhRecoMCPhi[kmcChHadron][(matched)]->Fill(phi,phiMC);
1313     if(eMC > 0) fhRecoMCRatioE  [kmcChHadron][(matched)]->Fill(e,e/eMC);
1314     fhRecoMCDeltaE  [kmcChHadron][(matched)]->Fill(e,eMC-e);
1315     fhRecoMCDeltaPhi[kmcChHadron][(matched)]->Fill(e,phiMC-phi);
1316     fhRecoMCDeltaEta[kmcChHadron][(matched)]->Fill(e,etaMC-eta);
1317     fhHaVxyz     ->Fill(vxMC,vyMC);//,vz);
1318     fhHaR        ->Fill(e,rVMC);
1319   }
1320
1321   if(primary || aodprimary) return kTRUE ;
1322   else                      return kFALSE;
1323
1324 }
1325
1326 //________________________________________________________________________________________________
1327 void AliAnaCalorimeterQA::ClusterMatchedWithTrackHistograms(AliVCluster *clus, TLorentzVector mom,
1328                                                             Bool_t okPrimary, Int_t pdg)
1329 {
1330   //Histograms for clusters matched with tracks
1331
1332   Float_t e   = mom.E();
1333   Float_t pt  = mom.Pt();
1334   Float_t eta = mom.Eta();
1335   Float_t phi = mom.Phi();
1336   if(phi < 0) phi +=TMath::TwoPi();
1337
1338   if(fFillAllTH12)
1339   {
1340     fhECharged      ->Fill(e);
1341     fhPtCharged     ->Fill(pt);
1342     fhPhiCharged    ->Fill(phi);
1343     fhEtaCharged    ->Fill(eta);
1344   }
1345
1346   //Study the track and matched cluster if track exists.
1347
1348   AliVTrack *track = GetCaloUtils()->GetMatchedTrack(clus, GetReader()->GetInputEvent());
1349
1350   if(!track) return ;
1351
1352   Double_t tpt   = track->Pt();
1353   Double_t tmom  = track->P();
1354   Double_t dedx  = track->GetTPCsignal();
1355   Int_t    nITS  = track->GetNcls(0);
1356   Int_t    nTPC  = track->GetNcls(1);
1357   Bool_t positive = kFALSE;
1358   if(track) positive = (track->Charge()>0);
1359
1360   // Residuals
1361   Float_t deta  = clus->GetTrackDz();
1362   Float_t dphi  = clus->GetTrackDx();
1363   Double_t  dR  = TMath::Sqrt(dphi*dphi + deta*deta);
1364
1365   if(TMath::Abs(dphi) < 999)
1366   {
1367     fhTrackMatchedDEta->Fill(e,deta);
1368     fhTrackMatchedDPhi->Fill(e,dphi);
1369     if(e > 0.5) fhTrackMatchedDEtaDPhi->Fill(deta,dphi);
1370
1371     if(track && positive)
1372     {
1373
1374       fhTrackMatchedDEtaPos->Fill(e,deta);
1375       fhTrackMatchedDPhiPos->Fill(e,dphi);
1376       if(e > 0.5) fhTrackMatchedDEtaDPhiPos->Fill(deta,dphi);
1377     }
1378   }
1379
1380   Double_t eOverP = e/tmom;
1381   fh1EOverP->Fill(tpt, eOverP);
1382   if(dR < 0.02)
1383   {
1384     fh1EOverPR02->Fill(tpt,eOverP);
1385     if(dedx > 60 && dedx < 100) fh1EleEOverP->Fill(tpt,eOverP);
1386   }
1387
1388   fh2dR->Fill(e,dR);
1389   fh2MatchdEdx->Fill(tmom,dedx);
1390
1391   if(IsDataMC() && okPrimary)
1392   {
1393     Double_t  charge = TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
1394
1395     if(TMath::Abs(pdg) == 11)
1396     {
1397       fhMCEle1EOverP->Fill(tpt,eOverP);
1398       fhMCEle1dR->Fill(dR);
1399       fhMCEle2MatchdEdx->Fill(tmom,dedx);
1400       if(dR < 0.02)
1401       {
1402         fhMCEle1EOverPR02->Fill(tpt,eOverP);
1403         if(dedx > 60 && dedx < 100) fhMCEle1EleEOverP->Fill(tpt,eOverP);
1404       }
1405     }
1406     else if(charge!=0)
1407     {
1408       fhMCChHad1EOverP->Fill(tpt,eOverP);
1409       fhMCChHad1dR->Fill(dR);
1410       fhMCChHad2MatchdEdx->Fill(tmom,dedx);
1411       if(dR < 0.02)
1412       {
1413         fhMCChHad1EOverPR02->Fill(tpt,eOverP);
1414         if(dedx > 60 && dedx < 100) fhMCChHad1EleEOverP->Fill(tpt,eOverP);
1415       }
1416     }
1417     else if(charge == 0)
1418     {
1419       fhMCNeutral1EOverP->Fill(tpt,eOverP);
1420       fhMCNeutral1dR->Fill(dR);
1421       fhMCNeutral2MatchdEdx->Fill(tmom,dedx);
1422       if(dR < 0.02)
1423       {
1424         fhMCNeutral1EOverPR02->Fill(tpt,eOverP);
1425         if(dedx > 60 && dedx < 100) fhMCNeutral1EleEOverP->Fill(tpt,eOverP);
1426       }
1427     }
1428   }//DataMC
1429
1430   if(dR < 0.02 && eOverP > 0.6 && eOverP< 1.2
1431      && clus->GetNCells() > 1 && nITS > 3 && nTPC > 20)
1432   {
1433     fh2EledEdx->Fill(tmom,dedx);
1434   }
1435
1436 }
1437
1438 //___________________________________
1439 void AliAnaCalorimeterQA::Correlate()
1440 {
1441   // Correlate information from PHOS and EMCAL and with V0 and track multiplicity
1442
1443   //Clusters
1444   TObjArray * caloClustersEMCAL = GetEMCALClusters();
1445   TObjArray * caloClustersPHOS  = GetPHOSClusters();
1446
1447   Int_t nclEMCAL = caloClustersEMCAL->GetEntriesFast();
1448   Int_t nclPHOS  = caloClustersPHOS ->GetEntriesFast();
1449
1450   Float_t cen = GetEventCentrality();
1451   Float_t ep  = GetEventPlaneAngle();
1452
1453   Float_t sumClusterEnergyEMCAL = 0;
1454   Float_t sumClusterEnergyPHOS  = 0;
1455   Int_t iclus = 0;
1456   for(iclus = 0 ; iclus <  caloClustersEMCAL->GetEntriesFast() ; iclus++)
1457     sumClusterEnergyEMCAL += ((AliVCluster*)caloClustersEMCAL->At(iclus))->E();
1458   for(iclus = 0 ; iclus <  caloClustersPHOS->GetEntriesFast(); iclus++)
1459     sumClusterEnergyPHOS += ((AliVCluster*)caloClustersPHOS->At(iclus))->E();
1460
1461   //Cells
1462
1463   AliVCaloCells * cellsEMCAL = GetEMCALCells();
1464   AliVCaloCells * cellsPHOS  = GetPHOSCells();
1465
1466   Int_t ncellsEMCAL = cellsEMCAL->GetNumberOfCells();
1467   Int_t ncellsPHOS  = cellsPHOS ->GetNumberOfCells();
1468
1469   Float_t sumCellEnergyEMCAL = 0;
1470   Float_t sumCellEnergyPHOS  = 0;
1471   Int_t icell = 0;
1472   for(icell = 0 ; icell < cellsEMCAL->GetNumberOfCells()  ; icell++)
1473     sumCellEnergyEMCAL += cellsEMCAL->GetAmplitude(icell);
1474   for(icell = 0 ; icell <  cellsPHOS->GetNumberOfCells(); icell++)
1475     sumCellEnergyPHOS += cellsPHOS->GetAmplitude(icell);
1476
1477
1478   //Fill Histograms
1479   fhCaloCorrNClusters->Fill(nclEMCAL,nclPHOS);
1480   fhCaloCorrEClusters->Fill(sumClusterEnergyEMCAL,sumClusterEnergyPHOS);
1481   fhCaloCorrNCells   ->Fill(ncellsEMCAL,ncellsPHOS);
1482   fhCaloCorrECells   ->Fill(sumCellEnergyEMCAL,sumCellEnergyPHOS);
1483
1484   Int_t v0S = GetV0Signal(0)+GetV0Signal(1);
1485   Int_t v0M = GetV0Multiplicity(0)+GetV0Multiplicity(1);
1486   Int_t trM = GetTrackMultiplicity();
1487   if(fCalorimeter=="PHOS")
1488   {
1489     fhCaloV0MCorrNClusters   ->Fill(v0M,nclPHOS);
1490     fhCaloV0MCorrEClusters   ->Fill(v0M,sumClusterEnergyPHOS);
1491     fhCaloV0MCorrNCells      ->Fill(v0M,ncellsPHOS);
1492     fhCaloV0MCorrECells      ->Fill(v0M,sumCellEnergyPHOS);
1493
1494     fhCaloV0SCorrNClusters   ->Fill(v0S,nclPHOS);
1495     fhCaloV0SCorrEClusters   ->Fill(v0S,sumClusterEnergyPHOS);
1496     fhCaloV0SCorrNCells      ->Fill(v0S,ncellsPHOS);
1497     fhCaloV0SCorrECells      ->Fill(v0S,sumCellEnergyPHOS);
1498
1499     fhCaloTrackMCorrNClusters->Fill(trM,nclPHOS);
1500     fhCaloTrackMCorrEClusters->Fill(trM,sumClusterEnergyPHOS);
1501     fhCaloTrackMCorrNCells   ->Fill(trM,ncellsPHOS);
1502     fhCaloTrackMCorrECells   ->Fill(trM,sumCellEnergyPHOS);
1503
1504     fhCaloCenNClusters       ->Fill(cen,nclPHOS);
1505     fhCaloCenEClusters       ->Fill(cen,sumClusterEnergyPHOS);
1506     fhCaloCenNCells          ->Fill(cen,ncellsPHOS);
1507     fhCaloCenECells          ->Fill(cen,sumCellEnergyPHOS);
1508
1509     fhCaloEvPNClusters       ->Fill(ep ,nclPHOS);
1510     fhCaloEvPEClusters       ->Fill(ep ,sumClusterEnergyPHOS);
1511     fhCaloEvPNCells          ->Fill(ep ,ncellsPHOS);
1512     fhCaloEvPECells          ->Fill(ep ,sumCellEnergyPHOS);
1513   }
1514   else
1515   {
1516     fhCaloV0MCorrNClusters   ->Fill(v0M,nclEMCAL);
1517     fhCaloV0MCorrEClusters   ->Fill(v0M,sumClusterEnergyEMCAL);
1518     fhCaloV0MCorrNCells      ->Fill(v0M,ncellsEMCAL);
1519     fhCaloV0MCorrECells      ->Fill(v0M,sumCellEnergyEMCAL);
1520
1521     fhCaloV0SCorrNClusters   ->Fill(v0S,nclEMCAL);
1522     fhCaloV0SCorrEClusters   ->Fill(v0S,sumClusterEnergyEMCAL);
1523     fhCaloV0SCorrNCells      ->Fill(v0S,ncellsEMCAL);
1524     fhCaloV0SCorrECells      ->Fill(v0S,sumCellEnergyEMCAL);
1525
1526     fhCaloTrackMCorrNClusters->Fill(trM,nclEMCAL);
1527     fhCaloTrackMCorrEClusters->Fill(trM,sumClusterEnergyEMCAL);
1528     fhCaloTrackMCorrNCells   ->Fill(trM,ncellsEMCAL);
1529     fhCaloTrackMCorrECells   ->Fill(trM,sumCellEnergyEMCAL);
1530
1531     fhCaloCenNClusters       ->Fill(cen,nclEMCAL);
1532     fhCaloCenEClusters       ->Fill(cen,sumClusterEnergyEMCAL);
1533     fhCaloCenNCells          ->Fill(cen,ncellsEMCAL);
1534     fhCaloCenECells          ->Fill(cen,sumCellEnergyEMCAL);
1535
1536     fhCaloEvPNClusters       ->Fill(ep ,nclEMCAL);
1537     fhCaloEvPEClusters       ->Fill(ep ,sumClusterEnergyEMCAL);
1538     fhCaloEvPNCells          ->Fill(ep ,ncellsEMCAL);
1539     fhCaloEvPECells          ->Fill(ep ,sumCellEnergyEMCAL);
1540   }
1541
1542   if(GetDebug() > 0 )
1543   {
1544     printf("AliAnaCalorimeterQA::Correlate(): \n");
1545     printf("\t EMCAL: N cells %d, N clusters  %d, summed E cells %f, summed E clusters %f \n",
1546            ncellsEMCAL,nclEMCAL, sumCellEnergyEMCAL,sumClusterEnergyEMCAL);
1547     printf("\t PHOS : N cells %d, N clusters  %d, summed E cells %f, summed E clusters %f \n",
1548            ncellsPHOS,nclPHOS,sumCellEnergyPHOS,sumClusterEnergyPHOS);
1549     printf("\t V0 : Signal %d, Multiplicity  %d, Track Multiplicity %d \n", v0S,v0M,trM);
1550     printf("\t centrality : %f, Event plane angle %f \n", cen,ep);
1551   }
1552
1553 }
1554
1555 //__________________________________________________
1556 TObjString * AliAnaCalorimeterQA::GetAnalysisCuts()
1557 {
1558   //Save parameters used for analysis
1559   TString parList ; //this will be list of parameters used for this analysis.
1560   const Int_t buffersize = 255;
1561   char onePar[buffersize] ;
1562
1563   snprintf(onePar,buffersize,"--- AliAnaCalorimeterQA ---\n") ;
1564   parList+=onePar ;
1565   snprintf(onePar,buffersize,"Calorimeter: %s\n",fCalorimeter.Data()) ;
1566   parList+=onePar ;
1567   snprintf(onePar,buffersize,"Time Cut : %2.2f < T < %2.2f ns  \n",fTimeCutMin, fTimeCutMax) ;
1568   parList+=onePar ;
1569   snprintf(onePar,buffersize,"PHOS Cell Amplitude > %2.2f GeV, EMCAL Cell Amplitude > %2.2f GeV  \n",fPHOSCellAmpMin, fEMCALCellAmpMin) ;
1570   parList+=onePar ;
1571   //Get parameters set in base class.
1572   //parList += GetBaseParametersList() ;
1573
1574   //Get parameters set in FiducialCut class (not available yet)
1575   //parlist += GetFidCut()->GetFidCutParametersList()
1576
1577   return new TObjString(parList) ;
1578 }
1579
1580 //_________________________________________________________________________________
1581 void AliAnaCalorimeterQA::ExoticHistograms(Int_t absIdMax, Float_t ampMax,
1582                                            AliVCluster *clus, AliVCaloCells* cells)
1583 {
1584   // Calculate weights
1585
1586   if(ampMax < 0.01)
1587   {
1588     printf("AliAnaCalorimeterQA::ExoticHistograms()- Low amplitude energy %f\n",ampMax);
1589     return;
1590   }
1591
1592   Float_t  l0   = clus->GetM02();
1593   Float_t  l1   = clus->GetM20();
1594   Float_t  en   = clus->E();
1595   Int_t    nc   = clus->GetNCells();
1596   Double_t tmax = clus->GetTOF()*1.e9; // recalibrated elsewhere
1597
1598   Float_t eCrossFrac = 1-GetECross(absIdMax,cells, 10000000)/ampMax;
1599
1600   if(en > 5)
1601   {
1602     fhExoL0ECross->Fill(eCrossFrac,l0);
1603     fhExoL1ECross->Fill(eCrossFrac,l1);
1604   }
1605
1606   for(Int_t ie = 0; ie < fExoNECrossCuts; ie++)
1607   {
1608     for(Int_t idt = 0; idt < fExoNDTimeCuts; idt++)
1609     {
1610       eCrossFrac = 1-GetECross(absIdMax,cells, fExoDTimeCuts[idt])/ampMax;
1611
1612       if(eCrossFrac > fExoECrossCuts[ie])
1613       {
1614         //Exotic
1615         fhExoL0    [ie][idt]->Fill(en,l0  );
1616         fhExoL1    [ie][idt]->Fill(en,l1  );
1617         fhExoTime  [ie][idt]->Fill(en,tmax);
1618
1619         if(en > 5)
1620         {
1621           fhExoL0NCell[ie][idt]->Fill(nc,l0);
1622           fhExoL1NCell[ie][idt]->Fill(nc,l1);
1623         }
1624
1625         // Diff time, do for one cut in e cross
1626         if(ie == 0)
1627         {
1628           for (Int_t icell = 0; icell < clus->GetNCells(); icell++)
1629           {
1630             Int_t absId  = clus->GetCellsAbsId()[icell];
1631             Double_t time  = cells->GetCellTime(absId);
1632             GetCaloUtils()->RecalibrateCellTime(time, fCalorimeter, absId,GetReader()->GetInputEvent()->GetBunchCrossNumber());
1633
1634             Float_t diff = (tmax-time)*1e9;
1635             fhExoDTime[idt]->Fill(en, diff);
1636           }
1637         }
1638       }
1639       else
1640       {
1641         fhExoECross[ie][idt]->Fill(en,eCrossFrac);
1642         fhExoNCell [ie][idt]->Fill(en,nc);
1643       }
1644     } // D time cut loop
1645   } // e cross cut loop
1646 }
1647
1648 //____________________________________________________
1649 TList * AliAnaCalorimeterQA::GetCreateOutputObjects()
1650 {
1651   // Create histograms to be saved in output file and
1652   // store them in outputContainer
1653
1654   TList * outputContainer = new TList() ;
1655   outputContainer->SetName("QAHistos") ;
1656
1657   //Histograms
1658   Int_t nptbins     = GetHistogramRanges()->GetHistoPtBins();           Float_t ptmax     = GetHistogramRanges()->GetHistoPtMax();           Float_t ptmin     = GetHistogramRanges()->GetHistoPtMin();
1659   Int_t nfineptbins = GetHistogramRanges()->GetHistoFinePtBins();           Float_t ptfinemax = GetHistogramRanges()->GetHistoFinePtMax();       Float_t ptfinemin = GetHistogramRanges()->GetHistoFinePtMin();
1660   Int_t nphibins    = GetHistogramRanges()->GetHistoPhiBins();              Float_t phimax    = GetHistogramRanges()->GetHistoPhiMax();          Float_t phimin    = GetHistogramRanges()->GetHistoPhiMin();
1661   Int_t netabins    = GetHistogramRanges()->GetHistoEtaBins();          Float_t etamax    = GetHistogramRanges()->GetHistoEtaMax();          Float_t etamin    = GetHistogramRanges()->GetHistoEtaMin();
1662   Int_t nmassbins   = GetHistogramRanges()->GetHistoMassBins();         Float_t massmax   = GetHistogramRanges()->GetHistoMassMax();           Float_t massmin   = GetHistogramRanges()->GetHistoMassMin();
1663   Int_t nasymbins   = GetHistogramRanges()->GetHistoAsymmetryBins();    Float_t asymmax   = GetHistogramRanges()->GetHistoAsymmetryMax();    Float_t asymmin   = GetHistogramRanges()->GetHistoAsymmetryMin();
1664   Int_t nPoverEbins = GetHistogramRanges()->GetHistoPOverEBins();       Float_t eOverPmax = GetHistogramRanges()->GetHistoPOverEMax();       Float_t eOverPmin = GetHistogramRanges()->GetHistoPOverEMin();
1665   Int_t ndedxbins   = GetHistogramRanges()->GetHistodEdxBins();         Float_t dedxmax   = GetHistogramRanges()->GetHistodEdxMax();         Float_t dedxmin   = GetHistogramRanges()->GetHistodEdxMin();
1666   Int_t ndRbins     = GetHistogramRanges()->GetHistodRBins();           Float_t dRmax     = GetHistogramRanges()->GetHistodRMax();           Float_t dRmin     = GetHistogramRanges()->GetHistodRMin();
1667   Int_t ntimebins   = GetHistogramRanges()->GetHistoTimeBins();         Float_t timemax   = GetHistogramRanges()->GetHistoTimeMax();         Float_t timemin   = GetHistogramRanges()->GetHistoTimeMin();
1668   Int_t nclbins     = GetHistogramRanges()->GetHistoNClustersBins();    Int_t   nclmax    = GetHistogramRanges()->GetHistoNClustersMax();    Int_t   nclmin    = GetHistogramRanges()->GetHistoNClustersMin();
1669   Int_t ncebins     = GetHistogramRanges()->GetHistoNCellsBins();       Int_t   ncemax    = GetHistogramRanges()->GetHistoNCellsMax();       Int_t   ncemin    = GetHistogramRanges()->GetHistoNCellsMin();
1670   Int_t nceclbins   = GetHistogramRanges()->GetHistoNClusterCellBins(); Int_t   nceclmax  = GetHistogramRanges()->GetHistoNClusterCellMax(); Int_t   nceclmin  = GetHistogramRanges()->GetHistoNClusterCellMin();
1671   Int_t nvdistbins  = GetHistogramRanges()->GetHistoVertexDistBins();   Float_t vdistmax  = GetHistogramRanges()->GetHistoVertexDistMax();   Float_t vdistmin  = GetHistogramRanges()->GetHistoVertexDistMin();
1672   Int_t rbins       = GetHistogramRanges()->GetHistoRBins();            Float_t rmax      = GetHistogramRanges()->GetHistoRMax();            Float_t rmin      = GetHistogramRanges()->GetHistoRMin();
1673   Int_t xbins       = GetHistogramRanges()->GetHistoXBins();            Float_t xmax      = GetHistogramRanges()->GetHistoXMax();            Float_t xmin      = GetHistogramRanges()->GetHistoXMin();
1674   Int_t ybins       = GetHistogramRanges()->GetHistoYBins();            Float_t ymax      = GetHistogramRanges()->GetHistoYMax();            Float_t ymin      = GetHistogramRanges()->GetHistoYMin();
1675   Int_t zbins       = GetHistogramRanges()->GetHistoZBins();            Float_t zmax      = GetHistogramRanges()->GetHistoZMax();            Float_t zmin      = GetHistogramRanges()->GetHistoZMin();
1676   Int_t ssbins      = GetHistogramRanges()->GetHistoShowerShapeBins();  Float_t ssmax     = GetHistogramRanges()->GetHistoShowerShapeMax();  Float_t ssmin     = GetHistogramRanges()->GetHistoShowerShapeMin();
1677   Int_t tdbins      = GetHistogramRanges()->GetHistoDiffTimeBins() ;    Float_t tdmax     = GetHistogramRanges()->GetHistoDiffTimeMax();     Float_t tdmin     = GetHistogramRanges()->GetHistoDiffTimeMin();
1678
1679   Int_t nv0sbins    = GetHistogramRanges()->GetHistoV0SignalBins();          Int_t nv0smax = GetHistogramRanges()->GetHistoV0SignalMax();          Int_t nv0smin = GetHistogramRanges()->GetHistoV0SignalMin();
1680   Int_t nv0mbins    = GetHistogramRanges()->GetHistoV0MultiplicityBins();    Int_t nv0mmax = GetHistogramRanges()->GetHistoV0MultiplicityMax();    Int_t nv0mmin = GetHistogramRanges()->GetHistoV0MultiplicityMin();
1681   Int_t ntrmbins    = GetHistogramRanges()->GetHistoTrackMultiplicityBins(); Int_t ntrmmax = GetHistogramRanges()->GetHistoTrackMultiplicityMax(); Int_t ntrmmin = GetHistogramRanges()->GetHistoTrackMultiplicityMin();
1682
1683   //EMCAL
1684   fNMaxCols = 48;
1685   fNMaxRows = 24;
1686   fNRCU     = 2 ;
1687   //PHOS
1688   if(fCalorimeter=="PHOS")
1689   {
1690     fNMaxCols = 56;
1691     fNMaxRows = 64;
1692     fNRCU     = 4 ;
1693   }
1694
1695   fhE  = new TH1F ("hE","E reconstructed clusters ", nptbins*5,ptmin,ptmax*5);
1696   fhE->SetXTitle("E (GeV)");
1697   outputContainer->Add(fhE);
1698
1699   if(fFillAllTH12)
1700   {
1701     fhPt  = new TH1F ("hPt","p_{T} reconstructed clusters", nptbins,ptmin,ptmax);
1702     fhPt->SetXTitle("p_{T} (GeV/c)");
1703     outputContainer->Add(fhPt);
1704
1705     fhPhi  = new TH1F ("hPhi","#phi reconstructed clusters ",nphibins,phimin,phimax);
1706     fhPhi->SetXTitle("#phi (rad)");
1707     outputContainer->Add(fhPhi);
1708
1709     fhEta  = new TH1F ("hEta","#eta reconstructed clusters ",netabins,etamin,etamax);
1710     fhEta->SetXTitle("#eta ");
1711     outputContainer->Add(fhEta);
1712   }
1713
1714   if(fFillAllTH3)
1715   {
1716     fhEtaPhiE  = new TH3F ("hEtaPhiE","#eta vs #phi vs energy, reconstructed clusters",
1717                            netabins,etamin,etamax,nphibins,phimin,phimax,nptbins,ptmin,ptmax);
1718     fhEtaPhiE->SetXTitle("#eta ");
1719     fhEtaPhiE->SetYTitle("#phi (rad)");
1720     fhEtaPhiE->SetZTitle("E (GeV) ");
1721     outputContainer->Add(fhEtaPhiE);
1722   }
1723
1724   fhClusterTimeEnergy  = new TH2F ("hClusterTimeEnergy","energy vs TOF, reconstructed clusters",
1725                                    nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
1726   fhClusterTimeEnergy->SetXTitle("E (GeV) ");
1727   fhClusterTimeEnergy->SetYTitle("TOF (ns)");
1728   outputContainer->Add(fhClusterTimeEnergy);
1729
1730   fhClusterPairDiffTimeE = new TH2F("hClusterPairDiffTimeE","cluster pair time difference vs E, only good clusters",
1731                                     nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
1732   fhClusterPairDiffTimeE->SetXTitle("E_{cluster} (GeV)");
1733   fhClusterPairDiffTimeE->SetYTitle("#Delta t (ns)");
1734   outputContainer->Add(fhClusterPairDiffTimeE);
1735
1736   fhLambda0  = new TH2F ("hLambda0","shower shape, #lambda^{2}_{0} vs E",
1737                          nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1738   fhLambda0->SetXTitle("E_{cluster}");
1739   fhLambda0->SetYTitle("#lambda^{2}_{0}");
1740   outputContainer->Add(fhLambda0);
1741
1742   fhLambda1  = new TH2F ("hLambda1","shower shape, #lambda^{2}_{1} vs E for bad cluster ",
1743                          nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1744   fhLambda1->SetXTitle("E_{cluster}");
1745   fhLambda1->SetYTitle("#lambda^{2}_{1}");
1746   outputContainer->Add(fhLambda1);
1747
1748   fhDispersion  = new TH2F ("hDispersion","shower shape, Dispersion^{2} vs E for bad cluster ",
1749                             nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1750   fhDispersion->SetXTitle("E_{cluster}");
1751   fhDispersion->SetYTitle("Dispersion");
1752   outputContainer->Add(fhDispersion);
1753
1754   fhClusterMaxCellCloseCellRatio  = new TH2F ("hClusterMaxCellCloseCellRatio","energy vs ratio of max cell / neighbour cell, reconstructed clusters",
1755                                               nptbins,ptmin,ptmax, 100,0,1.);
1756   fhClusterMaxCellCloseCellRatio->SetXTitle("E_{cluster} (GeV) ");
1757   fhClusterMaxCellCloseCellRatio->SetYTitle("E_{cell i}/E_{cell max}");
1758   outputContainer->Add(fhClusterMaxCellCloseCellRatio);
1759
1760   fhClusterMaxCellCloseCellDiff  = new TH2F ("hClusterMaxCellCloseCellDiff","energy vs ratio of max cell / neighbour cell, reconstructed clusters",
1761                                              nptbins,ptmin,ptmax, 500,0,100.);
1762   fhClusterMaxCellCloseCellDiff->SetXTitle("E_{cluster} (GeV) ");
1763   fhClusterMaxCellCloseCellDiff->SetYTitle("E_{cell max}-E_{cell i} (GeV)");
1764   outputContainer->Add(fhClusterMaxCellCloseCellDiff);
1765
1766   fhClusterMaxCellDiff  = new TH2F ("hClusterMaxCellDiff","energy vs difference of cluster energy - max cell energy / cluster energy, good clusters",
1767                                     nptbins,ptmin,ptmax, 500,0,1.);
1768   fhClusterMaxCellDiff->SetXTitle("E_{cluster} (GeV) ");
1769   fhClusterMaxCellDiff->SetYTitle("(E_{cluster} - E_{cell max})/ E_{cluster}");
1770   outputContainer->Add(fhClusterMaxCellDiff);
1771
1772   fhClusterMaxCellDiffNoCut  = new TH2F ("hClusterMaxCellDiffNoCut","energy vs difference of cluster energy - max cell energy / cluster energy",
1773                                          nptbins,ptmin,ptmax, 500,0,1.);
1774   fhClusterMaxCellDiffNoCut->SetXTitle("E_{cluster} (GeV) ");
1775   fhClusterMaxCellDiffNoCut->SetYTitle("(E_{cluster} - E_{cell max})/ E_{cluster}");
1776   outputContainer->Add(fhClusterMaxCellDiffNoCut);
1777
1778   fhClusterMaxCellECross  = new TH2F ("hClusterMaxCellECross","1 - Energy in cross around max energy cell / max energy cell vs cluster energy, good clusters",
1779                                        nptbins,ptmin,ptmax, 400,-1,1.);
1780   fhClusterMaxCellECross->SetXTitle("E_{cluster} (GeV) ");
1781   fhClusterMaxCellECross->SetYTitle("1- E_{cross}/E_{cell max}");
1782   outputContainer->Add(fhClusterMaxCellECross);
1783
1784   fhNCellsPerClusterNoCut  = new TH2F ("hNCellsPerClusterNoCut","# cells per cluster vs energy, no bad clusters cut",
1785                                        nptbins,ptmin,ptmax, nceclbins,nceclmin,nceclmax);
1786   fhNCellsPerClusterNoCut->SetXTitle("E (GeV)");
1787   fhNCellsPerClusterNoCut->SetYTitle("n cells");
1788   outputContainer->Add(fhNCellsPerClusterNoCut);
1789
1790   fhNCellsPerCluster  = new TH2F ("hNCellsPerCluster","# cells per cluster vs energy",nptbins,ptmin,ptmax, nceclbins,nceclmin,nceclmax);
1791   fhNCellsPerCluster->SetXTitle("E (GeV)");
1792   fhNCellsPerCluster->SetYTitle("n cells");
1793   outputContainer->Add(fhNCellsPerCluster);
1794
1795   fhNClusters  = new TH1F ("hNClusters","# clusters", nclbins,nclmin,nclmax);
1796   fhNClusters->SetXTitle("number of clusters");
1797   outputContainer->Add(fhNClusters);
1798
1799   if(fStudyBadClusters)
1800   {
1801     fhBadClusterEnergy  = new TH1F ("hBadClusterEnergy","Bad cluster energy", nptbins,ptmin,ptmax);
1802     fhBadClusterEnergy->SetXTitle("E_{cluster} (GeV) ");
1803     outputContainer->Add(fhBadClusterEnergy);
1804
1805     fhBadClusterMaxCellCloseCellRatio  = new TH2F ("hBadClusterMaxCellCloseCellRatio","energy vs ratio of max cell / neighbour cell constributing cell, reconstructed bad clusters",
1806                                                    nptbins,ptmin,ptmax, 100,0,1.);
1807     fhBadClusterMaxCellCloseCellRatio->SetXTitle("E_{cluster} (GeV) ");
1808     fhBadClusterMaxCellCloseCellRatio->SetYTitle("ratio");
1809     outputContainer->Add(fhBadClusterMaxCellCloseCellRatio);
1810
1811     fhBadClusterMaxCellCloseCellDiff  = new TH2F ("hBadClusterMaxCellCloseCellDiff","energy vs ratio of max cell - neighbour cell constributing cell, reconstructed bad clusters",
1812                                                   nptbins,ptmin,ptmax, 500,0,100);
1813     fhBadClusterMaxCellCloseCellDiff->SetXTitle("E_{cluster} (GeV) ");
1814     fhBadClusterMaxCellCloseCellDiff->SetYTitle("E_{cell max} - E_{cell i} (GeV)");
1815     outputContainer->Add(fhBadClusterMaxCellCloseCellDiff);
1816
1817     fhBadClusterMaxCellDiff  = new TH2F ("hBadClusterMaxCellDiff","energy vs difference of cluster energy - max cell energy / cluster energy for bad clusters",
1818                                          nptbins,ptmin,ptmax, 500,0,1.);
1819     fhBadClusterMaxCellDiff->SetXTitle("E_{cluster} (GeV) ");
1820     fhBadClusterMaxCellDiff->SetYTitle("(E_{cluster} - E_{cell max}) / E_{cluster}");
1821     outputContainer->Add(fhBadClusterMaxCellDiff);
1822
1823     fhBadClusterTimeEnergy  = new TH2F ("hBadClusterTimeEnergy","energy vs TOF of reconstructed bad clusters",
1824                                         nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
1825     fhBadClusterTimeEnergy->SetXTitle("E_{cluster} (GeV) ");
1826     fhBadClusterTimeEnergy->SetYTitle("TOF (ns)");
1827     outputContainer->Add(fhBadClusterTimeEnergy);
1828
1829     fhBadClusterPairDiffTimeE = new TH2F("hBadClusterPairDiffTimeE","cluster pair time difference (bad - good) vs E from bad cluster",nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
1830     fhBadClusterPairDiffTimeE->SetXTitle("E_{bad cluster} (GeV)");
1831     fhBadClusterPairDiffTimeE->SetYTitle("#Delta t (ns)");
1832     outputContainer->Add(fhBadClusterPairDiffTimeE);
1833
1834     fhBadClusterMaxCellECross  = new TH2F ("hBadClusterMaxCellECross","1 - Energy in cross around max energy cell / max energy cell vs cluster energy, bad clusters",
1835                                         nptbins,ptmin,ptmax, 400,-1,1.);
1836     fhBadClusterMaxCellECross->SetXTitle("E_{cluster} (GeV) ");
1837     fhBadClusterMaxCellECross->SetYTitle("1- E_{cross}/E_{cell max}");
1838     outputContainer->Add(fhBadClusterMaxCellECross);
1839
1840     if(fFillAllCellTimeHisto)
1841     {
1842       fhBadCellTimeSpreadRespectToCellMax = new TH2F ("hBadCellTimeSpreadRespectToCellMax","t_{cell max}-t_{cell i} from bad cluster", nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
1843       fhBadCellTimeSpreadRespectToCellMax->SetXTitle("E (GeV)");
1844       fhBadCellTimeSpreadRespectToCellMax->SetYTitle("#Delta t_{cell max - i} (ns)");
1845       outputContainer->Add(fhBadCellTimeSpreadRespectToCellMax);
1846
1847       fhBadClusterMaxCellDiffAverageTime = new TH2F ("hBadClusterMaxCellDiffAverageTime","t_{cell max}-t_{average} from bad cluster", nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
1848       fhBadClusterMaxCellDiffAverageTime->SetXTitle("E (GeV)");
1849       fhBadClusterMaxCellDiffAverageTime->SetYTitle("#Delta t_{cell max - average} (ns)");
1850       outputContainer->Add(fhBadClusterMaxCellDiffAverageTime);
1851
1852       fhBadClusterMaxCellDiffWeightedTime = new TH2F ("hBadClusterMaxCellDiffWeightedTime","t_{cell max}-t_{weighted} from bad cluster", nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
1853       fhBadClusterMaxCellDiffWeightedTime->SetXTitle("E (GeV)");
1854       fhBadClusterMaxCellDiffWeightedTime->SetYTitle("#Delta t_{cell max - weighted} (ns)");
1855       outputContainer->Add(fhBadClusterMaxCellDiffWeightedTime);
1856
1857     }
1858
1859   }
1860
1861   if(fStudyExotic)
1862   {
1863     fhExoL0ECross  = new TH2F("hExoL0_ECross",
1864                                "#lambda^{2}_{0} vs 1-E_{+}/E_{max} for E > 5 GeV",
1865                                400,0,1,ssbins,ssmin,ssmax);
1866     fhExoL0ECross ->SetXTitle("1-E_{+}/E_{cell max}");
1867     fhExoL0ECross ->SetYTitle("#lambda^{2}_{0}");
1868     outputContainer->Add(fhExoL0ECross) ;
1869
1870     fhExoL1ECross  = new TH2F("hExoL1_ECross",
1871                               "#lambda^{2}_{1} vs 1-E_{+}/E_{max} for E > 5 GeV",
1872                               400,0,1,ssbins,ssmin,ssmax);
1873     fhExoL1ECross ->SetXTitle("1-E_{+}/E_{cell max}");
1874     fhExoL1ECross ->SetYTitle("#lambda^{2}_{1}");
1875     outputContainer->Add(fhExoL1ECross) ;
1876
1877     for(Int_t ie = 0; ie <fExoNECrossCuts; ie++)
1878     {
1879
1880       fhExoDTime[ie]  = new TH2F(Form("hExoDTime_ECross%d",ie),
1881                                  Form("#Delta time = t_{max}-t_{cells} vs E_{cluster} for exotic, 1-E_{+}/E_{max} < %2.2f",fExoECrossCuts[ie]),
1882                                  nptbins,ptmin,ptmax,tdbins,tdmin,tdmax);
1883       fhExoDTime[ie] ->SetYTitle("#Delta t (ns)");
1884       fhExoDTime[ie] ->SetXTitle("E (GeV)");
1885       outputContainer->Add(fhExoDTime[ie]) ;
1886
1887       for(Int_t idt = 0; idt < fExoNDTimeCuts; idt++)
1888       {
1889         fhExoNCell[ie][idt]  = new TH2F(Form("hExoNCell_ECross%d_DT%d",ie,idt),
1890                                      Form("N cells per cluster vs E cluster, 1-E_{+}/E_{max} < %2.2f, #Delta t < %2.0f",fExoECrossCuts[ie],fExoDTimeCuts[idt]),
1891                                      nptbins,ptmin,ptmax,nceclbins,nceclmin,nceclmax);
1892         fhExoNCell[ie][idt] ->SetYTitle("N cells");
1893         fhExoNCell[ie][idt] ->SetXTitle("E (GeV)");
1894         outputContainer->Add(fhExoNCell[ie][idt]) ;
1895
1896         fhExoL0   [ie][idt]  = new TH2F(Form("hExoL0_ECross%d_DT%d",ie,idt),
1897                                      Form("#lambda^{2}_{0} vs E cluster for exotic, 1-E_{+}/E_{max} < %2.2f, #Delta t = %2.0f",fExoECrossCuts[ie],fExoDTimeCuts[idt]),
1898                                      nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1899         fhExoL0   [ie][idt] ->SetYTitle("#lambda^{2}_{0}");
1900         fhExoL0   [ie][idt] ->SetXTitle("E (GeV)");
1901         outputContainer->Add(fhExoL0[ie][idt]) ;
1902
1903         fhExoL1   [ie][idt]  = new TH2F(Form("hExoL1_ECross%d_DT%d",ie,idt),
1904                                         Form("#lambda^{2}_{1} vs E cluster for exotic, 1-E_{+}/E_{max} < %2.2f, #Delta t = %2.0f",fExoECrossCuts[ie],fExoDTimeCuts[idt]),
1905                                         nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1906         fhExoL1   [ie][idt] ->SetYTitle("#lambda^{2}_{1}");
1907         fhExoL1   [ie][idt] ->SetXTitle("E (GeV)");
1908         outputContainer->Add(fhExoL1[ie][idt]) ;
1909
1910         fhExoECross[ie][idt]  = new TH2F(Form("hExoECross_ECross%d_DT%d",ie,idt),
1911                                       Form("E cross for cells vs E cell, 1-E_{+}/E_{max} < %2.2f, #Delta t < %2.0f",fExoECrossCuts[ie],fExoDTimeCuts[idt]),
1912                                       nptbins,ptmin,ptmax,400,0,1);
1913         fhExoECross[ie][idt] ->SetYTitle("1-E_{+}/E_{cell max}");
1914         fhExoECross[ie][idt] ->SetXTitle("E_{cell} (GeV)");
1915         outputContainer->Add(fhExoECross[ie][idt]) ;
1916
1917         fhExoTime  [ie][idt]  = new TH2F(Form("hExoTime_ECross%d_DT%d",ie,idt),
1918                                         Form("Time of cluster (max cell) vs E cluster for exotic, 1-E_{+}/E_{max} < %2.2f, #Delta t = %2.0f",fExoECrossCuts[ie],fExoDTimeCuts[idt]),
1919                                         nptbins,ptmin,ptmax,ntimebins,timemin,timemax);
1920         fhExoTime  [ie][idt] ->SetYTitle("time_{max} (ns)");
1921         fhExoTime  [ie][idt] ->SetXTitle("E (GeV)");
1922         outputContainer->Add(fhExoTime[ie][idt]) ;
1923
1924         fhExoL0NCell[ie][idt]  = new TH2F(Form("hExoL0_NCell%d_DT%d",ie,idt),
1925                                           Form("#lambda^{2}_{0} vs N cells per clusters for E > 5 GeV, 1-E_{+}/E_{max} < %2.2f, #Delta t = %2.0f",fExoECrossCuts[ie],fExoDTimeCuts[idt]),
1926                                           nptbins,ptmin,ptmax,ntimebins,timemin,timemax);
1927         fhExoL0NCell[ie][idt] ->SetYTitle("N cells");
1928         fhExoL0NCell[ie][idt] ->SetXTitle("#lambda^{2}_{0}");
1929         outputContainer->Add(fhExoL0NCell[ie][idt]) ;
1930
1931         fhExoL1NCell[ie][idt]  = new TH2F(Form("hExoL1_NCell%d_DT%d",ie,idt),
1932                                           Form("#lambda^{2}_{1} vs N cells per clusters for E > 5 GeV, 1-E_{+}/E_{max} < %2.2f, #Delta t = %2.0f",fExoECrossCuts[ie],fExoDTimeCuts[idt]),
1933                                           nptbins,ptmin,ptmax,ntimebins,timemin,timemax);
1934         fhExoL1NCell[ie][idt] ->SetYTitle("N cells");
1935         fhExoL1NCell[ie][idt] ->SetXTitle("#lambda^{2}_{1}");
1936         outputContainer->Add(fhExoL1NCell[ie][idt]) ;
1937
1938       }
1939     }
1940   }
1941
1942   // Cluster size in terms of cells
1943   if(fStudyClustersAsymmetry)
1944   {
1945     fhDeltaIEtaDeltaIPhiE0[0]  = new TH2F ("hDeltaIEtaDeltaIPhiE0"," Cluster size in columns vs rows for E < 2 GeV, n cells > 3",
1946                                            50,0,50,50,0,50);
1947     fhDeltaIEtaDeltaIPhiE0[0]->SetXTitle("#Delta Column");
1948     fhDeltaIEtaDeltaIPhiE0[0]->SetYTitle("#Delta Row");
1949     outputContainer->Add(fhDeltaIEtaDeltaIPhiE0[0]);
1950
1951     fhDeltaIEtaDeltaIPhiE2[0]  = new TH2F ("hDeltaIEtaDeltaIPhiE2"," Cluster size in columns vs rows for 2 <E < 6 GeV, n cells > 3",
1952                                            50,0,50,50,0,50);
1953     fhDeltaIEtaDeltaIPhiE2[0]->SetXTitle("#Delta Column");
1954     fhDeltaIEtaDeltaIPhiE2[0]->SetYTitle("#Delta Row");
1955     outputContainer->Add(fhDeltaIEtaDeltaIPhiE2[0]);
1956
1957     fhDeltaIEtaDeltaIPhiE6[0]  = new TH2F ("hDeltaIEtaDeltaIPhiE6"," Cluster size in columns vs rows for E > 6 GeV, n cells > 3",
1958                                            50,0,50,50,0,50);
1959     fhDeltaIEtaDeltaIPhiE6[0]->SetXTitle("#Delta Column");
1960     fhDeltaIEtaDeltaIPhiE6[0]->SetYTitle("#Delta Row");
1961     outputContainer->Add(fhDeltaIEtaDeltaIPhiE6[0]);
1962
1963     fhDeltaIA[0]  = new TH2F ("hDeltaIA"," Cluster *asymmetry* in cell units vs E",
1964                               nptbins,ptmin,ptmax,21,-1.05,1.05);
1965     fhDeltaIA[0]->SetXTitle("E_{cluster}");
1966     fhDeltaIA[0]->SetYTitle("A_{cell in cluster}");
1967     outputContainer->Add(fhDeltaIA[0]);
1968
1969     fhDeltaIAL0[0]  = new TH2F ("hDeltaIAL0"," Cluster *asymmetry* in cell units vs #lambda^{2}_{0}",
1970                                 ssbins,ssmin,ssmax,21,-1.05,1.05);
1971     fhDeltaIAL0[0]->SetXTitle("#lambda^{2}_{0}");
1972     fhDeltaIAL0[0]->SetYTitle("A_{cell in cluster}");
1973     outputContainer->Add(fhDeltaIAL0[0]);
1974
1975     fhDeltaIAL1[0]  = new TH2F ("hDeltaIAL1"," Cluster *asymmetry* in cell units vs #lambda^{2}_{1}",
1976                                 ssbins,ssmin,ssmax,21,-1.05,1.05);
1977     fhDeltaIAL1[0]->SetXTitle("#lambda^{2}_{1}");
1978     fhDeltaIAL1[0]->SetYTitle("A_{cell in cluster}");
1979     outputContainer->Add(fhDeltaIAL1[0]);
1980
1981     fhDeltaIANCells[0]  = new TH2F ("hDeltaIANCells"," Cluster *asymmetry* in cell units vs N cells in cluster",
1982                                     nceclbins,nceclmin,nceclmax,21,-1.05,1.05);
1983     fhDeltaIANCells[0]->SetXTitle("N_{cell in cluster}");
1984     fhDeltaIANCells[0]->SetYTitle("A_{cell in cluster}");
1985     outputContainer->Add(fhDeltaIANCells[0]);
1986
1987
1988     fhDeltaIEtaDeltaIPhiE0[1]  = new TH2F ("hDeltaIEtaDeltaIPhiE0Charged"," Cluster size in columns vs rows for E < 2 GeV, n cells > 3, matched with track",
1989                                            50,0,50,50,0,50);
1990     fhDeltaIEtaDeltaIPhiE0[1]->SetXTitle("#Delta Column");
1991     fhDeltaIEtaDeltaIPhiE0[1]->SetYTitle("#Delta Row");
1992     outputContainer->Add(fhDeltaIEtaDeltaIPhiE0[1]);
1993
1994     fhDeltaIEtaDeltaIPhiE2[1]  = new TH2F ("hDeltaIEtaDeltaIPhiE2Charged"," Cluster size in columns vs rows for 2 <E < 6 GeV, n cells > 3, matched with track",
1995                                            50,0,50,50,0,50);
1996     fhDeltaIEtaDeltaIPhiE2[1]->SetXTitle("#Delta Column");
1997     fhDeltaIEtaDeltaIPhiE2[1]->SetYTitle("#Delta Row");
1998     outputContainer->Add(fhDeltaIEtaDeltaIPhiE2[1]);
1999
2000     fhDeltaIEtaDeltaIPhiE6[1]  = new TH2F ("hDeltaIEtaDeltaIPhiE6Charged"," Cluster size in columns vs rows for E > 6 GeV, n cells > 3, matched with track",
2001                                            50,0,50,50,0,50);
2002     fhDeltaIEtaDeltaIPhiE6[1]->SetXTitle("#Delta Column");
2003     fhDeltaIEtaDeltaIPhiE6[1]->SetYTitle("#Delta Row");
2004     outputContainer->Add(fhDeltaIEtaDeltaIPhiE6[1]);
2005
2006     fhDeltaIA[1]  = new TH2F ("hDeltaIACharged"," Cluster *asymmetry* in cell units vs E, matched with track",
2007                               nptbins,ptmin,ptmax,21,-1.05,1.05);
2008     fhDeltaIA[1]->SetXTitle("E_{cluster}");
2009     fhDeltaIA[1]->SetYTitle("A_{cell in cluster}");
2010     outputContainer->Add(fhDeltaIA[1]);
2011
2012     fhDeltaIAL0[1]  = new TH2F ("hDeltaIAL0Charged"," Cluster *asymmetry* in cell units vs #lambda^{2}_{0}, matched with track",
2013                                 ssbins,ssmin,ssmax,21,-1.05,1.05);
2014     fhDeltaIAL0[1]->SetXTitle("#lambda^{2}_{0}");
2015     fhDeltaIAL0[1]->SetYTitle("A_{cell in cluster}");
2016     outputContainer->Add(fhDeltaIAL0[1]);
2017
2018     fhDeltaIAL1[1]  = new TH2F ("hDeltaIAL1Charged"," Cluster *asymmetry* in cell units vs #lambda^{2}_{1}, matched with track",
2019                                 ssbins,ssmin,ssmax,21,-1.05,1.05);
2020     fhDeltaIAL1[1]->SetXTitle("#lambda^{2}_{1}");
2021     fhDeltaIAL1[1]->SetYTitle("A_{cell in cluster}");
2022     outputContainer->Add(fhDeltaIAL1[1]);
2023
2024     fhDeltaIANCells[1]  = new TH2F ("hDeltaIANCellsCharged"," Cluster *asymmetry* in cell units vs N cells in cluster, matched with track",
2025                                     nceclbins,nceclmin,nceclmax,21,-1.05,1.05);
2026     fhDeltaIANCells[1]->SetXTitle("N_{cell in cluster}");
2027     fhDeltaIANCells[1]->SetYTitle("A_{cell in cluster}");
2028     outputContainer->Add(fhDeltaIANCells[1]);
2029
2030     if(IsDataMC()){
2031       TString particle[]={"Photon","Electron","Conversion","Hadron"};
2032       for (Int_t iPart = 0; iPart < 4; iPart++) {
2033
2034         fhDeltaIAMC[iPart]  = new TH2F (Form("hDeltaIA_MC%s",particle[iPart].Data()),Form(" Cluster *asymmetry* in cell units vs E, from %s",particle[iPart].Data()),
2035                                         nptbins,ptmin,ptmax,21,-1.05,1.05);
2036         fhDeltaIAMC[iPart]->SetXTitle("E_{cluster}");
2037         fhDeltaIAMC[iPart]->SetYTitle("A_{cell in cluster}");
2038         outputContainer->Add(fhDeltaIAMC[iPart]);
2039       }
2040     }
2041
2042     if(fStudyBadClusters)
2043     {
2044       fhBadClusterDeltaIEtaDeltaIPhiE0  = new TH2F ("hBadClusterDeltaIEtaDeltaIPhiE0"," Cluster size in columns vs rows for E < 2 GeV, n cells > 3",
2045                                                     50,0,50,50,0,50);
2046       fhBadClusterDeltaIEtaDeltaIPhiE0->SetXTitle("#Delta Column");
2047       fhBadClusterDeltaIEtaDeltaIPhiE0->SetYTitle("#Delta Row");
2048       outputContainer->Add(fhBadClusterDeltaIEtaDeltaIPhiE0);
2049
2050       fhBadClusterDeltaIEtaDeltaIPhiE2  = new TH2F ("hBadClusterDeltaIEtaDeltaIPhiE2"," Cluster size in columns vs rows for 2 <E < 6 GeV, n cells > 3",
2051                                                     50,0,50,50,0,50);
2052       fhBadClusterDeltaIEtaDeltaIPhiE2->SetXTitle("#Delta Column");
2053       fhBadClusterDeltaIEtaDeltaIPhiE2->SetYTitle("#Delta Row");
2054       outputContainer->Add(fhBadClusterDeltaIEtaDeltaIPhiE2);
2055
2056       fhBadClusterDeltaIEtaDeltaIPhiE6  = new TH2F ("hBadClusterDeltaIEtaDeltaIPhiE6"," Cluster size in columns vs rows for E > 6 GeV, n cells > 3",
2057                                                     50,0,50,50,0,50);
2058       fhBadClusterDeltaIEtaDeltaIPhiE6->SetXTitle("#Delta Column");
2059       fhBadClusterDeltaIEtaDeltaIPhiE6->SetYTitle("#Delta Row");
2060       outputContainer->Add(fhBadClusterDeltaIEtaDeltaIPhiE6);
2061
2062       fhBadClusterDeltaIA  = new TH2F ("hBadClusterDeltaIA"," Cluster *asymmetry* in cell units vs E",
2063                                        nptbins,ptmin,ptmax,21,-1.05,1.05);
2064       fhBadClusterDeltaIA->SetXTitle("E_{cluster}");
2065       fhBadClusterDeltaIA->SetYTitle("A_{cell in cluster}");
2066       outputContainer->Add(fhBadClusterDeltaIA);
2067     }
2068   }
2069
2070   if(fStudyWeight)
2071   {
2072     fhECellClusterRatio  = new TH2F ("hECellClusterRatio"," cell energy / cluster energy vs cluster energy",
2073                                      nptbins,ptmin,ptmax, 100,0,1.);
2074     fhECellClusterRatio->SetXTitle("E_{cluster} (GeV) ");
2075     fhECellClusterRatio->SetYTitle("E_{cell i}/E_{cluster}");
2076     outputContainer->Add(fhECellClusterRatio);
2077
2078     fhECellClusterLogRatio  = new TH2F ("hECellClusterLogRatio"," Log(cell energy / cluster energy) vs cluster energy",
2079                                         nptbins,ptmin,ptmax, 100,-10,0);
2080     fhECellClusterLogRatio->SetXTitle("E_{cluster} (GeV) ");
2081     fhECellClusterLogRatio->SetYTitle("Log(E_{cell i}/E_{cluster})");
2082     outputContainer->Add(fhECellClusterLogRatio);
2083
2084     fhEMaxCellClusterRatio  = new TH2F ("hEMaxCellClusterRatio"," max cell energy / cluster energy vs cluster energy",
2085                                         nptbins,ptmin,ptmax, 100,0,1.);
2086     fhEMaxCellClusterRatio->SetXTitle("E_{cluster} (GeV) ");
2087     fhEMaxCellClusterRatio->SetYTitle("E_{max cell}/E_{cluster}");
2088     outputContainer->Add(fhEMaxCellClusterRatio);
2089
2090     fhEMaxCellClusterLogRatio  = new TH2F ("hEMaxCellClusterLogRatio"," Log(max cell energy / cluster energy) vs cluster energy",
2091                                            nptbins,ptmin,ptmax, 100,-10,0);
2092     fhEMaxCellClusterLogRatio->SetXTitle("E_{cluster} (GeV) ");
2093     fhEMaxCellClusterLogRatio->SetYTitle("Log (E_{max cell}/E_{cluster})");
2094     outputContainer->Add(fhEMaxCellClusterLogRatio);
2095
2096     for(Int_t iw = 0; iw < 14; iw++){
2097       fhLambda0ForW0[iw]  = new TH2F (Form("hLambda0ForW0%d",iw),Form("shower shape, #lambda^{2}_{0} vs E, w0 = %1.1f",1+0.5*iw),
2098                                       nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2099       fhLambda0ForW0[iw]->SetXTitle("E_{cluster}");
2100       fhLambda0ForW0[iw]->SetYTitle("#lambda^{2}_{0}");
2101       outputContainer->Add(fhLambda0ForW0[iw]);
2102
2103 //      fhLambda1ForW0[iw]  = new TH2F (Form("hLambda1ForW0%d",iw),Form("shower shape, #lambda^{2}_{1} vs E, w0 = %1.1f",1+0.5*iw),
2104 //                                      nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2105 //      fhLambda1ForW0[iw]->SetXTitle("E_{cluster}");
2106 //      fhLambda1ForW0[iw]->SetYTitle("#lambda^{2}_{1}");
2107 //      outputContainer->Add(fhLambda1ForW0[iw]);
2108
2109       if(IsDataMC()){
2110         TString mcnames[] = {"Photon", "Electron","Conversion","Pi0","Hadron"};
2111         for(Int_t imc = 0; imc < 5; imc++){
2112           fhLambda0ForW0MC[iw][imc]  = new TH2F (Form("hLambda0ForW0%d_MC%s",iw,mcnames[imc].Data()),
2113                                                  Form("shower shape, #lambda^{2}_{0} vs E, w0 = %1.1f, for MC %s",1+0.5*iw,mcnames[imc].Data()),
2114                                                  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2115           fhLambda0ForW0MC[iw][imc]->SetXTitle("E_{cluster}");
2116           fhLambda0ForW0MC[iw][imc]->SetYTitle("#lambda^{2}_{0}");
2117           outputContainer->Add(fhLambda0ForW0MC[iw][imc]);
2118
2119 //          fhLambda1ForW0MC[iw][imc]  = new TH2F (Form("hLambda1ForW0%d_MC%s",iw,mcnames[imc].Data()),
2120 //                                                 Form("shower shape, #lambda^{2}_{1} vs E, w0 = %1.1f, for MC %s",1+0.5*iw,mcnames[imc].Data()),
2121 //                                                 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2122 //          fhLambda1ForW0MC[iw][imc]->SetXTitle("E_{cluster}");
2123 //          fhLambda1ForW0MC[iw][imc]->SetYTitle("#lambda^{2}_{1}");
2124 //          outputContainer->Add(fhLambda1ForW0MC[iw][imc]);
2125         }
2126       }
2127     }
2128   }
2129
2130   //Track Matching
2131   if(fFillAllTMHisto)
2132   {
2133     Int_t   nresetabins = GetHistogramRanges()->GetHistoTrackResidualEtaBins();
2134     Float_t resetamax   = GetHistogramRanges()->GetHistoTrackResidualEtaMax();
2135     Float_t resetamin   = GetHistogramRanges()->GetHistoTrackResidualEtaMin();
2136     Int_t   nresphibins = GetHistogramRanges()->GetHistoTrackResidualPhiBins();
2137     Float_t resphimax   = GetHistogramRanges()->GetHistoTrackResidualPhiMax();
2138     Float_t resphimin   = GetHistogramRanges()->GetHistoTrackResidualPhiMin();
2139
2140     fhTrackMatchedDEta  = new TH2F("hTrackMatchedDEta","d#eta of cluster-track vs cluster energy",
2141                                    nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
2142     fhTrackMatchedDEta->SetYTitle("d#eta");
2143     fhTrackMatchedDEta->SetXTitle("E_{cluster} (GeV)");
2144
2145     fhTrackMatchedDPhi  = new TH2F("hTrackMatchedDPhi","d#phi of cluster-track vs cluster energy",
2146                                    nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
2147     fhTrackMatchedDPhi->SetYTitle("d#phi (rad)");
2148     fhTrackMatchedDPhi->SetXTitle("E_{cluster} (GeV)");
2149
2150     fhTrackMatchedDEtaDPhi  = new TH2F("hTrackMatchedDEtaDPhi","d#eta vs d#phi of cluster-track vs cluster energy",
2151                                        nresetabins,resetamin,resetamax,nresphibins,resphimin,resphimax);
2152     fhTrackMatchedDEtaDPhi->SetYTitle("d#phi (rad)");
2153     fhTrackMatchedDEtaDPhi->SetXTitle("d#eta");
2154
2155     fhTrackMatchedDEtaPos  = new TH2F("hTrackMatchedDEtaPos","d#eta of cluster-track vs cluster energy",
2156                                       nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
2157     fhTrackMatchedDEtaPos->SetYTitle("d#eta");
2158     fhTrackMatchedDEtaPos->SetXTitle("E_{cluster} (GeV)");
2159
2160     fhTrackMatchedDPhiPos  = new TH2F("hTrackMatchedDPhiPos","d#phi of cluster-track vs cluster energy",
2161                                       nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
2162     fhTrackMatchedDPhiPos->SetYTitle("d#phi (rad)");
2163     fhTrackMatchedDPhiPos->SetXTitle("E_{cluster} (GeV)");
2164
2165     fhTrackMatchedDEtaDPhiPos  = new TH2F("hTrackMatchedDEtaDPhiPos","d#eta vs d#phi of cluster-track vs cluster energy",
2166                                           nresetabins,resetamin,resetamax,nresphibins,resphimin,resphimax);
2167     fhTrackMatchedDEtaDPhiPos->SetYTitle("d#phi (rad)");
2168     fhTrackMatchedDEtaDPhiPos->SetXTitle("d#eta");
2169
2170     outputContainer->Add(fhTrackMatchedDEta) ;
2171     outputContainer->Add(fhTrackMatchedDPhi) ;
2172     outputContainer->Add(fhTrackMatchedDEtaDPhi) ;
2173     outputContainer->Add(fhTrackMatchedDEtaPos) ;
2174     outputContainer->Add(fhTrackMatchedDPhiPos) ;
2175     outputContainer->Add(fhTrackMatchedDEtaDPhiPos) ;
2176
2177     if(fFillAllTH12)
2178     {
2179       fhECharged  = new TH1F ("hECharged","E reconstructed clusters, matched with track", nptbins,ptmin,ptmax);
2180       fhECharged->SetXTitle("E (GeV)");
2181       outputContainer->Add(fhECharged);
2182
2183       fhPtCharged  = new TH1F ("hPtCharged","p_{T} reconstructed clusters, matched with track", nptbins,ptmin,ptmax);
2184       fhPtCharged->SetXTitle("p_{T} (GeV/c)");
2185       outputContainer->Add(fhPtCharged);
2186
2187       fhPhiCharged  = new TH1F ("hPhiCharged","#phi reconstructed clusters, matched with track",nphibins,phimin,phimax);
2188       fhPhiCharged->SetXTitle("#phi (rad)");
2189       outputContainer->Add(fhPhiCharged);
2190
2191       fhEtaCharged  = new TH1F ("hEtaCharged","#eta reconstructed clusters, matched with track",netabins,etamin,etamax);
2192       fhEtaCharged->SetXTitle("#eta ");
2193       outputContainer->Add(fhEtaCharged);
2194     }
2195     if(fFillAllTH3){
2196       fhEtaPhiECharged  = new TH3F ("hEtaPhiECharged","#eta vs #phi, reconstructed clusters, matched with track",
2197                                     netabins,etamin,etamax,nphibins,phimin,phimax,nptbins,ptmin,ptmax);
2198       fhEtaPhiECharged->SetXTitle("#eta ");
2199       fhEtaPhiECharged->SetYTitle("#phi ");
2200       fhEtaPhiECharged->SetZTitle("E (GeV) ");
2201       outputContainer->Add(fhEtaPhiECharged);
2202     }
2203
2204     fh1EOverP = new TH2F("h1EOverP","TRACK matches E/p",nptbins,ptmin,ptmax, nPoverEbins,eOverPmin,eOverPmax);
2205     fh1EOverP->SetYTitle("E/p");
2206     fh1EOverP->SetXTitle("p_{T} (GeV/c)");
2207     outputContainer->Add(fh1EOverP);
2208
2209     fh2dR = new TH2F("h2dR","TRACK matches dR",nptbins,ptmin,ptmax,ndRbins,dRmin,dRmax);
2210     fh2dR->SetXTitle("#Delta R (rad)");
2211     fh2dR->SetXTitle("E cluster (GeV)");
2212     outputContainer->Add(fh2dR) ;
2213
2214     fh2MatchdEdx = new TH2F("h2MatchdEdx","dE/dx vs. p for all matches",nptbins,ptmin,ptmax,ndedxbins,dedxmin,dedxmax);
2215     fh2MatchdEdx->SetXTitle("p (GeV/c)");
2216     fh2MatchdEdx->SetYTitle("<dE/dx>");
2217     outputContainer->Add(fh2MatchdEdx);
2218
2219     fh2EledEdx = new TH2F("h2EledEdx","dE/dx vs. p for electrons",nptbins,ptmin,ptmax,ndedxbins,dedxmin,dedxmax);
2220     fh2EledEdx->SetXTitle("p (GeV/c)");
2221     fh2EledEdx->SetYTitle("<dE/dx>");
2222     outputContainer->Add(fh2EledEdx) ;
2223
2224     fh1EOverPR02 = new TH2F("h1EOverPR02","TRACK matches E/p, all",nptbins,ptmin,ptmax, nPoverEbins,eOverPmin,eOverPmax);
2225     fh1EOverPR02->SetYTitle("E/p");
2226     fh1EOverPR02->SetXTitle("p_{T} (GeV/c)");
2227     outputContainer->Add(fh1EOverPR02);
2228
2229     fh1EleEOverP = new TH2F("h1EleEOverP","Electron candidates E/p (60<dEdx<100)",nptbins,ptmin,ptmax, nPoverEbins,eOverPmin,eOverPmax);
2230     fh1EleEOverP->SetYTitle("E/p");
2231     fh1EleEOverP->SetXTitle("p_{T} (GeV/c)");
2232     outputContainer->Add(fh1EleEOverP);
2233   }
2234
2235   if(fFillAllPi0Histo)
2236   {
2237     fhIM  = new TH2F ("hIM","Cluster pairs Invariant mass vs reconstructed pair energy, ncell > 1",nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
2238     fhIM->SetXTitle("p_{T, cluster pairs} (GeV) ");
2239     fhIM->SetYTitle("M_{cluster pairs} (GeV/c^{2})");
2240     outputContainer->Add(fhIM);
2241
2242     fhAsym  = new TH2F ("hAssym","Cluster pairs Asymmetry vs reconstructed pair energy",nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax);
2243     fhAsym->SetXTitle("p_{T, cluster pairs} (GeV) ");
2244     fhAsym->SetYTitle("Asymmetry");
2245     outputContainer->Add(fhAsym);
2246   }
2247
2248
2249   if(fFillAllPosHisto2)
2250   {
2251     if(fFillAllTH3)
2252     {
2253       fhXYZ  = new TH3F ("hXYZ","Cluster: x vs y vs z",xbins,xmin,xmax,ybins,ymin,ymax,zbins,zmin,zmax);
2254       fhXYZ->SetXTitle("x (cm)");
2255       fhXYZ->SetYTitle("y (cm)");
2256       fhXYZ->SetZTitle("z (cm) ");
2257       outputContainer->Add(fhXYZ);
2258     }
2259
2260     fhXNCells  = new TH2F ("hXNCells","Cluster X position vs N Cells per Cluster",xbins,xmin,xmax,nceclbins,nceclmin,nceclmax);
2261     fhXNCells->SetXTitle("x (cm)");
2262     fhXNCells->SetYTitle("N cells per cluster");
2263     outputContainer->Add(fhXNCells);
2264
2265     fhZNCells  = new TH2F ("hZNCells","Cluster Z position vs N Cells per Cluster",zbins,zmin,zmax,nceclbins,nceclmin,nceclmax);
2266     fhZNCells->SetXTitle("z (cm)");
2267     fhZNCells->SetYTitle("N cells per cluster");
2268     outputContainer->Add(fhZNCells);
2269
2270     fhXE  = new TH2F ("hXE","Cluster X position vs cluster energy",xbins,xmin,xmax,nptbins,ptmin,ptmax);
2271     fhXE->SetXTitle("x (cm)");
2272     fhXE->SetYTitle("E (GeV)");
2273     outputContainer->Add(fhXE);
2274
2275     fhZE  = new TH2F ("hZE","Cluster Z position vs cluster energy",zbins,zmin,zmax,nptbins,ptmin,ptmax);
2276     fhZE->SetXTitle("z (cm)");
2277     fhZE->SetYTitle("E (GeV)");
2278     outputContainer->Add(fhZE);
2279
2280     fhRNCells  = new TH2F ("hRNCells","Cluster R position vs N Cells per Cluster",rbins,rmin,rmax,nceclbins,nceclmin,nceclmax);
2281     fhRNCells->SetXTitle("r = #sqrt{x^{2}+y^{2}} (cm)");
2282     fhRNCells->SetYTitle("N cells per cluster");
2283     outputContainer->Add(fhRNCells);
2284
2285
2286     fhYNCells  = new TH2F ("hYNCells","Cluster Y position vs N Cells per Cluster",ybins,ymin,ymax,nceclbins,nceclmin,nceclmax);
2287     fhYNCells->SetXTitle("y (cm)");
2288     fhYNCells->SetYTitle("N cells per cluster");
2289     outputContainer->Add(fhYNCells);
2290
2291     fhRE  = new TH2F ("hRE","Cluster R position vs cluster energy",rbins,rmin,rmax,nptbins,ptmin,ptmax);
2292     fhRE->SetXTitle("r = #sqrt{x^{2}+y^{2}} (cm)");
2293     fhRE->SetYTitle("E (GeV)");
2294     outputContainer->Add(fhRE);
2295
2296     fhYE  = new TH2F ("hYE","Cluster Y position vs cluster energy",ybins,ymin,ymax,nptbins,ptmin,ptmax);
2297     fhYE->SetXTitle("y (cm)");
2298     fhYE->SetYTitle("E (GeV)");
2299     outputContainer->Add(fhYE);
2300   }
2301
2302   if(fFillAllPosHisto)
2303   {
2304     fhRCellE  = new TH2F ("hRCellE","Cell R position vs cell energy",rbins,rmin,rmax,nptbins,ptmin,ptmax);
2305     fhRCellE->SetXTitle("r = #sqrt{x^{2}+y^{2}} (cm)");
2306     fhRCellE->SetYTitle("E (GeV)");
2307     outputContainer->Add(fhRCellE);
2308
2309     fhXCellE  = new TH2F ("hXCellE","Cell X position vs cell energy",xbins,xmin,xmax,nptbins,ptmin,ptmax);
2310     fhXCellE->SetXTitle("x (cm)");
2311     fhXCellE->SetYTitle("E (GeV)");
2312     outputContainer->Add(fhXCellE);
2313
2314     fhYCellE  = new TH2F ("hYCellE","Cell Y position vs cell energy",ybins,ymin,ymax,nptbins,ptmin,ptmax);
2315     fhYCellE->SetXTitle("y (cm)");
2316     fhYCellE->SetYTitle("E (GeV)");
2317     outputContainer->Add(fhYCellE);
2318
2319     fhZCellE  = new TH2F ("hZCellE","Cell Z position vs cell energy",zbins,zmin,zmax,nptbins,ptmin,ptmax);
2320     fhZCellE->SetXTitle("z (cm)");
2321     fhZCellE->SetYTitle("E (GeV)");
2322     outputContainer->Add(fhZCellE);
2323
2324     fhXYZCell  = new TH3F ("hXYZCell","Cell : x vs y vs z",xbins,xmin,xmax,ybins,ymin,ymax,zbins,zmin,zmax);
2325     fhXYZCell->SetXTitle("x (cm)");
2326     fhXYZCell->SetYTitle("y (cm)");
2327     fhXYZCell->SetZTitle("z (cm)");
2328     outputContainer->Add(fhXYZCell);
2329
2330
2331     Float_t dx = TMath::Abs(xmin)+TMath::Abs(xmax);
2332     Float_t dy = TMath::Abs(ymin)+TMath::Abs(ymax);
2333     Float_t dz = TMath::Abs(zmin)+TMath::Abs(zmax);
2334     Float_t dr = TMath::Abs(rmin)+TMath::Abs(rmax);
2335
2336     fhDeltaCellClusterRNCells  = new TH2F ("hDeltaCellClusterRNCells","Cluster-Cell R position vs N Cells per Cluster",rbins*2,-dr,dr,nceclbins,nceclmin,nceclmax);
2337     fhDeltaCellClusterRNCells->SetXTitle("r = #sqrt{x^{2}+y^{2}} (cm)");
2338     fhDeltaCellClusterRNCells->SetYTitle("N cells per cluster");
2339     outputContainer->Add(fhDeltaCellClusterRNCells);
2340
2341     fhDeltaCellClusterXNCells  = new TH2F ("hDeltaCellClusterXNCells","Cluster-Cell X position vs N Cells per Cluster",xbins*2,-dx,dx,nceclbins,nceclmin,nceclmax);
2342     fhDeltaCellClusterXNCells->SetXTitle("x (cm)");
2343     fhDeltaCellClusterXNCells->SetYTitle("N cells per cluster");
2344     outputContainer->Add(fhDeltaCellClusterXNCells);
2345
2346     fhDeltaCellClusterYNCells  = new TH2F ("hDeltaCellClusterYNCells","Cluster-Cell Y position vs N Cells per Cluster",ybins*2,-dy,dy,nceclbins,nceclmin,nceclmax);
2347     fhDeltaCellClusterYNCells->SetXTitle("y (cm)");
2348     fhDeltaCellClusterYNCells->SetYTitle("N cells per cluster");
2349     outputContainer->Add(fhDeltaCellClusterYNCells);
2350
2351     fhDeltaCellClusterZNCells  = new TH2F ("hDeltaCellClusterZNCells","Cluster-Cell Z position vs N Cells per Cluster",zbins*2,-dz,dz,nceclbins,nceclmin,nceclmax);
2352     fhDeltaCellClusterZNCells->SetXTitle("z (cm)");
2353     fhDeltaCellClusterZNCells->SetYTitle("N cells per cluster");
2354     outputContainer->Add(fhDeltaCellClusterZNCells);
2355
2356     fhDeltaCellClusterRE  = new TH2F ("hDeltaCellClusterRE","Cluster-Cell R position vs cluster energy",rbins*2,-dr,dr,nptbins,ptmin,ptmax);
2357     fhDeltaCellClusterRE->SetXTitle("r = #sqrt{x^{2}+y^{2}} (cm)");
2358     fhDeltaCellClusterRE->SetYTitle("E (GeV)");
2359     outputContainer->Add(fhDeltaCellClusterRE);
2360
2361     fhDeltaCellClusterXE  = new TH2F ("hDeltaCellClusterXE","Cluster-Cell X position vs cluster energy",xbins*2,-dx,dx,nptbins,ptmin,ptmax);
2362     fhDeltaCellClusterXE->SetXTitle("x (cm)");
2363     fhDeltaCellClusterXE->SetYTitle("E (GeV)");
2364     outputContainer->Add(fhDeltaCellClusterXE);
2365
2366     fhDeltaCellClusterYE  = new TH2F ("hDeltaCellClusterYE","Cluster-Cell Y position vs cluster energy",ybins*2,-dy,dy,nptbins,ptmin,ptmax);
2367     fhDeltaCellClusterYE->SetXTitle("y (cm)");
2368     fhDeltaCellClusterYE->SetYTitle("E (GeV)");
2369     outputContainer->Add(fhDeltaCellClusterYE);
2370
2371     fhDeltaCellClusterZE  = new TH2F ("hDeltaCellClusterZE","Cluster-Cell Z position vs cluster energy",zbins*2,-dz,dz,nptbins,ptmin,ptmax);
2372     fhDeltaCellClusterZE->SetXTitle("z (cm)");
2373     fhDeltaCellClusterZE->SetYTitle("E (GeV)");
2374     outputContainer->Add(fhDeltaCellClusterZE);
2375
2376     fhEtaPhiAmp  = new TH3F ("hEtaPhiAmp","Cell #eta vs cell #phi vs cell energy",netabins,etamin,etamax,nphibins,phimin,phimax,nptbins,ptmin,ptmax);
2377     fhEtaPhiAmp->SetXTitle("#eta ");
2378     fhEtaPhiAmp->SetYTitle("#phi (rad)");
2379     fhEtaPhiAmp->SetZTitle("E (GeV) ");
2380     outputContainer->Add(fhEtaPhiAmp);
2381
2382   }
2383
2384   //Calo cells
2385   fhNCells  = new TH1F ("hNCells","# cells", ncebins,ncemin+0.5,ncemax);
2386   fhNCells->SetXTitle("n cells");
2387   outputContainer->Add(fhNCells);
2388
2389   fhAmplitude  = new TH1F ("hAmplitude","Cell Energy", nptbins*2,ptmin,ptmax);
2390   fhAmplitude->SetXTitle("Cell Energy (GeV)");
2391   outputContainer->Add(fhAmplitude);
2392
2393   fhAmpId  = new TH2F ("hAmpId","Cell Energy", nfineptbins,ptfinemin,ptfinemax,fNMaxRows*fNMaxCols*fNModules,0,fNMaxRows*fNMaxCols*fNModules);
2394   fhAmpId->SetXTitle("Cell Energy (GeV)");
2395   outputContainer->Add(fhAmpId);
2396
2397   if(fFillAllCellTimeHisto)
2398   {
2399     fhCellTimeSpreadRespectToCellMax = new TH2F ("hCellTimeSpreadRespectToCellMax","t_{cell max}-t_{cell i} per cluster", nptbins,ptmin,ptmax,tdbins,tdmin,tdmax);
2400     fhCellTimeSpreadRespectToCellMax->SetXTitle("E (GeV)");
2401     fhCellTimeSpreadRespectToCellMax->SetYTitle("#Delta t_{cell max-i} (ns)");
2402     outputContainer->Add(fhCellTimeSpreadRespectToCellMax);
2403
2404     fhClusterMaxCellDiffAverageTime = new TH2F ("hClusterMaxCellDiffAverageTime","t_{cell max}-t_{average} per cluster", nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
2405     fhClusterMaxCellDiffAverageTime->SetXTitle("E (GeV)");
2406     fhClusterMaxCellDiffAverageTime->SetYTitle("#Delta t_{cell max - average} (ns)");
2407     outputContainer->Add(fhClusterMaxCellDiffAverageTime);
2408
2409     fhClusterMaxCellDiffWeightedTime = new TH2F ("hClusterMaxCellDiffWeightedTime","t_{cell max}-t_{weighted} per cluster", nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
2410     fhClusterMaxCellDiffWeightedTime->SetXTitle("E (GeV)");
2411     fhClusterMaxCellDiffWeightedTime->SetYTitle("#Delta t_{cell max - weighted} (ns)");
2412     outputContainer->Add(fhClusterMaxCellDiffWeightedTime);
2413
2414     fhCellIdCellLargeTimeSpread= new TH1F ("hCellIdCellLargeTimeSpread","Cells with time 100 ns larger than cell max in cluster ",
2415                                            fNMaxCols*fNMaxRows*fNModules,0,fNMaxCols*fNMaxRows*fNModules);
2416     fhCellIdCellLargeTimeSpread->SetXTitle("Absolute Cell Id");
2417     outputContainer->Add(fhCellIdCellLargeTimeSpread);
2418
2419     fhTime  = new TH1F ("hTime","Cell Time",ntimebins,timemin,timemax);
2420     fhTime->SetXTitle("Cell Time (ns)");
2421     outputContainer->Add(fhTime);
2422
2423     fhTimeVz  = new TH2F ("hTimeVz","Cell Time vs vertex, amplitude > 0.5 GeV",100, 0, 50,ntimebins,timemin,timemax);
2424     fhTimeVz->SetXTitle("|v_{z}| (cm)");
2425     fhTimeVz->SetYTitle("Cell Time (ns)");
2426     outputContainer->Add(fhTimeVz);
2427
2428     fhTimeId  = new TH2F ("hTimeId","Cell Time vs Absolute Id",
2429                           ntimebins,timemin,timemax,fNMaxRows*fNMaxCols*fNModules,0,fNMaxRows*fNMaxCols*fNModules);
2430     fhTimeId->SetXTitle("Cell Time (ns)");
2431     fhTimeId->SetYTitle("Cell Absolute Id");
2432     outputContainer->Add(fhTimeId);
2433
2434     fhTimeAmp  = new TH2F ("hTimeAmp","Cell Time vs Cell Energy",nptbins*2,ptmin,ptmax,ntimebins,timemin,timemax);
2435     fhTimeAmp->SetYTitle("Cell Time (ns)");
2436     fhTimeAmp->SetXTitle("Cell Energy (GeV)");
2437     outputContainer->Add(fhTimeAmp);
2438
2439   }
2440
2441   fhCellECross  = new TH2F ("hCellECross","1 - Energy in cross around cell /  cell energy",
2442                             nptbins,ptmin,ptmax, 400,-1,1.);
2443   fhCellECross->SetXTitle("E_{cell} (GeV) ");
2444   fhCellECross->SetYTitle("1- E_{cross}/E_{cell}");
2445   outputContainer->Add(fhCellECross);
2446
2447
2448   if(fCorrelate)
2449   {
2450     //PHOS vs EMCAL
2451     fhCaloCorrNClusters  = new TH2F ("hCaloCorrNClusters","# clusters in EMCAL vs PHOS", nclbins,nclmin,nclmax,nclbins,nclmin,nclmax);
2452     fhCaloCorrNClusters->SetXTitle("number of clusters in EMCAL");
2453     fhCaloCorrNClusters->SetYTitle("number of clusters in PHOS");
2454     outputContainer->Add(fhCaloCorrNClusters);
2455
2456     fhCaloCorrEClusters  = new TH2F ("hCaloCorrEClusters","summed energy of clusters in EMCAL vs PHOS", nptbins,ptmin,ptmax*2,nptbins,ptmin,ptmax*2);
2457     fhCaloCorrEClusters->SetXTitle("#Sigma E of clusters in EMCAL (GeV)");
2458     fhCaloCorrEClusters->SetYTitle("#Sigma E of clusters in PHOS (GeV)");
2459     outputContainer->Add(fhCaloCorrEClusters);
2460
2461     fhCaloCorrNCells  = new TH2F ("hCaloCorrNCells","# Cells in EMCAL vs PHOS", ncebins,ncemin,ncemax, ncebins,ncemin,ncemax);
2462     fhCaloCorrNCells->SetXTitle("number of Cells in EMCAL");
2463     fhCaloCorrNCells->SetYTitle("number of Cells in PHOS");
2464     outputContainer->Add(fhCaloCorrNCells);
2465
2466     fhCaloCorrECells  = new TH2F ("hCaloCorrECells","summed energy of Cells in EMCAL vs PHOS", nptbins*2,ptmin,ptmax*4,nptbins*2,ptmin,ptmax*4);
2467     fhCaloCorrECells->SetXTitle("#Sigma E of Cells in EMCAL (GeV)");
2468     fhCaloCorrECells->SetYTitle("#Sigma E of Cells in PHOS (GeV)");
2469     outputContainer->Add(fhCaloCorrECells);
2470
2471     //Calorimeter VS V0 signal
2472     fhCaloV0SCorrNClusters  = new TH2F ("hCaloV0SNClusters",Form("# clusters in %s vs V0 signal",fCalorimeter.Data()), nv0sbins,nv0smin,nv0smax,nclbins,nclmin,nclmax);
2473     fhCaloV0SCorrNClusters->SetXTitle("V0 signal");
2474     fhCaloV0SCorrNClusters->SetYTitle(Form("number of clusters in %s",fCalorimeter.Data()));
2475     outputContainer->Add(fhCaloV0SCorrNClusters);
2476
2477     fhCaloV0SCorrEClusters  = new TH2F ("hCaloV0SEClusters",Form("summed energy of clusters in %s vs V0 signal",fCalorimeter.Data()), nv0sbins,nv0smin,nv0smax,nptbins,ptmin,ptmax*2);
2478     fhCaloV0SCorrEClusters->SetXTitle("V0 signal");
2479     fhCaloV0SCorrEClusters->SetYTitle(Form("#Sigma E of clusters in %s (GeV)",fCalorimeter.Data()));
2480     outputContainer->Add(fhCaloV0SCorrEClusters);
2481
2482     fhCaloV0SCorrNCells  = new TH2F ("hCaloV0SNCells",Form("# Cells in %s vs V0 signal",fCalorimeter.Data()), nv0sbins,nv0smin,nv0smax, ncebins,ncemin,ncemax);
2483     fhCaloV0SCorrNCells->SetXTitle("V0 signal");
2484     fhCaloV0SCorrNCells->SetYTitle(Form("number of Cells in %s",fCalorimeter.Data()));
2485     outputContainer->Add(fhCaloV0SCorrNCells);
2486
2487     fhCaloV0SCorrECells  = new TH2F ("hCaloV0SECells",Form("summed energy of Cells in %s vs V0 signal",fCalorimeter.Data()), nv0sbins,nv0smin,nv0smax,nptbins,ptmin,ptmax*2);
2488     fhCaloV0SCorrECells->SetXTitle("V0 signal");
2489     fhCaloV0SCorrECells->SetYTitle(Form("#Sigma E of Cells in %s (GeV)",fCalorimeter.Data()));
2490     outputContainer->Add(fhCaloV0SCorrECells);
2491
2492     //Calorimeter VS V0 multiplicity
2493     fhCaloV0MCorrNClusters  = new TH2F ("hCaloV0MNClusters",Form("# clusters in %s vs V0 signal",fCalorimeter.Data()), nv0mbins,nv0mmin,nv0mmax,nclbins,nclmin,nclmax);
2494     fhCaloV0MCorrNClusters->SetXTitle("V0 signal");
2495     fhCaloV0MCorrNClusters->SetYTitle(Form("number of clusters in %s",fCalorimeter.Data()));
2496     outputContainer->Add(fhCaloV0MCorrNClusters);
2497
2498     fhCaloV0MCorrEClusters  = new TH2F ("hCaloV0MEClusters",Form("summed energy of clusters in %s vs V0 signal",fCalorimeter.Data()), nv0mbins,nv0mmin,nv0mmax,nptbins,ptmin,ptmax*2);
2499     fhCaloV0MCorrEClusters->SetXTitle("V0 signal");
2500     fhCaloV0MCorrEClusters->SetYTitle(Form("#Sigma E of clusters in %s (GeV)",fCalorimeter.Data()));
2501     outputContainer->Add(fhCaloV0MCorrEClusters);
2502
2503     fhCaloV0MCorrNCells  = new TH2F ("hCaloV0MNCells",Form("# Cells in %s vs V0 signal",fCalorimeter.Data()), nv0mbins,nv0mmin,nv0mmax, ncebins,ncemin,ncemax);
2504     fhCaloV0MCorrNCells->SetXTitle("V0 signal");
2505     fhCaloV0MCorrNCells->SetYTitle(Form("number of Cells in %s",fCalorimeter.Data()));
2506     outputContainer->Add(fhCaloV0MCorrNCells);
2507
2508     fhCaloV0MCorrECells  = new TH2F ("hCaloV0MECells",Form("summed energy of Cells in %s vs V0 signal",fCalorimeter.Data()), nv0mbins,nv0mmin,nv0mmax,nptbins,ptmin,ptmax*2);
2509     fhCaloV0MCorrECells->SetXTitle("V0 signal");
2510     fhCaloV0MCorrECells->SetYTitle(Form("#Sigma E of Cells in %s (GeV)",fCalorimeter.Data()));
2511     outputContainer->Add(fhCaloV0MCorrECells);
2512
2513     //Calorimeter VS Track multiplicity
2514     fhCaloTrackMCorrNClusters  = new TH2F ("hCaloTrackMNClusters",Form("# clusters in %s vs # tracks",fCalorimeter.Data()), ntrmbins,ntrmmin,ntrmmax,nclbins,nclmin,nclmax);
2515     fhCaloTrackMCorrNClusters->SetXTitle("# tracks");
2516     fhCaloTrackMCorrNClusters->SetYTitle(Form("number of clusters in %s",fCalorimeter.Data()));
2517     outputContainer->Add(fhCaloTrackMCorrNClusters);
2518
2519     fhCaloTrackMCorrEClusters  = new TH2F ("hCaloTrackMEClusters",Form("summed energy of clusters in %s vs # tracks",fCalorimeter.Data()), ntrmbins,ntrmmin,ntrmmax,nptbins,ptmin,ptmax*2);
2520     fhCaloTrackMCorrEClusters->SetXTitle("# tracks");
2521     fhCaloTrackMCorrEClusters->SetYTitle(Form("#Sigma E of clusters in %s (GeV)",fCalorimeter.Data()));
2522     outputContainer->Add(fhCaloTrackMCorrEClusters);
2523
2524     fhCaloTrackMCorrNCells  = new TH2F ("hCaloTrackMNCells",Form("# Cells in %s vs # tracks",fCalorimeter.Data()), ntrmbins,ntrmmin,ntrmmax, ncebins,ncemin,ncemax);
2525     fhCaloTrackMCorrNCells->SetXTitle("# tracks");
2526     fhCaloTrackMCorrNCells->SetYTitle(Form("number of Cells in %s",fCalorimeter.Data()));
2527     outputContainer->Add(fhCaloTrackMCorrNCells);
2528
2529     fhCaloTrackMCorrECells  = new TH2F ("hCaloTrackMECells",Form("summed energy of Cells in %s vs # tracks",fCalorimeter.Data()), ntrmbins,ntrmmin,ntrmmax,nptbins,ptmin,ptmax*2);
2530     fhCaloTrackMCorrECells->SetXTitle("# tracks");
2531     fhCaloTrackMCorrECells->SetYTitle(Form("#Sigma E of Cells in %s (GeV)",fCalorimeter.Data()));
2532     outputContainer->Add(fhCaloTrackMCorrECells);
2533
2534     fhCaloCenNClusters  = new TH2F ("hCaloCenNClusters","# clusters in calorimeter vs centrality",100,0,100,nclbins,nclmin,nclmax);
2535     fhCaloCenNClusters->SetYTitle("number of clusters in calorimeter");
2536     fhCaloCenNClusters->SetXTitle("Centrality");
2537     outputContainer->Add(fhCaloCenNClusters);
2538
2539     fhCaloCenEClusters  = new TH2F ("hCaloCenEClusters","summed energy of clusters in calorimeter vs centrality",100,0,100,nptbins,ptmin,ptmax*2);
2540     fhCaloCenEClusters->SetYTitle("#Sigma E of clusters in calorimeter (GeV)");
2541     fhCaloCenEClusters->SetXTitle("Centrality");
2542     outputContainer->Add(fhCaloCenEClusters);
2543
2544     fhCaloCenNCells  = new TH2F ("hCaloCenNCells","# Cells in calorimeter vs centrality",100,0,100,ncebins,ncemin,ncemax);
2545     fhCaloCenNCells->SetYTitle("number of Cells in calorimeter");
2546     fhCaloCenNCells->SetXTitle("Centrality");
2547     outputContainer->Add(fhCaloCenNCells);
2548
2549     fhCaloCenECells  = new TH2F ("hCaloCenECells","summed energy of Cells in calorimeter vs centrality",100,0,100,nptbins*2,ptmin,ptmax*4);
2550     fhCaloCenECells->SetYTitle("#Sigma E of Cells in calorimeter (GeV)");
2551     fhCaloCenECells->SetXTitle("Centrality");
2552     outputContainer->Add(fhCaloCenECells);
2553
2554     fhCaloEvPNClusters  = new TH2F ("hCaloEvPNClusters","# clusters in calorimeter vs event plane angle",100,0,TMath::Pi(),nclbins,nclmin,nclmax);
2555     fhCaloEvPNClusters->SetYTitle("number of clusters in calorimeter");
2556     fhCaloEvPNClusters->SetXTitle("Event plane angle (rad)");
2557     outputContainer->Add(fhCaloEvPNClusters);
2558
2559     fhCaloEvPEClusters  = new TH2F ("hCaloEvPEClusters","summed energy of clusters in calorimeter vs  event plane angle",100,0,TMath::Pi(),nptbins,ptmin,ptmax*2);
2560     fhCaloEvPEClusters->SetYTitle("#Sigma E of clusters in calorimeter (GeV)");
2561     fhCaloEvPEClusters->SetXTitle("Event plane angle (rad)");
2562     outputContainer->Add(fhCaloEvPEClusters);
2563
2564     fhCaloEvPNCells  = new TH2F ("hCaloEvPNCells","# Cells in calorimeter vs  event plane angle",100,0,TMath::Pi(),ncebins,ncemin,ncemax);
2565     fhCaloEvPNCells->SetYTitle("number of Cells in calorimeter");
2566     fhCaloEvPNCells->SetXTitle("Event plane angle (rad)");
2567     outputContainer->Add(fhCaloEvPNCells);
2568
2569     fhCaloEvPECells  = new TH2F ("hCaloEvPECells","summed energy of Cells in calorimeter vs  event plane angle",100,0,TMath::Pi(),nptbins*2,ptmin,ptmax*4);
2570     fhCaloEvPECells->SetYTitle("#Sigma E of Cells in calorimeter (GeV)");
2571     fhCaloEvPECells->SetXTitle("Event plane angle (rad)");
2572     outputContainer->Add(fhCaloEvPECells);
2573
2574
2575   }//correlate calorimeters
2576
2577   //Module histograms
2578
2579   fhEMod  = new TH2F ("hE_Mod","Cluster reconstructed Energy in each present Module",nptbins,ptmin,ptmax,fNModules,0,fNModules);
2580   fhEMod->SetXTitle("E (GeV)");
2581   fhEMod->SetYTitle("Module");
2582   outputContainer->Add(fhEMod);
2583
2584   fhAmpMod  = new TH2F ("hAmp_Mod","Cell energy in each present Module",nptbins,ptmin,ptmax,fNModules,0,fNModules);
2585   fhAmpMod->SetXTitle("E (GeV)");
2586   fhAmpMod->SetYTitle("Module");
2587   outputContainer->Add(fhAmpMod);
2588
2589   if(fFillAllCellTimeHisto)
2590   {
2591     fhTimeMod  = new TH2F ("hTime_Mod","Cell time in each present Module",ntimebins,timemin,timemax,fNModules,0,fNModules);
2592     fhTimeMod->SetXTitle("t (ns)");
2593     fhTimeMod->SetYTitle("Module");
2594     outputContainer->Add(fhTimeMod);
2595   }
2596
2597   fhNClustersMod  = new TH2F ("hNClusters_Mod","# clusters vs Module", nclbins,nclmin+0.5,nclmax,fNModules,0,fNModules);
2598   fhNClustersMod->SetXTitle("number of clusters");
2599   fhNClustersMod->SetYTitle("Module");
2600   outputContainer->Add(fhNClustersMod);
2601
2602   fhNCellsMod  = new TH2F ("hNCells_Mod","# cells vs Module", ncebins,ncemin+0.5,ncemax,fNModules,0,fNModules);
2603   fhNCellsMod->SetXTitle("n cells");
2604   fhNCellsMod->SetYTitle("Module");
2605   outputContainer->Add(fhNCellsMod);
2606
2607   Int_t colmaxs = fNMaxCols;
2608   Int_t rowmaxs = fNMaxRows;
2609   if(fCalorimeter=="EMCAL")
2610   {
2611     colmaxs=2*fNMaxCols;
2612     rowmaxs=Int_t(fNModules/2)*fNMaxRows;
2613   }
2614   else
2615   {
2616     rowmaxs=fNModules*fNMaxRows;
2617   }
2618
2619   fhGridCells  = new TH2F ("hGridCells",Form("Entries in grid of cells"),
2620                            colmaxs+2,-1.5,colmaxs+0.5, rowmaxs+2,-1.5,rowmaxs+0.5);
2621   fhGridCells->SetYTitle("row (phi direction)");
2622   fhGridCells->SetXTitle("column (eta direction)");
2623   outputContainer->Add(fhGridCells);
2624
2625   fhGridCellsE  = new TH2F ("hGridCellsE","Accumulated energy in grid of cells",
2626                             colmaxs+2,-1.5,colmaxs+0.5, rowmaxs+2,-1.5,rowmaxs+0.5);
2627   fhGridCellsE->SetYTitle("row (phi direction)");
2628   fhGridCellsE->SetXTitle("column (eta direction)");
2629   outputContainer->Add(fhGridCellsE);
2630
2631   if(fFillAllCellTimeHisto)
2632   {
2633     fhGridCellsTime  = new TH2F ("hGridCellsTime","Accumulated time in grid of cells",
2634                                  colmaxs+2,-1.5,colmaxs+0.5, rowmaxs+2,-1.5,rowmaxs+0.5);
2635     fhGridCellsTime->SetYTitle("row (phi direction)");
2636     fhGridCellsTime->SetXTitle("column (eta direction)");
2637     outputContainer->Add(fhGridCellsTime);
2638   }
2639
2640   fhNCellsPerClusterMod      = new TH2F*[fNModules];
2641   fhNCellsPerClusterModNoCut = new TH2F*[fNModules];
2642   fhIMMod                    = new TH2F*[fNModules];
2643   if(fFillAllCellTimeHisto) fhTimeAmpPerRCU = new TH2F*[fNModules*fNRCU];
2644
2645   for(Int_t imod = 0; imod < fNModules; imod++)
2646   {
2647     fhNCellsPerClusterMod[imod]  = new TH2F (Form("hNCellsPerCluster_Mod%d",imod),
2648                                              Form("# cells per cluster vs cluster energy in Module %d",imod),
2649                                              nptbins,ptmin,ptmax, nceclbins,nceclmin,nceclmax);
2650     fhNCellsPerClusterMod[imod]->SetXTitle("E (GeV)");
2651     fhNCellsPerClusterMod[imod]->SetYTitle("n cells");
2652     outputContainer->Add(fhNCellsPerClusterMod[imod]);
2653
2654     fhNCellsPerClusterModNoCut[imod]  = new TH2F (Form("hNCellsPerClusterNoCut_Mod%d",imod),
2655                                                   Form("# cells per cluster vs cluster energy in Module %d, no cut",imod),
2656                                                   nptbins,ptmin,ptmax, nceclbins,nceclmin,nceclmax);
2657     fhNCellsPerClusterModNoCut[imod]->SetXTitle("E (GeV)");
2658     fhNCellsPerClusterModNoCut[imod]->SetYTitle("n cells");
2659     outputContainer->Add(fhNCellsPerClusterModNoCut[imod]);
2660
2661     if(fFillAllCellTimeHisto)
2662     {
2663       for(Int_t ircu = 0; ircu < fNRCU; ircu++)
2664       {
2665         fhTimeAmpPerRCU[imod*fNRCU+ircu]  = new TH2F (Form("hTimeAmp_Mod%d_RCU%d",imod,ircu),
2666                                                       Form("Cell Energy vs Cell Time in Module %d, RCU %d ",imod,ircu),
2667                                                       nptbins,ptmin,ptmax,ntimebins,timemin,timemax);
2668         fhTimeAmpPerRCU[imod*fNRCU+ircu]->SetXTitle("E (GeV)");
2669         fhTimeAmpPerRCU[imod*fNRCU+ircu]->SetYTitle("time (ns)");
2670         outputContainer->Add(fhTimeAmpPerRCU[imod*fNRCU+ircu]);
2671
2672       }
2673     }
2674
2675     if(fFillAllPi0Histo)
2676     {
2677       fhIMMod[imod]  = new TH2F (Form("hIM_Mod%d",imod),
2678                                  Form("Cluster pairs Invariant mass vs reconstructed pair energy in Module %d, n cell > 1",imod),
2679                                  nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
2680       fhIMMod[imod]->SetXTitle("p_{T, cluster pairs} (GeV) ");
2681       fhIMMod[imod]->SetYTitle("M_{cluster pairs} (GeV/c^{2})");
2682       outputContainer->Add(fhIMMod[imod]);
2683
2684     }
2685   }
2686
2687   // Monte Carlo Histograms
2688
2689   TString particleName[] = { "Photon", "Pi0", "Eta", "Electron", "NeutralHadron", "ChargedHadron" };
2690
2691   if(IsDataMC())
2692   {
2693     for(Int_t iPart = 0; iPart < 6; iPart++)
2694     {
2695       for(Int_t iCh = 0; iCh < 2; iCh++)
2696       {
2697         fhRecoMCRatioE[iPart][iCh]  = new TH2F (Form("hRecoMCRatioE_%s_Match%d",particleName[iPart].Data(),iCh),
2698                                                 Form("Reconstructed/Generated E, %s, Matched %d",particleName[iPart].Data(),iCh),
2699                                                 nptbins, ptmin, ptmax, 200,0,2);
2700         fhRecoMCRatioE[iPart][iCh]->SetYTitle("E_{reconstructed}/E_{generated}");
2701         fhRecoMCRatioE[iPart][iCh]->SetXTitle("E_{reconstructed} (GeV)");
2702         outputContainer->Add(fhRecoMCRatioE[iPart][iCh]);
2703
2704
2705         fhRecoMCDeltaE[iPart][iCh]  = new TH2F (Form("hRecoMCDeltaE_%s_Match%d",particleName[iPart].Data(),iCh),
2706                                                 Form("Generated - Reconstructed E, %s, Matched %d",particleName[iPart].Data(),iCh),
2707                                                 nptbins, ptmin, ptmax, nptbins*2,-ptmax,ptmax);
2708         fhRecoMCDeltaE[iPart][iCh]->SetYTitle("#Delta E (GeV)");
2709         fhRecoMCDeltaE[iPart][iCh]->SetXTitle("E_{reconstructed} (GeV)");
2710         outputContainer->Add(fhRecoMCDeltaE[iPart][iCh]);
2711
2712         fhRecoMCDeltaPhi[iPart][iCh]  = new TH2F (Form("hRecoMCDeltaPhi_%s_Match%d",particleName[iPart].Data(),iCh),
2713                                                   Form("Generated - Reconstructed #phi, %s, Matched %d",particleName[iPart].Data(),iCh),
2714                                                   nptbins, ptmin, ptmax, nphibins*2,-phimax,phimax);
2715         fhRecoMCDeltaPhi[iPart][iCh]->SetYTitle("#Delta #phi (rad)");
2716         fhRecoMCDeltaPhi[iPart][iCh]->SetXTitle("E_{reconstructed} (GeV)");
2717         outputContainer->Add(fhRecoMCDeltaPhi[iPart][iCh]);
2718
2719         fhRecoMCDeltaEta[iPart][iCh]  = new TH2F (Form("hRecoMCDeltaEta_%s_Match%d",particleName[iPart].Data(),iCh),
2720                                                   Form("Generated - Reconstructed #eta, %s, Matched %d",particleName[iPart].Data(),iCh),
2721                                                   nptbins, ptmin, ptmax,netabins*2,-etamax,etamax);
2722         fhRecoMCDeltaEta[iPart][iCh]->SetYTitle("#Delta #eta ");
2723         fhRecoMCDeltaEta[iPart][iCh]->SetXTitle("E_{reconstructed} (GeV)");
2724         outputContainer->Add(fhRecoMCDeltaEta[iPart][iCh]);
2725
2726         fhRecoMCE[iPart][iCh]  = new TH2F (Form("hRecoMCE_%s_Match%d",particleName[iPart].Data(),iCh),
2727                                            Form("E distribution, reconstructed vs generated, %s, Matched %d",particleName[iPart].Data(),iCh),
2728                                            nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2729         fhRecoMCE[iPart][iCh]->SetXTitle("E_{rec} (GeV)");
2730         fhRecoMCE[iPart][iCh]->SetYTitle("E_{gen} (GeV)");
2731         outputContainer->Add(fhRecoMCE[iPart][iCh]);
2732
2733         fhRecoMCPhi[iPart][iCh]  = new TH2F (Form("hRecoMCPhi_%s_Match%d",particleName[iPart].Data(),iCh),
2734                                              Form("#phi distribution, reconstructed vs generated, %s, Matched %d",particleName[iPart].Data(),iCh),
2735                                              nphibins,phimin,phimax, nphibins,phimin,phimax);
2736         fhRecoMCPhi[iPart][iCh]->SetXTitle("#phi_{reconstructed} (rad)");
2737         fhRecoMCPhi[iPart][iCh]->SetYTitle("#phi_{generated} (rad)");
2738         outputContainer->Add(fhRecoMCPhi[iPart][iCh]);
2739
2740         fhRecoMCEta[iPart][iCh]  = new TH2F (Form("hRecoMCEta_%s_Match%d",particleName[iPart].Data(),iCh),
2741                                              Form("#eta distribution, reconstructed vs generated, %s, Matched %d",particleName[iPart].Data(),iCh),
2742                                              netabins,etamin,etamax,netabins,etamin,etamax);
2743         fhRecoMCEta[iPart][iCh]->SetXTitle("#eta_{reconstructed} ");
2744         fhRecoMCEta[iPart][iCh]->SetYTitle("#eta_{generated} ");
2745         outputContainer->Add(fhRecoMCEta[iPart][iCh]);
2746       }
2747     }
2748
2749     //Pure MC
2750     for(Int_t iPart = 0; iPart < 4; iPart++)
2751     {
2752       fhGenMCE[iPart]     = new TH1F(Form("hGenMCE_%s",particleName[iPart].Data()) ,
2753                                      Form("p_{T} of generated %s",particleName[iPart].Data()),
2754                                      nptbins,ptmin,ptmax);
2755       fhGenMCEtaPhi[iPart] = new TH2F(Form("hGenMCEtaPhi_%s",particleName[iPart].Data()),
2756                                       Form("Y vs #phi of generated %s",particleName[iPart].Data()),
2757                                       netabins,etamin,etamax,nphibins,phimin,phimax);
2758
2759       fhGenMCE[iPart]    ->SetXTitle("p_{T} (GeV/c)");
2760       fhGenMCEtaPhi[iPart]->SetXTitle("#eta");
2761       fhGenMCEtaPhi[iPart]->SetYTitle("#phi (rad)");
2762
2763       outputContainer->Add(fhGenMCE[iPart]);
2764       outputContainer->Add(fhGenMCEtaPhi[iPart]);
2765
2766
2767       fhGenMCAccE[iPart]     = new TH1F(Form("hGenMCAccE_%s",particleName[iPart].Data()) ,
2768                                         Form("p_{T} of generated %s",particleName[iPart].Data()),
2769                                         nptbins,ptmin,ptmax);
2770       fhGenMCAccEtaPhi[iPart] = new TH2F(Form("hGenMCAccEtaPhi_%s",particleName[iPart].Data()),
2771                                          Form("Y vs #phi of generated %s",particleName[iPart].Data()),
2772                                          netabins,etamin,etamax,nphibins,phimin,phimax);
2773
2774       fhGenMCAccE[iPart]    ->SetXTitle("p_{T} (GeV/c)");
2775       fhGenMCAccEtaPhi[iPart]->SetXTitle("#eta");
2776       fhGenMCAccEtaPhi[iPart]->SetYTitle("#phi (rad)");
2777
2778       outputContainer->Add(fhGenMCAccE[iPart]);
2779       outputContainer->Add(fhGenMCAccEtaPhi[iPart]);
2780
2781     }
2782
2783     //Vertex of generated particles
2784
2785     fhEMVxyz  = new TH2F ("hEMVxyz","Production vertex of reconstructed ElectroMagnetic particles",nvdistbins,vdistmin,vdistmax,nvdistbins,vdistmin,vdistmax);//,100,0,500);
2786     fhEMVxyz->SetXTitle("v_{x}");
2787     fhEMVxyz->SetYTitle("v_{y}");
2788     //fhEMVxyz->SetZTitle("v_{z}");
2789     outputContainer->Add(fhEMVxyz);
2790
2791     fhHaVxyz  = new TH2F ("hHaVxyz","Production vertex of reconstructed hadrons",nvdistbins,vdistmin,vdistmax,nvdistbins,vdistmin,vdistmax);//,100,0,500);
2792     fhHaVxyz->SetXTitle("v_{x}");
2793     fhHaVxyz->SetYTitle("v_{y}");
2794     //fhHaVxyz->SetZTitle("v_{z}");
2795     outputContainer->Add(fhHaVxyz);
2796
2797     fhEMR  = new TH2F ("hEMR","Distance to production vertex of reconstructed ElectroMagnetic particles vs E rec",nptbins,ptmin,ptmax,nvdistbins,vdistmin,vdistmax);
2798     fhEMR->SetXTitle("E (GeV)");
2799     fhEMR->SetYTitle("TMath::Sqrt(v_{x}^{2}+v_{y}^{2})");
2800     outputContainer->Add(fhEMR);
2801
2802     fhHaR  = new TH2F ("hHaR","Distance to production vertex of reconstructed Hadrons vs E rec",nptbins,ptmin,ptmax,nvdistbins,vdistmin,vdistmax);
2803     fhHaR->SetXTitle("E (GeV)");
2804     fhHaR->SetYTitle("TMath::Sqrt(v_{x}^{2}+v_{y}^{2})");
2805     outputContainer->Add(fhHaR);
2806
2807
2808     //Track Matching
2809
2810     fhMCEle1EOverP = new TH2F("hMCEle1EOverP","TRACK matches E/p, MC electrons",nptbins,ptmin,ptmax, nPoverEbins,eOverPmin,eOverPmax);
2811     fhMCEle1EOverP->SetYTitle("E/p");
2812     fhMCEle1EOverP->SetXTitle("p_{T} (GeV/c)");
2813     outputContainer->Add(fhMCEle1EOverP);
2814
2815     fhMCEle1dR = new TH1F("hMCEle1dR","TRACK matches dR, MC electrons",ndRbins,dRmin,dRmax);
2816     fhMCEle1dR->SetXTitle("#Delta R (rad)");
2817     outputContainer->Add(fhMCEle1dR) ;
2818
2819     fhMCEle2MatchdEdx = new TH2F("hMCEle2MatchdEdx","dE/dx vs. p for all matches, MC electrons",nptbins,ptmin,ptmax,ndedxbins,dedxmin,dedxmax);
2820     fhMCEle2MatchdEdx->SetXTitle("p (GeV/c)");
2821     fhMCEle2MatchdEdx->SetYTitle("<dE/dx>");
2822     outputContainer->Add(fhMCEle2MatchdEdx);
2823
2824     fhMCChHad1EOverP = new TH2F("hMCChHad1EOverP","TRACK matches E/p, MC charged hadrons",nptbins,ptmin,ptmax, nPoverEbins,eOverPmin,eOverPmax);
2825     fhMCChHad1EOverP->SetYTitle("E/p");
2826     fhMCChHad1EOverP->SetXTitle("p_{T} (GeV/c)");
2827     outputContainer->Add(fhMCChHad1EOverP);
2828
2829     fhMCChHad1dR = new TH1F("hMCChHad1dR","TRACK matches dR, MC charged hadrons",ndRbins,dRmin,dRmax);
2830     fhMCChHad1dR->SetXTitle("#Delta R (rad)");
2831     outputContainer->Add(fhMCChHad1dR) ;
2832
2833     fhMCChHad2MatchdEdx = new TH2F("hMCChHad2MatchdEdx","dE/dx vs. p for all matches, MC charged hadrons",nptbins,ptmin,ptmax,ndedxbins,dedxmin,dedxmax);
2834     fhMCChHad2MatchdEdx->SetXTitle("p (GeV/c)");
2835     fhMCChHad2MatchdEdx->SetYTitle("<dE/dx>");
2836     outputContainer->Add(fhMCChHad2MatchdEdx);
2837
2838     fhMCNeutral1EOverP = new TH2F("hMCNeutral1EOverP","TRACK matches E/p, MC neutrals",nptbins,ptmin,ptmax, nPoverEbins,eOverPmin,eOverPmax);
2839     fhMCNeutral1EOverP->SetYTitle("E/p");
2840     fhMCNeutral1EOverP->SetXTitle("p_{T} (GeV/c)");
2841     outputContainer->Add(fhMCNeutral1EOverP);
2842
2843     fhMCNeutral1dR = new TH1F("hMCNeutral1dR","TRACK matches dR, MC neutrals",ndRbins,dRmin,dRmax);
2844     fhMCNeutral1dR->SetXTitle("#Delta R (rad)");
2845     outputContainer->Add(fhMCNeutral1dR) ;
2846
2847     fhMCNeutral2MatchdEdx = new TH2F("hMCNeutral2MatchdEdx","dE/dx vs. p for all matches, MC neutrals",nptbins,ptmin,ptmax,ndedxbins,dedxmin,dedxmax);
2848     fhMCNeutral2MatchdEdx->SetXTitle("p (GeV/c)");
2849     fhMCNeutral2MatchdEdx->SetYTitle("<dE/dx>");
2850     outputContainer->Add(fhMCNeutral2MatchdEdx);
2851
2852     fhMCEle1EOverPR02 = new TH2F("hMCEle1EOverPR02","TRACK matches E/p, MC electrons",nptbins,ptmin,ptmax, nPoverEbins,eOverPmin,eOverPmax);
2853     fhMCEle1EOverPR02->SetYTitle("E/p");
2854     fhMCEle1EOverPR02->SetXTitle("p_{T} (GeV/c)");
2855     outputContainer->Add(fhMCEle1EOverPR02);
2856
2857     fhMCChHad1EOverPR02 = new TH2F("hMCChHad1EOverPR02","TRACK matches E/p, MC charged hadrons",nptbins,ptmin,ptmax, nPoverEbins,eOverPmin,eOverPmax);
2858     fhMCChHad1EOverPR02->SetYTitle("E/p");
2859     fhMCChHad1EOverPR02->SetXTitle("p_{T} (GeV/c)");
2860     outputContainer->Add(fhMCChHad1EOverPR02);
2861
2862     fhMCNeutral1EOverPR02 = new TH2F("hMCNeutral1EOverPR02","TRACK matches E/p, MC neutrals",nptbins,ptmin,ptmax, nPoverEbins,eOverPmin,eOverPmax);
2863     fhMCNeutral1EOverPR02->SetYTitle("E/p");
2864     fhMCNeutral1EOverPR02->SetXTitle("p_{T} (GeV/c)");
2865     outputContainer->Add(fhMCNeutral1EOverPR02);
2866
2867     fhMCEle1EleEOverP = new TH2F("hMCEle1EleEOverP","Electron candidates E/p (60<dEdx<100), MC electrons",nptbins,ptmin,ptmax, nPoverEbins,eOverPmin,eOverPmax);
2868     fhMCEle1EleEOverP->SetYTitle("E/p");
2869     fhMCEle1EleEOverP->SetXTitle("p_{T} (GeV/c)");
2870     outputContainer->Add(fhMCEle1EleEOverP);
2871
2872     fhMCChHad1EleEOverP = new TH2F("hMCEle1EleEOverP","Electron candidates E/p (60<dEdx<100), MC charged hadrons",nptbins,ptmin,ptmax, nPoverEbins,eOverPmin,eOverPmax);
2873     fhMCChHad1EleEOverP->SetYTitle("E/p");
2874     fhMCChHad1EleEOverP->SetXTitle("p_{T} (GeV/c)");
2875     outputContainer->Add(fhMCChHad1EleEOverP);
2876
2877     fhMCNeutral1EleEOverP = new TH2F("hMCNeutral1EleEOverP","Electron candidates E/p (60<dEdx<100), MC neutrals",nptbins,ptmin,ptmax, nPoverEbins,eOverPmin,eOverPmax);
2878     fhMCNeutral1EleEOverP->SetYTitle("E/p");
2879     fhMCNeutral1EleEOverP->SetXTitle("p_{T} (GeV/c)");
2880     outputContainer->Add(fhMCNeutral1EleEOverP);
2881
2882   }
2883
2884   //  for(Int_t i = 0; i < outputContainer->GetEntries() ; i++)
2885   //    printf("i=%d, name= %s\n",i,outputContainer->At(i)->GetName());
2886
2887   return outputContainer;
2888 }
2889
2890 //______________________________________________________________________________________
2891 Float_t AliAnaCalorimeterQA::GetECross(Int_t absID, AliVCaloCells* cells, Float_t dtcut)
2892 {
2893   // Get energy in cross axis around maximum cell, for EMCAL only
2894
2895   Int_t icol =-1, irow=-1,iRCU = -1;
2896   Int_t imod = GetModuleNumberCellIndexes(absID, fCalorimeter, icol, irow, iRCU);
2897
2898   if(fCalorimeter=="EMCAL")
2899   {
2900     //Get close cells index, energy and time, not in corners
2901
2902     Int_t absID1 = -1;
2903     Int_t absID2 = -1;
2904
2905     if( irow < AliEMCALGeoParams::fgkEMCALRows-1) absID1 = GetCaloUtils()->GetEMCALGeometry()->GetAbsCellIdFromCellIndexes(imod, irow+1, icol);
2906     if( irow > 0 )                                absID2 = GetCaloUtils()->GetEMCALGeometry()->GetAbsCellIdFromCellIndexes(imod, irow-1, icol);
2907
2908     // In case of cell in eta = 0 border, depending on SM shift the cross cell index
2909     Int_t absID3 = -1;
2910     Int_t absID4 = -1;
2911
2912     if     ( icol == AliEMCALGeoParams::fgkEMCALCols - 1 && !(imod%2) )
2913     {
2914       absID3 = GetCaloUtils()->GetEMCALGeometry()-> GetAbsCellIdFromCellIndexes(imod+1, irow, 0);
2915       absID4 = GetCaloUtils()->GetEMCALGeometry()-> GetAbsCellIdFromCellIndexes(imod  , irow, icol-1);
2916     }
2917     else if( icol == 0 && imod%2 )
2918     {
2919       absID3 = GetCaloUtils()->GetEMCALGeometry()-> GetAbsCellIdFromCellIndexes(imod  , irow, icol+1);
2920       absID4 = GetCaloUtils()->GetEMCALGeometry()-> GetAbsCellIdFromCellIndexes(imod-1, irow, AliEMCALGeoParams::fgkEMCALCols-1);
2921     }
2922     else
2923     {
2924       if( icol < AliEMCALGeoParams::fgkEMCALCols-1 )
2925         absID3 = GetCaloUtils()->GetEMCALGeometry()-> GetAbsCellIdFromCellIndexes(imod, irow, icol+1);
2926       if( icol > 0 )
2927         absID4 = GetCaloUtils()->GetEMCALGeometry()-> GetAbsCellIdFromCellIndexes(imod, irow, icol-1);
2928     }
2929
2930     //Recalibrate cell energy if needed
2931     //Float_t  ecell = cells->GetCellAmplitude(absID);
2932     //GetCaloUtils()->RecalibrateCellAmplitude(ecell,fCalorimeter, absID);
2933     Double_t tcell = cells->GetCellTime(absID);
2934     GetCaloUtils()->RecalibrateCellTime(tcell, fCalorimeter, absID,GetReader()->GetInputEvent()->GetBunchCrossNumber());
2935
2936     Float_t  ecell1  = 0, ecell2  = 0, ecell3  = 0, ecell4  = 0;
2937     Double_t tcell1  = 0, tcell2  = 0, tcell3  = 0, tcell4  = 0;
2938
2939     if(absID1 >0 )
2940     {
2941       ecell1 = cells->GetCellAmplitude(absID1);
2942       GetCaloUtils()->RecalibrateCellAmplitude(ecell1, fCalorimeter, absID1);
2943       tcell1 = cells->GetCellTime(absID1);
2944       GetCaloUtils()->RecalibrateCellTime     (tcell1, fCalorimeter, absID1,GetReader()->GetInputEvent()->GetBunchCrossNumber());
2945     }
2946     if(absID2 >0 )
2947     {
2948       ecell2 = cells->GetCellAmplitude(absID2);
2949       GetCaloUtils()->RecalibrateCellAmplitude(ecell2, fCalorimeter, absID2);
2950       tcell2 = cells->GetCellTime(absID2);
2951       GetCaloUtils()->RecalibrateCellTime     (tcell2, fCalorimeter, absID2, GetReader()->GetInputEvent()->GetBunchCrossNumber());
2952     }
2953     if(absID3 >0 )
2954     {
2955       ecell3 = cells->GetCellAmplitude(absID3);
2956       GetCaloUtils()->RecalibrateCellAmplitude(ecell3, fCalorimeter, absID3);
2957       tcell3 = cells->GetCellTime(absID3);
2958       GetCaloUtils()->RecalibrateCellTime     (tcell3, fCalorimeter, absID3, GetReader()->GetInputEvent()->GetBunchCrossNumber());
2959     }
2960     if(absID4 >0 )
2961     {
2962       ecell4 = cells->GetCellAmplitude(absID4);
2963       GetCaloUtils()->RecalibrateCellAmplitude(ecell4, fCalorimeter, absID4);
2964       tcell4 = cells->GetCellTime(absID4);
2965       GetCaloUtils()->RecalibrateCellTime     (tcell4, fCalorimeter, absID4, GetReader()->GetInputEvent()->GetBunchCrossNumber());
2966     }
2967
2968     if(TMath::Abs(tcell-tcell1)*1.e9 > dtcut) ecell1 = 0 ;
2969     if(TMath::Abs(tcell-tcell2)*1.e9 > dtcut) ecell2 = 0 ;
2970     if(TMath::Abs(tcell-tcell3)*1.e9 > dtcut) ecell3 = 0 ;
2971     if(TMath::Abs(tcell-tcell4)*1.e9 > dtcut) ecell4 = 0 ;
2972
2973     return ecell1+ecell2+ecell3+ecell4;
2974   }
2975   else //PHOS
2976   {
2977
2978     Int_t absId1 = -1, absId2 = -1, absId3 = -1, absId4 = -1;
2979
2980     Int_t relId1[] = { imod+1, 0, irow+1, icol   };
2981     Int_t relId2[] = { imod+1, 0, irow-1, icol   };
2982     Int_t relId3[] = { imod+1, 0, irow  , icol+1 };
2983     Int_t relId4[] = { imod+1, 0, irow  , icol-1 };
2984
2985     GetCaloUtils()->GetPHOSGeometry()->RelToAbsNumbering(relId1, absId1);
2986     GetCaloUtils()->GetPHOSGeometry()->RelToAbsNumbering(relId2, absId2);
2987     GetCaloUtils()->GetPHOSGeometry()->RelToAbsNumbering(relId3, absId3);
2988     GetCaloUtils()->GetPHOSGeometry()->RelToAbsNumbering(relId4, absId4);
2989
2990     Float_t  ecell1  = 0, ecell2  = 0, ecell3  = 0, ecell4  = 0;
2991
2992     if(absId1 > 0 ) ecell1 = cells->GetCellAmplitude(absId1);
2993     if(absId2 > 0 ) ecell2 = cells->GetCellAmplitude(absId2);
2994     if(absId3 > 0 ) ecell3 = cells->GetCellAmplitude(absId3);
2995     if(absId4 > 0 ) ecell4 = cells->GetCellAmplitude(absId4);
2996
2997     return ecell1+ecell2+ecell3+ecell4;
2998
2999   }
3000
3001 }
3002
3003 //__________________________________________________________________________________________________
3004 void AliAnaCalorimeterQA::InvariantMassHistograms(Int_t iclus,   TLorentzVector mom,
3005                                                   Int_t nModule, const TObjArray* caloClusters,
3006                                                   AliVCaloCells * cells)
3007 {
3008   // Fill Invariant mass histograms
3009
3010   if(GetDebug()>1) printf("AliAnaCalorimeterQA::InvariantMassHistograms() - Start \n");
3011
3012   //Get vertex for photon momentum calculation and event selection
3013   Double_t v[3] = {0,0,0}; //vertex ;
3014   //GetReader()->GetVertex(v);
3015
3016   Int_t nModule2      = -1;
3017   TLorentzVector mom2 ;
3018   Int_t nCaloClusters = caloClusters->GetEntriesFast();
3019
3020   for(Int_t jclus = iclus + 1 ; jclus < nCaloClusters ; jclus++)
3021   {
3022     AliVCluster* clus2 =  (AliVCluster*)caloClusters->At(jclus);
3023
3024     Float_t maxCellFraction = 0.;
3025     Int_t absIdMax = GetCaloUtils()->GetMaxEnergyCell(cells, clus2,maxCellFraction);
3026
3027     // Try to rediuce background with a mild shower shape cut and no more than 1 maxima
3028     // in cluster and remove low energy clusters
3029     if( clus2->GetNCells() <= 1 || !IsGoodCluster(absIdMax,cells) ||
3030        GetCaloUtils()->GetNumberOfLocalMaxima(clus2,cells) > 1 ||
3031        clus2->GetM02() > 0.5 || clus2->E() < 0.3) continue;
3032
3033     //Get cluster kinematics
3034     clus2->GetMomentum(mom2,v);
3035
3036     //Check only certain regions
3037     Bool_t in2 = kTRUE;
3038     if(IsFiducialCutOn()) in2 =  GetFiducialCut()->IsInFiducialCut(mom2,fCalorimeter) ;
3039     if(!in2) continue;
3040
3041     //Get module of cluster
3042     nModule2 = GetModuleNumber(clus2);
3043
3044     //Fill histograms
3045
3046     //All modules
3047     fhIM  ->Fill((mom+mom2).Pt(),(mom+mom2).M());
3048
3049     //Single module
3050     if(nModule == nModule2 && nModule >=0 && nModule < fNModules)
3051       fhIMMod[nModule]->Fill((mom+mom2).Pt(),(mom+mom2).M());
3052
3053
3054     //Asymetry histograms
3055     fhAsym->Fill((mom+mom2).Pt(),TMath::Abs((mom.E()-mom2.E())/(mom.E()+mom2.E())));
3056
3057   }// 2nd cluster loop
3058
3059 }
3060
3061 //______________________________
3062 void AliAnaCalorimeterQA::Init()
3063 {
3064   //Check if the data or settings are ok
3065
3066   if(fCalorimeter != "PHOS" && fCalorimeter !="EMCAL")
3067     AliFatal(Form("Wrong calorimeter name <%s>", fCalorimeter.Data()));
3068
3069   //if(GetReader()->GetDataType()== AliCaloTrackReader::kMC)
3070   //  AliFatal("Analysis of reconstructed data, MC reader not aplicable");
3071
3072 }
3073
3074 //________________________________________
3075 void AliAnaCalorimeterQA::InitParameters()
3076 {
3077   //Initialize the parameters of the analysis.
3078   AddToHistogramsName("AnaCaloQA_");
3079
3080   fCalorimeter     = "EMCAL"; //or PHOS
3081   fNModules        = 12; // set maximum to maximum number of EMCAL modules
3082   fNRCU            = 2;  // set maximum number of RCU in EMCAL per SM
3083   fTimeCutMin      = -9999999;
3084   fTimeCutMax      =  9999999;
3085   fEMCALCellAmpMin = 0.2;
3086   fPHOSCellAmpMin  = 0.2;
3087
3088   // Exotic studies
3089   fExoNECrossCuts  = 10 ;
3090   fExoNDTimeCuts   = 4  ;
3091
3092   fExoDTimeCuts [0] = 1.e4 ; fExoDTimeCuts [1] = 50.0 ; fExoDTimeCuts [2] = 25.0 ; fExoDTimeCuts [3] = 10.0 ;
3093   fExoECrossCuts[0] = 0.80 ; fExoECrossCuts[1] = 0.85 ; fExoECrossCuts[2] = 0.90 ; fExoECrossCuts[3] = 0.92 ; fExoECrossCuts[4] = 0.94 ;
3094   fExoECrossCuts[5] = 0.95 ; fExoECrossCuts[6] = 0.96 ; fExoECrossCuts[7] = 0.97 ; fExoECrossCuts[8] = 0.98 ; fExoECrossCuts[9] = 0.99 ;
3095
3096 }
3097
3098 //_____________________________________________________________________________
3099 Bool_t AliAnaCalorimeterQA::IsGoodCluster(Int_t absIdMax, AliVCaloCells* cells)
3100 {
3101   //Identify cluster as exotic or not
3102
3103   if(!fStudyBadClusters) return kTRUE;
3104
3105   if(fCalorimeter=="EMCAL")
3106   {
3107     if(!GetCaloUtils()->GetEMCALRecoUtils()->IsRejectExoticCluster())
3108     {
3109       return !( GetCaloUtils()->GetEMCALRecoUtils()->IsExoticCell(absIdMax,cells,(GetReader()->GetInputEvent())->GetBunchCrossNumber()) );
3110     }
3111     else
3112     {
3113       return kTRUE;
3114     }
3115   }
3116   else // PHOS
3117   {
3118     Float_t ampMax = cells->GetCellAmplitude(absIdMax);
3119     GetCaloUtils()->RecalibrateCellAmplitude(ampMax, fCalorimeter, absIdMax);
3120
3121     if(ampMax < 0.01) return kFALSE;
3122
3123     if(1-GetECross(absIdMax,cells)/ampMax > 0.95) return kFALSE;
3124     else                                          return kTRUE;
3125   }
3126
3127 }
3128
3129 //_________________________________________________________
3130 void AliAnaCalorimeterQA::Print(const Option_t * opt) const
3131 {
3132   //Print some relevant parameters set for the analysis
3133   if(! opt)
3134     return;
3135
3136   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
3137   AliAnaCaloTrackCorrBaseClass::Print(" ");
3138
3139   printf("Select Calorimeter %s \n",fCalorimeter.Data());
3140   printf("Time Cut: %3.1f < TOF  < %3.1f\n", fTimeCutMin, fTimeCutMax);
3141   printf("EMCAL Min Amplitude   : %2.1f GeV/c\n", fEMCALCellAmpMin) ;
3142   printf("PHOS Min Amplitude    : %2.1f GeV/c\n", fPHOSCellAmpMin) ;
3143
3144 }
3145
3146 //_____________________________________________________
3147 void  AliAnaCalorimeterQA::MakeAnalysisFillHistograms()
3148 {
3149   //Fill Calorimeter QA histograms
3150
3151   //Play with the MC stack if available
3152   if(IsDataMC()) MCHistograms();
3153
3154   //Get List with CaloClusters
3155   TObjArray * caloClusters = NULL;
3156   if      (fCalorimeter == "PHOS")  caloClusters = GetPHOSClusters();
3157   else if (fCalorimeter == "EMCAL") caloClusters = GetEMCALClusters();
3158   else
3159     AliFatal(Form("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - Wrong calorimeter name <%s>, END\n", fCalorimeter.Data()));
3160
3161   // Do not do anything if there are no clusters
3162   if(caloClusters->GetEntriesFast() == 0) return;
3163
3164   //Get List with CaloCells
3165   AliVCaloCells * cells = 0x0;
3166   if(fCalorimeter == "PHOS") cells =  GetPHOSCells();
3167   else                                   cells =  GetEMCALCells();
3168
3169   if(!caloClusters || !cells) {
3170     AliFatal(Form("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - No CaloClusters or CaloCells available\n"));
3171     return; // trick coverity
3172   }
3173
3174   //printf("QA: N cells %d, N clusters %d \n",cells->GetNumberOfCells(),caloClusters->GetEntriesFast());
3175
3176   // Correlate Calorimeters and V0 and track Multiplicity
3177   if(fCorrelate)        Correlate();
3178
3179   // Clusters
3180   ClusterLoopHistograms(caloClusters,cells);
3181
3182   // Cells
3183   CellHistograms(cells);
3184
3185   if(GetDebug() > 0)
3186     printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - End \n");
3187
3188 }
3189
3190 //______________________________________
3191 void AliAnaCalorimeterQA::MCHistograms()
3192 {
3193   //Get the MC arrays and do some checks before filling MC histograms
3194
3195   TLorentzVector mom  ;
3196
3197   if(GetReader()->ReadStack()){
3198
3199     if(!GetMCStack())
3200       AliFatal("Stack not available, is the MC handler called?\n");
3201
3202     //Fill some pure MC histograms, only primaries.
3203     for(Int_t i=0 ; i<GetMCStack()->GetNprimary(); i++)
3204     {//Only primary particles, for all MC transport put GetNtrack()
3205       TParticle *primary = GetMCStack()->Particle(i) ;
3206
3207       if (primary->GetStatusCode() > 11) continue; //Working for PYTHIA and simple generators, check for HERWIG
3208       primary->Momentum(mom);
3209       MCHistograms(mom,TMath::Abs(primary->GetPdgCode()));
3210     } //primary loop
3211   }
3212   else if(GetReader()->ReadAODMCParticles()){
3213
3214     if(!GetReader()->GetAODMCParticles())
3215       AliFatal("AODMCParticles not available!");
3216
3217     //Fill some pure MC histograms, only primaries.
3218     for(Int_t i=0 ; i < (GetReader()->GetAODMCParticles())->GetEntriesFast(); i++)
3219     {
3220       AliAODMCParticle *aodprimary = (AliAODMCParticle*) (GetReader()->GetAODMCParticles())->At(i) ;
3221
3222       if (!aodprimary->IsPrimary()) continue; //accept all which is not MC transport generated. Don't know how to avoid partons
3223
3224       mom.SetPxPyPzE(aodprimary->Px(),aodprimary->Py(),aodprimary->Pz(),aodprimary->E());
3225       MCHistograms(mom,TMath::Abs(aodprimary->GetPdgCode()));
3226     } //primary loop
3227
3228   }
3229 }
3230
3231 //_______________________________________________________________________________
3232 void AliAnaCalorimeterQA::MCHistograms(TLorentzVector mom, Int_t pdg)
3233 {
3234   //Fill pure monte carlo related histograms
3235
3236   Float_t eMC    = mom.E();
3237   Float_t phiMC  = mom.Phi();
3238   if(phiMC < 0)
3239     phiMC  += TMath::TwoPi();
3240   Float_t etaMC  = mom.Eta();
3241
3242   if (TMath::Abs(etaMC) > 1) return;
3243
3244   Bool_t in = kFALSE;
3245
3246   //Rough stimate of acceptance for pi0, Eta and electrons
3247   if(fCalorimeter == "PHOS")
3248   {
3249     if(GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter))
3250       in = kTRUE ;
3251     if(GetDebug() > 2) printf("AliAnaCalorimeterQA::MCHistograms() - In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),in);
3252
3253   }
3254   else if(fCalorimeter == "EMCAL" && GetCaloUtils()->IsEMCALGeoMatrixSet())
3255   {
3256     if(GetEMCALGeometry())
3257     {
3258       Int_t absID=0;
3259       GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(mom.Eta(),mom.Phi(),absID);
3260
3261       if( absID >= 0)
3262         in = kTRUE;
3263
3264       if(GetDebug() > 2) printf("AliAnaCalorimeterQA::MCHistograms() - In %s Real acceptance? %d\n",fCalorimeter.Data(),in);
3265     }
3266     else
3267     {
3268       if(GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter))
3269         in = kTRUE ;
3270       if(GetDebug() > 2) printf("AliAnaCalorimeterQA::MCHistograms() - In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),in);
3271     }
3272   }
3273
3274   if (pdg==22)
3275   {
3276     fhGenMCE[kmcPhoton]    ->Fill(eMC);
3277     if(eMC > 0.5) fhGenMCEtaPhi[kmcPhoton]->Fill(etaMC,phiMC);
3278     if(in)
3279     {
3280       fhGenMCAccE[kmcPhoton]    ->Fill(eMC);
3281       if(eMC > 0.5) fhGenMCAccEtaPhi[kmcPhoton]->Fill(etaMC,phiMC);
3282     }
3283   }
3284   else if (pdg==111)
3285   {
3286     fhGenMCE[kmcPi0]    ->Fill(eMC);
3287     if(eMC > 0.5) fhGenMCEtaPhi[kmcPi0]->Fill(etaMC,phiMC);
3288     if(in)
3289     {
3290       fhGenMCAccE[kmcPi0]    ->Fill(eMC);
3291       if(eMC > 0.5) fhGenMCAccEtaPhi[kmcPi0]->Fill(etaMC,phiMC);
3292     }
3293   }
3294   else if (pdg==221)
3295   {
3296     fhGenMCE[kmcEta]    ->Fill(eMC);
3297     if(eMC > 0.5) fhGenMCEtaPhi[kmcEta]->Fill(etaMC,phiMC);
3298     if(in)
3299     {
3300       fhGenMCAccE[kmcEta]    ->Fill(eMC);
3301       if(eMC > 0.5) fhGenMCAccEtaPhi[kmcEta]->Fill(etaMC,phiMC);
3302     }
3303   }
3304   else if (TMath::Abs(pdg)==11)
3305   {
3306     fhGenMCE[kmcElectron]    ->Fill(eMC);
3307     if(eMC > 0.5) fhGenMCEtaPhi[kmcElectron]->Fill(etaMC,phiMC);
3308     if(in)
3309     {
3310       fhGenMCAccE[kmcElectron]    ->Fill(eMC);
3311       if(eMC > 0.5) fhGenMCAccEtaPhi[kmcElectron]->Fill(etaMC,phiMC);
3312     }
3313   }
3314 }
3315
3316 //_________________________________________________________________________________
3317 void AliAnaCalorimeterQA::WeightHistograms(AliVCluster *clus, AliVCaloCells* cells)
3318 {
3319   // Calculate weights
3320
3321   // First recalculate energy in case non linearity was applied
3322   Float_t  energy = 0;
3323   Float_t  ampMax = 0;
3324   for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++)
3325   {
3326     Int_t id       = clus->GetCellsAbsId()[ipos];
3327
3328     //Recalibrate cell energy if needed
3329     Float_t amp = cells->GetCellAmplitude(id);
3330     GetCaloUtils()->RecalibrateCellAmplitude(amp, fCalorimeter, id);
3331
3332     energy    += amp;
3333
3334     if(amp> ampMax)
3335       ampMax = amp;
3336
3337   } // energy loop
3338
3339   if(energy <=0 )
3340   {
3341     printf("AliAnaCalorimeterQA::WeightHistograms()- Wrong calculated energy %f\n",energy);
3342     return;
3343   }
3344
3345   fhEMaxCellClusterRatio   ->Fill(energy,ampMax/energy);
3346   fhEMaxCellClusterLogRatio->Fill(energy,TMath::Log(ampMax/energy));
3347
3348   //Get the ratio and log ratio to all cells in cluster
3349   for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++)
3350   {
3351     Int_t id       = clus->GetCellsAbsId()[ipos];
3352
3353     //Recalibrate cell energy if needed
3354     Float_t amp = cells->GetCellAmplitude(id);
3355     GetCaloUtils()->RecalibrateCellAmplitude(amp, fCalorimeter, id);
3356
3357     fhECellClusterRatio   ->Fill(energy,amp/energy);
3358     fhECellClusterLogRatio->Fill(energy,TMath::Log(amp/energy));
3359   }
3360
3361   //Recalculate shower shape for different W0
3362   if(fCalorimeter=="EMCAL")
3363   {
3364     Float_t l0org = clus->GetM02();
3365     Float_t l1org = clus->GetM20();
3366     Float_t dorg  = clus->GetDispersion();
3367
3368     for(Int_t iw = 0; iw < 14; iw++){
3369       GetCaloUtils()->GetEMCALRecoUtils()->SetW0(1+iw*0.5);
3370       GetCaloUtils()->GetEMCALRecoUtils()->RecalculateClusterShowerShapeParameters(GetEMCALGeometry(), cells, clus);
3371
3372       fhLambda0ForW0[iw]->Fill(energy,clus->GetM02());
3373       //fhLambda1ForW0[iw]->Fill(energy,clus->GetM20());
3374
3375       if(IsDataMC()){
3376
3377         Int_t tag = GetMCAnalysisUtils()->CheckOrigin(clus->GetLabels(),clus->GetNLabels(), GetReader());
3378
3379         if(   GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPhoton) &&
3380            !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPi0)      &&
3381            !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCEta)      &&
3382            !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion)        ){
3383           fhLambda0ForW0MC[iw][0]->Fill(energy,clus->GetM02());
3384           //fhLambda1ForW0MC[iw][0]->Fill(energy,clus->GetM20());
3385         } // Pure Photon
3386         else if( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCElectron) &&
3387                 !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion)        ){
3388           fhLambda0ForW0MC[iw][1]->Fill(energy,clus->GetM02());
3389           //fhLambda1ForW0MC[iw][1]->Fill(energy,clus->GetM20());
3390         } // Electron
3391         else if( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion)        ){
3392           fhLambda0ForW0MC[iw][2]->Fill(energy,clus->GetM02());
3393           //fhLambda1ForW0MC[iw][2]->Fill(energy,clus->GetM20());
3394         } // Conversion
3395         else if( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPi0) ){
3396           fhLambda0ForW0MC[iw][3]->Fill(energy,clus->GetM02());
3397           //fhLambda1ForW0MC[iw][3]->Fill(energy,clus->GetM20());
3398         }// Pi0
3399         else if(!GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCEta) &&
3400                 !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPhoton) ){
3401           fhLambda0ForW0MC[iw][4]->Fill(energy,clus->GetM02());
3402           //fhLambda1ForW0MC[iw][4]->Fill(energy,clus->GetM20());
3403         }// Hadron
3404
3405       }// Is MC
3406     } // w0 loop
3407
3408     // Set the original values back
3409     clus->SetM02(l0org);
3410     clus->SetM20(l1org);
3411     clus->SetDispersion(dorg);
3412
3413   }// EMCAL
3414
3415 }
3416
3417
3418