PWG4/PartCorrDep/AliAnaPhoton.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 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: AliAnaPhoton.cxx 28688 2008-09-11 15:04:07Z gconesab $ */
  16
  17 //_________________________________________________________________________
  18 //
  19 // Class for the photon identification.
  20 // Clusters from calorimeters are identified as photons
  21 // and kept in the AOD. Few histograms produced.
  22 // Produces input for other analysis classes like AliAnaPi0,
  23 // AliAnaParticleHadronCorrelation ...
  24 //
  25 // -- Author: Gustavo Conesa (LNF-INFN)
  26 //////////////////////////////////////////////////////////////////////////////
  27
  28
  29 // --- ROOT system ---
  30 #include <TH2F.h>
  31 #include <TH3D.h>
  32 #include <TClonesArray.h>
  33 #include <TObjString.h>
  34 //#include <Riostream.h>
  35 #include "TParticle.h"
  36 #include "TDatabasePDG.h"
  37
  38 // --- Analysis system ---
  39 #include "AliAnaPhoton.h"
  40 #include "AliCaloTrackReader.h"
  41 #include "AliStack.h"
  42 #include "AliCaloPID.h"
  43 #include "AliMCAnalysisUtils.h"
  44 #include "AliFiducialCut.h"
  45 #include "AliVCluster.h"
  46 #include "AliAODMCParticle.h"
  47 #include "AliMixedEvent.h"
  48
  49
  50 ClassImp(AliAnaPhoton)
  51
  52 //____________________________________________________________________________
  53   AliAnaPhoton::AliAnaPhoton() :
  54     AliAnaPartCorrBaseClass(), fCalorimeter(""),
  55     fMinDist(0.),fMinDist2(0.),fMinDist3(0.),fRejectTrackMatch(0),
  56     fCheckConversion(kFALSE),fAddConvertedPairsToAOD(kFALSE), fMassCut(0),
  57     fTimeCutMin(-1), fTimeCutMax(9999999), fNCellsCut(0),
  58     fConvDEtaCut(2.),fConvDPhiCut(6.),
  59     //fhVertex(0),
  60     fhNtraNclu(0), fhNCellsPt(0),
  61     fhPtPhoton(0),fhPhiPhoton(0),fhEtaPhoton(0), fhEtaPhiPhoton(0), fhEtaPhi05Photon(0),
  62     fhPtPhotonConv(0), fhEtaPhiPhotonConv(0),fhEtaPhi05PhotonConv(0),
  63     fhConvDeltaEta(0),fhConvDeltaPhi(0), fhConvDeltaEtaPhi(0), fhConvAsym(0), fhConvPt(0),
  64     //MC
  65     fhDeltaE(0), fhDeltaPt(0),fhRatioE(0), fhRatioPt(0),fh2E(0),fh2Pt(0),
  66     fhPtMCPhoton(0),fhPhiMCPhoton(0),fhEtaMCPhoton(0),
  67     fhPtPrompt(0),fhPhiPrompt(0),fhEtaPrompt(0),
  68     fhPtFragmentation(0),fhPhiFragmentation(0),fhEtaFragmentation(0),
  69     fhPtISR(0),fhPhiISR(0),fhEtaISR(0),
  70     fhPtPi0Decay(0),fhPhiPi0Decay(0),fhEtaPi0Decay(0),
  71     fhPtOtherDecay(0),  fhPhiOtherDecay(0),  fhEtaOtherDecay(0),
  72     fhPtConversion(0),  fhPhiConversion(0),  fhEtaConversion(0),fhPtConversionTagged(0),fhEtaPhiConversion(0),fhEtaPhi05Conversion(0),
  73     fhPtAntiNeutron(0), fhPhiAntiNeutron(0), fhEtaAntiNeutron(0),
  74     fhPtAntiProton(0),  fhPhiAntiProton(0),  fhEtaAntiProton(0),
  75     fhPtUnknown(0),     fhPhiUnknown(0),     fhEtaUnknown(0),
  76     fhConvDeltaEtaMCConversion(0),  fhConvDeltaPhiMCConversion(0),  fhConvDeltaEtaPhiMCConversion(0),  fhConvAsymMCConversion(0),  fhConvPtMCConversion(0),  fhConvDispersionMCConversion(0), fhConvM02MCConversion(0),
  77     fhConvDeltaEtaMCAntiNeutron(0), fhConvDeltaPhiMCAntiNeutron(0), fhConvDeltaEtaPhiMCAntiNeutron(0), fhConvAsymMCAntiNeutron(0), fhConvPtMCAntiNeutron(0), fhConvDispersionMCAntiNeutron(0),fhConvM02MCAntiNeutron(0),
  78     fhConvDeltaEtaMCAntiProton(0),  fhConvDeltaPhiMCAntiProton(0),  fhConvDeltaEtaPhiMCAntiProton(0),  fhConvAsymMCAntiProton(0),  fhConvPtMCAntiProton(0),  fhConvDispersionMCAntiProton(0), fhConvM02MCAntiProton(0),
  79     fhConvDeltaEtaMCString(0),      fhConvDeltaPhiMCString(0),      fhConvDeltaEtaPhiMCString(0),      fhConvAsymMCString(0),      fhConvPtMCString(0),      fhConvDispersionMCString(0),     fhConvM02MCString(0)
  80 {
  81   //default ctor
  82
  83   //Initialize parameters
  84   InitParameters();
  85
  86 }//____________________________________________________________________________
  87 AliAnaPhoton::~AliAnaPhoton()
  88 {
  89   //dtor
  90
  91 }
  92
  93 //________________________________________________________________________
  94 TObjString *  AliAnaPhoton::GetAnalysisCuts()
  95 {
  96   //Save parameters used for analysis
  97   TString parList ; //this will be list of parameters used for this analysis.
  98   const Int_t buffersize = 255;
  99   char onePar[buffersize] ;
 100
 101   snprintf(onePar,buffersize,"--- AliAnaPhoton ---\n") ;
 102   parList+=onePar ;
 103   snprintf(onePar,buffersize,"Calorimeter: %s\n",fCalorimeter.Data()) ;
 104   parList+=onePar ;
 105   snprintf(onePar,buffersize,"fMinDist =%2.2f (Minimal distance to bad channel to accept cluster) \n",fMinDist) ;
 106   parList+=onePar ;
 107   snprintf(onePar,buffersize,"fMinDist2=%2.2f (Cuts on Minimal distance to study acceptance evaluation) \n",fMinDist2) ;
 108   parList+=onePar ;
 109   snprintf(onePar,buffersize,"fMinDist3=%2.2f (One more cut on distance used for acceptance-efficiency study) \n",fMinDist3) ;
 110   parList+=onePar ;
 111   snprintf(onePar,buffersize,"fRejectTrackMatch: %d\n",fRejectTrackMatch) ;
 112   parList+=onePar ;
 113   snprintf(onePar,buffersize,"Conversion Selection: fConvDEtaCut %f fConvDPhiCut %f\n",fConvDEtaCut, fConvDPhiCut) ;
 114   parList+=onePar ;
 115
 116   //Get parameters set in base class.
 117   parList += GetBaseParametersList() ;
 118
 119   //Get parameters set in PID class.
 120   parList += GetCaloPID()->GetPIDParametersList() ;
 121
 122   //Get parameters set in FiducialCut class (not available yet)
 123   //parlist += GetFidCut()->GetFidCutParametersList()
 124
 125   return new TObjString(parList) ;
 126 }
 127
 128
 129 //________________________________________________________________________
 130 TList *  AliAnaPhoton::GetCreateOutputObjects()
 131 {
 132   // Create histograms to be saved in output file and
 133   // store them in outputContainer
 134   TList * outputContainer = new TList() ;
 135   outputContainer->SetName("PhotonHistos") ;
 136
 137   Int_t nptbins  = GetHistoPtBins();
 138   Int_t nphibins = GetHistoPhiBins();
 139   Int_t netabins = GetHistoEtaBins();
 140   Float_t ptmax  = GetHistoPtMax();
 141   Float_t phimax = GetHistoPhiMax();
 142   Float_t etamax = GetHistoEtaMax();
 143   Float_t ptmin  = GetHistoPtMin();
 144   Float_t phimin = GetHistoPhiMin();
 145   Float_t etamin = GetHistoEtaMin();
 146
 147   //Histograms of highest Photon identified in Event
 148 //  fhVertex  = new TH3D ("Vertex","vertex position", 20,-10.,10., 20,-10.,10., 80,-40.,40.);
 149 //  fhVertex->SetXTitle("X");
 150 //  fhVertex->SetYTitle("Y");
 151 //  fhVertex->SetZTitle("Z");
 152 //  outputContainer->Add(fhVertex);
 153
 154   fhNtraNclu  = new TH2F ("hNtracksNcluster","# of tracks vs # of clusters", 500,0,500, 500,0,500);
 155   fhNtraNclu->SetXTitle("# of tracks");
 156   fhNtraNclu->SetYTitle("# of clusters");
 157   outputContainer->Add(fhNtraNclu);
 158
 159   fhNCellsPt  = new TH2F ("hNCellsPt","# of tracks vs # of clusters", nptbins,ptmin, ptmax, 10,0,10);
 160   fhNCellsPt->SetXTitle("p_{T}");
 161   fhNCellsPt->SetYTitle("# of cells in cluster");
 162   outputContainer->Add(fhNCellsPt);
 163
 164   fhPtPhoton  = new TH1F("hPtPhoton","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
 165   fhPtPhoton->SetYTitle("N");
 166   fhPtPhoton->SetXTitle("p_{T #gamma}(GeV/c)");
 167   outputContainer->Add(fhPtPhoton) ;
 168
 169   fhPhiPhoton  = new TH2F
 170     ("hPhiPhoton","#phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 171   fhPhiPhoton->SetYTitle("#phi (rad)");
 172   fhPhiPhoton->SetXTitle("p_{T #gamma} (GeV/c)");
 173   outputContainer->Add(fhPhiPhoton) ;
 174
 175   fhEtaPhoton  = new TH2F
 176     ("hEtaPhoton","#eta_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 177   fhEtaPhoton->SetYTitle("#eta");
 178   fhEtaPhoton->SetXTitle("p_{T #gamma} (GeV/c)");
 179   outputContainer->Add(fhEtaPhoton) ;
 180
 181   fhEtaPhiPhoton  = new TH2F
 182   ("hEtaPhiPhoton","#eta vs #phi",netabins,etamin,etamax,nphibins,phimin,phimax);
 183   fhEtaPhiPhoton->SetYTitle("#phi (rad)");
 184   fhEtaPhiPhoton->SetXTitle("#eta");
 185   outputContainer->Add(fhEtaPhiPhoton) ;
 186
 187   fhEtaPhi05Photon  = new TH2F
 188   ("hEtaPhi05Photon","#eta vs #phi",netabins,etamin,etamax,nphibins,phimin,phimax);
 189   fhEtaPhi05Photon->SetYTitle("#phi (rad)");
 190   fhEtaPhi05Photon->SetXTitle("#eta");
 191   outputContainer->Add(fhEtaPhi05Photon) ;
 192
 193
 194   //Conversion
 195
 196   fhPtPhotonConv  = new TH1F("hPtPhotonConv","Number of #gamma over calorimeter, conversion",nptbins,ptmin,ptmax);
 197   fhPtPhotonConv->SetYTitle("N");
 198   fhPtPhotonConv->SetXTitle("p_{T #gamma}(GeV/c)");
 199   outputContainer->Add(fhPtPhotonConv) ;
 200
 201   fhEtaPhiPhotonConv  = new TH2F
 202   ("hEtaPhiPhotonConv","#eta vs #phi",netabins,etamin,etamax,nphibins,phimin,phimax);
 203   fhEtaPhiPhotonConv->SetYTitle("#phi (rad)");
 204   fhEtaPhiPhotonConv->SetXTitle("#eta");
 205   outputContainer->Add(fhEtaPhiPhotonConv) ;
 206
 207   fhEtaPhi05PhotonConv  = new TH2F
 208   ("hEtaPhi05PhotonConv","#eta vs #phi",netabins,etamin,etamax,nphibins,phimin,phimax);
 209   fhEtaPhi05PhotonConv->SetYTitle("#phi (rad)");
 210   fhEtaPhi05PhotonConv->SetXTitle("#eta");
 211   outputContainer->Add(fhEtaPhi05PhotonConv) ;
 212
 213   fhConvDeltaEta  = new TH2F
 214   ("hConvDeltaEta","#Delta #eta of selected conversion pairs",100,0,fMassCut,netabins,-0.5,0.5);
 215   fhConvDeltaEta->SetYTitle("#Delta #eta");
 216   fhConvDeltaEta->SetXTitle("Pair Mass (GeV/c^2)");
 217   outputContainer->Add(fhConvDeltaEta) ;
 218
 219   fhConvDeltaPhi  = new TH2F
 220   ("hConvDeltaPhi","#Delta #phi of selected conversion pairs",100,0,fMassCut,nphibins,-0.5,0.5);
 221   fhConvDeltaPhi->SetYTitle("#Delta #phi");
 222   fhConvDeltaPhi->SetXTitle("Pair Mass (GeV/c^2)");
 223   outputContainer->Add(fhConvDeltaPhi) ;
 224
 225   fhConvDeltaEtaPhi  = new TH2F
 226   ("hConvDeltaEtaPhi","#Delta #eta vs #Delta #phi of selected conversion pairs",netabins,-0.5,0.5,nphibins,-0.5,0.5);
 227   fhConvDeltaEtaPhi->SetYTitle("#Delta #phi");
 228   fhConvDeltaEtaPhi->SetXTitle("#Delta #eta");
 229   outputContainer->Add(fhConvDeltaEtaPhi) ;
 230
 231   fhConvAsym  = new TH2F
 232   ("hConvAsym","Asymmetry of selected conversion pairs",100,0,fMassCut,100,0,1);
 233   fhConvAsym->SetYTitle("Asymmetry");
 234   fhConvAsym->SetXTitle("Pair Mass (GeV/c^2)");
 235   outputContainer->Add(fhConvAsym) ;
 236
 237   fhConvPt  = new TH2F
 238   ("hConvPt","p_{T} of selected conversion pairs",100,0,fMassCut,100,0.,10.);
 239   fhConvPt->SetYTitle("Pair p_{T} (GeV/c)");
 240   fhConvPt->SetXTitle("Pair Mass (GeV/c^2)");
 241   outputContainer->Add(fhConvPt) ;
 242
 243   if(IsDataMC()){
 244     fhDeltaE  = new TH1F ("hDeltaE","MC - Reco E ", 200,-50,50);
 245     fhDeltaE->SetXTitle("#Delta E (GeV)");
 246     outputContainer->Add(fhDeltaE);
 247
 248     fhDeltaPt  = new TH1F ("hDeltaPt","MC - Reco p_{T} ", 200,-50,50);
 249     fhDeltaPt->SetXTitle("#Delta p_{T} (GeV/c)");
 250     outputContainer->Add(fhDeltaPt);
 251
 252     fhRatioE  = new TH1F ("hRatioE","Reco/MC E ", 200,0,2);
 253     fhRatioE->SetXTitle("E_{reco}/E_{gen}");
 254     outputContainer->Add(fhRatioE);
 255
 256     fhRatioPt  = new TH1F ("hRatioPt","Reco/MC p_{T} ", 200,0,2);
 257     fhRatioPt->SetXTitle("p_{T, reco}/p_{T, gen}");
 258     outputContainer->Add(fhRatioPt);
 259
 260     fh2E  = new TH2F ("h2E","E distribution, reconstructed vs generated", nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
 261     fh2E->SetXTitle("E_{rec} (GeV)");
 262     fh2E->SetYTitle("E_{gen} (GeV)");
 263     outputContainer->Add(fh2E);
 264
 265     fh2Pt  = new TH2F ("h2Pt","p_T distribution, reconstructed vs generated", nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
 266     fh2Pt->SetXTitle("p_{T,rec} (GeV/c)");
 267     fh2Pt->SetYTitle("p_{T,gen} (GeV/c)");
 268     outputContainer->Add(fh2Pt);
 269
 270     fhPtMCPhoton  = new TH1F("hPtMCPhoton","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
 271     fhPtMCPhoton->SetYTitle("N");
 272     fhPtMCPhoton->SetXTitle("p_{T #gamma}(GeV/c)");
 273     outputContainer->Add(fhPtMCPhoton) ;
 274
 275     fhPhiMCPhoton  = new TH2F
 276       ("hPhiMCPhoton","#phi_{#gamma}, #gamma in MC",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 277     fhPhiMCPhoton->SetYTitle("#phi");
 278     fhPhiMCPhoton->SetXTitle("p_{T #gamma} (GeV/c)");
 279     outputContainer->Add(fhPhiMCPhoton) ;
 280
 281     fhEtaMCPhoton  = new TH2F
 282       ("hEtaMCPhoton","#eta_{#gamma}, #gamma in MC",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 283     fhEtaMCPhoton->SetYTitle("#eta");
 284     fhEtaMCPhoton->SetXTitle("p_{T #gamma} (GeV/c)");
 285     outputContainer->Add(fhEtaMCPhoton) ;
 286
 287     fhPtPrompt  = new TH1F("hPtMCPrompt","Number of prompt #gamma over calorimeter",nptbins,ptmin,ptmax);
 288     fhPtPrompt->SetYTitle("N");
 289     fhPtPrompt->SetXTitle("p_{T #gamma}(GeV/c)");
 290     outputContainer->Add(fhPtPrompt) ;
 291
 292     fhPhiPrompt  = new TH2F
 293       ("hPhiMCPrompt","#phi_{#gamma}, prompt #gamma in MC",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 294     fhPhiPrompt->SetYTitle("#phi");
 295     fhPhiPrompt->SetXTitle("p_{T #gamma} (GeV/c)");
 296     outputContainer->Add(fhPhiPrompt) ;
 297
 298     fhEtaPrompt  = new TH2F
 299       ("hEtaMCPrompt","#eta_{#gamma}, prompt #gamma in MC",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 300     fhEtaPrompt->SetYTitle("#eta");
 301     fhEtaPrompt->SetXTitle("p_{T #gamma} (GeV/c)");
 302     outputContainer->Add(fhEtaPrompt) ;
 303
 304     fhPtFragmentation  = new TH1F("hPtMCFragmentation","Number of fragmentation #gamma over calorimeter",nptbins,ptmin,ptmax);
 305     fhPtFragmentation->SetYTitle("N");
 306     fhPtFragmentation->SetXTitle("p_{T #gamma}(GeV/c)");
 307     outputContainer->Add(fhPtFragmentation) ;
 308
 309     fhPhiFragmentation  = new TH2F
 310       ("hPhiMCFragmentation","#phi_{#gamma}, fragmentation #gamma in MC",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 311     fhPhiFragmentation->SetYTitle("#phi");
 312     fhPhiFragmentation->SetXTitle("p_{T #gamma} (GeV/c)");
 313     outputContainer->Add(fhPhiFragmentation) ;
 314
 315     fhEtaFragmentation  = new TH2F
 316       ("hEtaMCFragmentation","#eta_{#gamma}, fragmentation #gamma in MC",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 317     fhEtaFragmentation->SetYTitle("#eta");
 318     fhEtaFragmentation->SetXTitle("p_{T #gamma} (GeV/c)");
 319     outputContainer->Add(fhEtaFragmentation) ;
 320
 321     fhPtISR  = new TH1F("hPtMCISR","Number of initial state radiation #gamma over calorimeter",nptbins,ptmin,ptmax);
 322     fhPtISR->SetYTitle("N");
 323     fhPtISR->SetXTitle("p_{T #gamma}(GeV/c)");
 324     outputContainer->Add(fhPtISR) ;
 325
 326     fhPhiISR  = new TH2F
 327       ("hPhiMCISR","#phi_{#gamma} initial state radiation",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 328     fhPhiISR->SetYTitle("#phi");
 329     fhPhiISR->SetXTitle("p_{T #gamma} (GeV/c)");
 330     outputContainer->Add(fhPhiISR) ;
 331
 332     fhEtaISR  = new TH2F
 333       ("hEtaMCISR","#eta_{#gamma} initial state radiation",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 334     fhEtaISR->SetYTitle("#eta");
 335     fhEtaISR->SetXTitle("p_{T #gamma} (GeV/c)");
 336     outputContainer->Add(fhEtaISR) ;
 337
 338     fhPtPi0Decay  = new TH1F("hPtMCPi0Decay","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
 339     fhPtPi0Decay->SetYTitle("N");
 340     fhPtPi0Decay->SetXTitle("p_{T #gamma}(GeV/c)");
 341     outputContainer->Add(fhPtPi0Decay) ;
 342
 343     fhPhiPi0Decay  = new TH2F
 344       ("hPhiMCPi0Decay","#phi_{#gamma}, #pi^{0} decay #gamma in MC",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 345     fhPhiPi0Decay->SetYTitle("#phi");
 346     fhPhiPi0Decay->SetXTitle("p_{T #gamma} (GeV/c)");
 347     outputContainer->Add(fhPhiPi0Decay) ;
 348
 349     fhEtaPi0Decay  = new TH2F
 350       ("hEtaMCPi0Decay","#eta_{#gamma}, #pi^{0} #gamma in MC",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 351     fhEtaPi0Decay->SetYTitle("#eta");
 352     fhEtaPi0Decay->SetXTitle("p_{T #gamma} (GeV/c)");
 353     outputContainer->Add(fhEtaPi0Decay) ;
 354
 355     fhPtOtherDecay  = new TH1F("hPtMCOtherDecay","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
 356     fhPtOtherDecay->SetYTitle("N");
 357     fhPtOtherDecay->SetXTitle("p_{T #gamma}(GeV/c)");
 358     outputContainer->Add(fhPtOtherDecay) ;
 359
 360     fhPhiOtherDecay  = new TH2F
 361       ("hPhiMCOtherDecay","#phi_{#gamma}, other decay #gamma in MC",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 362     fhPhiOtherDecay->SetYTitle("#phi");
 363     fhPhiOtherDecay->SetXTitle("p_{T #gamma} (GeV/c)");
 364     outputContainer->Add(fhPhiOtherDecay) ;
 365
 366     fhEtaOtherDecay  = new TH2F
 367       ("hEtaMCOtherDecay","#eta_{#gamma}, other decay #gamma in MC",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 368     fhEtaOtherDecay->SetYTitle("#eta");
 369     fhEtaOtherDecay->SetXTitle("p_{T #gamma} (GeV/c)");
 370     outputContainer->Add(fhEtaOtherDecay) ;
 371
 372     fhPtConversion  = new TH1F("hPtMCConversion","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
 373     fhPtConversion->SetYTitle("N");
 374     fhPtConversion->SetXTitle("p_{T #gamma}(GeV/c)");
 375     outputContainer->Add(fhPtConversion) ;
 376
 377     fhPhiConversion  = new TH2F
 378       ("hPhiMCConversion","#phi_{#gamma}, conversion #gamma in MC",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 379     fhPhiConversion->SetYTitle("#phi");
 380     fhPhiConversion->SetXTitle("p_{T #gamma} (GeV/c)");
 381     outputContainer->Add(fhPhiConversion) ;
 382
 383     fhEtaConversion  = new TH2F
 384       ("hEtaMCConversion","#eta_{#gamma}, conversion #gamma in MC",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 385     fhEtaConversion->SetYTitle("#eta");
 386     fhEtaConversion->SetXTitle("p_{T #gamma} (GeV/c)");
 387     outputContainer->Add(fhEtaConversion) ;
 388
 389     fhPtConversionTagged  = new TH1F("hPtMCConversionTagged","Number of converted #gamma over calorimeter, tagged as converted",nptbins,ptmin,ptmax);
 390     fhPtConversionTagged->SetYTitle("N");
 391     fhPtConversionTagged->SetXTitle("p_{T #gamma}(GeV/c)");
 392     outputContainer->Add(fhPtConversionTagged) ;
 393
 394     fhEtaPhiConversion  = new TH2F
 395     ("hEtaPhiConversion","#eta vs #phi",netabins,etamin,etamax,nphibins,phimin,phimax);
 396     fhEtaPhiConversion->SetYTitle("#phi (rad)");
 397     fhEtaPhiConversion->SetXTitle("#eta");
 398     outputContainer->Add(fhEtaPhiConversion) ;
 399
 400     fhEtaPhi05Conversion  = new TH2F
 401     ("hEtaPhi05Conversion","#eta vs #phi",netabins,etamin,etamax,nphibins,phimin,phimax);
 402     fhEtaPhi05Conversion->SetYTitle("#phi (rad)");
 403     fhEtaPhi05Conversion->SetXTitle("#eta");
 404     outputContainer->Add(fhEtaPhi05Conversion) ;
 405
 406     fhPtAntiNeutron  = new TH1F("hPtMCAntiNeutron","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
 407     fhPtAntiNeutron->SetYTitle("N");
 408     fhPtAntiNeutron->SetXTitle("p_{T #gamma}(GeV/c)");
 409     outputContainer->Add(fhPtAntiNeutron) ;
 410
 411     fhPhiAntiNeutron  = new TH2F
 412     ("hPhiMCAntiNeutron","#phi_{#gamma}, unknown origin",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 413     fhPhiAntiNeutron->SetYTitle("#phi");
 414     fhPhiAntiNeutron->SetXTitle("p_{T #gamma} (GeV/c)");
 415     outputContainer->Add(fhPhiAntiNeutron) ;
 416
 417     fhEtaAntiNeutron  = new TH2F
 418     ("hEtaMCAntiNeutron","#eta_{#gamma}, unknown origin",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 419     fhEtaAntiNeutron->SetYTitle("#eta");
 420     fhEtaAntiNeutron->SetXTitle("p_{T #gamma} (GeV/c)");
 421     outputContainer->Add(fhEtaAntiNeutron) ;
 422
 423     fhPtAntiProton  = new TH1F("hPtMCAntiProtonConv","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
 424     fhPtAntiProton->SetYTitle("N");
 425     fhPtAntiProton->SetXTitle("p_{T #gamma}(GeV/c)");
 426     outputContainer->Add(fhPtAntiProton) ;
 427
 428     fhPhiAntiProton  = new TH2F
 429     ("hPhiMCAntiProton","#phi_{#gamma}, unknown origin",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 430     fhPhiAntiProton->SetYTitle("#phi");
 431     fhPhiAntiProton->SetXTitle("p_{T #gamma} (GeV/c)");
 432     outputContainer->Add(fhPhiAntiProton) ;
 433
 434     fhEtaAntiProton  = new TH2F
 435     ("hEtaMCAntiProton","#eta_{#gamma}, unknown origin",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 436     fhEtaAntiProton->SetYTitle("#eta");
 437     fhEtaAntiProton->SetXTitle("p_{T #gamma} (GeV/c)");
 438     outputContainer->Add(fhEtaAntiProton) ;
 439
 440     fhPtUnknown  = new TH1F("hPtMCUnknown","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
 441     fhPtUnknown->SetYTitle("N");
 442     fhPtUnknown->SetXTitle("p_{T #gamma}(GeV/c)");
 443     outputContainer->Add(fhPtUnknown) ;
 444
 445     fhPhiUnknown  = new TH2F
 446       ("hPhiMCUnknown","#phi_{#gamma}, unknown origin",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 447     fhPhiUnknown->SetYTitle("#phi");
 448     fhPhiUnknown->SetXTitle("p_{T #gamma} (GeV/c)");
 449     outputContainer->Add(fhPhiUnknown) ;
 450
 451     fhEtaUnknown  = new TH2F
 452       ("hEtaMCUnknown","#eta_{#gamma}, unknown origin",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 453     fhEtaUnknown->SetYTitle("#eta");
 454     fhEtaUnknown->SetXTitle("p_{T #gamma} (GeV/c)");
 455     outputContainer->Add(fhEtaUnknown) ;
 456
 457     fhConvDeltaEtaMCConversion  = new TH2F
 458     ("hConvDeltaEtaMCConversion","#Delta #eta of selected conversion pairs from real conversions",100,0,fMassCut,netabins,-0.5,0.5);
 459     fhConvDeltaEtaMCConversion->SetYTitle("#Delta #eta");
 460     fhConvDeltaEtaMCConversion->SetXTitle("Pair Mass (GeV/c^2)");
 461     outputContainer->Add(fhConvDeltaEtaMCConversion) ;
 462
 463     fhConvDeltaPhiMCConversion  = new TH2F
 464     ("hConvDeltaPhiMCConversion","#Delta #phi of selected conversion pairs from real conversions",100,0,fMassCut,nphibins,-0.5,0.5);
 465     fhConvDeltaPhiMCConversion->SetYTitle("#Delta #phi");
 466     fhConvDeltaPhiMCConversion->SetXTitle("Pair Mass (GeV/c^2)");
 467     outputContainer->Add(fhConvDeltaPhiMCConversion) ;
 468
 469     fhConvDeltaEtaPhiMCConversion  = new TH2F
 470     ("hConvDeltaEtaPhiMCConversion","#Delta #eta vs #Delta #phi of selected conversion pairs, from real conversions",netabins,-0.5,0.5,nphibins,-0.5,0.5);
 471     fhConvDeltaEtaPhiMCConversion->SetYTitle("#Delta #phi");
 472     fhConvDeltaEtaPhiMCConversion->SetXTitle("#Delta #eta");
 473     outputContainer->Add(fhConvDeltaEtaPhiMCConversion) ;
 474
 475     fhConvAsymMCConversion  = new TH2F
 476     ("hConvAsymMCConversion","Asymmetry of selected conversion pairs from real conversions",100,0,fMassCut,100,0,1);
 477     fhConvAsymMCConversion->SetYTitle("Asymmetry");
 478     fhConvAsymMCConversion->SetXTitle("Pair Mass (GeV/c^2)");
 479     outputContainer->Add(fhConvAsymMCConversion) ;
 480
 481     fhConvPtMCConversion  = new TH2F
 482     ("hConvPtMCConversion","p_{T} of selected conversion pairs from real conversions",100,0,fMassCut,100,0.,10.);
 483     fhConvPtMCConversion->SetYTitle("Pair p_{T} (GeV/c)");
 484     fhConvPtMCConversion->SetXTitle("Pair Mass (GeV/c^2)");
 485     outputContainer->Add(fhConvPtMCConversion) ;
 486
 487     fhConvDispersionMCConversion  = new TH2F
 488     ("hConvDispersionMCConversion","p_{T} of selected conversion pairs from real conversions",100,0.,1.,100,0.,1.);
 489     fhConvDispersionMCConversion->SetYTitle("Dispersion cluster 1");
 490     fhConvDispersionMCConversion->SetXTitle("Dispersion cluster 2");
 491     outputContainer->Add(fhConvDispersionMCConversion) ;
 492
 493     fhConvM02MCConversion  = new TH2F
 494     ("hConvM02MCConversion","p_{T} of selected conversion pairs from string",100,0.,1.,100,0.,1.);
 495     fhConvM02MCConversion->SetYTitle("M02 cluster 1");
 496     fhConvM02MCConversion->SetXTitle("M02 cluster 2");
 497     outputContainer->Add(fhConvM02MCConversion) ;
 498
 499     fhConvDeltaEtaMCAntiNeutron  = new TH2F
 500     ("hConvDeltaEtaMCAntiNeutron","#Delta #eta of selected conversion pairs from anti-neutrons",100,0,fMassCut,netabins,-0.5,0.5);
 501     fhConvDeltaEtaMCAntiNeutron->SetYTitle("#Delta #eta");
 502     fhConvDeltaEtaMCAntiNeutron->SetXTitle("Pair Mass (GeV/c^2)");
 503     outputContainer->Add(fhConvDeltaEtaMCAntiNeutron) ;
 504
 505     fhConvDeltaPhiMCAntiNeutron  = new TH2F
 506     ("hConvDeltaPhiMCAntiNeutron","#Delta #phi of selected conversion pairs from anti-neutrons",100,0,fMassCut,nphibins,-0.5,0.5);
 507     fhConvDeltaPhiMCAntiNeutron->SetYTitle("#Delta #phi");
 508     fhConvDeltaPhiMCAntiNeutron->SetXTitle("Pair Mass (GeV/c^2)");
 509     outputContainer->Add(fhConvDeltaPhiMCAntiNeutron) ;
 510
 511     fhConvDeltaEtaPhiMCAntiNeutron  = new TH2F
 512     ("hConvDeltaEtaPhiMCAntiNeutron","#Delta #eta vs #Delta #phi of selected conversion pairs from anti-neutrons",netabins,-0.5,0.5,nphibins,-0.5,0.5);
 513     fhConvDeltaEtaPhiMCAntiNeutron->SetYTitle("#Delta #phi");
 514     fhConvDeltaEtaPhiMCAntiNeutron->SetXTitle("#Delta #eta");
 515     outputContainer->Add(fhConvDeltaEtaPhiMCAntiNeutron) ;
 516
 517     fhConvAsymMCAntiNeutron  = new TH2F
 518     ("hConvAsymMCAntiNeutron","Asymmetry of selected conversion pairs from anti-neutrons",100,0,fMassCut,100,0,1);
 519     fhConvAsymMCAntiNeutron->SetYTitle("Asymmetry");
 520     fhConvAsymMCAntiNeutron->SetXTitle("Pair Mass (GeV/c^2)");
 521     outputContainer->Add(fhConvAsymMCAntiNeutron) ;
 522
 523     fhConvPtMCAntiNeutron  = new TH2F
 524     ("hConvPtMCAntiNeutron","p_{T} of selected conversion pairs from anti-neutrons",100,0,fMassCut,100,0.,10.);
 525     fhConvPtMCAntiNeutron->SetYTitle("Pair p_{T} (GeV/c)");
 526     fhConvPtMCAntiNeutron->SetXTitle("Pair Mass (GeV/c^2)");
 527     outputContainer->Add(fhConvPtMCAntiNeutron) ;
 528
 529     fhConvDispersionMCAntiNeutron  = new TH2F
 530     ("hConvDispersionMCAntiNeutron","p_{T} of selected conversion pairs from anti-neutrons",100,0.,1.,100,0.,1.);
 531     fhConvDispersionMCAntiNeutron->SetYTitle("Dispersion cluster 1");
 532     fhConvDispersionMCAntiNeutron->SetXTitle("Dispersion cluster 2");
 533     outputContainer->Add(fhConvDispersionMCAntiNeutron) ;
 534
 535     fhConvM02MCAntiNeutron  = new TH2F
 536     ("hConvM02MCAntiNeutron","p_{T} of selected conversion pairs from string",100,0.,1.,100,0.,1.);
 537     fhConvM02MCAntiNeutron->SetYTitle("M02 cluster 1");
 538     fhConvM02MCAntiNeutron->SetXTitle("M02 cluster 2");
 539     outputContainer->Add(fhConvM02MCAntiNeutron) ;
 540
 541     fhConvDeltaEtaMCAntiProton  = new TH2F
 542     ("hConvDeltaEtaMCAntiProton","#Delta #eta of selected conversion pairs from anti-protons",100,0,fMassCut,netabins,-0.5,0.5);
 543     fhConvDeltaEtaMCAntiProton->SetYTitle("#Delta #eta");
 544     fhConvDeltaEtaMCAntiProton->SetXTitle("Pair Mass (GeV/c^2)");
 545     outputContainer->Add(fhConvDeltaEtaMCAntiProton) ;
 546
 547     fhConvDeltaPhiMCAntiProton  = new TH2F
 548     ("hConvDeltaPhiMCAntiProton","#Delta #phi of selected conversion pairs from anti-protons",100,0,fMassCut,nphibins,-0.5,0.5);
 549     fhConvDeltaPhiMCAntiProton->SetYTitle("#Delta #phi");
 550     fhConvDeltaPhiMCAntiProton->SetXTitle("Pair Mass (GeV/c^2)");
 551     outputContainer->Add(fhConvDeltaPhiMCAntiProton) ;
 552
 553     fhConvDeltaEtaPhiMCAntiProton  = new TH2F
 554     ("hConvDeltaEtaPhiMCAntiProton","#Delta #eta vs #Delta #phi of selected conversion pairs from anti-protons",netabins,-0.5,0.5,nphibins,-0.5,0.5);
 555     fhConvDeltaEtaPhiMCAntiProton->SetYTitle("#Delta #phi");
 556     fhConvDeltaEtaPhiMCAntiProton->SetXTitle("#Delta #eta");
 557     outputContainer->Add(fhConvDeltaEtaPhiMCAntiProton) ;
 558
 559     fhConvAsymMCAntiProton  = new TH2F
 560     ("hConvAsymMCAntiProton","Asymmetry of selected conversion pairs from anti-protons",100,0,fMassCut,100,0,1);
 561     fhConvAsymMCAntiProton->SetYTitle("Asymmetry");
 562     fhConvAsymMCAntiProton->SetXTitle("Pair Mass (GeV/c^2)");
 563     outputContainer->Add(fhConvAsymMCAntiProton) ;
 564
 565     fhConvPtMCAntiProton  = new TH2F
 566     ("hConvPtMCAntiProton","p_{T} of selected conversion pairs from anti-protons",100,0,fMassCut,100,0.,10.);
 567     fhConvPtMCAntiProton->SetYTitle("Pair p_{T} (GeV/c)");
 568     fhConvPtMCAntiProton->SetXTitle("Pair Mass (GeV/c^2)");
 569     outputContainer->Add(fhConvPtMCAntiProton) ;
 570
 571     fhConvDispersionMCAntiProton  = new TH2F
 572     ("hConvDispersionMCAntiProton","p_{T} of selected conversion pairs from anti-protons",100,0.,1.,100,0.,1.);
 573     fhConvDispersionMCAntiProton->SetYTitle("Dispersion cluster 1");
 574     fhConvDispersionMCAntiProton->SetXTitle("Dispersion cluster 2");
 575     outputContainer->Add(fhConvDispersionMCAntiProton) ;
 576
 577     fhConvM02MCAntiProton  = new TH2F
 578     ("hConvM02MCAntiProton","p_{T} of selected conversion pairs from string",100,0.,1.,100,0.,1.);
 579     fhConvM02MCAntiProton->SetYTitle("M02 cluster 1");
 580     fhConvM02MCAntiProton->SetXTitle("M02 cluster 2");
 581     outputContainer->Add(fhConvM02MCAntiProton) ;
 582
 583     fhConvDeltaEtaMCString  = new TH2F
 584     ("hConvDeltaEtaMCString","#Delta #eta of selected conversion pairs from string",100,0,fMassCut,netabins,-0.5,0.5);
 585     fhConvDeltaEtaMCString->SetYTitle("#Delta #eta");
 586     fhConvDeltaEtaMCString->SetXTitle("Pair Mass (GeV/c^2)");
 587     outputContainer->Add(fhConvDeltaEtaMCString) ;
 588
 589     fhConvDeltaPhiMCString  = new TH2F
 590     ("hConvDeltaPhiMCString","#Delta #phi of selected conversion pairs from string",100,0,fMassCut,nphibins,-0.5,0.5);
 591     fhConvDeltaPhiMCString->SetYTitle("#Delta #phi");
 592     fhConvDeltaPhiMCString->SetXTitle("Pair Mass (GeV/c^2)");
 593     outputContainer->Add(fhConvDeltaPhiMCString) ;
 594
 595     fhConvDeltaEtaPhiMCString  = new TH2F
 596     ("hConvDeltaEtaPhiMCString","#Delta #eta vs #Delta #phi of selected conversion pairs from string",netabins,-0.5,0.5,nphibins,-0.5,0.5);
 597     fhConvDeltaEtaPhiMCString->SetYTitle("#Delta #phi");
 598     fhConvDeltaEtaPhiMCString->SetXTitle("#Delta #eta");
 599     outputContainer->Add(fhConvDeltaEtaPhiMCString) ;
 600
 601     fhConvAsymMCString  = new TH2F
 602     ("hConvAsymMCString","Asymmetry of selected conversion pairs from string",100,0,fMassCut,100,0,1);
 603     fhConvAsymMCString->SetYTitle("Asymmetry");
 604     fhConvAsymMCString->SetXTitle("Pair Mass (GeV/c^2)");
 605     outputContainer->Add(fhConvAsymMCString) ;
 606
 607     fhConvPtMCString  = new TH2F
 608     ("hConvPtMCString","p_{T} of selected conversion pairs from string",100,0,fMassCut,100,0.,10.);
 609     fhConvPtMCString->SetYTitle("Pair p_{T} (GeV/c)");
 610     fhConvPtMCString->SetXTitle("Pair Mass (GeV/c^2)");
 611     outputContainer->Add(fhConvPtMCString) ;
 612
 613     fhConvDispersionMCString  = new TH2F
 614     ("hConvDispersionMCString","p_{T} of selected conversion pairs from string",100,0.,1.,100,0.,1.);
 615     fhConvDispersionMCString->SetYTitle("Dispersion cluster 1");
 616     fhConvDispersionMCString->SetXTitle("Dispersion cluster 2");
 617     outputContainer->Add(fhConvDispersionMCString) ;
 618
 619     fhConvM02MCString  = new TH2F
 620     ("hConvM02MCString","p_{T} of selected conversion pairs from string",100,0.,1.,100,0.,1.);
 621     fhConvM02MCString->SetYTitle("M02 cluster 1");
 622     fhConvM02MCString->SetXTitle("M02 cluster 2");
 623     outputContainer->Add(fhConvM02MCString) ;
 624
 625
 626   }//Histos with MC
 627
 628   return outputContainer ;
 629
 630 }
 631
 632 //____________________________________________________________________________
 633 void AliAnaPhoton::Init()
 634 {
 635
 636   //Init
 637   //Do some checks
 638   if(fCalorimeter == "PHOS" && !GetReader()->IsPHOSSwitchedOn() && NewOutputAOD()){
 639     printf("AliAnaPhoton::Init() - !!STOP: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!\n");
 640     abort();
 641   }
 642   else  if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn() && NewOutputAOD()){
 643     printf("AliAnaPhoton::Init() - !!STOP: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!\n");
 644     abort();
 645   }
 646
 647 }
 648
 649
 650 //____________________________________________________________________________
 651 void AliAnaPhoton::InitParameters()
 652 {
 653
 654   //Initialize the parameters of the analysis.
 655   AddToHistogramsName("AnaPhoton_");
 656
 657   fCalorimeter = "EMCAL" ;
 658   fMinDist     = 2.;
 659   fMinDist2    = 4.;
 660   fMinDist3    = 5.;
 661   fMassCut     = 0.03; //30 MeV
 662
 663   fTimeCutMin  = -1;
 664   fTimeCutMax  = 9999999;
 665   fNCellsCut   = 0;
 666
 667   fRejectTrackMatch       = kTRUE ;
 668   fCheckConversion        = kFALSE;
 669   fAddConvertedPairsToAOD = kFALSE;
 670
 671 }
 672
 673 //__________________________________________________________________
 674 void  AliAnaPhoton::MakeAnalysisFillAOD()
 675 {
 676   //Do photon analysis and fill aods
 677
 678   //Get the vertex
 679   Double_t v[3] = {0,0,0}; //vertex ;
 680   GetReader()->GetVertex(v);
 681
 682   //Select the Calorimeter of the photon
 683   TObjArray * pl = 0x0;
 684   if(fCalorimeter == "PHOS")
 685     pl = GetPHOSClusters();
 686   else if (fCalorimeter == "EMCAL")
 687     pl = GetEMCALClusters();
 688
 689   if(!pl) {
 690     Info("MakeAnalysisFillAOD","TObjArray with %s clusters is NULL!\n",fCalorimeter.Data());
 691     return;
 692   }
 693
 694   //Init arrays, variables, get number of clusters
 695   TLorentzVector mom, mom2 ;
 696   Int_t nCaloClusters = pl->GetEntriesFast();
 697   //List to be used in conversion analysis, to tag the cluster as candidate for conversion
 698   Bool_t * indexConverted = new Bool_t[nCaloClusters];
 699   for (Int_t i = 0; i < nCaloClusters; i++)
 700     indexConverted[i] = kFALSE;
 701
 702   if(GetDebug() > 0) printf("AliAnaPhoton::MakeAnalysisFillAOD() - input %s cluster entries %d\n", fCalorimeter.Data(), nCaloClusters);
 703
 704   //----------------------------------------------------
 705   // Fill AOD with PHOS/EMCAL AliAODPWG4Particle objects
 706   //----------------------------------------------------
 707   // Loop on clusters
 708   for(Int_t icalo = 0; icalo < nCaloClusters; icalo++){
 709
 710           AliVCluster * calo =  (AliVCluster*) (pl->At(icalo));
 711     //printf("calo %d, %f\n",icalo,calo->E());
 712
 713     //Get the index where the cluster comes, to retrieve the corresponding vertex
 714     Int_t evtIndex = 0 ;
 715     if (GetMixedEvent()) {
 716       evtIndex=GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ;
 717       //Get the vertex and check it is not too large in z
 718       if(TMath::Abs(GetVertex(evtIndex)[2])> GetZvertexCut()) continue;
 719     }
 720
 721     //Cluster selection, not charged, with photon id and in fiducial cut
 722
 723     //Input from second AOD?
 724     //Int_t input = 0;
 725     //    if (fCalorimeter == "EMCAL" && GetReader()->GetEMCALClustersNormalInputEntries() <= icalo)
 726     //      input = 1 ;
 727     //    else if(fCalorimeter == "PHOS"  && GetReader()->GetPHOSClustersNormalInputEntries()  <= icalo)
 728     //      input = 1;
 729
 730     //Get Momentum vector,
 731     //if (input == 0)
 732     if(GetReader()->GetDataType() != AliCaloTrackReader::kMC){
 733       calo->GetMomentum(mom,GetVertex(evtIndex)) ;}//Assume that come from vertex in straight line
 734     else{
 735       Double_t vertex[]={0,0,0};
 736       calo->GetMomentum(mom,vertex) ;
 737     }
 738
 739     //    else if(input == 1)
 740     //      calo->GetMomentum(mom,vertex2);//Assume that come from vertex in straight line
 741
 742     //--------------------------------------
 743     // Cluster selection
 744     //--------------------------------------
 745     if(GetDebug() > 2)
 746       printf("AliAnaPhoton::MakeAnalysisFillAOD() Current Event %d; Before selection : E %2.2f, pT %2.2f, Ecl %2.2f, phi %2.2f, eta %2.2f\n",
 747              GetReader()->GetEventNumber(),
 748              mom.E(), mom.Pt(),calo->E(),mom.Phi()*TMath::RadToDeg(),mom.Eta());
 749
 750     //.......................................
 751     //If too small or big pt, skip it
 752     if(mom.Pt() < GetMinPt() || mom.Pt() > GetMaxPt() ) continue ;
 753     if(GetDebug() > 2) printf("\t Cluster %d Pass Pt Cut \n",icalo);
 754
 755     //.......................................
 756     // TOF cut, BE CAREFUL WITH THIS CUT
 757     Double_t tof = calo->GetTOF()*1e9;
 758     if(tof < fTimeCutMin || tof > fTimeCutMax) continue;
 759           if(GetDebug() > 2)  printf("\t Cluster %d Pass Time Cut \n",icalo);
 760
 761     //.......................................
 762     if(calo->GetNCells() <= fNCellsCut && GetReader()->GetDataType() != AliCaloTrackReader::kMC) continue;
 763     if(GetDebug() > 2) printf("\t Cluster %d Pass NCell Cut \n",icalo);
 764
 765     //.......................................
 766     //Check acceptance selection
 767     if(IsFiducialCutOn()){
 768       Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter) ;
 769       if(! in ) continue ;
 770     }
 771     if(GetDebug() > 2) printf("Fiducial cut passed \n");
 772
 773     //.......................................
 774     //Skip matched clusters with tracks
 775     if(fRejectTrackMatch){
 776       if(IsTrackMatched(calo)) {
 777         if(GetDebug() > 2) printf("\t Reject matched clusters\n");
 778         continue ;
 779       }
 780       else
 781         if(GetDebug() > 2)  printf(" matching cut passed cut passed \n");
 782     }// reject matched clusters
 783
 784     //.......................................
 785     //Check Distance to Bad channel, set bit.
 786     Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
 787     if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;
 788     if(distBad < fMinDist) {//In bad channel (PHOS cristal size 2.2x2.2 cm), EMCAL ( cell units )
 789       continue ;
 790     }
 791     else if(GetDebug() > 2) printf("\t Bad channel cut passed %4.2f > %2.2f \n",distBad, fMinDist);
 792
 793     if(GetDebug() > 0)
 794       printf("AliAnaPhoton::MakeAnalysisFillAOD() Current Event %d; After  selection : E %2.2f, pT %2.2f, Ecl %2.2f, phi %2.2f, eta %2.2f\n",
 795              GetReader()->GetEventNumber(),
 796              mom.E(), mom.Pt(),calo->E(),mom.Phi()*TMath::RadToDeg(),mom.Eta());
 797
 798
 799     //----------------------------
 800     //Create AOD for analysis
 801     //----------------------------
 802     AliAODPWG4Particle aodph = AliAODPWG4Particle(mom);
 803
 804     //...............................................
 805     //Set the indeces of the original caloclusters (MC, ID), and calorimeter
 806     Int_t label = calo->GetLabel();
 807     aodph.SetLabel(label);
 808     //aodph.SetInputFileIndex(input);
 809     aodph.SetCaloLabel(calo->GetID(),-1);
 810     aodph.SetDetector(fCalorimeter);
 811     //printf("Index %d, Id %d\n",icalo, calo->GetID());
 812
 813     //printf("Cluster %d Pass Bad Dist Cut \n",icalo);
 814
 815     //...............................................
 816     //Set bad channel distance bit
 817     if     (distBad > fMinDist3) aodph.SetDistToBad(2) ;
 818     else if(distBad > fMinDist2) aodph.SetDistToBad(1) ;
 819     else                         aodph.SetDistToBad(0) ;
 820     //printf("DistBad %f Bit %d\n",distBad, aodph.DistToBad());
 821
 822     //...............................................
 823     //Set number of cells in this cluster
 824     //Temporary patch FIXME
 825     aodph.SetBtag(calo->GetNCells());
 826     // MEFIX
 827
 828     //-------------------------------------
 829     //PID selection or bit setting
 830     //-------------------------------------
 831     // MC
 832     if(GetReader()->GetDataType() == AliCaloTrackReader::kMC){
 833       //Get most probable PID, check PID weights (in MC this option is mandatory)
 834       aodph.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,calo->GetPID(),mom.E()));//PID with weights
 835       if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PDG of identified particle %d\n",aodph.GetPdg());
 836       //If primary is not photon, skip it.
 837       if(aodph.GetPdg() != AliCaloPID::kPhoton) continue ;
 838     }
 839     //...............................................
 840     // Data, PID check on
 841     else if(IsCaloPIDOn()){
 842       //Get most probable PID, 2 options check PID weights
 843       //or redo PID, recommended option for EMCal.
 844       if(!IsCaloPIDRecalculationOn())
 845         aodph.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,calo->GetPID(),mom.E()));//PID with weights
 846       else
 847         aodph.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,mom,calo));//PID recalculated
 848
 849       if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PDG of identified particle %d\n",aodph.GetPdg());
 850
 851       //If cluster does not pass pid, not photon, skip it.
 852       if(aodph.GetPdg() != AliCaloPID::kPhoton) continue ;
 853
 854     }
 855     //...............................................
 856     // Data, PID check off
 857     else{
 858       //Set PID bits for later selection (AliAnaPi0 for example)
 859       //GetPDG already called in SetPIDBits.
 860       GetCaloPID()->SetPIDBits(fCalorimeter,calo,&aodph, GetCaloUtils());
 861       if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PID Bits set \n");
 862     }
 863
 864     if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - Photon selection cuts passed: pT %3.2f, pdg %d\n",aodph.Pt(), aodph.GetPdg());
 865
 866     //--------------------------------------------------------------------------------------
 867     //Play with the MC stack if available
 868     //--------------------------------------------------------------------------------------
 869
 870     //Check origin of the candidates
 871     if(IsDataMC()){
 872       aodph.SetTag(GetMCAnalysisUtils()->CheckOrigin(calo->GetLabels(),calo->GetNLabels(),GetReader(), aodph.GetInputFileIndex()));
 873       if(GetDebug() > 0) printf("AliAnaPhoton::MakeAnalysisFillAOD() - Origin of candidate, bit map %d\n",aodph.GetTag());
 874     }//Work with stack also
 875
 876     //--------------------------------------------------------------------------------------
 877     // Conversions pairs analysis
 878     // Check if cluster comes from a conversion in the material in front of the calorimeter
 879     // Do invariant mass of all pairs, if mass is close to 0, then it is conversion.
 880     //--------------------------------------------------------------------------------------
 881
 882     // Do analysis only if there are more than one cluster
 883     if( nCaloClusters > 1){
 884       Bool_t bConverted = kFALSE;
 885       Int_t id2 = -1;
 886
 887       //Check if set previously as converted couple, if so skip its use.
 888       if (fCheckConversion && indexConverted[icalo]) continue;
 889
 890       // Second cluster loop
 891       for(Int_t jcalo = icalo + 1 ; jcalo < nCaloClusters ; jcalo++) {
 892         //Check if set previously as converted couple, if so skip its use.
 893         if (indexConverted[jcalo]) continue;
 894         //printf("Check Conversion indeces %d and %d\n",icalo,jcalo);
 895         AliVCluster * calo2 =  (AliVCluster*) (pl->At(jcalo));              //Get cluster kinematics
 896
 897         //Mixed event, get index of event
 898         Int_t evtIndex2 = 0 ;
 899         if (GetMixedEvent()) {
 900           evtIndex2=GetMixedEvent()->EventIndexForCaloCluster(calo2->GetID()) ;
 901
 902         }
 903
 904         //Get kinematics of second cluster
 905         if(GetReader()->GetDataType() != AliCaloTrackReader::kMC){
 906           calo2->GetMomentum(mom2,GetVertex(evtIndex2)) ;}//Assume that come from vertex in straight line
 907         else{
 908           Double_t vertex[]={0,0,0};
 909           calo2->GetMomentum(mom2,vertex) ;
 910         }
 911
 912         //Check only certain regions
 913         Bool_t in2 = kTRUE;
 914         if(IsFiducialCutOn()) in2 =  GetFiducialCut()->IsInFiducialCut(mom2,fCalorimeter) ;
 915         if(!in2) continue;
 916
 917         //................................................
 918         //Get mass of pair, if small, take this pair.
 919         //printf("\t both in calo, mass %f, cut %f\n",(mom+mom2).M(),fMassCut);
 920         if((mom+mom2).M() < fMassCut){
 921
 922           id2 = calo2->GetID();
 923           indexConverted[icalo]=kTRUE;
 924           indexConverted[jcalo]=kTRUE;
 925
 926           if(GetDebug() > 2)
 927             printf("AliAnaPhoton::MakeAnalysisFillAOD(): Pair with mass %f < %2.4f; \n    cluster1 id %d, e %2.3f  SM %d, eta %2.3f, phi %2.3f ; \n    cluster2 id %d, e %2.3f, SM %d,eta %2.3f, phi %2.3f\n",
 928                    (mom+mom2).M(),fMassCut,
 929                    calo->GetID(),calo->E(),GetCaloUtils()->GetModuleNumber(calo), mom.Eta(), mom.Phi(),
 930                    id2, calo2->E(), GetCaloUtils()->GetModuleNumber(calo2),mom2.Eta(), mom2.Phi());
 931
 932           //...............................................
 933           //Fill few histograms with kinematics of the pair
 934           //FIXME, move all this to MakeAnalysisFillHistograms ...
 935           fhConvDeltaEta   ->Fill( (mom+mom2).M(), mom.Eta()-mom2.Eta() );
 936           fhConvDeltaPhi   ->Fill( (mom+mom2).M(), mom.Phi()-mom2.Phi() );
 937           fhConvAsym       ->Fill( (mom+mom2).M(), TMath::Abs(mom.E()-mom2.E())/(mom.E()+mom2.E()) );
 938           fhConvDeltaEtaPhi->Fill( mom.Eta()-mom2.Eta(), mom.Phi()-mom2.Phi() );
 939           fhConvPt         ->Fill( (mom+mom2).M(), (mom+mom2).Pt());
 940
 941           //...............................................
 942           //Select pairs in a eta-phi window
 943           if(TMath::Abs(mom.Eta()-mom2.Eta()) < fConvDEtaCut && TMath::Abs(mom.Phi()-mom2.Phi()) < fConvDPhiCut )
 944             bConverted = kTRUE;
 945
 946           //...........................................
 947           //Fill more histograms, simulated data
 948           //FIXME, move all this to MakeAnalysisFillHistograms ...
 949           if(IsDataMC()){
 950
 951             //Check the origin of the pair, look for conversion, antinucleons or jet correlations (strings)
 952             Int_t ancPDG    = 0;
 953             Int_t ancStatus = 0;
 954             TLorentzVector momentum;
 955             Int_t ancLabel  = GetMCAnalysisUtils()->CheckCommonAncestor(calo->GetLabel(), calo2->GetLabel(),
 956                                                                         GetReader(), ancPDG, ancStatus, momentum);
 957
 958             // printf("AliAnaPhoton::MakeAnalysisFillHistograms() - Common ancestor label %d, pdg %d, name %s, status %d; \n",
 959             //                          ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus);
 960
 961             Int_t tag2 = GetMCAnalysisUtils()->CheckOrigin(calo2->GetLabels(),calo2->GetNLabels(),GetReader(), 0);
 962             if(GetMCAnalysisUtils()->CheckTagBit(aodph.GetTag(),AliMCAnalysisUtils::kMCConversion)){
 963               if(GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCConversion) && (ancPDG==22 || TMath::Abs(ancPDG)==11) && ancLabel > -1){
 964                 fhConvDeltaEtaMCConversion   ->Fill( (mom+mom2).M(), mom.Eta()-mom2.Eta());
 965                 fhConvDeltaPhiMCConversion   ->Fill( (mom+mom2).M(), mom.Phi()-mom2.Phi());
 966                 fhConvAsymMCConversion       ->Fill( (mom+mom2).M(), TMath::Abs(mom.E()-mom2.E())/(mom.E()+mom2.E()) );
 967                 fhConvDeltaEtaPhiMCConversion->Fill( mom.Eta()-mom2.Eta(), mom.Phi()-mom2.Phi() );
 968                 fhConvPtMCConversion         ->Fill( (mom+mom2).M(), (mom+mom2).Pt());
 969                 fhConvDispersionMCConversion ->Fill( calo->GetDispersion(), calo2->GetDispersion());
 970                 fhConvM02MCConversion        ->Fill( calo->GetM02(), calo2->GetM02());
 971
 972               }
 973             }
 974             else if(GetMCAnalysisUtils()->CheckTagBit(aodph.GetTag(),AliMCAnalysisUtils::kMCAntiNeutron)){
 975               if(GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCAntiNeutron) && ancPDG==-2112 && ancLabel > -1){
 976                 fhConvDeltaEtaMCAntiNeutron    ->Fill( (mom+mom2).M(), mom.Eta()-mom2.Eta());
 977                 fhConvDeltaPhiMCAntiNeutron    ->Fill( (mom+mom2).M(), mom.Phi()-mom2.Phi());
 978                 fhConvAsymMCAntiNeutron        ->Fill( (mom+mom2).M(), TMath::Abs(mom.E()-mom2.E())/(mom.E()+mom2.E()) );
 979                 fhConvDeltaEtaPhiMCAntiNeutron ->Fill( mom.Eta()-mom2.Eta(), mom.Phi()-mom2.Phi() );
 980                 fhConvPtMCAntiNeutron          ->Fill( (mom+mom2).M(), (mom+mom2).Pt());
 981                 fhConvDispersionMCAntiNeutron  ->Fill( calo->GetDispersion(), calo2->GetDispersion());
 982                 fhConvM02MCAntiNeutron         ->Fill( calo->GetM02(), calo2->GetM02());
 983               }
 984             }
 985             else if(GetMCAnalysisUtils()->CheckTagBit(aodph.GetTag(),AliMCAnalysisUtils::kMCAntiProton)){
 986               if(GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCAntiProton) && ancPDG==-2212 && ancLabel > -1){
 987                 fhConvDeltaEtaMCAntiProton    ->Fill( (mom+mom2).M(), mom.Eta()-mom2.Eta());
 988                 fhConvDeltaPhiMCAntiProton    ->Fill( (mom+mom2).M(), mom.Phi()-mom2.Phi());
 989                 fhConvAsymMCAntiProton        ->Fill( (mom+mom2).M(), TMath::Abs(mom.E()-mom2.E())/(mom.E()+mom2.E()) );
 990                 fhConvDeltaEtaPhiMCAntiProton ->Fill( mom.Eta()-mom2.Eta(), mom.Phi()-mom2.Phi() );
 991                 fhConvPtMCAntiProton          ->Fill( (mom+mom2).M(), (mom+mom2).Pt());
 992                 fhConvDispersionMCAntiProton  ->Fill( calo->GetDispersion(), calo2->GetDispersion());
 993                 fhConvM02MCAntiProton         ->Fill( calo->GetM02(), calo2->GetM02());
 994               }
 995             }
 996
 997             //Pairs coming from fragmenting pairs.
 998             if(ancPDG < 22 && ancLabel > 7 && (ancStatus == 11 || ancStatus == 12) ){
 999               fhConvDeltaEtaMCString    ->Fill( (mom+mom2).M(), mom.Eta()-mom2.Eta());
1000               fhConvDeltaPhiMCString    ->Fill( (mom+mom2).M(), mom.Phi()-mom2.Phi());
1001               fhConvAsymMCString        ->Fill( (mom+mom2).M(), TMath::Abs(mom.E()-mom2.E())/(mom.E()+mom2.E()) );
1002               fhConvDeltaEtaPhiMCString ->Fill( mom.Eta()-mom2.Eta(), mom.Phi()-mom2.Phi() );
1003               fhConvPtMCString          ->Fill( (mom+mom2).M(), (mom+mom2).Pt());
1004               fhConvDispersionMCString  ->Fill( calo->GetDispersion(), calo2->GetDispersion());
1005               fhConvM02MCString         ->Fill( calo->GetM02(), calo2->GetM02());
1006             }
1007
1008           }// Data MC
1009
1010           break;
1011         }
1012
1013       }//Mass loop
1014
1015       //..........................................................................................................
1016       //Pair selected as converted, remove both clusters or recombine them into a photon and put them in the AOD
1017       if(bConverted){
1018         //Add to AOD
1019         if(fAddConvertedPairsToAOD){
1020           //Create AOD of pair analysis
1021           TLorentzVector mpair = mom+mom2;
1022           AliAODPWG4Particle aodpair = AliAODPWG4Particle(mpair);
1023           aodpair.SetLabel(aodph.GetLabel());
1024           //aodpair.SetInputFileIndex(input);
1025
1026           //printf("Index %d, Id %d\n",icalo, calo->GetID());
1027           //Set the indeces of the original caloclusters
1028           aodpair.SetCaloLabel(calo->GetID(),id2);
1029           aodpair.SetDetector(fCalorimeter);
1030           aodpair.SetPdg(aodph.GetPdg());
1031           aodpair.SetTag(aodph.GetTag());
1032           aodpair.SetTagged(kTRUE);
1033           //Add AOD with pair object to aod branch
1034           AddAODParticle(aodpair);
1035           //printf("\t \t both added pair\n");
1036         }
1037
1038         //Do not add the current calocluster
1039         if(fCheckConversion) continue;
1040         else {
1041           //printf("TAGGED\n");
1042           //Tag this cluster as likely conversion
1043           aodph.SetTagged(kTRUE);
1044         }
1045       }//converted pair
1046     }//check conversion
1047     //printf("\t \t added single cluster %d\n",icalo);
1048
1049     //FIXME, this to MakeAnalysisFillHistograms ...
1050     fhNCellsPt->Fill(aodph.Pt(),calo->GetNCells());
1051
1052     //Add AOD with photon object to aod branch
1053     AddAODParticle(aodph);
1054
1055   }//loop
1056
1057   delete [] indexConverted;
1058
1059   if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD()  End fill AODs, with %d entries \n",GetOutputAODBranch()->GetEntriesFast());
1060
1061 }
1062
1063 //__________________________________________________________________
1064 void  AliAnaPhoton::MakeAnalysisFillHistograms()
1065 {
1066   //Fill histograms
1067
1068   //-------------------------------------------------------------------
1069         // Access MC information in stack if requested, check that it exists.
1070         AliStack * stack = 0x0;
1071         TParticle * primary = 0x0;
1072         TClonesArray * mcparticles0 = 0x0;
1073         //TClonesArray * mcparticles1 = 0x0;
1074         AliAODMCParticle * aodprimary = 0x0;
1075         if(IsDataMC()){
1076
1077                 if(GetReader()->ReadStack()){
1078                         stack =  GetMCStack() ;
1079                         if(!stack) {
1080                                 printf("AliAnaPhoton::MakeAnalysisFillHistograms() - Stack not available, is the MC handler called? STOP\n");
1081                                 abort();
1082                         }
1083
1084                 }
1085                 else if(GetReader()->ReadAODMCParticles()){
1086
1087                         //Get the list of MC particles
1088                         mcparticles0 = GetReader()->GetAODMCParticles(0);
1089                         if(!mcparticles0 && GetDebug() > 0)     {
1090                                 printf("AliAnaPhoton::MakeAnalysisFillHistograms() -  Standard MCParticles not available!\n");
1091                         }
1092       //                        if(GetReader()->GetSecondInputAODTree()){
1093       //                                mcparticles1 = GetReader()->GetAODMCParticles(1);
1094       //                                if(!mcparticles1 && GetDebug() > 0)     {
1095       //                                        printf("AliAnaPhoton::MakeAnalysisFillHistograms() -  Second input MCParticles not available!\n");
1096       //                                }
1097       //                        }
1098
1099                 }
1100         }// is data and MC
1101
1102
1103   // Get vertex
1104   Double_t v[3] = {0,0,0}; //vertex ;
1105   GetReader()->GetVertex(v);
1106   //fhVertex->Fill(v[0],v[1],v[2]);
1107   if(TMath::Abs(v[2]) > GetZvertexCut()) return ; // done elsewhere for Single Event analysis, but there for mixed event
1108
1109   //----------------------------------
1110         //Loop on stored AOD photons
1111         Int_t naod = GetOutputAODBranch()->GetEntriesFast();
1112   fhNtraNclu->Fill(GetReader()->GetTrackMultiplicity(), naod);
1113         if(GetDebug() > 0) printf("AliAnaPhoton::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod);
1114
1115         for(Int_t iaod = 0; iaod < naod ; iaod++){
1116           AliAODPWG4Particle* ph =  (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
1117           Int_t pdg = ph->GetPdg();
1118
1119           if(GetDebug() > 3)
1120             printf("AliAnaPhoton::MakeAnalysisFillHistograms() - PDG %d, MC TAG %d, Calorimeter %s\n", ph->GetPdg(),ph->GetTag(), (ph->GetDetector()).Data()) ;
1121
1122           //If PID used, fill histos with photons in Calorimeter fCalorimeter
1123           if(IsCaloPIDOn() && pdg != AliCaloPID::kPhoton) continue;
1124           if(ph->GetDetector() != fCalorimeter) continue;
1125
1126           if(GetDebug() > 2)
1127             printf("AliAnaPhoton::MakeAnalysisFillHistograms() - ID Photon: pt %f, phi %f, eta %f\n", ph->Pt(),ph->Phi(),ph->Eta()) ;
1128
1129     //................................
1130           //Fill photon histograms
1131           Float_t ptcluster  = ph->Pt();
1132           Float_t phicluster = ph->Phi();
1133           Float_t etacluster = ph->Eta();
1134           Float_t ecluster   = ph->E();
1135
1136           fhPtPhoton  ->Fill(ptcluster);
1137     if(ph->IsTagged())fhPtPhotonConv->Fill(ptcluster);
1138           fhPhiPhoton ->Fill(ptcluster,phicluster);
1139           fhEtaPhoton ->Fill(ptcluster,etacluster);
1140     if(ptcluster > 0.5){
1141       fhEtaPhiPhoton ->Fill(etacluster, phicluster);
1142       if(ph->IsTagged())fhEtaPhiPhotonConv->Fill(etacluster, phicluster);
1143     }
1144     else {
1145       fhEtaPhi05Photon ->Fill(etacluster, phicluster);
1146       if(ph->IsTagged())fhEtaPhi05PhotonConv->Fill(etacluster, phicluster);
1147     }
1148
1149     //.......................................
1150           //Play with the MC data if available
1151           if(IsDataMC()){
1152
1153             Int_t tag =ph->GetTag();
1154
1155             if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton))
1156       {
1157         fhPtMCPhoton  ->Fill(ptcluster);
1158         fhPhiMCPhoton ->Fill(ptcluster,phicluster);
1159         fhEtaMCPhoton ->Fill(ptcluster,etacluster);
1160
1161         if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion))
1162         {
1163           fhPtConversion  ->Fill(ptcluster);
1164           fhPhiConversion ->Fill(ptcluster,phicluster);
1165           fhEtaConversion ->Fill(ptcluster,etacluster);
1166           if(ph->IsTagged()) fhPtConversionTagged ->Fill(ptcluster);
1167           if(ptcluster > 0.5)fhEtaPhiConversion   ->Fill(etacluster,phicluster);
1168           else               fhEtaPhi05Conversion ->Fill(etacluster,phicluster);
1169         }
1170
1171         if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt)){
1172           fhPtPrompt  ->Fill(ptcluster);
1173           fhPhiPrompt ->Fill(ptcluster,phicluster);
1174           fhEtaPrompt ->Fill(ptcluster,etacluster);
1175         }
1176         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation))
1177         {
1178           fhPtFragmentation  ->Fill(ptcluster);
1179           fhPhiFragmentation ->Fill(ptcluster,phicluster);
1180           fhEtaFragmentation ->Fill(ptcluster,etacluster);
1181         }
1182         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCISR))
1183         {
1184           fhPtISR  ->Fill(ptcluster);
1185           fhPhiISR ->Fill(ptcluster,phicluster);
1186           fhEtaISR ->Fill(ptcluster,etacluster);
1187         }
1188         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay))
1189         {
1190           fhPtPi0Decay  ->Fill(ptcluster);
1191           fhPhiPi0Decay ->Fill(ptcluster,phicluster);
1192           fhEtaPi0Decay ->Fill(ptcluster,etacluster);
1193         }
1194         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay) || GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay))
1195         {
1196           fhPtOtherDecay  ->Fill(ptcluster);
1197           fhPhiOtherDecay ->Fill(ptcluster,phicluster);
1198           fhEtaOtherDecay ->Fill(ptcluster,etacluster);
1199         }
1200       }
1201       else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiNeutron))
1202       {
1203         fhPtAntiNeutron  ->Fill(ptcluster);
1204         fhPhiAntiNeutron ->Fill(ptcluster,phicluster);
1205         fhEtaAntiNeutron ->Fill(ptcluster,etacluster);
1206       }
1207       else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiProton))
1208       {
1209         fhPtAntiProton  ->Fill(ptcluster);
1210         fhPhiAntiProton ->Fill(ptcluster,phicluster);
1211         fhEtaAntiProton ->Fill(ptcluster,etacluster);
1212       }
1213             else{
1214               fhPtUnknown  ->Fill(ptcluster);
1215               fhPhiUnknown ->Fill(ptcluster,phicluster);
1216               fhEtaUnknown ->Fill(ptcluster,etacluster);
1217
1218         //               printf(" AliAnaPhoton::MakeAnalysisFillHistograms() - Label %d, pT %2.3f Unknown, bits set: ",
1219         //                                      ph->GetLabel(),ph->Pt());
1220         //                for(Int_t i = 0; i < 20; i++) {
1221         //                        if(GetMCAnalysisUtils()->CheckTagBit(tag,i)) printf(" %d, ",i);
1222         //                }
1223         //                printf("\n");
1224
1225             }
1226
1227             //....................................................................
1228             // Access MC information in stack if requested, check that it exists.
1229             Int_t label =ph->GetLabel();
1230             if(label < 0) {
1231               printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** bad label ***:  label %d \n", label);
1232               continue;
1233             }
1234
1235             Float_t eprim   = 0;
1236             Float_t ptprim  = 0;
1237             if(GetReader()->ReadStack()){
1238
1239               if(label >=  stack->GetNtrack()) {
1240           if(GetDebug() > 2)  printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** large label ***:  label %d, n tracks %d \n", label, stack->GetNtrack());
1241           continue ;
1242               }
1243
1244               primary = stack->Particle(label);
1245               if(!primary){
1246           printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** no primary ***:  label %d \n", label);
1247           continue;
1248               }
1249               eprim   = primary->Energy();
1250               ptprim  = primary->Pt();
1251
1252             }
1253             else if(GetReader()->ReadAODMCParticles()){
1254               //Check which is the input
1255               if(ph->GetInputFileIndex() == 0){
1256           if(!mcparticles0) continue;
1257           if(label >=  mcparticles0->GetEntriesFast()) {
1258             if(GetDebug() > 2)  printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** large label ***:  label %d, n tracks %d \n",
1259                                        label, mcparticles0->GetEntriesFast());
1260             continue ;
1261           }
1262           //Get the particle
1263           aodprimary = (AliAODMCParticle*) mcparticles0->At(label);
1264
1265               }
1266 //            else {//Second input
1267 //          if(!mcparticles1) continue;
1268 //          if(label >=  mcparticles1->GetEntriesFast()) {
1269 //            if(GetDebug() > 2)  printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** large label ***:  label %d, n tracks %d \n",
1270 //                                       label, mcparticles1->GetEntriesFast());
1271 //            continue ;
1272 //          }
1273 //          //Get the particle
1274 //          aodprimary = (AliAODMCParticle*) mcparticles1->At(label);
1275 //
1276 //            }//second input
1277
1278               if(!aodprimary){
1279           printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** no primary ***:  label %d \n", label);
1280           continue;
1281               }
1282
1283               eprim   = aodprimary->E();
1284               ptprim  = aodprimary->Pt();
1285
1286             }
1287
1288             fh2E     ->Fill(ecluster, eprim);
1289             fh2Pt    ->Fill(ptcluster, ptprim);
1290             fhDeltaE ->Fill(eprim-ecluster);
1291             fhDeltaPt->Fill(ptprim-ptcluster);
1292             if(eprim > 0)  fhRatioE  ->Fill(ecluster/eprim);
1293             if(ptprim > 0) fhRatioPt ->Fill(ptcluster/ptprim);
1294
1295           }//Histograms with MC
1296
1297         }// aod loop
1298
1299 }
1300
1301
1302 //__________________________________________________________________
1303 void AliAnaPhoton::Print(const Option_t * opt) const
1304 {
1305   //Print some relevant parameters set for the analysis
1306
1307   if(! opt)
1308     return;
1309
1310   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
1311   AliAnaPartCorrBaseClass::Print(" ");
1312
1313   printf("Calorimeter            =     %s\n", fCalorimeter.Data()) ;
1314   printf("Min Distance to Bad Channel   = %2.1f\n",fMinDist);
1315   printf("Min Distance to Bad Channel 2 = %2.1f\n",fMinDist2);
1316   printf("Min Distance to Bad Channel 3 = %2.1f\n",fMinDist3);
1317   printf("Reject clusters with a track matched = %d\n",fRejectTrackMatch);
1318   printf("Check Pair Conversion                = %d\n",fCheckConversion);
1319   printf("Add conversion pair to AOD           = %d\n",fAddConvertedPairsToAOD);
1320   printf("Conversion pair mass cut             = %f\n",fMassCut);
1321   printf("Time Cut: %3.1f < TOF  < %3.1f\n", fTimeCutMin, fTimeCutMax);
1322   printf("Number of cells in cluster is        > %d \n", fNCellsCut);
1323   printf("    \n") ;
1324
1325 }