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 fMomentum(),                       fProdVertex(),
  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     {
 477       //Check if set previously as converted couple, if so skip its use.
 478       if (indexConverted[jaod]) continue;
 479       //printf("Check Conversion indeces %d and %d\n",iaod,jaod);
 480       calo2 =  (AliAODPWG4Particle*) (GetOutputAODBranch()->At(jaod));
 481
 482       //................................................
 483       //Get mass of pair, if small, take this pair.
 484       Float_t pairM     = calo->GetPairMass(calo2);
 485       //printf("\t both in calo, mass %f, cut %f\n",pairM,fMassCut);
 486       if(pairM < fMassCut){
 487         calo->SetTagged(kFALSE);
 488         id2 = calo2->GetCaloLabel(0);
 489         Float_t asymmetry = TMath::Abs(calo->E()-calo2->E())/(calo->E()+calo2->E());
 490         Float_t dPhi      = (calo->Momentum())->Phi()-(calo2->Momentum())->Phi();
 491         Float_t dEta      = (calo->Momentum())->Eta()-(calo2->Momentum())->Eta();
 492
 493         //...............................................
 494         //Fill few histograms with kinematics of the pair
 495         //FIXME, move all this to MakeAnalysisFillHistograms ...
 496
 497         fhConvDeltaEta   ->Fill( pairM, dPhi      );
 498         fhConvDeltaPhi   ->Fill( pairM, dEta      );
 499         fhConvAsym       ->Fill( pairM, asymmetry );
 500         fhConvDeltaEtaPhi->Fill( dEta , dPhi      );
 501         fhConvPt         ->Fill( pairM, (calo->Momentum())->Pt()+(calo2->Momentum())->Pt());
 502
 503         //Estimate conversion distance, T. Awes, M. Ivanov
 504         //Under the assumption that the pair has zero mass, and that each electron
 505         //of the pair has the same momentum, they will each have the same bend radius
 506         //given by R=p/(qB) = p / (300 B) with p in [MeV/c], B in [Tesla] and R in [m].
 507         //With nominal ALICE magnet current of 30kA B=0.5T, and so with E_cluster=p,
 508         //R = E/1.5 [cm].  Under these assumptions, the distance from the conversion
 509         //point to the MCEal can be related to the separation distance, L=2y, on the MCEal
 510         //as d = sqrt(R^2 -(R-y)^2) = sqrt(2Ry - y^2). And since R>>y we can write as
 511         //d = sqrt(E*L/1.5) where E is the cluster energy and L is the distance in cm between
 512         //the clusters.
 513
 514         TObjArray * clusters    = 0;
 515         if(calo->GetDetectorTag() == kEMCAL) clusters = GetEMCALClusters();
 516         else                                 clusters = GetPHOSClusters();
 517
 518
 519         Int_t iclus = -1;
 520         AliVCluster *cluster1 = FindCluster(clusters,calo ->GetCaloLabel(0),iclus);
 521         AliVCluster *cluster2 = FindCluster(clusters,calo2->GetCaloLabel(0),iclus);
 522
 523         Float_t pos1[3];
 524         cluster1->GetPosition(pos1);
 525         Float_t pos2[3];
 526         cluster2->GetPosition(pos2);
 527         Float_t clustDist = TMath::Sqrt((pos1[0]-pos2[0])*(pos1[0]-pos2[0])+
 528                                         (pos1[1]-pos2[1])*(pos1[1]-pos2[1])+
 529                                         (pos1[2]-pos2[2])*(pos1[2]-pos2[2]));
 530
 531         Float_t convDist  = TMath::Sqrt(calo->E() *clustDist*0.01/0.15);
 532         Float_t convDist2 = TMath::Sqrt(calo2->E()*clustDist*0.01/0.15);
 533         //printf("l = %f, e1 = %f, d1=%f, e2 = %f, d2=%f\n",clustDist,calo->E(),convDist,calo2->E(),convDist2);
 534         if(GetDebug() > 2)
 535           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",
 536                  pairM,fMassCut,fConvDPhiMinCut, dPhi, fConvDPhiMaxCut, dEta, fConvDEtaCut, asymmetry, fConvAsymCut,
 537                  calo->GetCaloLabel(0),calo->E(),GetCaloUtils()->GetModuleNumber(calo,GetReader()->GetInputEvent()), calo->Eta(), calo->Phi(),
 538                  id2, calo2->E(), GetCaloUtils()->GetModuleNumber(calo2,GetReader()->GetInputEvent()),calo2->Eta(), calo2->Phi());
 539
 540         fhConvDistEta ->Fill(calo ->Eta(),convDist );
 541         fhConvDistEta ->Fill(calo2->Eta(),convDist2);
 542         fhConvDistEn  ->Fill(calo ->E(),  convDist );
 543         fhConvDistEn  ->Fill(calo2->E(),  convDist2);
 544         fhConvDistMass->Fill(pairM, convDist );
 545         //dEta cut
 546         if(dEta<0.05){
 547           fhConvDistEtaCutEta ->Fill(calo->Eta(), convDist );
 548           fhConvDistEtaCutEta ->Fill(calo2->Eta(),convDist2);
 549           fhConvDistEnCutEta  ->Fill(calo->E(),   convDist );
 550           fhConvDistEnCutEta  ->Fill(calo2->E(),  convDist2);
 551           fhConvDistMassCutEta->Fill(pairM, convDist );
 552           //mass cut
 553           if(pairM<0.01){//10 MeV
 554             fhConvDistEtaCutMass ->Fill(calo ->Eta(), convDist );
 555             fhConvDistEtaCutMass ->Fill(calo2->Eta(), convDist2);
 556             fhConvDistEnCutMass  ->Fill(calo ->E(),   convDist );
 557             fhConvDistEnCutMass  ->Fill(calo2->E(),   convDist2);
 558             // asymmetry cut
 559             if(asymmetry<0.1){
 560               fhConvDistEtaCutAsy ->Fill(calo ->Eta(), convDist );
 561               fhConvDistEtaCutAsy ->Fill(calo2->Eta(), convDist2);
 562               fhConvDistEnCutAsy  ->Fill(calo ->E(),   convDist );
 563               fhConvDistEnCutAsy  ->Fill(calo2->E(),   convDist2);
 564             }//asymmetry cut
 565           }//mass cut
 566         }//dEta cut
 567
 568         //...............................................
 569         //Select pairs in a eta-phi window
 570         if(TMath::Abs(dEta) < fConvDEtaCut    &&
 571            TMath::Abs(dPhi) < fConvDPhiMaxCut &&
 572            TMath::Abs(dPhi) > fConvDPhiMinCut &&
 573            asymmetry        < fConvAsymCut       ){
 574           indexConverted[iaod] = kTRUE;
 575           indexConverted[jaod] = kTRUE;
 576           bConverted           = kTRUE;
 577         }
 578         //printf("Accepted? %d\n",bConverted);
 579         //...........................................
 580         //Fill more histograms, simulated data
 581         //FIXME, move all this to MakeAnalysisFillHistograms ...
 582         if(IsDataMC()){
 583
 584           //Check the origin of the pair, look for conversion, antinucleons or jet correlations (strings)
 585           Int_t ancPDG    = 0;
 586           Int_t ancStatus = 0;
 587           Int_t ancLabel  = GetMCAnalysisUtils()->CheckCommonAncestor(cluster1->GetLabel(), cluster2->GetLabel(),
 588                                                                       GetReader(), ancPDG, ancStatus, fMomentum, fProdVertex);
 589
 590           // printf("AliAnaPhotonConvInCalo::MakeAnalysisFillHistograms() - Common ancestor label %d, pdg %d, name %s, status %d; \n",
 591           //                          ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus);
 592
 593           Int_t tag1 = calo ->GetTag();
 594           Int_t tag2 = calo2->GetTag();
 595           if(GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCConversion)){
 596             if(GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCConversion) && (ancPDG==22 || TMath::Abs(ancPDG)==11) && ancLabel > -1){
 597               fhConvDeltaEtaMCConversion   ->Fill( pairM, dEta      );
 598               fhConvDeltaPhiMCConversion   ->Fill( pairM, dPhi      );
 599               fhConvAsymMCConversion       ->Fill( pairM, asymmetry );
 600               fhConvDeltaEtaPhiMCConversion->Fill( dEta , dPhi      );
 601               fhConvPtMCConversion         ->Fill( pairM, calo->Pt()+calo2->Pt());
 602               fhConvDispersionMCConversion ->Fill( cluster1->GetDispersion(), cluster2->GetDispersion());
 603               fhConvM02MCConversion        ->Fill( cluster1->GetM02(), cluster2->GetM02());
 604               fhConvDistMCConversion       ->Fill( convDist , fProdVertex.Mag() );
 605               fhConvDistMCConversion       ->Fill( convDist2, fProdVertex.Mag() );
 606
 607               if(dEta<0.05 && pairM<0.01 && asymmetry<0.1){
 608                 fhConvDistMCConversionCuts->Fill( convDist , fProdVertex.Mag() );
 609                 fhConvDistMCConversionCuts->Fill( convDist2, fProdVertex.Mag() );
 610               }
 611
 612             }
 613           }
 614           else if(GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCAntiNeutron)){
 615             if(GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCAntiNeutron) && ancPDG==-2112 && ancLabel > -1){
 616               fhConvDeltaEtaMCAntiNeutron    ->Fill( pairM, dEta      );
 617               fhConvDeltaPhiMCAntiNeutron    ->Fill( pairM, dPhi      );
 618               fhConvAsymMCAntiNeutron        ->Fill( pairM, asymmetry );
 619               fhConvDeltaEtaPhiMCAntiNeutron ->Fill( dEta , dPhi      );
 620               fhConvPtMCAntiNeutron          ->Fill( pairM, calo->Pt()+calo2->Pt());
 621               fhConvDispersionMCAntiNeutron  ->Fill( cluster1->GetDispersion(), cluster2->GetDispersion());
 622               fhConvM02MCAntiNeutron         ->Fill( cluster1->GetM02(), cluster2->GetM02());
 623             }
 624           }
 625           else if(GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCAntiProton)){
 626             if(GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCAntiProton) && ancPDG==-2212 && ancLabel > -1){
 627               fhConvDeltaEtaMCAntiProton    ->Fill( pairM, dEta      );
 628               fhConvDeltaPhiMCAntiProton    ->Fill( pairM, dPhi      );
 629               fhConvAsymMCAntiProton        ->Fill( pairM, asymmetry );
 630               fhConvDeltaEtaPhiMCAntiProton ->Fill( dEta , dPhi      );
 631               fhConvPtMCAntiProton          ->Fill( pairM, calo->Pt()+calo2->Pt());
 632               fhConvDispersionMCAntiProton  ->Fill( cluster1->GetDispersion(), cluster2->GetDispersion());
 633               fhConvM02MCAntiProton         ->Fill( cluster1->GetM02(), cluster2->GetM02());
 634             }
 635           }
 636
 637           //Pairs coming from fragmenting pairs.
 638           if(ancPDG < 22 && ancLabel > 7 && (ancStatus == 11 || ancStatus == 12) ){
 639             fhConvDeltaEtaMCString    ->Fill( pairM, dPhi);
 640             fhConvDeltaPhiMCString    ->Fill( pairM, dPhi);
 641             fhConvAsymMCString        ->Fill( pairM, TMath::Abs(calo->E()-calo2->E())/(calo->E()+calo2->E()) );
 642             fhConvDeltaEtaPhiMCString ->Fill( dPhi,  dPhi);
 643             fhConvPtMCString          ->Fill( pairM, calo->Pt()+calo2->Pt());
 644             fhConvDispersionMCString  ->Fill( cluster1->GetDispersion(), cluster2->GetDispersion());
 645             fhConvM02MCString         ->Fill( cluster1->GetM02(), cluster2->GetM02());
 646           }
 647
 648         }// Data MC
 649
 650         break;
 651       }
 652
 653     }//Mass loop
 654
 655     //..........................................................................................................
 656     //Pair selected as converted, remove both clusters or recombine them into a photon and put them in the AOD
 657     if(bConverted)
 658     {
 659       //Add to AOD
 660       if(fAddConvertedPairsToAOD){
 661         //Create AOD of pair analysis
 662         fMomentum = *(calo->Momentum())+*(calo2->Momentum());
 663         AliAODPWG4Particle aodpair = AliAODPWG4Particle(fMomentum);
 664         aodpair.SetLabel(calo->GetLabel());
 665
 666         //printf("Index %d, Id %d\n",iaod, calo->GetID());
 667         //Set the indeces of the original caloclusters
 668         aodpair.SetCaloLabel(calo->GetCaloLabel(0),id2);
 669         aodpair.SetDetectorTag(calo->GetDetectorTag());
 670         aodpair.SetIdentifiedParticleType(calo->GetIdentifiedParticleType());
 671         aodpair.SetTag(calo ->GetTag());
 672         aodpair.SetTagged(kTRUE);
 673         //Add AOD with pair object to aod branch
 674         AddAODParticle(aodpair);
 675         //printf("\t \t both added pair\n");
 676       }
 677
 678       //Do not add the current calocluster
 679       if(!fRemoveConvertedPair)
 680       {
 681         //printf("TAGGED\n");
 682         //Tag this cluster as likely conversion
 683         calo->SetTagged(kTRUE);
 684       }
 685     }//converted pair
 686
 687   }// main loop
 688
 689   // Remove entries identified as conversion electrons
 690   // Revise if this is OK
 691   if(fRemoveConvertedPair || fAddConvertedPairsToAOD){
 692     for(Int_t iaod = 0; iaod < naod ; iaod++)
 693       if(indexConverted[iaod])GetOutputAODBranch()->RemoveAt(iaod);
 694     GetOutputAODBranch()->Compress();
 695   }
 696
 697   delete [] indexConverted;
 698
 699   if(GetDebug() > 1) printf("AliAnaPhotonConvInCalo::MakeAnalysisFillAOD()  End fill AODs, with %d entries \n",GetOutputAODBranch()->GetEntriesFast());
 700
 701 }
 702
 703 //____________________________________________________
 704 void  AliAnaPhotonConvInCalo::MakeAnalysisFillHistograms()
 705 {
 706   //Fill histograms
 707
 708   //-------------------------------------------------------------------
 709   // Access MC information in stack if requested, check that it exists.
 710   AliStack         * stack       = 0x0;
 711   TParticle        * primary     = 0x0;
 712   TClonesArray     * mcparticles = 0x0;
 713   AliAODMCParticle * aodprimary  = 0x0;
 714
 715   if(IsDataMC()){
 716
 717     if(GetReader()->ReadStack()){
 718       stack =  GetMCStack() ;
 719       if(!stack) {
 720         printf("AliAnaPhotonConvInCalo::MakeAnalysisFillHistograms() - Stack not available, is the MC handler called? STOP\n");
 721         abort();
 722       }
 723
 724     }
 725     else if(GetReader()->ReadAODMCParticles()){
 726
 727       //Get the list of MC particles
 728       mcparticles = GetReader()->GetAODMCParticles();
 729       if(!mcparticles && GetDebug() > 0)
 730       {
 731         printf("AliAnaPhotonConvInCalo::MakeAnalysisFillHistograms() -  Standard MCParticles not available!\n");
 732       }
 733     }
 734   }// is data and MC
 735
 736   //----------------------------------
 737   //Loop on stored AOD photons
 738   Int_t naod = GetOutputAODBranch()->GetEntriesFast();
 739   if(GetDebug() > 0) printf("AliAnaPhotonConvInCalo::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod);
 740
 741   for(Int_t iaod = 0; iaod < naod ; iaod++){
 742     AliAODPWG4Particle* ph =  (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
 743     //Int_t pdg = ph->GetIdentifiedParticleType();
 744
 745     if(ph->IsTagged()){
 746
 747       if(GetDebug() > 2)
 748         printf("AliAnaPhotonConvInCalo::MakeAnalysisFillHistograms() - ID Photon: pt %f, phi %f, eta %f\n", ph->Pt(),ph->Phi(),ph->Eta()) ;
 749       //................................
 750       //Fill photon histograms
 751       Float_t ptcluster  = ph->Pt();
 752       Float_t phicluster = ph->Phi();
 753       Float_t etacluster = ph->Eta();
 754       Float_t ecluster   = ph->E();
 755
 756       fhPtPhotonConv->Fill(ptcluster);
 757       if(ecluster > 0.5)        fhEtaPhiPhotonConv  ->Fill(etacluster, phicluster);
 758       else if(GetMinPt() < 0.5) fhEtaPhi05PhotonConv->Fill(etacluster, phicluster);
 759
 760
 761       //.......................................
 762       //Play with the MC data if available
 763       if(IsDataMC()){
 764
 765
 766         //....................................................................
 767         // Access MC information in stack if requested, check that it exists.
 768         Int_t label =ph->GetLabel();
 769         if(label < 0) {
 770           if(GetDebug() > 1) printf("AliAnaPhotonConvInCalo::MakeAnalysisFillHistograms() *** bad label ***:  label %d \n", label);
 771           continue;
 772         }
 773
 774         //Float_t eprim   = 0;
 775         //Float_t ptprim  = 0;
 776         if(GetReader()->ReadStack()){
 777
 778           if(label >=  stack->GetNtrack()) {
 779             if(GetDebug() > 2)  printf("AliAnaPhotonConvInCalo::MakeAnalysisFillHistograms() *** large label ***:  label %d, n tracks %d \n", label, stack->GetNtrack());
 780             continue ;
 781           }
 782
 783           primary = stack->Particle(label);
 784           if(!primary){
 785             printf("AliAnaPhotonConvInCalo::MakeAnalysisFillHistograms() *** no primary ***:  label %d \n", label);
 786             continue;
 787           }
 788           //eprim   = primary->Energy();
 789           //ptprim  = primary->Pt();
 790
 791         }
 792         else if(GetReader()->ReadAODMCParticles()){
 793           //Check which is the input
 794           if(ph->GetInputFileIndex() == 0){
 795             if(!mcparticles) continue;
 796             if(label >=  mcparticles->GetEntriesFast()) {
 797               if(GetDebug() > 2)  printf("AliAnaPhotonConvInCalo::MakeAnalysisFillHistograms() *** large label ***:  label %d, n tracks %d \n",
 798                                          label, mcparticles->GetEntriesFast());
 799               continue ;
 800             }
 801             //Get the particle
 802             aodprimary = (AliAODMCParticle*) mcparticles->At(label);
 803
 804           }
 805
 806           if(!aodprimary){
 807             printf("AliAnaPhotonConvInCalo::MakeAnalysisFillHistograms() *** no primary ***:  label %d \n", label);
 808             continue;
 809           }
 810
 811           //eprim   = aodprimary->E();
 812           //ptprim  = aodprimary->Pt();
 813
 814         }
 815
 816         Int_t tag =ph->GetTag();
 817
 818         if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton))
 819         {
 820
 821           if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion))
 822           {
 823
 824             fhPtConversionTagged ->Fill(ptcluster);
 825
 826           }
 827         }
 828         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiNeutron))
 829         {
 830
 831           fhPtAntiNeutronTagged ->Fill(ptcluster);
 832
 833         }
 834         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiProton))
 835         {
 836           fhPtAntiProtonTagged ->Fill(ptcluster);
 837
 838         }
 839
 840         else {
 841           fhPtUnknownTagged ->Fill(ptcluster);
 842
 843         }
 844
 845       }//Histograms with MC
 846     }// tagged by conversion
 847   }// aod loop
 848
 849 }
 850
 851
 852 //________________________________________________________
 853 void AliAnaPhotonConvInCalo::Print(const Option_t * opt) const
 854 {
 855   //Print some relevant parameters set for the analysis
 856
 857   if(! opt)
 858     return;
 859
 860   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
 861   AliAnaCaloTrackCorrBaseClass::Print(" ");
 862
 863   printf("Add conversion pair to AOD           = %d\n",fAddConvertedPairsToAOD);
 864   printf("Conversion pair mass cut             = %f\n",fMassCut);
 865   printf("Conversion selection cut : A < %1.2f; %1.3f < Dphi < %1.3f; Deta < %1.3f\n",
 866          fConvAsymCut,fConvDPhiMinCut, fConvDPhiMaxCut, fConvDEtaCut);
 867
 868   printf("    \n") ;
 869
 870 }