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