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, 100,0,100);
 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   //Store calo PID histograms
 655   if(fRejectTrackMatch){
 656     TList * caloPIDHistos = GetCaloPID()->GetCreateOutputObjects() ;
 657     for(Int_t i = 0; i < caloPIDHistos->GetEntries(); i++) {
 658       outputContainer->Add(caloPIDHistos->At(i)) ;
 659     }
 660     delete caloPIDHistos;
 661   }
 662
 663   return outputContainer ;
 664
 665 }
 666
 667 //____________________________________________________________________________
 668 void AliAnaPhoton::Init()
 669 {
 670
 671   //Init
 672   //Do some checks
 673   if(fCalorimeter == "PHOS" && !GetReader()->IsPHOSSwitchedOn() && NewOutputAOD()){
 674     printf("AliAnaPhoton::Init() - !!STOP: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!\n");
 675     abort();
 676   }
 677   else  if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn() && NewOutputAOD()){
 678     printf("AliAnaPhoton::Init() - !!STOP: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!\n");
 679     abort();
 680   }
 681
 682 }
 683
 684
 685 //____________________________________________________________________________
 686 void AliAnaPhoton::InitParameters()
 687 {
 688
 689   //Initialize the parameters of the analysis.
 690   AddToHistogramsName("AnaPhoton_");
 691
 692   fCalorimeter = "EMCAL" ;
 693   fMinDist     = 2.;
 694   fMinDist2    = 4.;
 695   fMinDist3    = 5.;
 696   fMassCut     = 0.03; //30 MeV
 697
 698   fTimeCutMin  = -1;
 699   fTimeCutMax  = 9999999;
 700   fNCellsCut   = 0;
 701
 702   fRejectTrackMatch       = kTRUE ;
 703   fCheckConversion        = kFALSE;
 704   fRemoveConvertedPair    = kFALSE;
 705   fAddConvertedPairsToAOD = kFALSE;
 706
 707 }
 708
 709 //__________________________________________________________________
 710 void  AliAnaPhoton::MakeAnalysisFillAOD()
 711 {
 712   //Do photon analysis and fill aods
 713
 714   //Get the vertex
 715   Double_t v[3] = {0,0,0}; //vertex ;
 716   GetReader()->GetVertex(v);
 717
 718   //Select the Calorimeter of the photon
 719   TObjArray * pl = 0x0;
 720   if(fCalorimeter == "PHOS")
 721     pl = GetPHOSClusters();
 722   else if (fCalorimeter == "EMCAL")
 723     pl = GetEMCALClusters();
 724
 725   if(!pl) {
 726     Info("MakeAnalysisFillAOD","TObjArray with %s clusters is NULL!\n",fCalorimeter.Data());
 727     return;
 728   }
 729
 730   //Init arrays, variables, get number of clusters
 731   TLorentzVector mom, mom2 ;
 732   Int_t nCaloClusters = pl->GetEntriesFast();
 733   //List to be used in conversion analysis, to tag the cluster as candidate for conversion
 734   Bool_t * indexConverted = 0x0;
 735   if(fCheckConversion){
 736     indexConverted = new Bool_t[nCaloClusters];
 737     for (Int_t i = 0; i < nCaloClusters; i++)
 738       indexConverted[i] = kFALSE;
 739         }
 740
 741   if(GetDebug() > 0) printf("AliAnaPhoton::MakeAnalysisFillAOD() - input %s cluster entries %d\n", fCalorimeter.Data(), nCaloClusters);
 742
 743   //----------------------------------------------------
 744   // Fill AOD with PHOS/EMCAL AliAODPWG4Particle objects
 745   //----------------------------------------------------
 746   // Loop on clusters
 747   for(Int_t icalo = 0; icalo < nCaloClusters; icalo++){
 748
 749           AliVCluster * calo =  (AliVCluster*) (pl->At(icalo));
 750     //printf("calo %d, %f\n",icalo,calo->E());
 751
 752     //Get the index where the cluster comes, to retrieve the corresponding vertex
 753     Int_t evtIndex = 0 ;
 754     if (GetMixedEvent()) {
 755       evtIndex=GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ;
 756       //Get the vertex and check it is not too large in z
 757       if(TMath::Abs(GetVertex(evtIndex)[2])> GetZvertexCut()) continue;
 758     }
 759
 760     //Cluster selection, not charged, with photon id and in fiducial cut
 761
 762     //Input from second AOD?
 763     //Int_t input = 0;
 764     //    if (fCalorimeter == "EMCAL" && GetReader()->GetEMCALClustersNormalInputEntries() <= icalo)
 765     //      input = 1 ;
 766     //    else if(fCalorimeter == "PHOS"  && GetReader()->GetPHOSClustersNormalInputEntries()  <= icalo)
 767     //      input = 1;
 768
 769     //Get Momentum vector,
 770     //if (input == 0)
 771     if(GetReader()->GetDataType() != AliCaloTrackReader::kMC){
 772       calo->GetMomentum(mom,GetVertex(evtIndex)) ;}//Assume that come from vertex in straight line
 773     else{
 774       Double_t vertex[]={0,0,0};
 775       calo->GetMomentum(mom,vertex) ;
 776     }
 777
 778     //    else if(input == 1)
 779     //      calo->GetMomentum(mom,vertex2);//Assume that come from vertex in straight line
 780
 781     //--------------------------------------
 782     // Cluster selection
 783     //--------------------------------------
 784     if(GetDebug() > 2)
 785       printf("AliAnaPhoton::MakeAnalysisFillAOD() Current Event %d; Before selection : E %2.2f, pT %2.2f, Ecl %2.2f, phi %2.2f, eta %2.2f\n",
 786              GetReader()->GetEventNumber(),
 787              mom.E(), mom.Pt(),calo->E(),mom.Phi()*TMath::RadToDeg(),mom.Eta());
 788
 789     //.......................................
 790     //If too small or big pt, skip it
 791     if(mom.E() < GetMinPt() || mom.E() > GetMaxPt() ) continue ;
 792     if(GetDebug() > 2) printf("\t Cluster %d Pass E Cut \n",icalo);
 793
 794     //.......................................
 795     // TOF cut, BE CAREFUL WITH THIS CUT
 796     Double_t tof = calo->GetTOF()*1e9;
 797     if(tof < fTimeCutMin || tof > fTimeCutMax) continue;
 798           if(GetDebug() > 2)  printf("\t Cluster %d Pass Time Cut \n",icalo);
 799
 800     //.......................................
 801     if(calo->GetNCells() <= fNCellsCut && GetReader()->GetDataType() != AliCaloTrackReader::kMC) continue;
 802     if(GetDebug() > 2) printf("\t Cluster %d Pass NCell Cut \n",icalo);
 803
 804     //.......................................
 805     //Check acceptance selection
 806     if(IsFiducialCutOn()){
 807       Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter) ;
 808       if(! in ) continue ;
 809     }
 810     if(GetDebug() > 2) printf("Fiducial cut passed \n");
 811
 812     //.......................................
 813     //Skip matched clusters with tracks
 814     if(fRejectTrackMatch){
 815       if(IsTrackMatched(calo)) {
 816         if(GetDebug() > 2) printf("\t Reject track-matched clusters\n");
 817         continue ;
 818       }
 819       else
 820         if(GetDebug() > 2)  printf(" Track-matching cut passed \n");
 821     }// reject matched clusters
 822
 823     //.......................................
 824     //Check Distance to Bad channel, set bit.
 825     Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
 826     if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;
 827     if(distBad < fMinDist) {//In bad channel (PHOS cristal size 2.2x2.2 cm), EMCAL ( cell units )
 828       continue ;
 829     }
 830     else if(GetDebug() > 2) printf("\t Bad channel cut passed %4.2f > %2.2f \n",distBad, fMinDist);
 831
 832     if(GetDebug() > 0)
 833       printf("AliAnaPhoton::MakeAnalysisFillAOD() Current Event %d; After  selection : E %2.2f, pT %2.2f, Ecl %2.2f, phi %2.2f, eta %2.2f\n",
 834              GetReader()->GetEventNumber(),
 835              mom.E(), mom.Pt(),calo->E(),mom.Phi()*TMath::RadToDeg(),mom.Eta());
 836
 837
 838     //----------------------------
 839     //Create AOD for analysis
 840     //----------------------------
 841     AliAODPWG4Particle aodph = AliAODPWG4Particle(mom);
 842
 843     //...............................................
 844     //Set the indeces of the original caloclusters (MC, ID), and calorimeter
 845     Int_t label = calo->GetLabel();
 846     aodph.SetLabel(label);
 847     //aodph.SetInputFileIndex(input);
 848     aodph.SetCaloLabel(calo->GetID(),-1);
 849     aodph.SetDetector(fCalorimeter);
 850     //printf("Index %d, Id %d\n",icalo, calo->GetID());
 851
 852     //printf("Cluster %d Pass Bad Dist Cut \n",icalo);
 853
 854     //...............................................
 855     //Set bad channel distance bit
 856     if     (distBad > fMinDist3) aodph.SetDistToBad(2) ;
 857     else if(distBad > fMinDist2) aodph.SetDistToBad(1) ;
 858     else                         aodph.SetDistToBad(0) ;
 859     //printf("DistBad %f Bit %d\n",distBad, aodph.DistToBad());
 860
 861     //...............................................
 862     //Set number of cells in this cluster
 863     //Temporary patch FIXME
 864     aodph.SetBtag(calo->GetNCells());
 865     // MEFIX
 866
 867     //-------------------------------------
 868     //PID selection or bit setting
 869     //-------------------------------------
 870     // MC
 871     if(GetReader()->GetDataType() == AliCaloTrackReader::kMC){
 872       //Get most probable PID, check PID weights (in MC this option is mandatory)
 873       aodph.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,calo->GetPID(),mom.E()));//PID with weights
 874       if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PDG of identified particle %d\n",aodph.GetPdg());
 875       //If primary is not photon, skip it.
 876       if(aodph.GetPdg() != AliCaloPID::kPhoton) continue ;
 877     }
 878     //...............................................
 879     // Data, PID check on
 880     else if(IsCaloPIDOn()){
 881       //Get most probable PID, 2 options check PID weights
 882       //or redo PID, recommended option for EMCal.
 883       if(!IsCaloPIDRecalculationOn())
 884         aodph.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,calo->GetPID(),mom.E()));//PID with weights
 885       else
 886         aodph.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,mom,calo));//PID recalculated
 887
 888       if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PDG of identified particle %d\n",aodph.GetPdg());
 889
 890       //If cluster does not pass pid, not photon, skip it.
 891       if(aodph.GetPdg() != AliCaloPID::kPhoton) continue ;
 892
 893     }
 894     //...............................................
 895     // Data, PID check off
 896     else{
 897       //Set PID bits for later selection (AliAnaPi0 for example)
 898       //GetPDG already called in SetPIDBits.
 899       GetCaloPID()->SetPIDBits(fCalorimeter,calo,&aodph, GetCaloUtils());
 900       if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PID Bits set \n");
 901     }
 902
 903     if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - Photon selection cuts passed: pT %3.2f, pdg %d\n",aodph.Pt(), aodph.GetPdg());
 904
 905     //--------------------------------------------------------------------------------------
 906     //Play with the MC stack if available
 907     //--------------------------------------------------------------------------------------
 908
 909     //Check origin of the candidates
 910     if(IsDataMC()){
 911       aodph.SetTag(GetMCAnalysisUtils()->CheckOrigin(calo->GetLabels(),calo->GetNLabels(),GetReader(), aodph.GetInputFileIndex()));
 912       if(GetDebug() > 0)
 913         printf("AliAnaPhoton::MakeAnalysisFillAOD() - Origin of candidate, bit map %d\n",aodph.GetTag());
 914     }//Work with stack also
 915
 916     //--------------------------------------------------------------------------------------
 917     // Conversions pairs analysis
 918     // Check if cluster comes from a conversion in the material in front of the calorimeter
 919     // Do invariant mass of all pairs, if mass is close to 0, then it is conversion.
 920     //--------------------------------------------------------------------------------------
 921
 922     // Do analysis only if there are more than one cluster
 923     if( nCaloClusters > 1 && fCheckConversion){
 924       Bool_t bConverted = kFALSE;
 925       Int_t id2 = -1;
 926
 927       //Check if set previously as converted couple, if so skip its use.
 928       if (indexConverted[icalo]) continue;
 929
 930       // Second cluster loop
 931       for(Int_t jcalo = icalo + 1 ; jcalo < nCaloClusters ; jcalo++) {
 932         //Check if set previously as converted couple, if so skip its use.
 933         if (indexConverted[jcalo]) continue;
 934         //printf("Check Conversion indeces %d and %d\n",icalo,jcalo);
 935         AliVCluster * calo2 =  (AliVCluster*) (pl->At(jcalo));              //Get cluster kinematics
 936
 937         //Mixed event, get index of event
 938         Int_t evtIndex2 = 0 ;
 939         if (GetMixedEvent()) {
 940           evtIndex2=GetMixedEvent()->EventIndexForCaloCluster(calo2->GetID()) ;
 941
 942         }
 943
 944         //Get kinematics of second cluster
 945         if(GetReader()->GetDataType() != AliCaloTrackReader::kMC){
 946           calo2->GetMomentum(mom2,GetVertex(evtIndex2)) ;}//Assume that come from vertex in straight line
 947         else{
 948           Double_t vertex[]={0,0,0};
 949           calo2->GetMomentum(mom2,vertex) ;
 950         }
 951
 952         //Check only certain regions
 953         Bool_t in2 = kTRUE;
 954         if(IsFiducialCutOn()) in2 =  GetFiducialCut()->IsInFiducialCut(mom2,fCalorimeter) ;
 955         if(!in2) continue;
 956
 957         //................................................
 958         //Get mass of pair, if small, take this pair.
 959         Float_t pairM     = (mom+mom2).M();
 960         //printf("\t both in calo, mass %f, cut %f\n",pairM,fMassCut);
 961         if(pairM < fMassCut){
 962           aodph.SetTagged(kFALSE);
 963           id2 = calo2->GetID();
 964           indexConverted[icalo]=kTRUE;
 965           indexConverted[jcalo]=kTRUE;
 966
 967           Float_t asymmetry = TMath::Abs(mom.E()-mom2.E())/(mom.E()+mom2.E());
 968           Float_t dPhi      = mom.Phi()-mom2.Phi();
 969           Float_t dEta      = mom.Eta()-mom2.Eta();
 970
 971           //Estimate conversion distance, T. Awes, M. Ivanov
 972           //Under the assumption that the pair has zero mass, and that each electron
 973           //of the pair has the same momentum, they will each have the same bend radius
 974           //given by R=p/(qB) = p / (300 B) with p in [MeV/c], B in [Tesla] and R in [m].
 975           //With nominal ALICE magnet current of 30kA B=0.5T, and so with E_cluster=p,
 976           //R = E/1.5 [cm].  Under these assumptions, the distance from the conversion
 977           //point to the EMCal can be related to the separation distance, L=2y, on the EMCal
 978           //as d = sqrt(R^2 -(R-y)^2) = sqrt(2Ry - y^2). And since R>>y we can write as
 979           //d = sqrt(E*L/1.5) where E is the cluster energy and L is the distance in cm between
 980           //the clusters.
 981 //          Float_t pos1[3];
 982 //          calo->GetPosition(pos1);
 983 //          Float_t pos2[3];
 984 //          calo2->GetPosition(pos2);
 985 //          Float_t l = TMath::Sqrt((pos1[0]-pos2[0])*(pos1[0]-pos2[0])+
 986 //                                  (pos1[1]-pos2[1])*(pos1[1]-pos2[1])+
 987 //                                  (pos1[2]-pos2[2])*(pos1[2]-pos2[2]));
 988 //
 989 //          Float_t convDist  = TMath::Sqrt(mom.E() *l/1.5);
 990 //          Float_t convDist2 = TMath::Sqrt(mom2.E()*l/1.5);
 991 //          printf("l = %f, e1 = %f, d1=%f, e2 = %f, d2=%f\n",l,mom.E(),convDist,mom2.E(),convDist2);
 992           if(GetDebug() > 2)
 993             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",
 994                    pairM,fMassCut,fConvDPhiMinCut, dPhi, fConvDPhiMaxCut, dEta, fConvDEtaCut, asymmetry, fConvAsymCut,
 995                    calo->GetID(),calo->E(),GetCaloUtils()->GetModuleNumber(calo), mom.Eta(), mom.Phi(),
 996                    id2, calo2->E(), GetCaloUtils()->GetModuleNumber(calo2),mom2.Eta(), mom2.Phi());
 997
 998           //...............................................
 999           //Fill few histograms with kinematics of the pair
1000           //FIXME, move all this to MakeAnalysisFillHistograms ...
1001
1002           fhConvDeltaEta   ->Fill( pairM, dPhi      );
1003           fhConvDeltaPhi   ->Fill( pairM, dEta      );
1004           fhConvAsym       ->Fill( pairM, asymmetry );
1005           fhConvDeltaEtaPhi->Fill( dEta , dPhi      );
1006           fhConvPt         ->Fill( pairM, (mom+mom2).Pt());
1007
1008           //...............................................
1009           //Select pairs in a eta-phi window
1010           if(TMath::Abs(dEta) < fConvDEtaCut    &&
1011              TMath::Abs(dPhi) < fConvDPhiMaxCut &&
1012              TMath::Abs(dPhi) > fConvDPhiMinCut &&
1013              asymmetry        < fConvAsymCut       ){
1014               bConverted = kTRUE;
1015           }
1016           //printf("Accepted? %d\n",bConverted);
1017           //...........................................
1018           //Fill more histograms, simulated data
1019           //FIXME, move all this to MakeAnalysisFillHistograms ...
1020           if(IsDataMC()){
1021
1022             //Check the origin of the pair, look for conversion, antinucleons or jet correlations (strings)
1023             Int_t ancPDG    = 0;
1024             Int_t ancStatus = 0;
1025             TLorentzVector momentum;
1026             Int_t ancLabel  = GetMCAnalysisUtils()->CheckCommonAncestor(calo->GetLabel(), calo2->GetLabel(),
1027                                                                         GetReader(), ancPDG, ancStatus, momentum);
1028
1029             // printf("AliAnaPhoton::MakeAnalysisFillHistograms() - Common ancestor label %d, pdg %d, name %s, status %d; \n",
1030             //                          ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus);
1031
1032             Int_t tag2 = GetMCAnalysisUtils()->CheckOrigin(calo2->GetLabels(),calo2->GetNLabels(),GetReader(), 0);
1033             if(GetMCAnalysisUtils()->CheckTagBit(aodph.GetTag(),AliMCAnalysisUtils::kMCConversion)){
1034               if(GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCConversion) && (ancPDG==22 || TMath::Abs(ancPDG)==11) && ancLabel > -1){
1035                 fhConvDeltaEtaMCConversion   ->Fill( pairM, dEta      );
1036                 fhConvDeltaPhiMCConversion   ->Fill( pairM, dPhi      );
1037                 fhConvAsymMCConversion       ->Fill( pairM, asymmetry );
1038                 fhConvDeltaEtaPhiMCConversion->Fill( dEta , dPhi      );
1039                 fhConvPtMCConversion         ->Fill( pairM, (mom+mom2).Pt());
1040                 fhConvDispersionMCConversion ->Fill( calo->GetDispersion(), calo2->GetDispersion());
1041                 fhConvM02MCConversion        ->Fill( calo->GetM02(), calo2->GetM02());
1042
1043               }
1044             }
1045             else if(GetMCAnalysisUtils()->CheckTagBit(aodph.GetTag(),AliMCAnalysisUtils::kMCAntiNeutron)){
1046               if(GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCAntiNeutron) && ancPDG==-2112 && ancLabel > -1){
1047                 fhConvDeltaEtaMCAntiNeutron    ->Fill( pairM, dEta      );
1048                 fhConvDeltaPhiMCAntiNeutron    ->Fill( pairM, dPhi      );
1049                 fhConvAsymMCAntiNeutron        ->Fill( pairM, asymmetry );
1050                 fhConvDeltaEtaPhiMCAntiNeutron ->Fill( dEta , dPhi      );
1051                 fhConvPtMCAntiNeutron          ->Fill( pairM, (mom+mom2).Pt());
1052                 fhConvDispersionMCAntiNeutron  ->Fill( calo->GetDispersion(), calo2->GetDispersion());
1053                 fhConvM02MCAntiNeutron         ->Fill( calo->GetM02(), calo2->GetM02());
1054               }
1055             }
1056             else if(GetMCAnalysisUtils()->CheckTagBit(aodph.GetTag(),AliMCAnalysisUtils::kMCAntiProton)){
1057               if(GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCAntiProton) && ancPDG==-2212 && ancLabel > -1){
1058                 fhConvDeltaEtaMCAntiProton    ->Fill( pairM, dEta      );
1059                 fhConvDeltaPhiMCAntiProton    ->Fill( pairM, dPhi      );
1060                 fhConvAsymMCAntiProton        ->Fill( pairM, asymmetry );
1061                 fhConvDeltaEtaPhiMCAntiProton ->Fill( dEta , dPhi      );
1062                 fhConvPtMCAntiProton          ->Fill( pairM, (mom+mom2).Pt());
1063                 fhConvDispersionMCAntiProton  ->Fill( calo->GetDispersion(), calo2->GetDispersion());
1064                 fhConvM02MCAntiProton         ->Fill( calo->GetM02(), calo2->GetM02());
1065               }
1066             }
1067
1068             //Pairs coming from fragmenting pairs.
1069             if(ancPDG < 22 && ancLabel > 7 && (ancStatus == 11 || ancStatus == 12) ){
1070               fhConvDeltaEtaMCString    ->Fill( pairM, dPhi);
1071               fhConvDeltaPhiMCString    ->Fill( pairM, dPhi);
1072               fhConvAsymMCString        ->Fill( pairM, TMath::Abs(mom.E()-mom2.E())/(mom.E()+mom2.E()) );
1073               fhConvDeltaEtaPhiMCString ->Fill( dPhi, dPhi );
1074               fhConvPtMCString          ->Fill( pairM, (mom+mom2).Pt());
1075               fhConvDispersionMCString  ->Fill( calo->GetDispersion(), calo2->GetDispersion());
1076               fhConvM02MCString         ->Fill( calo->GetM02(), calo2->GetM02());
1077             }
1078
1079           }// Data MC
1080
1081           break;
1082         }
1083
1084       }//Mass loop
1085
1086       //..........................................................................................................
1087       //Pair selected as converted, remove both clusters or recombine them into a photon and put them in the AOD
1088       if(bConverted){
1089         //Add to AOD
1090         if(fAddConvertedPairsToAOD){
1091           //Create AOD of pair analysis
1092           TLorentzVector mpair = mom+mom2;
1093           AliAODPWG4Particle aodpair = AliAODPWG4Particle(mpair);
1094           aodpair.SetLabel(aodph.GetLabel());
1095           //aodpair.SetInputFileIndex(input);
1096
1097           //printf("Index %d, Id %d\n",icalo, calo->GetID());
1098           //Set the indeces of the original caloclusters
1099           aodpair.SetCaloLabel(calo->GetID(),id2);
1100           aodpair.SetDetector(fCalorimeter);
1101           aodpair.SetPdg(aodph.GetPdg());
1102           aodpair.SetTag(aodph.GetTag());
1103           aodpair.SetTagged(kTRUE);
1104           //Add AOD with pair object to aod branch
1105           AddAODParticle(aodpair);
1106           //printf("\t \t both added pair\n");
1107         }
1108
1109         //Do not add the current calocluster
1110         if(fRemoveConvertedPair) continue;
1111         else {
1112           //printf("TAGGED\n");
1113           //Tag this cluster as likely conversion
1114           aodph.SetTagged(kTRUE);
1115         }
1116       }//converted pair
1117     }//check conversion
1118     //printf("\t \t added single cluster %d\n",icalo);
1119
1120     //FIXME, this to MakeAnalysisFillHistograms ...
1121     fhNCellsPt->Fill(aodph.Pt(),calo->GetNCells());
1122
1123     //Add AOD with photon object to aod branch
1124     AddAODParticle(aodph);
1125
1126   }//loop
1127
1128   delete [] indexConverted;
1129
1130   if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD()  End fill AODs, with %d entries \n",GetOutputAODBranch()->GetEntriesFast());
1131
1132 }
1133
1134 //__________________________________________________________________
1135 void  AliAnaPhoton::MakeAnalysisFillHistograms()
1136 {
1137   //Fill histograms
1138
1139   //-------------------------------------------------------------------
1140         // Access MC information in stack if requested, check that it exists.
1141         AliStack * stack = 0x0;
1142         TParticle * primary = 0x0;
1143         TClonesArray * mcparticles0 = 0x0;
1144         //TClonesArray * mcparticles1 = 0x0;
1145         AliAODMCParticle * aodprimary = 0x0;
1146         if(IsDataMC()){
1147
1148                 if(GetReader()->ReadStack()){
1149                         stack =  GetMCStack() ;
1150                         if(!stack) {
1151                                 printf("AliAnaPhoton::MakeAnalysisFillHistograms() - Stack not available, is the MC handler called? STOP\n");
1152                                 abort();
1153                         }
1154
1155                 }
1156                 else if(GetReader()->ReadAODMCParticles()){
1157
1158                         //Get the list of MC particles
1159                         mcparticles0 = GetReader()->GetAODMCParticles(0);
1160                         if(!mcparticles0 && GetDebug() > 0)     {
1161                                 printf("AliAnaPhoton::MakeAnalysisFillHistograms() -  Standard MCParticles not available!\n");
1162                         }
1163       //                        if(GetReader()->GetSecondInputAODTree()){
1164       //                                mcparticles1 = GetReader()->GetAODMCParticles(1);
1165       //                                if(!mcparticles1 && GetDebug() > 0)     {
1166       //                                        printf("AliAnaPhoton::MakeAnalysisFillHistograms() -  Second input MCParticles not available!\n");
1167       //                                }
1168       //                        }
1169
1170                 }
1171         }// is data and MC
1172
1173
1174   // Get vertex
1175   Double_t v[3] = {0,0,0}; //vertex ;
1176   GetReader()->GetVertex(v);
1177   //fhVertex->Fill(v[0],v[1],v[2]);
1178   if(TMath::Abs(v[2]) > GetZvertexCut()) return ; // done elsewhere for Single Event analysis, but there for mixed event
1179
1180   //----------------------------------
1181         //Loop on stored AOD photons
1182         Int_t naod = GetOutputAODBranch()->GetEntriesFast();
1183   fhNtraNclu->Fill(GetReader()->GetTrackMultiplicity(), naod);
1184         if(GetDebug() > 0) printf("AliAnaPhoton::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod);
1185
1186         for(Int_t iaod = 0; iaod < naod ; iaod++){
1187           AliAODPWG4Particle* ph =  (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
1188           Int_t pdg = ph->GetPdg();
1189
1190           if(GetDebug() > 3)
1191             printf("AliAnaPhoton::MakeAnalysisFillHistograms() - PDG %d, MC TAG %d, Calorimeter %s\n", ph->GetPdg(),ph->GetTag(), (ph->GetDetector()).Data()) ;
1192
1193           //If PID used, fill histos with photons in Calorimeter fCalorimeter
1194           if(IsCaloPIDOn() && pdg != AliCaloPID::kPhoton) continue;
1195           if(ph->GetDetector() != fCalorimeter) continue;
1196
1197           if(GetDebug() > 2)
1198             printf("AliAnaPhoton::MakeAnalysisFillHistograms() - ID Photon: pt %f, phi %f, eta %f\n", ph->Pt(),ph->Phi(),ph->Eta()) ;
1199
1200     //................................
1201           //Fill photon histograms
1202           Float_t ptcluster  = ph->Pt();
1203           Float_t phicluster = ph->Phi();
1204           Float_t etacluster = ph->Eta();
1205           Float_t ecluster   = ph->E();
1206
1207     fhEPhoton   ->Fill(ecluster);
1208           fhPtPhoton  ->Fill(ptcluster);
1209           fhPhiPhoton ->Fill(ptcluster,phicluster);
1210           fhEtaPhoton ->Fill(ptcluster,etacluster);
1211     if(ecluster > 0.5)         fhEtaPhiPhoton ->Fill(etacluster, phicluster);
1212     else if(GetMinPt() < 0.5) fhEtaPhi05Photon->Fill(etacluster, phicluster);
1213
1214     if(fCheckConversion &&ph->IsTagged()){
1215       fhPtPhotonConv->Fill(ptcluster);
1216       if(ecluster > 0.5)        fhEtaPhiPhotonConv  ->Fill(etacluster, phicluster);
1217       else if(GetMinPt() < 0.5) fhEtaPhi05PhotonConv->Fill(etacluster, phicluster);
1218     }
1219
1220     //.......................................
1221           //Play with the MC data if available
1222           if(IsDataMC()){
1223
1224             Int_t tag =ph->GetTag();
1225
1226             if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton))
1227       {
1228         fhPtMCPhoton  ->Fill(ptcluster);
1229         fhPhiMCPhoton ->Fill(ptcluster,phicluster);
1230         fhEtaMCPhoton ->Fill(ptcluster,etacluster);
1231
1232         if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion))
1233         {
1234           fhPtConversion  ->Fill(ptcluster);
1235           fhPhiConversion ->Fill(ptcluster,phicluster);
1236           fhEtaConversion ->Fill(ptcluster,etacluster);
1237           if(ph->IsTagged()) fhPtConversionTagged ->Fill(ptcluster);
1238           if(ptcluster > 0.5)fhEtaPhiConversion   ->Fill(etacluster,phicluster);
1239           else               fhEtaPhi05Conversion ->Fill(etacluster,phicluster);
1240         }
1241
1242         if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt)){
1243           fhPtPrompt  ->Fill(ptcluster);
1244           fhPhiPrompt ->Fill(ptcluster,phicluster);
1245           fhEtaPrompt ->Fill(ptcluster,etacluster);
1246         }
1247         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation))
1248         {
1249           fhPtFragmentation  ->Fill(ptcluster);
1250           fhPhiFragmentation ->Fill(ptcluster,phicluster);
1251           fhEtaFragmentation ->Fill(ptcluster,etacluster);
1252         }
1253         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCISR))
1254         {
1255           fhPtISR  ->Fill(ptcluster);
1256           fhPhiISR ->Fill(ptcluster,phicluster);
1257           fhEtaISR ->Fill(ptcluster,etacluster);
1258         }
1259         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay))
1260         {
1261           fhPtPi0Decay  ->Fill(ptcluster);
1262           fhPhiPi0Decay ->Fill(ptcluster,phicluster);
1263           fhEtaPi0Decay ->Fill(ptcluster,etacluster);
1264         }
1265         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay) || GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay))
1266         {
1267           fhPtOtherDecay  ->Fill(ptcluster);
1268           fhPhiOtherDecay ->Fill(ptcluster,phicluster);
1269           fhEtaOtherDecay ->Fill(ptcluster,etacluster);
1270         }
1271       }
1272       else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiNeutron))
1273       {
1274         fhPtAntiNeutron  ->Fill(ptcluster);
1275         fhPhiAntiNeutron ->Fill(ptcluster,phicluster);
1276         fhEtaAntiNeutron ->Fill(ptcluster,etacluster);
1277         if(ph->IsTagged() && fCheckConversion) fhPtAntiNeutronTagged ->Fill(ptcluster);
1278
1279       }
1280       else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiProton))
1281       {
1282         fhPtAntiProton  ->Fill(ptcluster);
1283         fhPhiAntiProton ->Fill(ptcluster,phicluster);
1284         fhEtaAntiProton ->Fill(ptcluster,etacluster);
1285         if(ph->IsTagged() && fCheckConversion) fhPtAntiProtonTagged ->Fill(ptcluster);
1286
1287       }
1288             else{
1289               fhPtUnknown  ->Fill(ptcluster);
1290               fhPhiUnknown ->Fill(ptcluster,phicluster);
1291               fhEtaUnknown ->Fill(ptcluster,etacluster);
1292         if(ph->IsTagged() && fCheckConversion) fhPtUnknownTagged ->Fill(ptcluster);
1293
1294
1295         //               printf(" AliAnaPhoton::MakeAnalysisFillHistograms() - Label %d, pT %2.3f Unknown, bits set: ",
1296         //                                      ph->GetLabel(),ph->Pt());
1297         //                for(Int_t i = 0; i < 20; i++) {
1298         //                        if(GetMCAnalysisUtils()->CheckTagBit(tag,i)) printf(" %d, ",i);
1299         //                }
1300         //                printf("\n");
1301
1302             }
1303
1304             //....................................................................
1305             // Access MC information in stack if requested, check that it exists.
1306             Int_t label =ph->GetLabel();
1307             if(label < 0) {
1308               printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** bad label ***:  label %d \n", label);
1309               continue;
1310             }
1311
1312             Float_t eprim   = 0;
1313             Float_t ptprim  = 0;
1314             if(GetReader()->ReadStack()){
1315
1316               if(label >=  stack->GetNtrack()) {
1317           if(GetDebug() > 2)  printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** large label ***:  label %d, n tracks %d \n", label, stack->GetNtrack());
1318           continue ;
1319               }
1320
1321               primary = stack->Particle(label);
1322               if(!primary){
1323           printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** no primary ***:  label %d \n", label);
1324           continue;
1325               }
1326               eprim   = primary->Energy();
1327               ptprim  = primary->Pt();
1328
1329             }
1330             else if(GetReader()->ReadAODMCParticles()){
1331               //Check which is the input
1332               if(ph->GetInputFileIndex() == 0){
1333           if(!mcparticles0) continue;
1334           if(label >=  mcparticles0->GetEntriesFast()) {
1335             if(GetDebug() > 2)  printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** large label ***:  label %d, n tracks %d \n",
1336                                        label, mcparticles0->GetEntriesFast());
1337             continue ;
1338           }
1339           //Get the particle
1340           aodprimary = (AliAODMCParticle*) mcparticles0->At(label);
1341
1342               }
1343 //            else {//Second input
1344 //          if(!mcparticles1) continue;
1345 //          if(label >=  mcparticles1->GetEntriesFast()) {
1346 //            if(GetDebug() > 2)  printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** large label ***:  label %d, n tracks %d \n",
1347 //                                       label, mcparticles1->GetEntriesFast());
1348 //            continue ;
1349 //          }
1350 //          //Get the particle
1351 //          aodprimary = (AliAODMCParticle*) mcparticles1->At(label);
1352 //
1353 //            }//second input
1354
1355               if(!aodprimary){
1356           printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** no primary ***:  label %d \n", label);
1357           continue;
1358               }
1359
1360               eprim   = aodprimary->E();
1361               ptprim  = aodprimary->Pt();
1362
1363             }
1364
1365             fh2E     ->Fill(ecluster, eprim);
1366             fh2Pt    ->Fill(ptcluster, ptprim);
1367             fhDeltaE ->Fill(eprim-ecluster);
1368             fhDeltaPt->Fill(ptprim-ptcluster);
1369             if(eprim > 0)  fhRatioE  ->Fill(ecluster/eprim);
1370             if(ptprim > 0) fhRatioPt ->Fill(ptcluster/ptprim);
1371
1372           }//Histograms with MC
1373
1374         }// aod loop
1375
1376 }
1377
1378
1379 //__________________________________________________________________
1380 void AliAnaPhoton::Print(const Option_t * opt) const
1381 {
1382   //Print some relevant parameters set for the analysis
1383
1384   if(! opt)
1385     return;
1386
1387   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
1388   AliAnaPartCorrBaseClass::Print(" ");
1389
1390   printf("Calorimeter            =     %s\n", fCalorimeter.Data()) ;
1391   printf("Min Distance to Bad Channel   = %2.1f\n",fMinDist);
1392   printf("Min Distance to Bad Channel 2 = %2.1f\n",fMinDist2);
1393   printf("Min Distance to Bad Channel 3 = %2.1f\n",fMinDist3);
1394   printf("Reject clusters with a track matched = %d\n",fRejectTrackMatch);
1395   printf("Check Pair Conversion                = %d\n",fCheckConversion);
1396   printf("Add conversion pair to AOD           = %d\n",fAddConvertedPairsToAOD);
1397   printf("Conversion pair mass cut             = %f\n",fMassCut);
1398   printf("Conversion selection cut : A < %1.2f; %1.3f < Dphi < %1.3f; Deta < %1.3f\n",
1399          fConvAsymCut,fConvDPhiMinCut, fConvDPhiMaxCut, fConvDEtaCut);
1400   printf("Time Cut: %3.1f < TOF  < %3.1f\n", fTimeCutMin, fTimeCutMax);
1401   printf("Number of cells in cluster is        > %d \n", fNCellsCut);
1402   printf("    \n") ;
1403
1404 }