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