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