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