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