--- /dev/null
+/**************************************************************************\r
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *\r
+ * *\r
+ * Author: The ALICE Off-line Project. *\r
+ * Contributors are mentioned in the code where appropriate. *\r
+ * *\r
+ * Permission to use, copy, modify and distribute this software and its *\r
+ * documentation strictly for non-commercial purposes is hereby granted *\r
+ * without fee, provided that the above copyright notice appears in all *\r
+ * copies and that both the copyright notice and this permission notice *\r
+ * appear in the supporting documentation. The authors make no claims *\r
+ * about the suitability of this software for any purpose. It is *\r
+ * provided "as is" without express or implied warranty. *\r
+ **************************************************************************/\r
+\r
+/* $Id: */\r
+\r
+//_________________________________________________________________________\r
+// Analysis for Tagged Photons\r
+// Prepares all necessary histograms for calculation of \r
+// the yield of pi0 decay photon in calorimeter:\r
+// \r
+//\r
+//*-- Dmitry Blau \r
+//////////////////////////////////////////////////////////////////////////////\r
+\r
+#include <TCanvas.h>\r
+#include <TH1.h>\r
+#include <TH2.h>\r
+#include <TH3.h>\r
+#include <TROOT.h>\r
+#include <TLegend.h> \r
+#include <TNtuple.h>\r
+#include <TVector3.h> \r
+#include <TFile.h>\r
+\r
+#include "AliAnalysisTaskTaggedPhotons.h" \r
+#include "AliAnalysisManager.h"\r
+#include "AliESDEvent.h" \r
+#include "AliESDCaloCluster.h" \r
+#include "AliAODPWG4Particle.h"\r
+#include "AliLog.h"\r
+#include "AliESDVertex.h"\r
+#include "AliPHOSGeoUtils.h"\r
+#include "TGeoManager.h"\r
+#include "AliMCAnalysisUtils.h"
+#include "AliMCEventHandler.h"\r
+#include "AliAnalysisManager.h"\r
+#include "AliMCEvent.h"\r
+#include "AliStack.h"\r
+\r
+//______________________________________________________________________________\r
+AliAnalysisTaskTaggedPhotons::AliAnalysisTaskTaggedPhotons() : AliAnalysisTaskSE(),
+ fgeom(0x0),fStack(0x0),fDebug(0),fPHOS(1),
+ fPhotonId(1.0),fMinEnergyCut(0.4),
+ fPi0Mean_p0(0.1377),fPi0Mean_p1(-0.002566),fPi0Mean_p2(0.001216),fPi0Mean_p3(-0.0001256),
+ fPi0Sigma_p0(0.004508),fPi0Sigma_p1(0.005497),fPi0Sigma_p2(0.00000006),
+ fZmax(0.),fZmin(0.),fPhimax(0.),fPhimin(0.),
+ fOutputList(0x0),fEventList(0x0),
+ fhRecAll(0x0),fhRecAllArea1(0x0),fhRecAllArea2(0x0),fhRecAllArea3(0x0),
+ fhRecPhoton(0x0),fhRecOther(0x0),fhRecPhotPi0(0x0),fhRecPhotEta(0x0),
+ fhRecPhotOmega(0x0),fhRecPhotEtapr(0x0),fhRecPhotConv(0x0),fhRecPhotHadron(0x0),
+ fhRecPhotDirect(0x0),fhRecPhotOther(0x0),
+ fhDecWMCPartner(0x0),fhDecWMissedPartnerNotPhoton(0x0),fhDecWMissedPartnerAll(0x0),
+ fhDecWMissedPartnerEmin(0x0),fhDecWMissedPartnerConv(0x0),fhDecWMissedPartnerGeom0(0x0),
+ fhDecWMissedPartnerGeom1(0x0),fhDecWMissedPartnerGeom2(0x0),fhDecWMissedPartnerGeom3(0x0),
+ fhPartnerMCReg(0x0),fhPartnerMissedEmin(0x0),fhPartnerMissedConv(0x0),fhPartnerMissedGeo(0x0),
+ fhTaggedAll(0x0),fhTaggedArea1(0x0),fhTaggedArea2(0x0),fhTaggedArea3(0x0),fhTaggedMult(0x0),
+ fhTaggedMCTrue(0x0),fhMCMissedTagging(0x0),fhMCFakeTagged(0x0),
+ fhInvMassReal(0x0),fhInvMassMixed(0x0),fhMCMissedTaggingMass(0x0),
+ fhConversionRadius(0x0),fhInteractionRadius(0x0),fhEvents(0x0)
+{\r
+ //Default constructor\r
+}\r
+//______________________________________________________________________________\r
+AliAnalysisTaskTaggedPhotons::AliAnalysisTaskTaggedPhotons(const char *name) : \r
+ AliAnalysisTaskSE(name),\r
+ fgeom(0x0),fStack(0x0),fDebug(0),fPHOS(1),
+ fPhotonId(1.0),fMinEnergyCut(0.4),
+ fPi0Mean_p0(0.1377),fPi0Mean_p1(-0.002566),fPi0Mean_p2(0.001216),fPi0Mean_p3(-0.0001256),
+ fPi0Sigma_p0(0.004508),fPi0Sigma_p1(0.005497),fPi0Sigma_p2(0.00000006),
+ fZmax(0.),fZmin(0.),fPhimax(0.),fPhimin(0.),
+ fOutputList(0x0),fEventList(0x0),
+ fhRecAll(0x0),fhRecAllArea1(0x0),fhRecAllArea2(0x0),fhRecAllArea3(0x0),
+ fhRecPhoton(0x0),fhRecOther(0x0),fhRecPhotPi0(0x0),fhRecPhotEta(0x0),
+ fhRecPhotOmega(0x0),fhRecPhotEtapr(0x0),fhRecPhotConv(0x0),fhRecPhotHadron(0x0),
+ fhRecPhotDirect(0x0),fhRecPhotOther(0x0),
+ fhDecWMCPartner(0x0),fhDecWMissedPartnerNotPhoton(0x0),fhDecWMissedPartnerAll(0x0),
+ fhDecWMissedPartnerEmin(0x0),fhDecWMissedPartnerConv(0x0),fhDecWMissedPartnerGeom0(0x0),
+ fhDecWMissedPartnerGeom1(0x0),fhDecWMissedPartnerGeom2(0x0),fhDecWMissedPartnerGeom3(0x0),
+ fhPartnerMCReg(0x0),fhPartnerMissedEmin(0x0),fhPartnerMissedConv(0x0),fhPartnerMissedGeo(0x0),
+ fhTaggedAll(0x0),fhTaggedArea1(0x0),fhTaggedArea2(0x0),fhTaggedArea3(0x0),fhTaggedMult(0x0),
+ fhTaggedMCTrue(0x0),fhMCMissedTagging(0x0),fhMCFakeTagged(0x0),
+ fhInvMassReal(0x0),fhInvMassMixed(0x0),fhMCMissedTaggingMass(0x0),
+ fhConversionRadius(0x0),fhInteractionRadius(0x0),fhEvents(0x0)
+{\r
+ // Constructor.\r
+\r
+ // Output slots \r
+ DefineOutput(1, TList::Class()) ; \r
+}\r
+\r
+//____________________________________________________________________________\r
+AliAnalysisTaskTaggedPhotons::AliAnalysisTaskTaggedPhotons(const AliAnalysisTaskTaggedPhotons& ap) :\r
+ AliAnalysisTaskSE(ap.GetName()), \r
+ fgeom(0x0),fStack(0x0),fDebug(ap.fDebug),fPHOS(ap.fPHOS),
+ fPhotonId(ap.fPhotonId),fMinEnergyCut(ap.fMinEnergyCut),
+ fPi0Mean_p0(ap.fPi0Mean_p0),fPi0Mean_p1(ap.fPi0Mean_p1),fPi0Mean_p2(ap.fPi0Mean_p2),fPi0Mean_p3(ap.fPi0Mean_p3),
+ fPi0Sigma_p0(ap.fPi0Sigma_p0),fPi0Sigma_p1(ap.fPi0Sigma_p1),fPi0Sigma_p2(ap.fPi0Sigma_p2),
+ fZmax(ap.fZmax),fZmin(ap.fZmin),fPhimax(ap.fPhimax),fPhimin(ap.fPhimin),
+ fOutputList(0x0),fEventList(0x0),
+ fhRecAll(0x0),fhRecAllArea1(0x0),fhRecAllArea2(0x0),fhRecAllArea3(0x0),
+ fhRecPhoton(0x0),fhRecOther(0x0),fhRecPhotPi0(0x0),fhRecPhotEta(0x0),
+ fhRecPhotOmega(0x0),fhRecPhotEtapr(0x0),fhRecPhotConv(0x0),fhRecPhotHadron(0x0),
+ fhRecPhotDirect(0x0),fhRecPhotOther(0x0),
+ fhDecWMCPartner(0x0),fhDecWMissedPartnerNotPhoton(0x0),fhDecWMissedPartnerAll(0x0),
+ fhDecWMissedPartnerEmin(0x0),fhDecWMissedPartnerConv(0x0),fhDecWMissedPartnerGeom0(0x0),
+ fhDecWMissedPartnerGeom1(0x0),fhDecWMissedPartnerGeom2(0x0),fhDecWMissedPartnerGeom3(0x0),
+ fhPartnerMCReg(0x0),fhPartnerMissedEmin(0x0),fhPartnerMissedConv(0x0),fhPartnerMissedGeo(0x0),
+ fhTaggedAll(0x0),fhTaggedArea1(0x0),fhTaggedArea2(0x0),fhTaggedArea3(0x0),fhTaggedMult(0x0),
+ fhTaggedMCTrue(0x0),fhMCMissedTagging(0x0),fhMCFakeTagged(0x0),
+ fhInvMassReal(0x0),fhInvMassMixed(0x0),fhMCMissedTaggingMass(0x0),
+ fhConversionRadius(0x0),fhInteractionRadius(0x0),fhEvents(0x0)
+{
+ // cpy ctor\r
+}\r
+\r
+//_____________________________________________________________________________\r
+AliAnalysisTaskTaggedPhotons& AliAnalysisTaskTaggedPhotons::operator = (const AliAnalysisTaskTaggedPhotons& ap)\r
+{\r
+// assignment operator\r
+\r
+ this->~AliAnalysisTaskTaggedPhotons();\r
+ new(this) AliAnalysisTaskTaggedPhotons(ap);\r
+ return *this;\r
+}\r
+\r
+//______________________________________________________________________________\r
+AliAnalysisTaskTaggedPhotons::~AliAnalysisTaskTaggedPhotons()\r
+{\r
+ // dtor\r
+ if(fOutputList) {\r
+ fOutputList->Clear() ; \r
+ delete fOutputList ;\r
+ }\r
+}\r
+\r
+\r
+//________________________________________________________________________\r
+void AliAnalysisTaskTaggedPhotons::UserCreateOutputObjects()\r
+{ \r
+\r
+ //Load geometry\r
+ //if file "geometry.root" exists, load it\r
+ //otherwise use misaligned matrixes stored in ESD\r
+ TFile *geoFile = new TFile("geometry.root","read");\r
+ if(geoFile->IsZombie()){ //no file, load geo matrixes from ESD\r
+ //todo\r
+ AliInfo("Can not find file geometry.root, reading misalignment matrixes from AliESDs") ;
+ AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent()) ;
+ if(!esd)
+ AliFatal("Can not read geometry even from ESD. Note, that AOD does not contain PHOS/EMCAL geometry") ;
+ if(fPHOS){//reading PHOS matrixes
+ fgeom = new AliPHOSGeoUtils("IHEP","");\r
+ for(Int_t mod=0; mod<5; mod++){
+ const TGeoHMatrix* m=esd->GetPHOSMatrix(mod) ;
+ fgeom->SetMisalMatrix(m, mod) ;
+ }
+ }
+ }\r
+ else{\r
+ gGeoManager = (TGeoManager*)geoFile->Get("Geometry");\r
+ fgeom = new AliPHOSGeoUtils("IHEP","");\r
+ }\r
+\r
+\r
+ if(fgeom==NULL){\r
+ AliError("Error loading Geometry\n");\r
+ }\r
+ else\r
+ AliInfo("Geometry loaded... OK\n");\r
+\r
+ //Evaluate active PHOS/EMCAL area\r
+ //To be fixed todo\r
+ fZmax= 2.25*56/2. ;\r
+ fZmin=-2.25*56/2. ;\r
+ fPhimax=220./180.*TMath::Pi() ;\r
+ fPhimin=320./180.*TMath::Pi() ;\r
+\r
+\r
+ // Create the outputs containers\r
+\r
+ OpenFile(1) ; \r
+ const Int_t nPtBins=51 ;\r
+ Double_t ptBins[nPtBins+1] ;\r
+ for(Int_t i=0;i<=20;i++)ptBins[i]=0.1*i ; //0-2 GeV: 0.1 GeV/bin\r
+ for(Int_t i=21;i<=30;i++)ptBins[i]=2.+0.2*(i-20) ; //2-4 GeV: 0.2 GeV/bin \r
+ for(Int_t i=31;i<=38;i++)ptBins[i]=4.+0.5*(i-30) ; //4-8 GeV: 0.5 GeV/bin \r
+ for(Int_t i=39;i<=42;i++)ptBins[i]=8.+1.0*(i-38) ; //8-12 GeV: 1. GeV/bin \r
+ for(Int_t i=43;i<=52;i++)ptBins[i]=12.+2.0*(i-42) ; //12-30 GeV: 2. GeV/bin \r
+\r
+\r
+ //Reconstructed spectra\r
+ fhRecAll = new TH1D("fhRecAll","Spectrum of all reconstructed particles", nPtBins, ptBins ) ;\r
+ fhRecAllArea1 = new TH1D("fhRecAllArea1","Spectrum of rec particles in Fid. Area 1", nPtBins, ptBins ) ;\r
+ fhRecAllArea2 = new TH1D("fhRecAllArea2","Spectrum of rec particles in Fid. Area 2", nPtBins, ptBins ) ;\r
+ fhRecAllArea3 = new TH1D("fhRecAllArea3","Spectrum of rec particles in Fid. Area 3", nPtBins, ptBins ) ;\r
+\r
+ //Sort registered particles spectra according MC information\r
+ fhRecPhoton = new TH1D("fhRecPhoton","Spectrum of rec. with primary==22 and no PID criteria", nPtBins, ptBins ) ;\r
+ fhRecOther = new TH1D("fhRecOther"," Spectrum of rec. with primary!=22 and no PID criteria", nPtBins, ptBins ) ;\r
+ fhRecPhotonPID[0]= new TH1D("fhRecPhotonPID0","Spectrum of rec. with primary==22 and PID=0", nPtBins, ptBins ) ;\r
+ fhRecPhotonPID[1]= new TH1D("fhRecPhotonPID1","Spectrum of rec. with primary==22 and PID=1", nPtBins, ptBins ) ;\r
+ fhRecPhotonPID[2]= new TH1D("fhRecPhotonPID2","Spectrum of rec. with primary==22 and PID=2", nPtBins, ptBins ) ;\r
+ fhRecPhotonPID[3]= new TH1D("fhRecPhotonPID3","Spectrum of rec. with primary==22 and PID=3", nPtBins, ptBins ) ;\r
+ fhRecOtherPID[0]= new TH1D("fhRecOtherPID0","Spectrum of rec. with primary!=22 and PID=0", nPtBins, ptBins ) ;\r
+ fhRecOtherPID[1]= new TH1D("fhRecOtherPID1","Spectrum of rec. with primary!=22 and PID=1", nPtBins, ptBins ) ;\r
+ fhRecOtherPID[2]= new TH1D("fhRecOtherPID2","Spectrum of rec. with primary!=22 and PID=2", nPtBins, ptBins ) ;\r
+ fhRecOtherPID[3]= new TH1D("fhRecOtherPID3","Spectrum of rec. with primary!=22 and PID=3", nPtBins, ptBins ) ;\r
+ fhRecPhotPi0 = new TH1D("fhRecPhotPi0","Spectrum of rec. photons from pi0 decays", nPtBins, ptBins ) ;\r
+ fhRecPhotEta = new TH1D("fhRecPhotEta","Spectrum of rec. photons from eta decays", nPtBins, ptBins ) ;\r
+ fhRecPhotOmega = new TH1D("fhRecPhotOmega","Spectrum of rec. photons from omega decays", nPtBins, ptBins ) ;\r
+ fhRecPhotEtapr = new TH1D("fhRecPhotEtapr","Spectrum of rec. photons from eta prime decays", nPtBins, ptBins ) ;\r
+ fhRecPhotConv = new TH1D("fhRecPhotConv"," Spectrum of rec. photons from conversion", nPtBins, ptBins ) ;\r
+ fhRecPhotHadron = new TH1D("fhRecPhotHadron","Spectrum of rec. photons from hadron-matter interactions", nPtBins, ptBins ) ;\r
+ fhRecPhotDirect = new TH1D("fhRecPhotDirect","Spectrum of rec. photons direct or no primary", nPtBins, ptBins ) ;\r
+ fhRecPhotOther = new TH1D("fhRecPhotOther","Spectrum of rec. photons from other hadron decays", nPtBins, ptBins ) ;\r
+\r
+ //MC tagging: reasons of partner loss etc.\r
+ fhDecWMCPartner = new TH1D("fhDecWMCPartner","pi0 decay photon which partner should be registered according to MC", nPtBins, ptBins ) ;\r
+ fhDecWMissedPartnerNotPhoton = new TH1D("fhDecWMissedPartnerNotPhoton","Decay photon with missed non-photon partner", nPtBins, ptBins ) ;\r
+ fhDecWMissedPartnerAll = new TH1D("fhDecWMissedPartnerAll","Decay photons with partner missed due to some reason", nPtBins, ptBins ) ;\r
+ fhDecWMissedPartnerEmin = new TH1D("fhDecWMissedPartnerEmin","Decay photons with partner missed due to low energy", nPtBins, ptBins ) ;\r
+ fhDecWMissedPartnerConv = new TH1D("fhDecWMissedPartnerConv","Decay photons with partner missed due to conversion", nPtBins, ptBins ) ;\r
+ fhDecWMissedPartnerGeom0 = new TH1D("fhDecWMissedPartnerGeom0","Decay photons with partner missed due geometry", nPtBins, ptBins ) ;\r
+ fhDecWMissedPartnerGeom1 = new TH1D("fhDecWMissedPartnerGeom1","Decay photons with partner missed due geometry Fid. area. 1", nPtBins, ptBins ) ;\r
+ fhDecWMissedPartnerGeom2 = new TH1D("fhDecWMissedPartnerGeom2","Decay photons with partner missed due geometry Fid. area. 2", nPtBins, ptBins ) ;\r
+ fhDecWMissedPartnerGeom3 = new TH1D("fhDecWMissedPartnerGeom3","Decay photons with partner missed due geometry Fid. area. 3", nPtBins, ptBins ) ;\r
+\r
+ //MC tagging: Decay partners spectra\r
+ fhPartnerMCReg = new TH1D("fhPartnerMCReg","Spectrum of decay partners which should be registered (MC)", nPtBins, ptBins ) ;\r
+ fhPartnerMissedEmin = new TH1D("fhPartnerMissedEmin","Spectrum of decay partners lost due to Emin cut", nPtBins, ptBins ) ;\r
+ fhPartnerMissedConv = new TH1D("fhPartnerMissedConv","Spectrum of decay partners lost due to conversion", nPtBins, ptBins ) ;\r
+ fhPartnerMissedGeo = new TH1D("fhPartnerMissedGeo","Spectrum of decay partners lost due to acceptance", nPtBins, ptBins ) ;\r
+\r
+ //Tagging\r
+ fhTaggedAll = new TH1D("fhTaggedAll","Spectrum of all tagged photons", nPtBins, ptBins ) ;\r
+ fhTaggedArea1 = new TH1D("fhTaggedArea1","Spectrum of all tagged photons Fid. area1", nPtBins, ptBins ) ;\r
+ fhTaggedArea2 = new TH1D("fhTaggedArea2","Spectrum of all tagged photons Fid. area2", nPtBins, ptBins ) ;\r
+ fhTaggedArea3 = new TH1D("fhTaggedArea3","Spectrum of all tagged photons Fid. area3", nPtBins, ptBins ) ;\r
+ fhTaggedPID[0] = new TH1D("fhTaggedPID0","Spectrum of tagged photons for PID=0", nPtBins, ptBins ) ;\r
+ fhTaggedPID[1] = new TH1D("fhTaggedPID1","Spectrum of tagged photons for PID=1", nPtBins, ptBins ) ;\r
+ fhTaggedPID[2] = new TH1D("fhTaggedPID2","Spectrum of tagged photons for PID=2", nPtBins, ptBins ) ;\r
+ fhTaggedPID[3] = new TH1D("fhTaggedPID3","Spectrum of tagged photons for PID=3", nPtBins, ptBins ) ;\r
+ fhTaggedMult = new TH1D("fhTaggedMult","Spectrum of multiply tagged photons", nPtBins, ptBins ) ;\r
+\r
+ //Tagging: use MC information if available\r
+ fhTaggedMCTrue = new TH1D("fhTaggedMCTrue","Spectrum of correctly tagged pi0 decay photons", nPtBins, ptBins ) ;\r
+ fhMCMissedTagging = new TH1D("fhMCMissedTagging","Spectrum of pi0 decay photons missed tagging due to wrong pair mass", nPtBins, ptBins ) ;\r
+ fhMCFakeTagged = new TH1D("fhMCFakeTagged","Spectrum of photons wrongly tagged according to MC", nPtBins, ptBins ) ;\r
+\r
+ //Invariant mass distributions for fake corrections\r
+ const Int_t nmass=200 ;\r
+ Double_t masses[nmass+1] ;\r
+ for(Int_t i=0;i<=nmass;i++)masses[i]=0.005*i ;\r
+ fhInvMassReal = new TH2D("fhInvMassReal","Two-photon inv. mass vs first photon pt",nmass,masses,nPtBins, ptBins ) ;\r
+ fhInvMassMixed = new TH2D("fhInvMassMixed","Two-photon inv. mass vs first photon pt",nmass,masses,nPtBins, ptBins ) ;\r
+ fhMCMissedTaggingMass= new TH2D("fhMCMissedTaggingMass","Inv mass of pairs missed tagging",nmass,masses,nPtBins, ptBins ) ;\r
+\r
+ //Conversion and annihilation radius distributions\r
+ fhConversionRadius = new TH1D("fhConversionRadius","Radis of photon production (conversion)",100,0.,500.) ;\r
+ fhInteractionRadius = new TH1D("fhInteractionRadius","Radis of photon production (hadron interaction)",100,0.,500.);\r
+\r
+ fhEvents = new TH1D("hEvents", "Number of Events processed", 1, 0, 1);\r
+\r
+// create output container\r
+ \r
+ fOutputList = new TList() ;\r
+ fEventList = new TList() ; \r
+ fOutputList->SetName(GetName()) ; \r
+\r
+ fOutputList->AddAt(fhRecAll, 0) ;\r
+ fOutputList->AddAt(fhRecAllArea1, 1) ;\r
+ fOutputList->AddAt(fhRecAllArea2, 2) ;\r
+ fOutputList->AddAt(fhRecAllArea3, 3) ;\r
+\r
+ fOutputList->AddAt(fhRecPhoton, 4) ;\r
+ fOutputList->AddAt(fhRecOther, 5) ;\r
+ fOutputList->AddAt(fhRecPhotonPID[0], 6) ;\r
+ fOutputList->AddAt(fhRecPhotonPID[1], 7) ;\r
+ fOutputList->AddAt(fhRecPhotonPID[2], 8) ;\r
+ fOutputList->AddAt(fhRecPhotonPID[3], 9) ;\r
+ fOutputList->AddAt(fhRecOtherPID[0], 10) ;\r
+ fOutputList->AddAt(fhRecOtherPID[1], 11) ;\r
+ fOutputList->AddAt(fhRecOtherPID[2], 12) ;\r
+ fOutputList->AddAt(fhRecOtherPID[3], 13) ;\r
+ fOutputList->AddAt(fhRecPhotPi0, 14) ;\r
+ fOutputList->AddAt(fhRecPhotEta, 15) ;\r
+ fOutputList->AddAt(fhRecPhotOmega, 16) ;\r
+ fOutputList->AddAt(fhRecPhotEtapr, 17) ;\r
+ fOutputList->AddAt(fhRecPhotConv, 18) ;\r
+ fOutputList->AddAt(fhRecPhotHadron, 19) ;\r
+ fOutputList->AddAt(fhRecPhotDirect, 20) ;\r
+ fOutputList->AddAt(fhRecPhotOther, 21) ;\r
+\r
+ fOutputList->AddAt(fhDecWMCPartner, 22) ;\r
+ fOutputList->AddAt(fhDecWMissedPartnerNotPhoton, 23) ;\r
+ fOutputList->AddAt(fhDecWMissedPartnerAll, 24) ;\r
+ fOutputList->AddAt(fhDecWMissedPartnerEmin, 25) ;\r
+ fOutputList->AddAt(fhDecWMissedPartnerConv, 26) ;\r
+ fOutputList->AddAt(fhDecWMissedPartnerGeom0, 27) ;\r
+ fOutputList->AddAt(fhDecWMissedPartnerGeom1, 28) ;\r
+ fOutputList->AddAt(fhDecWMissedPartnerGeom2, 29) ;\r
+ fOutputList->AddAt(fhDecWMissedPartnerGeom3, 30) ;\r
+\r
+ fOutputList->AddAt(fhPartnerMCReg, 31) ;\r
+ fOutputList->AddAt(fhPartnerMissedEmin, 32) ;\r
+ fOutputList->AddAt(fhPartnerMissedConv, 33) ;\r
+ fOutputList->AddAt(fhPartnerMissedGeo, 34) ;\r
+\r
+ fOutputList->AddAt(fhTaggedAll, 35) ;\r
+ fOutputList->AddAt(fhTaggedArea1, 36) ;\r
+ fOutputList->AddAt(fhTaggedArea2, 37) ;\r
+ fOutputList->AddAt(fhTaggedArea3, 38) ;\r
+ fOutputList->AddAt(fhTaggedPID[0], 39) ;\r
+ fOutputList->AddAt(fhTaggedPID[1], 40) ;\r
+ fOutputList->AddAt(fhTaggedPID[2], 41) ;\r
+ fOutputList->AddAt(fhTaggedPID[3], 42) ;\r
+ fOutputList->AddAt(fhTaggedMult, 43) ;\r
+\r
+ fOutputList->AddAt(fhTaggedMCTrue, 44) ;\r
+ fOutputList->AddAt(fhMCMissedTagging, 45) ;\r
+ fOutputList->AddAt(fhMCFakeTagged, 46) ;\r
+\r
+ fOutputList->AddAt(fhInvMassReal, 47) ;\r
+ fOutputList->AddAt(fhInvMassMixed, 48) ;\r
+ fOutputList->AddAt(fhMCMissedTaggingMass, 49) ;\r
+\r
+ fOutputList->AddAt(fhConversionRadius, 50) ;\r
+ fOutputList->AddAt(fhInteractionRadius, 51) ;\r
+\r
+ fOutputList->AddAt(fhEvents, 52) ;\r
+\r
+\r
+/*\r
+ fOutputList->AddAt(fhPHOSPos, 0) ; \r
+ fOutputList->AddAt(fhPHOS, 1) ; \r
+ fOutputList->AddAt(fhAllPhotons, 2) ;\r
+ fOutputList->AddAt(fhNotPhotons, 3) ;\r
+ fOutputList->AddAt(fhAllPhotonsPrimary, 4) ;\r
+ fOutputList->AddAt(fhNotPhotonsPrimary, 5) ;\r
+ fOutputList->AddAt(fhfakeNotPhotons, 6) ;\r
+\r
+ fOutputList->AddAt(fhTaggedPhotons, 7) ; \r
+ fOutputList->AddAt(fhfakeTaggedPhotons, 8) ;\r
+ fOutputList->AddAt(fhDecayNotTaggedPhotons, 9) ;\r
+ fOutputList->AddAt(fhPi0DecayPhotonsPrimary, 10);\r
+ fOutputList->AddAt(fhEtaDecayPhotonsPrimary, 11);\r
+ fOutputList->AddAt(fhOmegaDecayPhotonsPrimary, 12);\r
+ fOutputList->AddAt(fhEtaSDecayPhotonsPrimary, 13);\r
+ fOutputList->AddAt(fhOtherDecayPhotonsPrimary, 14);\r
+ fOutputList->AddAt(fhDecayPhotonsPrimary, 15);\r
+\r
+ fOutputList->AddAt(fhConvertedPhotonsPrimary, 16);\r
+ fOutputList->AddAt(fhConvertedPhotonsPrimaryHadronsDecays, 17);\r
+ fOutputList->AddAt(fhCoordsConversion, 18);\r
+ fOutputList->AddAt(fhCoordsConversion2, 19);\r
+\r
+ fOutputList->AddAt(fhPHOSInvariantMassReal, 20);\r
+ fOutputList->AddAt(fhPHOSInvariantMassMixed, 21);\r
+\r
+ fOutputList->AddAt(fhPHOSPi0, 22);\r
+\r
+ fOutputList->AddAt(fhPi0DecayPhotonsGeomfake, 23);\r
+ fOutputList->AddAt(fhPi0DecayPhotonsTaggedPrimary, 24);\r
+ fOutputList->AddAt(fhPi0DecayPhotonsTaggedPrimaryPair, 25);\r
+ fOutputList->AddAt(fhPi0DecayPhotonsBigDecay, 26);\r
+ fOutputList->AddAt(fhPi0DecayPhotonsPConv, 27);\r
+ fOutputList->AddAt(fhPi0DecayPhotonsPGeo, 28);\r
+ fOutputList->AddAt(fhPi0DecayPhotonsPReg, 29);\r
+\r
+ fOutputList->AddAt(fhfakeTaggedPhotonsConv, 30) ;\r
+ fOutputList->AddAt(fhfakeTaggedPhotonsPID, 31) ;\r
+ fOutputList->AddAt(fhstrangeNotTaggedPhotons, 32) ;\r
+ fOutputList->AddAt(fhstrangeNotTaggedPhotonsPair, 33) ;\r
+ fOutputList->AddAt(fhstrangeNotTaggedPhotonsRegCut, 34) ;\r
+ fOutputList->AddAt(fhPi0DecayPhotonsTaggedPrimaryRegCut, 35);\r
+ fOutputList->AddAt(fhstrangeNotTaggedPhotonsPairRegCut, 36) ;\r
+ fOutputList->AddAt(fhPi0DecayPhotonsTaggedPrimaryPairRegCut, 37);\r
+ fOutputList->AddAt(fhTrackRefCoords, 38);\r
+*/\r
+ fOutputList->AddAt(fhEvents, 39);\r
+}\r
+\r
+//______________________________________________________________________________\r
+void AliAnalysisTaskTaggedPhotons::UserExec(Option_t *) \r
+{\r
+ fhEvents->Fill(0.);\r
+\r
+ // Processing of one event\r
+ if(fDebug>1)\r
+ AliInfo(Form("\n\n Processing event # %lld", Entry())) ; \r
+ AliESDEvent* esd = (AliESDEvent*)InputEvent();\r
+\r
+ if(fDebug>2){ //DP: Check these histograms <------- \r
+// printf("Tagged: %f ",fhTaggedPhotons->GetEntries());\r
+// printf("fakeTagged: %f ",fhfakeTaggedPhotons->GetEntries());\r
+// printf("fakeNotTagged: %f ",fhfakeNotTaggedPhotons->GetEntries());\r
+// printf("strangeNotTagged: %f ",fhstrangeNotTaggedPhotons->GetEntries());\r
+// printf("trueTagged: %f ",fhPi0DecayPhotonsTrueTagged->GetEntries());\r
+// printf("fakeTaggedConv: %f ",fhfakeTaggedPhotonsConv->GetEntries());\r
+// printf("\nDIFF: %f\n",fhTaggedPhotons->GetEntries()-fhfakeTaggedPhotons->GetEntries()+fhfakeNotTaggedPhotons->GetEntries()+fhstrangeNotTaggedPhotons->GetEntries()-fhPi0DecayPhotonsTrueTagged->GetEntries()-fhfakeTaggedPhotonsConv->GetEntries());\r
+ }\r
+\r
+ //MC stack init\r
+ AliMCEventHandler* mctruth = (AliMCEventHandler*)((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler());\r
+ fStack = mctruth->MCEvent()->Stack();\r
+\r
+ if(!fStack && gDebug>1)\r
+ AliInfo("No stack! \n");\r
+\r
+ //************************ PHOS *************************************\r
+ TRefArray * caloClustersArr = new TRefArray(); \r
+ esd->GetPHOSClusters(caloClustersArr);\r
+ const Int_t kNumberOfPhosClusters = caloClustersArr->GetEntries() ; \r
+\r
+ TClonesArray * fCaloPhotonsArr = new TClonesArray("AliAODPWG4Particle",kNumberOfPhosClusters);\r
+ Int_t inList = 0; //counter of caloClusters\r
+\r
+ Int_t phosCluster ; \r
+\r
+ // loop over PHOS Clusters\r
+ for(phosCluster = 0 ; phosCluster < kNumberOfPhosClusters ; phosCluster++) {\r
+ AliESDCaloCluster * caloCluster = static_cast<AliESDCaloCluster*>(caloClustersArr->At(phosCluster)) ;\r
+ \r
+ if((fPHOS && !caloCluster->IsPHOS()) ||\r
+ (!fPHOS && caloCluster->IsPHOS()))\r
+ continue ; \r
+\r
+ Double_t v[3] ; //vertex ;\r
+ esd->GetVertex()->GetXYZ(v) ;\r
+\r
+ TLorentzVector momentum ;\r
+ caloCluster->GetMomentum(momentum, v);\r
+ new ((*fCaloPhotonsArr)[inList]) AliAODPWG4Particle(momentum.Px(),momentum.Py(),momentum.Pz(),caloCluster->E() );\r
+ AliAODPWG4Particle *p = static_cast<AliAODPWG4Particle*>(fCaloPhotonsArr->At(inList));\r
+ inList++;\r
+\r
+ p->SetCaloLabel(phosCluster,-1); //This and partner cluster
+ p->SetDistToBad(Int_t(caloCluster->GetDistanceToBadChannel()));\r
+
+ p->SetTag(AliMCAnalysisUtils::kMCUnknown);\r
+ p->SetTagged(kFALSE); //Reconstructed pairs found\r
+ p->SetLabel(caloCluster->GetLabel());\r
+ Float_t pos[3] ;
+ caloCluster->GetPosition(pos) ;
+ p->SetFiducialArea(GetFiducialArea(pos)) ;
+\r
+ fhRecAll->Fill( p->Pt() ) ; //All recontructed particles\r
+ Int_t iFidArea = p->GetFiducialArea(); \r
+ if(iFidArea>0){\r
+ fhRecAllArea1->Fill(p->Pt() ) ; \r
+ if(iFidArea>1){\r
+ fhRecAllArea2->Fill(p->Pt() ) ; \r
+ if(iFidArea>2){\r
+ fhRecAllArea3->Fill(p->Pt() ) ; \r
+ }\r
+ }\r
+ }\r
+\r
+ if(fStack){\r
+ TParticle * prim = fStack->Particle(caloCluster->GetLabel()) ;\r
+ if(fDebug>2)\r
+ printf("Pdgcode = %d\n",prim->GetPdgCode());\r
+\r
+ if(prim->GetPdgCode()!=22){ //not photon\r
+//<--DP p->SetPhoton(kFALSE);\r
+ fhRecOther->Fill(p->Pt()); //not photons real spectra\r
+ for(Int_t iPID=0; iPID<4; iPID++){\r
+ if(p->IsPIDOK(iPID,22))\r
+ fhRecOtherPID[iPID]->Fill(p->Pt());\r
+ }\r
+ }\r
+ else{ //Primary photon (as in MC)\r
+//<--DP p->SetPhoton(kTRUE);\r
+ fhRecPhoton->Fill(p->Pt()); //Reconstructed with primary photon\r
+ for(Int_t iPID=0; iPID<4; iPID++){\r
+ if(p->IsPIDOK(iPID,22))\r
+ fhRecPhotonPID[iPID]->Fill(p->Pt());\r
+ }\r
+ Int_t pi0i=prim->GetFirstMother();\r
+ Int_t grandpaPDG=-1 ;\r
+ TParticle * pi0p = 0;\r
+ if(pi0i>=0){\r
+ pi0p=fStack->Particle(pi0i);\r
+ grandpaPDG=pi0p->GetPdgCode() ;\r
+ }\r
+ switch(grandpaPDG){\r
+ case 111: //Pi0 decay\r
+ //Primary decay photon (as in MC)\r
+ fhRecPhotPi0->Fill(p->Pt());\r
+ break ;\r
+ case 11:\r
+ case -11: //electron/positron conversion\r
+//<--DP p->SetConverted(1);\r
+ fhRecPhotConv->Fill(p->Pt()); //Reconstructed with photon from conversion primary\r
+ fhConversionRadius->Fill(prim->R());\r
+ break ;\r
+ case -2212:\r
+ case -2112: //antineutron & antiproton conversion\r
+ fhRecPhotHadron->Fill(p->Pt()); //Reconstructed with photon from antibaryon annihilation\r
+ fhInteractionRadius->Fill(prim->R());\r
+ break ;\r
+ \r
+ case 221: //eta decay\r
+ fhRecPhotEta->Fill(p->Pt());\r
+ break ; \r
+ \r
+ case 223: //omega meson decay\r
+ fhRecPhotOmega->Fill(p->Pt());\r
+ break ;\r
+ \r
+ case 331: //eta' decay\r
+ fhRecPhotEtapr->Fill(p->Pt());\r
+ break ;\r
+ \r
+ case -1: //direct photon or no primary\r
+ fhRecPhotDirect->Fill(p->Pt());\r
+ break ;\r
+ \r
+ default: \r
+ fhRecPhotOther->Fill(p->Pt());\r
+ break ;\r
+ } \r
+\r
+ //Now classify pi0 decay photon\r
+ if(grandpaPDG==111){\r
+//<--DP p->Pi0Decay(kTRUE); //Mark this photon as primary decayed\r
+//<--DP p->Pi0Id(pi0i); //remember id of the parent pi0\r
+\r
+ //Now check if second (partner) photon from pi0 decay hits PHOS or not\r
+ //i.e. both photons can be tagged or it's the systematic error\r
+//<--DP p->SetPartnerPt(0.); \r
+ Int_t indexdecay1,indexdecay2;\r
+\r
+ indexdecay1=pi0p->GetFirstDaughter();\r
+ indexdecay2=pi0p->GetLastDaughter();\r
+ Int_t indexdecay=-1;\r
+ if(fDebug>2)\r
+ printf("checking pi0 decay...index1=%d, index2=%d, index_pi0=%d, index_ph_prim=%d\n", indexdecay1,indexdecay2,pi0i,caloCluster->GetLabel());\r
+ \r
+ if(indexdecay1!=caloCluster->GetLabel()) \r
+ indexdecay=indexdecay1;\r
+ if(indexdecay2!=caloCluster->GetLabel()) \r
+ indexdecay=indexdecay2;\r
+ if(indexdecay==-1){\r
+ if(fDebug>2){\r
+ printf("Probably the other photon is not in the stack!\n");\r
+ printf("Number of daughters: %d\n",pi0p->GetNDaughters());\r
+ }\r
+ } \r
+ else{\r
+ TParticle *partner = fStack->Particle(indexdecay);\r
+ Int_t modulenum;\r
+ Double_t xtmp,ztmp;\r
+//<--DP p->SetPartnerPt(partner->Pt());\r
+ if(partner->GetPdgCode()==22){ \r
+ Bool_t isPartnerLost=kFALSE; //If partner is lost for some reason\r
+ if(partner->Energy()<fMinEnergyCut){ //energy is not enough to be registered by PHOS\r
+ if(fDebug>2)\r
+ printf("P_Reg, E=%f\n",partner->Energy());\r
+ fhPartnerMissedEmin->Fill(partner->Pt()); //Spectrum of missed partners\r
+ fhDecWMissedPartnerEmin->Fill(p->Pt()) ; //Spectrum of tagged with missed partner\r
+ isPartnerLost=kTRUE;\r
+ }\r
+ if(partner->GetNDaughters()!=0){ //this photon is converted before it is registered by some detector\r
+ if(fDebug>2)\r
+ printf("P_Conv, daughters=%d\n",partner->GetNDaughters());\r
+//<--DP p->SetConvertedPartner(1);\r
+ fhPartnerMissedConv->Fill(partner->Pt());\r
+ fhDecWMissedPartnerConv->Fill(p->Pt()) ; //Spectrum of tagged with missed partner\r
+ isPartnerLost=kTRUE;\r
+ }\r
+ if(!fgeom->ImpactOnEmc(partner,modulenum,ztmp,xtmp)){ //this photon cannot hit PHOS\r
+ if(fDebug>2)\r
+ printf("P_Geo\n");\r
+ fhPartnerMissedGeo->Fill(partner->Pt());\r
+ fhDecWMissedPartnerGeom0->Fill(p->Pt()) ; //Spectrum of tagged with missed partner\r
+ if(iFidArea>0){\r
+ fhDecWMissedPartnerGeom1->Fill(p->Pt()) ; //Spectrum of tagged with missed partner\r
+ if(iFidArea>1){\r
+ fhDecWMissedPartnerGeom2->Fill(p->Pt()) ; //Spectrum of tagged with missed partner\r
+ if(iFidArea>2){\r
+ fhDecWMissedPartnerGeom3->Fill(p->Pt()) ; //Spectrum of tagged with missed partner\r
+ }\r
+ }\r
+ }\r
+ isPartnerLost=kTRUE;\r
+ }\r
+ if(!isPartnerLost){\r
+// p->SetMCTagged(1); //set this photon as primary tagged\r
+ fhDecWMCPartner->Fill(p->Pt());\r
+ fhPartnerMCReg->Fill(partner->Pt());\r
+ if(fDebug>2){\r
+ 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", caloCluster->E(),partner->Pt(),partner->Energy(),partner->Px(),partner->Py(),partner->Pz(),partner->GetNDaughters());\r
+ printf(", module: %d, x tower: %f, z tower: %f\n", modulenum,xtmp,ztmp);\r
+ }\r
+ }\r
+ else{\r
+ fhDecWMissedPartnerAll->Fill(p->Pt());\r
+ }\r
+ }//Partner - photon\r
+ else{//partner not photon\r
+ fhDecWMissedPartnerNotPhoton->Fill(p->Pt()); \r
+ }\r
+ }\r
+ }\r
+ }\r
+ }\r
+ } //PHOS clusters\r
+ \r
+ if(fDebug>1) \r
+ printf("number of clusters: %d\n",inList);\r
+\r
+ //Invariant Mass analysis\r
+ for(Int_t phosPhoton1 = 0 ; phosPhoton1 < inList-1 ; phosPhoton1++) {\r
+ AliAODPWG4Particle * p1 = static_cast<AliAODPWG4Particle*>(fCaloPhotonsArr->At(phosPhoton1)); \r
+ for(Int_t phosPhoton2 = phosPhoton1 + 1 ; phosPhoton2 < inList ; phosPhoton2++) {\r
+ AliAODPWG4Particle * p2 = static_cast<AliAODPWG4Particle*>(fCaloPhotonsArr->At(phosPhoton2));\r
+\r
+ Double_t invMass = p1->GetPairMass(p2);\r
+ fhInvMassReal->Fill(invMass,p1->Pt());\r
+ fhInvMassReal->Fill(invMass,p2->Pt());\r
+ if(fDebug>2)\r
+ printf("Pair i=%d,j=%d, M=%f\n",phosPhoton1,phosPhoton2,invMass);\r
+\r
+ Bool_t makePi01=IsInPi0Band(invMass,p1->Pt());\r
+ Bool_t makePi02=IsInPi0Band(invMass,p2->Pt());\r
+\r
+\r
+ if(makePi01 && p1->IsTagged()){//Multiple tagging\r
+ fhTaggedMult->Fill(p1->Pt());\r
+ } \r
+ if(makePi01 && !p1->IsTagged()){//Each photon should enter histogram once, even if tagged several times\r
+ fhTaggedAll->Fill(p1->Pt());\r
+ Int_t iFidArea = p1->GetFiducialArea(); \r
+ if(iFidArea>0){\r
+ fhTaggedArea1->Fill(p1->Pt() ) ; \r
+ if(iFidArea>1){\r
+ fhTaggedArea2->Fill(p1->Pt() ) ; \r
+ if(iFidArea>2){\r
+ fhTaggedArea3->Fill(p1->Pt() ) ; \r
+ }\r
+ }\r
+ }\r
+\r
+ for(Int_t iPID=0; iPID<4; iPID++){\r
+ if(p1->IsPIDOK(iPID,22))\r
+ fhTaggedPID[iPID]->Fill(p1->Pt());\r
+ }\r
+\r
+ p1->SetTagged(kTRUE) ;\r
+ } \r
+ if(makePi02 && p2->IsTagged()){//Multiple tagging\r
+ fhTaggedMult->Fill(p2->Pt());\r
+ } \r
+ if(makePi02 && !p2->IsTagged()){//How should be account for multiply tagged photons?\r
+ fhTaggedAll->Fill(p2->Pt());\r
+ p2->SetTagged(kTRUE) ;\r
+ }\r
+ \r
+ //Now get use MC information\r
+ //First chesk if this is true pi0 pair\r
+ if(IsSamePi0(p1,p2)){ //Correctly tagged - from the same pi0\r
+// p1->SetTrueTagged(1);\r
+// p2->SetTrueTagged(1);\r
+ if(makePi01)//Correctly tagged photons\r
+ fhTaggedMCTrue->Fill(p1->Pt());\r
+ else{ //Decay pair missed tagging \r
+ fhMCMissedTagging->Fill(p1->Pt());\r
+ fhMCMissedTaggingMass->Fill(invMass,p1->Pt()) ;\r
+ //Clussify why missed tagging (todo)\r
+ //Converted\r
+ //Partner not a photon\r
+ //Tagged not a photon\r
+ //Just wrong inv.mass \r
+ } \r
+ if(makePi02)\r
+ fhTaggedMCTrue->Fill(p2->Pt());\r
+ else{ \r
+ fhMCMissedTagging->Fill(p2->Pt());\r
+ fhMCMissedTaggingMass->Fill(invMass,p2->Pt()) ;\r
+ //Clussify why missed tagging (todo)\r
+ //Converted\r
+ //Partner not a photon\r
+ //Tagged not a photon\r
+ //Just wrong inv.mass \r
+ } \r
+ }\r
+ else{//Fake tagged - not from the same pi0\r
+ if(makePi01)//Fake pair\r
+ fhMCFakeTagged->Fill(p1->Pt());\r
+ if(makePi02)\r
+ fhMCFakeTagged->Fill(p2->Pt());\r
+ }\r
+ }\r
+ }\r
+\r
+ //Fill Mixed InvMass distributions:\r
+ TIter nextEv(fEventList) ;\r
+ while(TClonesArray * event2 = static_cast<TClonesArray*>(nextEv())){\r
+ Int_t nPhotons2 = event2->GetEntriesFast() ;\r
+ for(Int_t i=0; i < inList ; i++){\r
+ AliAODPWG4Particle * p1 = static_cast<AliAODPWG4Particle*>(fCaloPhotonsArr->At(i)) ;\r
+ for(Int_t j=0; j < nPhotons2 ; j++){\r
+ AliAODPWG4Particle * p2 = static_cast<AliAODPWG4Particle*>(event2->At(j)) ;\r
+ Double_t invMass = p1->GetPairMass(p2);\r
+ fhInvMassMixed->Fill(invMass,p1->Pt());\r
+ fhInvMassMixed->Fill(invMass,p2->Pt());\r
+ }\r
+ }\r
+ }\r
+\r
+/* if(inList==1){ //We have only one photon in PHOS, but still there can be missing partner for it\r
+ AliAODPWG4Particle * p = static_cast<AliAODPWG4Particle*>(fCaloPhotonsArr->At(0));\r
+// if(p->IsMCTagged()){\r
+ //should have pair but partner was not registered in PHOS\r
+// double pairpt=p->GetPairPt();\r
+// if(fDebug>2)\r
+// printf("pair not found for phot: E=%f, pair Pt=%f\n",p->Energy(),pairpt);\r
+// fhMCPartnerMissed->Fill(pairpt);\r
+// fhstra->Fill(p->Pt());\r
+// }\r
+ }\r
+*/\r
+\r
+ //Remove old events\r
+ fEventList->AddFirst(fCaloPhotonsArr);\r
+ if(fEventList->GetSize() > 10){\r
+ TClonesArray *tmp = static_cast <TClonesArray*> (fEventList->Last());\r
+ fEventList->Remove(tmp);\r
+ delete tmp;\r
+ }\r
+\r
+ PostData(1, fOutputList);\r
+}\r
+\r
+\r
+//______________________________________________________________________________\r
+void AliAnalysisTaskTaggedPhotons::Init()\r
+{\r
+ // Intialisation of parameters\r
+ AliInfo("Doing initialisation") ; \r
+ SetPhotonId(0.9) ; \r
+ SetMinEnergyCut(0.4);\r
+ SetPi0MeanParameters(0.1377,-0.002566,0.001216,-0.0001256);\r
+ SetPi0SigmaParameters(0.004508,0.005497,0.00000006);\r
+}\r
+\r
+//______________________________________________________________________________\r
+void AliAnalysisTaskTaggedPhotons::Terminate(Option_t *)\r
+{\r
+ // Processing when the event loop is ended\r
+\r
+ //Write everything to the file\r
+ char outname[55];\r
+ if(fPHOS)\r
+ sprintf(outname,"Tagging_PHOS.root") ;\r
+ else \r
+ sprintf(outname,"Tagging_EMCAL.root") ;\r
+ TFile *outfile = new TFile (outname,"recreate");\r
+\r
+fhRecAll->Write();\r
+fhRecAllArea1->Write();\r
+fhRecAllArea2->Write();\r
+fhRecAllArea3->Write();\r
+fhRecPhoton->Write();\r
+fhRecOther->Write();\r
+fhRecPhotonPID[0]->Write();\r
+fhRecPhotonPID[1]->Write();\r
+fhRecPhotonPID[2]->Write();\r
+fhRecPhotonPID[3]->Write();\r
+fhRecOtherPID[0]->Write();\r
+fhRecOtherPID[1]->Write();\r
+fhRecOtherPID[2]->Write();\r
+fhRecOtherPID[3]->Write();\r
+fhRecPhotPi0->Write();\r
+fhRecPhotEta->Write();\r
+fhRecPhotOmega->Write();\r
+fhRecPhotEtapr->Write();\r
+fhRecPhotConv->Write();\r
+fhRecPhotHadron->Write();\r
+fhRecPhotDirect->Write();\r
+fhRecPhotOther->Write();\r
+fhDecWMCPartner->Write();\r
+fhDecWMissedPartnerNotPhoton->Write();\r
+fhDecWMissedPartnerAll->Write();\r
+fhDecWMissedPartnerEmin->Write();\r
+fhDecWMissedPartnerConv->Write();\r
+fhDecWMissedPartnerGeom0->Write();\r
+fhDecWMissedPartnerGeom1->Write();\r
+fhDecWMissedPartnerGeom2->Write();\r
+fhDecWMissedPartnerGeom3->Write();\r
+fhPartnerMCReg->Write();\r
+fhPartnerMissedEmin->Write();\r
+fhPartnerMissedConv->Write();\r
+fhPartnerMissedGeo->Write();\r
+fhTaggedAll->Write();\r
+fhTaggedArea1->Write();\r
+fhTaggedArea2->Write();\r
+fhTaggedArea3->Write();\r
+\r
+fhTaggedPID[0]->Write();\r
+fhTaggedPID[1]->Write();\r
+fhTaggedPID[2]->Write();\r
+fhTaggedPID[3]->Write();\r
+fhTaggedMult->Write();\r
+fhTaggedMCTrue->Write();\r
+fhMCMissedTagging->Write();\r
+fhMCFakeTagged->Write();\r
+fhInvMassReal->Write();\r
+fhInvMassMixed->Write();\r
+fhMCMissedTaggingMass->Write();\r
+fhConversionRadius->Write();\r
+fhInteractionRadius->Write();\r
+fhEvents->Write();\r
+\r
+/*\r
+ fhAllPhotons->Write() ;\r
+ fhNotPhotons->Write() ;\r
+ fhAllPhotonsPrimary->Write() ;\r
+ fhNotPhotonsPrimary->Write() ;\r
+ fhfakeNotPhotons->Write() ;\r
+\r
+ fhTaggedPhotons->Write();\r
+ fhfakeTaggedPhotons->Write();\r
+ fhDecayNotTaggedPhotons->Write();\r
+ fhstrangeNotTaggedPhotons->Write();\r
+ fhstrangeNotTaggedPhotonsPair->Write();\r
+ fhstrangeNotTaggedPhotonsRegCut->Write();\r
+ fhstrangeNotTaggedPhotonsPairRegCut->Write();\r
+\r
+ fhPi0DecayPhotonsPrimary->Write();\r
+ fhEtaDecayPhotonsPrimary->Write();\r
+ fhOmegaDecayPhotonsPrimary->Write();\r
+ fhEtaSDecayPhotonsPrimary->Write();\r
+ fhOtherDecayPhotonsPrimary->Write();\r
+ fhDecayPhotonsPrimary->Write();\r
+ fhConvertedPhotonsPrimary->Write();\r
+ fhConvertedPhotonsPrimaryHadronsDecays->Write();\r
+ fhCoordsConvertion->Write();\r
+ fhCoordsConvertion2->Write();\r
+\r
+ fhPHOSInvariantMassReal->Write();\r
+ fhPHOSInvariantMassMixed->Write();\r
+\r
+ fhPHOSPi0->Write();\r
+\r
+ fhPi0DecayPhotonsGeomfake->Write();\r
+ fhPi0DecayPhotonsTaggedPrimary->Write();\r
+ fhPi0DecayPhotonsTaggedPrimaryPair->Write();\r
+ fhPi0DecayPhotonsTaggedPrimaryRegCut->Write();\r
+ fhPi0DecayPhotonsTaggedPrimaryPairRegCut->Write();\r
+\r
+ fhPi0DecayPhotonsBigDecay->Write();\r
+ fhPi0DecayPhotonsPConv->Write();\r
+ fhPi0DecayPhotonsPGeo->Write();\r
+ fhPi0DecayPhotonsPReg->Write();\r
+\r
+ fhfakeTaggedPhotonsConv->Write();\r
+ fhfakeTaggedPhotonsPID->Write();\r
+\r
+ fhTrackRefCoords->Write();\r
+ fhEvents->Write();\r
+*/\r
+\r
+outfile->Close();\r
+\r
+}\r
+//______________________________________________________________________________\r
+Bool_t AliAnalysisTaskTaggedPhotons::IsInPi0Band(Double_t m, Double_t pt)\r
+{\r
+ //Parameterization of the pi0 peak region\r
+ Double_t mpi0_mean = fPi0Mean_p0 + fPi0Mean_p1 * pt + fPi0Mean_p2 * pt*pt + fPi0Mean_p3 * pt*pt*pt;\r
+ Double_t mpi0_sigma = TMath::Sqrt(fPi0Sigma_p0 * fPi0Sigma_p0 / pt + fPi0Sigma_p1 * fPi0Sigma_p1 + fPi0Sigma_p2 * fPi0Sigma_p2 / pt / pt);\r
+ \r
+ return (m>mpi0_mean-2*mpi0_sigma && m>mpi0_mean+2*mpi0_sigma) ;\r
+}\r
+//______________________________________________________________________________\r
+Bool_t AliAnalysisTaskTaggedPhotons::IsSamePi0(AliAODPWG4Particle *p1, AliAODPWG4Particle *p2){\r
+ //Looks through parents and finds if there was commont pi0 among ancestors\r
+\r
+ if(!fStack)\r
+ return kFALSE ; //can not say anything\r
+\r
+ Int_t prim1 = p1->GetLabel();\r
+ while(prim1!=-1){ \r
+ Int_t prim2 = p2->GetLabel();\r
+ while(prim2!=-1){ \r
+ if(prim1==prim2){\r
+ if(fStack->Particle(prim1)->GetPdgCode()==111)\r
+ return kTRUE ;\r
+ else\r
+ return kFALSE ;\r
+ }\r
+ prim2=fStack->Particle(prim2)->GetFirstMother() ;\r
+ }\r
+ prim1=fStack->Particle(prim1)->GetFirstMother() ;\r
+ }\r
+ return kFALSE ;\r
+}\r
+//______________________________________________________________________________\r
+Int_t AliAnalysisTaskTaggedPhotons::GetFiducialArea(Float_t * pos){\r
+ //calculates in which kind of fiducial area photon hit\r
+ if(fPHOS){\r
+ Double_t phi=TMath::ATan2(pos[1],pos[0]) ;\r
+ Double_t z=pos[2] ;
+ while(phi>TMath::TwoPi())phi-=TMath::TwoPi() ;\r
+ while(phi<0.)phi+=TMath::TwoPi() ;\r
+ //From active PHOS areat remove bands in 10 cm\r
+ const Double_t kphi=TMath::ATan(10./460.) ; //angular band width\r
+ Double_t dzMax=TMath::Ceil((fZmax-z)/10.) ;\r
+ Double_t dzMin=TMath::Ceil((z-fZmin)/10.) ;\r
+ Double_t dphiMax=TMath::Ceil((fPhimax-phi)/kphi);\r
+ Double_t dphiMin=TMath::Ceil((phi-fPhimin)/kphi);\r
+ return (Int_t)TMath::Min(TMath::Min(dzMax,dzMin),TMath::Min(dphiMax,dphiMin)); \r
+ }\r
+ else{//EMCAL\r
+ //For the moment whole EMCAL is considered as area 1\r
+ return 1 ;\r
+ }\r
+}\r