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