PWGGA/CaloTrackCorrelations/AliAnaPhotonConvInCalo.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
  16 //_________________________________________________________________________
  17 //
  18 // Conversions pairs analysis
  19 // Check if cluster comes from a conversion in the material in front of the calorimeter
  20 // Do invariant mass of all pairs, if mass is close to 0, then it is conversion.
  21 // Input are selected clusters with AliAnaPhoton
  22 //
  23 //
  24 //-- Author: Gustavo Conesa (LPSC-IN2P3-CNRS)
  25 //////////////////////////////////////////////////////////////////////////////
  26
  27
  28 // --- ROOT system ---
  29 #include <TH2F.h>
  30 #include <TH3D.h>
  31 #include <TClonesArray.h>
  32 #include <TObjString.h>
  33 #include "TParticle.h"
  34 #include "TDatabasePDG.h"
  35
  36 // --- Analysis system ---
  37 #include "AliAnaPhotonConvInCalo.h"
  38 #include "AliCaloTrackReader.h"
  39 #include "AliStack.h"
  40 #include "AliCaloPID.h"
  41 #include "AliMCAnalysisUtils.h"
  42 #include "AliFiducialCut.h"
  43 #include "AliVCluster.h"
  44 #include "AliAODMCParticle.h"
  45
  46 ClassImp(AliAnaPhotonConvInCalo)
  47
  48 //________________________________________
  49 AliAnaPhotonConvInCalo::AliAnaPhotonConvInCalo() :
  50 AliAnaCaloTrackCorrBaseClass(),
  51 fRemoveConvertedPair(kFALSE),
  52 fAddConvertedPairsToAOD(kFALSE),
  53 fMassCut(0),
  54 fConvAsymCut(1.),                  fConvDEtaCut(2.),
  55 fConvDPhiMinCut(-1.),              fConvDPhiMaxCut(7.),
  56
  57 // Histograms
  58 fhPtPhotonConv(0),                 fhEtaPhiPhotonConv(0),          fhEtaPhi05PhotonConv(0),
  59 fhConvDeltaEta(0),                 fhConvDeltaPhi(0),              fhConvDeltaEtaPhi(0),
  60 fhConvAsym(0),                     fhConvPt(0),
  61 fhConvDistEta(0),                  fhConvDistEn(0),                fhConvDistMass(0),
  62 fhConvDistEtaCutEta(0),            fhConvDistEnCutEta(0),          fhConvDistMassCutEta(0),
  63 fhConvDistEtaCutMass(0),           fhConvDistEnCutMass(0),
  64 fhConvDistEtaCutAsy(0),            fhConvDistEnCutAsy(0),
  65
  66 // MC histograms
  67 fhPtConversionTagged(0),           fhPtAntiNeutronTagged(0),
  68 fhPtAntiProtonTagged(0),           fhPtUnknownTagged(0),
  69
  70 fhConvDeltaEtaMCConversion(0),     fhConvDeltaPhiMCConversion(0),  fhConvDeltaEtaPhiMCConversion(0),
  71 fhConvAsymMCConversion(0),         fhConvPtMCConversion(0),
  72 fhConvDispersionMCConversion(0),   fhConvM02MCConversion(0),
  73
  74 fhConvDeltaEtaMCAntiNeutron(0),    fhConvDeltaPhiMCAntiNeutron(0), fhConvDeltaEtaPhiMCAntiNeutron(0),
  75 fhConvAsymMCAntiNeutron(0),        fhConvPtMCAntiNeutron(0),
  76 fhConvDispersionMCAntiNeutron(0),  fhConvM02MCAntiNeutron(0),
  77 fhConvDeltaEtaMCAntiProton(0),     fhConvDeltaPhiMCAntiProton(0),  fhConvDeltaEtaPhiMCAntiProton(0),
  78 fhConvAsymMCAntiProton(0),         fhConvPtMCAntiProton(0),
  79 fhConvDispersionMCAntiProton(0),   fhConvM02MCAntiProton(0),
  80 fhConvDeltaEtaMCString(0),         fhConvDeltaPhiMCString(0),      fhConvDeltaEtaPhiMCString(0),
  81 fhConvAsymMCString(0),             fhConvPtMCString(0),
  82 fhConvDispersionMCString(0),       fhConvM02MCString(0),
  83 fhConvDistMCConversion(0),         fhConvDistMCConversionCuts(0)
  84 {
  85   //default ctor
  86
  87   //Initialize parameters
  88   InitParameters();
  89
  90 }
  91
  92 //_________________________________________________
  93 TObjString *  AliAnaPhotonConvInCalo::GetAnalysisCuts()
  94 {
  95   //Save parameters used for analysis
  96   TString parList ; //this will be list of parameters used for this analysis.
  97   const Int_t buffersize = 255;
  98   char onePar[buffersize] ;
  99
 100   snprintf(onePar,buffersize,"--- AliAnaPhotonConvInCalo---\n") ;
 101   parList+=onePar ;
 102   snprintf(onePar,buffersize,"Conversion Selection: fConvAsymCut %1.2f, fConvDEtaCut %1.2f fConvDPhiCut (%1.2f,%1.2f)\n",
 103            fConvAsymCut, fConvDEtaCut, fConvDPhiMinCut, fConvDPhiMaxCut) ;
 104   parList+=onePar ;
 105
 106   return new TObjString(parList) ;
 107 }
 108
 109 //___________________________________________________
 110 TList *  AliAnaPhotonConvInCalo::GetCreateOutputObjects()
 111 {
 112   // Create histograms to be saved in output file and
 113   // store them in outputContainer
 114   TList * outputContainer = new TList() ;
 115   outputContainer->SetName("PhotonConvInCaloHistos") ;
 116
 117   Int_t nptbins  = GetHistogramRanges()->GetHistoPtBins();  Float_t ptmax  = GetHistogramRanges()->GetHistoPtMax();  Float_t ptmin  = GetHistogramRanges()->GetHistoPtMin();
 118   Int_t nphibins = GetHistogramRanges()->GetHistoPhiBins(); Float_t phimax = GetHistogramRanges()->GetHistoPhiMax(); Float_t phimin = GetHistogramRanges()->GetHistoPhiMin();
 119   Int_t netabins = GetHistogramRanges()->GetHistoEtaBins(); Float_t etamax = GetHistogramRanges()->GetHistoEtaMax(); Float_t etamin = GetHistogramRanges()->GetHistoEtaMin();
 120
 121   fhPtPhotonConv  = new TH1F("hPtPhotonConv","Number of #gamma over calorimeter, conversion",nptbins,ptmin,ptmax);
 122   fhPtPhotonConv->SetYTitle("N");
 123   fhPtPhotonConv->SetXTitle("p_{T #gamma}(GeV/c)");
 124   outputContainer->Add(fhPtPhotonConv) ;
 125
 126   fhEtaPhiPhotonConv  = new TH2F
 127   ("hEtaPhiPhotonConv","#eta vs #phi",netabins,etamin,etamax,nphibins,phimin,phimax);
 128   fhEtaPhiPhotonConv->SetYTitle("#phi (rad)");
 129   fhEtaPhiPhotonConv->SetXTitle("#eta");
 130   outputContainer->Add(fhEtaPhiPhotonConv) ;
 131   if(GetMinPt() < 0.5){
 132     fhEtaPhi05PhotonConv  = new TH2F
 133     ("hEtaPhi05PhotonConv","#eta vs #phi, E > 0.5",netabins,etamin,etamax,nphibins,phimin,phimax);
 134     fhEtaPhi05PhotonConv->SetYTitle("#phi (rad)");
 135     fhEtaPhi05PhotonConv->SetXTitle("#eta");
 136     outputContainer->Add(fhEtaPhi05PhotonConv) ;
 137   }
 138
 139   fhConvDeltaEta  = new TH2F
 140   ("hConvDeltaEta","#Delta #eta of selected conversion pairs",100,0,fMassCut,netabins*2,-0.5,0.5);
 141   fhConvDeltaEta->SetYTitle("#Delta #eta");
 142   fhConvDeltaEta->SetXTitle("Pair Mass (GeV/c^2)");
 143   outputContainer->Add(fhConvDeltaEta) ;
 144
 145   fhConvDeltaPhi  = new TH2F
 146   ("hConvDeltaPhi","#Delta #phi of selected conversion pairs",100,0,fMassCut,nphibins*2,-0.5,0.5);
 147   fhConvDeltaPhi->SetYTitle("#Delta #phi");
 148   fhConvDeltaPhi->SetXTitle("Pair Mass (GeV/c^2)");
 149   outputContainer->Add(fhConvDeltaPhi) ;
 150
 151   fhConvDeltaEtaPhi  = new TH2F
 152   ("hConvDeltaEtaPhi","#Delta #eta vs #Delta #phi of selected conversion pairs",netabins,-0.5,0.5,nphibins,-0.5,0.5);
 153   fhConvDeltaEtaPhi->SetYTitle("#Delta #phi");
 154   fhConvDeltaEtaPhi->SetXTitle("#Delta #eta");
 155   outputContainer->Add(fhConvDeltaEtaPhi) ;
 156
 157   fhConvAsym  = new TH2F
 158   ("hConvAsym","Asymmetry of selected conversion pairs",100,0,fMassCut,100,0,1);
 159   fhConvAsym->SetYTitle("Asymmetry");
 160   fhConvAsym->SetXTitle("Pair Mass (GeV/c^2)");
 161   outputContainer->Add(fhConvAsym) ;
 162
 163   fhConvPt  = new TH2F
 164   ("hConvPt","p_{T} of selected conversion pairs",100,0,fMassCut,100,0.,10.);
 165   fhConvPt->SetYTitle("Pair p_{T} (GeV/c)");
 166   fhConvPt->SetXTitle("Pair Mass (GeV/c^2)");
 167   outputContainer->Add(fhConvPt) ;
 168
 169   fhConvDistEta  = new TH2F
 170   ("hConvDistEta","distance to conversion vertex",100,-0.7,0.7,100,0.,5.);
 171   fhConvDistEta->SetXTitle("#eta");
 172   fhConvDistEta->SetYTitle(" distance (m)");
 173   outputContainer->Add(fhConvDistEta) ;
 174
 175   fhConvDistEn  = new TH2F
 176   ("hConvDistEn","distance to conversion vertex",nptbins,ptmin,ptmax,100,0.,5.);
 177   fhConvDistEn->SetXTitle("E (GeV)");
 178   fhConvDistEn->SetYTitle(" distance (m)");
 179   outputContainer->Add(fhConvDistEn) ;
 180
 181   fhConvDistMass  = new TH2F
 182   ("hConvDistMass","distance to conversion vertex",100,0,fMassCut,100,0.,5.);
 183   fhConvDistMass->SetXTitle("m (GeV/c^2)");
 184   fhConvDistMass->SetYTitle(" distance (m)");
 185   outputContainer->Add(fhConvDistMass) ;
 186
 187   fhConvDistEtaCutEta  = new TH2F
 188   ("hConvDistEtaCutEta","distance to conversion vertex, dEta < 0.05",100,-0.7,0.7,100,0.,5.);
 189   fhConvDistEtaCutEta->SetXTitle("#eta");
 190   fhConvDistEtaCutEta->SetYTitle(" distance (m)");
 191   outputContainer->Add(fhConvDistEtaCutEta) ;
 192
 193   fhConvDistEnCutEta  = new TH2F
 194   ("hConvDistEnCutEta","distance to conversion vertex, dEta < 0.05",nptbins,ptmin,ptmax,100,0.,5.);
 195   fhConvDistEnCutEta->SetXTitle("E (GeV)");
 196   fhConvDistEnCutEta->SetYTitle(" distance (m)");
 197   outputContainer->Add(fhConvDistEnCutEta) ;
 198
 199   fhConvDistMassCutEta  = new TH2F
 200   ("hConvDistMassCutEta","distance to conversion vertex, dEta < 0.05",100,0,fMassCut,100,0.,5.);
 201   fhConvDistMassCutEta->SetXTitle("m (GeV/c^2)");
 202   fhConvDistMassCutEta->SetYTitle(" distance (m)");
 203   outputContainer->Add(fhConvDistMassCutEta) ;
 204
 205   fhConvDistEtaCutMass  = new TH2F
 206   ("hConvDistEtaCutMass","distance to conversion vertex, dEta < 0.05, m < 10 MeV",100,-0.7,0.7,100,0.,5.);
 207   fhConvDistEtaCutMass->SetXTitle("#eta");
 208   fhConvDistEtaCutMass->SetYTitle(" distance (m)");
 209   outputContainer->Add(fhConvDistEtaCutMass) ;
 210
 211   fhConvDistEnCutMass  = new TH2F
 212   ("hConvDistEnCutMass","distance to conversion vertex, dEta < 0.05, m < 10 MeV",nptbins,ptmin,ptmax,100,0.,5.);
 213   fhConvDistEnCutMass->SetXTitle("E (GeV)");
 214   fhConvDistEnCutMass->SetYTitle(" distance (m)");
 215   outputContainer->Add(fhConvDistEnCutMass) ;
 216
 217   fhConvDistEtaCutAsy  = new TH2F
 218   ("hConvDistEtaCutAsy","distance to conversion vertex, dEta < 0.05, m < 10 MeV, A < 0.1",100,-0.7,0.7,100,0.,5.);
 219   fhConvDistEtaCutAsy->SetXTitle("#eta");
 220   fhConvDistEtaCutAsy->SetYTitle(" distance (m)");
 221   outputContainer->Add(fhConvDistEtaCutAsy) ;
 222
 223   fhConvDistEnCutAsy  = new TH2F
 224   ("hConvDistEnCutAsy","distance to conversion vertex, dEta < 0.05, m < 10 MeV, A < 0.1",nptbins,ptmin,ptmax,100,0.,5.);
 225   fhConvDistEnCutAsy->SetXTitle("E (GeV)");
 226   fhConvDistEnCutAsy->SetYTitle(" distance (m)");
 227   outputContainer->Add(fhConvDistEnCutAsy) ;
 228
 229   if(IsDataMC()){
 230
 231     fhPtConversionTagged  = new TH1F("hPtMCConversionTagged","Number of converted #gamma over calorimeter, tagged as converted",nptbins,ptmin,ptmax);
 232     fhPtConversionTagged->SetYTitle("N");
 233     fhPtConversionTagged->SetXTitle("p_{T #gamma}(GeV/c)");
 234     outputContainer->Add(fhPtConversionTagged) ;
 235
 236
 237     fhPtAntiNeutronTagged  = new TH1F("hPtMCAntiNeutronTagged","Number of AntiNeutron id as Photon over calorimeter, tagged as converted",nptbins,ptmin,ptmax);
 238     fhPtAntiNeutronTagged->SetYTitle("N");
 239     fhPtAntiNeutronTagged->SetXTitle("p_{T #gamma}(GeV/c)");
 240     outputContainer->Add(fhPtAntiNeutronTagged) ;
 241
 242     fhPtAntiProtonTagged  = new TH1F("hPtMCAntiProtonTagged","Number of AntiProton id as Photon over calorimeter, tagged as converted",nptbins,ptmin,ptmax);
 243     fhPtAntiProtonTagged->SetYTitle("N");
 244     fhPtAntiProtonTagged->SetXTitle("p_{T #gamma}(GeV/c)");
 245     outputContainer->Add(fhPtAntiProtonTagged) ;
 246
 247     fhPtUnknownTagged  = new TH1F("hPtMCUnknownTagged","Number of Unknown id as Photon over calorimeter, tagged as converted",nptbins,ptmin,ptmax);
 248     fhPtUnknownTagged->SetYTitle("N");
 249     fhPtUnknownTagged->SetXTitle("p_{T #gamma}(GeV/c)");
 250     outputContainer->Add(fhPtUnknownTagged) ;
 251
 252     fhConvDeltaEtaMCConversion  = new TH2F
 253     ("hConvDeltaEtaMCConversion","#Delta #eta of selected conversion pairs from real conversions",100,0,fMassCut,netabins,-0.5,0.5);
 254     fhConvDeltaEtaMCConversion->SetYTitle("#Delta #eta");
 255     fhConvDeltaEtaMCConversion->SetXTitle("Pair Mass (GeV/c^2)");
 256     outputContainer->Add(fhConvDeltaEtaMCConversion) ;
 257
 258     fhConvDeltaPhiMCConversion  = new TH2F
 259     ("hConvDeltaPhiMCConversion","#Delta #phi of selected conversion pairs from real conversions",100,0,fMassCut,nphibins,-0.5,0.5);
 260     fhConvDeltaPhiMCConversion->SetYTitle("#Delta #phi");
 261     fhConvDeltaPhiMCConversion->SetXTitle("Pair Mass (GeV/c^2)");
 262     outputContainer->Add(fhConvDeltaPhiMCConversion) ;
 263
 264     fhConvDeltaEtaPhiMCConversion  = new TH2F
 265     ("hConvDeltaEtaPhiMCConversion","#Delta #eta vs #Delta #phi of selected conversion pairs, from real conversions",netabins,-0.5,0.5,nphibins,-0.5,0.5);
 266     fhConvDeltaEtaPhiMCConversion->SetYTitle("#Delta #phi");
 267     fhConvDeltaEtaPhiMCConversion->SetXTitle("#Delta #eta");
 268     outputContainer->Add(fhConvDeltaEtaPhiMCConversion) ;
 269
 270     fhConvAsymMCConversion  = new TH2F
 271     ("hConvAsymMCConversion","Asymmetry of selected conversion pairs from real conversions",100,0,fMassCut,100,0,1);
 272     fhConvAsymMCConversion->SetYTitle("Asymmetry");
 273     fhConvAsymMCConversion->SetXTitle("Pair Mass (GeV/c^2)");
 274     outputContainer->Add(fhConvAsymMCConversion) ;
 275
 276     fhConvPtMCConversion  = new TH2F
 277     ("hConvPtMCConversion","p_{T} of selected conversion pairs from real conversions",100,0,fMassCut,100,0.,10.);
 278     fhConvPtMCConversion->SetYTitle("Pair p_{T} (GeV/c)");
 279     fhConvPtMCConversion->SetXTitle("Pair Mass (GeV/c^2)");
 280     outputContainer->Add(fhConvPtMCConversion) ;
 281
 282     fhConvDispersionMCConversion  = new TH2F
 283     ("hConvDispersionMCConversion","p_{T} of selected conversion pairs from real conversions",100,0.,1.,100,0.,1.);
 284     fhConvDispersionMCConversion->SetYTitle("Dispersion cluster 1");
 285     fhConvDispersionMCConversion->SetXTitle("Dispersion cluster 2");
 286     outputContainer->Add(fhConvDispersionMCConversion) ;
 287
 288     fhConvM02MCConversion  = new TH2F
 289     ("hConvM02MCConversion","p_{T} of selected conversion pairs from string",100,0.,1.,100,0.,1.);
 290     fhConvM02MCConversion->SetYTitle("M02 cluster 1");
 291     fhConvM02MCConversion->SetXTitle("M02 cluster 2");
 292     outputContainer->Add(fhConvM02MCConversion) ;
 293
 294     fhConvDeltaEtaMCAntiNeutron  = new TH2F
 295     ("hConvDeltaEtaMCAntiNeutron","#Delta #eta of selected conversion pairs from anti-neutrons",100,0,fMassCut,netabins,-0.5,0.5);
 296     fhConvDeltaEtaMCAntiNeutron->SetYTitle("#Delta #eta");
 297     fhConvDeltaEtaMCAntiNeutron->SetXTitle("Pair Mass (GeV/c^2)");
 298     outputContainer->Add(fhConvDeltaEtaMCAntiNeutron) ;
 299
 300     fhConvDeltaPhiMCAntiNeutron  = new TH2F
 301     ("hConvDeltaPhiMCAntiNeutron","#Delta #phi of selected conversion pairs from anti-neutrons",100,0,fMassCut,nphibins,-0.5,0.5);
 302     fhConvDeltaPhiMCAntiNeutron->SetYTitle("#Delta #phi");
 303     fhConvDeltaPhiMCAntiNeutron->SetXTitle("Pair Mass (GeV/c^2)");
 304     outputContainer->Add(fhConvDeltaPhiMCAntiNeutron) ;
 305
 306     fhConvDeltaEtaPhiMCAntiNeutron  = new TH2F
 307     ("hConvDeltaEtaPhiMCAntiNeutron","#Delta #eta vs #Delta #phi of selected conversion pairs from anti-neutrons",netabins,-0.5,0.5,nphibins,-0.5,0.5);
 308     fhConvDeltaEtaPhiMCAntiNeutron->SetYTitle("#Delta #phi");
 309     fhConvDeltaEtaPhiMCAntiNeutron->SetXTitle("#Delta #eta");
 310     outputContainer->Add(fhConvDeltaEtaPhiMCAntiNeutron) ;
 311
 312     fhConvAsymMCAntiNeutron  = new TH2F
 313     ("hConvAsymMCAntiNeutron","Asymmetry of selected conversion pairs from anti-neutrons",100,0,fMassCut,100,0,1);
 314     fhConvAsymMCAntiNeutron->SetYTitle("Asymmetry");
 315     fhConvAsymMCAntiNeutron->SetXTitle("Pair Mass (GeV/c^2)");
 316     outputContainer->Add(fhConvAsymMCAntiNeutron) ;
 317
 318     fhConvPtMCAntiNeutron  = new TH2F
 319     ("hConvPtMCAntiNeutron","p_{T} of selected conversion pairs from anti-neutrons",100,0,fMassCut,100,0.,10.);
 320     fhConvPtMCAntiNeutron->SetYTitle("Pair p_{T} (GeV/c)");
 321     fhConvPtMCAntiNeutron->SetXTitle("Pair Mass (GeV/c^2)");
 322     outputContainer->Add(fhConvPtMCAntiNeutron) ;
 323
 324     fhConvDispersionMCAntiNeutron  = new TH2F
 325     ("hConvDispersionMCAntiNeutron","p_{T} of selected conversion pairs from anti-neutrons",100,0.,1.,100,0.,1.);
 326     fhConvDispersionMCAntiNeutron->SetYTitle("Dispersion cluster 1");
 327     fhConvDispersionMCAntiNeutron->SetXTitle("Dispersion cluster 2");
 328     outputContainer->Add(fhConvDispersionMCAntiNeutron) ;
 329
 330     fhConvM02MCAntiNeutron  = new TH2F
 331     ("hConvM02MCAntiNeutron","p_{T} of selected conversion pairs from string",100,0.,1.,100,0.,1.);
 332     fhConvM02MCAntiNeutron->SetYTitle("M02 cluster 1");
 333     fhConvM02MCAntiNeutron->SetXTitle("M02 cluster 2");
 334     outputContainer->Add(fhConvM02MCAntiNeutron) ;
 335
 336     fhConvDeltaEtaMCAntiProton  = new TH2F
 337     ("hConvDeltaEtaMCAntiProton","#Delta #eta of selected conversion pairs from anti-protons",100,0,fMassCut,netabins,-0.5,0.5);
 338     fhConvDeltaEtaMCAntiProton->SetYTitle("#Delta #eta");
 339     fhConvDeltaEtaMCAntiProton->SetXTitle("Pair Mass (GeV/c^2)");
 340     outputContainer->Add(fhConvDeltaEtaMCAntiProton) ;
 341
 342     fhConvDeltaPhiMCAntiProton  = new TH2F
 343     ("hConvDeltaPhiMCAntiProton","#Delta #phi of selected conversion pairs from anti-protons",100,0,fMassCut,nphibins,-0.5,0.5);
 344     fhConvDeltaPhiMCAntiProton->SetYTitle("#Delta #phi");
 345     fhConvDeltaPhiMCAntiProton->SetXTitle("Pair Mass (GeV/c^2)");
 346     outputContainer->Add(fhConvDeltaPhiMCAntiProton) ;
 347
 348     fhConvDeltaEtaPhiMCAntiProton  = new TH2F
 349     ("hConvDeltaEtaPhiMCAntiProton","#Delta #eta vs #Delta #phi of selected conversion pairs from anti-protons",netabins,-0.5,0.5,nphibins,-0.5,0.5);
 350     fhConvDeltaEtaPhiMCAntiProton->SetYTitle("#Delta #phi");
 351     fhConvDeltaEtaPhiMCAntiProton->SetXTitle("#Delta #eta");
 352     outputContainer->Add(fhConvDeltaEtaPhiMCAntiProton) ;
 353
 354     fhConvAsymMCAntiProton  = new TH2F
 355     ("hConvAsymMCAntiProton","Asymmetry of selected conversion pairs from anti-protons",100,0,fMassCut,100,0,1);
 356     fhConvAsymMCAntiProton->SetYTitle("Asymmetry");
 357     fhConvAsymMCAntiProton->SetXTitle("Pair Mass (GeV/c^2)");
 358     outputContainer->Add(fhConvAsymMCAntiProton) ;
 359
 360     fhConvPtMCAntiProton  = new TH2F
 361     ("hConvPtMCAntiProton","p_{T} of selected conversion pairs from anti-protons",100,0,fMassCut,100,0.,10.);
 362     fhConvPtMCAntiProton->SetYTitle("Pair p_{T} (GeV/c)");
 363     fhConvPtMCAntiProton->SetXTitle("Pair Mass (GeV/c^2)");
 364     outputContainer->Add(fhConvPtMCAntiProton) ;
 365
 366     fhConvDispersionMCAntiProton  = new TH2F
 367     ("hConvDispersionMCAntiProton","p_{T} of selected conversion pairs from anti-protons",100,0.,1.,100,0.,1.);
 368     fhConvDispersionMCAntiProton->SetYTitle("Dispersion cluster 1");
 369     fhConvDispersionMCAntiProton->SetXTitle("Dispersion cluster 2");
 370     outputContainer->Add(fhConvDispersionMCAntiProton) ;
 371
 372     fhConvM02MCAntiProton  = new TH2F
 373     ("hConvM02MCAntiProton","p_{T} of selected conversion pairs from string",100,0.,1.,100,0.,1.);
 374     fhConvM02MCAntiProton->SetYTitle("M02 cluster 1");
 375     fhConvM02MCAntiProton->SetXTitle("M02 cluster 2");
 376     outputContainer->Add(fhConvM02MCAntiProton) ;
 377
 378     fhConvDeltaEtaMCString  = new TH2F
 379     ("hConvDeltaEtaMCString","#Delta #eta of selected conversion pairs from string",100,0,fMassCut,netabins,-0.5,0.5);
 380     fhConvDeltaEtaMCString->SetYTitle("#Delta #eta");
 381     fhConvDeltaEtaMCString->SetXTitle("Pair Mass (GeV/c^2)");
 382     outputContainer->Add(fhConvDeltaEtaMCString) ;
 383
 384     fhConvDeltaPhiMCString  = new TH2F
 385     ("hConvDeltaPhiMCString","#Delta #phi of selected conversion pairs from string",100,0,fMassCut,nphibins,-0.5,0.5);
 386     fhConvDeltaPhiMCString->SetYTitle("#Delta #phi");
 387     fhConvDeltaPhiMCString->SetXTitle("Pair Mass (GeV/c^2)");
 388     outputContainer->Add(fhConvDeltaPhiMCString) ;
 389
 390     fhConvDeltaEtaPhiMCString  = new TH2F
 391     ("hConvDeltaEtaPhiMCString","#Delta #eta vs #Delta #phi of selected conversion pairs from string",netabins,-0.5,0.5,nphibins,-0.5,0.5);
 392     fhConvDeltaEtaPhiMCString->SetYTitle("#Delta #phi");
 393     fhConvDeltaEtaPhiMCString->SetXTitle("#Delta #eta");
 394     outputContainer->Add(fhConvDeltaEtaPhiMCString) ;
 395
 396     fhConvAsymMCString  = new TH2F
 397     ("hConvAsymMCString","Asymmetry of selected conversion pairs from string",100,0,fMassCut,100,0,1);
 398     fhConvAsymMCString->SetYTitle("Asymmetry");
 399     fhConvAsymMCString->SetXTitle("Pair Mass (GeV/c^2)");
 400     outputContainer->Add(fhConvAsymMCString) ;
 401
 402     fhConvPtMCString  = new TH2F
 403     ("hConvPtMCString","p_{T} of selected conversion pairs from string",100,0,fMassCut,100,0.,10.);
 404     fhConvPtMCString->SetYTitle("Pair p_{T} (GeV/c)");
 405     fhConvPtMCString->SetXTitle("Pair Mass (GeV/c^2)");
 406     outputContainer->Add(fhConvPtMCString) ;
 407
 408     fhConvDispersionMCString  = new TH2F
 409     ("hConvDispersionMCString","p_{T} of selected conversion pairs from string",100,0.,1.,100,0.,1.);
 410     fhConvDispersionMCString->SetYTitle("Dispersion cluster 1");
 411     fhConvDispersionMCString->SetXTitle("Dispersion cluster 2");
 412     outputContainer->Add(fhConvDispersionMCString) ;
 413
 414     fhConvM02MCString  = new TH2F
 415     ("hConvM02MCString","p_{T} of selected conversion pairs from string",100,0.,1.,100,0.,1.);
 416     fhConvM02MCString->SetYTitle("M02 cluster 1");
 417     fhConvM02MCString->SetXTitle("M02 cluster 2");
 418     outputContainer->Add(fhConvM02MCString) ;
 419
 420     fhConvDistMCConversion  = new TH2F
 421     ("hConvDistMCConversion","calculated conversion distance vs real vertes for MC conversion",100,0.,5.,100,0.,5.);
 422     fhConvDistMCConversion->SetYTitle("distance");
 423     fhConvDistMCConversion->SetXTitle("vertex R");
 424     outputContainer->Add(fhConvDistMCConversion) ;
 425
 426     fhConvDistMCConversionCuts  = new TH2F
 427     ("hConvDistMCConversionCuts","calculated conversion distance vs real vertes for MC conversion, deta < 0.05, m < 10 MeV, asym < 0.1",100,0.,5.,100,0.,5.);
 428     fhConvDistMCConversionCuts->SetYTitle("distance");
 429     fhConvDistMCConversionCuts->SetXTitle("vertex R");
 430     outputContainer->Add(fhConvDistMCConversionCuts) ;
 431
 432   }
 433
 434   return outputContainer ;
 435
 436 }
 437
 438 //_______________________________________
 439 void AliAnaPhotonConvInCalo::InitParameters()
 440 {
 441
 442   //Initialize the parameters of the analysis.
 443   AddToHistogramsName("AnaPhotonConvInCalo_");
 444
 445   fMassCut                = 0.03; //30 MeV
 446   fRemoveConvertedPair    = kFALSE;
 447   fAddConvertedPairsToAOD = kFALSE;
 448
 449 }
 450
 451 //_____________________________________________
 452 void  AliAnaPhotonConvInCalo::MakeAnalysisFillAOD()
 453 {
 454   //Do conversion photon analysis and fill aods
 455
 456   //Loop on stored AOD photons
 457   Int_t naod = GetOutputAODBranch()->GetEntriesFast();
 458   if(GetDebug() > 0) printf("AliAnaPhotonConvInCalo::MakeAnalysisFillAOD() - aod branch entries %d\n", naod);
 459
 460   //List to be used in conversion analysis, to tag the cluster as candidate for conversion
 461   Bool_t * indexConverted = new Bool_t[naod];
 462   for (Int_t i = 0; i < naod; i++) indexConverted[i] = kFALSE;
 463
 464   for(Int_t iaod = 0; iaod < naod ; iaod++){
 465     AliAODPWG4Particle* calo =  (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
 466
 467     Bool_t bConverted = kFALSE;
 468     Int_t id2 = -1;
 469
 470     //Check if set previously as converted couple, if so skip its use.
 471     if (indexConverted[iaod]) continue;
 472
 473     // Second cluster loop
 474     AliAODPWG4Particle* calo2 = 0;
 475     for(Int_t jaod = iaod + 1 ; jaod < naod ; jaod++) {
 476       //Check if set previously as converted couple, if so skip its use.
 477       if (indexConverted[jaod]) continue;
 478       //printf("Check Conversion indeces %d and %d\n",iaod,jaod);
 479       calo2 =  (AliAODPWG4Particle*) (GetOutputAODBranch()->At(jaod));
 480
 481       //................................................
 482       //Get mass of pair, if small, take this pair.
 483       Float_t pairM     = calo->GetPairMass(calo2);
 484       //printf("\t both in calo, mass %f, cut %f\n",pairM,fMassCut);
 485       if(pairM < fMassCut){
 486         calo->SetTagged(kFALSE);
 487         id2 = calo2->GetCaloLabel(0);
 488         Float_t asymmetry = TMath::Abs(calo->E()-calo2->E())/(calo->E()+calo2->E());
 489         Float_t dPhi      = (calo->Momentum())->Phi()-(calo2->Momentum())->Phi();
 490         Float_t dEta      = (calo->Momentum())->Eta()-(calo2->Momentum())->Eta();
 491
 492         //...............................................
 493         //Fill few histograms with kinematics of the pair
 494         //FIXME, move all this to MakeAnalysisFillHistograms ...
 495
 496         fhConvDeltaEta   ->Fill( pairM, dPhi      );
 497         fhConvDeltaPhi   ->Fill( pairM, dEta      );
 498         fhConvAsym       ->Fill( pairM, asymmetry );
 499         fhConvDeltaEtaPhi->Fill( dEta , dPhi      );
 500         fhConvPt         ->Fill( pairM, (calo->Momentum())->Pt()+(calo2->Momentum())->Pt());
 501
 502         //Estimate conversion distance, T. Awes, M. Ivanov
 503         //Under the assumption that the pair has zero mass, and that each electron
 504         //of the pair has the same momentum, they will each have the same bend radius
 505         //given by R=p/(qB) = p / (300 B) with p in [MeV/c], B in [Tesla] and R in [m].
 506         //With nominal ALICE magnet current of 30kA B=0.5T, and so with E_cluster=p,
 507         //R = E/1.5 [cm].  Under these assumptions, the distance from the conversion
 508         //point to the MCEal can be related to the separation distance, L=2y, on the MCEal
 509         //as d = sqrt(R^2 -(R-y)^2) = sqrt(2Ry - y^2). And since R>>y we can write as
 510         //d = sqrt(E*L/1.5) where E is the cluster energy and L is the distance in cm between
 511         //the clusters.
 512
 513         TObjArray * clusters    = 0;
 514         if(calo->GetDetector() == "EMCAL"){
 515           clusters = GetEMCALClusters();
 516         }
 517         else{
 518           clusters = GetPHOSClusters();
 519         }
 520
 521         Int_t iclus = -1;
 522         AliVCluster *cluster1 = FindCluster(clusters,calo ->GetCaloLabel(0),iclus);
 523         AliVCluster *cluster2 = FindCluster(clusters,calo2->GetCaloLabel(0),iclus);
 524
 525         Float_t pos1[3];
 526         cluster1->GetPosition(pos1);
 527         Float_t pos2[3];
 528         cluster2->GetPosition(pos2);
 529         Float_t clustDist = TMath::Sqrt((pos1[0]-pos2[0])*(pos1[0]-pos2[0])+
 530                                         (pos1[1]-pos2[1])*(pos1[1]-pos2[1])+
 531                                         (pos1[2]-pos2[2])*(pos1[2]-pos2[2]));
 532
 533         Float_t convDist  = TMath::Sqrt(calo->E() *clustDist*0.01/0.15);
 534         Float_t convDist2 = TMath::Sqrt(calo2->E()*clustDist*0.01/0.15);
 535         //printf("l = %f, e1 = %f, d1=%f, e2 = %f, d2=%f\n",clustDist,calo->E(),convDist,calo2->E(),convDist2);
 536         if(GetDebug() > 2)
 537           printf("AliAnaPhotonConvInCalo::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",
 538                  pairM,fMassCut,fConvDPhiMinCut, dPhi, fConvDPhiMaxCut, dEta, fConvDEtaCut, asymmetry, fConvAsymCut,
 539                  calo->GetCaloLabel(0),calo->E(),GetCaloUtils()->GetModuleNumber(calo,GetReader()->GetInputEvent()), calo->Eta(), calo->Phi(),
 540                  id2, calo2->E(), GetCaloUtils()->GetModuleNumber(calo2,GetReader()->GetInputEvent()),calo2->Eta(), calo2->Phi());
 541
 542         fhConvDistEta ->Fill(calo ->Eta(),convDist );
 543         fhConvDistEta ->Fill(calo2->Eta(),convDist2);
 544         fhConvDistEn  ->Fill(calo ->E(),  convDist );
 545         fhConvDistEn  ->Fill(calo2->E(),  convDist2);
 546         fhConvDistMass->Fill(pairM, convDist );
 547         //dEta cut
 548         if(dEta<0.05){
 549           fhConvDistEtaCutEta ->Fill(calo->Eta(), convDist );
 550           fhConvDistEtaCutEta ->Fill(calo2->Eta(),convDist2);
 551           fhConvDistEnCutEta  ->Fill(calo->E(),   convDist );
 552           fhConvDistEnCutEta  ->Fill(calo2->E(),  convDist2);
 553           fhConvDistMassCutEta->Fill(pairM, convDist );
 554           //mass cut
 555           if(pairM<0.01){//10 MeV
 556             fhConvDistEtaCutMass ->Fill(calo ->Eta(), convDist );
 557             fhConvDistEtaCutMass ->Fill(calo2->Eta(), convDist2);
 558             fhConvDistEnCutMass  ->Fill(calo ->E(),   convDist );
 559             fhConvDistEnCutMass  ->Fill(calo2->E(),   convDist2);
 560             // asymmetry cut
 561             if(asymmetry<0.1){
 562               fhConvDistEtaCutAsy ->Fill(calo ->Eta(), convDist );
 563               fhConvDistEtaCutAsy ->Fill(calo2->Eta(), convDist2);
 564               fhConvDistEnCutAsy  ->Fill(calo ->E(),   convDist );
 565               fhConvDistEnCutAsy  ->Fill(calo2->E(),   convDist2);
 566             }//asymmetry cut
 567           }//mass cut
 568         }//dEta cut
 569
 570         //...............................................
 571         //Select pairs in a eta-phi window
 572         if(TMath::Abs(dEta) < fConvDEtaCut    &&
 573            TMath::Abs(dPhi) < fConvDPhiMaxCut &&
 574            TMath::Abs(dPhi) > fConvDPhiMinCut &&
 575            asymmetry        < fConvAsymCut       ){
 576           indexConverted[iaod] = kTRUE;
 577           indexConverted[jaod] = kTRUE;
 578           bConverted           = kTRUE;
 579         }
 580         //printf("Accepted? %d\n",bConverted);
 581         //...........................................
 582         //Fill more histograms, simulated data
 583         //FIXME, move all this to MakeAnalysisFillHistograms ...
 584         if(IsDataMC()){
 585
 586           //Check the origin of the pair, look for conversion, antinucleons or jet correlations (strings)
 587           Int_t ancPDG    = 0;
 588           Int_t ancStatus = 0;
 589           TLorentzVector momentum;
 590           TVector3 prodVertex;
 591           Int_t ancLabel  = GetMCAnalysisUtils()->CheckCommonAncestor(cluster1->GetLabel(), cluster2->GetLabel(),
 592                                                                       GetReader(), ancPDG, ancStatus, momentum, prodVertex);
 593
 594           // printf("AliAnaPhotonConvInCalo::MakeAnalysisFillHistograms() - Common ancestor label %d, pdg %d, name %s, status %d; \n",
 595           //                          ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus);
 596
 597           Int_t tag1 = calo ->GetTag();
 598           Int_t tag2 = calo2->GetTag();
 599           if(GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCConversion)){
 600             if(GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCConversion) && (ancPDG==22 || TMath::Abs(ancPDG)==11) && ancLabel > -1){
 601               fhConvDeltaEtaMCConversion   ->Fill( pairM, dEta      );
 602               fhConvDeltaPhiMCConversion   ->Fill( pairM, dPhi      );
 603               fhConvAsymMCConversion       ->Fill( pairM, asymmetry );
 604               fhConvDeltaEtaPhiMCConversion->Fill( dEta , dPhi      );
 605               fhConvPtMCConversion         ->Fill( pairM, calo->Pt()+calo2->Pt());
 606               fhConvDispersionMCConversion ->Fill( cluster1->GetDispersion(), cluster2->GetDispersion());
 607               fhConvM02MCConversion        ->Fill( cluster1->GetM02(), cluster2->GetM02());
 608               fhConvDistMCConversion       ->Fill( convDist , prodVertex.Mag() );
 609               fhConvDistMCConversion       ->Fill( convDist2, prodVertex.Mag() );
 610
 611               if(dEta<0.05 && pairM<0.01 && asymmetry<0.1){
 612                 fhConvDistMCConversionCuts->Fill( convDist , prodVertex.Mag() );
 613                 fhConvDistMCConversionCuts->Fill( convDist2, prodVertex.Mag() );
 614               }
 615
 616             }
 617           }
 618           else if(GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCAntiNeutron)){
 619             if(GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCAntiNeutron) && ancPDG==-2112 && ancLabel > -1){
 620               fhConvDeltaEtaMCAntiNeutron    ->Fill( pairM, dEta      );
 621               fhConvDeltaPhiMCAntiNeutron    ->Fill( pairM, dPhi      );
 622               fhConvAsymMCAntiNeutron        ->Fill( pairM, asymmetry );
 623               fhConvDeltaEtaPhiMCAntiNeutron ->Fill( dEta , dPhi      );
 624               fhConvPtMCAntiNeutron          ->Fill( pairM, calo->Pt()+calo2->Pt());
 625               fhConvDispersionMCAntiNeutron  ->Fill( cluster1->GetDispersion(), cluster2->GetDispersion());
 626               fhConvM02MCAntiNeutron         ->Fill( cluster1->GetM02(), cluster2->GetM02());
 627             }
 628           }
 629           else if(GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCAntiProton)){
 630             if(GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCAntiProton) && ancPDG==-2212 && ancLabel > -1){
 631               fhConvDeltaEtaMCAntiProton    ->Fill( pairM, dEta      );
 632               fhConvDeltaPhiMCAntiProton    ->Fill( pairM, dPhi      );
 633               fhConvAsymMCAntiProton        ->Fill( pairM, asymmetry );
 634               fhConvDeltaEtaPhiMCAntiProton ->Fill( dEta , dPhi      );
 635               fhConvPtMCAntiProton          ->Fill( pairM, calo->Pt()+calo2->Pt());
 636               fhConvDispersionMCAntiProton  ->Fill( cluster1->GetDispersion(), cluster2->GetDispersion());
 637               fhConvM02MCAntiProton         ->Fill( cluster1->GetM02(), cluster2->GetM02());
 638             }
 639           }
 640
 641           //Pairs coming from fragmenting pairs.
 642           if(ancPDG < 22 && ancLabel > 7 && (ancStatus == 11 || ancStatus == 12) ){
 643             fhConvDeltaEtaMCString    ->Fill( pairM, dPhi);
 644             fhConvDeltaPhiMCString    ->Fill( pairM, dPhi);
 645             fhConvAsymMCString        ->Fill( pairM, TMath::Abs(calo->E()-calo2->E())/(calo->E()+calo2->E()) );
 646             fhConvDeltaEtaPhiMCString ->Fill( dPhi,  dPhi);
 647             fhConvPtMCString          ->Fill( pairM, calo->Pt()+calo2->Pt());
 648             fhConvDispersionMCString  ->Fill( cluster1->GetDispersion(), cluster2->GetDispersion());
 649             fhConvM02MCString         ->Fill( cluster1->GetM02(), cluster2->GetM02());
 650           }
 651
 652         }// Data MC
 653
 654         break;
 655       }
 656
 657     }//Mass loop
 658
 659     //..........................................................................................................
 660     //Pair selected as converted, remove both clusters or recombine them into a photon and put them in the AOD
 661     if(bConverted){
 662       //Add to AOD
 663       if(fAddConvertedPairsToAOD){
 664         //Create AOD of pair analysis
 665         TLorentzVector mpair = *(calo->Momentum())+*(calo2->Momentum());
 666         AliAODPWG4Particle aodpair = AliAODPWG4Particle(mpair);
 667         aodpair.SetLabel(calo->GetLabel());
 668
 669         //printf("Index %d, Id %d\n",iaod, calo->GetID());
 670         //Set the indeces of the original caloclusters
 671         aodpair.SetCaloLabel(calo->GetCaloLabel(0),id2);
 672         aodpair.SetDetector(calo->GetDetector());
 673         aodpair.SetIdentifiedParticleType(calo->GetIdentifiedParticleType());
 674         aodpair.SetTag(calo ->GetTag());
 675         aodpair.SetTagged(kTRUE);
 676         //Add AOD with pair object to aod branch
 677         AddAODParticle(aodpair);
 678         //printf("\t \t both added pair\n");
 679       }
 680
 681       //Do not add the current calocluster
 682       if(!fRemoveConvertedPair)
 683       {
 684         //printf("TAGGED\n");
 685         //Tag this cluster as likely conversion
 686         calo->SetTagged(kTRUE);
 687       }
 688     }//converted pair
 689
 690   }// main loop
 691
 692   // Remove entries identified as conversion electrons
 693   // Revise if this is OK
 694   if(fRemoveConvertedPair || fAddConvertedPairsToAOD){
 695     for(Int_t iaod = 0; iaod < naod ; iaod++)
 696       if(indexConverted[iaod])GetOutputAODBranch()->RemoveAt(iaod);
 697     GetOutputAODBranch()->Compress();
 698   }
 699
 700   delete [] indexConverted;
 701
 702   if(GetDebug() > 1) printf("AliAnaPhotonConvInCalo::MakeAnalysisFillAOD()  End fill AODs, with %d entries \n",GetOutputAODBranch()->GetEntriesFast());
 703
 704 }
 705
 706 //____________________________________________________
 707 void  AliAnaPhotonConvInCalo::MakeAnalysisFillHistograms()
 708 {
 709   //Fill histograms
 710
 711   //-------------------------------------------------------------------
 712   // Access MC information in stack if requested, check that it exists.
 713   AliStack         * stack       = 0x0;
 714   TParticle        * primary     = 0x0;
 715   TClonesArray     * mcparticles = 0x0;
 716   AliAODMCParticle * aodprimary  = 0x0;
 717
 718   if(IsDataMC()){
 719
 720     if(GetReader()->ReadStack()){
 721       stack =  GetMCStack() ;
 722       if(!stack) {
 723         printf("AliAnaPhotonConvInCalo::MakeAnalysisFillHistograms() - Stack not available, is the MC handler called? STOP\n");
 724         abort();
 725       }
 726
 727     }
 728     else if(GetReader()->ReadAODMCParticles()){
 729
 730       //Get the list of MC particles
 731       mcparticles = GetReader()->GetAODMCParticles();
 732       if(!mcparticles && GetDebug() > 0)
 733       {
 734         printf("AliAnaPhotonConvInCalo::MakeAnalysisFillHistograms() -  Standard MCParticles not available!\n");
 735       }
 736     }
 737   }// is data and MC
 738
 739   //----------------------------------
 740   //Loop on stored AOD photons
 741   Int_t naod = GetOutputAODBranch()->GetEntriesFast();
 742   if(GetDebug() > 0) printf("AliAnaPhotonConvInCalo::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod);
 743
 744   for(Int_t iaod = 0; iaod < naod ; iaod++){
 745     AliAODPWG4Particle* ph =  (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
 746     //Int_t pdg = ph->GetIdentifiedParticleType();
 747
 748     if(ph->IsTagged()){
 749
 750       if(GetDebug() > 2)
 751         printf("AliAnaPhotonConvInCalo::MakeAnalysisFillHistograms() - ID Photon: pt %f, phi %f, eta %f\n", ph->Pt(),ph->Phi(),ph->Eta()) ;
 752       //................................
 753       //Fill photon histograms
 754       Float_t ptcluster  = ph->Pt();
 755       Float_t phicluster = ph->Phi();
 756       Float_t etacluster = ph->Eta();
 757       Float_t ecluster   = ph->E();
 758
 759       fhPtPhotonConv->Fill(ptcluster);
 760       if(ecluster > 0.5)        fhEtaPhiPhotonConv  ->Fill(etacluster, phicluster);
 761       else if(GetMinPt() < 0.5) fhEtaPhi05PhotonConv->Fill(etacluster, phicluster);
 762
 763
 764       //.......................................
 765       //Play with the MC data if available
 766       if(IsDataMC()){
 767
 768
 769         //....................................................................
 770         // Access MC information in stack if requested, check that it exists.
 771         Int_t label =ph->GetLabel();
 772         if(label < 0) {
 773           if(GetDebug() > 1) printf("AliAnaPhotonConvInCalo::MakeAnalysisFillHistograms() *** bad label ***:  label %d \n", label);
 774           continue;
 775         }
 776
 777         //Float_t eprim   = 0;
 778         //Float_t ptprim  = 0;
 779         if(GetReader()->ReadStack()){
 780
 781           if(label >=  stack->GetNtrack()) {
 782             if(GetDebug() > 2)  printf("AliAnaPhotonConvInCalo::MakeAnalysisFillHistograms() *** large label ***:  label %d, n tracks %d \n", label, stack->GetNtrack());
 783             continue ;
 784           }
 785
 786           primary = stack->Particle(label);
 787           if(!primary){
 788             printf("AliAnaPhotonConvInCalo::MakeAnalysisFillHistograms() *** no primary ***:  label %d \n", label);
 789             continue;
 790           }
 791           //eprim   = primary->Energy();
 792           //ptprim  = primary->Pt();
 793
 794         }
 795         else if(GetReader()->ReadAODMCParticles()){
 796           //Check which is the input
 797           if(ph->GetInputFileIndex() == 0){
 798             if(!mcparticles) continue;
 799             if(label >=  mcparticles->GetEntriesFast()) {
 800               if(GetDebug() > 2)  printf("AliAnaPhotonConvInCalo::MakeAnalysisFillHistograms() *** large label ***:  label %d, n tracks %d \n",
 801                                          label, mcparticles->GetEntriesFast());
 802               continue ;
 803             }
 804             //Get the particle
 805             aodprimary = (AliAODMCParticle*) mcparticles->At(label);
 806
 807           }
 808
 809           if(!aodprimary){
 810             printf("AliAnaPhotonConvInCalo::MakeAnalysisFillHistograms() *** no primary ***:  label %d \n", label);
 811             continue;
 812           }
 813
 814           //eprim   = aodprimary->E();
 815           //ptprim  = aodprimary->Pt();
 816
 817         }
 818
 819         Int_t tag =ph->GetTag();
 820
 821         if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton))
 822         {
 823
 824           if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion))
 825           {
 826
 827             fhPtConversionTagged ->Fill(ptcluster);
 828
 829           }
 830         }
 831         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiNeutron))
 832         {
 833
 834           fhPtAntiNeutronTagged ->Fill(ptcluster);
 835
 836         }
 837         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiProton))
 838         {
 839           fhPtAntiProtonTagged ->Fill(ptcluster);
 840
 841         }
 842
 843         else {
 844           fhPtUnknownTagged ->Fill(ptcluster);
 845
 846         }
 847
 848       }//Histograms with MC
 849     }// tagged by conversion
 850   }// aod loop
 851
 852 }
 853
 854
 855 //________________________________________________________
 856 void AliAnaPhotonConvInCalo::Print(const Option_t * opt) const
 857 {
 858   //Print some relevant parameters set for the analysis
 859
 860   if(! opt)
 861     return;
 862
 863   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
 864   AliAnaCaloTrackCorrBaseClass::Print(" ");
 865
 866   printf("Add conversion pair to AOD           = %d\n",fAddConvertedPairsToAOD);
 867   printf("Conversion pair mass cut             = %f\n",fMassCut);
 868   printf("Conversion selection cut : A < %1.2f; %1.3f < Dphi < %1.3f; Deta < %1.3f\n",
 869          fConvAsymCut,fConvDPhiMinCut, fConvDPhiMaxCut, fConvDEtaCut);
 870
 871   printf("    \n") ;
 872
 873 }