PWG4/PartCorrDep/AliAnalysisTaskTaggedPhotons.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 is 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 /* $Id: */
  17
  18 //_________________________________________________________________________
  19 // Analysis for Tagged Photons
  20 // Prepares all necessary histograms for calculation of
  21 // the yield of pi0 decay photon in calorimeter:
  22 // Marks photons which makes pi0 with some other and
  23 // fill invariant mass distributions for estimate background below pi0
  24 // peak so that estimate proportion of fake pairs.
  25 // Fills as well controll MC histograms with clasification of the photon origin
  26 // and check of the proportion of truly tagged photons.
  27 //
  28 //
  29 //*-- Dmitry Blau
  30 //////////////////////////////////////////////////////////////////////////////
  31
  32 #include <TH1.h>
  33 #include <TH2.h>
  34 #include <TH3.h>
  35 #include <TFile.h>
  36 #include <TROOT.h>
  37
  38 #include "AliAnalysisTaskTaggedPhotons.h"
  39 #include "AliAnalysisManager.h"
  40 #include "AliESDVertex.h"
  41 #include "AliESDEvent.h"
  42 #include "AliAODEvent.h"
  43 #include "AliESDCaloCluster.h"
  44 #include "AliAODPWG4Particle.h"
  45 #include "AliAnalysisManager.h"
  46 #include "AliLog.h"
  47 #include "TGeoManager.h"
  48 #include "AliMCAnalysisUtils.h"
  49 #include "AliMCEventHandler.h"
  50 #include "AliMCEvent.h"
  51 #include "AliStack.h"
  52 #include "AliPHOSGeoUtils.h"
  53 #include "AliEMCALGeoUtils.h"
  54
  55
  56
  57 //______________________________________________________________________________
  58 AliAnalysisTaskTaggedPhotons::AliAnalysisTaskTaggedPhotons() : AliAnalysisTaskSE(),
  59   fPHOSgeom(0x0),
  60   fEMCALgeom(0x0),
  61   fStack(0x0),fPHOS(1),
  62   fPhotonId(1.0),fMinEnergyCut(0.4),
  63   fPi0MeanP0(0.136),fPi0MeanP1(0.0),fPi0MeanP2(0.0),fPi0MeanP3(0.0),
  64   fPi0SigmaP0(0.004508),fPi0SigmaP1(0.005497),fPi0SigmaP2(0.00000006),
  65   fZmax(0.),fZmin(0.),fPhimax(0.),fPhimin(0.),
  66   fOutputList(0x0),fEventList(0x0),
  67   fhRecAllArea1(0x0),fhRecAllArea2(0x0),fhRecAllArea3(0x0),
  68   fhRecPhoton(0x0),fhRecOther(0x0),fhRecPhotPi0(0x0),fhRecPhotEta(0x0),
  69   fhRecPhotOmega(0x0),fhRecPhotEtapr(0x0),fhRecPhotConv(0x0),fhRecPhotHadron(0x0),
  70   fhRecPhotDirect(0x0),fhRecPhotOther(0x0),
  71   fhDecWMCPartner(0x0),fhDecWMissedPartnerNotPhoton(0x0),fhDecWMissedPartnerAll(0x0),
  72   fhDecWMissedPartnerEmin(0x0),fhDecWMissedPartnerConv(0x0),fhDecWMissedPartnerStack(0x0),
  73   fhDecWMissedPartnerGeom0(0x0),
  74   fhDecWMissedPartnerGeom1(0x0),fhDecWMissedPartnerGeom2(0x0),fhDecWMissedPartnerGeom3(0x0),
  75   fhPartnerMCReg(0x0),fhPartnerMissedEmin(0x0),fhPartnerMissedConv(0x0),fhPartnerMissedGeo(0x0),
  76   fhTaggedAll(0x0),fhTaggedArea1(0x0),fhTaggedArea2(0x0),fhTaggedArea3(0x0),fhTaggedMult(0x0),
  77   fhTaggedMCTrue(0x0),fhMCMissedTagging(0x0),fhMCFakeTagged(0x0),
  78   fhMCMissedTaggingMass(0x0),
  79   fhConversionRadius(0x0),fhInteractionRadius(0x0),fhEvents(0x0)
  80 {
  81   //Default constructor
  82   for (Int_t i = 0; i < 4; i++)
  83     {
  84       fhRecAll[i] = 0;
  85       fhRecPhotonPID[i] = 0;
  86       fhRecOtherPID[i] = 0;
  87       fhTaggedPID[i] = 0;
  88       fhInvMassReal[i] = 0;
  89       fhInvMassMixed[i] = 0;
  90     }
  91
  92 }
  93 //______________________________________________________________________________
  94 AliAnalysisTaskTaggedPhotons::AliAnalysisTaskTaggedPhotons(const char *name) :
  95   AliAnalysisTaskSE(name),
  96   fPHOSgeom(0x0),
  97   fEMCALgeom(0x0),
  98   fStack(0x0),fPHOS(1),
  99   fPhotonId(1.0),fMinEnergyCut(0.4),
 100   fPi0MeanP0(0.136),fPi0MeanP1(0.0),fPi0MeanP2(0.0),fPi0MeanP3(0.0),
 101   fPi0SigmaP0(0.004508),fPi0SigmaP1(0.005497),fPi0SigmaP2(0.00000006),
 102   fZmax(0.),fZmin(0.),fPhimax(0.),fPhimin(0.),
 103   fOutputList(0x0),fEventList(0x0),
 104   fhRecAllArea1(0x0),fhRecAllArea2(0x0),fhRecAllArea3(0x0),
 105   fhRecPhoton(0x0),fhRecOther(0x0),fhRecPhotPi0(0x0),fhRecPhotEta(0x0),
 106   fhRecPhotOmega(0x0),fhRecPhotEtapr(0x0),fhRecPhotConv(0x0),fhRecPhotHadron(0x0),
 107   fhRecPhotDirect(0x0),fhRecPhotOther(0x0),
 108   fhDecWMCPartner(0x0),fhDecWMissedPartnerNotPhoton(0x0),fhDecWMissedPartnerAll(0x0),
 109   fhDecWMissedPartnerEmin(0x0),fhDecWMissedPartnerConv(0x0),fhDecWMissedPartnerStack(0x0),
 110   fhDecWMissedPartnerGeom0(0x0),
 111   fhDecWMissedPartnerGeom1(0x0),fhDecWMissedPartnerGeom2(0x0),fhDecWMissedPartnerGeom3(0x0),
 112   fhPartnerMCReg(0x0),fhPartnerMissedEmin(0x0),fhPartnerMissedConv(0x0),fhPartnerMissedGeo(0x0),
 113   fhTaggedAll(0x0),fhTaggedArea1(0x0),fhTaggedArea2(0x0),fhTaggedArea3(0x0),fhTaggedMult(0x0),
 114   fhTaggedMCTrue(0x0),fhMCMissedTagging(0x0),fhMCFakeTagged(0x0),
 115   fhMCMissedTaggingMass(0x0),
 116   fhConversionRadius(0x0),fhInteractionRadius(0x0),fhEvents(0x0)
 117 {
 118   // Constructor.
 119   for (Int_t i = 0; i < 4; i++)
 120     {
 121       fhRecAll[i] = 0;
 122       fhRecPhotonPID[i] = 0;
 123       fhRecOtherPID[i] = 0;
 124       fhTaggedPID[i] = 0;
 125       fhInvMassReal[i] = 0;
 126       fhInvMassMixed[i] = 0;
 127     }
 128
 129   // Output slots
 130   DefineOutput(1,  TList::Class()) ;
 131 }
 132
 133 //____________________________________________________________________________
 134 AliAnalysisTaskTaggedPhotons::AliAnalysisTaskTaggedPhotons(const AliAnalysisTaskTaggedPhotons& ap) :
 135   AliAnalysisTaskSE(ap.GetName()),
 136   fPHOSgeom(0x0),
 137   fEMCALgeom(0x0),
 138   fStack(0x0),fPHOS(ap.fPHOS),
 139   fPhotonId(ap.fPhotonId),fMinEnergyCut(ap.fMinEnergyCut),
 140   fPi0MeanP0(ap.fPi0MeanP0),fPi0MeanP1(ap.fPi0MeanP1),fPi0MeanP2(ap.fPi0MeanP2),fPi0MeanP3(ap.fPi0MeanP3),
 141   fPi0SigmaP0(ap.fPi0SigmaP0),fPi0SigmaP1(ap.fPi0SigmaP1),fPi0SigmaP2(ap.fPi0SigmaP2),
 142   fZmax(ap.fZmax),fZmin(ap.fZmin),fPhimax(ap.fPhimax),fPhimin(ap.fPhimin),
 143   fOutputList(0x0),fEventList(0x0),
 144   fhRecAllArea1(0x0),fhRecAllArea2(0x0),fhRecAllArea3(0x0),
 145   fhRecPhoton(0x0),fhRecOther(0x0),fhRecPhotPi0(0x0),fhRecPhotEta(0x0),
 146   fhRecPhotOmega(0x0),fhRecPhotEtapr(0x0),fhRecPhotConv(0x0),fhRecPhotHadron(0x0),
 147   fhRecPhotDirect(0x0),fhRecPhotOther(0x0),
 148   fhDecWMCPartner(0x0),fhDecWMissedPartnerNotPhoton(0x0),fhDecWMissedPartnerAll(0x0),
 149   fhDecWMissedPartnerEmin(0x0),fhDecWMissedPartnerConv(0x0),fhDecWMissedPartnerStack(0x0),
 150   fhDecWMissedPartnerGeom0(0x0),
 151   fhDecWMissedPartnerGeom1(0x0),fhDecWMissedPartnerGeom2(0x0),fhDecWMissedPartnerGeom3(0x0),
 152   fhPartnerMCReg(0x0),fhPartnerMissedEmin(0x0),fhPartnerMissedConv(0x0),
 153   fhPartnerMissedGeo(0x0),
 154   fhTaggedAll(0x0),fhTaggedArea1(0x0),fhTaggedArea2(0x0),fhTaggedArea3(0x0),fhTaggedMult(0x0),
 155   fhTaggedMCTrue(0x0),fhMCMissedTagging(0x0),fhMCFakeTagged(0x0),
 156   fhMCMissedTaggingMass(0x0),
 157   fhConversionRadius(0x0),fhInteractionRadius(0x0),fhEvents(0x0)
 158 {
 159   // cpy ctor
 160   for (Int_t i = 0; i < 4; i++)
 161     {
 162       fhRecAll[i] = 0;
 163       fhRecPhotonPID[i] = 0;
 164       fhRecOtherPID[i] = 0;
 165       fhTaggedPID[i] = 0;
 166       fhInvMassReal[i] = 0;
 167       fhInvMassMixed[i] = 0;
 168     }
 169
 170 }
 171
 172 //_____________________________________________________________________________
 173 AliAnalysisTaskTaggedPhotons& AliAnalysisTaskTaggedPhotons::operator = (const AliAnalysisTaskTaggedPhotons& ap)
 174 {
 175 // assignment operator
 176
 177   this->~AliAnalysisTaskTaggedPhotons();
 178   new(this) AliAnalysisTaskTaggedPhotons(ap);
 179   return *this;
 180 }
 181
 182 //______________________________________________________________________________
 183 AliAnalysisTaskTaggedPhotons::~AliAnalysisTaskTaggedPhotons()
 184 {
 185   // dtor
 186   if(fOutputList) {
 187     fOutputList->Clear() ;
 188     delete fOutputList ;
 189   }
 190 }
 191
 192
 193 //________________________________________________________________________
 194 void AliAnalysisTaskTaggedPhotons::UserCreateOutputObjects()
 195 {
 196
 197
 198   //Load geometry
 199   //if file "geometry.root" exists, load it
 200   //otherwise use misaligned matrixes stored in ESD
 201
 202   // Create the outputs containers
 203
 204   OpenFile(1) ;
 205   const Int_t nPtBins=52 ;
 206   Double_t ptBins[nPtBins+1] ;
 207   for(Int_t i=0;i<=20;i++)ptBins[i]=0.1*i ; //0-2 GeV:  0.1 GeV/bin
 208   for(Int_t i=21;i<=30;i++)ptBins[i]=2.+0.2*(i-20) ; //2-4 GeV:  0.2 GeV/bin
 209   for(Int_t i=31;i<=38;i++)ptBins[i]=4.+0.5*(i-30) ; //4-8 GeV:  0.5 GeV/bin
 210   for(Int_t i=39;i<=42;i++)ptBins[i]=8.+1.0*(i-38) ; //8-12 GeV: 1. GeV/bin
 211   for(Int_t i=43;i<=52;i++)ptBins[i]=12.+2.0*(i-42) ; //12-30 GeV: 2. GeV/bin
 212
 213
 214   //Reconstructed spectra
 215     fhRecAll[0]   = new TH1D("fhRecAll0","Spectrum of all reconstructed particles, no PID", nPtBins, ptBins ) ;
 216     fhRecAll[1]   = new TH1D("fhRecAll1","Spectrum of all reconstructed particles, PID=1", nPtBins, ptBins ) ;
 217     fhRecAll[2]   = new TH1D("fhRecAll2","Spectrum of all reconstructed particles, PID=2", nPtBins, ptBins ) ;
 218     fhRecAll[3]   = new TH1D("fhRecAll3","Spectrum of all reconstructed particles, PID=3", nPtBins, ptBins ) ;
 219
 220     fhRecAllArea1 = new TH1D("fhRecAllArea1","Spectrum of rec particles in Fid. Area 1", nPtBins, ptBins ) ;
 221     fhRecAllArea2 = new TH1D("fhRecAllArea2","Spectrum of rec particles in Fid. Area 2", nPtBins, ptBins ) ;
 222     fhRecAllArea3 = new TH1D("fhRecAllArea3","Spectrum of rec particles in Fid. Area 3", nPtBins, ptBins ) ;
 223
 224     //Sort registered particles spectra according MC information
 225     fhRecPhoton   = new TH1D("fhRecPhoton","Spectrum of rec. with primary==22 and no PID criteria", nPtBins, ptBins ) ;
 226     fhRecOther    = new TH1D("fhRecOther"," Spectrum of rec. with primary!=22 and no PID criteria", nPtBins, ptBins ) ;
 227     fhRecPhotonPID[0]= new TH1D("fhRecPhotonPID0","Spectrum of rec. with primary==22 and PID=0", nPtBins, ptBins ) ;
 228     fhRecPhotonPID[1]= new TH1D("fhRecPhotonPID1","Spectrum of rec. with primary==22 and PID=1", nPtBins, ptBins ) ;
 229     fhRecPhotonPID[2]= new TH1D("fhRecPhotonPID2","Spectrum of rec. with primary==22 and PID=2", nPtBins, ptBins ) ;
 230     fhRecPhotonPID[3]= new TH1D("fhRecPhotonPID3","Spectrum of rec. with primary==22 and PID=3", nPtBins, ptBins ) ;
 231     fhRecOtherPID[0]= new TH1D("fhRecOtherPID0","Spectrum of rec. with primary!=22 and PID=0", nPtBins, ptBins ) ;
 232     fhRecOtherPID[1]= new TH1D("fhRecOtherPID1","Spectrum of rec. with primary!=22 and PID=1", nPtBins, ptBins ) ;
 233     fhRecOtherPID[2]= new TH1D("fhRecOtherPID2","Spectrum of rec. with primary!=22 and PID=2", nPtBins, ptBins ) ;
 234     fhRecOtherPID[3]= new TH1D("fhRecOtherPID3","Spectrum of rec. with primary!=22 and PID=3", nPtBins, ptBins ) ;
 235     fhRecPhotPi0    = new TH1D("fhRecPhotPi0","Spectrum of rec. photons from pi0 decays", nPtBins, ptBins ) ;
 236     fhRecPhotEta    = new TH1D("fhRecPhotEta","Spectrum of rec. photons from eta decays", nPtBins, ptBins ) ;
 237     fhRecPhotOmega  = new TH1D("fhRecPhotOmega","Spectrum of rec. photons from omega decays", nPtBins, ptBins ) ;
 238     fhRecPhotEtapr  = new TH1D("fhRecPhotEtapr","Spectrum of rec. photons from eta prime decays", nPtBins, ptBins ) ;
 239     fhRecPhotConv   = new TH1D("fhRecPhotConv"," Spectrum of rec. photons from conversion", nPtBins, ptBins ) ;
 240     fhRecPhotHadron = new TH1D("fhRecPhotHadron","Spectrum of rec. photons from hadron-matter interactions", nPtBins, ptBins ) ;
 241     fhRecPhotDirect = new TH1D("fhRecPhotDirect","Spectrum of rec. photons direct or no primary", nPtBins, ptBins ) ;
 242     fhRecPhotOther  = new TH1D("fhRecPhotOther","Spectrum of rec. photons from other hadron decays", nPtBins, ptBins ) ;
 243
 244     //MC tagging: reasons of partner loss etc.
 245     fhDecWMCPartner        = new TH1D("fhDecWMCPartner","pi0 decay photon which partner should be registered according to MC", nPtBins, ptBins ) ;
 246     fhDecWMissedPartnerNotPhoton = new TH1D("fhDecWMissedPartnerNotPhoton","Decay photon with missed non-photon partner", nPtBins, ptBins ) ;
 247     fhDecWMissedPartnerAll  = new TH1D("fhDecWMissedPartnerAll","Decay photons with partner missed due to some reason", nPtBins, ptBins ) ;
 248     fhDecWMissedPartnerEmin = new TH1D("fhDecWMissedPartnerEmin","Decay photons with partner missed due to low energy", nPtBins, ptBins ) ;
 249     fhDecWMissedPartnerConv = new TH1D("fhDecWMissedPartnerConv","Decay photons with partner missed due to conversion", nPtBins, ptBins ) ;
 250     fhDecWMissedPartnerStack= new TH1D("fhDecWMissedPartnerStack","Decay photons with partner not in Stack", nPtBins, ptBins ) ;
 251     fhDecWMissedPartnerGeom0 = new TH1D("fhDecWMissedPartnerGeom0","Decay photons with partner missed due geometry", nPtBins, ptBins ) ;
 252     fhDecWMissedPartnerGeom1 = new TH1D("fhDecWMissedPartnerGeom1","Decay photons with partner missed due geometry Fid. area. 1", nPtBins, ptBins ) ;
 253     fhDecWMissedPartnerGeom2 = new TH1D("fhDecWMissedPartnerGeom2","Decay photons with partner missed due geometry Fid. area. 2", nPtBins, ptBins ) ;
 254     fhDecWMissedPartnerGeom3 = new TH1D("fhDecWMissedPartnerGeom3","Decay photons with partner missed due geometry Fid. area. 3", nPtBins, ptBins ) ;
 255
 256     //MC tagging: Decay partners spectra
 257     fhPartnerMCReg      = new TH1D("fhPartnerMCReg","Spectrum of decay partners which should be registered (MC)", nPtBins, ptBins ) ;
 258     fhPartnerMissedEmin = new TH1D("fhPartnerMissedEmin","Spectrum of decay partners lost due to Emin cut", nPtBins, ptBins ) ;
 259     fhPartnerMissedConv = new TH1D("fhPartnerMissedConv","Spectrum of decay partners lost due to conversion", nPtBins, ptBins ) ;
 260     fhPartnerMissedGeo  = new TH1D("fhPartnerMissedGeo","Spectrum of decay partners lost due to acceptance", nPtBins, ptBins ) ;
 261
 262     //Tagging
 263     fhTaggedAll    = new TH1D("fhTaggedAll","Spectrum of all tagged photons", nPtBins, ptBins ) ;
 264     fhTaggedArea1  = new TH1D("fhTaggedArea1","Spectrum of all tagged photons Fid. area1", nPtBins, ptBins ) ;
 265     fhTaggedArea2  = new TH1D("fhTaggedArea2","Spectrum of all tagged photons Fid. area2", nPtBins, ptBins ) ;
 266     fhTaggedArea3  = new TH1D("fhTaggedArea3","Spectrum of all tagged photons Fid. area3", nPtBins, ptBins ) ;
 267     fhTaggedPID[0] = new TH1D("fhTaggedPID0","Spectrum of tagged photons for PID=0", nPtBins, ptBins ) ;
 268     fhTaggedPID[1] = new TH1D("fhTaggedPID1","Spectrum of tagged photons for PID=1", nPtBins, ptBins ) ;
 269     fhTaggedPID[2] = new TH1D("fhTaggedPID2","Spectrum of tagged photons for PID=2", nPtBins, ptBins ) ;
 270     fhTaggedPID[3] = new TH1D("fhTaggedPID3","Spectrum of tagged photons for PID=3", nPtBins, ptBins ) ;
 271     fhTaggedMult   = new TH1D("fhTaggedMult","Spectrum of multiply tagged photons", nPtBins, ptBins ) ;
 272
 273     //Tagging: use MC information if available
 274     fhTaggedMCTrue    = new TH1D("fhTaggedMCTrue","Spectrum of correctly tagged pi0 decay photons", nPtBins, ptBins ) ;
 275     fhMCMissedTagging = new TH1D("fhMCMissedTagging","Spectrum of pi0 decay photons missed tagging due to wrong pair mass", nPtBins, ptBins ) ;
 276     fhMCFakeTagged    = new TH1D("fhMCFakeTagged","Spectrum of photons wrongly tagged according to MC", nPtBins, ptBins ) ;
 277
 278     //Invariant mass distributions for fake corrections
 279     const Int_t nmass=500 ;
 280     Double_t masses[nmass+1] ;
 281     for(Int_t i=0;i<=nmass;i++)masses[i]=0.002*i ;
 282     fhInvMassReal[0] = new TH2D("hInvMassRealPID0","Two-photon inv. mass vs first photon pt",nmass,masses,nPtBins, ptBins ) ;
 283     fhInvMassReal[1] = new TH2D("hInvMassRealPID1","Two-photon inv. mass vs first photon pt",nmass,masses,nPtBins, ptBins ) ;
 284     fhInvMassReal[2] = new TH2D("hInvMassRealPID2","Two-photon inv. mass vs first photon pt",nmass,masses,nPtBins, ptBins ) ;
 285     fhInvMassReal[3] = new TH2D("hInvMassRealPID3","Two-photon inv. mass vs first photon pt",nmass,masses,nPtBins, ptBins ) ;
 286     fhInvMassMixed[0]  = new TH2D("hInvMassMixedPID0","Two-photon inv. mass vs first photon pt",nmass,masses,nPtBins, ptBins ) ;
 287     fhInvMassMixed[1]  = new TH2D("hInvMassMixedPID1","Two-photon inv. mass vs first photon pt",nmass,masses,nPtBins, ptBins ) ;
 288     fhInvMassMixed[2]  = new TH2D("hInvMassMixedPID2","Two-photon inv. mass vs first photon pt",nmass,masses,nPtBins, ptBins ) ;
 289     fhInvMassMixed[3]  = new TH2D("hInvMassMixedPID3","Two-photon inv. mass vs first photon pt",nmass,masses,nPtBins, ptBins ) ;
 290     fhMCMissedTaggingMass= new TH2D("hMCMissedTaggingMass","Inv mass of pairs missed tagging",nmass,masses,nPtBins, ptBins ) ;
 291
 292     //Conversion and annihilation radius distributions
 293     fhConversionRadius  = new TH1D("fhConversionRadius","Radii of photon production (conversion)",100,0.,500.) ;
 294     fhInteractionRadius = new TH1D("fhInteractionRadius","Radii of photon production (hadron interaction)",100,0.,500.);
 295
 296     fhEvents               = new TH1D("hEvents", "Number of Events processed", 1, 0, 1);
 297
 298 // create output container
 299
 300   fOutputList = new TList() ;
 301   fEventList = new TList() ;
 302   fOutputList->SetName(GetName()) ;
 303
 304   fOutputList->AddAt(fhRecAll[0],                            0) ;
 305   fOutputList->AddAt(fhRecAll[1],                            1) ;
 306   fOutputList->AddAt(fhRecAll[2],                            2) ;
 307   fOutputList->AddAt(fhRecAll[3],                            3) ;
 308   fOutputList->AddAt(fhRecAllArea1,                          4) ;
 309   fOutputList->AddAt(fhRecAllArea2,                          5) ;
 310   fOutputList->AddAt(fhRecAllArea3,                          6) ;
 311
 312   fOutputList->AddAt(fhRecPhoton,                            7) ;
 313   fOutputList->AddAt(fhRecOther,                             8) ;
 314   fOutputList->AddAt(fhRecPhotonPID[0],                      9) ;
 315   fOutputList->AddAt(fhRecPhotonPID[1],                     10) ;
 316   fOutputList->AddAt(fhRecPhotonPID[2],                     11) ;
 317   fOutputList->AddAt(fhRecPhotonPID[3],                     12) ;
 318   fOutputList->AddAt(fhRecOtherPID[0],                      13) ;
 319   fOutputList->AddAt(fhRecOtherPID[1],                      14) ;
 320   fOutputList->AddAt(fhRecOtherPID[2],                      15) ;
 321   fOutputList->AddAt(fhRecOtherPID[3],                      16) ;
 322   fOutputList->AddAt(fhRecPhotPi0,                          17) ;
 323   fOutputList->AddAt(fhRecPhotEta,                          18) ;
 324   fOutputList->AddAt(fhRecPhotOmega,                        19) ;
 325   fOutputList->AddAt(fhRecPhotEtapr,                        20) ;
 326   fOutputList->AddAt(fhRecPhotConv,                         21) ;
 327   fOutputList->AddAt(fhRecPhotHadron,                       22) ;
 328   fOutputList->AddAt(fhRecPhotDirect,                       23) ;
 329   fOutputList->AddAt(fhRecPhotOther,                        24) ;
 330
 331   fOutputList->AddAt(fhDecWMCPartner,                       25) ;
 332   fOutputList->AddAt(fhDecWMissedPartnerNotPhoton,          26) ;
 333   fOutputList->AddAt(fhDecWMissedPartnerAll,                27) ;
 334   fOutputList->AddAt(fhDecWMissedPartnerEmin,               28) ;
 335   fOutputList->AddAt(fhDecWMissedPartnerConv,               29) ;
 336   fOutputList->AddAt(fhDecWMissedPartnerStack,              30) ;
 337   fOutputList->AddAt(fhDecWMissedPartnerGeom0,              31) ;
 338   fOutputList->AddAt(fhDecWMissedPartnerGeom1,              32) ;
 339   fOutputList->AddAt(fhDecWMissedPartnerGeom2,              33) ;
 340   fOutputList->AddAt(fhDecWMissedPartnerGeom3,              34) ;
 341
 342   fOutputList->AddAt(fhPartnerMCReg,                        35) ;
 343   fOutputList->AddAt(fhPartnerMissedEmin,                   36) ;
 344   fOutputList->AddAt(fhPartnerMissedConv,                   37) ;
 345   fOutputList->AddAt(fhPartnerMissedGeo,                    38) ;
 346
 347   fOutputList->AddAt(fhTaggedAll,                           39) ;
 348   fOutputList->AddAt(fhTaggedArea1,                         40) ;
 349   fOutputList->AddAt(fhTaggedArea2,                         41) ;
 350   fOutputList->AddAt(fhTaggedArea3,                         42) ;
 351   fOutputList->AddAt(fhTaggedPID[0],                        43) ;
 352   fOutputList->AddAt(fhTaggedPID[1],                        44) ;
 353   fOutputList->AddAt(fhTaggedPID[2],                        45) ;
 354   fOutputList->AddAt(fhTaggedPID[3],                        46) ;
 355   fOutputList->AddAt(fhTaggedMult,                          47) ;
 356
 357   fOutputList->AddAt(fhTaggedMCTrue,                        48) ;
 358   fOutputList->AddAt(fhMCMissedTagging,                     49) ;
 359   fOutputList->AddAt(fhMCFakeTagged,                        50) ;
 360
 361   fOutputList->AddAt(fhInvMassReal[0],                      51) ;
 362   fOutputList->AddAt(fhInvMassReal[1],                      52) ;
 363   fOutputList->AddAt(fhInvMassReal[2],                      53) ;
 364   fOutputList->AddAt(fhInvMassReal[3],                      54) ;
 365   fOutputList->AddAt(fhInvMassMixed[0],                     55) ;
 366   fOutputList->AddAt(fhInvMassMixed[1],                     56) ;
 367   fOutputList->AddAt(fhInvMassMixed[2],                     57) ;
 368   fOutputList->AddAt(fhInvMassMixed[3],                     58) ;
 369   fOutputList->AddAt(fhMCMissedTaggingMass,                 59) ;
 370
 371   fOutputList->AddAt(fhConversionRadius,                    60) ;
 372   fOutputList->AddAt(fhInteractionRadius,                   61) ;
 373
 374   fOutputList->AddAt(fhEvents,                              62) ;
 375
 376 }
 377
 378 //______________________________________________________________________________
 379 void AliAnalysisTaskTaggedPhotons::UserExec(Option_t *)
 380 {
 381   //Fill all histograms
 382
 383
 384   // We prepared 4 PID histograms, choose which PID combinations to use
 385   // see AliAODPWG4Particle for more choises
 386   //    0=all,2=Disp,4=Ch,6=Ch+Disp
 387   const Int_t pidMap[4]={0,2,4,6} ;
 388
 389   // Processing of one event
 390   if(fDebug>1)
 391     AliInfo(Form("\n\n Processing event # %lld",  Entry())) ;
 392   AliESDEvent* esd = (AliESDEvent*)InputEvent();
 393
 394
 395   //read geometry if not read yet
 396   if((fPHOS && fPHOSgeom==0) || (!fPHOS && fEMCALgeom==0))
 397     InitGeometry() ;
 398   if(fPHOS && fPHOSgeom!=0){
 399   }
 400   if(!fPHOS && fEMCALgeom!=0){ //EMCAL geometry initialized, but still need to set matrixes
 401     AliAODEvent * aod = 0x0 ;
 402     if(!esd)
 403       aod=dynamic_cast<AliAODEvent*>(InputEvent()) ;
 404     if(!esd && !aod)
 405       AliFatal("Can not read geometry matrixes from ESD/AOD: NO ESD") ;
 406     for(Int_t mod=0; mod < (fEMCALgeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
 407       if(esd){
 408         const TGeoHMatrix* m=esd->GetEMCALMatrix(mod) ;
 409         fEMCALgeom->SetMisalMatrix(m, mod) ;
 410       }
 411       else{
 412         const TGeoHMatrix* m=aod->GetHeader()->GetEMCALMatrix(mod) ;
 413         fEMCALgeom->SetMisalMatrix(m, mod) ;
 414       }
 415     }
 416   }
 417   //MC stack init
 418   fStack=0x0 ;
 419   if(AliAnalysisManager::GetAnalysisManager()){
 420     AliMCEventHandler* mctruth = (AliMCEventHandler*)((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler());
 421     if(mctruth)
 422       fStack = mctruth->MCEvent()->Stack();
 423   }
 424
 425   if(!fStack && gDebug>1)
 426     AliInfo("No stack! \n");
 427
 428
 429 //  Here one can choose trigger. But according to framework it should be chosen
 430 //  by external filters???
 431 //  TString trigClasses = esd->GetFiredTriggerClasses();
 432 //  if (!fStack && !trigClasses.Contains("CINT1B-ABCE-NOPF-ALL") ){
 433 //    Printf("Skip event with trigger class \"%s\"",trigClasses.Data());
 434 //    return;
 435 //  }
 436
 437   fhEvents->Fill(0.);
 438
 439   //************************  PHOS/EMCAL *************************************
 440   TRefArray * caloClustersArr  = new TRefArray();
 441   if(fPHOS)
 442     esd->GetPHOSClusters(caloClustersArr);
 443   else
 444     esd->GetEMCALClusters(caloClustersArr);
 445   const Int_t kNumberOfClusters = caloClustersArr->GetEntries() ;
 446
 447   TClonesArray * fCaloPhotonsArr   = new TClonesArray("AliAODPWG4Particle",kNumberOfClusters);
 448   Int_t inList = 0; //counter of caloClusters
 449
 450   Int_t cluster ;
 451
 452   // loop over Clusters
 453   for(cluster = 0 ; cluster < kNumberOfClusters ; cluster++) {
 454     AliESDCaloCluster * caloCluster = static_cast<AliESDCaloCluster*>(caloClustersArr->At(cluster)) ;
 455
 456     if((fPHOS && !caloCluster->IsPHOS()) ||
 457        (!fPHOS && caloCluster->IsPHOS()))
 458       continue ;
 459
 460     Double_t v[3] ; //vertex ;
 461     esd->GetVertex()->GetXYZ(v) ;
 462
 463     TLorentzVector momentum ;
 464     caloCluster->GetMomentum(momentum, v);
 465     new ((*fCaloPhotonsArr)[inList]) AliAODPWG4Particle(momentum.Px(),momentum.Py(),momentum.Pz(),caloCluster->E() );
 466     AliAODPWG4Particle *p = static_cast<AliAODPWG4Particle*>(fCaloPhotonsArr->At(inList));
 467     inList++;
 468
 469     p->SetCaloLabel(cluster,-1); //This and partner cluster
 470     p->SetDistToBad(Int_t(caloCluster->GetDistanceToBadChannel()));
 471
 472     p->SetTag(AliMCAnalysisUtils::kMCUnknown);
 473     p->SetTagged(kFALSE);   //Reconstructed pairs found
 474     p->SetLabel(caloCluster->GetLabel());
 475     Float_t pos[3] ;
 476     caloCluster->GetPosition(pos) ;
 477     p->SetFiducialArea(GetFiducialArea(pos)) ;
 478
 479     //PID criteria
 480     p->SetDispBit(TestDisp(caloCluster->GetM02(),caloCluster->GetM20(),caloCluster->E())) ;
 481     p->SetTOFBit(TestTOF(caloCluster->GetTOF(),caloCluster->E())) ;
 482     p->SetChargedBit(TestCharged(caloCluster->GetEmcCpvDistance(),caloCluster->E())) ;
 483
 484     for(Int_t iPID=0; iPID<4; iPID++){
 485       if(p->IsPIDOK(pidMap[iPID],22))
 486         fhRecAll[iPID]->Fill( p->Pt() ) ; //All recontructed particles
 487     }
 488     Int_t iFidArea = p->GetFiducialArea();
 489     if(iFidArea>0){
 490       fhRecAllArea1->Fill(p->Pt() ) ;
 491       if(iFidArea>1){
 492         fhRecAllArea2->Fill(p->Pt() ) ;
 493         if(iFidArea>2){
 494           fhRecAllArea3->Fill(p->Pt() ) ;
 495         }
 496       }
 497     }
 498
 499     if(fStack){
 500       TParticle * prim = fStack->Particle(caloCluster->GetLabel()) ;
 501       if(fDebug>2)
 502         printf("Pdgcode = %d\n",prim->GetPdgCode());
 503
 504       if(prim->GetPdgCode()!=22){ //not photon
 505 //<--DP        p->SetPhoton(kFALSE);
 506         fhRecOther->Fill(p->Pt()); //not photons real spectra
 507         for(Int_t iPID=0; iPID<4; iPID++){
 508           if(p->IsPIDOK(pidMap[iPID],22))
 509             fhRecOtherPID[iPID]->Fill(p->Pt());
 510         }
 511       }
 512       else{ //Primary photon (as in MC)
 513 //<--DP        p->SetPhoton(kTRUE);
 514         fhRecPhoton->Fill(p->Pt()); //Reconstructed with primary photon
 515         for(Int_t iPID=0; iPID<4; iPID++){
 516           if(p->IsPIDOK(pidMap[iPID],22))
 517             fhRecPhotonPID[iPID]->Fill(p->Pt());
 518         }
 519         Int_t pi0i=prim->GetFirstMother();
 520         Int_t grandpaPDG=-1 ;
 521         TParticle * pi0p = 0;
 522         if(pi0i>=0){
 523           pi0p=fStack->Particle(pi0i);
 524           grandpaPDG=pi0p->GetPdgCode() ;
 525         }
 526         switch(grandpaPDG){
 527         case 111: //Pi0 decay
 528                 //Primary decay photon (as in MC)
 529                 fhRecPhotPi0->Fill(p->Pt());
 530                 break ;
 531         case  11:
 532         case -11: //electron/positron conversion
 533 //<--DP                p->SetConverted(1);
 534                 fhRecPhotConv->Fill(p->Pt());  //Reconstructed with photon from conversion primary
 535                 fhConversionRadius->Fill(prim->R());
 536                 break ;
 537         case -2212:
 538         case -2112: //antineutron & antiproton conversion
 539                 fhRecPhotHadron->Fill(p->Pt());  //Reconstructed with photon from antibaryon annihilation
 540                 fhInteractionRadius->Fill(prim->R());
 541                 break ;
 542
 543         case 221: //eta decay
 544                 fhRecPhotEta->Fill(p->Pt());
 545                 break ;
 546
 547         case 223: //omega meson decay
 548                 fhRecPhotOmega->Fill(p->Pt());
 549                 break ;
 550
 551         case 331: //eta' decay
 552                 fhRecPhotEtapr->Fill(p->Pt());
 553                 break ;
 554
 555         case -1: //direct photon or no primary
 556                 fhRecPhotDirect->Fill(p->Pt());
 557                 break ;
 558
 559         default:
 560                 fhRecPhotOther->Fill(p->Pt());
 561                 break ;
 562         }
 563
 564         //Now classify pi0 decay photon
 565         if(grandpaPDG==111){
 566 //<--DP          p->Pi0Decay(kTRUE); //Mark this photon as primary decayed
 567 //<--DP          p->Pi0Id(pi0i); //remember id of the parent pi0
 568
 569           //Now check if second (partner) photon from pi0 decay hits PHOS or not
 570           //i.e. both photons can be tagged or it's the systematic error
 571 //<--DP          p->SetPartnerPt(0.);
 572           Int_t indexdecay1,indexdecay2;
 573
 574           indexdecay1=pi0p->GetFirstDaughter();
 575           indexdecay2=pi0p->GetLastDaughter();
 576           Int_t indexdecay=-1;
 577           if(fDebug>2)
 578                 printf("checking pi0 decay...index1=%d, index2=%d, index_pi0=%d, index_ph_prim=%d\n", indexdecay1,indexdecay2,pi0i,caloCluster->GetLabel());
 579
 580           if(indexdecay1!=caloCluster->GetLabel())
 581             indexdecay=indexdecay1;
 582           if(indexdecay2!=caloCluster->GetLabel())
 583             indexdecay=indexdecay2;
 584           if(indexdecay==-1){
 585             if(fDebug>2){
 586               printf("Probably the other photon is not in the stack!\n");
 587               printf("Number of daughters: %d\n",pi0p->GetNDaughters());
 588             }
 589             fhDecWMissedPartnerStack->Fill(p->Pt()) ;
 590           }
 591           else{
 592             TParticle *partner = fStack->Particle(indexdecay);
 593 //<--DP            p->SetPartnerPt(partner->Pt());
 594             if(partner->GetPdgCode()==22){
 595               Bool_t isPartnerLost=kFALSE; //If partner is lost for some reason
 596               if(partner->GetNDaughters()!=0){ //this photon is converted before it is registered by some detector
 597                 if(fDebug>2)
 598                   printf("P_Conv, daughters=%d\n",partner->GetNDaughters());
 599 //<--DP                p->SetConvertedPartner(1);
 600                 fhPartnerMissedConv->Fill(partner->Pt());
 601                 fhDecWMissedPartnerConv->Fill(p->Pt()) ;  //Spectrum of tagged with missed partner
 602                 isPartnerLost=kTRUE;
 603               }
 604               Bool_t impact = kFALSE ;
 605               if(fPHOS){
 606                 Int_t modulenum ;
 607                 Double_t ztmp,xtmp ;
 608                 impact=fPHOSgeom->ImpactOnEmc(partner,modulenum,ztmp,xtmp) ;
 609                 if(fDebug>2){
 610                   printf("Impact on PHOS: module: %d, x tower: %f, z tower: %f\n", modulenum,xtmp,ztmp);
 611                 }
 612               }
 613               else{
 614                 impact = fEMCALgeom->Impact(partner) ;
 615               }
 616               if(!impact){ //this photon cannot hit PHOS
 617                 if(fDebug>2)
 618                   printf("P_Geo\n");
 619                 fhPartnerMissedGeo->Fill(partner->Pt());
 620                 fhDecWMissedPartnerGeom0->Fill(p->Pt()) ;  //Spectrum of tagged with missed partner
 621                 if(iFidArea>0){
 622                   fhDecWMissedPartnerGeom1->Fill(p->Pt()) ;  //Spectrum of tagged with missed partner
 623                   if(iFidArea>1){
 624                     fhDecWMissedPartnerGeom2->Fill(p->Pt()) ;  //Spectrum of tagged with missed partner
 625                     if(iFidArea>2){
 626                       fhDecWMissedPartnerGeom3->Fill(p->Pt()) ;  //Spectrum of tagged with missed partner
 627                     }
 628                   }
 629                 }
 630                 isPartnerLost=kTRUE;
 631               }
 632               if(!isPartnerLost && partner->Energy()<fMinEnergyCut){ //energy is not enough to be registered by PHOS
 633                 if(fDebug>2)
 634                   printf("P_Reg, E=%f\n",partner->Energy());
 635                 fhPartnerMissedEmin->Fill(partner->Pt());  //Spectrum of missed partners
 636                 fhDecWMissedPartnerEmin->Fill(p->Pt()) ;  //Spectrum of tagged with missed partner
 637                 isPartnerLost=kTRUE;
 638               }
 639               if(!isPartnerLost){
 640 //                p->SetMCTagged(1); //set this photon as primary tagged
 641                 fhDecWMCPartner->Fill(p->Pt());
 642                 fhPartnerMCReg->Fill(partner->Pt());
 643                 if(fDebug>2){
 644                   printf("both photons are inside PHOS. Energy: %f, Pt of pair photon: %f, E of pair photon: %f, Px: %f Py: %f Pz: %f, num of daughters: %d \n", caloCluster->E(),partner->Pt(),partner->Energy(),partner->Px(),partner->Py(),partner->Pz(),partner->GetNDaughters());
 645                 }
 646               }
 647               else{
 648                 fhDecWMissedPartnerAll->Fill(p->Pt());
 649               }
 650             }//Partner - photon
 651             else{//partner not photon
 652               fhDecWMissedPartnerNotPhoton->Fill(p->Pt());
 653             }
 654           }
 655         }
 656       }
 657     }
 658   } //PHOS/EMCAL clusters
 659
 660   if(fDebug>1)
 661     printf("number of clusters: %d\n",inList);
 662
 663   //Invariant Mass analysis
 664   for(Int_t phosPhoton1 = 0 ; phosPhoton1 < inList ; phosPhoton1++) {
 665     AliAODPWG4Particle * p1 = static_cast<AliAODPWG4Particle*>(fCaloPhotonsArr->At(phosPhoton1));
 666     for(Int_t phosPhoton2 = 0 ; phosPhoton2 < inList ; phosPhoton2++) {
 667       if(phosPhoton1==phosPhoton2)
 668         continue ;
 669       AliAODPWG4Particle * p2 = static_cast<AliAODPWG4Particle*>(fCaloPhotonsArr->At(phosPhoton2));
 670
 671       Double_t invMass = p1->GetPairMass(p2);
 672       for(Int_t iPID=0; iPID<4; iPID++){
 673         if(p1->IsPIDOK(pidMap[iPID],22))
 674           fhInvMassReal[iPID]->Fill(invMass,p1->Pt());
 675         if(p2->IsPIDOK(pidMap[iPID],22))
 676           fhInvMassReal[iPID]->Fill(invMass,p2->Pt());
 677       }
 678       if(fDebug>2)
 679           printf("Pair i=%d,j=%d, M=%f\n",phosPhoton1,phosPhoton2,invMass);
 680
 681       Bool_t makePi0=IsInPi0Band(invMass,p1->Pt());
 682
 683       if(makePi0 && p1->IsTagged()){//Multiple tagging
 684         fhTaggedMult->Fill(p1->Pt());
 685       }
 686       if(makePi0 && !p1->IsTagged()){//Each photon should enter histogram once, even if tagged several times
 687         fhTaggedAll->Fill(p1->Pt());
 688         Int_t iFidArea = p1->GetFiducialArea();
 689         if(iFidArea>0){
 690           fhTaggedArea1->Fill(p1->Pt() ) ;
 691           if(iFidArea>1){
 692             fhTaggedArea2->Fill(p1->Pt() ) ;
 693             if(iFidArea>2){
 694               fhTaggedArea3->Fill(p1->Pt() ) ;
 695             }
 696           }
 697         }
 698
 699         for(Int_t iPID=0; iPID<4; iPID++){
 700           if(p1->IsPIDOK(pidMap[iPID],22))
 701             fhTaggedPID[iPID]->Fill(p1->Pt());
 702         }
 703
 704         p1->SetTagged(kTRUE) ;
 705       }
 706
 707       //First chesk if this is true pi0 pair
 708       if(IsSamePi0(p1,p2)){ //Correctly tagged - from the same pi0
 709 //        p1->SetTrueTagged(1);
 710 //        p2->SetTrueTagged(1);
 711         if(makePi0)//Correctly tagged photons
 712           fhTaggedMCTrue->Fill(p1->Pt());
 713         else{ //Decay pair missed tagging
 714           fhMCMissedTagging->Fill(p1->Pt());
 715           fhMCMissedTaggingMass->Fill(invMass,p1->Pt()) ;
 716           //Clussify why missed tagging (todo)
 717           //Converted
 718           //Partner not a photon
 719           //Tagged not a photon
 720           //Just wrong inv.mass
 721         }
 722       }
 723       else{//Fake tagged - not from the same pi0
 724         if(makePi0)//Fake pair
 725           fhMCFakeTagged->Fill(p1->Pt());
 726       }
 727     }
 728   }
 729
 730   //Fill Mixed InvMass distributions:
 731   TIter nextEv(fEventList) ;
 732   while(TClonesArray * event2 = static_cast<TClonesArray*>(nextEv())){
 733     Int_t nPhotons2 = event2->GetEntriesFast() ;
 734     for(Int_t i=0; i < inList ; i++){
 735       AliAODPWG4Particle * p1 = static_cast<AliAODPWG4Particle*>(fCaloPhotonsArr->At(i)) ;
 736       for(Int_t j=0; j < nPhotons2 ; j++){
 737         AliAODPWG4Particle * p2 = static_cast<AliAODPWG4Particle*>(event2->At(j)) ;
 738         Double_t invMass = p1->GetPairMass(p2);
 739         for(Int_t iPID=0; iPID<4; iPID++){
 740           if(p1->IsPIDOK(pidMap[iPID],22))
 741             fhInvMassMixed[iPID]->Fill(invMass,p1->Pt());
 742           if(p2->IsPIDOK(pidMap[iPID],22))
 743             fhInvMassMixed[iPID]->Fill(invMass,p2->Pt());
 744         }
 745       }
 746     }
 747   }
 748
 749   //Remove old events
 750   fEventList->AddFirst(fCaloPhotonsArr);
 751   if(fEventList->GetSize() > 10){
 752     TClonesArray *tmp = static_cast <TClonesArray*> (fEventList->Last());
 753     fEventList->Remove(tmp);
 754     delete tmp;
 755   }
 756
 757   delete caloClustersArr;
 758
 759   PostData(1, fOutputList);
 760 }
 761
 762
 763 //______________________________________________________________________________
 764 void AliAnalysisTaskTaggedPhotons::Init()
 765 {
 766   // Intialisation of parameters
 767   AliInfo("Doing initialisation") ;
 768   SetPhotonId(0.9) ;
 769   SetMinEnergyCut(0.4);
 770   SetPi0MeanParameters(0.136,0.,0.0,0.0);
 771 //  SetPi0MeanParameters(0.1377,-0.002566,0.001216,-0.0001256);
 772   SetPi0SigmaParameters(0.004508,0.005497,0.00000006);
 773 }
 774
 775 //______________________________________________________________________________
 776 void AliAnalysisTaskTaggedPhotons::Terminate(Option_t *)
 777 {
 778   // Processing when the event loop is ended
 779   if (fDebug > 1) Printf("Terminate()");
 780 }
 781 //______________________________________________________________________________
 782 Bool_t AliAnalysisTaskTaggedPhotons::IsInPi0Band(Double_t m, Double_t pt)const
 783 {
 784   //Parameterization of the pi0 peak region
 785   Double_t mpi0mean = fPi0MeanP0 + fPi0MeanP1 * pt + fPi0MeanP2 * pt*pt + fPi0MeanP3 * pt*pt*pt;
 786   Double_t mpi0sigma = TMath::Sqrt(fPi0SigmaP0 * fPi0SigmaP0 / pt + fPi0SigmaP1 * fPi0SigmaP1 + fPi0SigmaP2 * fPi0SigmaP2 / pt / pt);
 787
 788   return (m>mpi0mean-2*mpi0sigma && m<mpi0mean+2*mpi0sigma) ;
 789 }
 790 //______________________________________________________________________________
 791 Bool_t AliAnalysisTaskTaggedPhotons::IsSamePi0(const AliAODPWG4Particle *p1, const AliAODPWG4Particle *p2)const{
 792   //Looks through parents and finds if there was commont pi0 among ancestors
 793
 794   if(!fStack)
 795     return kFALSE ; //can not say anything
 796
 797   Int_t prim1 = p1->GetLabel();
 798   while(prim1!=-1){
 799     Int_t prim2 = p2->GetLabel();
 800     while(prim2!=-1){
 801       if(prim1==prim2){
 802         if(fStack->Particle(prim1)->GetPdgCode()==111)
 803           return kTRUE ;
 804         else
 805           return kFALSE ;
 806       }
 807       prim2=fStack->Particle(prim2)->GetFirstMother() ;
 808     }
 809     prim1=fStack->Particle(prim1)->GetFirstMother() ;
 810   }
 811   return kFALSE ;
 812 }
 813 //______________________________________________________________________________
 814 Int_t AliAnalysisTaskTaggedPhotons::GetFiducialArea(Float_t * pos)const{
 815   //calculates in which kind of fiducial area photon hit
 816   Double_t phi=TMath::ATan2(pos[1],pos[0]) ;
 817   Double_t z=pos[2] ;
 818   while(phi>TMath::TwoPi())phi-=TMath::TwoPi() ;
 819   while(phi<0.)phi+=TMath::TwoPi() ;
 820   if(fDebug>2){
 821     printf("FiducialArea: phi=%f, z=%f \n",phi,z) ;
 822         printf("   fZmax=%f, fZmin=%f, fPhimax=%f, fPhimin=%f \n",fZmax,fZmin,fPhimax,fPhimin) ;
 823   }
 824   if(fPHOS){
 825     //From active PHOS area remove bands in 10 cm
 826     const Double_t kphi=TMath::ATan(10./460.) ; //angular band width
 827     Double_t dzMax=TMath::Ceil((fZmax-z)/10.) ;
 828     Double_t dzMin=TMath::Ceil((z-fZmin)/10.) ;
 829     Double_t dphiMax=TMath::Ceil((fPhimax-phi)/kphi);
 830     Double_t dphiMin=TMath::Ceil((phi-fPhimin)/kphi);
 831     if(fDebug>2){
 832      printf("In PHOS \n") ;
 833     printf("    dzMax=%f, dzMin=%f, dphiMax=%f, dphiMin=%f ret=%d\n",dzMax,dzMin,dphiMax,dphiMin,(Int_t)TMath::Min(TMath::Min(dzMax,dzMin),TMath::Min(dphiMax,dphiMin))) ;
 834     }
 835     return (Int_t)TMath::Min(TMath::Min(dzMax,dzMin),TMath::Min(dphiMax,dphiMin));
 836   }
 837   else{//EMCAL
 838     //From active EMCAL area remove bands in 20 cm
 839     const Double_t kphi=TMath::ATan(20./428.) ; //angular band width
 840     Double_t dzMax=TMath::Ceil((fZmax-z)/20.) ;
 841     Double_t dzMin=TMath::Ceil((z-fZmin)/20.) ;
 842     Double_t dphiMax=TMath::Ceil((fPhimax-phi)/kphi);
 843     Double_t dphiMin=TMath::Ceil((phi-fPhimin)/kphi);
 844     return (Int_t)TMath::Min(TMath::Min(dzMax,dzMin),TMath::Min(dphiMax,dphiMin));
 845   }
 846 }
 847 //______________________________________________________________________________
 848 Bool_t  AliAnalysisTaskTaggedPhotons::TestDisp(Double_t l0, Double_t l1, Double_t /*e*/)const{
 849   //test if dispersion corresponds to those of photon
 850   if(fPHOS){
 851 //    Double_t l0mean=1.38736+0.490405*TMath::Exp(-e*0.286170) ;
 852 //    Double_t l1mean=1.09786-0.323469*TMath::Exp(-e*0.918719) ;
 853 //    Double_t l0sigma=0.159905+0.829831/e-0.158067/e/e ;
 854 //    Double_t l1sigma=0.133170+0.404387/e-0.0426302/e/e ;
 855 //    Double_t c =-0.382233 ;
 856 //    return ((l0-l0mean)*(l0-l0mean)/l0sigma/l0sigma + (l1-l1mean)*(l1-l1mean)/l1sigma/l1sigma+c*(l0-l0mean)*(l1-l1mean)/l0sigma/l1sigma)<1. ;
 857
 858 if(l1>= 0   && l0>= 0   && l1 < 0.1 && l0 < 0.1) return kFALSE;
 859 if(l1>= 0   && l0 > 0.5 && l1 < 0.1 && l0 < 1.5) return kTRUE;
 860 if(l1>= 0   && l0 > 2.0 && l1 < 0.1 && l0 < 2.7) return kFALSE;
 861 if(l1>= 0   && l0 > 2.7 && l1 < 0.1 && l0 < 4.0) return kFALSE;
 862 if(l1 > 0.1 && l1 < 0.7 && l0 > 0.7 && l0 < 2.1) return kTRUE;
 863 if(l1 > 0.1 && l1 < 0.3 && l0 > 3.0 && l0 < 5.0) return kFALSE;
 864 if(l1 > 0.3 && l1 < 0.7 && l0 > 2.5 && l0 < 4.0) return kFALSE;
 865 if(l1 > 0.7 && l1 < 1.3 && l0 > 1.0 && l0 < 1.6) return kTRUE;
 866 if(l1 > 0.7 && l1 < 1.3 && l0 > 1.6 && l0 < 3.5) return kTRUE;
 867 if(l1 > 1.3 && l1 < 3.5 && l0 > 1.3 && l0 < 3.5) return kTRUE;
 868    return kFALSE ;
 869
 870   }
 871   else{ //EMCAL: not ready yet
 872    return kTRUE ;
 873
 874   }
 875
 876 }
 877 //______________________________________________________________________________^M
 878 Bool_t  AliAnalysisTaskTaggedPhotons::TestCharged(Double_t dr,Double_t /*en*/)const{
 879
 880   if(dr<15.) return kFALSE ;
 881   return kTRUE ;
 882
 883 }
 884 //______________________________________________________________________________^M
 885 void  AliAnalysisTaskTaggedPhotons::InitGeometry(){
 886   //Read rotation matrixes from ESD
 887
 888   if(fDebug>1)printf("Init geometry \n") ;
 889   AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent()) ;
 890   AliAODEvent * aod = 0x0 ;
 891   if(!esd)
 892     aod=dynamic_cast<AliAODEvent*>(InputEvent()) ;
 893   if(!esd && !aod)
 894     AliFatal("Can not read geometry matrixes from ESD/AOD: NO ESD") ;
 895   if(fPHOS){//reading PHOS matrixes
 896     fPHOSgeom = new AliPHOSGeoUtils("IHEP","");
 897     Bool_t activeMod[5]={0} ;
 898     for(Int_t mod=0; mod<5; mod++){
 899       if(esd){
 900         const TGeoHMatrix* m=esd->GetPHOSMatrix(mod) ;
 901         fPHOSgeom->SetMisalMatrix(m, mod) ;
 902         if(m!=0)
 903          activeMod[mod]=kTRUE ;
 904       }
 905       else{
 906         const TGeoHMatrix* m=aod->GetHeader()->GetPHOSMatrix(mod) ;
 907         fPHOSgeom->SetMisalMatrix(m, mod) ;
 908         if(m!=0)
 909          activeMod[mod]=kTRUE ;
 910       }
 911     }
 912
 913     if(fZmax>fZmin) //already  set manually
 914       return ;
 915     fZmax= 999. ;
 916     fZmin=-999. ;
 917     fPhimax=-999. ;
 918     fPhimin= 999. ;
 919     for(Int_t imod=1; imod<=5; imod++){
 920       if(!activeMod[imod-1])
 921         continue ;
 922     //Find exact coordinates of PHOS corners
 923       Int_t relId[4]={imod,0,1,1} ;
 924       Float_t x,z ;
 925       fPHOSgeom->RelPosInModule(relId,x,z) ;
 926       TVector3 pos ;
 927       fPHOSgeom->Local2Global(imod,x,z,pos) ;
 928       Double_t phi=pos.Phi() ;
 929       while(phi<0.)phi+=TMath::TwoPi() ;
 930       while(phi>TMath::TwoPi())phi-=TMath::TwoPi() ;
 931       fPhimin=TMath::Min(fPhimin,float(phi)) ;
 932       fZmin=TMath::Max(fZmin,float(pos.Z())) ;
 933       relId[2]=64 ;
 934       relId[3]=56 ;
 935       fPHOSgeom->RelPosInModule(relId,x,z) ;
 936       fPHOSgeom->Local2Global(imod,x,z,pos) ;
 937       phi=pos.Phi() ;
 938       while(phi<0.)phi+=TMath::TwoPi() ;
 939       while(phi>TMath::TwoPi())phi-=TMath::TwoPi() ;
 940       fPhimax=TMath::Max(fPhimax,float(phi)) ;
 941       fZmax=TMath::Min(fZmax,float(pos.Z())) ;
 942     }
 943   }
 944   else{ //EMCAL
 945     fEMCALgeom = new AliEMCALGeoUtils("EMCAL_COMPLETE");
 946     for(Int_t mod=0; mod < (fEMCALgeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){ //<---Gustavo, could you check???
 947       if(esd){
 948         const TGeoHMatrix* m=esd->GetEMCALMatrix(mod) ;
 949         fEMCALgeom->SetMisalMatrix(m, mod) ;
 950       }
 951       else{
 952         const TGeoHMatrix* m=aod->GetHeader()->GetEMCALMatrix(mod) ;
 953         fEMCALgeom->SetMisalMatrix(m, mod) ;
 954       }
 955     }
 956     if(fZmax>fZmin) //already  set manually
 957       return ;
 958
 959     fZmax= 999. ;
 960     fZmin=-999. ;
 961     fPhimax=-999. ;
 962     fPhimin= 999. ;
 963     for(Int_t imod=0; imod<12; imod++){
 964     //Find exact coordinates of SM corners
 965       Int_t absId = fEMCALgeom->GetAbsCellIdFromCellIndexes(imod, 0, 0);
 966       TVector3 pos ;
 967       //Get the position of this tower.
 968       fEMCALgeom->RelPosCellInSModule(absId,pos);
 969       Double_t phi=pos.Phi() ;
 970       while(phi<0.)phi+=TMath::TwoPi() ;
 971       while(phi>TMath::TwoPi())phi-=TMath::TwoPi() ;
 972       fPhimin=TMath::Min(fPhimin,float(phi)) ;
 973       fZmin=TMath::Max(fZmin,float(pos.Z())) ;
 974       absId = fEMCALgeom->GetAbsCellIdFromCellIndexes(imod, 24, 48);
 975       fEMCALgeom->RelPosCellInSModule(absId,pos);
 976       phi=pos.Phi() ;
 977       while(phi<0.)phi+=TMath::TwoPi() ;
 978       while(phi>TMath::TwoPi())phi-=TMath::TwoPi() ;
 979       fPhimax=TMath::Max(fPhimax,float(phi)) ;
 980       fZmax=TMath::Min(fZmax,float(pos.Z())) ;
 981     }
 982   }
 983 }