PWG4/PartCorrDep/AliAnaPi0.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 collect two-photon invariant mass distributions for
  19 // extracting raw pi0 yield.
  20 // Input is produced by AliAnaPhoton (or any other analysis producing output AliAODPWG4Particles),
  21 // it will do nothing if executed alone
  22 //
  23 //-- Author: Dmitri Peressounko (RRC "KI")
  24 //-- Adapted to PartCorr frame by Lamia Benhabib (SUBATECH)
  25 //-- and Gustavo Conesa (INFN-Frascati)
  26 //_________________________________________________________________________
  27
  28
  29 // --- ROOT system ---
  30 #include "TH3.h"
  31 #include "TH2D.h"
  32 //#include "Riostream.h"
  33 #include "TCanvas.h"
  34 #include "TPad.h"
  35 #include "TROOT.h"
  36 #include "TClonesArray.h"
  37 #include "TObjString.h"
  38 #include "TDatabasePDG.h"
  39
  40 //---- AliRoot system ----
  41 #include "AliAnaPi0.h"
  42 #include "AliCaloTrackReader.h"
  43 #include "AliCaloPID.h"
  44 #include "AliStack.h"
  45 #include "AliFiducialCut.h"
  46 #include "TParticle.h"
  47 #include "AliVEvent.h"
  48 #include "AliESDCaloCluster.h"
  49 #include "AliESDEvent.h"
  50 #include "AliAODEvent.h"
  51 #include "AliNeutralMesonSelection.h"
  52 #include "AliMixedEvent.h"
  53 #include "AliAODMCParticle.h"
  54
  55 ClassImp(AliAnaPi0)
  56
  57 //________________________________________________________________________________________________________________________________________________
  58 AliAnaPi0::AliAnaPi0() : AliAnaPartCorrBaseClass(),
  59 fDoOwnMix(kFALSE),fNCentrBin(0),//fNZvertBin(0),fNrpBin(0),
  60 fNmaxMixEv(0), fCalorimeter(""),
  61 fNModules(12), fUseAngleCut(kFALSE), fUseAngleEDepCut(kFALSE),fAngleCut(0), fAngleMaxCut(7.),fEventsList(0x0), fMultiCutAna(kFALSE), fMultiCutAnaSim(kFALSE),
  62 fNPtCuts(0),fNAsymCuts(0), fNCellNCuts(0),fNPIDBits(0),  fMakeInvPtPlots(kFALSE), fSameSM(kFALSE),
  63 fUseTrackMultBins(kFALSE),fUsePhotonMultBins(kFALSE),fUseAverCellEBins(kFALSE), fUseAverClusterEBins(kFALSE),
  64 fUseAverClusterEDenBins(0), //fUseAverClusterPairRBins(0), fUseAverClusterPairRWeightBins(0), fUseEMaxBins(0),
  65 fFillBadDistHisto(kFALSE),
  66 fhAverTotECluster(0), fhAverTotECell(0), fhAverTotECellvsCluster(0),
  67 fhEDensityCluster(0), fhEDensityCell(0), fhEDensityCellvsCluster(0),
  68 //fhClusterPairDist(0), fhClusterPairDistWeight(0),  fhAverClusterPairDist(0), fhAverClusterPairDistWeight(0),
  69 //fhAverClusterPairDistvsAverE(0), fhAverClusterPairDistWeightvsAverE(0),fhAverClusterPairDistvsN(0), fhAverClusterPairDistWeightvsN(0),
  70 //fhMaxEvsClustMult(0), fhMaxEvsClustEDen(0),
  71 fhReMod(0x0),    fhReSameSideEMCALMod(0x0), fhReSameSectorEMCALMod(0x0), fhReDiffPHOSMod(0x0),
  72 fhMiMod(0x0),    fhMiSameSideEMCALMod(0x0), fhMiSameSectorEMCALMod(0x0), fhMiDiffPHOSMod(0x0),
  73 fhReConv(0x0),   fhMiConv(0x0),    fhReConv2(0x0),  fhMiConv2(0x0),
  74 fhRe1(0x0),      fhMi1(0x0),       fhRe2(0x0),      fhMi2(0x0),      fhRe3(0x0),      fhMi3(0x0),
  75 fhReInvPt1(0x0), fhMiInvPt1(0x0),  fhReInvPt2(0x0), fhMiInvPt2(0x0), fhReInvPt3(0x0), fhMiInvPt3(0x0),
  76 fhRePtNCellAsymCuts(0x0), fhRePtNCellAsymCutsSM0(0x0), fhRePtNCellAsymCutsSM1(0x0), fhRePtNCellAsymCutsSM2(0x0), fhRePtNCellAsymCutsSM3(0x0), fhMiPtNCellAsymCuts(0x0),
  77 fhRePIDBits(0x0),fhRePtMult(0x0), fhRePtAsym(0x0), fhRePtAsymPi0(0x0),fhRePtAsymEta(0x0),
  78 fhEvents(0x0),   fhCentrality(0x0),
  79 fhRealOpeningAngle(0x0),fhRealCosOpeningAngle(0x0), fhMixedOpeningAngle(0x0),fhMixedCosOpeningAngle(0x0),
  80 fhPrimPi0Pt(0x0), fhPrimPi0AccPt(0x0), fhPrimPi0Y(0x0), fhPrimPi0AccY(0x0), fhPrimPi0Phi(0x0), fhPrimPi0AccPhi(0x0),
  81 fhPrimPi0OpeningAngle(0x0), fhPrimPi0CosOpeningAngle(0x0),
  82 fhPrimEtaPt(0x0), fhPrimEtaAccPt(0x0), fhPrimEtaY(0x0), fhPrimEtaAccY(0x0), fhPrimEtaPhi(0x0), fhPrimEtaAccPhi(0x0),
  83 fhPrimPi0PtOrigin(0x0),  fhPrimEtaPtOrigin(0x0),
  84 fhMCOrgMass(),fhMCOrgAsym(),  fhMCOrgDeltaEta(),fhMCOrgDeltaPhi(),
  85 fhMCPi0MassPtRec(), fhMCPi0MassPtTrue(), fhMCPi0PtTruePtRec(), fhMCEtaMassPtRec(), fhMCEtaMassPtTrue(), fhMCEtaPtTruePtRec(),
  86 fhMCPi0PtOrigin(0x0),  fhMCEtaPtOrigin(0x0)
  87 {
  88 //Default Ctor
  89  InitParameters();
  90
  91 }
  92
  93 //________________________________________________________________________________________________________________________________________________
  94 AliAnaPi0::~AliAnaPi0() {
  95   // Remove event containers
  96
  97   if(fDoOwnMix && fEventsList){
  98     for(Int_t ic=0; ic<fNCentrBin; ic++){
  99       for(Int_t iz=0; iz<GetNZvertBin(); iz++){
 100         for(Int_t irp=0; irp<GetNRPBin(); irp++){
 101           fEventsList[ic*GetNZvertBin()*GetNRPBin()+iz*GetNRPBin()+irp]->Delete() ;
 102           delete fEventsList[ic*GetNZvertBin()*GetNRPBin()+iz*GetNRPBin()+irp] ;
 103         }
 104       }
 105     }
 106     delete[] fEventsList;
 107     fEventsList=0 ;
 108   }
 109
 110 }
 111
 112 //________________________________________________________________________________________________________________________________________________
 113 void AliAnaPi0::InitParameters()
 114 {
 115 //Init parameters when first called the analysis
 116 //Set default parameters
 117   SetInputAODName("PWG4Particle");
 118
 119   AddToHistogramsName("AnaPi0_");
 120   fNModules = 12; // set maximum to maximum number of EMCAL modules
 121   fNCentrBin = 1;
 122 //  fNZvertBin = 1;
 123 //  fNrpBin    = 1;
 124   fNmaxMixEv = 10;
 125
 126   fCalorimeter  = "PHOS";
 127   fUseAngleCut = kFALSE;
 128   fUseAngleEDepCut = kFALSE;
 129   fAngleCut    = 0.;
 130   fAngleMaxCut = TMath::Pi();
 131
 132   fMultiCutAna = kFALSE;
 133
 134   fNPtCuts = 3;
 135   fPtCuts[0] = 0.; fPtCuts[1] = 0.3;   fPtCuts[2] = 0.5;
 136   for(Int_t i = fNPtCuts; i < 10; i++)fPtCuts[i] = 0.;
 137
 138   fNAsymCuts = 4;
 139   fAsymCuts[0] = 1.;  fAsymCuts[1] = 0.8; fAsymCuts[2] = 0.6;   fAsymCuts[3] = 0.1;
 140   for(Int_t i = fNAsymCuts; i < 10; i++)fAsymCuts[i] = 0.;
 141
 142   fNCellNCuts = 3;
 143   fCellNCuts[0] = 0; fCellNCuts[1] = 1;   fCellNCuts[2] = 2;
 144   for(Int_t i = fNCellNCuts; i < 10; i++)fCellNCuts[i]  = 0;
 145
 146   fNPIDBits = 2;
 147   fPIDBits[0] = 0;   fPIDBits[1] = 2; //  fPIDBits[2] = 4; fPIDBits[3] = 6;// check, no cut,  dispersion, neutral, dispersion&&neutral
 148   for(Int_t i = fNPIDBits; i < 10; i++)fPIDBits[i] = 0;
 149
 150 }
 151
 152
 153 //________________________________________________________________________________________________________________________________________________
 154 TObjString * AliAnaPi0::GetAnalysisCuts()
 155 {
 156   //Save parameters used for analysis
 157   TString parList ; //this will be list of parameters used for this analysis.
 158   const Int_t buffersize = 255;
 159   char onePar[buffersize] ;
 160   snprintf(onePar,buffersize,"--- AliAnaPi0 ---\n") ;
 161   parList+=onePar ;
 162   snprintf(onePar,buffersize,"Number of bins in Centrality:  %d \n",fNCentrBin) ;
 163   parList+=onePar ;
 164   snprintf(onePar,buffersize,"Number of bins in Z vert. pos: %d \n",GetNZvertBin()) ;
 165   parList+=onePar ;
 166   snprintf(onePar,buffersize,"Number of bins in Reac. Plain: %d \n",GetNRPBin()) ;
 167   parList+=onePar ;
 168   snprintf(onePar,buffersize,"Depth of event buffer: %d \n",fNmaxMixEv) ;
 169   parList+=onePar ;
 170   snprintf(onePar,buffersize,"Pair in same Module: %d ; TrackMult as centrality: %d; PhotonMult as centrality: %d; cluster E as centrality: %d; cell as centrality: %d; Fill InvPt histos %d\n",
 171            fSameSM, fUseTrackMultBins, fUsePhotonMultBins, fUseAverClusterEBins, fUseAverCellEBins, fMakeInvPtPlots) ;
 172   parList+=onePar ;
 173   snprintf(onePar,buffersize,"Select pairs with their angle: %d, edep %d, min angle %2.3f, max angle %2.3f,\n",fUseAngleCut, fUseAngleEDepCut,fAngleCut,fAngleMaxCut) ;
 174   parList+=onePar ;
 175   snprintf(onePar,buffersize," Asymmetry cuts: n = %d, asymmetry < ",fNAsymCuts) ;
 176   for(Int_t i = 0; i < fNAsymCuts; i++) snprintf(onePar,buffersize,"%s %2.2f;",onePar,fAsymCuts[i]);
 177   parList+=onePar ;
 178   snprintf(onePar,buffersize," PID selection bits: n = %d, PID bit =\n",fNPIDBits) ;
 179   for(Int_t i = 0; i < fNPIDBits; i++) snprintf(onePar,buffersize,"%s %d;",onePar,fPIDBits[i]);
 180   parList+=onePar ;
 181   snprintf(onePar,buffersize,"Cuts: \n") ;
 182   parList+=onePar ;
 183   snprintf(onePar,buffersize,"Z vertex position: -%f < z < %f \n",GetZvertexCut(),GetZvertexCut()) ;
 184   parList+=onePar ;
 185   snprintf(onePar,buffersize,"Calorimeter: %s \n",fCalorimeter.Data()) ;
 186   parList+=onePar ;
 187   snprintf(onePar,buffersize,"Number of modules: %d \n",fNModules) ;
 188   parList+=onePar ;
 189   if(fMultiCutAna){
 190     snprintf(onePar, buffersize," pT cuts: n = %d, pt > ",fNPtCuts) ;
 191     for(Int_t i = 0; i < fNPtCuts; i++) snprintf(onePar,buffersize,"%s %2.2f;",onePar,fPtCuts[i]);
 192     parList+=onePar ;
 193     snprintf(onePar,buffersize, " N cell in cluster cuts: n = %d, nCell > ",fNCellNCuts) ;
 194     for(Int_t i = 0; i < fNCellNCuts; i++) snprintf(onePar,buffersize,"%s %d;",onePar,fCellNCuts[i]);
 195     parList+=onePar ;
 196   }
 197
 198   return new TObjString(parList) ;
 199 }
 200
 201 //________________________________________________________________________________________________________________________________________________
 202 TList * AliAnaPi0::GetCreateOutputObjects()
 203 {
 204   // Create histograms to be saved in output file and
 205   // store them in fOutputContainer
 206
 207   //create event containers
 208   fEventsList = new TList*[fNCentrBin*GetNZvertBin()*GetNRPBin()] ;
 209
 210   for(Int_t ic=0; ic<fNCentrBin; ic++){
 211     for(Int_t iz=0; iz<GetNZvertBin(); iz++){
 212       for(Int_t irp=0; irp<GetNRPBin(); irp++){
 213         fEventsList[ic*GetNZvertBin()*GetNRPBin()+iz*GetNRPBin()+irp] = new TList() ;
 214         fEventsList[ic*GetNZvertBin()*GetNRPBin()+iz*GetNRPBin()+irp]->SetOwner(kFALSE);
 215       }
 216     }
 217   }
 218
 219   TList * outputContainer = new TList() ;
 220   outputContainer->SetName(GetName());
 221
 222   fhReMod                = new TH2D*[fNModules]   ;
 223   fhMiMod                = new TH2D*[fNModules]   ;
 224
 225   if(fCalorimeter == "PHOS"){
 226     fhReDiffPHOSMod        = new TH2D*[fNModules]   ;
 227     fhMiDiffPHOSMod        = new TH2D*[fNModules]   ;
 228   }
 229   else{
 230     fhReSameSectorEMCALMod = new TH2D*[fNModules/2] ;
 231     fhReSameSideEMCALMod   = new TH2D*[fNModules-2] ;
 232     fhMiSameSectorEMCALMod = new TH2D*[fNModules/2] ;
 233     fhMiSameSideEMCALMod   = new TH2D*[fNModules-2] ;
 234   }
 235
 236
 237   fhRe1 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
 238   fhMi1 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
 239   if(fFillBadDistHisto){
 240     fhRe2 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
 241     fhRe3 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
 242     fhMi2 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
 243     fhMi3 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
 244   }
 245   if(fMakeInvPtPlots) {
 246     fhReInvPt1 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
 247     fhMiInvPt1 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
 248     if(fFillBadDistHisto){
 249       fhReInvPt2 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
 250       fhReInvPt3 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
 251       fhMiInvPt2 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
 252       fhMiInvPt3 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
 253     }
 254   }
 255
 256   const Int_t buffersize = 255;
 257   char key[buffersize] ;
 258   char title[buffersize] ;
 259
 260   Int_t nptbins   = GetHistoPtBins();
 261   Int_t nphibins  = GetHistoPhiBins();
 262   Int_t netabins  = GetHistoEtaBins();
 263   Float_t ptmax   = GetHistoPtMax();
 264   Float_t phimax  = GetHistoPhiMax();
 265   Float_t etamax  = GetHistoEtaMax();
 266   Float_t ptmin   = GetHistoPtMin();
 267   Float_t phimin  = GetHistoPhiMin();
 268   Float_t etamin  = GetHistoEtaMin();
 269
 270   Int_t nmassbins = GetHistoMassBins();
 271   Int_t nasymbins = GetHistoAsymmetryBins();
 272   Float_t massmax = GetHistoMassMax();
 273   Float_t asymmax = GetHistoAsymmetryMax();
 274   Float_t massmin = GetHistoMassMin();
 275   Float_t asymmin = GetHistoAsymmetryMin();
 276   Int_t ntrmbins  = GetHistoTrackMultiplicityBins();
 277   Int_t ntrmmax   = GetHistoTrackMultiplicityMax();
 278   Int_t ntrmmin   = GetHistoTrackMultiplicityMin();
 279
 280   fhAverTotECluster = new TH1F("hAverTotECluster","hAverTotECluster",200,0,50) ;
 281   fhAverTotECluster->SetXTitle("E_{cluster, aver. SM} (GeV)");
 282   outputContainer->Add(fhAverTotECluster) ;
 283
 284   fhAverTotECell    = new TH1F("hAverTotECell","hAverTotECell",200,0,50) ;
 285   fhAverTotECell->SetXTitle("E_{cell, aver. SM} (GeV)");
 286   outputContainer->Add(fhAverTotECell) ;
 287
 288   fhAverTotECellvsCluster    = new TH2F("hAverTotECellvsCluster","hAverTotECellvsCluster",200,0,50,200,0,50) ;
 289   fhAverTotECellvsCluster->SetYTitle("E_{cell, aver. SM} (GeV)");
 290   fhAverTotECellvsCluster->SetXTitle("E_{cluster, aver. SM} (GeV)");
 291   outputContainer->Add(fhAverTotECellvsCluster) ;
 292
 293   fhEDensityCluster = new TH1F("hEDensityCluster","hEDensityCluster",200,0,50) ;
 294   fhEDensityCluster->SetXTitle("#Sigma E_{cluster} / N_{cluster} (GeV)");
 295   outputContainer->Add(fhEDensityCluster) ;
 296
 297   fhEDensityCell    = new TH1F("hEDensityCell","hEDensityCell",200,0,50) ;
 298   fhEDensityCell->SetXTitle("#Sigma E_{cell} / N_{cell} (GeV)");
 299   outputContainer->Add(fhEDensityCell) ;
 300
 301   fhEDensityCellvsCluster    = new TH2F("hEDensityCellvsCluster","hEDensityCellvsCluster",200,0,50,200,0,50) ;
 302   fhEDensityCellvsCluster->SetYTitle("#Sigma E_{cell} / N_{cell} (GeV)");
 303   fhEDensityCellvsCluster->SetXTitle("#Sigma E_{cluster} / N_{cluster} (GeV)");
 304   outputContainer->Add(fhEDensityCellvsCluster) ;
 305
 306 //  fhClusterPairDist = new TH1F("hClusterPairDist","Distance between clusters",250,0,750) ;
 307 //  fhClusterPairDist->SetXTitle("#sqrt{(x_{1}-x_{2})^2+(z_{1}-z_{2})^2} (cm)");
 308 //  outputContainer->Add(fhClusterPairDist) ;
 309 //
 310 //  fhClusterPairDistWeight = new TH1F("hClusterPairDistWeighted","Distance between clusters, weighted by pair energy",200,0,400) ;
 311 //  fhClusterPairDistWeight->SetXTitle("#sqrt{(x_{1}E_{1}-x_{2}E_{2})^{2}+(z_{1}E_{1}-z_{2}E_{2})^{2}}/ (E_{1}+E_{2}) (cm)");
 312 //  outputContainer->Add(fhClusterPairDistWeight) ;
 313 //
 314 //  fhAverClusterPairDist = new TH1F("hAverClusterPairDist","Average distance between clusters",250,0,750) ;
 315 //  fhAverClusterPairDist->SetXTitle("#Sigma (#sqrt{(x_{1}-x_{2})^{2}+(z_{1}-z_{2})^{2}}) / N_{pairs} (cm)");
 316 //  outputContainer->Add(fhAverClusterPairDist) ;
 317 //
 318 //  fhAverClusterPairDistWeight = new TH1F("hAverClusterPairDistWeighted","Average distance between clusters, weighted by pair energy",200,0,400) ;
 319 //  fhAverClusterPairDistWeight->SetXTitle("#Sigma (#sqrt{(x_{1}E_{1}-x_{2}E_{2})^{2}+(z_{1}E_{1}-z_{2}E_{2})^{2}}/ (E_{1}+E_{2})) / N_{pairs} (cm)");
 320 //  outputContainer->Add(fhAverClusterPairDistWeight) ;
 321 //
 322 //  fhAverClusterPairDistvsAverE = new TH2F("hAverClusterPairDistvsAverE","Average distance between clusters",250,0,750,200,0,50) ;
 323 //  fhAverClusterPairDistvsAverE->SetXTitle("#Sigma (#sqrt{(x_{1}-x_{2})^{2}+(z_{1}-z_{2})^{2}}) / N_{pairs} (cm)");
 324 //  fhAverClusterPairDistvsAverE->SetYTitle("#Sigma E_{cluster} / N_{cluster} (GeV)");
 325 //  outputContainer->Add(fhAverClusterPairDistvsAverE) ;
 326 //
 327 //  fhAverClusterPairDistWeightvsAverE = new TH2F("hAverClusterPairDistWeightedvsAverE","Average distance between clusters, weighted by pair energy",200,0,400,200,0,50) ;
 328 //  fhAverClusterPairDistWeightvsAverE->SetXTitle("#Sigma (#sqrt{(x_{1}E_{1}-x_{2}E_{2})^2+(z_{1}E_{1}-z_{2}E_{2})^2}/ (E_{1}+E_{2})) / N_{pairs} (cm/GeV)");
 329 //  fhAverClusterPairDistWeightvsAverE->SetYTitle("#Sigma E_{cluster} / N_{cluster} (GeV)");
 330 //  outputContainer->Add(fhAverClusterPairDistWeightvsAverE) ;
 331
 332 //  fhAverClusterPairDistvsN = new TH2F("hAverClusterPairDistvsN","Average distance between clusters",250,0,750,200,0,50) ;
 333 //  fhAverClusterPairDistvsN->SetXTitle("#Sigma (#sqrt{(x_{1}-x_{2})^{2}+(z_{1}-z_{2})^{2}}) / N_{pairs} (cm)");
 334 //  fhAverClusterPairDistvsN->SetYTitle("N_{cluster}");
 335 //  outputContainer->Add(fhAverClusterPairDistvsN) ;
 336 //
 337 //  fhAverClusterPairDistWeightvsN = new TH2F("hAverClusterPairDistWeightedvsN","Average distance between clusters, weighted by pair energy",200,0,400,200,0,50) ;
 338 //  fhAverClusterPairDistWeightvsN->SetXTitle("#Sigma (#sqrt{(x_{1}E_{1}-x_{2}E_{2})^{2}+(z_{1}E_{1}-z_{2}E_{2})^{2}}/ (E_{1}+E_{2})) / N_{pairs} (cm)");
 339 //  fhAverClusterPairDistWeightvsN->SetYTitle("N_{cluster}");
 340 //  outputContainer->Add(fhAverClusterPairDistWeightvsN) ;
 341
 342 //  fhMaxEvsClustMult = new TH2F("hMaxEvsClustMult","",nptbins,ptmin,ptmax,50,0,50) ;
 343 //  fhMaxEvsClustMult->SetXTitle("E_{max}");
 344 //  fhMaxEvsClustMult->SetYTitle("N_{cluster}");
 345 //  outputContainer->Add(fhMaxEvsClustMult) ;
 346 //
 347 //  fhMaxEvsClustEDen = new TH2F("hMaxEvsClustEDen","",nptbins,ptmin,ptmax,200,0,50) ;
 348 //  fhMaxEvsClustEDen->SetXTitle("E_{max}");
 349 //  fhMaxEvsClustEDen->SetYTitle("#Sigma E_{cluster} / N_{cluster} (GeV)");
 350 //  outputContainer->Add(fhMaxEvsClustEDen) ;
 351
 352   fhReConv = new TH2D("hReConv","Real Pair with one recombined conversion ",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 353   fhReConv->SetXTitle("p_{T} (GeV/c)");
 354   fhReConv->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 355   outputContainer->Add(fhReConv) ;
 356
 357   fhReConv2 = new TH2D("hReConv2","Real Pair with 2 recombined conversion ",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 358   fhReConv2->SetXTitle("p_{T} (GeV/c)");
 359   fhReConv2->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 360   outputContainer->Add(fhReConv2) ;
 361
 362   if(fDoOwnMix){
 363     fhMiConv = new TH2D("hMiConv","Mixed Pair with one recombined conversion ",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 364     fhMiConv->SetXTitle("p_{T} (GeV/c)");
 365     fhMiConv->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 366     outputContainer->Add(fhMiConv) ;
 367
 368     fhMiConv2 = new TH2D("hMiConv2","Mixed Pair with 2 recombined conversion ",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 369     fhMiConv2->SetXTitle("p_{T} (GeV/c)");
 370     fhMiConv2->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 371     outputContainer->Add(fhMiConv2) ;
 372   }
 373
 374   for(Int_t ic=0; ic<fNCentrBin; ic++){
 375       for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
 376         for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
 377           Int_t index = ((ic*fNPIDBits)+ipid)*fNAsymCuts + iasym;
 378           //printf("cen %d, pid %d, asy %d, Index %d\n",ic,ipid,iasym,index);
 379           //Distance to bad module 1
 380           snprintf(key, buffersize,"hRe_cen%d_pidbit%d_asy%d_dist1",ic,ipid,iasym) ;
 381           snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 1",
 382                    ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
 383           fhRe1[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 384           fhRe1[index]->SetXTitle("p_{T} (GeV/c)");
 385           fhRe1[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 386           //printf("name: %s\n ",fhRe1[index]->GetName());
 387           outputContainer->Add(fhRe1[index]) ;
 388
 389           if(fFillBadDistHisto){
 390             //Distance to bad module 2
 391             snprintf(key, buffersize,"hRe_cen%d_pidbit%d_asy%d_dist2",ic,ipid,iasym) ;
 392             snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 2",
 393                      ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
 394             fhRe2[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 395             fhRe2[index]->SetXTitle("p_{T} (GeV/c)");
 396             fhRe2[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 397             outputContainer->Add(fhRe2[index]) ;
 398
 399             //Distance to bad module 3
 400             snprintf(key, buffersize,"hRe_cen%d_pidbit%d_asy%d_dist3",ic,ipid,iasym) ;
 401             snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 3",
 402                      ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
 403             fhRe3[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 404             fhRe3[index]->SetXTitle("p_{T} (GeV/c)");
 405             fhRe3[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 406             outputContainer->Add(fhRe3[index]) ;
 407           }
 408
 409           //Inverse pT
 410           if(fMakeInvPtPlots){
 411             //Distance to bad module 1
 412             snprintf(key, buffersize,"hReInvPt_cen%d_pidbit%d_asy%d_dist1",ic,ipid,iasym) ;
 413             snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 1",
 414                      ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
 415             fhReInvPt1[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 416             fhReInvPt1[index]->SetXTitle("p_{T} (GeV/c)");
 417             fhReInvPt1[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 418             outputContainer->Add(fhReInvPt1[index]) ;
 419
 420             if(fFillBadDistHisto){
 421               //Distance to bad module 2
 422               snprintf(key, buffersize,"hReInvPt_cen%d_pidbit%d_asy%d_dist2",ic,ipid,iasym) ;
 423               snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 2",
 424                        ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
 425               fhReInvPt2[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 426               fhReInvPt2[index]->SetXTitle("p_{T} (GeV/c)");
 427               fhReInvPt2[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 428               outputContainer->Add(fhReInvPt2[index]) ;
 429
 430               //Distance to bad module 3
 431               snprintf(key, buffersize,"hReInvPt_cen%d_pidbit%d_asy%d_dist3",ic,ipid,iasym) ;
 432               snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 3",
 433                        ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
 434               fhReInvPt3[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 435               fhReInvPt3[index]->SetXTitle("p_{T} (GeV/c)");
 436               fhReInvPt3[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 437               outputContainer->Add(fhReInvPt3[index]) ;
 438             }
 439           }
 440           if(fDoOwnMix){
 441             //Distance to bad module 1
 442             snprintf(key, buffersize,"hMi_cen%d_pidbit%d_asy%d_dist1",ic,ipid,iasym) ;
 443             snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 1",
 444                      ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
 445             fhMi1[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 446             fhMi1[index]->SetXTitle("p_{T} (GeV/c)");
 447             fhMi1[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 448             outputContainer->Add(fhMi1[index]) ;
 449             if(fFillBadDistHisto){
 450               //Distance to bad module 2
 451               snprintf(key, buffersize,"hMi_cen%d_pidbit%d_asy%d_dist2",ic,ipid,iasym) ;
 452               snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 2",
 453                        ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
 454               fhMi2[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 455               fhMi2[index]->SetXTitle("p_{T} (GeV/c)");
 456               fhMi2[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 457               outputContainer->Add(fhMi2[index]) ;
 458
 459               //Distance to bad module 3
 460               snprintf(key, buffersize,"hMi_cen%d_pidbit%d_asy%d_dist3",ic,ipid,iasym) ;
 461               snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 3",
 462                        ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
 463               fhMi3[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 464               fhMi3[index]->SetXTitle("p_{T} (GeV/c)");
 465               fhMi3[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 466               outputContainer->Add(fhMi3[index]) ;
 467             }
 468             //Inverse pT
 469             if(fMakeInvPtPlots){
 470               //Distance to bad module 1
 471               snprintf(key, buffersize,"hMiInvPt_cen%d_pidbit%d_asy%d_dist1",ic,ipid,iasym) ;
 472               snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 1",
 473                        ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
 474               fhMiInvPt1[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 475               fhMiInvPt1[index]->SetXTitle("p_{T} (GeV/c)");
 476               fhMiInvPt1[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 477               outputContainer->Add(fhMiInvPt1[index]) ;
 478               if(fFillBadDistHisto){
 479                 //Distance to bad module 2
 480                 snprintf(key, buffersize,"hMiInvPt_cen%d_pidbit%d_asy%d_dist2",ic,ipid,iasym) ;
 481                 snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 2",
 482                          ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
 483                 fhMiInvPt2[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 484                 fhMiInvPt2[index]->SetXTitle("p_{T} (GeV/c)");
 485                 fhMiInvPt2[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 486                 outputContainer->Add(fhMiInvPt2[index]) ;
 487
 488                 //Distance to bad module 3
 489                 snprintf(key, buffersize,"hMiInvPt_cen%d_pidbit%d_asy%d_dist3",ic,ipid,iasym) ;
 490                 snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f,dist bad 3",
 491                          ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
 492                 fhMiInvPt3[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 493                 fhMiInvPt3[index]->SetXTitle("p_{T} (GeV/c)");
 494                 fhMiInvPt3[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 495                 outputContainer->Add(fhMiInvPt3[index]) ;
 496               }
 497             }
 498           }
 499         }
 500       }
 501     }
 502
 503   fhRePtAsym = new TH2D("hRePtAsym","Asymmetry vs pt, for pairs",nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax) ;
 504   fhRePtAsym->SetXTitle("p_{T} (GeV/c)");
 505   fhRePtAsym->SetYTitle("Asymmetry");
 506   outputContainer->Add(fhRePtAsym);
 507
 508   fhRePtAsymPi0 = new TH2D("hRePtAsymPi0","Asymmetry vs pt, for pairs close to #pi^{0} mass",nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax) ;
 509   fhRePtAsymPi0->SetXTitle("p_{T} (GeV/c)");
 510   fhRePtAsymPi0->SetYTitle("Asymmetry");
 511   outputContainer->Add(fhRePtAsymPi0);
 512
 513   fhRePtAsymEta = new TH2D("hRePtAsymEta","Asymmetry vs pt, for pairs close to #eta mass",nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax) ;
 514   fhRePtAsymEta->SetXTitle("p_{T} (GeV/c)");
 515   fhRePtAsymEta->SetYTitle("Asymmetry");
 516   outputContainer->Add(fhRePtAsymEta);
 517
 518   if(fMultiCutAna){
 519
 520     fhRePIDBits         = new TH2D*[fNPIDBits];
 521     for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
 522       snprintf(key,   buffersize,"hRe_pidbit%d",ipid) ;
 523       snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for PIDBit=%d",fPIDBits[ipid]) ;
 524       fhRePIDBits[ipid] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 525       fhRePIDBits[ipid]->SetXTitle("p_{T} (GeV/c)");
 526       fhRePIDBits[ipid]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 527       outputContainer->Add(fhRePIDBits[ipid]) ;
 528     }// pid bit loop
 529
 530     fhRePtNCellAsymCuts    = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
 531     fhRePtNCellAsymCutsSM0 = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
 532     fhRePtNCellAsymCutsSM1 = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
 533     fhRePtNCellAsymCutsSM2 = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
 534     fhRePtNCellAsymCutsSM3 = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
 535     fhMiPtNCellAsymCuts    = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
 536     for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
 537       for(Int_t icell=0; icell<fNCellNCuts; icell++){
 538         for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
 539           snprintf(key,   buffersize,"hRe_pt%d_cell%d_asym%d",ipt,icell,iasym) ;
 540           snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for pt >%2.2f, ncell>%d and asym >%1.2f ",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]) ;
 541           Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
 542           //printf("ipt %d, icell %d, iassym %d, index %d\n",ipt, icell, iasym, index);
 543           fhRePtNCellAsymCuts[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 544           fhRePtNCellAsymCuts[index]->SetXTitle("p_{T} (GeV/c)");
 545           fhRePtNCellAsymCuts[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 546           outputContainer->Add(fhRePtNCellAsymCuts[index]) ;
 547
 548           snprintf(key,   buffersize,"hRe_pt%d_cell%d_asym%d_SM0",ipt,icell,iasym) ;
 549           snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for pt >%2.2f, ncell>%d and asym >%1.2f, SM 0 ",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]) ;
 550           fhRePtNCellAsymCutsSM0[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 551           fhRePtNCellAsymCutsSM0[index]->SetXTitle("p_{T} (GeV/c)");
 552           fhRePtNCellAsymCutsSM0[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 553           outputContainer->Add(fhRePtNCellAsymCutsSM0[index]) ;
 554
 555           snprintf(key,   buffersize,"hRe_pt%d_cell%d_asym%d_SM1",ipt,icell,iasym) ;
 556           snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for pt >%2.2f, ncell>%d and asym >%1.2f, SM 1 ",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]) ;
 557           fhRePtNCellAsymCutsSM1[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 558           fhRePtNCellAsymCutsSM1[index]->SetXTitle("p_{T} (GeV/c)");
 559           fhRePtNCellAsymCutsSM1[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 560           outputContainer->Add(fhRePtNCellAsymCutsSM1[index]) ;
 561
 562           snprintf(key,   buffersize,"hRe_pt%d_cell%d_asym%d_SM2",ipt,icell,iasym) ;
 563           snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for pt >%2.2f, ncell>%d and asym >%1.2f, SM 2 ",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]) ;
 564           fhRePtNCellAsymCutsSM2[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 565           fhRePtNCellAsymCutsSM2[index]->SetXTitle("p_{T} (GeV/c)");
 566           fhRePtNCellAsymCutsSM2[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 567           outputContainer->Add(fhRePtNCellAsymCutsSM2[index]) ;
 568
 569           snprintf(key,   buffersize,"hRe_pt%d_cell%d_asym%d_SM3",ipt,icell,iasym) ;
 570           snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for pt >%2.2f, ncell>%d and asym >%1.2f, SM 3 ",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]) ;
 571           fhRePtNCellAsymCutsSM3[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 572           fhRePtNCellAsymCutsSM3[index]->SetXTitle("p_{T} (GeV/c)");
 573           fhRePtNCellAsymCutsSM3[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 574           outputContainer->Add(fhRePtNCellAsymCutsSM3[index]) ;
 575
 576           snprintf(key,   buffersize,"hMi_pt%d_cell%d_asym%d",ipt,icell,iasym) ;
 577           snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for pt >%2.2f, ncell>%d and asym >%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]) ;
 578           fhMiPtNCellAsymCuts[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 579           fhMiPtNCellAsymCuts[index]->SetXTitle("p_{T} (GeV/c)");
 580           fhMiPtNCellAsymCuts[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 581           outputContainer->Add(fhMiPtNCellAsymCuts[index]) ;
 582
 583         }
 584       }
 585     }
 586
 587     fhRePtMult = new TH3D*[fNAsymCuts] ;
 588     for(Int_t iasym = 0; iasym<fNAsymCuts; iasym++){
 589       fhRePtMult[iasym] = new TH3D(Form("hRePtMult_asym%d",iasym),Form("(p_{T},C,M)_{#gamma#gamma}, A<%1.2f",fAsymCuts[iasym]),
 590                                    nptbins,ptmin,ptmax,ntrmbins,ntrmmin,ntrmmax,nmassbins,massmin,massmax);
 591       fhRePtMult[iasym]->SetXTitle("p_{T} (GeV/c)");
 592       fhRePtMult[iasym]->SetYTitle("Track multiplicity");
 593       fhRePtMult[iasym]->SetZTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 594       outputContainer->Add(fhRePtMult[iasym]) ;
 595     }
 596
 597   }// multi cuts analysis
 598
 599   fhEvents=new TH3D("hEvents","Number of events",fNCentrBin,0.,1.*fNCentrBin,
 600                     GetNZvertBin(),0.,1.*GetNZvertBin(),GetNRPBin(),0.,1.*GetNRPBin()) ;
 601
 602   fhEvents->SetXTitle("Centrality bin");
 603   fhEvents->SetYTitle("Z vertex bin bin");
 604   fhEvents->SetZTitle("RP bin");
 605   outputContainer->Add(fhEvents) ;
 606
 607   fhCentrality=new TH1D("hCentralityBin","Number of events in centrality bin",fNCentrBin*10,0.,1.*fNCentrBin) ;
 608   fhCentrality->SetXTitle("Centrality bin");
 609   outputContainer->Add(fhCentrality) ;
 610
 611   fhRealOpeningAngle  = new TH2D
 612   ("hRealOpeningAngle","Angle between all #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,300,0,TMath::Pi());
 613   fhRealOpeningAngle->SetYTitle("#theta(rad)");
 614   fhRealOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
 615   outputContainer->Add(fhRealOpeningAngle) ;
 616
 617   fhRealCosOpeningAngle  = new TH2D
 618   ("hRealCosOpeningAngle","Cosinus of angle between all #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,100,0,1);
 619   fhRealCosOpeningAngle->SetYTitle("cos (#theta) ");
 620   fhRealCosOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
 621   outputContainer->Add(fhRealCosOpeningAngle) ;
 622
 623   if(fDoOwnMix){
 624
 625     fhMixedOpeningAngle  = new TH2D
 626     ("hMixedOpeningAngle","Angle between all #gamma pair vs E_{#pi^{0}}, Mixed pairs",nptbins,ptmin,ptmax,300,0,TMath::Pi());
 627     fhMixedOpeningAngle->SetYTitle("#theta(rad)");
 628     fhMixedOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
 629     outputContainer->Add(fhMixedOpeningAngle) ;
 630
 631     fhMixedCosOpeningAngle  = new TH2D
 632     ("hMixedCosOpeningAngle","Cosinus of angle between all #gamma pair vs E_{#pi^{0}}, Mixed pairs",nptbins,ptmin,ptmax,100,0,1);
 633     fhMixedCosOpeningAngle->SetYTitle("cos (#theta) ");
 634     fhMixedCosOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
 635     outputContainer->Add(fhMixedCosOpeningAngle) ;
 636
 637   }
 638
 639   //Histograms filled only if MC data is requested
 640   if(IsDataMC()){
 641     //Pi0
 642     fhPrimPi0Pt     = new TH1D("hPrimPi0Pt","Primary pi0 pt",nptbins,ptmin,ptmax) ;
 643     fhPrimPi0AccPt  = new TH1D("hPrimPi0AccPt","Primary pi0 pt with both photons in acceptance",nptbins,ptmin,ptmax) ;
 644     fhPrimPi0Pt   ->SetXTitle("p_{T} (GeV/c)");
 645     fhPrimPi0AccPt->SetXTitle("p_{T} (GeV/c)");
 646     outputContainer->Add(fhPrimPi0Pt) ;
 647     outputContainer->Add(fhPrimPi0AccPt) ;
 648
 649     fhPrimPi0Y      = new TH2D("hPrimPi0Rapidity","Rapidity of primary pi0",nptbins,ptmin,ptmax, netabins,etamin,etamax) ;
 650     fhPrimPi0Y   ->SetYTitle("Rapidity");
 651     fhPrimPi0Y   ->SetXTitle("p_{T} (GeV/c)");
 652     outputContainer->Add(fhPrimPi0Y) ;
 653
 654     fhPrimPi0AccY   = new TH2D("hPrimPi0AccRapidity","Rapidity of primary pi0",nptbins,ptmin,ptmax, netabins,etamin,etamax) ;
 655     fhPrimPi0AccY->SetYTitle("Rapidity");
 656     fhPrimPi0AccY->SetXTitle("p_{T} (GeV/c)");
 657     outputContainer->Add(fhPrimPi0AccY) ;
 658
 659     fhPrimPi0Phi    = new TH2D("hPrimPi0Phi","Azimuthal of primary pi0",nptbins,ptmin,ptmax, nphibins,phimin*TMath::RadToDeg(),phimax*TMath::RadToDeg()) ;
 660     fhPrimPi0Phi->SetYTitle("#phi (deg)");
 661     fhPrimPi0Phi->SetXTitle("p_{T} (GeV/c)");
 662     outputContainer->Add(fhPrimPi0Phi) ;
 663
 664     fhPrimPi0AccPhi = new TH2D("hPrimPi0AccPhi","Azimuthal of primary pi0 with accepted daughters",nptbins,ptmin,ptmax, nphibins,phimin*TMath::RadToDeg(),phimax*TMath::RadToDeg()) ;
 665     fhPrimPi0AccPhi->SetYTitle("#phi (deg)");
 666     fhPrimPi0AccPhi->SetXTitle("p_{T} (GeV/c)");
 667     outputContainer->Add(fhPrimPi0AccPhi) ;
 668
 669     //Eta
 670     fhPrimEtaPt     = new TH1D("hPrimEtaPt","Primary eta pt",nptbins,ptmin,ptmax) ;
 671     fhPrimEtaAccPt  = new TH1D("hPrimEtaAccPt","Primary eta pt with both photons in acceptance",nptbins,ptmin,ptmax) ;
 672     fhPrimEtaPt   ->SetXTitle("p_{T} (GeV/c)");
 673     fhPrimEtaAccPt->SetXTitle("p_{T} (GeV/c)");
 674     outputContainer->Add(fhPrimEtaPt) ;
 675     outputContainer->Add(fhPrimEtaAccPt) ;
 676
 677     fhPrimEtaY      = new TH2D("hPrimEtaRapidity","Rapidity of primary eta",nptbins,ptmin,ptmax,netabins,etamin,etamax) ;
 678     fhPrimEtaY->SetYTitle("Rapidity");
 679     fhPrimEtaY->SetXTitle("p_{T} (GeV/c)");
 680     outputContainer->Add(fhPrimEtaY) ;
 681
 682     fhPrimEtaAccY   = new TH2D("hPrimEtaAccRapidity","Rapidity of primary eta",nptbins,ptmin,ptmax, netabins,etamin,etamax) ;
 683     fhPrimEtaAccY->SetYTitle("Rapidity");
 684     fhPrimEtaAccY->SetXTitle("p_{T} (GeV/c)");
 685     outputContainer->Add(fhPrimEtaAccY) ;
 686
 687     fhPrimEtaPhi    = new TH2D("hPrimEtaPhi","Azimuthal of primary eta",nptbins,ptmin,ptmax, nphibins,phimin*TMath::RadToDeg(),phimax*TMath::RadToDeg()) ;
 688     fhPrimEtaPhi->SetYTitle("#phi (deg)");
 689     fhPrimEtaPhi->SetXTitle("p_{T} (GeV/c)");
 690     outputContainer->Add(fhPrimEtaPhi) ;
 691
 692     fhPrimEtaAccPhi = new TH2D("hPrimEtaAccPhi","Azimuthal of primary eta with accepted daughters",nptbins,ptmin,ptmax, nphibins,phimin*TMath::RadToDeg(),phimax*TMath::RadToDeg()) ;
 693     fhPrimEtaAccPhi->SetYTitle("#phi (deg)");
 694     fhPrimEtaAccPhi->SetXTitle("p_{T} (GeV/c)");
 695     outputContainer->Add(fhPrimEtaAccPhi) ;
 696
 697
 698     //Prim origin
 699     //Pi0
 700     fhPrimPi0PtOrigin     = new TH2D("hPrimPi0PtOrigin","Primary pi0 pt vs origin",nptbins,ptmin,ptmax,11,0,11) ;
 701     fhPrimPi0PtOrigin->SetXTitle("p_{T} (GeV/c)");
 702     fhPrimPi0PtOrigin->SetYTitle("Origin");
 703     fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(1 ,"Status 21");
 704     fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(2 ,"Quark");
 705     fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(3 ,"qq Resonances ");
 706     fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(4 ,"Resonances");
 707     fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(5 ,"#rho");
 708     fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(6 ,"#omega");
 709     fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(7 ,"K");
 710     fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(8 ,"Other");
 711     fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(9 ,"#eta");
 712     fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(10 ,"#eta prime");
 713     outputContainer->Add(fhPrimPi0PtOrigin) ;
 714
 715     fhMCPi0PtOrigin     = new TH2D("hMCPi0PtOrigin","Reconstructed pair from generated pi0 pt vs origin",nptbins,ptmin,ptmax,11,0,11) ;
 716     fhMCPi0PtOrigin->SetXTitle("p_{T} (GeV/c)");
 717     fhMCPi0PtOrigin->SetYTitle("Origin");
 718     fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(1 ,"Status 21");
 719     fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(2 ,"Quark");
 720     fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(3 ,"qq Resonances");
 721     fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(4 ,"Resonances");
 722     fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(5 ,"#rho");
 723     fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(6 ,"#omega");
 724     fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(7 ,"K");
 725     fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(8 ,"Other");
 726     fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(9 ,"#eta");
 727     fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(10 ,"#eta prime");
 728     outputContainer->Add(fhMCPi0PtOrigin) ;
 729
 730     //Eta
 731     fhPrimEtaPtOrigin     = new TH2D("hPrimEtaPtOrigin","Primary pi0 pt vs origin",nptbins,ptmin,ptmax,7,0,7) ;
 732     fhPrimEtaPtOrigin->SetXTitle("p_{T} (GeV/c)");
 733     fhPrimEtaPtOrigin->SetYTitle("Origin");
 734     fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(1 ,"Status 21");
 735     fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(2 ,"Quark");
 736     fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(3 ,"qq Resonances");
 737     fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(4 ,"Resonances");
 738     fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(5 ,"Other");
 739     fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(6 ,"#eta prime ");
 740
 741     outputContainer->Add(fhPrimEtaPtOrigin) ;
 742
 743     fhMCEtaPtOrigin     = new TH2D("hMCEtaPtOrigin","Reconstructed pair from generated pi0 pt vs origin",nptbins,ptmin,ptmax,7,0,7) ;
 744     fhMCEtaPtOrigin->SetXTitle("p_{T} (GeV/c)");
 745     fhMCEtaPtOrigin->SetYTitle("Origin");
 746     fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(1 ,"Status 21");
 747     fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(2 ,"Quark");
 748     fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(3 ,"qq Resonances");
 749     fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(4 ,"Resonances");
 750     fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(5 ,"Other");
 751     fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(6 ,"#eta prime");
 752
 753     outputContainer->Add(fhMCEtaPtOrigin) ;
 754
 755
 756     fhPrimPi0OpeningAngle  = new TH2D
 757     ("hPrimPi0OpeningAngle","Angle between all primary #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,100,0,0.5);
 758     fhPrimPi0OpeningAngle->SetYTitle("#theta(rad)");
 759     fhPrimPi0OpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
 760     outputContainer->Add(fhPrimPi0OpeningAngle) ;
 761
 762     fhPrimPi0CosOpeningAngle  = new TH2D
 763     ("hPrimPi0CosOpeningAngle","Cosinus of angle between all primary #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,100,-1,1);
 764     fhPrimPi0CosOpeningAngle->SetYTitle("cos (#theta) ");
 765     fhPrimPi0CosOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
 766     outputContainer->Add(fhPrimPi0CosOpeningAngle) ;
 767
 768     for(Int_t i = 0; i<13; i++){
 769       fhMCOrgMass[i] = new TH2D(Form("hMCOrgMass_%d",i),Form("mass vs pt, origin %d",i),nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 770       fhMCOrgMass[i]->SetXTitle("p_{T} (GeV/c)");
 771       fhMCOrgMass[i]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 772       outputContainer->Add(fhMCOrgMass[i]) ;
 773
 774       fhMCOrgAsym[i]= new TH2D(Form("hMCOrgAsym_%d",i),Form("asymmetry vs pt, origin %d",i),nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax) ;
 775       fhMCOrgAsym[i]->SetXTitle("p_{T} (GeV/c)");
 776       fhMCOrgAsym[i]->SetYTitle("A");
 777       outputContainer->Add(fhMCOrgAsym[i]) ;
 778
 779       fhMCOrgDeltaEta[i] = new TH2D(Form("hMCOrgDeltaEta_%d",i),Form("#Delta #eta of pair vs pt, origin %d",i),nptbins,ptmin,ptmax,netabins,-1.4,1.4) ;
 780       fhMCOrgDeltaEta[i]->SetXTitle("p_{T} (GeV/c)");
 781       fhMCOrgDeltaEta[i]->SetYTitle("#Delta #eta");
 782       outputContainer->Add(fhMCOrgDeltaEta[i]) ;
 783
 784       fhMCOrgDeltaPhi[i]= new TH2D(Form("hMCOrgDeltaPhi_%d",i),Form("#Delta #phi of pair vs p_{T}, origin %d",i),nptbins,ptmin,ptmax,nphibins,-0.7,0.7) ;
 785       fhMCOrgDeltaPhi[i]->SetXTitle("p_{T} (GeV/c)");
 786       fhMCOrgDeltaPhi[i]->SetYTitle("#Delta #phi (rad)");
 787       outputContainer->Add(fhMCOrgDeltaPhi[i]) ;
 788
 789     }
 790
 791     if(fMultiCutAnaSim){
 792       fhMCPi0MassPtTrue  = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
 793       fhMCPi0MassPtRec   = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
 794       fhMCPi0PtTruePtRec = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
 795       fhMCEtaMassPtRec   = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
 796       fhMCEtaMassPtTrue  = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
 797       fhMCEtaPtTruePtRec = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
 798       for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
 799         for(Int_t icell=0; icell<fNCellNCuts; icell++){
 800           for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
 801             Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
 802
 803             fhMCPi0MassPtRec[index] = new TH2D(Form("hMCPi0MassPtRec_pt%d_cell%d_asym%d",ipt,icell,iasym),
 804                                                Form("Reconstructed Mass vs reconstructed p_T of true #pi^{0} cluster pairs for pt >%2.2f, ncell>%d and asym >%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]),
 805                                                nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 806             fhMCPi0MassPtRec[index]->SetXTitle("p_{T, reconstructed} (GeV/c)");
 807             fhMCPi0MassPtRec[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 808             outputContainer->Add(fhMCPi0MassPtRec[index]) ;
 809
 810             fhMCPi0MassPtTrue[index] = new TH2D(Form("hMCPi0MassPtTrue_pt%d_cell%d_asym%d",ipt,icell,iasym),
 811                                                 Form("Reconstructed Mass vs generated p_T of true #pi^{0} cluster pairs for pt >%2.2f, ncell>%d and asym >%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]),
 812                                                 nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 813             fhMCPi0MassPtTrue[index]->SetXTitle("p_{T, generated} (GeV/c)");
 814             fhMCPi0MassPtTrue[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 815             outputContainer->Add(fhMCPi0MassPtTrue[index]) ;
 816
 817             fhMCPi0PtTruePtRec[index] = new TH2D(Form("hMCPi0PtTruePtRec_pt%d_cell%d_asym%d",ipt,icell,iasym),
 818                                                  Form("Generated vs reconstructed p_T of true #pi^{0} cluster pairs, 0.01 < rec. mass < 0.17 MeV/c^{2} for pt >%2.2f, ncell>%d and asym >%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]),
 819                                                  nptbins,ptmin,ptmax,nptbins,ptmin,ptmax) ;
 820             fhMCPi0PtTruePtRec[index]->SetXTitle("p_{T, generated} (GeV/c)");
 821             fhMCPi0PtTruePtRec[index]->SetYTitle("p_{T, reconstructed} (GeV/c)");
 822             outputContainer->Add(fhMCPi0PtTruePtRec[index]) ;
 823
 824             fhMCEtaMassPtRec[index] = new TH2D(Form("hMCEtaMassPtRec_pt%d_cell%d_asym%d",ipt,icell,iasym),
 825                                                Form("Reconstructed Mass vs reconstructed p_T of true #eta cluster pairs for pt >%2.2f, ncell>%d and asym >%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]),
 826                                                nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 827             fhMCEtaMassPtRec[index]->SetXTitle("p_{T, generated} (GeV/c)");
 828             fhMCEtaMassPtRec[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 829             outputContainer->Add(fhMCEtaMassPtRec[index]) ;
 830
 831             fhMCEtaMassPtTrue[index] = new TH2D(Form("hMCEtaMassPtTrue_pt%d_cell%d_asym%d",ipt,icell,iasym),
 832                                                 Form("Reconstructed Mass vs generated p_T of true #eta cluster pairs for pt >%2.2f, ncell>%d and asym >%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]),
 833                                                 nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 834             fhMCEtaMassPtTrue[index]->SetXTitle("p_{T, generated} (GeV/c)");
 835             fhMCEtaMassPtTrue[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 836             outputContainer->Add(fhMCEtaMassPtTrue[index]) ;
 837
 838             fhMCEtaPtTruePtRec[index] = new TH2D(Form("hMCEtaPtTruePtRec_pt%d_cell%d_asym%d",ipt,icell,iasym),
 839                                                  Form("Generated vs reconstructed p_T of true #eta cluster pairs, 0.01 < rec. mass < 0.17 MeV/c^{2} for pt >%2.2f, ncell>%d and asym >%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]),
 840                                                  nptbins,ptmin,ptmax,nptbins,ptmin,ptmax) ;
 841             fhMCEtaPtTruePtRec[index]->SetXTitle("p_{T, generated} (GeV/c)");
 842             fhMCEtaPtTruePtRec[index]->SetYTitle("p_{T, reconstructed} (GeV/c)");
 843             outputContainer->Add(fhMCEtaPtTruePtRec[index]) ;
 844           }
 845         }
 846       }
 847     }//multi cut ana
 848     else {
 849       fhMCPi0MassPtTrue  = new TH2D*[1];
 850       fhMCPi0PtTruePtRec = new TH2D*[1];
 851       fhMCEtaMassPtTrue  = new TH2D*[1];
 852       fhMCEtaPtTruePtRec = new TH2D*[1];
 853
 854       fhMCPi0MassPtTrue[0] = new TH2D("hMCPi0MassPtTrue","Reconstructed Mass vs generated p_T of true #pi^{0} cluster pairs",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 855       fhMCPi0MassPtTrue[0]->SetXTitle("p_{T, generated} (GeV/c)");
 856       fhMCPi0MassPtTrue[0]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 857       outputContainer->Add(fhMCPi0MassPtTrue[0]) ;
 858
 859       fhMCPi0PtTruePtRec[0]= new TH2D("hMCPi0PtTruePtRec","Generated vs reconstructed p_T of true #pi^{0} cluster pairs, 0.01 < rec. mass < 0.17 MeV/c^{2}",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax) ;
 860       fhMCPi0PtTruePtRec[0]->SetXTitle("p_{T, generated} (GeV/c)");
 861       fhMCPi0PtTruePtRec[0]->SetYTitle("p_{T, reconstructed} (GeV/c)");
 862       outputContainer->Add(fhMCPi0PtTruePtRec[0]) ;
 863
 864       fhMCEtaMassPtTrue[0] = new TH2D("hMCEtaMassPtTrue","Reconstructed Mass vs generated p_T of true #eta cluster pairs",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 865       fhMCEtaMassPtTrue[0]->SetXTitle("p_{T, generated} (GeV/c)");
 866       fhMCEtaMassPtTrue[0]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 867       outputContainer->Add(fhMCEtaMassPtTrue[0]) ;
 868
 869       fhMCEtaPtTruePtRec[0]= new TH2D("hMCEtaPtTruePtRec","Generated vs reconstructed p_T of true #eta cluster pairs, 0.01 < rec. mass < 0.17 MeV/c^{2}",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax) ;
 870       fhMCEtaPtTruePtRec[0]->SetXTitle("p_{T, generated} (GeV/c)");
 871       fhMCEtaPtTruePtRec[0]->SetYTitle("p_{T, reconstructed} (GeV/c)");
 872       outputContainer->Add(fhMCEtaPtTruePtRec[0]) ;
 873     }
 874   }
 875
 876   TString pairnamePHOS[] = {"(0-1)","(0-2)","(1-2)","(0-3)","(0-4)","(1-3)","(1-4)","(2-3)","(2-4)","(3-4)"};
 877   for(Int_t imod=0; imod<fNModules; imod++){
 878     //Module dependent invariant mass
 879     snprintf(key, buffersize,"hReMod_%d",imod) ;
 880     snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for Module %d",imod) ;
 881     fhReMod[imod]  = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 882     fhReMod[imod]->SetXTitle("p_{T} (GeV/c)");
 883     fhReMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 884     outputContainer->Add(fhReMod[imod]) ;
 885     if(fCalorimeter=="PHOS"){
 886       snprintf(key, buffersize,"hReDiffPHOSMod_%d",imod) ;
 887       snprintf(title, buffersize,"Real pairs PHOS, clusters in different Modules: %s",(pairnamePHOS[imod]).Data()) ;
 888       fhReDiffPHOSMod[imod]  = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 889       fhReDiffPHOSMod[imod]->SetXTitle("p_{T} (GeV/c)");
 890       fhReDiffPHOSMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 891       outputContainer->Add(fhReDiffPHOSMod[imod]) ;
 892     }
 893     else{//EMCAL
 894       if(imod<fNModules/2){
 895         snprintf(key, buffersize,"hReSameSectorEMCAL_%d",imod) ;
 896         snprintf(title, buffersize,"Real pairs EMCAL, clusters in same sector, SM(%d,%d)",imod*2,imod*2+1) ;
 897         fhReSameSectorEMCALMod[imod]  = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 898         fhReSameSectorEMCALMod[imod]->SetXTitle("p_{T} (GeV/c)");
 899         fhReSameSectorEMCALMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 900         outputContainer->Add(fhReSameSectorEMCALMod[imod]) ;
 901       }
 902       if(imod<fNModules-2){
 903         snprintf(key, buffersize,"hReSameSideEMCAL_%d",imod) ;
 904         snprintf(title, buffersize,"Real pairs EMCAL, clusters in same side SM(%d,%d)",imod, imod+2) ;
 905         fhReSameSideEMCALMod[imod]  = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 906         fhReSameSideEMCALMod[imod]->SetXTitle("p_{T} (GeV/c)");
 907         fhReSameSideEMCALMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 908         outputContainer->Add(fhReSameSideEMCALMod[imod]) ;
 909       }
 910     }//EMCAL
 911
 912     if(fDoOwnMix){
 913       snprintf(key, buffersize,"hMiMod_%d",imod) ;
 914       snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for Module %d",imod) ;
 915       fhMiMod[imod]  = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 916       fhMiMod[imod]->SetXTitle("p_{T} (GeV/c)");
 917       fhMiMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 918       outputContainer->Add(fhMiMod[imod]) ;
 919
 920       if(fCalorimeter=="PHOS"){
 921         snprintf(key, buffersize,"hMiDiffPHOSMod_%d",imod) ;
 922         snprintf(title, buffersize,"Mixed pairs PHOS, clusters in different Modules: %s",(pairnamePHOS[imod]).Data()) ;
 923         fhMiDiffPHOSMod[imod]  = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 924         fhMiDiffPHOSMod[imod]->SetXTitle("p_{T} (GeV/c)");
 925         fhMiDiffPHOSMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 926         outputContainer->Add(fhMiDiffPHOSMod[imod]) ;
 927       }//PHOS
 928       else{//EMCAL
 929         if(imod<fNModules/2){
 930           snprintf(key, buffersize,"hMiSameSectorEMCALMod_%d",imod) ;
 931           snprintf(title, buffersize,"Mixed pairs EMCAL, clusters in same sector, SM(%d,%d)",imod*2,imod*2+1) ;
 932           fhMiSameSectorEMCALMod[imod]  = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 933           fhMiSameSectorEMCALMod[imod]->SetXTitle("p_{T} (GeV/c)");
 934           fhMiSameSectorEMCALMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 935           outputContainer->Add(fhMiSameSectorEMCALMod[imod]) ;
 936         }
 937         if(imod<fNModules-2){
 938           snprintf(key, buffersize,"hMiSameSideEMCALMod_%d",imod) ;
 939           snprintf(title, buffersize,"Mixed pairs EMCAL, clusters in same side SM(%d,%d)",imod, imod+2) ;
 940           fhMiSameSideEMCALMod[imod]  = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 941           fhMiSameSideEMCALMod[imod]->SetXTitle("p_{T} (GeV/c)");
 942           fhMiSameSideEMCALMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 943           outputContainer->Add(fhMiSameSideEMCALMod[imod]) ;
 944         }
 945       }//EMCAL
 946      }
 947   }
 948
 949 //  for(Int_t i = 0; i < outputContainer->GetEntries() ; i++){
 950 //
 951 //    printf("Histogram %d, name: %s\n ",i, outputContainer->At(i)->GetName());
 952 //
 953 //  }
 954
 955   return outputContainer;
 956 }
 957
 958 //_________________________________________________________________________________________________________________________________________________
 959 void AliAnaPi0::Print(const Option_t * /*opt*/) const
 960 {
 961   //Print some relevant parameters set for the analysis
 962   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
 963   AliAnaPartCorrBaseClass::Print(" ");
 964
 965   printf("Number of bins in Centrality:  %d \n",fNCentrBin) ;
 966   printf("Number of bins in Z vert. pos: %d \n",GetNZvertBin()) ;
 967   printf("Number of bins in Reac. Plain: %d \n",GetNRPBin()) ;
 968   printf("Depth of event buffer: %d \n",fNmaxMixEv) ;
 969   printf("Pair in same Module: %d \n",fSameSM) ;
 970   printf("Cuts: \n") ;
 971 // printf("Z vertex position: -%2.3f < z < %2.3f \n",GetZvertexCut(),GetZvertexCut()) ; //It crashes here, why?
 972   printf("Number of modules:             %d \n",fNModules) ;
 973   printf("Select pairs with their angle: %d, edep %d, min angle %2.3f, max angle %2.3f \n",fUseAngleCut, fUseAngleEDepCut, fAngleCut, fAngleMaxCut) ;
 974   printf("Asymmetry cuts: n = %d, \n",fNAsymCuts) ;
 975   printf("\tasymmetry < ");
 976   for(Int_t i = 0; i < fNAsymCuts; i++) printf("%2.2f ",fAsymCuts[i]);
 977   printf("\n");
 978
 979   printf("PID selection bits: n = %d, \n",fNPIDBits) ;
 980   printf("\tPID bit = ");
 981   for(Int_t i = 0; i < fNPIDBits; i++) printf("%d ",fPIDBits[i]);
 982   printf("\n");
 983
 984   if(fMultiCutAna){
 985     printf("pT cuts: n = %d, \n",fNPtCuts) ;
 986     printf("\tpT > ");
 987     for(Int_t i = 0; i < fNPtCuts; i++) printf("%2.2f ",fPtCuts[i]);
 988     printf("GeV/c\n");
 989
 990     printf("N cell in cluster cuts: n = %d, \n",fNCellNCuts) ;
 991     printf("\tnCell > ");
 992     for(Int_t i = 0; i < fNCellNCuts; i++) printf("%d ",fCellNCuts[i]);
 993     printf("\n");
 994
 995   }
 996   printf("------------------------------------------------------\n") ;
 997 }
 998
 999 //_____________________________________________________________
1000 void AliAnaPi0::FillAcceptanceHistograms(){
1001   //Fill acceptance histograms if MC data is available
1002
1003   if(GetReader()->ReadStack()){
1004     AliStack * stack = GetMCStack();
1005     if(stack){
1006       for(Int_t i=0 ; i<stack->GetNprimary(); i++){
1007         TParticle * prim = stack->Particle(i) ;
1008         Int_t pdg = prim->GetPdgCode();
1009         if( pdg == 111 || pdg == 221){
1010           Double_t pi0Pt = prim->Pt() ;
1011           //printf("pi0, pt %2.2f\n",pi0Pt);
1012           if(prim->Energy() == TMath::Abs(prim->Pz()))  continue ; //Protection against floating point exception
1013           Double_t pi0Y  = 0.5*TMath::Log((prim->Energy()-prim->Pz())/(prim->Energy()+prim->Pz())) ;
1014           Double_t phi   = TMath::RadToDeg()*prim->Phi() ;
1015           if(pdg == 111){
1016             if(TMath::Abs(pi0Y) < 1.0){
1017               fhPrimPi0Pt->Fill(pi0Pt) ;
1018             }
1019             fhPrimPi0Y  ->Fill(pi0Pt, pi0Y) ;
1020             fhPrimPi0Phi->Fill(pi0Pt, phi) ;
1021           }
1022           else if(pdg == 221){
1023             if(TMath::Abs(pi0Y) < 1.0){
1024               fhPrimEtaPt->Fill(pi0Pt) ;
1025             }
1026             fhPrimEtaY  ->Fill(pi0Pt, pi0Y) ;
1027             fhPrimEtaPhi->Fill(pi0Pt, phi) ;
1028           }
1029
1030           //Origin of meson
1031           Int_t momindex  = prim->GetFirstMother();
1032           if(momindex < 0) continue;
1033           TParticle* mother = stack->Particle(momindex);
1034           Int_t mompdg    = TMath::Abs(mother->GetPdgCode());
1035           Int_t momstatus = mother->GetStatusCode();
1036           if(pdg == 111){
1037             if     (momstatus  == 21)fhPrimPi0PtOrigin->Fill(pi0Pt,0.5);//parton
1038             else if(mompdg     < 22 ) fhPrimPi0PtOrigin->Fill(pi0Pt,1.5);//quark
1039             else if(mompdg     > 2100  && mompdg   < 2210) fhPrimPi0PtOrigin->Fill(pi0Pt,2.5);// resonances
1040             else if(mompdg    == 221) fhPrimPi0PtOrigin->Fill(pi0Pt,8.5);//eta
1041             else if(mompdg    == 331) fhPrimPi0PtOrigin->Fill(pi0Pt,9.5);//eta prime
1042             else if(mompdg    == 213) fhPrimPi0PtOrigin->Fill(pi0Pt,4.5);//rho
1043             else if(mompdg    == 223) fhPrimPi0PtOrigin->Fill(pi0Pt,5.5);//omega
1044             else if(mompdg    >= 310   && mompdg    <= 323) fhPrimPi0PtOrigin->Fill(pi0Pt,6.5);//k0S, k+-,k*
1045             else if(mompdg    == 130) fhPrimPi0PtOrigin->Fill(pi0Pt,6.5);//k0L
1046             else if(momstatus == 11 || momstatus  == 12 ) fhPrimPi0PtOrigin->Fill(pi0Pt,3.5);//resonances
1047             else                      fhPrimPi0PtOrigin->Fill(pi0Pt,7.5);//other?
1048           }//pi0
1049           else {
1050             if     (momstatus == 21 ) fhPrimEtaPtOrigin->Fill(pi0Pt,0.5);//parton
1051             else if(mompdg    < 22  ) fhPrimEtaPtOrigin->Fill(pi0Pt,1.5);//quark
1052             else if(mompdg    > 2100  && mompdg   < 2210) fhPrimEtaPtOrigin->Fill(pi0Pt,2.5);//qq resonances
1053             else if(mompdg    == 331) fhPrimEtaPtOrigin->Fill(pi0Pt,5.5);//eta prime
1054             else if(momstatus == 11 || momstatus  == 12 ) fhPrimEtaPtOrigin->Fill(pi0Pt,3.5);//resonances
1055             else fhPrimEtaPtOrigin->Fill(pi0Pt,4.5);//stable, conversions?
1056             //printf("Other Meson pdg %d, Mother %s, pdg %d, status %d\n",pdg, TDatabasePDG::Instance()->GetParticle(mompdg)->GetName(),mompdg, momstatus );
1057           }
1058
1059
1060           //Check if both photons hit Calorimeter
1061           if(prim->GetNDaughters()!=2) return; //Only interested in 2 gamma decay
1062           Int_t iphot1=prim->GetFirstDaughter() ;
1063           Int_t iphot2=prim->GetLastDaughter() ;
1064           if(iphot1>-1 && iphot1<stack->GetNtrack() && iphot2>-1 && iphot2<stack->GetNtrack()){
1065             TParticle * phot1 = stack->Particle(iphot1) ;
1066             TParticle * phot2 = stack->Particle(iphot2) ;
1067             if(phot1 && phot2 && phot1->GetPdgCode()==22 && phot2->GetPdgCode()==22){
1068               //printf("2 photons: photon 1: pt %2.2f, phi %3.2f, eta %1.2f; photon 2: pt %2.2f, phi %3.2f, eta %1.2f\n",
1069               //        phot1->Pt(), phot1->Phi()*180./3.1415, phot1->Eta(), phot2->Pt(), phot2->Phi()*180./3.1415, phot2->Eta());
1070
1071               TLorentzVector lv1, lv2;
1072               phot1->Momentum(lv1);
1073               phot2->Momentum(lv2);
1074
1075               Bool_t inacceptance = kFALSE;
1076               if(fCalorimeter == "PHOS"){
1077                 if(GetPHOSGeometry() && GetCaloUtils()->IsPHOSGeoMatrixSet()){
1078                   Int_t mod ;
1079                   Double_t x,z ;
1080                   if(GetPHOSGeometry()->ImpactOnEmc(phot1,mod,z,x) && GetPHOSGeometry()->ImpactOnEmc(phot2,mod,z,x))
1081                     inacceptance = kTRUE;
1082                   if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
1083                 }
1084                 else{
1085
1086                   if(GetFiducialCut()->IsInFiducialCut(lv1,fCalorimeter) && GetFiducialCut()->IsInFiducialCut(lv2,fCalorimeter))
1087                     inacceptance = kTRUE ;
1088                   if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
1089                 }
1090
1091               }
1092               else if(fCalorimeter == "EMCAL" && GetCaloUtils()->IsEMCALGeoMatrixSet()){
1093                 if(GetEMCALGeometry()){
1094
1095                   Int_t absID1=0;
1096                   Int_t absID2=0;
1097
1098                   GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(phot1->Eta(),phot1->Phi(),absID1);
1099                   GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(phot2->Eta(),phot2->Phi(),absID2);
1100
1101                   if( absID1 >= 0 && absID2 >= 0)
1102                     inacceptance = kTRUE;
1103
1104 //                  if(GetEMCALGeometry()->Impact(phot1) && GetEMCALGeometry()->Impact(phot2))
1105 //                    inacceptance = kTRUE;
1106                   if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
1107                 }
1108                 else{
1109                   if(GetFiducialCut()->IsInFiducialCut(lv1,fCalorimeter) && GetFiducialCut()->IsInFiducialCut(lv2,fCalorimeter))
1110                     inacceptance = kTRUE ;
1111                   if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
1112                 }
1113               }
1114
1115               if(inacceptance){
1116                 if(pdg==111){
1117                   fhPrimPi0AccPt ->Fill(pi0Pt) ;
1118                   fhPrimPi0AccPhi->Fill(pi0Pt, phi) ;
1119                   fhPrimPi0AccY  ->Fill(pi0Pt, pi0Y) ;
1120                   Double_t angle  = lv1.Angle(lv2.Vect());
1121                   fhPrimPi0OpeningAngle   ->Fill(pi0Pt,angle);
1122                   fhPrimPi0CosOpeningAngle->Fill(pi0Pt,TMath::Cos(angle));
1123                 }
1124                 else if(pdg==221){
1125                   fhPrimEtaAccPt ->Fill(pi0Pt) ;
1126                   fhPrimEtaAccPhi->Fill(pi0Pt, phi) ;
1127                   fhPrimEtaAccY  ->Fill(pi0Pt, pi0Y) ;
1128                 }
1129               }//Accepted
1130             }// 2 photons
1131           }//Check daughters exist
1132         }// Primary pi0 or eta
1133       }//loop on primaries
1134     }//stack exists and data is MC
1135   }//read stack
1136   else if(GetReader()->ReadAODMCParticles()){
1137
1138     TClonesArray * mcparticles = GetReader()->GetAODMCParticles(0);
1139     if(mcparticles){
1140       Int_t nprim = mcparticles->GetEntriesFast();
1141       for(Int_t i=0; i < nprim; i++)
1142       {
1143         AliAODMCParticle * prim = (AliAODMCParticle *) mcparticles->At(i);
1144         // Only generator particles
1145         if( prim->GetStatus() == 0) break;
1146
1147         Int_t pdg = prim->GetPdgCode();
1148         if( pdg == 111 || pdg == 221){
1149           Double_t pi0Pt = prim->Pt() ;
1150           //printf("pi0, pt %2.2f\n",pi0Pt);
1151           if(prim->E() == TMath::Abs(prim->Pz()))  continue ; //Protection against floating point exception
1152           Double_t pi0Y  = 0.5*TMath::Log((prim->E()-prim->Pz())/(prim->E()+prim->Pz())) ;
1153           Double_t phi   = TMath::RadToDeg()*prim->Phi() ;
1154           if(pdg == 111){
1155             if(TMath::Abs(pi0Y) < 0.5){
1156               fhPrimPi0Pt->Fill(pi0Pt) ;
1157             }
1158             fhPrimPi0Y  ->Fill(pi0Pt, pi0Y) ;
1159             fhPrimPi0Phi->Fill(pi0Pt, phi) ;
1160           }
1161           else if(pdg == 221){
1162             if(TMath::Abs(pi0Y) < 0.5){
1163               fhPrimEtaPt->Fill(pi0Pt) ;
1164             }
1165             fhPrimEtaY  ->Fill(pi0Pt, pi0Y) ;
1166             fhPrimEtaPhi->Fill(pi0Pt, phi) ;
1167           }
1168
1169           //Origin of meson
1170           Int_t momindex  = prim->GetMother();
1171           if(momindex < 0) continue;
1172           AliAODMCParticle* mother = (AliAODMCParticle *) mcparticles->At(momindex);
1173           Int_t mompdg    = TMath::Abs(mother->GetPdgCode());
1174           Int_t momstatus = mother->GetStatus();
1175             if(pdg == 111){
1176               if     (momstatus  == 21) fhPrimPi0PtOrigin->Fill(pi0Pt,0.5);//parton
1177               else if(mompdg     < 22 ) fhPrimPi0PtOrigin->Fill(pi0Pt,1.5);//quark
1178               else if(mompdg     > 2100  && mompdg   < 2210) fhPrimPi0PtOrigin->Fill(pi0Pt,2.5);// resonances
1179               else if(mompdg    == 221) fhPrimPi0PtOrigin->Fill(pi0Pt,8.5);//eta
1180               else if(mompdg    == 331) fhPrimPi0PtOrigin->Fill(pi0Pt,9.5);//eta prime
1181               else if(mompdg    == 213) fhPrimPi0PtOrigin->Fill(pi0Pt,4.5);//rho
1182               else if(mompdg    == 223) fhPrimPi0PtOrigin->Fill(pi0Pt,5.5);//omega
1183               else if(mompdg    >= 310   && mompdg    <= 323) fhPrimPi0PtOrigin->Fill(pi0Pt,6.5);//k0S, k+-,k*
1184               else if(mompdg    == 130) fhPrimPi0PtOrigin->Fill(pi0Pt,6.5);//k0L
1185               else if(momstatus == 11 || momstatus  == 12 ) fhPrimPi0PtOrigin->Fill(pi0Pt,3.5);//resonances
1186               else                      fhPrimPi0PtOrigin->Fill(pi0Pt,7.5);//other?
1187             }//pi0
1188             else {
1189               if     (momstatus == 21 ) fhPrimEtaPtOrigin->Fill(pi0Pt,0.5);//parton
1190               else if(mompdg    < 22  ) fhPrimEtaPtOrigin->Fill(pi0Pt,1.5);//quark
1191               else if(mompdg    > 2100  && mompdg   < 2210) fhPrimEtaPtOrigin->Fill(pi0Pt,2.5);//qq resonances
1192               else if(mompdg    == 331) fhPrimEtaPtOrigin->Fill(pi0Pt,5.5);//eta prime
1193               else if(momstatus == 11 || momstatus  == 12 ) fhPrimEtaPtOrigin->Fill(pi0Pt,3.5);//resonances
1194               else fhPrimEtaPtOrigin->Fill(pi0Pt,4.5);//stable, conversions?
1195               //printf("Other Meson pdg %d, Mother %s, pdg %d, status %d\n",pdg, TDatabasePDG::Instance()->GetParticle(mompdg)->GetName(),mompdg, momstatus );
1196             }
1197
1198           //Check if both photons hit Calorimeter
1199           if(prim->GetNDaughters()!=2) return; //Only interested in 2 gamma decay
1200           Int_t iphot1=prim->GetDaughter(0) ;
1201           Int_t iphot2=prim->GetDaughter(1) ;
1202           if(iphot1>-1 && iphot1<nprim && iphot2>-1 && iphot2<nprim){
1203             AliAODMCParticle * phot1 = (AliAODMCParticle *) mcparticles->At(iphot1);
1204             AliAODMCParticle * phot2 = (AliAODMCParticle *) mcparticles->At(iphot2);
1205             if(phot1 && phot2 && phot1->GetPdgCode()==22 && phot2->GetPdgCode()==22){
1206               TLorentzVector lv1, lv2;
1207               lv1.SetPxPyPzE(phot1->Px(),phot1->Py(),phot1->Pz(),phot1->E());
1208               lv2.SetPxPyPzE(phot2->Px(),phot2->Py(),phot2->Pz(),phot2->E());
1209
1210               Bool_t inacceptance = kFALSE;
1211               if(fCalorimeter == "PHOS"){
1212                 if(GetPHOSGeometry() && GetCaloUtils()->IsPHOSGeoMatrixSet()){
1213                   Int_t mod ;
1214                   Double_t x,z ;
1215                   Double_t vtx []={phot1->Xv(),phot1->Yv(),phot1->Zv()};
1216                   Double_t vtx2[]={phot2->Xv(),phot2->Yv(),phot2->Zv()};
1217                   if(GetPHOSGeometry()->ImpactOnEmc(vtx, phot1->Theta(),phot1->Phi(),mod,z,x) &&
1218                      GetPHOSGeometry()->ImpactOnEmc(vtx2,phot2->Theta(),phot2->Phi(),mod,z,x))
1219                     inacceptance = kTRUE;
1220                   if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
1221                 }
1222                 else{
1223
1224                   if(GetFiducialCut()->IsInFiducialCut(lv1,fCalorimeter) && GetFiducialCut()->IsInFiducialCut(lv2,fCalorimeter))
1225                     inacceptance = kTRUE ;
1226                   if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
1227                 }
1228
1229               }
1230               else if(fCalorimeter == "EMCAL" && GetCaloUtils()->IsEMCALGeoMatrixSet()){
1231                 if(GetEMCALGeometry()){
1232
1233                   Int_t absID1=0;
1234                   Int_t absID2=0;
1235
1236                   //TVector3 vtx(phot1->Xv(),phot1->Yv(),phot1->Zv());
1237                   //TVector3 vimpact(0,0,0);
1238
1239                   //GetEMCALGeometry()->ImpactOnEmcal(vtx,phot1->Theta(),phot1->Phi(),absID1,vimpact);
1240                   //TVector3 vtx2(phot2->Xv(),phot2->Yv(),phot2->Zv());
1241                   //TVector3 vimpact2(0,0,0);
1242                   //GetEMCALGeometry()->ImpactOnEmcal(vtx2,phot2->Theta(),phot2->Phi(),absID2,vimpact2);
1243
1244                   GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(phot1->Eta(),phot1->Phi(),absID1);
1245                   GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(phot2->Eta(),phot2->Phi(),absID2);
1246
1247 //                  if(TMath::Abs(phot1->Eta()) < 0.7 && phot1->Phi() > 80*TMath::DegToRad() && phot1->Phi() < 120*TMath::DegToRad() )
1248 //                    printf("Phot1 ccepted? %d\n",absID1);
1249 //                  if(TMath::Abs(phot2->Eta()) < 0.7 && phot2->Phi() > 80*TMath::DegToRad() && phot2->Phi() < 120*TMath::DegToRad() )
1250 //                    printf("Phot2 accepted? %d\n",absID2);
1251
1252                   if( absID1 >= 0 && absID2 >= 0)
1253                     inacceptance = kTRUE;
1254
1255 //                  if(pdg==111 && inacceptance) printf("2 photons: photon 1: absId %d, pt %2.2f, phi %3.2f, eta %1.2f; photon 2: absId %d, pt %2.2f, phi %3.2f, eta %1.2f\n",
1256 //                                                       absID1,phot1->Pt(), phot1->Phi()*TMath::RadToDeg(), phot1->Eta(),
1257 //                                                       absID2,phot2->Pt(), phot2->Phi()*TMath::RadToDeg(), phot2->Eta());
1258
1259
1260
1261                   if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
1262                 }
1263                 else{
1264                   if(GetFiducialCut()->IsInFiducialCut(lv1,fCalorimeter) && GetFiducialCut()->IsInFiducialCut(lv2,fCalorimeter))
1265                     inacceptance = kTRUE ;
1266                   if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
1267                 }
1268               }
1269
1270               if(inacceptance){
1271                 if(pdg==111){
1272   //                printf("ACCEPTED pi0: pt %2.2f, phi %3.2f, eta %1.2f\n",pi0Pt,phi,pi0Y);
1273                   fhPrimPi0AccPt ->Fill(pi0Pt) ;
1274                   fhPrimPi0AccPhi->Fill(pi0Pt, phi) ;
1275                   fhPrimPi0AccY  ->Fill(pi0Pt, pi0Y) ;
1276                   Double_t angle  = lv1.Angle(lv2.Vect());
1277                   fhPrimPi0OpeningAngle   ->Fill(pi0Pt,angle);
1278                   fhPrimPi0CosOpeningAngle->Fill(pi0Pt,TMath::Cos(angle));
1279                 }
1280                 else if(pdg==221){
1281                   fhPrimEtaAccPt ->Fill(pi0Pt) ;
1282                   fhPrimEtaAccPhi->Fill(pi0Pt, phi) ;
1283                   fhPrimEtaAccY  ->Fill(pi0Pt, pi0Y) ;
1284                 }
1285               }//Accepted
1286             }// 2 photons
1287           }//Check daughters exist
1288         }// Primary pi0 or eta
1289       }//loop on primaries
1290     }//stack exists and data is MC
1291
1292
1293   }     // read AOD MC
1294 }
1295
1296 //_____________________________________________________________
1297 void AliAnaPi0::FillMCVersusRecDataHistograms(const Int_t index1,  const Int_t index2,
1298                                               const Float_t pt1,   const Float_t pt2,
1299                                               const Int_t ncell1,  const Int_t ncell2,
1300                                               const Double_t mass, const Double_t pt, const Double_t asym,
1301                                               const Double_t deta, const Double_t dphi){
1302   //Do some MC checks on the origin of the pair, is there any common ancestor and if there is one, who?
1303   //Adjusted for Pythia, need to see what to do for other generators.
1304   //Array of histograms ordered as follows: 0-Photon, 1-electron, 2-pi0, 3-eta, 4-a-proton, 5-a-neutron, 6-stable particles,
1305   // 7-other decays, 8-string, 9-final parton, 10-initial parton, intermediate, 11-colliding proton, 12-unrelated
1306
1307   Int_t ancPDG    = 0;
1308   Int_t ancStatus = 0;
1309   TLorentzVector ancMomentum;
1310   Int_t ancLabel  = GetMCAnalysisUtils()->CheckCommonAncestor(index1, index2,
1311                                                               GetReader(), ancPDG, ancStatus,ancMomentum);
1312
1313   Int_t momindex  = -1;
1314   Int_t mompdg    = -1;
1315   Int_t momstatus = -1;
1316   if(GetDebug() > 1) printf("AliAnaPi0::FillMCVersusRecDataHistograms() - Common ancestor label %d, pdg %d, name %s, status %d; \n",
1317                             ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus);
1318
1319     if(ancLabel > -1){
1320       if(ancPDG==22){//gamma
1321         fhMCOrgMass[0]->Fill(pt,mass);
1322         fhMCOrgAsym[0]->Fill(pt,asym);
1323         fhMCOrgDeltaEta[0]->Fill(pt,deta);
1324         fhMCOrgDeltaPhi[0]->Fill(pt,dphi);
1325       }
1326       else if(TMath::Abs(ancPDG)==11){//e
1327         fhMCOrgMass[1]->Fill(pt,mass);
1328         fhMCOrgAsym[1]->Fill(pt,asym);
1329         fhMCOrgDeltaEta[1]->Fill(pt,deta);
1330         fhMCOrgDeltaPhi[1]->Fill(pt,dphi);
1331       }
1332       else if(ancPDG==111){//Pi0
1333         fhMCOrgMass[2]->Fill(pt,mass);
1334         fhMCOrgAsym[2]->Fill(pt,asym);
1335         fhMCOrgDeltaEta[2]->Fill(pt,deta);
1336         fhMCOrgDeltaPhi[2]->Fill(pt,dphi);
1337         if(fMultiCutAnaSim){
1338           for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
1339             for(Int_t icell=0; icell<fNCellNCuts; icell++){
1340               for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
1341                 Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
1342                 if(pt1    >  fPtCuts[ipt]      && pt2    >  fPtCuts[ipt]        &&
1343                    asym   <  fAsymCuts[iasym]                                   &&
1344                    ncell1 >= fCellNCuts[icell] && ncell2 >= fCellNCuts[icell]){
1345                   fhMCPi0MassPtRec [index]->Fill(pt,mass);
1346                   fhMCPi0MassPtTrue[index]->Fill(ancMomentum.Pt(),mass);
1347                   if(mass < 0.17 && mass > 0.1) fhMCPi0PtTruePtRec[index]->Fill(ancMomentum.Pt(),pt);
1348                 }//pass the different cuts
1349               }// pid bit cut loop
1350             }// icell loop
1351           }// pt cut loop
1352         }//Multi cut ana sim
1353         else {
1354           fhMCPi0MassPtTrue[0]->Fill(ancMomentum.Pt(),mass);
1355           if(mass < 0.17 && mass > 0.1) {
1356             fhMCPi0PtTruePtRec[0]->Fill(ancMomentum.Pt(),pt);
1357
1358             if(GetReader()->ReadStack()){
1359               TParticle* ancestor = GetMCStack()->Particle(ancLabel);
1360               momindex  = ancestor->GetFirstMother();
1361               if(momindex < 0) return;
1362               TParticle* mother = GetMCStack()->Particle(momindex);
1363               mompdg    = TMath::Abs(mother->GetPdgCode());
1364               momstatus = mother->GetStatusCode();
1365             }
1366             else {
1367               TClonesArray * mcparticles = GetReader()->GetAODMCParticles(0);
1368               AliAODMCParticle* ancestor = (AliAODMCParticle *) mcparticles->At(ancLabel);
1369               momindex  = ancestor->GetMother();
1370               if(momindex < 0) return;
1371               AliAODMCParticle* mother = (AliAODMCParticle *) mcparticles->At(momindex);
1372               mompdg    = TMath::Abs(mother->GetPdgCode());
1373               momstatus = mother->GetStatus();
1374             }
1375
1376             if     (momstatus  == 21) fhMCPi0PtOrigin->Fill(pt,0.5);//parton
1377             else if(mompdg     < 22 ) fhMCPi0PtOrigin->Fill(pt,1.5);//quark
1378             else if(mompdg     > 2100  && mompdg   < 2210) fhMCPi0PtOrigin->Fill(pt,2.5);// resonances
1379             else if(mompdg    == 221) fhMCPi0PtOrigin->Fill(pt,8.5);//eta
1380             else if(mompdg    == 331) fhMCPi0PtOrigin->Fill(pt,9.5);//eta prime
1381             else if(mompdg    == 213) fhMCPi0PtOrigin->Fill(pt,4.5);//rho
1382             else if(mompdg    == 223) fhMCPi0PtOrigin->Fill(pt,5.5);//omega
1383             else if(mompdg    >= 310   && mompdg    <= 323) fhMCPi0PtOrigin->Fill(pt,6.5);//k0S, k+-,k*
1384             else if(mompdg    == 130) fhMCPi0PtOrigin->Fill(pt,6.5);//k0L
1385             else if(momstatus == 11 || momstatus  == 12 ) fhMCPi0PtOrigin->Fill(pt,3.5);//resonances
1386             else                      fhMCPi0PtOrigin->Fill(pt,7.5);//other?
1387
1388           }//pi0 mass region
1389
1390         }
1391       }
1392       else if(ancPDG==221){//Eta
1393         fhMCOrgMass[3]->Fill(pt,mass);
1394         fhMCOrgAsym[3]->Fill(pt,asym);
1395         fhMCOrgDeltaEta[3]->Fill(pt,deta);
1396         fhMCOrgDeltaPhi[3]->Fill(pt,dphi);
1397         if(fMultiCutAnaSim){
1398           for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
1399             for(Int_t icell=0; icell<fNCellNCuts; icell++){
1400               for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
1401                 Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
1402                 if(pt1    >  fPtCuts[ipt]      && pt2    >  fPtCuts[ipt]        &&
1403                    asym   <  fAsymCuts[iasym]                                   &&
1404                    ncell1 >= fCellNCuts[icell] && ncell2 >= fCellNCuts[icell]){
1405                   fhMCEtaMassPtRec [index]->Fill(pt,mass);
1406                   fhMCEtaMassPtTrue[index]->Fill(ancMomentum.Pt(),mass);
1407                   if(mass < 0.65 && mass > 0.45) fhMCEtaPtTruePtRec[index]->Fill(ancMomentum.Pt(),pt);
1408                 }//pass the different cuts
1409               }// pid bit cut loop
1410             }// icell loop
1411           }// pt cut loop
1412         } //Multi cut ana sim
1413         else {
1414           fhMCEtaMassPtTrue[0]->Fill(ancMomentum.Pt(),mass);
1415           if(mass < 0.65 && mass > 0.45) fhMCEtaPtTruePtRec[0]->Fill(ancMomentum.Pt(),pt);
1416
1417           if(GetReader()->ReadStack()){
1418             TParticle* ancestor = GetMCStack()->Particle(ancLabel);
1419             momindex  = ancestor->GetFirstMother();
1420             if(momindex < 0) return;
1421             TParticle* mother = GetMCStack()->Particle(momindex);
1422             mompdg    = TMath::Abs(mother->GetPdgCode());
1423             momstatus = mother->GetStatusCode();
1424           }
1425           else {
1426             TClonesArray * mcparticles = GetReader()->GetAODMCParticles(0);
1427             AliAODMCParticle* ancestor = (AliAODMCParticle *) mcparticles->At(ancLabel);
1428             momindex  = ancestor->GetMother();
1429             if(momindex < 0) return;
1430             AliAODMCParticle* mother = (AliAODMCParticle *) mcparticles->At(momindex);
1431             mompdg    = TMath::Abs(mother->GetPdgCode());
1432             momstatus = mother->GetStatus();
1433           }
1434
1435           if     (momstatus == 21 ) fhMCEtaPtOrigin->Fill(pt,0.5);//parton
1436           else if(mompdg    < 22  ) fhMCEtaPtOrigin->Fill(pt,1.5);//quark
1437           else if(mompdg    > 2100  && mompdg   < 2210) fhMCEtaPtOrigin->Fill(pt,2.5);//qq resonances
1438           else if(mompdg    == 331) fhMCEtaPtOrigin->Fill(pt,5.5);//eta prime
1439           else if(momstatus == 11 || momstatus  == 12 ) fhMCEtaPtOrigin->Fill(pt,3.5);//resonances
1440           else fhMCEtaPtOrigin->Fill(pt,4.5);//stable, conversions?
1441           //printf("Other Meson pdg %d, Mother %s, pdg %d, status %d\n",pdg, TDatabasePDG::Instance()->GetParticle(mompdg)->GetName(),mompdg, momstatus );
1442         }// eta mass region
1443       }
1444       else if(ancPDG==-2212){//AProton
1445         fhMCOrgMass[4]->Fill(pt,mass);
1446         fhMCOrgAsym[4]->Fill(pt,asym);
1447         fhMCOrgDeltaEta[4]->Fill(pt,deta);
1448         fhMCOrgDeltaPhi[4]->Fill(pt,dphi);
1449       }
1450       else if(ancPDG==-2112){//ANeutron
1451         fhMCOrgMass[5]->Fill(pt,mass);
1452         fhMCOrgAsym[5]->Fill(pt,asym);
1453         fhMCOrgDeltaEta[5]->Fill(pt,deta);
1454         fhMCOrgDeltaPhi[5]->Fill(pt,dphi);
1455       }
1456       else if(TMath::Abs(ancPDG)==13){//muons
1457         fhMCOrgMass[6]->Fill(pt,mass);
1458         fhMCOrgAsym[6]->Fill(pt,asym);
1459         fhMCOrgDeltaEta[6]->Fill(pt,deta);
1460         fhMCOrgDeltaPhi[6]->Fill(pt,dphi);
1461       }
1462       else if (TMath::Abs(ancPDG) > 100 && ancLabel > 7) {
1463         if(ancStatus==1){//Stable particles, converted? not decayed resonances
1464           fhMCOrgMass[6]->Fill(pt,mass);
1465           fhMCOrgAsym[6]->Fill(pt,asym);
1466           fhMCOrgDeltaEta[6]->Fill(pt,deta);
1467           fhMCOrgDeltaPhi[6]->Fill(pt,dphi);
1468         }
1469         else{//resonances and other decays, more hadron conversions?
1470           fhMCOrgMass[7]->Fill(pt,mass);
1471           fhMCOrgAsym[7]->Fill(pt,asym);
1472           fhMCOrgDeltaEta[7]->Fill(pt,deta);
1473           fhMCOrgDeltaPhi[7]->Fill(pt,dphi);
1474         }
1475       }
1476       else {//Partons, colliding protons, strings, intermediate corrections
1477         if(ancStatus==11 || ancStatus==12){//String fragmentation
1478           fhMCOrgMass[8]->Fill(pt,mass);
1479           fhMCOrgAsym[8]->Fill(pt,asym);
1480           fhMCOrgDeltaEta[8]->Fill(pt,deta);
1481           fhMCOrgDeltaPhi[8]->Fill(pt,dphi);
1482         }
1483         else if (ancStatus==21){
1484           if(ancLabel < 2) {//Colliding protons
1485             fhMCOrgMass[11]->Fill(pt,mass);
1486             fhMCOrgAsym[11]->Fill(pt,asym);
1487             fhMCOrgDeltaEta[11]->Fill(pt,deta);
1488             fhMCOrgDeltaPhi[11]->Fill(pt,dphi);
1489           }//colliding protons
1490           else if(ancLabel < 6){//partonic initial states interactions
1491             fhMCOrgMass[9]->Fill(pt,mass);
1492             fhMCOrgAsym[9]->Fill(pt,asym);
1493             fhMCOrgDeltaEta[9]->Fill(pt,deta);
1494             fhMCOrgDeltaPhi[9]->Fill(pt,dphi);
1495           }
1496           else if(ancLabel < 8){//Final state partons radiations?
1497             fhMCOrgMass[10]->Fill(pt,mass);
1498             fhMCOrgAsym[10]->Fill(pt,asym);
1499             fhMCOrgDeltaEta[10]->Fill(pt,deta);
1500             fhMCOrgDeltaPhi[10]->Fill(pt,dphi);
1501           }
1502           else {
1503             printf("AliAnaPi0::FillMCVersusRecDataHistograms() - Check ** Common ancestor label %d, pdg %d, name %s, status %d; \n",
1504                    ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus);
1505           }
1506         }//status 21
1507         else {
1508           printf("AliAnaPi0::FillMCVersusRecDataHistograms() - Check *** Common ancestor label %d, pdg %d, name %s, status %d; \n",
1509                  ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus);
1510         }
1511       }////Partons, colliding protons, strings, intermediate corrections
1512     }//ancLabel > -1
1513     else { //ancLabel <= -1
1514       //printf("Not related at all label = %d\n",ancLabel);
1515       fhMCOrgMass[12]->Fill(pt,mass);
1516       fhMCOrgAsym[12]->Fill(pt,asym);
1517       fhMCOrgDeltaEta[12]->Fill(pt,deta);
1518       fhMCOrgDeltaPhi[12]->Fill(pt,dphi);
1519     }
1520 }
1521
1522 //____________________________________________________________________________________________________________________________________________________
1523 void AliAnaPi0::MakeAnalysisFillHistograms()
1524 {
1525   //Process one event and extract photons from AOD branch
1526   // filled with AliAnaPhoton and fill histos with invariant mass
1527
1528   //In case of simulated data, fill acceptance histograms
1529   if(IsDataMC())FillAcceptanceHistograms();
1530
1531   //if (GetReader()->GetEventNumber()%10000 == 0)
1532   // printf("--- Event %d ---\n",GetReader()->GetEventNumber());
1533
1534   //Init some variables
1535 //Int_t   iRun     = (GetReader()->GetInputEvent())->GetRunNumber() ;
1536   Int_t   nPhot    = GetInputAODBranch()->GetEntriesFast() ;
1537   Int_t   nClus    = 0;
1538   Int_t   nCell    = 0;
1539   Float_t eClusTot = 0;
1540   Float_t eCellTot = 0;
1541   Float_t eDenClus = 0;
1542   Float_t eDenCell = 0;
1543 //  Int_t   ncomb    = 0;
1544 //  Float_t rtmp     = 0;
1545 //  Float_t rtmpw    = 0;
1546 //  Float_t rxz      = 0;
1547 //  Float_t rxzw     = 0;
1548 //  Float_t pos1[3];
1549 //  Float_t pos2[3];
1550 //  Float_t emax     = 0;
1551
1552   if(GetDebug() > 1)
1553     printf("AliAnaPi0::MakeAnalysisFillHistograms() - Photon entries %d\n", nPhot);
1554
1555   //If less than photon 2 entries in the list, skip this event
1556   if(nPhot < 2 ) return ;
1557
1558   // Count the number of clusters and cells, in case multiplicity bins dependent on such numbers
1559   // are requested
1560   if(fCalorimeter=="EMCAL"){
1561     nClus = GetEMCALClusters()  ->GetEntriesFast();
1562     nCell = GetEMCALCells()->GetNumberOfCells();
1563     for(Int_t icl=0; icl < nClus; icl++) {
1564       Float_t e1 = ((AliVCluster*)GetEMCALClusters()->At(icl))->E();
1565       eClusTot +=  e1;
1566 //      if(e1 > emax) emax = e1;
1567 //      ((AliVCluster*)GetEMCALClusters()->At(icl))->GetPosition(pos1);
1568 //      for(Int_t icl2=icl+1; icl2 < nClus; icl2++) {
1569 //        Float_t e2 = ((AliVCluster*)GetEMCALClusters()->At(icl2))->E();
1570 //        ((AliVCluster*)GetEMCALClusters()->At(icl2))->GetPosition(pos2);
1571 //        rtmp  =  TMath::Sqrt((pos1[0]-pos2[0])*(pos1[0]-pos2[0]) + (pos1[2]-pos2[2])*(pos1[2]-pos2[2]));
1572 //        rtmpw =  TMath::Sqrt((pos1[0]*e1-pos2[0]*e2)*(pos1[0]*e1-pos2[0]*e2) + (pos1[2]*e1-pos2[2]*e2)*(pos1[2]*e1-pos2[2]*e2))/(e1+e2);
1573 //        rxz  += rtmp;
1574 //        rxzw += rtmpw;
1575 //        ncomb++;
1576 //        fhClusterPairDist      ->Fill(rtmp);
1577 //        fhClusterPairDistWeight->Fill(rtmpw);
1578 //        //printf("Distance: %f; weighted  %f\n ",rtmp,rtmp/(e1+((AliVCluster*)GetEMCALClusters()->At(icl2))->E()));
1579 //
1580 //      }// second cluster loop
1581     }// first cluster
1582
1583     for(Int_t jce=0; jce < nCell; jce++) eCellTot +=  GetEMCALCells()->GetAmplitude(jce);
1584   }
1585   else {
1586     nClus = GetPHOSClusters()->GetEntriesFast();
1587     nCell = GetPHOSCells()   ->GetNumberOfCells();
1588     for(Int_t icl=0; icl < nClus; icl++) {
1589       Float_t e1 = ((AliVCluster*)GetPHOSClusters()->At(icl))->E();
1590       eClusTot +=  e1;
1591 //      ((AliVCluster*)GetPHOSClusters()->At(icl))->GetPosition(pos1);
1592 //      for(Int_t icl2=icl+1; icl2 < nClus; icl2++) {
1593 //        Float_t e2 = ((AliVCluster*)GetPHOSClusters()->At(icl2))->E();
1594 //        ((AliVCluster*)GetPHOSClusters()->At(icl2))->GetPosition(pos2);
1595 //        rtmp  = TMath::Sqrt((pos1[0]-pos2[0])*(pos1[0]-pos2[0]) + (pos1[2]-pos2[2])*(pos1[2]-pos2[2]));
1596 //        rtmpw =  TMath::Sqrt((pos1[0]*e1-pos2[0]*e2)*(pos1[0]*e1-pos2[0]*e2) + (pos1[2]*e1-pos2[2]*e2)*(pos1[2]*e1-pos2[2]*e2))/(e1+e2);
1597 //        rxz  += rtmp;
1598 //        rxzw += rtmpw;
1599 //        ncomb++;
1600 //        fhClusterPairDist      ->Fill(rtmp);
1601 //        fhClusterPairDistWeight->Fill(rtmpw);
1602 //      }// second cluster loop
1603     }// first cluster
1604     for(Int_t jce=0; jce < nCell; jce++) eCellTot +=  GetPHOSCells()->GetAmplitude(jce);
1605   }
1606   if(GetDebug() > 1)
1607     printf("AliAnaPi0::MakeAnalysisFillHistograms() - # Clusters %d, sum cluster E per SM %f,# Cells %d, sum cell E per SM %f\n", nClus,eClusTot,nCell,eCellTot);
1608
1609   //Fill histograms with "energy density", ncell and nclust will be > 0 since there are at least 2 "photons"
1610   eDenClus = eClusTot/nClus;
1611   eDenCell = eCellTot/nCell;
1612   fhEDensityCluster      ->Fill(eDenClus);
1613   fhEDensityCell         ->Fill(eDenCell);
1614   fhEDensityCellvsCluster->Fill(eDenClus, eDenCell);
1615   //Fill the average number of cells or clusters per SM
1616   eClusTot /=fNModules;
1617   eCellTot /=fNModules;
1618   fhAverTotECluster      ->Fill(eClusTot);
1619   fhAverTotECell         ->Fill(eCellTot);
1620   fhAverTotECellvsCluster->Fill(eClusTot, eCellTot);
1621   //printf("Average Cluster: E %f, density %f;  Average Cell E %f, density  %f\n ",eClusTot,eDenClus,eCellTot,eDenCell);
1622
1623 //  //Average weighted pair distance
1624 //  rxz  /= ncomb;
1625 //  rxzw /= ncomb;
1626 //
1627 //  fhAverClusterPairDist             ->Fill(rxz );
1628 //  fhAverClusterPairDistWeight       ->Fill(rxzw);
1629 //  fhAverClusterPairDistvsAverE      ->Fill(rxz ,eDenClus);
1630 //  fhAverClusterPairDistWeightvsAverE->Fill(rxzw,eDenClus);
1631 //  fhAverClusterPairDistvsN          ->Fill(rxz ,nClus);
1632 //  fhAverClusterPairDistWeightvsN    ->Fill(rxzw,nClus);
1633 //
1634 //  //emax
1635 //  fhMaxEvsClustEDen->Fill(emax,eDenClus);
1636 //  fhMaxEvsClustMult->Fill(emax,nPhot);
1637
1638   //printf("Average Distance: %f; weighted  %f\n ",rxz,rxzw);
1639
1640
1641   //Init variables
1642   Int_t module1         = -1;
1643   Int_t module2         = -1;
1644   Double_t vert[]       = {0.0, 0.0, 0.0} ; //vertex
1645   Int_t evtIndex1       = 0 ;
1646   Int_t currentEvtIndex = -1;
1647   Int_t curCentrBin     = 0 ;
1648   Int_t curRPBin        = 0 ;
1649   Int_t curZvertBin     = 0 ;
1650
1651   //---------------------------------
1652   //First loop on photons/clusters
1653   //---------------------------------
1654   for(Int_t i1=0; i1<nPhot-1; i1++){
1655     AliAODPWG4Particle * p1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i1)) ;
1656     //printf("AliAnaPi0::MakeAnalysisFillHistograms() : cluster1 id %d\n",p1->GetCaloLabel(0));
1657
1658     // get the event index in the mixed buffer where the photon comes from
1659     // in case of mixing with analysis frame, not own mixing
1660     evtIndex1 = GetEventIndex(p1, vert) ;
1661     //printf("charge = %d\n", track->Charge());
1662     if ( evtIndex1 == -1 )
1663       return ;
1664     if ( evtIndex1 == -2 )
1665       continue ;
1666
1667     //printf("z vertex %f < %f\n",vert[2],GetZvertexCut());
1668     if(TMath::Abs(vert[2]) > GetZvertexCut()) continue ;   //vertex cut
1669
1670
1671     //----------------------------------------------------------------------------
1672     // Get the multiplicity bin. Different cases: centrality (PbPb),
1673     // average cluster multiplicity, average cell multiplicity, track multiplicity
1674     // default is centrality bins
1675     //----------------------------------------------------------------------------
1676     if (evtIndex1 != currentEvtIndex) {
1677       if(fUseTrackMultBins){ // Track multiplicity bins
1678         //printf("track  mult %d\n",GetTrackMultiplicity());
1679         curCentrBin = (GetTrackMultiplicity()-1)/5;
1680         if(curCentrBin > fNCentrBin-1) curCentrBin=fNCentrBin-1;
1681         //printf("track mult bin %d\n",curCentrBin);
1682       }
1683       else if(fUsePhotonMultBins){ // Photon multiplicity bins
1684         //printf("photon  mult %d cluster mult %d\n",nPhot, nClus);
1685         curRPBin = nPhot-2;
1686         if(curRPBin > GetNRPBin() -1) curRPBin=GetNRPBin()-1;
1687         //printf("photon mult bin %d\n",curRPBin);
1688       }
1689       else if(fUseAverClusterEBins){ // Cluster average energy bins
1690         //Bins for pp, if needed can be done in a more general way
1691         curCentrBin = (Int_t) eClusTot/10 * fNCentrBin;
1692         if(curCentrBin > fNCentrBin-1) curCentrBin=fNCentrBin-1;
1693         //printf("cluster E average %f, bin %d \n",eClusTot,curCentrBin);
1694       }
1695       else if(fUseAverCellEBins){ // Cell average energy bins
1696         //Bins for pp, if needed can be done in a more general way
1697         curCentrBin = (Int_t) eCellTot/10*fNCentrBin;
1698         if(curCentrBin > fNCentrBin-1) curCentrBin=fNCentrBin-1;
1699         //printf("cell E average %f, bin %d \n",eCellTot,curCentrBin);
1700       }
1701       else if(fUseAverClusterEDenBins){ // Energy density bins
1702         //Bins for pp, if needed can be done in a more general way
1703         curCentrBin = (Int_t) eDenClus/10*fNCentrBin;
1704         if(curCentrBin > fNCentrBin-1) curCentrBin=fNCentrBin-1;
1705         //printf("cluster Eden average %f, bin %d \n",eDenClus,curCentrBin);
1706       }
1707 //      else if(fUseAverClusterPairRBins){ // Cluster average distance bins
1708 //        //Bins for pp, if needed can be done in a more general way
1709 //        curCentrBin = rxz/650*fNCentrBin;
1710 //        if(curCentrBin > fNCentrBin-1) curCentrBin=fNCentrBin-1;
1711 //        //printf("cluster pair R average %f, bin %d \n",rxz,curCentrBin);
1712 //      }
1713 //      else if(fUseAverClusterPairRWeightBins){ // Cluster average distance bins
1714 //        //Bins for pp, if needed can be done in a more general way
1715 //        curCentrBin = rxzw/350*fNCentrBin;
1716 //        if(curCentrBin > fNCentrBin-1) curCentrBin=fNCentrBin-1;
1717 //        //printf("cluster pair rW average %f, bin %d \n",rxzw,curCentrBin);
1718 //      }
1719 //      else if(fUseEMaxBins){ // Cluster average distance bins
1720 //        //Bins for pp, if needed can be done in a more general way
1721 //        curCentrBin = emax/20*fNCentrBin;
1722 //        if(curCentrBin > fNCentrBin-1) curCentrBin=fNCentrBin-1;
1723 //        //printf("cluster pair rW average %f, bin %d \n",rxzw,curCentrBin);
1724 //      }
1725       else { //Event centrality
1726           curCentrBin = GetEventCentrality();
1727         //printf("curCentrBin %d\n",curCentrBin);
1728       }
1729
1730       if (curCentrBin < 0 || curCentrBin >= fNCentrBin){
1731         if(GetDebug() > 0)
1732           printf("AliAnaPi0::MakeAnalysisFillHistograms() - Centrality bin <%d> not expected, n bins <%d> , return\n",curCentrBin,fNCentrBin);
1733         return;
1734       }
1735
1736       //Reaction plane bin
1737       curRPBin    = 0 ;
1738
1739       //Get vertex z bin
1740       curZvertBin = (Int_t)(0.5*GetNZvertBin()*(vert[2]+GetZvertexCut())/GetZvertexCut()) ;
1741
1742       //Fill event bin info
1743       fhEvents    ->Fill(curCentrBin+0.5,curZvertBin+0.5,curRPBin+0.5) ;
1744       fhCentrality->Fill(curCentrBin);
1745       currentEvtIndex = evtIndex1 ;
1746       if(GetDebug() > 1)
1747         printf("AliAnaPi0::MakeAnalysisFillHistograms() - Centrality %d, Vertex Bin %d, RP bin %d \n",curCentrBin,curRPBin,curZvertBin);
1748     }
1749
1750     //printf("AliAnaPi0::MakeAnalysisFillHistograms(): Photon 1 Evt %d  Vertex : %f,%f,%f\n",evtIndex1, GetVertex(evtIndex1)[0] ,GetVertex(evtIndex1)[1],GetVertex(evtIndex1)[2]);
1751
1752     //Get the momentum of this cluster
1753     TLorentzVector photon1(p1->Px(),p1->Py(),p1->Pz(),p1->E());
1754
1755     //Get (Super)Module number of this cluster
1756     module1 = GetModuleNumber(p1);
1757
1758     //---------------------------------
1759     //Second loop on photons/clusters
1760     //---------------------------------
1761     for(Int_t i2=i1+1; i2<nPhot; i2++){
1762       AliAODPWG4Particle * p2 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i2)) ;
1763
1764       //In case of mixing frame, check we are not in the same event as the first cluster
1765       Int_t evtIndex2 = GetEventIndex(p2, vert) ;
1766       if ( evtIndex2 == -1 )
1767         return ;
1768       if ( evtIndex2 == -2 )
1769         continue ;
1770       if (GetMixedEvent() && (evtIndex1 == evtIndex2))
1771         continue ;
1772
1773       //printf("AliAnaPi0::MakeAnalysisFillHistograms(): Photon 2 Evt %d  Vertex : %f,%f,%f\n",evtIndex2, GetVertex(evtIndex2)[0] ,GetVertex(evtIndex2)[1],GetVertex(evtIndex2)[2]);
1774
1775       //Get the momentum of this cluster
1776       TLorentzVector photon2(p2->Px(),p2->Py(),p2->Pz(),p2->E());
1777       //Get module number
1778       module2       = GetModuleNumber(p2);
1779
1780       //---------------------------------
1781       // Get pair kinematics
1782       //---------------------------------
1783       Double_t m    = (photon1 + photon2).M() ;
1784       Double_t pt   = (photon1 + photon2).Pt();
1785       Double_t deta = photon1.Eta() - photon2.Eta();
1786       Double_t dphi = photon1.Phi() - photon2.Phi();
1787       Double_t a    = TMath::Abs(p1->E()-p2->E())/(p1->E()+p2->E()) ;
1788
1789       if(GetDebug() > 2)
1790         printf(" E: photon1 %f, photon2 %f; Pair: pT %f, mass %f, a %f\n", p1->E(), p2->E(), (photon1 + photon2).E(),m,a);
1791
1792       //--------------------------------
1793       // Opening angle selection
1794       //--------------------------------
1795       //Check if opening angle is too large or too small compared to what is expected
1796       Double_t angle   = photon1.Angle(photon2.Vect());
1797       if(fUseAngleEDepCut && !GetNeutralMesonSelection()->IsAngleInWindow((photon1+photon2).E(),angle+0.05)) {
1798         if(GetDebug() > 2)
1799           printf("AliAnaPi0::MakeAnalysisFillHistograms() -Real pair angle %f not in E %f window\n",angle, (photon1+photon2).E());
1800         continue;
1801       }
1802
1803       if(fUseAngleCut && (angle < fAngleCut || angle > fAngleMaxCut)) {
1804         if(GetDebug() > 2)
1805           printf("AliAnaPi0::MakeAnalysisFillHistograms() - Real pair cut %f < angle %f < cut %f\n",fAngleCut, angle, fAngleMaxCut);
1806         continue;
1807       }
1808
1809       //-------------------------------------------------------------------------------------------------
1810       //Fill module dependent histograms, put a cut on assymmetry on the first available cut in the array
1811       //-------------------------------------------------------------------------------------------------
1812       if(a < fAsymCuts[0]){
1813         if(module1==module2 && module1 >=0 && module1<fNModules)
1814           fhReMod[module1]->Fill(pt,m) ;
1815
1816         if(fCalorimeter=="EMCAL"){
1817
1818           // Same sector
1819           Int_t j=0;
1820           for(Int_t i = 0; i < fNModules/2; i++){
1821             j=2*i;
1822             if((module1==j && module2==j+1) || (module1==j+1 && module2==j)) fhReSameSectorEMCALMod[i]->Fill(pt,m) ;
1823           }
1824
1825           // Same side
1826           for(Int_t i = 0; i < fNModules-2; i++){
1827             if((module1==i && module2==i+2) || (module1==i+2 && module2==i)) fhReSameSideEMCALMod[i]->Fill(pt,m);
1828           }
1829         }//EMCAL
1830         else {//PHOS
1831           if((module1==0 && module2==1) || (module1==1 && module2==0)) fhReDiffPHOSMod[0]->Fill(pt,m) ;
1832           if((module1==0 && module2==2) || (module1==2 && module2==0)) fhReDiffPHOSMod[1]->Fill(pt,m) ;
1833           if((module1==1 && module2==2) || (module1==2 && module2==1)) fhReDiffPHOSMod[2]->Fill(pt,m) ;
1834         }//PHOS
1835       }
1836
1837       //In case we want only pairs in same (super) module, check their origin.
1838       Bool_t ok = kTRUE;
1839       if(fSameSM && module1!=module2) ok=kFALSE;
1840       if(ok){
1841
1842         //Check if one of the clusters comes from a conversion
1843         if     (p1->IsTagged() && p2->IsTagged()) fhReConv2->Fill(pt,m);
1844         else if(p1->IsTagged() || p2->IsTagged()) fhReConv ->Fill(pt,m);
1845
1846         //Fill histograms for different bad channel distance, centrality, assymmetry cut and pid bit
1847         for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
1848           if((p1->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton)) && (p2->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton))){
1849             for(Int_t iasym=0; iasym < fNAsymCuts; iasym++){
1850               if(a < fAsymCuts[iasym]){
1851                 Int_t index = ((curCentrBin*fNPIDBits)+ipid)*fNAsymCuts + iasym;
1852                 //printf("index %d :(cen %d * nPID %d + ipid %d)*nasym %d + iasym %d\n",index,curCentrBin,fNPIDBits,ipid,fNAsymCuts,iasym);
1853                 fhRe1     [index]->Fill(pt,m);
1854                 if(fMakeInvPtPlots)fhReInvPt1[index]->Fill(pt,m,1./pt) ;
1855                 if(fFillBadDistHisto){
1856                   if(p1->DistToBad()>0 && p2->DistToBad()>0){
1857                     fhRe2     [index]->Fill(pt,m) ;
1858                     if(fMakeInvPtPlots)fhReInvPt2[index]->Fill(pt,m,1./pt) ;
1859                     if(p1->DistToBad()>1 && p2->DistToBad()>1){
1860                       fhRe3     [index]->Fill(pt,m) ;
1861                       if(fMakeInvPtPlots)fhReInvPt3[index]->Fill(pt,m,1./pt) ;
1862                     }// bad 3
1863                   }// bad2
1864                 }// Fill bad dist histos
1865               }//assymetry cut
1866             }// asymmetry cut loop
1867           }// bad 1
1868         }// pid bit loop
1869
1870         //Fill histograms with opening angle
1871         fhRealOpeningAngle   ->Fill(pt,angle);
1872         fhRealCosOpeningAngle->Fill(pt,TMath::Cos(angle));
1873
1874         //Fill histograms with pair assymmetry
1875         fhRePtAsym->Fill(pt,a);
1876         if(m > 0.10 && m < 0.17) fhRePtAsymPi0->Fill(pt,a);
1877         if(m > 0.45 && m < 0.65) fhRePtAsymEta->Fill(pt,a);
1878
1879         //-------------------------------------------------------
1880         //Get the number of cells needed for multi cut analysis.
1881         //-------------------------------------------------------
1882         Int_t ncell1 = 0;
1883         Int_t ncell2 = 0;
1884         if(fMultiCutAna || (IsDataMC() && fMultiCutAnaSim)){
1885
1886           AliVEvent * event = GetReader()->GetInputEvent();
1887           if(event){
1888             for(Int_t iclus = 0; iclus < event->GetNumberOfCaloClusters(); iclus++){
1889               AliVCluster *cluster = event->GetCaloCluster(iclus);
1890
1891               Bool_t is = kFALSE;
1892               if     (fCalorimeter == "EMCAL" && GetReader()->IsEMCALCluster(cluster)) is = kTRUE;
1893               else if(fCalorimeter == "PHOS"  && GetReader()->IsPHOSCluster (cluster)) is = kTRUE;
1894
1895               if(is){
1896                 if      (p1->GetCaloLabel(0) == cluster->GetID()) ncell1 = cluster->GetNCells();
1897                 else if (p2->GetCaloLabel(0) == cluster->GetID()) ncell2 = cluster->GetNCells();
1898               } // PHOS or EMCAL cluster as requested in analysis
1899
1900               if(ncell2 > 0 &&  ncell1 > 0) break; // No need to continue the iteration
1901
1902             }
1903             //printf("e 1: %2.2f, e 2: %2.2f, ncells: n1 %d, n2 %d\n", p1->E(), p2->E(),ncell1,ncell2);
1904           }
1905         }
1906
1907         //---------
1908         // MC data
1909         //---------
1910         //Do some MC checks on the origin of the pair, is there any common ancestor and if there is one, who?
1911         if(IsDataMC()) FillMCVersusRecDataHistograms(p1->GetLabel(), p2->GetLabel(),p1->Pt(), p2->Pt(),ncell1, ncell2, m, pt, a,deta, dphi);
1912
1913         //-----------------------
1914         //Multi cuts analysis
1915         //-----------------------
1916         if(fMultiCutAna){
1917           //Histograms for different PID bits selection
1918           for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
1919
1920             if(p1->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton)    &&
1921                p2->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton))   fhRePIDBits[ipid]->Fill(pt,m) ;
1922
1923             //printf("ipt %d, ipid%d, name %s\n",ipt, ipid, fhRePtPIDCuts[ipt*fNPIDBitsBits+ipid]->GetName());
1924           } // pid bit cut loop
1925
1926           //Several pt,ncell and asymmetry cuts
1927           for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
1928             for(Int_t icell=0; icell<fNCellNCuts; icell++){
1929               for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
1930                 Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
1931                 if(p1->Pt() >   fPtCuts[ipt]      && p2->Pt() > fPtCuts[ipt]        &&
1932                    a        <   fAsymCuts[iasym]                                    &&
1933                    ncell1   >=  fCellNCuts[icell] && ncell2   >= fCellNCuts[icell]){
1934                     fhRePtNCellAsymCuts[index]->Fill(pt,m) ;
1935                   //printf("ipt %d, icell%d, iasym %d, name %s\n",ipt, icell, iasym,  fhRePtNCellAsymCuts[((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym]->GetName());
1936                   if(module1==module2){
1937                     if     (module1==0)  fhRePtNCellAsymCutsSM0[index]->Fill(pt,m) ;
1938                     else if(module1==1)  fhRePtNCellAsymCutsSM1[index]->Fill(pt,m) ;
1939                     else if(module1==2)  fhRePtNCellAsymCutsSM2[index]->Fill(pt,m) ;
1940                     else if(module1==3)  fhRePtNCellAsymCutsSM3[index]->Fill(pt,m) ;
1941                     else printf("AliAnaPi0::FillHistograms() - WRONG SM NUMBER\n");
1942                   }
1943                 }
1944               }// pid bit cut loop
1945             }// icell loop
1946           }// pt cut loop
1947           for(Int_t iasym = 0; iasym < fNAsymCuts; iasym++){
1948             if(a < fAsymCuts[iasym])fhRePtMult[iasym]->Fill(pt,GetTrackMultiplicity(),m) ;
1949           }
1950         }// multiple cuts analysis
1951       }// ok if same sm
1952     }// second same event particle
1953   }// first cluster
1954
1955   //-------------------------------------------------------------
1956   // Mixing
1957   //-------------------------------------------------------------
1958   if(fDoOwnMix){
1959     //printf("Cen bin %d, RP bin %d, e aver %f, mult %d\n",curCentrBin,curRPBin, eClusTot, nPhot);
1960     //Recover events in with same characteristics as the current event
1961     TList * evMixList=fEventsList[curCentrBin*GetNZvertBin()*GetNRPBin()+curZvertBin*GetNRPBin()+curRPBin] ;
1962     Int_t nMixed = evMixList->GetSize() ;
1963     for(Int_t ii=0; ii<nMixed; ii++){
1964       TClonesArray* ev2= (TClonesArray*) (evMixList->At(ii));
1965       Int_t nPhot2=ev2->GetEntriesFast() ;
1966       Double_t m = -999;
1967       if(GetDebug() > 1)
1968         printf("AliAnaPi0::MakeAnalysisFillHistograms() - Mixed event %d photon entries %d, centrality bin %d\n", ii, nPhot2, curCentrBin);
1969
1970       //---------------------------------
1971       //First loop on photons/clusters
1972       //---------------------------------
1973       for(Int_t i1=0; i1<nPhot; i1++){
1974         AliAODPWG4Particle * p1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i1)) ;
1975         if(fSameSM && GetModuleNumber(p1)!=module1) continue;
1976
1977         //Get kinematics of cluster and (super) module of this cluster
1978         TLorentzVector photon1(p1->Px(),p1->Py(),p1->Pz(),p1->E());
1979         module1 = GetModuleNumber(p1);
1980
1981         //---------------------------------
1982         //First loop on photons/clusters
1983         //---------------------------------
1984         for(Int_t i2=0; i2<nPhot2; i2++){
1985           AliAODPWG4Particle * p2 = (AliAODPWG4Particle*) (ev2->At(i2)) ;
1986
1987           //Get kinematics of second cluster and calculate those of the pair
1988           TLorentzVector photon2(p2->Px(),p2->Py(),p2->Pz(),p2->E());
1989           m           = (photon1+photon2).M() ;
1990           Double_t pt = (photon1 + photon2).Pt();
1991           Double_t a  = TMath::Abs(p1->E()-p2->E())/(p1->E()+p2->E()) ;
1992
1993           //Check if opening angle is too large or too small compared to what is expected
1994           Double_t angle   = photon1.Angle(photon2.Vect());
1995           if(fUseAngleEDepCut && !GetNeutralMesonSelection()->IsAngleInWindow((photon1+photon2).E(),angle+0.05)){
1996             if(GetDebug() > 2)
1997               printf("AliAnaPi0::MakeAnalysisFillHistograms() -Mix pair angle %f not in E %f window\n",angle, (photon1+photon2).E());
1998             continue;
1999           }
2000           if(fUseAngleCut && (angle < fAngleCut || angle > fAngleMaxCut)) {
2001             if(GetDebug() > 2)
2002               printf("AliAnaPi0::MakeAnalysisFillHistograms() -Mix pair angle %f < cut %f\n",angle,fAngleCut);
2003             continue;
2004
2005           }
2006
2007           if(GetDebug() > 2)
2008             printf("AliAnaPi0::MakeAnalysisFillHistograms() - Mixed Event: pT: photon1 %2.2f, photon2 %2.2f; Pair: pT %2.2f, mass %2.3f, a %f2.3\n",
2009                    p1->Pt(), p2->Pt(), pt,m,a);
2010
2011           //In case we want only pairs in same (super) module, check their origin.
2012           module2 = GetModuleNumber(p2);
2013
2014           //-------------------------------------------------------------------------------------------------
2015           //Fill module dependent histograms, put a cut on assymmetry on the first available cut in the array
2016           //-------------------------------------------------------------------------------------------------
2017           if(a < fAsymCuts[0]){
2018             if(module1==module2 && module1 >=0 && module1<fNModules)
2019               fhMiMod[module1]->Fill(pt,m) ;
2020
2021             if(fCalorimeter=="EMCAL"){
2022
2023               // Same sector
2024               Int_t j=0;
2025               for(Int_t i = 0; i < fNModules/2; i++){
2026                 j=2*i;
2027                 if((module1==j && module2==j+1) || (module1==j+1 && module2==j)) fhMiSameSectorEMCALMod[i]->Fill(pt,m) ;
2028               }
2029
2030               // Same side
2031               for(Int_t i = 0; i < fNModules-2; i++){
2032                 if((module1==i && module2==i+2) || (module1==i+2 && module2==i)) fhMiSameSideEMCALMod[i]->Fill(pt,m);
2033               }
2034             }//EMCAL
2035             else {//PHOS
2036               if((module1==0 && module2==1) || (module1==1 && module2==0)) fhMiDiffPHOSMod[0]->Fill(pt,m) ;
2037               if((module1==0 && module2==2) || (module1==2 && module2==0)) fhMiDiffPHOSMod[1]->Fill(pt,m) ;
2038               if((module1==1 && module2==2) || (module1==2 && module2==1)) fhMiDiffPHOSMod[2]->Fill(pt,m) ;
2039             }//PHOS
2040
2041
2042           }
2043
2044           Bool_t ok = kTRUE;
2045           if(fSameSM && module1!=module2) ok=kFALSE;
2046           if(ok){
2047
2048             //Check if one of the clusters comes from a conversion
2049             if     (p1->IsTagged() && p2->IsTagged()) fhMiConv2->Fill(pt,m);
2050             else if(p1->IsTagged() || p2->IsTagged()) fhMiConv ->Fill(pt,m);
2051
2052             //Fill histograms for different bad channel distance, centrality, assymmetry cut and pid bit
2053             for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
2054               if((p1->IsPIDOK(ipid,AliCaloPID::kPhoton)) && (p2->IsPIDOK(ipid,AliCaloPID::kPhoton))){
2055                 for(Int_t iasym=0; iasym < fNAsymCuts; iasym++){
2056                   if(a < fAsymCuts[iasym]){
2057                     Int_t index = ((curCentrBin*fNPIDBits)+ipid)*fNAsymCuts + iasym;
2058                     fhMi1     [index]->Fill(pt,m) ;
2059                     if(fMakeInvPtPlots)fhMiInvPt1[index]->Fill(pt,m,1./pt) ;
2060                     if(fFillBadDistHisto){
2061                       if(p1->DistToBad()>0 && p2->DistToBad()>0){
2062                         fhMi2     [index]->Fill(pt,m) ;
2063                         if(fMakeInvPtPlots)fhMiInvPt2[index]->Fill(pt,m,1./pt) ;
2064                         if(p1->DistToBad()>1 && p2->DistToBad()>1){
2065                           fhMi3     [index]->Fill(pt,m) ;
2066                           if(fMakeInvPtPlots)fhMiInvPt3[index]->Fill(pt,m,1./pt) ;
2067                         }
2068                       }
2069                     }// Fill bad dist histo
2070                   }//Asymmetry cut
2071                 }// Asymmetry loop
2072               }//PID cut
2073             }//loop for histograms
2074
2075             //-----------------------
2076             //Multi cuts analysis
2077             //-----------------------
2078             if(fMultiCutAna){
2079               //Several pt,ncell and asymmetry cuts
2080               for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
2081                 for(Int_t icell=0; icell<fNCellNCuts; icell++){
2082                   for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
2083                     Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
2084                     if(p1->Pt() >   fPtCuts[ipt]      && p2->Pt() > fPtCuts[ipt]        &&
2085                        a        <   fAsymCuts[iasym]                                    &&
2086                        p1->GetBtag() >=  fCellNCuts[icell] && p2->GetBtag() >= fCellNCuts[icell]){
2087                       fhMiPtNCellAsymCuts[index]->Fill(pt,m) ;
2088                       //printf("ipt %d, icell%d, iasym %d, name %s\n",ipt, icell, iasym,  fhRePtNCellAsymCuts[((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym]->GetName());
2089                     }
2090                   }// pid bit cut loop
2091                 }// icell loop
2092               }// pt cut loop
2093             } // Multi cut ana
2094
2095             //Fill histograms with opening angle
2096             fhMixedOpeningAngle   ->Fill(pt,angle);
2097             fhMixedCosOpeningAngle->Fill(pt,TMath::Cos(angle));
2098           }//ok
2099         }// second cluster loop
2100       }//first cluster loop
2101     }//loop on mixed events
2102
2103     //--------------------------------------------------------
2104     //Add the current event to the list of events for mixing
2105     //--------------------------------------------------------
2106     TClonesArray *currentEvent = new TClonesArray(*GetInputAODBranch());
2107     //Add current event to buffer and Remove redundant events
2108     if(currentEvent->GetEntriesFast()>0){
2109       evMixList->AddFirst(currentEvent) ;
2110       currentEvent=0 ; //Now list of particles belongs to buffer and it will be deleted with buffer
2111       if(evMixList->GetSize()>=fNmaxMixEv)
2112       {
2113         TClonesArray * tmp = (TClonesArray*) (evMixList->Last()) ;
2114         evMixList->RemoveLast() ;
2115         delete tmp ;
2116       }
2117     }
2118     else{ //empty event
2119       delete currentEvent ;
2120       currentEvent=0 ;
2121     }
2122   }// DoOwnMix
2123
2124 }
2125
2126 //________________________________________________________________________
2127 void AliAnaPi0::ReadHistograms(TList* outputList)
2128 {
2129   // Needed when Terminate is executed in distributed environment
2130   // Refill analysis histograms of this class with corresponding histograms in output list.
2131
2132   // Histograms of this analsys are kept in the same list as other analysis, recover the position of
2133   // the first one and then add the next.
2134   Int_t index = outputList->IndexOf(outputList->FindObject(GetAddedHistogramsStringToName()+"hRe_cen0_pid0_dist1"));
2135
2136   if(!fhRe1) fhRe1 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
2137   if(!fhRe2) fhRe2 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
2138   if(!fhRe3) fhRe3 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
2139   if(!fhMi1) fhMi1 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
2140   if(!fhMi2) fhMi2 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
2141   if(!fhMi3) fhMi3 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
2142   if(!fhReInvPt1) fhReInvPt1  = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
2143   if(!fhReInvPt2) fhReInvPt2  = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
2144   if(!fhReInvPt3) fhReInvPt3  = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
2145   if(!fhMiInvPt1) fhMiInvPt1  = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
2146   if(!fhMiInvPt2) fhMiInvPt2  = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
2147   if(!fhMiInvPt3) fhMiInvPt3  = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
2148   if(!fhReMod)    fhReMod     = new TH2D*[fNModules]   ;
2149   if(!fhReDiffPHOSMod)       fhReDiffPHOSMod        = new TH2D*[fNModules] ;
2150   if(!fhReSameSectorEMCALMod)fhReSameSectorEMCALMod = new TH2D*[fNModules/2] ;
2151   if(!fhReSameSideEMCALMod)  fhReSameSideEMCALMod   = new TH2D*[fNModules-2] ;
2152   if(!fhMiMod)    fhMiMod     = new TH2D*[fNModules]   ;
2153   if(!fhMiDiffPHOSMod)       fhMiDiffPHOSMod        = new TH2D*[fNModules] ;
2154   if(!fhMiSameSectorEMCALMod)fhMiSameSectorEMCALMod = new TH2D*[fNModules/2] ;
2155   if(!fhMiSameSideEMCALMod)  fhMiSameSideEMCALMod   = new TH2D*[fNModules-2] ;
2156
2157   fhReConv  = (TH2D*) outputList->At(index++);
2158   fhMiConv  = (TH2D*) outputList->At(index++);
2159   fhReConv2 = (TH2D*) outputList->At(index++);
2160   fhMiConv2 = (TH2D*) outputList->At(index++);
2161
2162   for(Int_t ic=0; ic<fNCentrBin; ic++){
2163     for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
2164       for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
2165         Int_t ihisto = ((ic*fNPIDBits)+ipid)*fNAsymCuts + iasym;
2166
2167         fhRe1[ihisto] = (TH2D*) outputList->At(index++);
2168         fhRe2[ihisto] = (TH2D*) outputList->At(index++);
2169         fhRe3[ihisto] = (TH2D*) outputList->At(index++);
2170
2171         fhReInvPt1[ihisto] = (TH2D*) outputList->At(index++);
2172         fhReInvPt2[ihisto] = (TH2D*) outputList->At(index++);
2173         fhReInvPt3[ihisto] = (TH2D*) outputList->At(index++);
2174
2175         if(fDoOwnMix){
2176           fhMi1[ihisto] = (TH2D*) outputList->At(index++);
2177           fhMi2[ihisto] = (TH2D*) outputList->At(index++);
2178           fhMi3[ihisto] = (TH2D*) outputList->At(index++);
2179
2180           fhMiInvPt1[ihisto] = (TH2D*) outputList->At(index++);
2181           fhMiInvPt2[ihisto] = (TH2D*) outputList->At(index++);
2182           fhMiInvPt3[ihisto] = (TH2D*) outputList->At(index++);
2183         }//Own mix
2184       }//asymmetry loop
2185     }// pid loop
2186   }// centrality loop
2187
2188   fhRePtAsym    = (TH2D*)outputList->At(index++);
2189   fhRePtAsymPi0 = (TH2D*)outputList->At(index++);
2190   fhRePtAsymEta = (TH2D*)outputList->At(index++);
2191
2192   if(fMultiCutAna){
2193
2194     if(!fhRePtNCellAsymCuts) fhRePtNCellAsymCuts = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
2195     if(!fhRePIDBits)         fhRePIDBits         = new TH2D*[fNPIDBits];
2196
2197     for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
2198       fhRePIDBits[ipid] = (TH2D*) outputList->At(index++);
2199     }// ipid loop
2200
2201     for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
2202       for(Int_t icell=0; icell<fNCellNCuts; icell++){
2203         for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
2204           fhRePtNCellAsymCuts[((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym] = (TH2D*) outputList->At(index++);
2205         }// iasym
2206       }// icell loop
2207     }// ipt loop
2208
2209     if(!fhRePtMult) fhRePtMult  = new TH3D*[fNAsymCuts]  ;
2210     for(Int_t iasym = 0; iasym < fNAsymCuts; iasym++)
2211       fhRePtMult[iasym] = (TH3D*) outputList->At(index++);
2212   }// multi cut analysis
2213
2214   fhEvents     = (TH3D *) outputList->At(index++);
2215   fhCentrality = (TH1D *) outputList->At(index++);
2216
2217   fhRealOpeningAngle     = (TH2D*)  outputList->At(index++);
2218   fhRealCosOpeningAngle  = (TH2D*)  outputList->At(index++);
2219   if(fDoOwnMix){
2220     fhMixedOpeningAngle     = (TH2D*)  outputList->At(index++);
2221     fhMixedCosOpeningAngle  = (TH2D*)  outputList->At(index++);
2222   }
2223
2224   //Histograms filled only if MC data is requested
2225   if(IsDataMC() || (GetReader()->GetDataType() == AliCaloTrackReader::kMC) ){
2226     fhPrimPi0Pt     = (TH1D*)  outputList->At(index++);
2227     fhPrimPi0AccPt  = (TH1D*)  outputList->At(index++);
2228     fhPrimPi0Y      = (TH2D*)  outputList->At(index++);
2229     fhPrimPi0AccY   = (TH2D*)  outputList->At(index++);
2230     fhPrimPi0Phi    = (TH2D*)  outputList->At(index++);
2231     fhPrimPi0AccPhi = (TH2D*)  outputList->At(index++);
2232     fhPrimEtaPt     = (TH1D*)  outputList->At(index++);
2233     fhPrimEtaAccPt  = (TH1D*)  outputList->At(index++);
2234     fhPrimEtaY      = (TH2D*)  outputList->At(index++);
2235     fhPrimEtaAccY   = (TH2D*)  outputList->At(index++);
2236     fhPrimEtaPhi    = (TH2D*)  outputList->At(index++);
2237     fhPrimEtaAccPhi = (TH2D*)  outputList->At(index++);
2238     for(Int_t i = 0; i<13; i++){
2239       fhMCOrgMass[i]     = (TH2D*)  outputList->At(index++);
2240       fhMCOrgAsym[i]     = (TH2D*)  outputList->At(index++);
2241       fhMCOrgDeltaEta[i] = (TH2D*)  outputList->At(index++);
2242       fhMCOrgDeltaPhi[i] = (TH2D*)  outputList->At(index++);
2243     }
2244
2245     if(fMultiCutAnaSim){
2246       fhMCPi0MassPtTrue  = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
2247       fhMCPi0MassPtRec   = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
2248       fhMCPi0PtTruePtRec = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
2249       fhMCEtaMassPtTrue  = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
2250       fhMCEtaMassPtRec   = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
2251       fhMCEtaPtTruePtRec = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
2252       for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
2253         for(Int_t icell=0; icell<fNCellNCuts; icell++){
2254           for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
2255             Int_t in = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
2256             fhMCPi0MassPtTrue[in]  = (TH2D*)  outputList->At(index++);
2257             fhMCPi0PtTruePtRec[in] = (TH2D*)  outputList->At(index++);
2258             fhMCEtaMassPtTrue[in]  = (TH2D*)  outputList->At(index++);
2259             fhMCEtaPtTruePtRec[in] = (TH2D*)  outputList->At(index++);
2260           }
2261         }
2262       }
2263     }
2264     else{
2265       fhMCPi0MassPtTrue  = new TH2D*[1];
2266       fhMCPi0PtTruePtRec = new TH2D*[1];
2267       fhMCEtaMassPtTrue  = new TH2D*[1];
2268       fhMCEtaPtTruePtRec = new TH2D*[1];
2269
2270       fhMCPi0MassPtTrue[0]  = (TH2D*)  outputList->At(index++);
2271       fhMCPi0PtTruePtRec[0] = (TH2D*)  outputList->At(index++);
2272       fhMCEtaMassPtTrue[0]  = (TH2D*)  outputList->At(index++);
2273       fhMCEtaPtTruePtRec[0] = (TH2D*)  outputList->At(index++);
2274     }
2275   }
2276
2277   for(Int_t imod=0; imod < fNModules; imod++){
2278     fhReMod[imod]                = (TH2D*) outputList->At(index++);
2279     if(fCalorimeter=="EMCAL"){
2280       if(imod < fNModules/2) fhReSameSectorEMCALMod[imod] = (TH2D*) outputList->At(index++);
2281       if(imod < fNModules-2) fhReSameSideEMCALMod[imod]   = (TH2D*) outputList->At(index++);
2282     }
2283     else     fhReDiffPHOSMod[imod]        = (TH2D*) outputList->At(index++);
2284
2285     if(fDoOwnMix){
2286       fhMiMod[imod]         = (TH2D*) outputList->At(index++);
2287       if(fCalorimeter=="EMCAL"){
2288         if(imod < fNModules/2) fhMiSameSectorEMCALMod[imod] = (TH2D*) outputList->At(index++);
2289         if(imod < fNModules-2) fhMiSameSideEMCALMod[imod]   = (TH2D*) outputList->At(index++);
2290       }
2291       else     fhMiDiffPHOSMod[imod]        = (TH2D*) outputList->At(index++);
2292     }
2293   }
2294
2295 }
2296
2297
2298 //____________________________________________________________________________________________________________________________________________________
2299 void AliAnaPi0::Terminate(TList* outputList)
2300 {
2301   //Do some calculations and plots from the final histograms.
2302
2303   printf(" *** %s Terminate:\n", GetName()) ;
2304
2305   //Recover histograms from output histograms list, needed for distributed analysis.
2306   ReadHistograms(outputList);
2307
2308   if (!fhRe1) {
2309     printf("AliAnaPi0::Terminate() - Error: Remote output histograms not imported in AliAnaPi0 object");
2310     return;
2311   }
2312
2313   printf("AliAnaPi0::Terminate()         Mgg Real        : %5.3f , RMS : %5.3f \n", fhRe1[0]->GetMean(),   fhRe1[0]->GetRMS() ) ;
2314
2315   const Int_t buffersize = 255;
2316
2317   char nameIM[buffersize];
2318   snprintf(nameIM, buffersize,"AliAnaPi0_%s_cPt",fCalorimeter.Data());
2319   TCanvas  * cIM = new TCanvas(nameIM, "", 400, 10, 600, 700) ;
2320   cIM->Divide(2, 2);
2321
2322   cIM->cd(1) ;
2323   //gPad->SetLogy();
2324   TH1D * hIMAllPt = (TH1D*) fhRe1[0]->ProjectionY(Form("IMPtAll_%s",fCalorimeter.Data()));
2325   hIMAllPt->SetLineColor(2);
2326   hIMAllPt->SetTitle("No cut on  p_{T, #gamma#gamma} ");
2327   hIMAllPt->Draw();
2328
2329   cIM->cd(2) ;
2330   TH1D * hIMPt5 = (TH1D*) fhRe1[0]->ProjectionY(Form("IMPt0-5_%s",fCalorimeter.Data()),0, fhRe1[0]->GetXaxis()->FindBin(5.));
2331 //  hRe1Pt5->GetXaxis()->SetRangeUser(0,5);
2332 //  TH1D * hIMPt5 = (TH1D*) hRe1Pt5->Project3D(Form("IMPt5_%s_pz",fCalorimeter.Data()));
2333   hIMPt5->SetLineColor(2);
2334   hIMPt5->SetTitle("0 < p_{T, #gamma#gamma} < 5 GeV/c");
2335   hIMPt5->Draw();
2336
2337   cIM->cd(3) ;
2338   TH1D * hIMPt10 =  (TH1D*) fhRe1[0]->ProjectionY(Form("IMPt5-10_%s",fCalorimeter.Data()), fhRe1[0]->GetXaxis()->FindBin(5.),fhRe1[0]->GetXaxis()->FindBin(10.));
2339 //  hRe1Pt10->GetXaxis()->SetRangeUser(5,10);
2340 //  TH1D * hIMPt10 = (TH1D*) hRe1Pt10->Project3D(Form("IMPt10_%s_pz",fCalorimeter.Data()));
2341   hIMPt10->SetLineColor(2);
2342   hIMPt10->SetTitle("5 < p_{T, #gamma#gamma} < 10 GeV/c");
2343   hIMPt10->Draw();
2344
2345   cIM->cd(4) ;
2346   TH1D * hIMPt20 =  (TH1D*) fhRe1[0]->ProjectionY(Form("IMPt10-20_%s",fCalorimeter.Data()), fhRe1[0]->GetXaxis()->FindBin(10.),fhRe1[0]->GetXaxis()->FindBin(20.));
2347  // TH3F * hRe1Pt20 =  (TH3F*)fhRe1[0]->Clone(Form("IMPt20_%s",fCalorimeter.Data()));
2348 //  hRe1Pt20->GetXaxis()->SetRangeUser(10,20);
2349 //  TH1D * hIMPt20 = (TH1D*) hRe1Pt20->Project3D(Form("IMPt20_%s_pz",fCalorimeter.Data()));
2350   hIMPt20->SetLineColor(2);
2351   hIMPt20->SetTitle("10 < p_{T, #gamma#gamma} < 20 GeV/c");
2352   hIMPt20->Draw();
2353
2354   char nameIMF[buffersize];
2355   snprintf(nameIMF,buffersize,"AliAnaPi0_%s_Mgg.eps",fCalorimeter.Data());
2356   cIM->Print(nameIMF);
2357
2358   char namePt[buffersize];
2359   snprintf(namePt,buffersize,"AliAnaPi0_%s_cPt",fCalorimeter.Data());
2360   TCanvas  * cPt = new TCanvas(namePt, "", 400, 10, 600, 700) ;
2361   cPt->Divide(2, 2);
2362
2363   cPt->cd(1) ;
2364   //gPad->SetLogy();
2365   TH1D * hPt = (TH1D*) fhRe1[0]->ProjectionX(Form("Pt0_%s",fCalorimeter.Data()),-1,-1);
2366   hPt->SetLineColor(2);
2367   hPt->SetTitle("No cut on  M_{#gamma#gamma} ");
2368   hPt->Draw();
2369
2370   cPt->cd(2) ;
2371   TH1D * hPtIM1 = (TH1D*)fhRe1[0]->ProjectionX(Form("Pt1_%s",fCalorimeter.Data()), fhRe1[0]->GetZaxis()->FindBin(0.05),fhRe1[0]->GetZaxis()->FindBin(0.21));
2372 //  TH3F * hRe1IM1 = (TH3F*)fhRe1[0]->Clone(Form("Pt1_%s",fCalorimeter.Data()));
2373 //  hRe1IM1->GetZaxis()->SetRangeUser(0.05,0.21);
2374 //  TH1D * hPtIM1 = (TH1D*) hRe1IM1->Project3D("x");
2375   hPtIM1->SetLineColor(2);
2376   hPtIM1->SetTitle("0.05 < M_{#gamma#gamma} < 0.21 GeV/c^{2}");
2377   hPtIM1->Draw();
2378
2379   cPt->cd(3) ;
2380   TH1D * hPtIM2 = (TH1D*)fhRe1[0]->ProjectionX(Form("Pt2_%s",fCalorimeter.Data()), fhRe1[0]->GetZaxis()->FindBin(0.09),fhRe1[0]->GetZaxis()->FindBin(0.17));
2381 //  TH3F * hRe1IM2 = (TH3F*)fhRe1[0]->Clone(Form("Pt2_%s",fCalorimeter.Data()));
2382 //  hRe1IM2->GetZaxis()->SetRangeUser(0.09,0.17);
2383 //  TH1D * hPtIM2 = (TH1D*) hRe1IM2->Project3D("x");
2384   hPtIM2->SetLineColor(2);
2385   hPtIM2->SetTitle("0.09 < M_{#gamma#gamma} < 0.17 GeV/c^{2}");
2386   hPtIM2->Draw();
2387
2388   cPt->cd(4) ;
2389   TH1D * hPtIM3 = (TH1D*)fhRe1[0]->ProjectionX(Form("Pt3_%s",fCalorimeter.Data()), fhRe1[0]->GetZaxis()->FindBin(0.11),fhRe1[0]->GetZaxis()->FindBin(0.15));
2390 //  TH3F * hRe1IM3 = (TH3F*)fhRe1[0]->Clone(Form("Pt3_%s",fCalorimeter.Data()));
2391 //  hRe1IM3->GetZaxis()->SetRangeUser(0.11,0.15);
2392 //  TH1D * hPtIM3 = (TH1D*) hRe1IM1->Project3D("x");
2393   hPtIM3->SetLineColor(2);
2394   hPtIM3->SetTitle("0.11 < M_{#gamma#gamma} < 0.15 GeV/c^{2}");
2395   hPtIM3->Draw();
2396
2397   char namePtF[buffersize];
2398   snprintf(namePtF,buffersize,"AliAnaPi0_%s_Pt.eps",fCalorimeter.Data());
2399   cPt->Print(namePtF);
2400
2401   char line[buffersize] ;
2402   snprintf(line,buffersize,".!tar -zcf %s_%s.tar.gz *.eps", GetName(),fCalorimeter.Data()) ;
2403   gROOT->ProcessLine(line);
2404   snprintf(line, buffersize,".!rm -fR AliAnaPi0_%s*.eps",fCalorimeter.Data());
2405   gROOT->ProcessLine(line);
2406
2407   printf(" AliAnaPi0::Terminate() - !! All the eps files are in %s_%s.tar.gz !!!\n", GetName(), fCalorimeter.Data());
2408
2409 }
2410   //____________________________________________________________________________________________________________________________________________________
2411 Int_t AliAnaPi0::GetEventIndex(AliAODPWG4Particle * part, Double_t * vert)
2412 {
2413   // retieves the event index and checks the vertex
2414   //    in the mixed buffer returns -2 if vertex NOK
2415   //    for normal events   returns 0 if vertex OK and -1 if vertex NOK
2416
2417   Int_t evtIndex = -1 ;
2418   if(GetReader()->GetDataType()!=AliCaloTrackReader::kMC){
2419
2420     if (GetMixedEvent()){
2421
2422       evtIndex = GetMixedEvent()->EventIndexForCaloCluster(part->GetCaloLabel(0)) ;
2423       GetVertex(vert,evtIndex);
2424
2425       if(TMath::Abs(vert[2])> GetZvertexCut())
2426         evtIndex = -2 ; //Event can not be used (vertex, centrality,... cuts not fulfilled)
2427     } else {// Single event
2428
2429       GetVertex(vert);
2430
2431       if(TMath::Abs(vert[2])> GetZvertexCut())
2432         evtIndex = -1 ; //Event can not be used (vertex, centrality,... cuts not fulfilled)
2433       else
2434         evtIndex = 0 ;
2435     }
2436   }//No MC reader
2437   else {
2438     evtIndex = 0;
2439     vert[0] = 0. ;
2440     vert[1] = 0. ;
2441     vert[2] = 0. ;
2442   }
2443
2444   return evtIndex ;
2445 }
2446