#include "AliAODEvent.h"
#include "AliAODEvent.h"
#include "AliVCluster.h"
-#include "AliAODPWG4Particle.h"
+#include "AliCaloPhoton.h"
#include "AliAODMCParticle.h"
#include "AliAnalysisManager.h"
#include "AliLog.h"
#include "AliEMCALGeometry.h"
#include "AliAnalysisUtils.h"
#include "AliOADBContainer.h"
+#include "AliAODMCHeader.h"
ClassImp(AliAnalysisTaskTaggedPhotons)
fZmin(0.),
fPhimax(0.),
fPhimin(0.),
+ fMinBCDistance(0.),
fCentrality(0.),
fCentBin(0),
- fIsMB(0)
+ fIsMB(0),
+ fIsMC(0),
+ fIsFastMC(0)
{
//Deafult constructor
//no memory allocations
fZmin(60.),
fPhimax(250.),
fPhimin(320.),
+ fMinBCDistance(0.),
fCentrality(0.),
fCentBin(0),
- fIsMB(0)
+ fIsMB(0),
+ fIsMC(0),
+ fIsFastMC(0)
{
// Constructor.
fZmin(60.),
fPhimax(250.),
fPhimin(320.),
+ fMinBCDistance(0.),
fCentrality(0.),
fCentBin(0),
- fIsMB(0)
+ fIsMB(0),
+ fIsMC(0),
+ fIsFastMC(0)
{
// cpy ctor
fZmax=ap.fZmax ;
fOutputContainer->Add(new TH1F("hCentrality","Ccentrality",100,0.,100.));
fOutputContainer->Add(new TH2F("hPHOSCentrality","PHOS vs centrality",100,0.,100.,100,0.,100.));
fOutputContainer->Add(new TH2F("hTOF","cluster TOF",200,0.,20.,300,-3.e-6,6.e-6));
- fOutputContainer->Add(new TH2F("hCenTrack","Centrality vs tracks", 100,0.,100.,100,0.,1500.)) ;
fOutputContainer->Add(new TH2F("hCluNXZM1","Clu (X,Z), M1" ,64,0.5,64.5, 56,0.5,56.5));
fOutputContainer->Add(new TH2F("hCluNXZM2","Clu (X,Z), M2" ,64,0.5,64.5, 56,0.5,56.5));
snprintf(cPID[3],5,"Both");
- const Int_t nPt=400 ;
- const Double_t ptMax=40. ;
+ const Int_t nPt=500 ;
+ const Double_t ptMax=50. ;
const Int_t nM=400 ;
const Double_t mMax=1. ;
+ //QA histograms
+ for(Int_t iPID=0; iPID<4; iPID++){
+ fOutputContainer->Add(new TH1F(Form("hPhotM1_%s",cPID[iPID]),"Spectrum of all reconstructed particles, M1",nPt,0.,ptMax)) ;
+ fOutputContainer->Add(new TH1F(Form("hPhotM2_%s",cPID[iPID]),"Spectrum of all reconstructed particles, M2",nPt,0.,ptMax)) ;
+ fOutputContainer->Add(new TH1F(Form("hPhotM3_%s",cPID[iPID]),"Spectrum of all reconstructed particles, M3",nPt,0.,ptMax)) ;
+ fOutputContainer->Add(new TH1F(Form("hPhotM1A2_%s",cPID[iPID]),"Spectrum of all reconstructed particles, M1",nPt,0.,ptMax)) ;
+ fOutputContainer->Add(new TH1F(Form("hPhotM2A2_%s",cPID[iPID]),"Spectrum of all reconstructed particles, M2",nPt,0.,ptMax)) ;
+ fOutputContainer->Add(new TH1F(Form("hPhotM3A2_%s",cPID[iPID]),"Spectrum of all reconstructed particles, M3",nPt,0.,ptMax)) ;
+ fOutputContainer->Add(new TH1F(Form("hPhotM1A3_%s",cPID[iPID]),"Spectrum of all reconstructed particles, M1",nPt,0.,ptMax)) ;
+ fOutputContainer->Add(new TH1F(Form("hPhotM2A3_%s",cPID[iPID]),"Spectrum of all reconstructed particles, M2",nPt,0.,ptMax)) ;
+ fOutputContainer->Add(new TH1F(Form("hPhotM3A3_%s",cPID[iPID]),"Spectrum of all reconstructed particles, M3",nPt,0.,ptMax)) ;
+ }
+
+
const Int_t nCenBin=5;
for(Int_t cen=0; cen<nCenBin; cen++){
fOutputContainer->Add(new TH1F(Form("hPhot_Area2_%s_cent%d",cPID[iPID],cen),"Spectrum of all reconstructed particles",nPt,0.,ptMax)) ;
fOutputContainer->Add(new TH1F(Form("hPhot_Area3_%s_cent%d",cPID[iPID],cen),"Spectrum of all reconstructed particles",nPt,0.,ptMax)) ;
- fOutputContainer->Add(new TH1F(Form("hPhot_nStrictTagged_Area1_%s_cent%d",cPID[iPID],cen),"Spectrum of all reconstructed particles, no PID",nPt,0.,ptMax)) ;
-
for(Int_t itag=0; itag<18; itag++){
- fOutputContainer->Add(new TH1F(Form("hPhot_nTagged%d_Area1_%s_cent%d",itag,cPID[iPID],cen),"Spectrum of all reconstructed particles, no PID",nPt,0.,ptMax)) ;
- fOutputContainer->Add(new TH1F(Form("hPhot_nTagged%d_Area2_%s_cent%d",itag,cPID[iPID],cen),"Spectrum of all reconstructed particles, no PID",nPt,0.,ptMax)) ;
- fOutputContainer->Add(new TH1F(Form("hPhot_nTagged%d_Area3_%s_cent%d",itag,cPID[iPID],cen),"Spectrum of all reconstructed particles, no PID",nPt,0.,ptMax)) ;
+ fOutputContainer->Add(new TH1F(Form("hPhot_TaggedMult%d_Area1_%s_cent%d",itag,cPID[iPID],cen),"Spectrum of multiply tagged photons",nPt,0.,ptMax)) ;
+ fOutputContainer->Add(new TH1F(Form("hPhot_TaggedMult%d_Area2_%s_cent%d",itag,cPID[iPID],cen),"Spectrum of multiply tagged photons",nPt,0.,ptMax)) ;
+ fOutputContainer->Add(new TH1F(Form("hPhot_TaggedMult%d_Area3_%s_cent%d",itag,cPID[iPID],cen),"Spectrum of multiply tagged photons",nPt,0.,ptMax)) ;
+ fOutputContainer->Add(new TH1F(Form("hPhot_Tagged%d_Area1_%s_cent%d",itag,cPID[iPID],cen),"Spectrum of all reconstructed particles, no PID",nPt,0.,ptMax)) ;
+ fOutputContainer->Add(new TH1F(Form("hPhot_Tagged%d_Area2_%s_cent%d",itag,cPID[iPID],cen),"Spectrum of all reconstructed particles, no PID",nPt,0.,ptMax)) ;
+ fOutputContainer->Add(new TH1F(Form("hPhot_Tagged%d_Area3_%s_cent%d",itag,cPID[iPID],cen),"Spectrum of all reconstructed particles, no PID",nPt,0.,ptMax)) ;
+
+ fOutputContainer->Add(new TH1F(Form("hPhot_TrueTagged%d_%s_cent%d",itag,cPID[iPID],cen),"Spectrum of all tagged particles",nPt,0.,ptMax)) ;
}
for(Int_t kind=1; kind<33; kind*=2){
fOutputContainer->Add(new TH1F(Form("hPhot_Isolation%d_%s_cent%d",kind,cPID[iPID],cen),"Spectrum of all reconstructed particles, no PID",nPt,0.,ptMax)) ;
fOutputContainer->Add(new TH1F(Form("hPhot_Isolation%d_Area1_%s_cent%d",kind,cPID[iPID],cen),"Spectrum of all reconstructed particles, no PID",nPt,0.,ptMax)) ;
- fOutputContainer->Add(new TH1F(Form("hPhot_nStrictTagged_Isolation%d_Area1_%s_cent%d",kind,cPID[iPID],cen),"Spectrum of all reconstructed particles, no PID",nPt,0.,ptMax)) ;
fOutputContainer->Add(new TH1F(Form("hPhot_nTagged_Isolation%d_Area1_%s_cent%d",kind,cPID[iPID],cen),"Spectrum of all reconstructed particles, no PID",nPt,0.,ptMax)) ;
}
}
-
+ for(Int_t kind=1; kind<33; kind*=2){
+ fOutputContainer->Add(new TH1F(Form("hPi_Isolation%d_cent%d",kind,cen),"Spectrum of all reconstructed particles, no PID",nPt,0.,ptMax)) ;
+ }
+
- fOutputContainer->Add(new TH1F(Form("hTaggedMult_cent%d",cen),"Spectrum of multiply tagged photons",nPt,0.,ptMax)) ;
//Invariant mass distributions for fake corrections
fOutputContainer->Add(new TH2F(Form("QA_Cone1_Tracks_cent%d",cen),"Two-photon inv. mass vs first photon pt",50,0.,50.,200,0.,100.)) ;
fOutputContainer->Add(new TH2F(Form("QA_Cone2_Tracks_cent%d",cen),"Two-photon inv. mass vs first photon pt",50,0.,50.,200,0.,100.)) ;
fOutputContainer->Add(new TH2F(Form("QA_Cone3_Tracks_cent%d",cen),"Two-photon inv. mass vs first photon pt",50,0.,50.,200,0.,100.)) ;
+ fOutputContainer->Add(new TH2F(Form("QA_PCone1_Tracks_cent%d",cen),"Two-photon inv. mass vs first photon pt",50,0.,50.,200,0.,100.)) ;
+ fOutputContainer->Add(new TH2F(Form("QA_PCone2_Tracks_cent%d",cen),"Two-photon inv. mass vs first photon pt",50,0.,50.,200,0.,100.)) ;
+ fOutputContainer->Add(new TH2F(Form("QA_PCone3_Tracks_cent%d",cen),"Two-photon inv. mass vs first photon pt",50,0.,50.,200,0.,100.)) ;
+ fOutputContainer->Add(new TH2F(Form("QA_Pi0Cone1_Tracks_cent%d",cen),"Two-photon inv. mass vs first photon pt",50,0.,50.,200,0.,100.)) ;
+ fOutputContainer->Add(new TH2F(Form("QA_Pi0Cone2_Tracks_cent%d",cen),"Two-photon inv. mass vs first photon pt",50,0.,50.,200,0.,100.)) ;
+ fOutputContainer->Add(new TH2F(Form("QA_Pi0Cone3_Tracks_cent%d",cen),"Two-photon inv. mass vs first photon pt",50,0.,50.,200,0.,100.)) ;
+ fOutputContainer->Add(new TH2F(Form("QA_Pi0PCone1_Tracks_cent%d",cen),"Two-photon inv. mass vs first photon pt",50,0.,50.,200,0.,100.)) ;
+ fOutputContainer->Add(new TH2F(Form("QA_Pi0PCone2_Tracks_cent%d",cen),"Two-photon inv. mass vs first photon pt",50,0.,50.,200,0.,100.)) ;
+ fOutputContainer->Add(new TH2F(Form("QA_Pi0PCone3_Tracks_cent%d",cen),"Two-photon inv. mass vs first photon pt",50,0.,50.,200,0.,100.)) ;
}//centrality
//MC
- char partName[4][10] ;
+ char partName[15][10] ;
snprintf(partName[0],10,"gamma") ;
snprintf(partName[1],10,"pi0");
snprintf(partName[2],10,"eta") ;
snprintf(partName[3],10,"omega");
+ snprintf(partName[4],10,"K0s");
+ snprintf(partName[5],10,"Kpm");
+ snprintf(partName[6],10,"pipm");
+ snprintf(partName[7],10,"n");
+ snprintf(partName[8],10,"nbar");
+ snprintf(partName[9],10,"p");
+ snprintf(partName[10],10,"pbar");
+
- if(AliAnalysisManager::GetAnalysisManager()){
- AliMCEventHandler* mctruth = (AliMCEventHandler*)((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler());
- if(mctruth){
+ if(fIsMC){
fOutputContainer->Add(new TH1F("hMCConversionRadius","Clusters without label",600,0.,600.)) ;
fOutputContainer->Add(new TH2F("hMCRecPi0Vtx","Secondary pi0s",100,0.,10.,600,0.,600.)) ;
fOutputContainer->Add(new TH2F("hMCRecK0lVtx","Secondary K0l",100,0.,10.,600,0.,600.)) ;
fOutputContainer->Add(new TH2F("hMCGammaPi0MisConvR","Converted photons",400,0.,40.,600,0.,600.)) ;
- for(Int_t cen=0; cen<nCenBin; cen++){
- for(Int_t ipart=0; ipart<4; ipart++){
- fOutputContainer->Add(new TH2F(Form("hMCpi0_ptrap_cen%d",cen),"Spectrum of primary photons",100,0.,10.,100,-2.,2.)) ;
- fOutputContainer->Add(new TH2F(Form("hMCpi0_ptphi_cen%d",cen),"Spectrum of primary photons",100,0.,10.,100,0.,TMath::TwoPi())) ;
- fOutputContainer->Add(new TH2F(Form("hMC_%s_vertex",partName[ipart]),"vertex",nPt,0.,ptMax,150,0.,600.)) ;
- fOutputContainer->Add(new TH1F(Form("hMC_all_%s",partName[ipart]),"Spectum (full rapifity)",nPt,0.,ptMax)) ;
- fOutputContainer->Add(new TH1F(Form("hMC_unitEta_%s",partName[ipart]),"Spectum, |y|<0.15",nPt,0.,ptMax)) ;
- fOutputContainer->Add(new TH1F(Form("hMC_rap_%s",partName[ipart]),"Rapidity",100,-5.,5.)) ;
- fOutputContainer->Add(new TH1F(Form("hMC_phi_%s",partName[ipart]),"Azimuthal angle",100,0.,TMath::TwoPi())) ;
+ for(Int_t ipart=0; ipart<11; ipart++){
+ fOutputContainer->Add(new TH2F(Form("hMC%s_ptrap",partName[ipart]),"Spectrum of primary photons",100,0.,10.,200,-1.,1.)) ;
+ fOutputContainer->Add(new TH2F(Form("hMC%s_ptphi",partName[ipart]),"Spectrum of primary photons",100,0.,10.,100,0.,TMath::TwoPi())) ;
+ fOutputContainer->Add(new TH2F(Form("hMC_%s_vertex",partName[ipart]),"vertex",nPt,0.,ptMax,150,0.,600.)) ;
+ for(Int_t cen=0; cen<nCenBin; cen++){
+ fOutputContainer->Add(new TH1F(Form("hMC_all_%s_cent%d",partName[ipart],cen),"Spectum (full rapifity)",nPt,0.,ptMax)) ;
+ fOutputContainer->Add(new TH1F(Form("hMC_unitEta_%s_cent%d",partName[ipart],cen),"Spectum, |y|<0.15",nPt,0.,ptMax)) ;
+ fOutputContainer->Add(new TH1F(Form("hMC_rap_%s_cent%d",partName[ipart],cen),"Rapidity",100,-5.,5.)) ;
+ fOutputContainer->Add(new TH1F(Form("hMC_phi_%s_cent%d",partName[ipart],cen),"Azimuthal angle",100,0.,TMath::TwoPi())) ;
}
-
+ }
+ for(Int_t cen=0; cen<nCenBin; cen++){
+ fOutputContainer->Add(new TH1F(Form("hMC_unitEta_GammaPi0_cent%d",cen),"Spectum, |y|<0.15",nPt,0.,ptMax)) ;
+ }
+
+ for(Int_t cen=0; cen<nCenBin; cen++){
fOutputContainer->Add(new TH2F(Form("LabelsNPrim_cent%d",cen),"vertex",nPt,0.,ptMax,20,0.,20.)) ;
fOutputContainer->Add(new TH1F(Form("LabelsGamma_cent%d",cen),"vertex",nPt,0.,ptMax)) ;
fOutputContainer->Add(new TH2F(Form("LabelsGammaE_cent%d",cen),"vertex",nPt,0.,ptMax,100,0.,2.)) ;
+
+
+ //Sort registered particles spectra according MC information
+ for(Int_t iPID=0; iPID<4; iPID++){
+ fOutputContainer->Add(new TH1F(Form("hMCRecPhoton_%s_cent%d",cPID[iPID],cen),"Spectrum of rec. photons", nPt,0.,ptMax )) ;
+ fOutputContainer->Add(new TH1F(Form("hMCRecPhoton_Area1_%s_cent%d",cPID[iPID],cen),"Spectrum of rec. photons", nPt,0.,ptMax )) ;
+ fOutputContainer->Add(new TH1F(Form("hMCRecPhoton_Area2_%s_cent%d",cPID[iPID],cen),"Spectrum of rec. photons", nPt,0.,ptMax )) ;
+ fOutputContainer->Add(new TH1F(Form("hMCRecPhoton_Area3_%s_cent%d",cPID[iPID],cen),"Spectrum of rec. photons", nPt,0.,ptMax )) ;
+ fOutputContainer->Add(new TH1F(Form("hMCRecE_%s_cent%d",cPID[iPID],cen),"Spectrum of rec. electrons", nPt,0.,ptMax )) ;
+ fOutputContainer->Add(new TH1F(Form("hMCRecPbar_%s_cent%d",cPID[iPID],cen),"Spectrum of rec. electrons", nPt,0.,ptMax )) ;
+ fOutputContainer->Add(new TH1F(Form("hMCRecNbar_%s_cent%d",cPID[iPID],cen),"Spectrum of rec. electrons", nPt,0.,ptMax )) ;
+ fOutputContainer->Add(new TH1F(Form("hMCRecP_%s_cent%d",cPID[iPID],cen),"Spectrum of rec. electrons", nPt,0.,ptMax )) ;
+ fOutputContainer->Add(new TH1F(Form("hMCRecPipm_%s_cent%d",cPID[iPID],cen),"Spectrum of rec. electrons", nPt,0.,ptMax )) ;
+ fOutputContainer->Add(new TH1F(Form("hMCRecKpm_%s_cent%d",cPID[iPID],cen),"Spectrum of rec. electrons", nPt,0.,ptMax )) ;
+ fOutputContainer->Add(new TH1F(Form("hMCRecN_%s_cent%d",cPID[iPID],cen),"Spectrum of rec. electrons", nPt,0.,ptMax )) ;
+ fOutputContainer->Add(new TH1F(Form("hMCRecCharg_%s_cent%d",cPID[iPID],cen),"Spectrum of rec. electrons", nPt,0.,ptMax )) ;
+ fOutputContainer->Add(new TH1F(Form("hMCRecNeutral_%s_cent%d",cPID[iPID],cen),"Spectrum of rec. electrons", nPt,0.,ptMax )) ;
+ fOutputContainer->Add(new TH1F(Form("hMCRecK0s_%s_cent%d",cPID[iPID],cen),"Spectrum of rec. electrons", nPt,0.,ptMax )) ;
+ fOutputContainer->Add(new TH1F(Form("hMCRecNoPRim_%s_cent%d",cPID[iPID],cen),"Spectrum of rec. electrons", nPt,0.,ptMax )) ;
+
+ //Decay photons
+ fOutputContainer->Add(new TH1F(Form("hMCRecPhotPi0_%s_cent%d",cPID[iPID],cen),"Spectrum of rec. photons", nPt,0.,ptMax )) ;
+ fOutputContainer->Add(new TH1F(Form("hMCRecPhotEta_%s_cent%d",cPID[iPID],cen),"Spectrum of rec. photons", nPt,0.,ptMax )) ;
+ fOutputContainer->Add(new TH1F(Form("hMCRecPhotOmega_%s_cent%d",cPID[iPID],cen),"Spectrum of rec. photons", nPt,0.,ptMax )) ;
+ fOutputContainer->Add(new TH1F(Form("hMCRecPhotOther_%s_cent%d",cPID[iPID],cen),"Spectrum of rec. photons", nPt,0.,ptMax )) ;
+ fOutputContainer->Add(new TH1F(Form("hMCRecPhotNoPrim_%s_cent%d",cPID[iPID],cen),"Spectrum of rec. photons", nPt,0.,ptMax )) ;
+
+
+ //MC tagging: reasons of partner loss etc.
+ fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnStack_%s_cent%d",cPID[iPID],cen),"Decay photons with partner not in Stack", nPt,0.,ptMax )) ;
+ for(Int_t iType=0; iType<9; iType++)
+ fOutputContainer->Add(new TH1F(Form("hMCDecWithFoundPartnType%d_%s_cent%d",iType,cPID[iPID],cen),"Decay photon with found partner", nPt,0.,ptMax )) ;
+ fOutputContainer->Add(new TH1F(Form("hMCDecWithWrongMass_%s_cent%d",cPID[iPID],cen),"Decay photon with wrong mass", nPt,0.,ptMax )) ;
+ fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnAccept_%s_cent%d",cPID[iPID],cen),"Decay photon with parttner not in PHOS acc", nPt,0.,ptMax )) ;
+ fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnAcceptFA1_%s_cent%d",cPID[iPID],cen),"Decay photons with partner missed due geometry Fid. area. 1", nPt,0.,ptMax )) ;
+ fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnAcceptFA2_%s_cent%d",cPID[iPID],cen),"Decay photons with partner missed due geometry Fid. area. 2", nPt,0.,ptMax)) ;
+ fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnAcceptFA3_%s_cent%d",cPID[iPID],cen),"Decay photons with partner missed due geometry Fid. area. 3", nPt,0.,ptMax )) ;
+ fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnConv_%s_cent%d",cPID[iPID],cen),"Decay photons with partner missed due to conversion", nPt,0.,ptMax )) ;
+ fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnEmin_%s_cent%d",cPID[iPID],cen),"Decay photons with partner missed due to low energy", nPt,0.,ptMax )) ;
+ fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnOther_%s_cent%d",cPID[iPID],cen),"Decay photons with partner missed due unknown reason", nPt,0.,ptMax )) ;
+ fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnAll_%s_cent%d",cPID[iPID],cen),"Decay photons with partner missed due to any reason", nPt,0.,ptMax )) ;
+ fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnNPhot_%s_cent%d",cPID[iPID],cen),"pi0 decay photon with non-photon partner", nPt,0.,ptMax )) ;
+
+ fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnCutEmin_%s_cent%d",cPID[iPID],cen),"Decay photons with rec. partner but failed Emin cut", nPt,0.,ptMax )) ;
+ fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnCutNcell_%s_cent%d",cPID[iPID],cen),"Decay photons with rec. partner but failed Ncell cut", nPt,0.,ptMax )) ;
+ fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnCutEcross_%s_cent%d",cPID[iPID],cen),"Decay photons with rec. partner but failed Ecross cut", nPt,0.,ptMax )) ;
+ fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnCutM02_%s_cent%d",cPID[iPID],cen),"Decay photons with rec. partner but failed M02 cut", nPt,0.,ptMax )) ;
+ fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnDefCuts_%s_cent%d",cPID[iPID],cen),"Decay photons with rec. partner but failed default cuts", nPt,0.,ptMax )) ;
+ fOutputContainer->Add(new TH1F(Form("hMCDecWRecPartn_%s_cent%d",cPID[iPID],cen),"Decay photons with rec partner", nPt,0.,ptMax )) ;
+ fOutputContainer->Add(new TH1F(Form("hMCDecWRecUniqPartn_%s_cent%d",cPID[iPID],cen),"Decay photons with rec partner", nPt,0.,ptMax )) ;
+
+ fOutputContainer->Add(new TH1F(Form("hMCDecMerged_%s_cent%d",cPID[iPID],cen),"Decay photons with rec partner", nPt,0.,ptMax )) ;
+ fOutputContainer->Add(new TH1F(Form("hMCDecUnfolded_%s_cent%d",cPID[iPID],cen),"Decay photons with rec partner", nPt,0.,ptMax )) ;
+
+
+ fOutputContainer->Add(new TH2F(Form("hMC_InvM_Re_%s_cent%d",cPID[iPID],cen),"Two-photon inv. mass vs first photon pt",nM,0.,mMax,nPt,0.,ptMax)) ;
+ fOutputContainer->Add(new TH2F(Form("hMC_InvM_Re_Strict_%s_cent%d",cPID[iPID],cen),"Two-photon inv. mass vs first photon pt",nM,0.,mMax,nPt,0.,ptMax)) ;
+ }
+ fOutputContainer->Add(new TH2F(Form("hMCmass_cent%d",cen),"Mass with reconstructed decay partner",nM,0.,mMax,nPt,0.,ptMax )) ;
+ }
+
+/*
fOutputContainer->Add(new TH1F(Form("hMCRecNoLabel_cent%d",cen),"Clusters without label",nPt,0.,ptMax)) ;
fOutputContainer->Add(new TH1F(Form("hMCConversionRadius_cent%d",cen),"Clusters without label",600,0.,600.)) ;
fOutputContainer->Add(new TH2F(Form("hMCRecPi0Vtx_cent%d",cen),"Secondary pi0s",100,0.,10.,600,0.,600.)) ;
//Pi0 decay photons
//MC tagging: reasons of partner loss etc.
- fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnStack_cen%d",cen),"Decay photons with partner not in Stack", nPt,0.,ptMax)) ;
- fOutputContainer->Add(new TH1F(Form("hMCDecWithFoundPartn_cen%d",cen),"Decay photon with found partner", nPt,0.,ptMax)) ;
- fOutputContainer->Add(new TH1F(Form("hMCDecWithWrongMass_cen%d",cen),"Decay photon with wrong mass", nPt,0.,ptMax)) ;
- fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnAccept_cen%d",cen),"Decay photon with parttner not in PHOS acc", nPt,0.,ptMax)) ;
- fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnAcceptFA1_cen%d",cen),"Decay photons with partner missed due geometry Fid. area. 1", nPt,0.,ptMax)) ;
- fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnAcceptFA2_cen%d",cen),"Decay photons with partner missed due geometry Fid. area. 2", nPt,0.,ptMax)) ;
- fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnAcceptFA3_cen%d",cen),"Decay photons with partner missed due geometry Fid. area. 3", nPt,0.,ptMax)) ;
- fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnConv_cen%d",cen),"Decay photons with partner missed due to conversion", nPt,0.,ptMax)) ;
- fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnEmin_cen%d",cen),"Decay photons with partner missed due to low energy", nPt,0.,ptMax)) ;
- fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnOther_cen%d",cen),"Decay photons with partner missed due unknown reason", nPt,0.,ptMax)) ;
- fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnAll_cen%d",cen),"Decay photons with partner missed due to any reason", nPt,0.,ptMax)) ;
- fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnNPhot_cen%d",cen),"pi0 decay photon with non-photon partner", nPt,0.,ptMax)) ;
+ fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnStack_cent%d",cen),"Decay photons with partner not in Stack", nPt,0.,ptMax)) ;
+ fOutputContainer->Add(new TH1F(Form("hMCDecWithFoundPartn_cent%d",cen),"Decay photon with found partner", nPt,0.,ptMax)) ;
+ fOutputContainer->Add(new TH1F(Form("hMCDecWithWrongMass_cent%d",cen),"Decay photon with wrong mass", nPt,0.,ptMax)) ;
+ fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnAccept_cent%d",cen),"Decay photon with parttner not in PHOS acc", nPt,0.,ptMax)) ;
+ fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnAcceptFA1_cent%d",cen),"Decay photons with partner missed due geometry Fid. area. 1", nPt,0.,ptMax)) ;
+ fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnAcceptFA2_cent%d",cen),"Decay photons with partner missed due geometry Fid. area. 2", nPt,0.,ptMax)) ;
+ fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnAcceptFA3_cent%d",cen),"Decay photons with partner missed due geometry Fid. area. 3", nPt,0.,ptMax)) ;
+ fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnConv_cent%d",cen),"Decay photons with partner missed due to conversion", nPt,0.,ptMax)) ;
+ fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnEmin_cent%d",cen),"Decay photons with partner missed due to low energy", nPt,0.,ptMax)) ;
+ fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnOther_cent%d",cen),"Decay photons with partner missed due unknown reason", nPt,0.,ptMax)) ;
+ fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnAll_cent%d",cen),"Decay photons with partner missed due to any reason", nPt,0.,ptMax)) ;
+ fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnNPhot_cent%d",cen),"pi0 decay photon with non-photon partner", nPt,0.,ptMax)) ;
- fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnCutEmin_cen%d",cen),"Decay photons with rec. partner but failed Emin cut", nPt,0.,ptMax)) ;
- fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnCutNcell_cen%d",cen),"Decay photons with rec. partner but failed Ncell cut", nPt,0.,ptMax)) ;
- fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnCutEcross_cen%d",cen),"Decay photons with rec. partner but failed Ecross cut", nPt,0.,ptMax)) ;
- fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnCutM02_cen%d",cen),"Decay photons with rec. partner but failed M02 cut", nPt,0.,ptMax)) ;
- fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnDefCuts_cen%d",cen),"Decay photons with rec. partner but failed default cuts", nPt,0.,ptMax)) ;
- fOutputContainer->Add(new TH1F(Form("hMCDecWRecPartn_cen%d",cen),"Decay photons with rec partner", nPt,0.,ptMax)) ;
+ fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnCutEmin_cent%d",cen),"Decay photons with rec. partner but failed Emin cut", nPt,0.,ptMax)) ;
+ fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnCutNcell_cent%d",cen),"Decay photons with rec. partner but failed Ncell cut", nPt,0.,ptMax)) ;
+ fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnCutEcross_cent%d",cen),"Decay photons with rec. partner but failed Ecross cut", nPt,0.,ptMax)) ;
+ fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnCutM02_cent%d",cen),"Decay photons with rec. partner but failed M02 cut", nPt,0.,ptMax)) ;
+ fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnDefCuts_cent%d",cen),"Decay photons with rec. partner but failed default cuts", nPt,0.,ptMax)) ;
+ fOutputContainer->Add(new TH1F(Form("hMCDecWRecPartn_cent%d",cen),"Decay photons with rec partner", nPt,0.,ptMax)) ;
- fOutputContainer->Add(new TH2F(Form("hMCmass_cen%d",cen),"Mass with reconstructed decay partner",nM,0.,mMax,nPt,0.,ptMax)) ;
+ fOutputContainer->Add(new TH2F(Form("hMCmass_cent%d",cen),"Mass with reconstructed decay partner",nM,0.,mMax,nPt,0.,ptMax)) ;
fOutputContainer->Add(new TH1F(Form("hMCRecGammaPi0Dalitz_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, from pi0",nPt,0.,ptMax)) ;
//all clusters fake tagged
fOutputContainer->Add(new TH1F(Form("hMCAllFakeTagged_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, from pi0",nPt,0.,ptMax)) ;
- fOutputContainer->Add(new TH2F(Form("hMC_InvM_Re_%s_cent%d",cPID[iPID],cen),"Two-photon inv. mass vs first photon pt",nM,0.,mMax,nPt,0.,ptMax)) ;
- fOutputContainer->Add(new TH2F(Form("hMC_InvM_Re_Strict_%s_cent%d",cPID[iPID],cen),"Two-photon inv. mass vs first photon pt",nM,0.,mMax,nPt,0.,ptMax)) ;
}
fOutputContainer->Add(new TH1F(Form("hMCGammaPi0MisPartner_cent%d",cen),"Spectrum of missed partners",nPt,0.,ptMax)) ;
fOutputContainer->Add(new TH2F(Form("hMCGammaPi0MisPartnerEtaPhi_cent%d",cen),"Spectrum of missed partners",100,-0.2,0.2,100,4.5,5.6)) ;
}
- }
- }//If MC handler exists...
+*/
+ }
+
+ //If we work with MC, need to set Sumw2 - we will use weights
+ if(fIsMC){
+ for(Int_t i=0; i<fOutputContainer->GetSize();i++){
+ ((TH1*)fOutputContainer->At(i))->Sumw2() ;
+ }
+ }
+
- for(Int_t i=0;i<1;i++)
+
+ for(Int_t i=0;i<10;i++)
for(Int_t j=0;j<5;j++)
fPHOSEvents[i][j]=0x0 ; //Container for PHOS photons
//Select photons
//Fill QA histograms
//Fill Tagging histogsms
-
const Double_t kEcrossCut=0.98 ;
const Double_t kTOFMaxCut= 100.e-9 ;
const Double_t kTOFMinCut=-100.e-9 ;
// Event selection flags
-
-// FillHistogram("hSelEvents",0) ;
+ // FillHistogram("hSelEvents",0) ;
AliVEvent* event = (AliVEvent*)InputEvent();
if(!event){
return;
}
FillHistogram("hSelEvents",1) ;
+
+ //MC stack init
+ fStack = (TClonesArray*)event->FindListObject(AliAODMCParticle::StdBranchName());
//read geometry if not read yet
- if(fPHOSgeom==0)
+ if(fPHOSgeom==0){
InitGeometry() ;
+ }
- //MC stack init
- fStack = (TClonesArray*)event->FindListObject(AliAODMCParticle::StdBranchName());
-
- if(!fUtils)
- fUtils = new AliAnalysisUtils();
+ if(!fIsFastMC){
+ if(!fUtils)
+ fUtils = new AliAnalysisUtils();
- Bool_t isMB = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kINT7) ;
- Bool_t isPHI7 = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kPHI7);
-
- if((fIsMB && !isMB) || (fIsMB && isPHI7) || (!fIsMB && !isPHI7)){
- PostData(1, fOutputContainer);
- return;
+ Bool_t isMB = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kINT7) ;
+ Bool_t isPHI7 = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kPHI7);
+
+ if((fIsMB && !isMB) || (!fIsMB && !isPHI7)){
+ PostData(1, fOutputContainer);
+ return;
+ }
}
FillHistogram("hSelEvents",2) ;
FillHistogram("hNvertexTracks",event->GetPrimaryVertex()->GetNContributors());
FillHistogram("hZvertex" ,vtx5[2]);
+ if(fIsFastMC){ //vertex from header
+ AliAODMCHeader *cHeaderAOD = dynamic_cast<AliAODMCHeader*>(event->FindListObject(AliAODMCHeader::StdBranchName()));
+ if(!cHeaderAOD){
+ PostData(1, fOutputContainer);
+ return ;
+ }
+ cHeaderAOD->GetVertex(vtx5);
+ }
if (TMath::Abs(vtx5[2]) > 10. ){
PostData(1, fOutputContainer);
return ;
}
-
FillHistogram("hSelEvents",3) ;
+ //Vtx class z-bin
+ Int_t zvtx = TMath::Min(9,Int_t((vtx5[2]+10.)/2.)) ;
+
+ if(!fIsFastMC){
+// if (event->IsPileupFromSPD()){
+// PostData(1, fOutputContainer);
+// return ;
+// }
-// if (event->IsPileupFromSPD()){
-// PostData(1, fOutputContainer);
-// return ;
-// }
+ if(!fUtils->IsVertexSelected2013pA(event)){
+ PostData(1, fOutputContainer);
+ return ;
+ }
- if(!fUtils->IsVertexSelected2013pA(event)){
- PostData(1, fOutputContainer);
- return ;
- }
- FillHistogram("hSelEvents",4) ;
+ FillHistogram("hSelEvents",4) ;
- if(fUtils->IsPileUpEvent(event)){
- PostData(1, fOutputContainer);
- return ;
+ if(fUtils->IsPileUpEvent(event)){
+ PostData(1, fOutputContainer);
+ return ;
+ }
+ FillHistogram("hSelEvents",5) ;
}
- FillHistogram("hSelEvents",5) ;
//centrality
- AliCentrality *centrality = event->GetCentrality();
- if( centrality )
- fCentrality=centrality->GetCentralityPercentile("V0M");
- else {
- AliError("Event has 0x0 centrality");
- fCentrality = -1.;
- }
- FillHistogram("hCentrality",fCentrality) ;
+ if(!fIsFastMC){
+ AliCentrality *centrality = event->GetCentrality();
+ if( centrality )
+ fCentrality=centrality->GetCentralityPercentile("V0M");
+ else {
+ AliError("Event has 0x0 centrality");
+ fCentrality = -1.;
+ }
+ FillHistogram("hCentrality",fCentrality) ;
- if(fCentrality<0. || fCentrality>=100.){
- PostData(1, fOutputContainer);
- return ;
+ if(fCentrality<0. || fCentrality>=100.){
+ PostData(1, fOutputContainer);
+ return ;
+ }
+ }
+ else{
+ fCentrality=1.;
}
fCentBin = (Int_t)(fCentrality/20.) ;
FillHistogram("hSelEvents",6) ;
- //Vtx class z-bin
- Int_t zvtx = 0 ;
-
//Calculate charged multiplicity
Int_t trackMult = 0;
if(fTrackEvent)
fTrackEvent->Clear() ;
else
- fTrackEvent = new TClonesArray("AliAODPWG4Particle",event->GetNumberOfTracks()) ;
+ fTrackEvent = new TClonesArray("AliCaloPhoton",event->GetNumberOfTracks()) ;
for (Int_t i=0;i<event->GetNumberOfTracks();++i) {
AliAODTrack *track = (AliAODTrack*)event->GetTrack(i) ;
if(TMath::Abs(track->Eta())< 0.9){
if(trackMult>=fTrackEvent->GetSize())
fTrackEvent->Expand(2*trackMult) ;
- new ((*fTrackEvent)[trackMult]) AliAODPWG4Particle(track->Px(),track->Py(),track->Pz(),track->P());
+ new ((*fTrackEvent)[trackMult]) AliCaloPhoton(track->Px(),track->Py(),track->Pz(),track->P());
trackMult++;
FillHistogram("hTrackEtaPhi",track->Eta(),track->Phi()) ;
FillHistogram("hTrackEtaPt",track->Eta(),track->Pt()) ;
//---------Select photons-------------------
Int_t multClust = event->GetNumberOfCaloClusters();
if(!fPHOSEvent)
- fPHOSEvent = new TClonesArray("AliAODPWG4Particle",multClust);
+ fPHOSEvent = new TClonesArray("AliCaloPhoton",multClust);
else
fPHOSEvent->Clear() ;
Int_t inList = 0; //counter of caloClusters
if(clu->GetMCEnergyFraction()>kEcrossCut) //Ecross cut, should be filled with Tender
continue ;
+
+ if(clu->GetDistanceToBadChannel()<fMinBCDistance)
+ continue ;
Float_t pos[3] ;
clu->GetPosition(pos) ;
Int_t mod = relId[0] ;
Int_t cellX = relId[2];
Int_t cellZ = relId[3] ;
-
+
FillHistogram("hTOF",clu->E(),clu->GetTOF()) ;
FillHistogram(Form("hTofM%d",mod),clu->GetTOF(),clu->E()) ;
- if(clu->GetTOF() < kTOFMinCut || clu->GetTOF() > kTOFMaxCut)
+ if((!fIsMC) && (clu->GetTOF() < kTOFMinCut || clu->GetTOF() > kTOFMaxCut))
continue ;
// if(clu->GetDistanceToBadChannel()<2.5)
TLorentzVector momentum ;
clu->GetMomentum(momentum, vtx5);
- AliAODPWG4Particle *p = new ((*fPHOSEvent)[inList]) AliAODPWG4Particle(momentum.Px(),momentum.Py(),momentum.Pz(),clu->E() );
+
+ AliCaloPhoton *p = new ((*fPHOSEvent)[inList]) AliCaloPhoton(momentum.Px(),momentum.Py(),momentum.Pz(),clu->E() );
+// momentum*= clu->GetCoreEnergy()/momentum.E() ;
+// AliCaloPhoton *p = new ((*fPHOSEvent)[inList]) AliCaloPhoton(momentum.Px(),momentum.Py(),momentum.Pz(),clu->GetCoreEnergy() );
inList++;
-
- Int_t isolation = EvalIsolation(&momentum) ;
- p->SetBtag(isolation) ;
- p->SetCaloLabel(i,-1); //This and partner cluster
- p->SetDistToBad((Int_t)(1.+clu->GetDistanceToBadChannel()/2.2));
+ Int_t isolation = EvalIsolation(&momentum,kTRUE) ;
+ p->SetIsolationTag(isolation) ;
- p->SetTag(0); //Strict PID pi0 partner not found
+ p->SetDistToBad((Int_t)(1.+clu->GetDistanceToBadChannel()/2.2));
+ p->SetBC(i) ; //reference to CaloCluster
+ p->SetTagInfo(0); //No pi0 partners found so far
p->SetTagged(kFALSE); //Reconstructed pairs found
- Bool_t sure = 0;
p->SetFiducialArea(fidArea) ;
- if(fStack){
+ if(fIsMC){
//This is primary entered PHOS
- Int_t primLabel=FindPrimary(clu,sure) ;
+ FillHistogram(Form("LabelsNPrim_cent%d",fCentBin),clu->E(),float(clu->GetNLabels())) ;
+ Int_t primLabel=clu->GetLabelAt(0) ; //FindPrimary(clu,sure) ;
//Look what particle left vertex
if(primLabel>-1){
AliAODMCParticle * prim = (AliAODMCParticle*)fStack->At(primLabel) ;
while((r2 > rcut*rcut) && (iparent>-1)){
iparent=parent->GetMother();
parent=(AliAODMCParticle*)fStack->At(iparent);
+ r2=parent->Xv()*parent->Xv()+parent->Yv()*parent->Yv() ;
}
- p->SetCaloLabel(primLabel,iparent); //This and partner cluster
+ p->SetPrimary(primLabel) ;
+ p->SetPrimaryAtVertex(iparent) ;
+ p->SetWeight(PrimaryParticleWeight(parent)) ;
}
else{
- p->SetCaloLabel(-1,-1); //This and partner cluster
+ p->SetPrimary(-1); //Primary index
+ p->SetPrimaryAtVertex(-1) ;
+ p->SetWeight(1.) ;
}
}
else{
- //This is primary at vertex(R<1cm)
- p->SetCaloLabel(-1,-1); //This and partner cluster
+ p->SetPrimary(-1); //Primary index
+ p->SetPrimaryAtVertex(-1) ;
+ p->SetWeight(1.) ;
}
- p->SetLabel(i); //Cluster index
//PID criteria
- p->SetDispBit(clu->Chi2()<2.5) ;
+// Cut on full Lambdas
+// p->SetDispBit(clu->Chi2()<2.5*2.5) ;
+// Cut on CoreLamdas
+ p->SetDispBit(clu->GetDispersion()<2.5*2.5) ;
p->SetTOFBit(TestTOF(clu->GetTOF(),clu->E())) ;
- p->SetChargedBit(clu->GetEmcCpvDistance()>2.5) ;
+ p->SetCPVBit(clu->GetEmcCpvDistance()>2.5) ;
+
+ FillHistogram(Form("hPhotM%d_All",mod),p->Pt()) ;
+ if(fidArea>1){
+ FillHistogram(Form("hPhotM%dA2_All",mod),p->Pt()) ;
+ if(fidArea>2){
+ FillHistogram(Form("hPhotM%dA3_All",mod),p->Pt()) ;
+ }
+ }
+ if(p->IsDispOK()){
+ FillHistogram(Form("hPhotM%d_Disp",mod),p->Pt()) ;
+ if(fidArea>1){
+ FillHistogram(Form("hPhotM%dA2_Disp",mod),p->Pt()) ;
+ if(fidArea>2){
+ FillHistogram(Form("hPhotM%dA3_Disp",mod),p->Pt()) ;
+ }
+ }
+ if(p->IsCPVOK()){
+ FillHistogram(Form("hPhotM%d_Both",mod),p->Pt()) ;
+ if(fidArea>1){
+ FillHistogram(Form("hPhotM%dA2_Both",mod),p->Pt()) ;
+ if(fidArea>2){
+ FillHistogram(Form("hPhotM%dA3_Both",mod),p->Pt()) ;
+ }
+ }
+ }
+ }
+ if(p->IsCPVOK()){
+ FillHistogram(Form("hPhotM%d_CPV",mod),p->Pt()) ;
+ if(fidArea>1){
+ FillHistogram(Form("hPhotM%dA2_CPV",mod),p->Pt()) ;
+ if(fidArea>2){
+ FillHistogram(Form("hPhotM%dA3_CPV",mod),p->Pt()) ;
+ }
+ }
+ }
}
FillHistogram("hPHOSCentrality",fCentrality,inList+0.5) ;
+
FillMCHistos() ;
FillTaggingHistos() ;
void AliAnalysisTaskTaggedPhotons::FillMCHistos(){
//MC info about this particle
- if(!fStack)
+ if(!fIsMC)
return ;
const Double_t rcut=1. ; //cut on vertex to consider particle as "primary"
const Double_t phiMin=260.*TMath::Pi()/180. ;
const Double_t phiMax=320.*TMath::Pi()/180. ;
+
+ AliVEvent* event = (AliVEvent*)InputEvent();
Int_t nPrim = fStack->GetEntriesFast() ;
//Fill Primary particl yields
if(abs(pdg) == 211)
snprintf(partName,30,"pipm") ;
else
- continue ;
-
+ if(abs(pdg) == 2212)
+ snprintf(partName,30,"p") ;
+ else
+ if(abs(pdg) ==-2212)
+ snprintf(partName,30,"pbar") ;
+ else
+ if(abs(pdg) == 2112)
+ snprintf(partName,30,"n") ;
+ else
+ if(abs(pdg) ==-2112)
+ snprintf(partName,30,"nbar") ;
+ else
+ continue ;
//Primary particle
Double_t phi=prim->Phi() ;
//Total number of pi0 with creation radius <1 cm
Double_t w = PrimaryParticleWeight(prim) ;
- FillHistogram(Form("hMC_all_%s",partName),pt,w) ;
+ FillHistogram(Form("hMC_all_%s_cent%d",partName,fCentBin),pt,w) ;
if(TMath::Abs(prim->Y())<0.13){
- FillHistogram(Form("hMC_phi_%s",partName),phi,w) ;
- if(phi>phiMin && phi<phiMax)
- FillHistogram(Form("hMC_%s_pt_cen%d",partName,fCentBin),pt,w) ;
+ FillHistogram(Form("hMC_phi_%s_cent%d",partName,fCentBin),phi,w) ;
+ if(phi>phiMin && phi<phiMax){
+ FillHistogram(Form("hMC_unitEta_%s_cent%d",partName,fCentBin),pt,w) ;
+ if(prim->GetMother()>-1){
+ AliAODMCParticle * primGM = (AliAODMCParticle*)fStack->At(prim->GetMother()) ;
+ if(primGM->GetPdgCode()==111)
+ FillHistogram(Form("hMC_unitEta_GammaPi0_cent%d",fCentBin),pt,w) ;
+ }
+ }
}
- FillHistogram(Form("hMC_rap_%s",partName),prim->Y(),w) ;
+ FillHistogram(Form("hMC_rap_%s_cent%d",partName,fCentBin),prim->Y(),w) ;
//Some additional QA
if(pdg == 111){
- FillHistogram(Form("hMCpi0_ptrap_cen%d",fCentBin),pt,prim->Y(),w) ;
- FillHistogram(Form("hMCpi0_ptphi_cen%d",fCentBin),pt,phi,w) ;
+ FillHistogram("hMCpi0_ptrap",pt,prim->Y(),w) ;
+ FillHistogram("hMCpi0_ptphi",pt,phi,w) ;
}
if(pdg == 22){
- FillHistogram(Form("hMCgamma_ptrap_cen%d",fCentBin),pt,prim->Y(),w) ;
- FillHistogram(Form("hMCgamma_ptphi_cen%d",fCentBin),pt,phi,w) ;
+ FillHistogram("hMCgamma_ptrap",pt,prim->Y(),w) ;
+ FillHistogram("hMCgamma_ptphi",pt,phi,w) ;
}
}
const Int_t n=fPHOSEvent->GetEntriesFast() ;
for(Int_t i=0;i<n;i++){
- AliAODPWG4Particle *p = static_cast<AliAODPWG4Particle*>(fPHOSEvent->At(i));
- Int_t label=p->GetLabel() ;
+ AliCaloPhoton *p = static_cast<AliCaloPhoton*>(fPHOSEvent->At(i));
+ Int_t label=p->GetPrimary() ;
if(label<0){ //No label!
- FillHistogram("hMCRecNoLabel",p->Pt());
+ FillHistogram("hMCRecNoLabel",p->Pt(),p->GetWeight());
continue ;
}
- AliAODMCParticle * prim = (AliAODMCParticle*)fStack->At(p->GetLabel()) ;
+ AliAODMCParticle * prim = (AliAODMCParticle*)fStack->At(p->GetPrimary()) ;
//Look what particle left virtex
- Int_t iparent=p->GetLabel();
+ Int_t iparent=p->GetPrimary();
AliAODMCParticle * parent = prim;
while(parent->Xv()*parent->Xv()+parent->Yv()*parent->Yv() > rcut*rcut){
iparent=parent->GetMother();
break ;
parent = (AliAODMCParticle*)fStack->At(iparent) ;
}
- Int_t parentPDG=parent->GetPdgCode() ;
+ Int_t parentPDG=parent->GetPdgCode() ;
+ Int_t iFidArea=p->GetFiducialArea();
switch(parentPDG){
case 22: //electron/positron conversion
FillPIDHistograms("hMCRecPhoton",p); //Reconstructed with photon from conversion primary
+ if(iFidArea>0){
+ FillPIDHistograms("hMCRecPhoton_Area1",p);
+ if(iFidArea>1){
+ FillPIDHistograms("hMCRecPhoton_Area2",p);
+ if(iFidArea>1){
+ FillPIDHistograms("hMCRecPhoton_Area3",p);
+ }
+ }
+ }
break ;
case 11:
case -11: //electron/positron conversion
break ;
case 211:
case -211:
+ FillPIDHistograms("hMCRecPipm",p); //Reconstructed with photon from antibaryon annihilation
+ break ;
case 2212:
+ FillPIDHistograms("hMCRecP",p); //Reconstructed with photon from antibaryon annihilation
+ break ;
case 321:
case -321:
- FillPIDHistograms("hMCRecCharg",p); //Reconstructed with photon from conversion primary
+ FillPIDHistograms("hMCRecKpm",p); //Reconstructed with photon from conversion primary
break ;
case 310:
FillPIDHistograms("hMCRecK0s",p); //Reconstructed with photon from conversion primary
break ;
case 2112: //antineutron & antiproton conversion
- case 130:
- FillPIDHistograms("hMCRecNeutral",p); //Reconstructed with photon from antibaryon annihilation
- break ;
+ FillPIDHistograms("hMCRecN",p); //Reconstructed with photon from antibaryon annihilation
+ break ;
case -1: //direct photon or no primary
FillPIDHistograms("hMCRecNoPRim",p);
break ;
default:
- printf("Unknown PDG: %d \n",parentPDG) ;
- FillPIDHistograms("hMCRecUnknown",p);
- break ;
+ if(parent->Charge()!=0)
+ FillPIDHistograms("hMCRecCharg",p); //Reconstructed with photon from antibaryon annihilation
+ else
+ FillPIDHistograms("hMCRecNeutral",p); //Reconstructed with photon from antibaryon annihilation
}
ipartner=-1 ;
}
}
-
//There is no partner in stack
if(ipartner==-1){
FillPIDHistograms("hMCDecWMisPartnStack",p) ;
AliAODMCParticle * partner = (AliAODMCParticle *)fStack->At(ipartner);
//Check if partner is registered and made correct mass
//If not - trace the reason
- AliAODPWG4Particle *pp = 0x0 ;
+ AliCaloPhoton *pp = 0x0 ;
for(Int_t ii=0;ii<n;ii++){
if(i==ii) continue;
- AliAODPWG4Particle * tmp = static_cast<AliAODPWG4Particle*>(fPHOSEvent->At(ii));
- Int_t ipartnPrim = tmp->GetLabel() ;
+ AliCaloPhoton * tmp = static_cast<AliCaloPhoton*>(fPHOSEvent->At(ii));
+ Int_t ipartnPrim = tmp->GetPrimary() ;
while(ipartnPrim>-1){
- if(ipartnPrim==ipartner)
+ if(ipartnPrim==ipartner){
break ;
+ }
ipartnPrim = ((AliAODMCParticle *)fStack->At(ipartnPrim))->GetMother();
}
if(ipartnPrim==ipartner){
break ;
}
}
+ if(pp){
+ //Partner reconstructed, but did not pass cuts
+ FillPIDHistograms("hMCDecWRecUniqPartn",p) ;
+ }
+ //Partner not found. Check if it is not dominant contributor?
+ if(!pp){
+ for(Int_t ii=0;(ii<n) && (!pp);ii++){
+ if(i==ii) continue;
+ AliCaloPhoton * tmp = static_cast<AliCaloPhoton*>(fPHOSEvent->At(ii));
+ Int_t iCaloCluster=tmp->GetBC();//index of AODCaloCluster
+ AliVCluster * clu = event->GetCaloCluster(iCaloCluster);
+ Int_t nCluPrimaries = clu->GetNLabels() ;
+ for(Int_t iAODLabel=0; (iAODLabel<nCluPrimaries) && (!pp); iAODLabel++){
+ Int_t ipartnPrim = clu->GetLabelAt(iAODLabel) ;
+ while(ipartnPrim>-1){
+ if(ipartnPrim==ipartner){
+ break ;
+ }
+ ipartnPrim = ((AliAODMCParticle *)fStack->At(ipartnPrim))->GetMother();
+ }
+ if(ipartnPrim==ipartner){
+ pp=tmp ;
+ break ;
+ }
+ }
+ }
+ }
+ //If partner still not found, check if it is in same cluster
+ if(!pp){
+ Int_t iCaloCluster=p->GetBC();//index of AODCaloCluster
+ AliVCluster * clu = event->GetCaloCluster(iCaloCluster);
+ Int_t nCluPrimaries = clu->GetNLabels() ;
+ for(Int_t iAODLabel=0; iAODLabel<nCluPrimaries ; iAODLabel++){
+ Int_t ipartnPrim = clu->GetLabelAt(iAODLabel) ;
+ while(ipartnPrim>-1){
+ if(ipartnPrim==ipartner){
+ break ;
+ }
+ ipartnPrim = ((AliAODMCParticle *)fStack->At(ipartnPrim))->GetMother();
+ }
+ if(ipartnPrim==ipartner){ //yes, this cluster contains both primary
+ if(clu->GetNExMax()<2){ //was not unfolded
+ FillPIDHistograms("hMCDecMerged",p) ;
+ }
+ else{
+ FillPIDHistograms("hMCDecUnfolded",p) ;
+ }
+ }
+ }
+ }
if(pp){
//Partner reconstructed, but did not pass cuts
- Bool_t cutEmin = (pp->GetInputFileIndex() &1) ;
- Bool_t cutNcell = (pp->GetInputFileIndex() &2) ;
- Bool_t cutEcross = (pp->GetInputFileIndex() &4) ;
- Bool_t cutM02 = (pp->GetInputFileIndex() &8) ;
-
- if(!cutEmin)FillPIDHistograms("hMCDecWMisPartnCutEmin",p) ;
- if(!cutNcell)FillPIDHistograms("hMCDecWMisPartnCutNcell",p) ;
- if(!cutEcross)FillPIDHistograms("hMCDecWMisPartnCutEcross",p) ;
- if(!cutM02)FillPIDHistograms("hMCDecWMisPartnCutM02",p) ;
-
- if(cutEmin && cutNcell && cutEcross && cutM02){
FillPIDHistograms("hMCDecWRecPartn",p) ;
- Double_t invMass=p->GetPairMass(pp) ;
- FillHistogram("hMCmass",invMass,p->Pt()) ;
- if(IsInPi0Band(invMass,p->Pt())){
- FillPIDHistograms("hMCDecWithFoundPartn",p) ;
+ Double_t invMass=(*p+ *pp).M() ;
+ FillHistogram(Form("hMCmass_cent%d",fCentBin),invMass,p->Pt(),p->GetWeight()) ;
+ Double_t nSigma=InPi0Band(invMass,p->Pt()) ;
+ // analog to Tag
+ for(Int_t eminType=0; eminType<3; eminType++){
+ if(pp->E()>0.1*(eminType+1)){
+ for(Int_t isigma=0; isigma<3; isigma++){
+ if(nSigma<1.+isigma){
+ Int_t iType=3*eminType+isigma ;
+ FillPIDHistograms(Form("hMCDecWithFoundPartnType%d",iType),p) ;
+ }
+ }
+ }
}
- else{
+ if(nSigma>3.){
FillPIDHistograms("hMCDecWithWrongMass",p) ;
}
- }
- else{
- FillPIDHistograms("hMCDecWMisPartnDefCuts",p) ;
- }
}
else{//Partner not reconstructed
if(partner->GetPdgCode()==22){
if(!impact){ //this photon cannot hit PHOS
FillPIDHistograms("hMCDecWMisPartnAccept",p) ; //Spectrum of tagged with missed partner
- Int_t iFidArea = p->GetFiducialArea();
if(iFidArea>0){
FillPIDHistograms("hMCDecWMisPartnAcceptFA1",p) ; //Spectrum of tagged with missed partner
if(iFidArea>1){
}
if(!isPartnerLost){ //Reason not found!!!!!
FillPIDHistograms("hMCDecWMisPartnOther",p);
+ Int_t multClust = event->GetNumberOfCaloClusters();
+ for (Int_t iclu=0; (iclu<multClust) && (!isPartnerLost); iclu++) {
+ AliVCluster * clu = event->GetCaloCluster(iclu);
+ if(!clu->IsPHOS())
+ continue ;
+ if(clu->E()==p->E()) //same cluster as current
+ continue ;
+ Int_t nCluPrimaries = clu->GetNLabels() ;
+ for(Int_t iAODLabel=0; (iAODLabel<nCluPrimaries) && (!isPartnerLost); iAODLabel++){
+ Int_t ipartnPrim = clu->GetLabelAt(iAODLabel) ;
+ while(ipartnPrim>-1){
+ if(ipartnPrim==ipartner)
+ break ;
+ ipartnPrim = ((AliAODMCParticle *)fStack->At(ipartnPrim))->GetMother();
+ }
+ if(ipartnPrim==ipartner){
+ isPartnerLost=kTRUE;
+ break ;
+ }
+ }
+ }
+ if(isPartnerLost){//Did not pass default cuts
+ FillPIDHistograms("hMCDecWMisPartnDefCuts",p);
+ }
}
else{//Sum of all missed partners
FillPIDHistograms("hMCDecWMisPartnAll",p);
//Invariant Mass analysis
const Int_t n=fPHOSEvent->GetEntriesFast() ;
for(Int_t i=0;i<n-1;i++){
- AliAODPWG4Particle *p1 = static_cast<AliAODPWG4Particle*>(fPHOSEvent->At(i));
+ AliCaloPhoton *p1 = static_cast<AliCaloPhoton*>(fPHOSEvent->At(i));
for(Int_t j = i+1 ; j < n ; j++) {
- AliAODPWG4Particle * p2 = static_cast<AliAODPWG4Particle*>(fPHOSEvent->At(j));
+ AliCaloPhoton * p2 = static_cast<AliCaloPhoton*>(fPHOSEvent->At(j));
- Double_t invMass = p1->GetPairMass(p2);
+ Double_t invMass = (*p1 + *p2).M();
+
+ Double_t nsigma1 = InPi0Band(invMass,p1->Pt()) ; //in band with n sigmas
+ Double_t nsigma2 = InPi0Band(invMass,p2->Pt()) ; //in band with n sigmas
if((p1->E()>0.1) && (p2->E()>0.1)){
- FillPIDHistograms("hInvM_Re_Emin1",p1,p2,invMass) ;
+ FillPIDHistograms("hInvM_Re_Emin1",p1,p2,invMass,kTRUE) ;
if((p1->E()>0.2) && (p2->E()>0.2)){
- FillPIDHistograms("hInvM_Re_Emin2",p1,p2,invMass) ;
+ FillPIDHistograms("hInvM_Re_Emin2",p1,p2,invMass,kTRUE) ;
if((p1->E()>0.3) && (p2->E()>0.3)){
- FillPIDHistograms("hInvM_Re_Emin3",p1,p2,invMass) ;
+ FillPIDHistograms("hInvM_Re_Emin3",p1,p2,invMass,kTRUE) ;
+
+ //Fill izolated pi0s
+ if(nsigma1<2 || nsigma2<2){ //2 sigma band
+ TLorentzVector pi0=*p1+*p2 ;
+ Int_t isolation=EvalIsolation(&pi0,0) ;
+ for(Int_t kind=1; kind<33; kind*=2){
+ if((isolation&kind)){
+ FillHistogram(Form("hPi_Isolation%d_cent%d",kind,fCentBin),pi0.Pt()) ;
+ }
+ }
+ }
+
}
}
}
}
}
}
- if(IsSamePi0(p1,p2)){
+ if(IsSameParent(p1,p2)==111){
FillPIDHistograms("hMC_InvM_Re",p1,invMass) ;
FillPIDHistograms("hMC_InvM_Re",p2,invMass) ;
if(TestPID(3, p2)){
//Tagging: 1,2,3 sigma
//Emin=100,200,300 Mev
- //IsInPi0Band(mass, Ptphoton, type Emin cut
+ //InPi0Band(mass, Ptphoton, type Emin cut
+ //Set Bits to 1 if tagged
Int_t tag1=0 ;
for(Int_t eminType=0; eminType<3; eminType++){
if(p2->E()>0.1*(eminType+1)){
//Set corresponding bit
- Double_t nsigma = IsInPi0Band(invMass,p1->Pt()) ; //in band with n sigmas
for(Int_t isigma=0; isigma<3; isigma++){
- if(nsigma<1+isigma){
+ if(nsigma1<1+isigma){
tag1|= (1 << (3*eminType+isigma)) ;
- if(TestPID(3, p2))
- tag1|= (1 << (3*eminType+isigma+9)) ;
+ if(p2->IsPIDOK(3))
+ tag1|= (1 << (3*eminType+isigma+9)) ;
}
}
}
}
- p1->SetTag(tag1) ;
+ Int_t oldTag1=p1->GetTagInfo() ;
+ for(Int_t ibit=0; ibit<18; ibit++){
+ if(((oldTag1 & (1<<ibit))!=0) && //Already tagged
+ ((tag1 & (1<<ibit))!=0)){//Multiple tagging
+ Int_t iFidArea = p1->GetFiducialArea();
+ if(iFidArea>0){
+ FillPIDHistograms(Form("hPhot_TaggedMult%d_Area1",ibit),p1) ;
+ if(iFidArea>1){
+ FillPIDHistograms(Form("hPhot_TaggedMult%d_Area2",ibit),p1) ;
+ if(iFidArea>2){
+ FillPIDHistograms(Form("hPhot_TaggedMult%d_Area3",ibit),p1) ;
+ }
+ }
+ }
+ }
+ }
+ tag1=tag1|oldTag1 ;
+ p1->SetTagInfo(tag1) ;
Int_t tag2=0 ;
for(Int_t eminType=0; eminType<3; eminType++){
if(p1->E()>0.1*(eminType+1)){
//Set corresponding bit
- Double_t nsigma = IsInPi0Band(invMass,p2->Pt()) ; //in band with n sigmas
for(Int_t isigma=0; isigma<3; isigma++){
- if(nsigma<1+isigma){
+ if(nsigma2<1+isigma){
tag2|= (1 << (3*eminType+isigma)) ;
- if(TestPID(3, p2))
- tag2|= (1 << (3*eminType+isigma+9)) ;
+ if(p1->IsPIDOK(3))
+ tag2|= (1 << (3*eminType+isigma+9)) ;
}
}
}
}
- p2->SetTag(tag2) ;
-
- if(tag1 & (1<<7)){ //2 sigma, Emin=0.3: default tagging
- if(p1->IsTagged()){//Multiple tagging
- FillHistogram(Form("hTaggedMult_cent%d",fCentBin),p1->Pt());
- }
- p1->SetTagged(kTRUE) ;
+ Int_t oldTag2=p2->GetTagInfo() ;
+ for(Int_t ibit=0; ibit<18; ibit++){
+ if(((oldTag2 & (1<<ibit))!=0) && //Already tagged
+ ((tag2 & (1<<ibit))!=0)){//Multiple tagging
+ Int_t iFidArea = p2->GetFiducialArea();
+ if(iFidArea>0){
+ FillPIDHistograms(Form("hPhot_TaggedMult%d_Area1",ibit),p2) ;
+ if(iFidArea>1){
+ FillPIDHistograms(Form("hPhot_TaggedMult%d_Area2",ibit),p2) ;
+ if(iFidArea>2){
+ FillPIDHistograms(Form("hPhot_TaggedMult%d_Area3",ibit),p2) ;
+ }
+ }
+ }
+ }
}
- if(tag2 & (1<<7)){ //2 sigma, Emin=0.3: default tagging
- if(p2->IsTagged()){//Multiple tagging
- FillHistogram(Form("hTaggedMult_cent%d",fCentBin),p2->Pt());
- }
- p2->SetTagged(kTRUE) ;
- }
+ tag2=tag2|oldTag2 ;
+ p2->SetTagInfo(tag2) ;
+
+ Bool_t trueTag=(IsSameParent(p1,p2)==111); //true tag
+ p1->SetUnfolded(p1->IsntUnfolded()|trueTag) ;
+ p2->SetUnfolded(p2->IsntUnfolded()|trueTag) ;
+
}
}
//Single particle histgams
for(Int_t i=0;i<n;i++){
- AliAODPWG4Particle *p = static_cast<AliAODPWG4Particle*>(fPHOSEvent->At(i));
+ AliCaloPhoton *p = static_cast<AliCaloPhoton*>(fPHOSEvent->At(i));
- Int_t isolation = p->GetBtag();
+ Int_t isolation = p->GetIsolationTag();
//Inclusive spectra
FillPIDHistograms("hPhot",p) ;
FillPIDHistograms(Form("hPhot_Isolation%d",kind),p) ;
}
}
-
+
Int_t iFidArea = p->GetFiducialArea();
if(iFidArea>0){
FillPIDHistograms("hPhot_Area1",p) ;
//3 Emin cuts
//Default Emin, 1,2,3 sigmas
//strict and loose PID cut on partner
- Int_t tag=p->GetTag() ;
+ Int_t tag=p->GetTagInfo() ;
for(Int_t ibit=0; ibit<18; ibit++){
- if((tag & (1<<ibit))==0){
- FillPIDHistograms(Form("hPhot_nTagged%d_Area1",ibit),p) ;
-// for(Int_t kind=1; kind<33; kind*=2){
-// if((isolation&kind)){
-// FillPIDHistograms(Form("hPhot_nStrictTagged_Isolation%d_Area1",kind),p) ;
-// }
-// }
+ if((tag & (1<<ibit))!=0){
+ FillPIDHistograms(Form("hPhot_Tagged%d_Area1",ibit),p) ;
+ if(p->IsntUnfolded()) //true tag
+ FillPIDHistograms(Form("hPhot_TrueTagged%d",ibit),p) ;
+
}
}
if(iFidArea>1){
FillPIDHistograms("hPhot_Area2",p) ;
for(Int_t ibit=0; ibit<18; ibit++){
- if((tag & (1<<ibit))==0){
- FillPIDHistograms(Form("hPhot_nTagged%d_Area2",ibit),p) ;
+ if((tag & (1<<ibit))!=0){
+ FillPIDHistograms(Form("hPhot_Tagged%d_Area2",ibit),p) ;
}
}
if(iFidArea>2){
FillPIDHistograms("hPhot_Area3",p) ;
for(Int_t ibit=0; ibit<18; ibit++){
- if((tag & (1<<ibit))==0){
- FillPIDHistograms(Form("hPhot_nTagged%d_Area3",ibit),p) ;
+ if((tag & (1<<ibit))!=0){
+ FillPIDHistograms(Form("hPhot_Tagged%d_Area3",ibit),p) ;
}
}
}
//Fill Mixed InvMass distributions:
TIter nextEv(fCurrentMixedList) ;
for(Int_t i=0;i<n;i++){
- AliAODPWG4Particle *p1 = static_cast<AliAODPWG4Particle*>(fPHOSEvent->At(i));
+ AliCaloPhoton *p1 = static_cast<AliCaloPhoton*>(fPHOSEvent->At(i));
while(TClonesArray * event2 = static_cast<TClonesArray*>(nextEv())){
Int_t nPhotons2 = event2->GetEntriesFast() ;
for(Int_t j=0; j < nPhotons2 ; j++){
- AliAODPWG4Particle * p2 = static_cast<AliAODPWG4Particle*>(event2->At(j)) ;
- Double_t invMass = p1->GetPairMass(p2);
+ AliCaloPhoton * p2 = static_cast<AliCaloPhoton*>(event2->At(j)) ;
+ Double_t invMass = (*p1 + *p2).M();
if((p1->E()>0.1) && (p2->E()>0.1)){
- FillPIDHistograms("hInvM_Mi_Emin1",p1,p2,invMass) ;
+ FillPIDHistograms("hInvM_Mi_Emin1",p1,p2,invMass,kFALSE) ;
if((p1->E())>0.2 && (p2->E()>0.2)){
- FillPIDHistograms("hInvM_Mi_Emin2",p1,p2,invMass) ;
+ FillPIDHistograms("hInvM_Mi_Emin2",p1,p2,invMass,kFALSE) ;
if((p1->E())>0.3 && (p2->E()>0.3)){
- FillPIDHistograms("hInvM_Mi_Emin3",p1,p2,invMass) ;
+ FillPIDHistograms("hInvM_Mi_Emin3",p1,p2,invMass,kFALSE) ;
}
}
}
if (fDebug > 1) Printf("Terminate()");
}
//______________________________________________________________________________
-Bool_t AliAnalysisTaskTaggedPhotons::IsInPi0Band(Double_t m, Double_t pt)const
+Double_t AliAnalysisTaskTaggedPhotons::InPi0Band(Double_t m, Double_t pt)const
{
//Parameterization of the pi0 peak region
//for LHC13bcdef
- Double_t mpi0mean = 0.13447 - 1.41259e-03 * TMath::Exp(-pt/1.30044) ;
+// Double_t mpi0mean = 0.13447 - 1.41259e-03 * TMath::Exp(-pt/1.30044) ;
+ //Parameterization of data 30.08.2014
+ Double_t mpi0mean = 0.135 ;
- Double_t mpi0sigma=TMath::Sqrt(5.22245e-03*5.22245e-03 +2.86851e-03*2.86851e-03/pt) + 9.09932e-05*pt ;
-
- return (m>mpi0mean-2*mpi0sigma && m<mpi0mean+2*mpi0sigma) ;
+ //Double_t mpi0sigma=TMath::Sqrt(5.22245e-03*5.22245e-03 +2.86851e-03*2.86851e-03/pt) + 9.09932e-05*pt ;
+ //Parameterization of data 30.08.2014
+ Double_t mpi0sigma=TMath::Sqrt(4.67491e-03*4.67491e-03 +3.42884e-03*3.42884e-03/pt) + 1.24324e-04*pt ;
+
+ return TMath::Abs(m-mpi0mean)/mpi0sigma ;
}
//______________________________________________________________________________
-Bool_t AliAnalysisTaskTaggedPhotons::IsSamePi0(const AliAODPWG4Particle *p1, const AliAODPWG4Particle *p2)const{
+Int_t AliAnalysisTaskTaggedPhotons::IsSameParent(const AliCaloPhoton *p1, const AliCaloPhoton *p2)const{
//Looks through parents and finds if there was commont pi0 among ancestors
- if(!fStack)
- return kFALSE ; //can not say anything
+ if(!fIsMC)
+ return 0 ; //can not say anything
- Int_t prim1 = p1->GetCaloLabel(0);
+ Int_t prim1 = p1->GetPrimary();
while(prim1!=-1){
- Int_t prim2 = p2->GetCaloLabel(0);
- while(prim2!=-1){
+ Int_t prim2 = p2->GetPrimary();
+
+ while(prim2!=-1){
if(prim1==prim2){
- if(((AliAODMCParticle*)fStack->At(prim1))->GetPdgCode()==111)
- return kTRUE ;
- else
- return kFALSE ;
+ return ((AliAODMCParticle*)fStack->At(prim1))->GetPdgCode() ;
}
prim2=((AliAODMCParticle*)fStack->At(prim2))->GetMother() ;
}
prim1=((AliAODMCParticle*)fStack->At(prim1))->GetMother() ;
}
- return kFALSE ;
+ return 0 ;
}
//______________________________________________________________________________
Int_t AliAnalysisTaskTaggedPhotons::GetFiducialArea(const Float_t * position)const{
}
}
//_____________________________________________________________________________
-void AliAnalysisTaskTaggedPhotons::FillPIDHistograms(const char * name, const AliAODPWG4Particle * p) const{
-
- FillHistogram(Form("%s_All_cent%d",name,fCentBin),p->Pt()) ;
- if(p->GetDispBit())
- FillHistogram(Form("%s_Disp_cent%d",name,fCentBin),p->Pt()) ;
- if(p->GetChargedBit())
- FillHistogram(Form("%s_CPV_cent%d",name,fCentBin),p->Pt()) ;
- if(p->GetDispBit() && p->GetChargedBit())
- FillHistogram(Form("%s_Both_cent%d",name,fCentBin),p->Pt()) ;
+void AliAnalysisTaskTaggedPhotons::FillPIDHistograms(const char * name, const AliCaloPhoton * p) const{
+
+ FillHistogram(Form("%s_All_cent%d",name,fCentBin),p->Pt(),p->GetWeight()) ;
+ if(p->IsDispOK())
+ FillHistogram(Form("%s_Disp_cent%d",name,fCentBin),p->Pt(),p->GetWeight()) ;
+ if(p->IsCPVOK())
+ FillHistogram(Form("%s_CPV_cent%d",name,fCentBin),p->Pt(),p->GetWeight()) ;
+ if(p->IsDispOK() && p->IsCPVOK())
+ FillHistogram(Form("%s_Both_cent%d",name,fCentBin),p->Pt(),p->GetWeight()) ;
}
//_____________________________________________________________________________
-void AliAnalysisTaskTaggedPhotons::FillPIDHistograms(const char * name, const AliAODPWG4Particle * p,Double_t x) const{
-
- FillHistogram(Form("%s_All_cent%d",name,fCentBin),x,p->Pt()) ;
- if(p->GetDispBit())
- FillHistogram(Form("%s_Disp_cent%d",name,fCentBin),x,p->Pt()) ;
- if(p->GetChargedBit())
- FillHistogram(Form("%s_CPV_cent%d",name,fCentBin),x,p->Pt()) ;
- if(p->GetDispBit() && p->GetChargedBit())
- FillHistogram(Form("%s_Both_cent%d",name,fCentBin),x,p->Pt()) ;
+void AliAnalysisTaskTaggedPhotons::FillPIDHistograms(const char * name, const AliCaloPhoton * p,Double_t x) const{
+
+ FillHistogram(Form("%s_All_cent%d",name,fCentBin),x,p->Pt(),p->GetWeight()) ;
+ if(p->IsDispOK())
+ FillHistogram(Form("%s_Disp_cent%d",name,fCentBin),x,p->Pt(),p->GetWeight()) ;
+ if(p->IsCPVOK())
+ FillHistogram(Form("%s_CPV_cent%d",name,fCentBin),x,p->Pt(),p->GetWeight()) ;
+ if(p->IsDispOK() && p->IsCPVOK())
+ FillHistogram(Form("%s_Both_cent%d",name,fCentBin),x,p->Pt(),p->GetWeight()) ;
}
//_____________________________________________________________________________
-void AliAnalysisTaskTaggedPhotons::FillPIDHistograms(const char * name, const AliAODPWG4Particle * p1,const AliAODPWG4Particle * p2,Double_t x) const{
-
- Double_t ptPi = (*(p1->Momentum()) + *(p2->Momentum())).Pt() ;
- FillHistogram(Form("%s_All_cent%d",name,fCentBin),x,ptPi) ;
- if(p1->GetDispBit() && p2->GetDispBit())
- FillHistogram(Form("%s_Disp_cent%d",name,fCentBin),x,ptPi) ;
- if(p1->GetChargedBit() && p2->GetChargedBit())
- FillHistogram(Form("%s_CPV_cent%d",name,fCentBin),x,ptPi) ;
- if(p1->GetDispBit() && p1->GetChargedBit() && p2->GetDispBit() && p2->GetChargedBit())
- FillHistogram(Form("%s_Both_cent%d",name,fCentBin),x,ptPi) ;
+void AliAnalysisTaskTaggedPhotons::FillPIDHistograms(const char * name, const AliCaloPhoton * p1,const AliCaloPhoton * p2,Double_t x, Bool_t /*isRe*/) const{
+
+ Double_t ptPi = (*p1 + *p2).Pt() ;
+ Double_t w=TMath::Sqrt(p1->GetWeight()*p2->GetWeight()) ;
+// if(isRe){
+// Int_t pdg=IsSameParent(p1,p2) ;
+// if(pdg!=0)
+// w=p1->GetWeight() ;
+// }
+ FillHistogram(Form("%s_All_cent%d",name,fCentBin),x,ptPi,w) ;
+ if(p1->IsDispOK() && p2->IsDispOK())
+ FillHistogram(Form("%s_Disp_cent%d",name,fCentBin),x,ptPi,w) ;
+ if(p1->IsCPVOK() && p2->IsCPVOK())
+ FillHistogram(Form("%s_CPV_cent%d",name,fCentBin),x,ptPi,w) ;
+ if(p1->IsDispOK() && p1->IsCPVOK() && p2->IsDispOK() && p2->IsCPVOK())
+ FillHistogram(Form("%s_Both_cent%d",name,fCentBin),x,ptPi,w) ;
}
//_____________________________________________________________________________
}
//_________________________________________________________________________________
-Int_t AliAnalysisTaskTaggedPhotons::EvalIsolation(TLorentzVector * ph){
+Int_t AliAnalysisTaskTaggedPhotons::EvalIsolation(TLorentzVector * ph, Bool_t isPhoton){
// Check if this particle is isolated.
//We use several cone radii and epsilons.
const Double_t coneR2=0.4 ;
const Double_t coneR3=0.5 ;
- const Double_t epsilon1=0.1 ;
- const Double_t epsilon2=0.05 ;
+ const Double_t cutEcone1=3. ;
+ const Double_t cutEcone2=4.5 ;
+ const Double_t cutEcone3=6.5 ;
if(!ph) return 0 ;
Double_t eCone1 = 0;
Double_t eCone2 = 0;
Double_t eCone3 = 0;
+ //Cross-check, energy in cone perpr to photon
+ Double_t eP1Cone1 = 0;
+ Double_t eP1Cone2 = 0;
+ Double_t eP1Cone3 = 0;
+ Double_t eP2Cone1 = 0;
+ Double_t eP2Cone2 = 0;
+ Double_t eP2Cone3 = 0;
+
Double_t phiTrig = ph->Phi();
Double_t etaTrig = ph->Eta();
Int_t n=fTrackEvent->GetEntriesFast() ;
for(Int_t itr=0; itr<n; itr++){
- AliAODPWG4Particle * track = (AliAODPWG4Particle*)fTrackEvent->At(itr) ;
+ AliCaloPhoton * track = (AliCaloPhoton*)fTrackEvent->At(itr) ;
Double_t deleta = etaTrig - track->Eta() ;
Double_t delphi = phiTrig - track->Phi() ;
while(delphi<-TMath::Pi()) delphi+=TMath::TwoPi() ;
while(delphi>TMath::Pi()) delphi-=TMath::TwoPi() ;
+ //Perp cones
+ Double_t delphiP1 = phiTrig +TMath::PiOver2()- track->Phi() ;
+ while(delphiP1<-TMath::Pi()) delphiP1+=TMath::TwoPi() ;
+ while(delphiP1>TMath::Pi()) delphiP1-=TMath::TwoPi() ;
+ Double_t delphiP2 = phiTrig -TMath::PiOver2()- track->Phi() ;
+ while(delphiP2<-TMath::Pi()) delphiP2+=TMath::TwoPi() ;
+ while(delphiP2>TMath::Pi()) delphiP2-=TMath::TwoPi() ;
+
+
Double_t dr = TMath::Sqrt(deleta * deleta + delphi * delphi);
+ Double_t drP1 = TMath::Sqrt(deleta * deleta + delphiP1 * delphiP1);
+ Double_t drP2 = TMath::Sqrt(deleta * deleta + delphiP2 * delphiP2);
if(dr<coneR3){
eCone3+=track->Pt() ;
}
}
}
+
+ if(drP1<coneR3){
+ eP1Cone3+=track->Pt() ;
+ if(drP1<coneR2){
+ eP1Cone2+=track->Pt() ;
+ if(drP1<coneR1){
+ eP1Cone1+=track->Pt() ;
+ }
+ }
+ }
+
+ if(drP2<coneR3){
+ eP2Cone3+=track->Pt() ;
+ if(drP2<coneR2){
+ eP2Cone2+=track->Pt() ;
+ if(drP2<coneR1){
+ eP2Cone1+=track->Pt() ;
+ }
+ }
+ }
}
//Fill QA histgams
Double_t ptTrig=ph->Pt() ;
+ if(isPhoton){
FillHistogram(Form("QA_Cone1_Tracks_cent%d",fCentBin),ptTrig,eCone1) ;
FillHistogram(Form("QA_Cone2_Tracks_cent%d",fCentBin),ptTrig,eCone2) ;
FillHistogram(Form("QA_Cone3_Tracks_cent%d",fCentBin),ptTrig,eCone3) ;
-
+ FillHistogram(Form("QA_PCone1_Tracks_cent%d",fCentBin),ptTrig,eP1Cone1) ;
+ FillHistogram(Form("QA_PCone2_Tracks_cent%d",fCentBin),ptTrig,eP1Cone2) ;
+ FillHistogram(Form("QA_PCone3_Tracks_cent%d",fCentBin),ptTrig,eP1Cone3) ;
+ FillHistogram(Form("QA_PCone1_Tracks_cent%d",fCentBin),ptTrig,eP2Cone1) ;
+ FillHistogram(Form("QA_PCone2_Tracks_cent%d",fCentBin),ptTrig,eP2Cone2) ;
+ FillHistogram(Form("QA_PCone3_Tracks_cent%d",fCentBin),ptTrig,eP2Cone3) ;
+ }
+ else{
+ FillHistogram(Form("QA_Pi0Cone1_Tracks_cent%d",fCentBin),ptTrig,eCone1) ;
+ FillHistogram(Form("QA_Pi0Cone2_Tracks_cent%d",fCentBin),ptTrig,eCone2) ;
+ FillHistogram(Form("QA_Pi0Cone3_Tracks_cent%d",fCentBin),ptTrig,eCone3) ;
+ FillHistogram(Form("QA_Pi0PCone1_Tracks_cent%d",fCentBin),ptTrig,eP1Cone1) ;
+ FillHistogram(Form("QA_Pi0PCone2_Tracks_cent%d",fCentBin),ptTrig,eP1Cone2) ;
+ FillHistogram(Form("QA_Pi0PCone3_Tracks_cent%d",fCentBin),ptTrig,eP1Cone3) ;
+ FillHistogram(Form("QA_Pi0PCone1_Tracks_cent%d",fCentBin),ptTrig,eP2Cone1) ;
+ FillHistogram(Form("QA_Pi0PCone2_Tracks_cent%d",fCentBin),ptTrig,eP2Cone2) ;
+ FillHistogram(Form("QA_Pi0PCone3_Tracks_cent%d",fCentBin),ptTrig,eP2Cone3) ;
+
+ }
//Fill Bits
- Int_t iCone1E1 = (epsilon1*ptTrig > eCone1) ;
- Int_t iCone2E1 = (epsilon1*ptTrig > eCone2) ;
- Int_t iCone3E1 = (epsilon1*ptTrig > eCone3) ;
+ Int_t iCone1E1 = (cutEcone1 > eCone1) ;
+ Int_t iCone2E1 = (cutEcone2 > eCone2) ;
+ Int_t iCone3E1 = (cutEcone3 > eCone3) ;
- Int_t iCone1E2 = (epsilon2*ptTrig > eCone1) ;
- Int_t iCone2E2 = (epsilon2*ptTrig > eCone2) ;
- Int_t iCone3E2 = (epsilon2*ptTrig > eCone3) ;
+ Int_t iCone1E2 = (1.5*cutEcone1 > eCone1) ;
+ Int_t iCone2E2 = (1.5*cutEcone2 > eCone2) ;
+ Int_t iCone3E2 = (1.5*cutEcone3 > eCone3) ;
Int_t isolation= iCone1E1+ 2*iCone2E1 +4*iCone3E1+
return isolation ;
}
//_________________________________________________________________________________
-Bool_t AliAnalysisTaskTaggedPhotons::TestPID(Int_t iPID, AliAODPWG4Particle* part){
+Bool_t AliAnalysisTaskTaggedPhotons::TestPID(Int_t iPID, AliCaloPhoton* part){
// //Checks PID of particle
if(!part) return kFALSE ;
if(iPID==0) return kTRUE ;
- if(iPID==1) return part->GetDispBit() ;
- if(iPID==2) return part->GetChargedBit() ;
- if(iPID==3) return part->GetDispBit() && part->GetChargedBit() ;
+ if(iPID==1) return part->IsDispOK() ;
+ if(iPID==2) return part->IsCPVOK() ;
+ if(iPID==3) return part->IsDispOK() && part->IsCPVOK() ;
return kFALSE ;
}
//_________________________________________________________________________________
-Double_t AliAnalysisTaskTaggedPhotons::PrimaryParticleWeight(AliAODMCParticle * /*particle*/){
- return 1;
+Double_t AliAnalysisTaskTaggedPhotons::PrimaryParticleWeight(AliAODMCParticle * particle){
+ if(!fIsMC) //This is real data
+ return 1;
+ //Classify parent at vertex
+ //Introduce for eta and pi0 weights
+
+ Double_t r2=particle->Xv()*particle->Xv()+particle->Yv()*particle->Yv() ;
+ Int_t mother = particle->GetMother() ;
+ while(mother>-1){
+ if(r2<1.)
+ break ;
+ particle = (AliAODMCParticle*) fStack->At(mother);
+ mother = particle->GetMother() ;
+ r2=particle->Xv()*particle->Xv()+particle->Yv()*particle->Yv() ;
+ }
+
+ //Particle within 1 cm from the virtex
+ Int_t pdg = particle->GetPdgCode() ;
+ Double_t x = particle->Pt() ;
+ mother = particle->GetMother() ;
+ while(TMath::Abs(pdg)<100){//gamma, electrons, muons
+ if(mother>-1){
+ AliAODMCParticle * tmpP=(AliAODMCParticle*)fStack->At(mother) ;
+ pdg=tmpP->GetPdgCode() ;
+ x = tmpP->Pt() ;
+ mother = tmpP->GetMother() ;
+ }
+ else{ //direct photon/electron....
+ return 1.;
+ }
+ }
+ if(pdg == 111){
+ //Pi0
+ if(x<1) return 1. ;
+ else return fWeightParamPi0[0]*(TMath::Power(x,fWeightParamPi0[1])*
+ (1.+fWeightParamPi0[2]*x+fWeightParamPi0[3]*x*x)/
+ (1.+fWeightParamPi0[4]*x+fWeightParamPi0[5]*x*x) +fWeightParamPi0[6]) ;
+ }
+ if(pdg == 221){
+ //Eta - same same shape, but yield 0.48 and Br(eta->2gamma)
+ Double_t norm=0.48 * 0.3943;
+ if(x<1) return norm ;
+ else return norm*fWeightParamPi0[0]*(TMath::Power(x,fWeightParamPi0[1])*
+ (1.+fWeightParamPi0[2]*x+fWeightParamPi0[3]*x*x)/
+ (1.+fWeightParamPi0[4]*x+fWeightParamPi0[5]*x*x) +fWeightParamPi0[6]) ;
+ }
+ return 1. ;
+}
+//_________________________________________________________________________________
+void AliAnalysisTaskTaggedPhotons::SetPi0WeightParameters(TArrayD * ar){
+
+ for(Int_t i=0; i<7; i++){ //Array range
+ if(ar->GetSize()>i) fWeightParamPi0[i]=ar->At(i) ;
+ else fWeightParamPi0[i]=0.;
+ }
+ //normalize to obtain smooth transition at 1 GeV
+ Double_t x=1. ;
+ fWeightParamPi0[0]=1./(TMath::Power(x,fWeightParamPi0[1])*
+ (1.+fWeightParamPi0[2]*x+fWeightParamPi0[3]*x*x)/
+ (1.+fWeightParamPi0[4]*x+fWeightParamPi0[5]*x*x) +fWeightParamPi0[6]) ;
+
}
//___________________________________________________________________________
//Check if this channel belogs to the good ones
if(mod>4 || mod<1){
- return kTRUE ;
+ return kFALSE ;
}
if(!fPHOSBadMap[mod]){
return kTRUE ;
else
return kTRUE ;
}
-
-
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff