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
 405     if(!esd && !aod){
 406       AliFatal("Can not read geometry matrixes from ESD/AOD: NO ESD") ;
 407     }
 408     else{
 409
 410
 411       for(Int_t mod=0; mod < (fEMCALgeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
 412         if(esd){
 413           const TGeoHMatrix* m=esd->GetEMCALMatrix(mod) ;
 414           fEMCALgeom->SetMisalMatrix(m, mod) ;
 415         }
 416         else{
 417           const TGeoHMatrix* m=aod->GetHeader()->GetEMCALMatrix(mod) ;
 418           fEMCALgeom->SetMisalMatrix(m, mod) ;
 419         }
 420       }
 421     }
 422     //MC stack init
 423     fStack=0x0 ;
 424     if(AliAnalysisManager::GetAnalysisManager()){
 425       AliMCEventHandler* mctruth = (AliMCEventHandler*)((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler());
 426       if(mctruth)
 427         fStack = mctruth->MCEvent()->Stack();
 428     }
 429
 430     if(!fStack && gDebug>1)
 431       AliInfo("No stack! \n");
 432
 433
 434     //  Here one can choose trigger. But according to framework it should be chosen
 435     //  by external filters???
 436     //  TString trigClasses = esd->GetFiredTriggerClasses();
 437     //  if (!fStack && !trigClasses.Contains("CINT1B-ABCE-NOPF-ALL") ){
 438     //    Printf("Skip event with trigger class \"%s\"",trigClasses.Data());
 439     //    return;
 440     //  }
 441
 442     fhEvents->Fill(0.);
 443
 444     //************************  PHOS/EMCAL *************************************
 445     TRefArray * caloClustersArr  = new TRefArray();
 446     if(fPHOS)
 447       esd->GetPHOSClusters(caloClustersArr);
 448     else
 449       esd->GetEMCALClusters(caloClustersArr);
 450     const Int_t kNumberOfClusters = caloClustersArr->GetEntries() ;
 451
 452     TClonesArray * fCaloPhotonsArr   = new TClonesArray("AliAODPWG4Particle",kNumberOfClusters);
 453     Int_t inList = 0; //counter of caloClusters
 454
 455     Int_t cluster ;
 456
 457     // loop over Clusters
 458     for(cluster = 0 ; cluster < kNumberOfClusters ; cluster++) {
 459       AliESDCaloCluster * caloCluster = static_cast<AliESDCaloCluster*>(caloClustersArr->At(cluster)) ;
 460
 461       if((fPHOS && !caloCluster->IsPHOS()) ||
 462          (!fPHOS && caloCluster->IsPHOS()))
 463         continue ;
 464
 465       Double_t v[3] ; //vertex ;
 466       esd->GetVertex()->GetXYZ(v) ;
 467
 468       TLorentzVector momentum ;
 469       caloCluster->GetMomentum(momentum, v);
 470       new ((*fCaloPhotonsArr)[inList]) AliAODPWG4Particle(momentum.Px(),momentum.Py(),momentum.Pz(),caloCluster->E() );
 471       AliAODPWG4Particle *p = static_cast<AliAODPWG4Particle*>(fCaloPhotonsArr->At(inList));
 472       inList++;
 473
 474       p->SetCaloLabel(cluster,-1); //This and partner cluster
 475       p->SetDistToBad(Int_t(caloCluster->GetDistanceToBadChannel()));
 476
 477       p->SetTag(AliMCAnalysisUtils::kMCUnknown);
 478       p->SetTagged(kFALSE);   //Reconstructed pairs found
 479       p->SetLabel(caloCluster->GetLabel());
 480       Float_t pos[3] ;
 481       caloCluster->GetPosition(pos) ;
 482       p->SetFiducialArea(GetFiducialArea(pos)) ;
 483
 484       //PID criteria
 485       p->SetDispBit(TestDisp(caloCluster->GetM02(),caloCluster->GetM20(),caloCluster->E())) ;
 486       p->SetTOFBit(TestTOF(caloCluster->GetTOF(),caloCluster->E())) ;
 487       p->SetChargedBit(TestCharged(caloCluster->GetEmcCpvDistance(),caloCluster->E())) ;
 488
 489       for(Int_t iPID=0; iPID<4; iPID++){
 490         if(p->IsPIDOK(pidMap[iPID],22))
 491           fhRecAll[iPID]->Fill( p->Pt() ) ; //All recontructed particles
 492       }
 493       Int_t iFidArea = p->GetFiducialArea();
 494       if(iFidArea>0){
 495         fhRecAllArea1->Fill(p->Pt() ) ;
 496         if(iFidArea>1){
 497           fhRecAllArea2->Fill(p->Pt() ) ;
 498           if(iFidArea>2){
 499             fhRecAllArea3->Fill(p->Pt() ) ;
 500           }
 501         }
 502       }
 503
 504       if(fStack){
 505         TParticle * prim = fStack->Particle(caloCluster->GetLabel()) ;
 506         if(fDebug>2)
 507           printf("Pdgcode = %d\n",prim->GetPdgCode());
 508
 509         if(prim->GetPdgCode()!=22){ //not photon
 510           //<--DP        p->SetPhoton(kFALSE);
 511           fhRecOther->Fill(p->Pt()); //not photons real spectra
 512           for(Int_t iPID=0; iPID<4; iPID++){
 513             if(p->IsPIDOK(pidMap[iPID],22))
 514               fhRecOtherPID[iPID]->Fill(p->Pt());
 515           }
 516         }
 517         else{ //Primary photon (as in MC)
 518           //<--DP        p->SetPhoton(kTRUE);
 519           fhRecPhoton->Fill(p->Pt()); //Reconstructed with primary photon
 520           for(Int_t iPID=0; iPID<4; iPID++){
 521             if(p->IsPIDOK(pidMap[iPID],22))
 522               fhRecPhotonPID[iPID]->Fill(p->Pt());
 523           }
 524           Int_t pi0i=prim->GetFirstMother();
 525           Int_t grandpaPDG=-1 ;
 526           TParticle * pi0p = 0;
 527           if(pi0i>=0){
 528             pi0p=fStack->Particle(pi0i);
 529             grandpaPDG=pi0p->GetPdgCode() ;
 530           }
 531           switch(grandpaPDG){
 532             case 111: //Pi0 decay
 533               //Primary decay photon (as in MC)
 534               fhRecPhotPi0->Fill(p->Pt());
 535               break ;
 536             case  11:
 537             case -11: //electron/positron conversion
 538               //<--DP                p->SetConverted(1);
 539               fhRecPhotConv->Fill(p->Pt());  //Reconstructed with photon from conversion primary
 540               fhConversionRadius->Fill(prim->R());
 541               break ;
 542             case -2212:
 543             case -2112: //antineutron & antiproton conversion
 544               fhRecPhotHadron->Fill(p->Pt());  //Reconstructed with photon from antibaryon annihilation
 545               fhInteractionRadius->Fill(prim->R());
 546               break ;
 547
 548             case 221: //eta decay
 549               fhRecPhotEta->Fill(p->Pt());
 550               break ;
 551
 552             case 223: //omega meson decay
 553               fhRecPhotOmega->Fill(p->Pt());
 554               break ;
 555
 556             case 331: //eta' decay
 557               fhRecPhotEtapr->Fill(p->Pt());
 558               break ;
 559
 560             case -1: //direct photon or no primary
 561               fhRecPhotDirect->Fill(p->Pt());
 562               break ;
 563
 564             default:
 565               fhRecPhotOther->Fill(p->Pt());
 566               break ;
 567           }
 568
 569           //Now classify pi0 decay photon
 570           if(grandpaPDG==111){
 571             //<--DP          p->Pi0Decay(kTRUE); //Mark this photon as primary decayed
 572             //<--DP          p->Pi0Id(pi0i); //remember id of the parent pi0
 573
 574             //Now check if second (partner) photon from pi0 decay hits PHOS or not
 575             //i.e. both photons can be tagged or it's the systematic error
 576             //<--DP          p->SetPartnerPt(0.);
 577             Int_t indexdecay1,indexdecay2;
 578
 579             indexdecay1=pi0p->GetFirstDaughter();
 580             indexdecay2=pi0p->GetLastDaughter();
 581             Int_t indexdecay=-1;
 582             if(fDebug>2)
 583               printf("checking pi0 decay...index1=%d, index2=%d, index_pi0=%d, index_ph_prim=%d\n", indexdecay1,indexdecay2,pi0i,caloCluster->GetLabel());
 584
 585             if(indexdecay1!=caloCluster->GetLabel())
 586               indexdecay=indexdecay1;
 587             if(indexdecay2!=caloCluster->GetLabel())
 588               indexdecay=indexdecay2;
 589             if(indexdecay==-1){
 590               if(fDebug>2){
 591                 printf("Probably the other photon is not in the stack!\n");
 592                 printf("Number of daughters: %d\n",pi0p->GetNDaughters());
 593               }
 594               fhDecWMissedPartnerStack->Fill(p->Pt()) ;
 595             }
 596             else{
 597               TParticle *partner = fStack->Particle(indexdecay);
 598               //<--DP            p->SetPartnerPt(partner->Pt());
 599               if(partner->GetPdgCode()==22){
 600                 Bool_t isPartnerLost=kFALSE; //If partner is lost for some reason
 601                 if(partner->GetNDaughters()!=0){ //this photon is converted before it is registered by some detector
 602                   if(fDebug>2)
 603                     printf("P_Conv, daughters=%d\n",partner->GetNDaughters());
 604                   //<--DP                p->SetConvertedPartner(1);
 605                   fhPartnerMissedConv->Fill(partner->Pt());
 606                   fhDecWMissedPartnerConv->Fill(p->Pt()) ;  //Spectrum of tagged with missed partner
 607                   isPartnerLost=kTRUE;
 608                 }
 609                 Bool_t impact = kFALSE ;
 610                 if(fPHOS){
 611                   Int_t modulenum ;
 612                   Double_t ztmp,xtmp ;
 613                   impact=fPHOSgeom->ImpactOnEmc(partner,modulenum,ztmp,xtmp) ;
 614                   if(fDebug>2){
 615                     printf("Impact on PHOS: module: %d, x tower: %f, z tower: %f\n", modulenum,xtmp,ztmp);
 616                   }
 617                 }
 618                 else{
 619                   impact = fEMCALgeom->Impact(partner) ;
 620                 }
 621                 if(!impact){ //this photon cannot hit PHOS
 622                   if(fDebug>2)
 623                     printf("P_Geo\n");
 624                   fhPartnerMissedGeo->Fill(partner->Pt());
 625                   fhDecWMissedPartnerGeom0->Fill(p->Pt()) ;  //Spectrum of tagged with missed partner
 626                   if(iFidArea>0){
 627                     fhDecWMissedPartnerGeom1->Fill(p->Pt()) ;  //Spectrum of tagged with missed partner
 628                     if(iFidArea>1){
 629                       fhDecWMissedPartnerGeom2->Fill(p->Pt()) ;  //Spectrum of tagged with missed partner
 630                       if(iFidArea>2){
 631                         fhDecWMissedPartnerGeom3->Fill(p->Pt()) ;  //Spectrum of tagged with missed partner
 632                       }
 633                     }
 634                   }
 635                   isPartnerLost=kTRUE;
 636                 }
 637                 if(!isPartnerLost && partner->Energy()<fMinEnergyCut){ //energy is not enough to be registered by PHOS
 638                   if(fDebug>2)
 639                     printf("P_Reg, E=%f\n",partner->Energy());
 640                   fhPartnerMissedEmin->Fill(partner->Pt());  //Spectrum of missed partners
 641                   fhDecWMissedPartnerEmin->Fill(p->Pt()) ;  //Spectrum of tagged with missed partner
 642                   isPartnerLost=kTRUE;
 643                 }
 644                 if(!isPartnerLost){
 645                   //                p->SetMCTagged(1); //set this photon as primary tagged
 646                   fhDecWMCPartner->Fill(p->Pt());
 647                   fhPartnerMCReg->Fill(partner->Pt());
 648                   if(fDebug>2){
 649                     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());
 650                   }
 651                 }
 652                 else{
 653                   fhDecWMissedPartnerAll->Fill(p->Pt());
 654                 }
 655               }//Partner - photon
 656               else{//partner not photon
 657                 fhDecWMissedPartnerNotPhoton->Fill(p->Pt());
 658               }
 659             }
 660           }
 661         }
 662       }
 663     } //PHOS/EMCAL clusters
 664
 665     if(fDebug>1)
 666       printf("number of clusters: %d\n",inList);
 667
 668     //Invariant Mass analysis
 669     for(Int_t phosPhoton1 = 0 ; phosPhoton1 < inList ; phosPhoton1++) {
 670       AliAODPWG4Particle * p1 = static_cast<AliAODPWG4Particle*>(fCaloPhotonsArr->At(phosPhoton1));
 671       for(Int_t phosPhoton2 = 0 ; phosPhoton2 < inList ; phosPhoton2++) {
 672         if(phosPhoton1==phosPhoton2)
 673           continue ;
 674         AliAODPWG4Particle * p2 = static_cast<AliAODPWG4Particle*>(fCaloPhotonsArr->At(phosPhoton2));
 675
 676         Double_t invMass = p1->GetPairMass(p2);
 677         for(Int_t iPID=0; iPID<4; iPID++){
 678           if(p1->IsPIDOK(pidMap[iPID],22))
 679             fhInvMassReal[iPID]->Fill(invMass,p1->Pt());
 680           if(p2->IsPIDOK(pidMap[iPID],22))
 681             fhInvMassReal[iPID]->Fill(invMass,p2->Pt());
 682         }
 683         if(fDebug>2)
 684           printf("Pair i=%d,j=%d, M=%f\n",phosPhoton1,phosPhoton2,invMass);
 685
 686         Bool_t makePi0=IsInPi0Band(invMass,p1->Pt());
 687
 688         if(makePi0 && p1->IsTagged()){//Multiple tagging
 689           fhTaggedMult->Fill(p1->Pt());
 690         }
 691         if(makePi0 && !p1->IsTagged()){//Each photon should enter histogram once, even if tagged several times
 692           fhTaggedAll->Fill(p1->Pt());
 693           Int_t iFidArea = p1->GetFiducialArea();
 694           if(iFidArea>0){
 695             fhTaggedArea1->Fill(p1->Pt() ) ;
 696             if(iFidArea>1){
 697               fhTaggedArea2->Fill(p1->Pt() ) ;
 698               if(iFidArea>2){
 699                 fhTaggedArea3->Fill(p1->Pt() ) ;
 700               }
 701             }
 702           }
 703
 704           for(Int_t iPID=0; iPID<4; iPID++){
 705             if(p1->IsPIDOK(pidMap[iPID],22))
 706               fhTaggedPID[iPID]->Fill(p1->Pt());
 707           }
 708
 709           p1->SetTagged(kTRUE) ;
 710         }
 711
 712         //First chesk if this is true pi0 pair
 713         if(IsSamePi0(p1,p2)){ //Correctly tagged - from the same pi0
 714           //        p1->SetTrueTagged(1);
 715           //        p2->SetTrueTagged(1);
 716           if(makePi0)//Correctly tagged photons
 717             fhTaggedMCTrue->Fill(p1->Pt());
 718           else{ //Decay pair missed tagging
 719             fhMCMissedTagging->Fill(p1->Pt());
 720             fhMCMissedTaggingMass->Fill(invMass,p1->Pt()) ;
 721             //Clussify why missed tagging (todo)
 722             //Converted
 723             //Partner not a photon
 724             //Tagged not a photon
 725             //Just wrong inv.mass
 726           }
 727         }
 728         else{//Fake tagged - not from the same pi0
 729           if(makePi0)//Fake pair
 730             fhMCFakeTagged->Fill(p1->Pt());
 731         }
 732       }
 733     }
 734
 735     //Fill Mixed InvMass distributions:
 736     TIter nextEv(fEventList) ;
 737     while(TClonesArray * event2 = static_cast<TClonesArray*>(nextEv())){
 738       Int_t nPhotons2 = event2->GetEntriesFast() ;
 739       for(Int_t i=0; i < inList ; i++){
 740         AliAODPWG4Particle * p1 = static_cast<AliAODPWG4Particle*>(fCaloPhotonsArr->At(i)) ;
 741         for(Int_t j=0; j < nPhotons2 ; j++){
 742           AliAODPWG4Particle * p2 = static_cast<AliAODPWG4Particle*>(event2->At(j)) ;
 743           Double_t invMass = p1->GetPairMass(p2);
 744           for(Int_t iPID=0; iPID<4; iPID++){
 745             if(p1->IsPIDOK(pidMap[iPID],22))
 746               fhInvMassMixed[iPID]->Fill(invMass,p1->Pt());
 747             if(p2->IsPIDOK(pidMap[iPID],22))
 748               fhInvMassMixed[iPID]->Fill(invMass,p2->Pt());
 749           }
 750         }
 751       }
 752     }
 753
 754     //Remove old events
 755     fEventList->AddFirst(fCaloPhotonsArr);
 756     if(fEventList->GetSize() > 10){
 757       TClonesArray *tmp = static_cast <TClonesArray*> (fEventList->Last());
 758       fEventList->Remove(tmp);
 759       delete tmp;
 760     }
 761
 762     delete caloClustersArr;
 763
 764     PostData(1, fOutputList);
 765   }// esd or aod exist
 766 }
 767
 768
 769 //______________________________________________________________________________
 770 void AliAnalysisTaskTaggedPhotons::Init()
 771 {
 772   // Intialisation of parameters
 773   AliInfo("Doing initialisation") ;
 774   SetPhotonId(0.9) ;
 775   SetMinEnergyCut(0.4);
 776   SetPi0MeanParameters(0.136,0.,0.0,0.0);
 777 //  SetPi0MeanParameters(0.1377,-0.002566,0.001216,-0.0001256);
 778   SetPi0SigmaParameters(0.004508,0.005497,0.00000006);
 779 }
 780
 781 //______________________________________________________________________________
 782 void AliAnalysisTaskTaggedPhotons::Terminate(Option_t *)
 783 {
 784   // Processing when the event loop is ended
 785   if (fDebug > 1) Printf("Terminate()");
 786 }
 787 //______________________________________________________________________________
 788 Bool_t AliAnalysisTaskTaggedPhotons::IsInPi0Band(Double_t m, Double_t pt)const
 789 {
 790   //Parameterization of the pi0 peak region
 791   Double_t mpi0mean = fPi0MeanP0 + fPi0MeanP1 * pt + fPi0MeanP2 * pt*pt + fPi0MeanP3 * pt*pt*pt;
 792   Double_t mpi0sigma = TMath::Sqrt(fPi0SigmaP0 * fPi0SigmaP0 / pt + fPi0SigmaP1 * fPi0SigmaP1 + fPi0SigmaP2 * fPi0SigmaP2 / pt / pt);
 793
 794   return (m>mpi0mean-2*mpi0sigma && m<mpi0mean+2*mpi0sigma) ;
 795 }
 796 //______________________________________________________________________________
 797 Bool_t AliAnalysisTaskTaggedPhotons::IsSamePi0(const AliAODPWG4Particle *p1, const AliAODPWG4Particle *p2)const{
 798   //Looks through parents and finds if there was commont pi0 among ancestors
 799
 800   if(!fStack)
 801     return kFALSE ; //can not say anything
 802
 803   Int_t prim1 = p1->GetLabel();
 804   while(prim1!=-1){
 805     Int_t prim2 = p2->GetLabel();
 806     while(prim2!=-1){
 807       if(prim1==prim2){
 808         if(fStack->Particle(prim1)->GetPdgCode()==111)
 809           return kTRUE ;
 810         else
 811           return kFALSE ;
 812       }
 813       prim2=fStack->Particle(prim2)->GetFirstMother() ;
 814     }
 815     prim1=fStack->Particle(prim1)->GetFirstMother() ;
 816   }
 817   return kFALSE ;
 818 }
 819 //______________________________________________________________________________
 820 Int_t AliAnalysisTaskTaggedPhotons::GetFiducialArea(Float_t * pos)const{
 821   //calculates in which kind of fiducial area photon hit
 822   Double_t phi=TMath::ATan2(pos[1],pos[0]) ;
 823   Double_t z=pos[2] ;
 824   while(phi>TMath::TwoPi())phi-=TMath::TwoPi() ;
 825   while(phi<0.)phi+=TMath::TwoPi() ;
 826   if(fDebug>2){
 827     printf("FiducialArea: phi=%f, z=%f \n",phi,z) ;
 828         printf("   fZmax=%f, fZmin=%f, fPhimax=%f, fPhimin=%f \n",fZmax,fZmin,fPhimax,fPhimin) ;
 829   }
 830   if(fPHOS){
 831     //From active PHOS area remove bands in 10 cm
 832     const Double_t kphi=TMath::ATan(10./460.) ; //angular band width
 833     Double_t dzMax=TMath::Ceil((fZmax-z)/10.) ;
 834     Double_t dzMin=TMath::Ceil((z-fZmin)/10.) ;
 835     Double_t dphiMax=TMath::Ceil((fPhimax-phi)/kphi);
 836     Double_t dphiMin=TMath::Ceil((phi-fPhimin)/kphi);
 837     if(fDebug>2){
 838      printf("In PHOS \n") ;
 839     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))) ;
 840     }
 841     return (Int_t)TMath::Min(TMath::Min(dzMax,dzMin),TMath::Min(dphiMax,dphiMin));
 842   }
 843   else{//EMCAL
 844     //From active EMCAL area remove bands in 20 cm
 845     const Double_t kphi=TMath::ATan(20./428.) ; //angular band width
 846     Double_t dzMax=TMath::Ceil((fZmax-z)/20.) ;
 847     Double_t dzMin=TMath::Ceil((z-fZmin)/20.) ;
 848     Double_t dphiMax=TMath::Ceil((fPhimax-phi)/kphi);
 849     Double_t dphiMin=TMath::Ceil((phi-fPhimin)/kphi);
 850     return (Int_t)TMath::Min(TMath::Min(dzMax,dzMin),TMath::Min(dphiMax,dphiMin));
 851   }
 852 }
 853 //______________________________________________________________________________
 854 Bool_t  AliAnalysisTaskTaggedPhotons::TestDisp(Double_t l0, Double_t l1, Double_t /*e*/)const{
 855   //test if dispersion corresponds to those of photon
 856   if(fPHOS){
 857 //    Double_t l0mean=1.38736+0.490405*TMath::Exp(-e*0.286170) ;
 858 //    Double_t l1mean=1.09786-0.323469*TMath::Exp(-e*0.918719) ;
 859 //    Double_t l0sigma=0.159905+0.829831/e-0.158067/e/e ;
 860 //    Double_t l1sigma=0.133170+0.404387/e-0.0426302/e/e ;
 861 //    Double_t c =-0.382233 ;
 862 //    return ((l0-l0mean)*(l0-l0mean)/l0sigma/l0sigma + (l1-l1mean)*(l1-l1mean)/l1sigma/l1sigma+c*(l0-l0mean)*(l1-l1mean)/l0sigma/l1sigma)<1. ;
 863
 864 if(l1>= 0   && l0>= 0   && l1 < 0.1 && l0 < 0.1) return kFALSE;
 865 if(l1>= 0   && l0 > 0.5 && l1 < 0.1 && l0 < 1.5) return kTRUE;
 866 if(l1>= 0   && l0 > 2.0 && l1 < 0.1 && l0 < 2.7) return kFALSE;
 867 if(l1>= 0   && l0 > 2.7 && l1 < 0.1 && l0 < 4.0) return kFALSE;
 868 if(l1 > 0.1 && l1 < 0.7 && l0 > 0.7 && l0 < 2.1) return kTRUE;
 869 if(l1 > 0.1 && l1 < 0.3 && l0 > 3.0 && l0 < 5.0) return kFALSE;
 870 if(l1 > 0.3 && l1 < 0.7 && l0 > 2.5 && l0 < 4.0) return kFALSE;
 871 if(l1 > 0.7 && l1 < 1.3 && l0 > 1.0 && l0 < 1.6) return kTRUE;
 872 if(l1 > 0.7 && l1 < 1.3 && l0 > 1.6 && l0 < 3.5) return kTRUE;
 873 if(l1 > 1.3 && l1 < 3.5 && l0 > 1.3 && l0 < 3.5) return kTRUE;
 874    return kFALSE ;
 875
 876   }
 877   else{ //EMCAL: not ready yet
 878    return kTRUE ;
 879
 880   }
 881
 882 }
 883 //______________________________________________________________________________^M
 884 Bool_t  AliAnalysisTaskTaggedPhotons::TestCharged(Double_t dr,Double_t /*en*/)const{
 885
 886   if(dr<15.) return kFALSE ;
 887   return kTRUE ;
 888
 889 }
 890 //______________________________________________________________________________^M
 891 void  AliAnalysisTaskTaggedPhotons::InitGeometry(){
 892   //Read rotation matrixes from ESD
 893
 894   if(fDebug>1)printf("Init geometry \n") ;
 895   AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent()) ;
 896   AliAODEvent * aod = 0x0 ;
 897   if(!esd)
 898     aod=dynamic_cast<AliAODEvent*>(InputEvent()) ;
 899   if(!esd && !aod){
 900     AliFatal("Can not read geometry matrixes from ESD/AOD: NO ESD") ;
 901   }
 902   else {
 903
 904     if(fPHOS){//reading PHOS matrixes
 905       fPHOSgeom = new AliPHOSGeoUtils("IHEP","");
 906       Bool_t activeMod[5]={0} ;
 907       for(Int_t mod=0; mod<5; mod++){
 908         if(esd){
 909           const TGeoHMatrix* m=esd->GetPHOSMatrix(mod) ;
 910           fPHOSgeom->SetMisalMatrix(m, mod) ;
 911           if(m!=0)
 912             activeMod[mod]=kTRUE ;
 913         }
 914         else{
 915           const TGeoHMatrix* m=aod->GetHeader()->GetPHOSMatrix(mod) ;
 916           fPHOSgeom->SetMisalMatrix(m, mod) ;
 917           if(m!=0)
 918             activeMod[mod]=kTRUE ;
 919         }
 920       }
 921
 922       if(fZmax>fZmin) //already  set manually
 923         return ;
 924       fZmax= 999. ;
 925       fZmin=-999. ;
 926       fPhimax=-999. ;
 927       fPhimin= 999. ;
 928       for(Int_t imod=1; imod<=5; imod++){
 929         if(!activeMod[imod-1])
 930           continue ;
 931         //Find exact coordinates of PHOS corners
 932         Int_t relId[4]={imod,0,1,1} ;
 933         Float_t x,z ;
 934         fPHOSgeom->RelPosInModule(relId,x,z) ;
 935         TVector3 pos ;
 936         fPHOSgeom->Local2Global(imod,x,z,pos) ;
 937         Double_t phi=pos.Phi() ;
 938         while(phi<0.)phi+=TMath::TwoPi() ;
 939         while(phi>TMath::TwoPi())phi-=TMath::TwoPi() ;
 940         fPhimin=TMath::Min(fPhimin,float(phi)) ;
 941         fZmin=TMath::Max(fZmin,float(pos.Z())) ;
 942         relId[2]=64 ;
 943         relId[3]=56 ;
 944         fPHOSgeom->RelPosInModule(relId,x,z) ;
 945         fPHOSgeom->Local2Global(imod,x,z,pos) ;
 946         phi=pos.Phi() ;
 947         while(phi<0.)phi+=TMath::TwoPi() ;
 948         while(phi>TMath::TwoPi())phi-=TMath::TwoPi() ;
 949         fPhimax=TMath::Max(fPhimax,float(phi)) ;
 950         fZmax=TMath::Min(fZmax,float(pos.Z())) ;
 951       }
 952     }
 953     else{ //EMCAL
 954       fEMCALgeom = new AliEMCALGeoUtils("EMCAL_FIRSTYEAR");
 955       for(Int_t mod=0; mod < (fEMCALgeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
 956         if(esd){
 957           const TGeoHMatrix* m=esd->GetEMCALMatrix(mod) ;
 958           fEMCALgeom->SetMisalMatrix(m, mod) ;
 959         }
 960         else{
 961           const TGeoHMatrix* m=aod->GetHeader()->GetEMCALMatrix(mod) ;
 962           fEMCALgeom->SetMisalMatrix(m, mod) ;
 963         }
 964       }
 965       if(fZmax>fZmin) //already  set manually
 966         return ;
 967
 968       fZmax= 999. ;
 969       fZmin=-999. ;
 970       fPhimax=-999. ;
 971       fPhimin= 999. ;
 972       for(Int_t imod=0; imod<12; imod++){
 973         //Find exact coordinates of SM corners
 974         Int_t absId = fEMCALgeom->GetAbsCellIdFromCellIndexes(imod, 0, 0);
 975         TVector3 pos ;
 976         //Get the position of this tower.
 977         fEMCALgeom->RelPosCellInSModule(absId,pos);
 978         Double_t phi=pos.Phi() ;
 979         while(phi<0.)phi+=TMath::TwoPi() ;
 980         while(phi>TMath::TwoPi())phi-=TMath::TwoPi() ;
 981         fPhimin=TMath::Min(fPhimin,float(phi)) ;
 982         fZmin=TMath::Max(fZmin,float(pos.Z())) ;
 983         absId = fEMCALgeom->GetAbsCellIdFromCellIndexes(imod, 24, 48);
 984         fEMCALgeom->RelPosCellInSModule(absId,pos);
 985         phi=pos.Phi() ;
 986         while(phi<0.)phi+=TMath::TwoPi() ;
 987         while(phi>TMath::TwoPi())phi-=TMath::TwoPi() ;
 988         fPhimax=TMath::Max(fPhimax,float(phi)) ;
 989         fZmax=TMath::Min(fZmax,float(pos.Z())) ;
 990       }
 991     }
 992   } // ESD or AOD available
 993 }