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