From 8bdca7f1e883b70f0567d11ce1677f035d06af60 Mon Sep 17 00:00:00 2001 From: fbock Date: Fri, 27 Jun 2014 15:20:32 +0200 Subject: [PATCH] - added new tasks for combination of EMCAL and PCM for pi0 and photons --- PWGGA/CMakelibPWGGAGammaConv.pkg | 4 +- .../AliAnalysisTaskGammaConvCalo.cxx | 2913 +++++++++++++++++ .../GammaConv/AliAnalysisTaskGammaConvCalo.h | 300 ++ PWGGA/GammaConv/AliCaloPhotonCuts.cxx | 1268 +++++++ PWGGA/GammaConv/AliCaloPhotonCuts.h | 183 ++ .../macros/AddTask_GammaConvCalo_PbPb.C | 205 ++ .../macros/AddTask_GammaConvCalo_pPb.C | 181 + .../macros/AddTask_GammaConvCalo_pp.C | 162 + .../macros/AddTask_GammaConvV1_pPb.C | 16 +- PWGGA/PWGGAGammaConvLinkDef.h | 25 +- 10 files changed, 5236 insertions(+), 21 deletions(-) create mode 100644 PWGGA/GammaConv/AliAnalysisTaskGammaConvCalo.cxx create mode 100644 PWGGA/GammaConv/AliAnalysisTaskGammaConvCalo.h create mode 100644 PWGGA/GammaConv/AliCaloPhotonCuts.cxx create mode 100644 PWGGA/GammaConv/AliCaloPhotonCuts.h create mode 100644 PWGGA/GammaConv/macros/AddTask_GammaConvCalo_PbPb.C create mode 100644 PWGGA/GammaConv/macros/AddTask_GammaConvCalo_pPb.C create mode 100644 PWGGA/GammaConv/macros/AddTask_GammaConvCalo_pp.C diff --git a/PWGGA/CMakelibPWGGAGammaConv.pkg b/PWGGA/CMakelibPWGGAGammaConv.pkg index 441f49aaa63..8ab2d123634 100644 --- a/PWGGA/CMakelibPWGGAGammaConv.pkg +++ b/PWGGA/CMakelibPWGGAGammaConv.pkg @@ -33,6 +33,7 @@ set ( SRCS GammaConv/AliAODConversionPhoton.cxx GammaConv/AliKFConversionPhoton.cxx GammaConv/AliKFConversionMother.cxx + GammaConv/AliCaloPhotonCuts.cxx GammaConv/AliConversionCuts.cxx GammaConv/AliConversionSelection.cxx GammaConv/AliConversionMesonCuts.cxx @@ -57,6 +58,7 @@ set ( SRCS GammaConv/AliPrimaryPionSelector.cxx GammaConv/AliPrimaryPionCuts.cxx GammaConv/AliAnalysisTaskEtaToPiPlPiMiGamma.cxx + GammaConv/AliAnalysisTaskGammaConvCalo.cxx ) string ( REPLACE ".cxx" ".h" HDRS "${SRCS}" ) @@ -65,4 +67,4 @@ set ( DHDR PWGGAGammaConvLinkDef.h) string ( REPLACE ".cxx" ".h" EXPORT "${SRCS}" ) -set ( EINCLUDE PWGGA/GammaConv PWG/TRD CORRFW STEER/AOD STEER/ESD STEER/STEERBase ANALYSIS) +set ( EINCLUDE PWGGA/PWGGAUtils PWGGA/GammaConv PWG/TRD PWG/EMCAL EMCAL OADB STEER/STEER CORRFW STEER/AOD STEER/ESD STEER/STEERBase ANALYSIS) diff --git a/PWGGA/GammaConv/AliAnalysisTaskGammaConvCalo.cxx b/PWGGA/GammaConv/AliAnalysisTaskGammaConvCalo.cxx new file mode 100644 index 00000000000..85fe4c8cef8 --- /dev/null +++ b/PWGGA/GammaConv/AliAnalysisTaskGammaConvCalo.cxx @@ -0,0 +1,2913 @@ +/************************************************************************** + * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * * + * Author: Baldo Sahlmueller, Friederike Bock * + * Version 1.0 * + * * + * * + * 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. * + **************************************************************************/ + +////////////////////////////////////////////////////////////////// +//---------------------------------------------------------------- +// Class used to do analysis on conversion photons + calo photons +//---------------------------------------------------------------- +////////////////////////////////////////////////////////////////// +#include "TChain.h" +#include "TTree.h" +#include "TBranch.h" +#include "TFile.h" +#include "TH1F.h" +#include "TH2F.h" +#include "TH3F.h" +#include "THnSparse.h" +#include "TCanvas.h" +#include "TNtuple.h" +#include "AliAnalysisTask.h" +#include "AliAnalysisManager.h" +#include "AliESDEvent.h" +#include "AliESDInputHandler.h" +#include "AliMCEventHandler.h" +#include "AliMCEvent.h" +#include "AliMCParticle.h" +#include "AliCentrality.h" +#include "AliESDVZERO.h" +#include "AliESDpid.h" +#include "AliAnalysisTaskGammaConvCalo.h" +#include "AliVParticle.h" +#include "AliESDtrack.h" +#include "AliESDtrackCuts.h" +#include "AliKFVertex.h" +#include "AliV0ReaderV1.h" +#include "AliGenCocktailEventHeader.h" +#include "AliConversionAODBGHandlerRP.h" +#include "AliAODMCParticle.h" +#include "AliAODMCHeader.h" +#include "AliEventplane.h" +#include "AliAnalysisTaskEMCALClusterizeFast.h" +#include "AliAODEvent.h" +#include "AliAODInputHandler.h" +#include "AliESDEvent.h" +#include "AliESDInputHandler.h" +#include "AliInputEventHandler.h" + +ClassImp(AliAnalysisTaskGammaConvCalo) + +//________________________________________________________________________ +AliAnalysisTaskGammaConvCalo::AliAnalysisTaskGammaConvCalo(): AliAnalysisTaskSE(), + fV0Reader(NULL), + fBGHandler(NULL), + fBGHandlerRP(NULL), + fBGClusHandler(NULL), + fBGClusHandlerRP(NULL), + fInputEvent(NULL), + fMCEvent(NULL), + fMCStack(NULL), + fCutFolder(NULL), + fESDList(NULL), + fBackList(NULL), + fMotherList(NULL), + fPhotonDCAList(NULL), + fMesonDCAList(NULL), + fTrueList(NULL), + fMCList(NULL), + fHeaderNameList(NULL), + fTagOutputList(NULL), + fOutputContainer(NULL), + fReaderGammas(NULL), + fGammaCandidates(NULL), + fClusterCandidates(NULL), + fCutArray(NULL), + fConversionCuts(NULL), + fClusterCutArray(NULL), + fCaloPhotonCuts(NULL), + fMesonCutArray(NULL), + fMesonCuts(NULL), + fHistoConvGammaPt(NULL), + fHistoConvGammaR(NULL), + fHistoConvGammaEta(NULL), + fTreeConvGammaPtDcazCat(NULL), + fPtGamma(0), + fDCAzPhoton(0), + fRConvPhoton(0), + fEtaPhoton(0), + fCharCatPhoton(0), + fCharPhotonMCInfo(0), + fHistoMotherInvMassPt(NULL), + fSparseMotherInvMassPtZM(NULL), + fHistoMotherBackInvMassPt(NULL), + fSparseMotherBackInvMassPtZM(NULL), + fHistoMotherInvMassEalpha(NULL), + fHistoMotherPi0PtY(NULL), + fHistoMotherEtaPtY(NULL), + fHistoMotherPi0PtAlpha(NULL), + fHistoMotherEtaPtAlpha(NULL), + fHistoMotherPi0PtOpenAngle(NULL), + fHistoMotherEtaPtOpenAngle(NULL), + fTreeMesonsInvMassPtDcazMinDcazMaxFlag(NULL), + fInvMass(0), + fPt(0), + fDCAzGammaMin(0), + fDCAzGammaMax(0), + fCharFlag(0), + fCharMesonMCInfo(0), + fHistoConvGammaUntagged(NULL), + fHistoConvGammaTagged(NULL), + fHistoConvGammaPi0Tagged(NULL), + fHistoConvGammaEtaTagged(NULL), + fHistoPhotonPairAll(NULL), + fHistoPhotonPairAllGam(NULL), + fHistoClusGammaPt(NULL), + fHistoMCHeaders(NULL), + fHistoMCAllGammaPt(NULL), + fHistoMCDecayGammaPi0Pt(NULL), + fHistoMCDecayGammaRhoPt(NULL), + fHistoMCDecayGammaEtaPt(NULL), + fHistoMCDecayGammaOmegaPt(NULL), + fHistoMCDecayGammaEtapPt(NULL), + fHistoMCDecayGammaPhiPt(NULL), + fHistoMCDecayGammaSigmaPt(NULL), + fHistoMCConvGammaPt(NULL), + fHistoMCConvGammaR(NULL), + fHistoMCConvGammaEta(NULL), + fHistoMCPi0Pt(NULL), + fHistoMCPi0WOWeightPt(NULL), + fHistoMCEtaPt(NULL), + fHistoMCEtaWOWeightPt(NULL), + fHistoMCPi0InAccPt(NULL), + fHistoMCEtaInAccPt(NULL), + fHistoMCPi0PtY(NULL), + fHistoMCEtaPtY(NULL), + fHistoMCK0sPt(NULL), + fHistoMCK0sWOWeightPt(NULL), + fHistoMCK0sPtY(NULL), + fHistoMCSecPi0PtvsSource(NULL), + fHistoMCSecPi0Source(NULL), + fHistoMCSecEtaPt(NULL), + fHistoMCSecEtaSource(NULL), + fHistoTrueMotherInvMassPt(NULL), + fHistoTruePrimaryMotherInvMassPt(NULL), + fHistoTruePrimaryMotherW0WeightingInvMassPt(NULL), + fProfileTruePrimaryMotherWeightsInvMassPt(NULL), + fHistoTruePrimaryPi0MCPtResolPt(NULL), + fHistoTruePrimaryEtaMCPtResolPt(NULL), + fHistoTrueSecondaryMotherInvMassPt(NULL), + fHistoTrueSecondaryMotherFromK0sInvMassPt(NULL), + fHistoTrueK0sWithPi0DaughterMCPt(NULL), + fHistoTrueSecondaryMotherFromEtaInvMassPt(NULL), + fHistoTrueEtaWithPi0DaughterMCPt(NULL), + fHistoTrueSecondaryMotherFromLambdaInvMassPt(NULL), + fHistoTrueLambdaWithPi0DaughterMCPt(NULL), + fHistoTrueBckGGInvMassPt(NULL), + fHistoTrueBckContInvMassPt(NULL), + fHistoTruePi0PtY(NULL), + fHistoTrueEtaPtY(NULL), + fHistoTruePi0PtAlpha(NULL), + fHistoTrueEtaPtAlpha(NULL), + fHistoTruePi0PtOpenAngle(NULL), + fHistoTrueEtaPtOpenAngle(NULL), + fHistoTrueMotherDalitzInvMassPt(NULL), + fHistoTrueConvGammaPt(NULL), + fHistoTrueConvPi0GammaPt(NULL), + fHistoTrueConvGammaEta(NULL), + fHistoCombinatorialPt(NULL), + fHistoTruePrimaryConvGammaPt(NULL), + fHistoTruePrimaryConvGammaESDPtMCPt(NULL), + fHistoTrueSecondaryConvGammaPt(NULL), + fHistoTrueSecondaryConvGammaFromXFromK0sPt(NULL), + fHistoTrueSecondaryConvGammaFromXFromLambdaPt(NULL), + fHistoTrueDalitzPsiPairDeltaPhi(NULL), + fHistoTrueGammaPsiPairDeltaPhi(NULL), + fHistoTrueClusGammaPt(NULL), + fHistoTruePrimaryClusGammaPt(NULL), + fHistoTruePrimaryClusGammaESDPtMCPt(NULL), + fHistoNEvents(NULL), + fHistoNGoodESDTracks(NULL), + fHistoNGammaCandidates(NULL), + fHistoNGoodESDTracksVsNGammaCanditates(NULL), + fHistoNV0Tracks(NULL), + fProfileEtaShift(NULL), + fEventPlaneAngle(-100), + fRandom(0), + fNGammaCandidates(0), + fUnsmearedPx(NULL), + fUnsmearedPy(NULL), + fUnsmearedPz(NULL), + fUnsmearedE(NULL), + fMCStackPos(NULL), + fMCStackNeg(NULL), + fESDArrayPos(NULL), + fESDArrayNeg(NULL), + fnCuts(0), + fiCut(0), + fMoveParticleAccordingToVertex(kTRUE), + fIsHeavyIon(0), + fDoMesonAnalysis(kTRUE), + fDoMesonQA(0), + fDoPhotonQA(0), + fIsFromMBHeader(kTRUE), + fIsMC(kFALSE), + fMinE(0.1), + fNminCells(2), + fEMCm02cut(0.5) +{ + +} + +//________________________________________________________________________ +AliAnalysisTaskGammaConvCalo::AliAnalysisTaskGammaConvCalo(const char *name): + AliAnalysisTaskSE(name), + fV0Reader(NULL), + fBGHandler(NULL), + fBGHandlerRP(NULL), + fBGClusHandler(NULL), + fBGClusHandlerRP(NULL), + fInputEvent(NULL), + fMCEvent(NULL), + fMCStack(NULL), + fCutFolder(NULL), + fESDList(NULL), + fBackList(NULL), + fMotherList(NULL), + fPhotonDCAList(NULL), + fMesonDCAList(NULL), + fTrueList(NULL), + fMCList(NULL), + fHeaderNameList(NULL), + fTagOutputList(NULL), + fOutputContainer(0), + fReaderGammas(NULL), + fGammaCandidates(NULL), + fClusterCandidates(NULL), + fCutArray(NULL), + fConversionCuts(NULL), + fClusterCutArray(NULL), + fCaloPhotonCuts(NULL), + fMesonCutArray(NULL), + fMesonCuts(NULL), + fHistoConvGammaPt(NULL), + fHistoConvGammaR(NULL), + fHistoConvGammaEta(NULL), + fTreeConvGammaPtDcazCat(NULL), + fPtGamma(0), + fDCAzPhoton(0), + fRConvPhoton(0), + fEtaPhoton(0), + fCharCatPhoton(0), + fCharPhotonMCInfo(0), + fHistoMotherInvMassPt(NULL), + fSparseMotherInvMassPtZM(NULL), + fHistoMotherBackInvMassPt(NULL), + fSparseMotherBackInvMassPtZM(NULL), + fHistoMotherInvMassEalpha(NULL), + fHistoMotherPi0PtY(NULL), + fHistoMotherEtaPtY(NULL), + fHistoMotherPi0PtAlpha(NULL), + fHistoMotherEtaPtAlpha(NULL), + fHistoMotherPi0PtOpenAngle(NULL), + fHistoMotherEtaPtOpenAngle(NULL), + fTreeMesonsInvMassPtDcazMinDcazMaxFlag(NULL), + fInvMass(0), + fPt(0), + fDCAzGammaMin(0), + fDCAzGammaMax(0), + fCharFlag(0), + fCharMesonMCInfo(0), + fHistoConvGammaUntagged(NULL), + fHistoConvGammaTagged(NULL), + fHistoConvGammaPi0Tagged(NULL), + fHistoConvGammaEtaTagged(NULL), + fHistoPhotonPairAll(NULL), + fHistoPhotonPairAllGam(NULL), + fHistoClusGammaPt(NULL), + fHistoMCHeaders(NULL), + fHistoMCAllGammaPt(NULL), + fHistoMCDecayGammaPi0Pt(NULL), + fHistoMCDecayGammaRhoPt(NULL), + fHistoMCDecayGammaEtaPt(NULL), + fHistoMCDecayGammaOmegaPt(NULL), + fHistoMCDecayGammaEtapPt(NULL), + fHistoMCDecayGammaPhiPt(NULL), + fHistoMCDecayGammaSigmaPt(NULL), + fHistoMCConvGammaPt(NULL), + fHistoMCConvGammaR(NULL), + fHistoMCConvGammaEta(NULL), + fHistoMCPi0Pt(NULL), + fHistoMCPi0WOWeightPt(NULL), + fHistoMCEtaPt(NULL), + fHistoMCEtaWOWeightPt(NULL), + fHistoMCPi0InAccPt(NULL), + fHistoMCEtaInAccPt(NULL), + fHistoMCPi0PtY(NULL), + fHistoMCEtaPtY(NULL), + fHistoMCK0sPt(NULL), + fHistoMCK0sWOWeightPt(NULL), + fHistoMCK0sPtY(NULL), + fHistoMCSecPi0PtvsSource(NULL), + fHistoMCSecPi0Source(NULL), + fHistoMCSecEtaPt(NULL), + fHistoMCSecEtaSource(NULL), + fHistoTrueMotherInvMassPt(NULL), + fHistoTruePrimaryMotherInvMassPt(NULL), + fHistoTruePrimaryMotherW0WeightingInvMassPt(NULL), + fProfileTruePrimaryMotherWeightsInvMassPt(NULL), + fHistoTruePrimaryPi0MCPtResolPt(NULL), + fHistoTruePrimaryEtaMCPtResolPt(NULL), + fHistoTrueSecondaryMotherInvMassPt(NULL), + fHistoTrueSecondaryMotherFromK0sInvMassPt(NULL), + fHistoTrueK0sWithPi0DaughterMCPt(NULL), + fHistoTrueSecondaryMotherFromEtaInvMassPt(NULL), + fHistoTrueEtaWithPi0DaughterMCPt(NULL), + fHistoTrueSecondaryMotherFromLambdaInvMassPt(NULL), + fHistoTrueLambdaWithPi0DaughterMCPt(NULL), + fHistoTrueBckGGInvMassPt(NULL), + fHistoTrueBckContInvMassPt(NULL), + fHistoTruePi0PtY(NULL), + fHistoTrueEtaPtY(NULL), + fHistoTruePi0PtAlpha(NULL), + fHistoTrueEtaPtAlpha(NULL), + fHistoTruePi0PtOpenAngle(NULL), + fHistoTrueEtaPtOpenAngle(NULL), + fHistoTrueMotherDalitzInvMassPt(NULL), + fHistoTrueConvGammaPt(NULL), + fHistoTrueConvPi0GammaPt(NULL), + fHistoTrueConvGammaEta(NULL), + fHistoCombinatorialPt(NULL), + fHistoTruePrimaryConvGammaPt(NULL), + fHistoTruePrimaryConvGammaESDPtMCPt(NULL), + fHistoTrueSecondaryConvGammaPt(NULL), + fHistoTrueSecondaryConvGammaFromXFromK0sPt(NULL), + fHistoTrueSecondaryConvGammaFromXFromLambdaPt(NULL), + fHistoTrueDalitzPsiPairDeltaPhi(NULL), + fHistoTrueGammaPsiPairDeltaPhi(NULL), + fHistoTrueClusGammaPt(NULL), + fHistoTruePrimaryClusGammaPt(NULL), + fHistoTruePrimaryClusGammaESDPtMCPt(NULL), + fHistoNEvents(NULL), + fHistoNGoodESDTracks(NULL), + fHistoNGammaCandidates(NULL), + fHistoNGoodESDTracksVsNGammaCanditates(NULL), + fHistoNV0Tracks(NULL), + fProfileEtaShift(NULL), + fEventPlaneAngle(-100), + fRandom(0), + fNGammaCandidates(0), + fUnsmearedPx(NULL), + fUnsmearedPy(NULL), + fUnsmearedPz(NULL), + fUnsmearedE(NULL), + fMCStackPos(NULL), + fMCStackNeg(NULL), + fESDArrayPos(NULL), + fESDArrayNeg(NULL), + fnCuts(0), + fiCut(0), + fMoveParticleAccordingToVertex(kTRUE), + fIsHeavyIon(0), + fDoMesonAnalysis(kTRUE), + fDoMesonQA(0), + fDoPhotonQA(0), + fIsFromMBHeader(kTRUE), + fIsMC(kFALSE), + fMinE(0.1), + fNminCells(2), + fEMCm02cut(0.5) +{ + // Define output slots here + DefineOutput(1, TList::Class()); +} + +AliAnalysisTaskGammaConvCalo::~AliAnalysisTaskGammaConvCalo() +{ + if(fGammaCandidates){ + delete fGammaCandidates; + fGammaCandidates = 0x0; + } + if(fClusterCandidates){ + delete fClusterCandidates; + fClusterCandidates = 0x0; + } + if(fBGHandler){ + delete[] fBGHandler; + fBGHandler = 0x0; + } + if(fBGHandlerRP){ + delete[] fBGHandlerRP; + fBGHandlerRP = 0x0; + } + if(fBGClusHandler){ + delete[] fBGClusHandler; + fBGClusHandler = 0x0; + } + if(fBGClusHandlerRP){ + delete[] fBGClusHandlerRP; + fBGClusHandlerRP = 0x0; + } +} +//___________________________________________________________ +void AliAnalysisTaskGammaConvCalo::InitBack(){ + + const Int_t nDim = 4; + Int_t nBins[nDim] = {800,250,7,4}; + Double_t xMin[nDim] = {0,0, 0,0}; + Double_t xMax[nDim] = {0.8,25,7,4}; + + fSparseMotherInvMassPtZM = new THnSparseF*[fnCuts]; + fSparseMotherBackInvMassPtZM = new THnSparseF*[fnCuts]; + + fBGHandler = new AliGammaConversionAODBGHandler*[fnCuts]; + fBGHandlerRP = new AliConversionAODBGHandlerRP*[fnCuts]; + + fBGClusHandler = new AliGammaConversionAODBGHandler*[fnCuts]; + fBGClusHandlerRP = new AliConversionAODBGHandlerRP*[fnCuts]; + + for(Int_t iCut = 0; iCutAt(iCut))->DoBGCalculation()){ + TString cutstring = ((AliConversionCuts*)fCutArray->At(iCut))->GetCutNumber(); + TString cutstringCalo = ((AliCaloPhotonCuts*)fClusterCutArray->At(iCut))->GetCutNumber(); + TString cutstringMeson = ((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->GetCutNumber(); + + Int_t collisionSystem = atoi((TString)(((AliConversionCuts*)fCutArray->At(iCut))->GetCutNumber())(0,1)); + Int_t centMin = atoi((TString)(((AliConversionCuts*)fCutArray->At(iCut))->GetCutNumber())(1,1)); + Int_t centMax = atoi((TString)(((AliConversionCuts*)fCutArray->At(iCut))->GetCutNumber())(2,1)); + + if(collisionSystem == 1 || collisionSystem == 2 || + collisionSystem == 5 || collisionSystem == 8 || + collisionSystem == 9){ + centMin = centMin*10; + centMax = centMax*10; + if(centMax ==0 && centMax!=centMin) centMax=100; + } else if(collisionSystem == 3 || collisionSystem == 6){ + centMin = centMin*5; + centMax = centMax*5; + } else if(collisionSystem == 4 || collisionSystem == 7){ + centMin = ((centMin*5)+45); + centMax = ((centMax*5)+45); + } + + fBackList[iCut] = new TList(); + fBackList[iCut]->SetName(Form("%s_%s_%s Back histograms",cutstring.Data(),cutstringCalo.Data(), cutstringMeson.Data())); + fBackList[iCut]->SetOwner(kTRUE); + fCutFolder[iCut]->Add(fBackList[iCut]); + + fSparseMotherBackInvMassPtZM[iCut] = new THnSparseF("Back_Back_InvMass_Pt_z_m","Back_Back_InvMass_Pt_z_m",nDim,nBins,xMin,xMax); + fBackList[iCut]->Add(fSparseMotherBackInvMassPtZM[iCut]); + + fMotherList[iCut] = new TList(); + fMotherList[iCut]->SetName(Form("%s_%s_%s Mother histograms",cutstring.Data(),cutstringCalo.Data(),cutstringMeson.Data())); + fMotherList[iCut]->SetOwner(kTRUE); + fCutFolder[iCut]->Add(fMotherList[iCut]); + + fSparseMotherInvMassPtZM[iCut] = new THnSparseF("Back_Mother_InvMass_Pt_z_m","Back_Mother_InvMass_Pt_z_m",nDim,nBins,xMin,xMax); + fMotherList[iCut]->Add(fSparseMotherInvMassPtZM[iCut]); + + if(((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->BackgroundHandlerType() == 0){ + fBGHandler[iCut] = new AliGammaConversionAODBGHandler( + collisionSystem,centMin,centMax, + ((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->GetNumberOfBGEvents(), + ((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->UseTrackMultiplicity()); + fBGClusHandler[iCut] = new AliGammaConversionAODBGHandler( + collisionSystem,centMin,centMax, + ((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->GetNumberOfBGEvents(), + ((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->UseTrackMultiplicity()); + fBGHandlerRP[iCut] = NULL; + } else{ + fBGHandlerRP[iCut] = new AliConversionAODBGHandlerRP( + ((AliConversionCuts*)fCutArray->At(fiCut))->IsHeavyIon(), + ((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->UseTrackMultiplicity(), + ((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->GetNumberOfBGEvents()); + fBGClusHandlerRP[iCut] = new AliConversionAODBGHandlerRP( + ((AliConversionCuts*)fCutArray->At(fiCut))->IsHeavyIon(), + ((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->UseTrackMultiplicity(), + ((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->GetNumberOfBGEvents()); + fBGHandler[iCut] = NULL; + } + } + } +} +//________________________________________________________________________ +void AliAnalysisTaskGammaConvCalo::UserCreateOutputObjects(){ + + // Create histograms + if(fOutputContainer != NULL){ + delete fOutputContainer; + fOutputContainer = NULL; + } + if(fOutputContainer == NULL){ + fOutputContainer = new TList(); + fOutputContainer->SetOwner(kTRUE); + } + + // Array of current cut's gammas + fGammaCandidates = new TList(); + fClusterCandidates = new TList(); + + fCutFolder = new TList*[fnCuts]; + fESDList = new TList*[fnCuts]; + fBackList = new TList*[fnCuts]; + fMotherList = new TList*[fnCuts]; + fHistoNEvents = new TH1I*[fnCuts]; + fHistoNGoodESDTracks = new TH1I*[fnCuts]; + fHistoNGammaCandidates = new TH1I*[fnCuts]; + fHistoNGoodESDTracksVsNGammaCanditates = new TH2F*[fnCuts]; + fHistoNV0Tracks = new TH1I*[fnCuts]; + fProfileEtaShift = new TProfile*[fnCuts]; + fHistoConvGammaPt = new TH1F*[fnCuts]; + + if (fDoPhotonQA == 2){ + fPhotonDCAList = new TList*[fnCuts]; + fTreeConvGammaPtDcazCat = new TTree*[fnCuts]; + } + if (fDoPhotonQA > 0){ + fHistoConvGammaR = new TH1F*[fnCuts]; + fHistoConvGammaEta = new TH1F*[fnCuts]; + } + + if(fDoMesonAnalysis){ + fHistoMotherInvMassPt = new TH2F*[fnCuts]; + fHistoMotherBackInvMassPt = new TH2F*[fnCuts]; + fHistoMotherInvMassEalpha = new TH2F*[fnCuts]; + if (fDoMesonQA == 2){ + fMesonDCAList = new TList*[fnCuts]; + fTreeMesonsInvMassPtDcazMinDcazMaxFlag = new TTree*[fnCuts]; + } + if (fDoMesonQA > 0){ + fHistoMotherPi0PtY = new TH2F*[fnCuts]; + fHistoMotherEtaPtY = new TH2F*[fnCuts]; + fHistoMotherPi0PtAlpha = new TH2F*[fnCuts]; + fHistoMotherEtaPtAlpha = new TH2F*[fnCuts]; + fHistoMotherPi0PtOpenAngle = new TH2F*[fnCuts]; + fHistoMotherEtaPtOpenAngle = new TH2F*[fnCuts]; + } + } + fTagOutputList = new TList*[fnCuts]; + + fHistoConvGammaUntagged = new TH1F*[fnCuts]; + fHistoConvGammaTagged = new TH1F*[fnCuts]; + fHistoConvGammaPi0Tagged = new TH1F*[fnCuts]; + fHistoConvGammaEtaTagged = new TH1F*[fnCuts]; + fHistoPhotonPairAll = new TH2F*[fnCuts]; + fHistoPhotonPairAllGam = new TH2F*[fnCuts]; + + fHistoClusGammaPt = new TH1F*[fnCuts]; + + for(Int_t iCut = 0; iCutAt(iCut))->GetCutNumber(); + TString cutstringCalo = ((AliCaloPhotonCuts*)fClusterCutArray->At(iCut))->GetCutNumber(); + TString cutstringMeson = "NoMesonCut"; + if(fDoMesonAnalysis)cutstringMeson = ((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->GetCutNumber(); + + fCutFolder[iCut] = new TList(); + fCutFolder[iCut]->SetName(Form("Cut Number %s_%s_%s",cutstring.Data(),cutstringCalo.Data(),cutstringMeson.Data())); + fCutFolder[iCut]->SetOwner(kTRUE); + fOutputContainer->Add(fCutFolder[iCut]); + fESDList[iCut] = new TList(); + fESDList[iCut]->SetName(Form("%s_%s_%s ESD histograms",cutstring.Data(),cutstringCalo.Data(),cutstringMeson.Data())); + fESDList[iCut]->SetOwner(kTRUE); + fCutFolder[iCut]->Add(fESDList[iCut]); + + fHistoNEvents[iCut] = new TH1I("NEvents","NEvents",9,-0.5,8.5); + fHistoNEvents[iCut]->GetXaxis()->SetBinLabel(1,"Accepted"); + fHistoNEvents[iCut]->GetXaxis()->SetBinLabel(2,"Centrality"); + fHistoNEvents[iCut]->GetXaxis()->SetBinLabel(3,"Missing MC"); + if (((AliConversionCuts*)fCutArray->At(iCut))->IsSpecialTrigger() == 4 ){ + TString TriggerNames = "Not Trigger: "; + TriggerNames = TriggerNames+ ( (AliConversionCuts*)fCutArray->At(iCut))->GetSpecialTriggerName(); + fHistoNEvents[iCut]->GetXaxis()->SetBinLabel(4,TriggerNames.Data()); + } else { + fHistoNEvents[iCut]->GetXaxis()->SetBinLabel(4,"Trigger"); + } + fHistoNEvents[iCut]->GetXaxis()->SetBinLabel(5,"Vertex Z"); + fHistoNEvents[iCut]->GetXaxis()->SetBinLabel(6,"Cont. Vertex"); + fHistoNEvents[iCut]->GetXaxis()->SetBinLabel(7,"Pile-Up"); + fHistoNEvents[iCut]->GetXaxis()->SetBinLabel(8,"no SDD"); + fHistoNEvents[iCut]->GetXaxis()->SetBinLabel(9,"no V0AND"); + fESDList[iCut]->Add(fHistoNEvents[iCut]); + + if(fIsHeavyIon == 1) fHistoNGoodESDTracks[iCut] = new TH1I("GoodESDTracks","GoodESDTracks",4000,0,4000); + else if(fIsHeavyIon == 2) fHistoNGoodESDTracks[iCut] = new TH1I("GoodESDTracks","GoodESDTracks",400,0,400); + else fHistoNGoodESDTracks[iCut] = new TH1I("GoodESDTracks","GoodESDTracks",200,0,200); + fESDList[iCut]->Add(fHistoNGoodESDTracks[iCut]); + if(fIsHeavyIon == 1) fHistoNGammaCandidates[iCut] = new TH1I("GammaCandidates","GammaCandidates",100,0,100); + else if(fIsHeavyIon == 2) fHistoNGammaCandidates[iCut] = new TH1I("GammaCandidates","GammaCandidates",50,0,50); + else fHistoNGammaCandidates[iCut] = new TH1I("GammaCandidates","GammaCandidates",50,0,50); + fESDList[iCut]->Add(fHistoNGammaCandidates[iCut]); + if(fIsHeavyIon == 1) fHistoNGoodESDTracksVsNGammaCanditates[iCut] = new TH2F("GoodESDTracksVsGammaCandidates","GoodESDTracksVsGammaCandidates",4000,0,4000,100,0,100); + else if(fIsHeavyIon == 2) fHistoNGoodESDTracksVsNGammaCanditates[iCut] = new TH2F("GoodESDTracksVsGammaCandidates","GoodESDTracksVsGammaCandidates",400,0,400,50,0,50); + else fHistoNGoodESDTracksVsNGammaCanditates[iCut] = new TH2F("GoodESDTracksVsGammaCandidates","GoodESDTracksVsGammaCandidates",200,0,200,50,0,50); + fESDList[iCut]->Add(fHistoNGoodESDTracksVsNGammaCanditates[iCut]); + + + if(fIsHeavyIon == 1) fHistoNV0Tracks[iCut] = new TH1I("V0 Multiplicity","V0 Multiplicity",30000,0,30000); + else if(fIsHeavyIon == 2) fHistoNV0Tracks[iCut] = new TH1I("V0 Multiplicity","V0 Multiplicity",2500,0,2500); + else fHistoNV0Tracks[iCut] = new TH1I("V0 Multiplicity","V0 Multiplicity",1500,0,1500); + fESDList[iCut]->Add(fHistoNV0Tracks[iCut]); + fProfileEtaShift[iCut] = new TProfile("Eta Shift","Eta Shift",1, -0.5,0.5); + fESDList[iCut]->Add(fProfileEtaShift[iCut]); + fHistoConvGammaPt[iCut] = new TH1F("ESD_ConvGamma_Pt","ESD_ConvGamma_Pt",250,0,25); + fESDList[iCut]->Add(fHistoConvGammaPt[iCut]); + + if (fDoPhotonQA == 2){ + fPhotonDCAList[iCut] = new TList(); + fPhotonDCAList[iCut]->SetName(Form("%s_%s_%s Photon DCA tree",cutstring.Data(),cutstringCalo.Data(),cutstringMeson.Data())); + fPhotonDCAList[iCut]->SetOwner(kTRUE); + fCutFolder[iCut]->Add(fPhotonDCAList[iCut]); + + fTreeConvGammaPtDcazCat[iCut] = new TTree("ESD_ConvGamma_Pt_Dcaz_R_Eta","ESD_ConvGamma_Pt_Dcaz_R_Eta_Cat"); + fTreeConvGammaPtDcazCat[iCut]->Branch("Pt",&fPtGamma,"fPtGamma/F"); + fTreeConvGammaPtDcazCat[iCut]->Branch("DcaZPhoton",&fDCAzPhoton,"fDCAzPhoton/F"); + // fTreeConvGammaPtDcazCat[iCut]->Branch("R",&fRConvPhoton,"fRConvPhoton/F"); + // fTreeConvGammaPtDcazCat[iCut]->Branch("Eta",&fEtaPhoton,"fEtaPhoton/F"); + + fTreeConvGammaPtDcazCat[iCut]->Branch("cat",&fCharCatPhoton,"fCharCatPhoton/b"); + if(fIsMC){ + fTreeConvGammaPtDcazCat[iCut]->Branch("photonMCInfo",&fCharPhotonMCInfo,"fCharPhotonMCInfo/b"); + } + fPhotonDCAList[iCut]->Add(fTreeConvGammaPtDcazCat[iCut]); + } + + if (fDoPhotonQA > 0){ + fHistoConvGammaR[iCut] = new TH1F("ESD_ConvGamma_R","ESD_ConvGamma_R",800,0,200); + fESDList[iCut]->Add(fHistoConvGammaR[iCut]); + fHistoConvGammaEta[iCut] = new TH1F("ESD_ConvGamma_Eta","ESD_ConvGamma_Eta",2000,-2,2); + fESDList[iCut]->Add(fHistoConvGammaEta[iCut]); + } + + fTagOutputList[iCut] = new TList(); + fTagOutputList[iCut]->SetName(Form("Tagging Output Cut Number %s_%s_%s",cutstring.Data(),cutstringCalo.Data(),cutstringMeson.Data())); + fTagOutputList[iCut]->SetOwner(1); + fCutFolder[iCut]->Add(fTagOutputList[iCut]); + + const Int_t nptbins = 200; + const Double_t ptmin = 0.; + const Double_t ptmax = 20.; + + const Int_t nmbins = 180; + const Double_t mmin = 0.; + const Double_t mmax = 0.9; + + // photon candidates + // this is maybe not necessary ... + + fHistoConvGammaUntagged[iCut] = new TH1F("ConvGammaUntagged","",nptbins,ptmin,ptmax); + fHistoConvGammaUntagged[iCut]->SetXTitle("p_{T} (GeV/c)"); + fTagOutputList[iCut]->Add(fHistoConvGammaUntagged[iCut]); + + fHistoConvGammaTagged[iCut] = new TH1F("ConvGammaTagged","",nptbins,ptmin,ptmax); + fHistoConvGammaTagged[iCut]->SetXTitle("p_{T} (GeV/c)"); + fTagOutputList[iCut]->Add(fHistoConvGammaTagged[iCut]); + + fHistoConvGammaPi0Tagged[iCut] = new TH1F("ConvGammaPi0Tagged","",nptbins,ptmin,ptmax); + fHistoConvGammaPi0Tagged[iCut]->SetXTitle("p_{T} (GeV/c)"); + fTagOutputList[iCut]->Add(fHistoConvGammaPi0Tagged[iCut]); + + fHistoConvGammaEtaTagged[iCut] = new TH1F("ConvGammaEtaTagged","",nptbins,ptmin,ptmax); + fHistoConvGammaEtaTagged[iCut]->SetXTitle("p_{T} (GeV/c)"); + fTagOutputList[iCut]->Add(fHistoConvGammaEtaTagged[iCut]); + + // pairs + fHistoPhotonPairAll[iCut] = new TH2F("PhotonPairAll","",nmbins,mmin,mmax,nptbins,ptmin,ptmax); + fHistoPhotonPairAll[iCut]->SetXTitle("M_{inv} (GeV/cc)"); + fHistoPhotonPairAll[iCut]->SetYTitle("p_{T} (GeV/c)"); + fTagOutputList[iCut]->Add(fHistoPhotonPairAll[iCut]); + + fHistoPhotonPairAllGam[iCut] = new TH2F("PhotonPairAllGammaConvPt","",nmbins,mmin,mmax,nptbins,ptmin,ptmax); + fHistoPhotonPairAllGam[iCut]->SetXTitle("M_{inv} (GeV/cc)"); + fHistoPhotonPairAllGam[iCut]->SetYTitle("#gamma^{conv} p_{T} (GeV/c)"); + fTagOutputList[iCut]->Add(fHistoPhotonPairAllGam[iCut]); + + fHistoClusGammaPt[iCut] = new TH1F("ClusGamma_Pt","ClusGamma_Pt",250,0,25); + fTagOutputList[iCut]->Add(fHistoClusGammaPt[iCut]); + + + if(fDoMesonAnalysis){ + fHistoMotherInvMassPt[iCut] = new TH2F("ESD_Mother_InvMass_Pt","ESD_Mother_InvMass_Pt",800,0,0.8,250,0,25); + fESDList[iCut]->Add(fHistoMotherInvMassPt[iCut]); + fHistoMotherBackInvMassPt[iCut] = new TH2F("ESD_Background_InvMass_Pt","ESD_Background_InvMass_Pt",800,0,0.8,250,0,25); + fESDList[iCut]->Add(fHistoMotherBackInvMassPt[iCut]); + fHistoMotherInvMassEalpha[iCut] = new TH2F("ESD_Mother_InvMass_vs_E_alpha","ESD_Mother_InvMass_vs_E_alpha",800,0,0.8,250,0,25); + fESDList[iCut]->Add(fHistoMotherInvMassEalpha[iCut]); + if (fDoMesonQA == 2){ + fMesonDCAList[iCut] = new TList(); + fMesonDCAList[iCut]->SetName(Form("%s_%s_%s Meson DCA tree",cutstring.Data(),cutstringCalo.Data(),cutstringMeson.Data())); + fMesonDCAList[iCut]->SetOwner(kTRUE); + fCutFolder[iCut]->Add(fMesonDCAList[iCut]); + + fTreeMesonsInvMassPtDcazMinDcazMaxFlag[iCut] = new TTree("ESD_Mesons_InvMass_Pt_DcazMin_DcazMax_Flag","ESD_Mesons_InvMass_Pt_DcazMin_DcazMax_Flag"); + fTreeMesonsInvMassPtDcazMinDcazMaxFlag[iCut]->Branch("InvMass",&fInvMass,"fInvMass/F"); + fTreeMesonsInvMassPtDcazMinDcazMaxFlag[iCut]->Branch("Pt",&fPt,"fPt/F"); + fTreeMesonsInvMassPtDcazMinDcazMaxFlag[iCut]->Branch("DcaZMin",&fDCAzGammaMin,"fDCAzGammaMin/F"); + fTreeMesonsInvMassPtDcazMinDcazMaxFlag[iCut]->Branch("DcaZMax",&fDCAzGammaMax,"fDCAzGammaMax/F"); + fTreeMesonsInvMassPtDcazMinDcazMaxFlag[iCut]->Branch("kind",&fCharFlag,"fCharFlag/b"); + if(fIsMC){ + fTreeMesonsInvMassPtDcazMinDcazMaxFlag[iCut]->Branch("mesonMCInfo",&fCharMesonMCInfo,"fCharMesonMCInfo/b"); + } + fMesonDCAList[iCut]->Add(fTreeMesonsInvMassPtDcazMinDcazMaxFlag[iCut]); + + } + if (fDoMesonQA > 0 ){ + fHistoMotherPi0PtY[iCut] = new TH2F("ESD_MotherPi0_Pt_Y","ESD_MotherPi0_Pt_Y",150,0.03,15.,150,-1.5,1.5); + SetLogBinningXTH2(fHistoMotherPi0PtY[iCut]); + fESDList[iCut]->Add(fHistoMotherPi0PtY[iCut]); + fHistoMotherEtaPtY[iCut] = new TH2F("ESD_MotherEta_Pt_Y","ESD_MotherEta_Pt_Y",150,0.03,15.,150,-1.5,1.5); + SetLogBinningXTH2(fHistoMotherEtaPtY[iCut]); + fESDList[iCut]->Add(fHistoMotherEtaPtY[iCut]); + fHistoMotherPi0PtAlpha[iCut] = new TH2F("ESD_MotherPi0_Pt_Alpha","ESD_MotherPi0_Pt_Alpha",150,0.03,15.,100,0,1); + SetLogBinningXTH2(fHistoMotherPi0PtAlpha[iCut]); + fESDList[iCut]->Add(fHistoMotherPi0PtAlpha[iCut]); + fHistoMotherEtaPtAlpha[iCut] = new TH2F("ESD_MotherEta_Pt_Alpha","ESD_MotherEta_Pt_Alpha",150,0.03,15.,100,0,1); + SetLogBinningXTH2(fHistoMotherEtaPtAlpha[iCut]); + fESDList[iCut]->Add(fHistoMotherEtaPtAlpha[iCut]); + fHistoMotherPi0PtOpenAngle[iCut] = new TH2F("ESD_MotherPi0_Pt_OpenAngle","ESD_MotherPi0_Pt_OpenAngle",150,0.03,15.,200,0,2*TMath::Pi()); + SetLogBinningXTH2(fHistoMotherPi0PtOpenAngle[iCut]); + fESDList[iCut]->Add(fHistoMotherPi0PtOpenAngle[iCut]); + fHistoMotherEtaPtOpenAngle[iCut] = new TH2F("ESD_MotherEta_Pt_OpenAngle","ESD_MotherEta_Pt_OpenAngle",150,0.03,15.,200,0,2*TMath::Pi()); + SetLogBinningXTH2(fHistoMotherEtaPtOpenAngle[iCut]); + fESDList[iCut]->Add(fHistoMotherEtaPtOpenAngle[iCut]); + } + } + } + if(fDoMesonAnalysis){ + InitBack(); // Init Background Handler + } + + if(fIsMC){ + // MC Histogramms + fMCList = new TList*[fnCuts]; + // True Histogramms + fTrueList = new TList*[fnCuts]; + // Selected Header List + fHeaderNameList = new TList*[fnCuts]; + fHistoMCHeaders = new TH1I*[fnCuts]; + fHistoMCAllGammaPt = new TH1F*[fnCuts]; + fHistoMCDecayGammaPi0Pt = new TH1F*[fnCuts]; + fHistoMCDecayGammaRhoPt = new TH1F*[fnCuts]; + fHistoMCDecayGammaEtaPt = new TH1F*[fnCuts]; + fHistoMCDecayGammaOmegaPt = new TH1F*[fnCuts]; + fHistoMCDecayGammaEtapPt = new TH1F*[fnCuts]; + fHistoMCDecayGammaPhiPt = new TH1F*[fnCuts]; + fHistoMCDecayGammaSigmaPt = new TH1F*[fnCuts]; + fHistoMCConvGammaPt = new TH1F*[fnCuts]; + fHistoTrueConvGammaPt = new TH1F*[fnCuts]; + fHistoTrueConvPi0GammaPt = new TH1F*[fnCuts]; + + fHistoCombinatorialPt = new TH2F*[fnCuts]; + fHistoTruePrimaryConvGammaPt = new TH1F*[fnCuts]; + fHistoTruePrimaryConvGammaESDPtMCPt = new TH2F*[fnCuts]; + fHistoTrueSecondaryConvGammaPt = new TH1F*[fnCuts]; + fHistoTrueSecondaryConvGammaFromXFromK0sPt = new TH1F*[fnCuts]; + fHistoTrueSecondaryConvGammaFromXFromLambdaPt = new TH1F*[fnCuts]; + + fHistoTrueDalitzPsiPairDeltaPhi= new TH2F*[fnCuts]; + fHistoTrueGammaPsiPairDeltaPhi= new TH2F*[fnCuts]; + + fHistoTrueClusGammaPt = new TH1F*[fnCuts]; + fHistoTruePrimaryClusGammaPt = new TH1F*[fnCuts]; + fHistoTruePrimaryClusGammaESDPtMCPt = new TH2F*[fnCuts]; + + if (fDoPhotonQA > 0){ + fHistoMCConvGammaR = new TH1F*[fnCuts]; + fHistoMCConvGammaEta = new TH1F*[fnCuts]; + fHistoTrueConvGammaEta = new TH1F*[fnCuts]; + } + + if(fDoMesonAnalysis){ + fHistoMCPi0Pt = new TH1F*[fnCuts]; + fHistoMCPi0WOWeightPt = new TH1F*[fnCuts]; + fHistoMCEtaPt = new TH1F*[fnCuts]; + fHistoMCEtaWOWeightPt = new TH1F*[fnCuts]; + fHistoMCPi0InAccPt = new TH1F*[fnCuts]; + fHistoMCEtaInAccPt = new TH1F*[fnCuts]; + + fHistoTrueMotherInvMassPt = new TH2F*[fnCuts]; + fHistoTruePrimaryMotherInvMassPt = new TH2F*[fnCuts]; + fHistoTruePrimaryMotherW0WeightingInvMassPt = new TH2F*[fnCuts]; + fProfileTruePrimaryMotherWeightsInvMassPt = new TProfile2D*[fnCuts]; + fHistoTrueSecondaryMotherInvMassPt = new TH2F*[fnCuts]; + fHistoTrueSecondaryMotherFromK0sInvMassPt = new TH2F*[fnCuts]; + fHistoTrueSecondaryMotherFromEtaInvMassPt = new TH2F*[fnCuts]; + fHistoTrueSecondaryMotherFromLambdaInvMassPt = new TH2F*[fnCuts]; + fHistoTrueMotherDalitzInvMassPt = new TH2F*[fnCuts]; + if (fDoMesonQA > 0){ + fHistoMCPi0PtY = new TH2F*[fnCuts]; + fHistoMCEtaPtY = new TH2F*[fnCuts]; + fHistoMCK0sPt = new TH1F*[fnCuts]; + fHistoMCK0sWOWeightPt = new TH1F*[fnCuts]; + fHistoMCK0sPtY = new TH2F*[fnCuts]; + fHistoMCSecPi0PtvsSource= new TH2F*[fnCuts]; + fHistoMCSecPi0Source = new TH1F*[fnCuts]; + fHistoMCSecEtaPt = new TH1F*[fnCuts]; + fHistoMCSecEtaSource = new TH1F*[fnCuts]; + fHistoTruePrimaryPi0MCPtResolPt = new TH2F*[fnCuts]; + fHistoTruePrimaryEtaMCPtResolPt = new TH2F*[fnCuts]; + fHistoTrueK0sWithPi0DaughterMCPt = new TH1F*[fnCuts]; + fHistoTrueEtaWithPi0DaughterMCPt = new TH1F*[fnCuts]; + fHistoTrueLambdaWithPi0DaughterMCPt = new TH1F*[fnCuts]; + fHistoTrueBckGGInvMassPt = new TH2F*[fnCuts]; + fHistoTrueBckContInvMassPt = new TH2F*[fnCuts]; + fHistoTruePi0PtY = new TH2F*[fnCuts]; + fHistoTrueEtaPtY = new TH2F*[fnCuts]; + fHistoTruePi0PtAlpha = new TH2F*[fnCuts]; + fHistoTrueEtaPtAlpha = new TH2F*[fnCuts]; + fHistoTruePi0PtOpenAngle = new TH2F*[fnCuts]; + fHistoTrueEtaPtOpenAngle = new TH2F*[fnCuts]; + } + } + + + + for(Int_t iCut = 0; iCutAt(iCut))->GetCutNumber(); + TString cutstringCalo = ((AliCaloPhotonCuts*)fClusterCutArray->At(iCut))->GetCutNumber(); + TString cutstringMeson = "NoMesonCut"; + if(fDoMesonAnalysis)cutstringMeson = ((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->GetCutNumber(); + + fMCList[iCut] = new TList(); + fMCList[iCut]->SetName(Form("%s_%s_%s MC histograms",cutstring.Data(),cutstringCalo.Data(),cutstringMeson.Data())); + fMCList[iCut]->SetOwner(kTRUE); + fCutFolder[iCut]->Add(fMCList[iCut]); + fHistoMCHeaders[iCut] = new TH1I("MC_Headers","MC_Headers",20,0,20); + fMCList[iCut]->Add(fHistoMCHeaders[iCut]); + fHistoMCAllGammaPt[iCut] = new TH1F("MC_AllGamma_Pt","MC_AllGamma_Pt",250,0,25); + fMCList[iCut]->Add(fHistoMCAllGammaPt[iCut]); + fHistoMCDecayGammaPi0Pt[iCut] = new TH1F("MC_DecayGammaPi0_Pt","MC_DecayGammaPi0_Pt",250,0,25); + fMCList[iCut]->Add(fHistoMCDecayGammaPi0Pt[iCut]); + fHistoMCDecayGammaRhoPt[iCut] = new TH1F("MC_DecayGammaRho_Pt","MC_DecayGammaRho_Pt",250,0,25); + fMCList[iCut]->Add(fHistoMCDecayGammaRhoPt[iCut]); + fHistoMCDecayGammaEtaPt[iCut] = new TH1F("MC_DecayGammaEta_Pt","MC_DecayGammaEta_Pt",250,0,25); + fMCList[iCut]->Add(fHistoMCDecayGammaEtaPt[iCut]); + fHistoMCDecayGammaOmegaPt[iCut] = new TH1F("MC_DecayGammaOmega_Pt","MC_DecayGammaOmmega_Pt",250,0,25); + fMCList[iCut]->Add(fHistoMCDecayGammaOmegaPt[iCut]); + fHistoMCDecayGammaEtapPt[iCut] = new TH1F("MC_DecayGammaEtap_Pt","MC_DecayGammaEtap_Pt",250,0,25); + fMCList[iCut]->Add(fHistoMCDecayGammaEtapPt[iCut]); + fHistoMCDecayGammaPhiPt[iCut] = new TH1F("MC_DecayGammaPhi_Pt","MC_DecayGammaPhi_Pt",250,0,25); + fMCList[iCut]->Add(fHistoMCDecayGammaPhiPt[iCut]); + fHistoMCDecayGammaSigmaPt[iCut] = new TH1F("MC_DecayGammaSigma_Pt","MC_DecayGammaSigma_Pt",250,0,25); + fMCList[iCut]->Add(fHistoMCDecayGammaSigmaPt[iCut]); + fHistoMCConvGammaPt[iCut] = new TH1F("MC_ConvGamma_Pt","MC_ConvGamma_Pt",250,0,25); + fMCList[iCut]->Add(fHistoMCConvGammaPt[iCut]); + + if (fDoPhotonQA > 0){ + fHistoMCConvGammaR[iCut] = new TH1F("MC_ConvGamma_R","MC_ConvGamma_R",800,0,200); + fMCList[iCut]->Add(fHistoMCConvGammaR[iCut]); + fHistoMCConvGammaEta[iCut] = new TH1F("MC_ConvGamma_Eta","MC_ConvGamma_Eta",2000,-2,2); + fMCList[iCut]->Add(fHistoMCConvGammaEta[iCut]); + } + + if(fDoMesonAnalysis){ + fHistoMCPi0Pt[iCut] = new TH1F("MC_Pi0_Pt","MC_Pi0_Pt",250,0,25); + fHistoMCPi0Pt[iCut]->Sumw2(); + fMCList[iCut]->Add(fHistoMCPi0Pt[iCut]); + fHistoMCPi0WOWeightPt[iCut] = new TH1F("MC_Pi0_WOWeights_Pt","MC_Pi0_WOWeights_Pt",250,0,25); + fHistoMCPi0WOWeightPt[iCut]->Sumw2(); + fMCList[iCut]->Add(fHistoMCPi0WOWeightPt[iCut]); + + fHistoMCEtaPt[iCut] = new TH1F("MC_Eta_Pt","MC_Eta_Pt",250,0,25); + fHistoMCEtaPt[iCut]->Sumw2(); + fMCList[iCut]->Add(fHistoMCEtaPt[iCut]); + fHistoMCEtaWOWeightPt[iCut] = new TH1F("MC_Eta_WOWeights_Pt","MC_Eta_WOWeights_Pt",250,0,25); + fHistoMCEtaWOWeightPt[iCut]->Sumw2(); + fMCList[iCut]->Add(fHistoMCEtaWOWeightPt[iCut]); + fHistoMCPi0InAccPt[iCut] = new TH1F("MC_Pi0InAcc_Pt","MC_Pi0InAcc_Pt",250,0,25); + fHistoMCPi0InAccPt[iCut]->Sumw2(); + fMCList[iCut]->Add(fHistoMCPi0InAccPt[iCut]); + fHistoMCEtaInAccPt[iCut] = new TH1F("MC_EtaInAcc_Pt","MC_EtaInAcc_Pt",250,0,25); + fHistoMCEtaInAccPt[iCut]->Sumw2(); + fMCList[iCut]->Add(fHistoMCEtaInAccPt[iCut]); + if (fDoMesonQA > 0){ + fHistoMCPi0PtY[iCut] = new TH2F("MC_Pi0_Pt_Y","MC_Pi0_Pt_Y",150,0.03,15.,150,-1.5,1.5); + fHistoMCPi0PtY[iCut]->Sumw2(); + SetLogBinningXTH2(fHistoMCPi0PtY[iCut]); + fMCList[iCut]->Add(fHistoMCPi0PtY[iCut]); + fHistoMCEtaPtY[iCut] = new TH2F("MC_Eta_Pt_Y","MC_Eta_Pt_Y",150,0.03,15.,150,-1.5,1.5); + fHistoMCEtaPtY[iCut]->Sumw2(); + SetLogBinningXTH2(fHistoMCEtaPtY[iCut]); + fMCList[iCut]->Add(fHistoMCEtaPtY[iCut]); + fHistoMCK0sPt[iCut] = new TH1F("MC_K0s_Pt","MC_K0s_Pt",150,0,15); + fHistoMCK0sPt[iCut]->Sumw2(); + fMCList[iCut]->Add(fHistoMCK0sPt[iCut]); + fHistoMCK0sWOWeightPt[iCut] = new TH1F("MC_K0s_WOWeights_Pt","MC_K0s_WOWeights_Pt",150,0,15); + fHistoMCK0sWOWeightPt[iCut]->Sumw2(); + fMCList[iCut]->Add(fHistoMCK0sWOWeightPt[iCut]); + fHistoMCK0sPtY[iCut] = new TH2F("MC_K0s_Pt_Y","MC_K0s_Pt_Y",150,0.03,15.,150,-1.5,1.5); + fHistoMCK0sPtY[iCut]->Sumw2(); + SetLogBinningXTH2(fHistoMCK0sPtY[iCut]); + fMCList[iCut]->Add(fHistoMCK0sPtY[iCut]); + + fHistoMCSecPi0Source[iCut] = new TH1F("MC_SecPi0_Source","MC_SecPi0_Source",5000,0.,5000); + fMCList[iCut]->Add(fHistoMCSecPi0Source[iCut]); + fHistoMCSecEtaSource[iCut] = new TH1F("MC_SecEta_Source","MC_SecEta_Source",5000,0,5000); + fMCList[iCut]->Add(fHistoMCSecEtaSource[iCut]); + fHistoMCSecPi0PtvsSource[iCut] = new TH2F("MC_SecPi0_Pt_Source","MC_SecPi0_Pt_Source",250,0.0,25.,16,-0.5,15.5); + fHistoMCSecPi0PtvsSource[iCut]->Sumw2(); + fMCList[iCut]->Add(fHistoMCSecPi0PtvsSource[iCut]); + fHistoMCSecEtaPt[iCut] = new TH1F("MC_SecEta_Pt","MC_SecEta_Pt",250,0,25); + fHistoMCSecEtaPt[iCut]->Sumw2(); + fMCList[iCut]->Add(fHistoMCSecEtaPt[iCut]); + } + + } + fTrueList[iCut] = new TList(); + fTrueList[iCut]->SetName(Form("%s_%s_%s True histograms",cutstring.Data(),cutstringCalo.Data(),cutstringMeson.Data())); + fTrueList[iCut]->SetOwner(kTRUE); + fCutFolder[iCut]->Add(fTrueList[iCut]); + + fHistoTrueConvGammaPt[iCut] = new TH1F("ESD_TrueConvGamma_Pt","ESD_TrueConvGamma_Pt",250,0,25); + fTrueList[iCut]->Add(fHistoTrueConvGammaPt[iCut]); + + fHistoTrueConvPi0GammaPt[iCut] = new TH1F("ESD_TrueConvGamma_Pt","ESD_TrueConvGamma_Pt",250,0,25); + fTrueList[iCut]->Add(fHistoTrueConvPi0GammaPt[iCut]); + + fHistoCombinatorialPt[iCut] = new TH2F("ESD_TrueCombinatorial_Pt","ESD_TrueCombinatorial_Pt",250,0,25,16,-0.5,15.5); + fHistoCombinatorialPt[iCut]->GetYaxis()->SetBinLabel( 1,"Elec+Elec"); + fHistoCombinatorialPt[iCut]->GetYaxis()->SetBinLabel( 2,"Elec+Pion"); + fHistoCombinatorialPt[iCut]->GetYaxis()->SetBinLabel( 3,"Elec+Kaon"); + fHistoCombinatorialPt[iCut]->GetYaxis()->SetBinLabel( 4,"Elec+Proton"); + fHistoCombinatorialPt[iCut]->GetYaxis()->SetBinLabel( 5,"Elec+Muon"); + fHistoCombinatorialPt[iCut]->GetYaxis()->SetBinLabel( 6,"Pion+Pion"); + fHistoCombinatorialPt[iCut]->GetYaxis()->SetBinLabel( 7,"Pion+Kaon"); + fHistoCombinatorialPt[iCut]->GetYaxis()->SetBinLabel( 8,"Pion+Proton"); + fHistoCombinatorialPt[iCut]->GetYaxis()->SetBinLabel( 9,"Pion+Muon"); + fHistoCombinatorialPt[iCut]->GetYaxis()->SetBinLabel(10,"Kaon+Kaon"); + fHistoCombinatorialPt[iCut]->GetYaxis()->SetBinLabel(11,"Kaon+Proton"); + fHistoCombinatorialPt[iCut]->GetYaxis()->SetBinLabel(12,"Kaon+Muon"); + fHistoCombinatorialPt[iCut]->GetYaxis()->SetBinLabel(13,"Proton+Proton"); + fHistoCombinatorialPt[iCut]->GetYaxis()->SetBinLabel(14,"Proton+Muon"); + fHistoCombinatorialPt[iCut]->GetYaxis()->SetBinLabel(15,"Muon+Muon"); + fHistoCombinatorialPt[iCut]->GetYaxis()->SetBinLabel(16,"Rest"); + fTrueList[iCut]->Add(fHistoCombinatorialPt[iCut]); + fHistoTruePrimaryConvGammaPt[iCut] = new TH1F("ESD_TruePrimaryConvGamma_Pt","ESD_TruePrimaryConvGamma_Pt",250,0,25); + fTrueList[iCut]->Add(fHistoTruePrimaryConvGammaPt[iCut]); + fHistoTrueSecondaryConvGammaPt[iCut] = new TH1F("ESD_TrueSecondaryConvGamma_Pt","ESD_TrueSecondaryConvGamma_Pt",250,0,25); + fTrueList[iCut]->Add(fHistoTrueSecondaryConvGammaPt[iCut]); + + fHistoTrueSecondaryConvGammaFromXFromK0sPt[iCut] = new TH1F("ESD_TrueSecondaryConvGammaFromXFromK0s_Pt", "ESD_TrueSecondaryConvGammaFromXFromK0s_Pt",250,0,25); + fTrueList[iCut]->Add(fHistoTrueSecondaryConvGammaFromXFromK0sPt[iCut]); + fHistoTrueSecondaryConvGammaFromXFromLambdaPt[iCut] = new TH1F("ESD_TrueSecondaryConvGammaFromXFromLambda_Pt", "ESD_TrueSecondaryConvGammaFromXFromLambda_Pt",250,0,25); + fTrueList[iCut]->Add(fHistoTrueSecondaryConvGammaFromXFromLambdaPt[iCut]); + + fHistoTrueDalitzPsiPairDeltaPhi[iCut] = new TH2F("ESD_TrueDalitzPsiPairDeltaPhi_Pt", "ESD_TrueDalitzPsiPairDeltaPhi_Pt",100,-0.5,2,100,-0.5,0.5); + fTrueList[iCut]->Add(fHistoTrueDalitzPsiPairDeltaPhi[iCut]); + + fHistoTrueGammaPsiPairDeltaPhi[iCut] = new TH2F("ESD_TrueGammaPsiPairDeltaPhi_Pt", "ESD_TrueGammaPsiPairDeltaPhi_Pt",100,-0.5,2,100,-0.5,0.5); + fTrueList[iCut]->Add(fHistoTrueGammaPsiPairDeltaPhi[iCut]); + + fHistoTruePrimaryConvGammaESDPtMCPt[iCut] = new TH2F("ESD_TruePrimaryConvGammaESD_PtMCPt", "ESD_TruePrimaryConvGammaESD_PtMCPt",250,0,25,250,0,25); + fTrueList[iCut]->Add(fHistoTruePrimaryConvGammaESDPtMCPt[iCut]); + + fHistoTrueClusGammaPt[iCut] = new TH1F("TrueClusGamma_Pt","ESD_TruePrimaryClusGamma_Pt",250,0,25); + fTagOutputList[iCut]->Add(fHistoTrueClusGammaPt[iCut]); + fHistoTruePrimaryClusGammaPt[iCut] = new TH1F("TruePrimaryClusGamma_Pt","ESD_TruePrimaryClusGamma_Pt",250,0,25); + fTagOutputList[iCut]->Add(fHistoTruePrimaryClusGammaPt[iCut]); + fHistoTruePrimaryClusGammaESDPtMCPt[iCut] = new TH2F("TruePrimaryClusGamma_Pt_MCPt","ESD_TruePrimaryClusGamma_MCPt",250,0,25,250,0,25); + fTagOutputList[iCut]->Add(fHistoTruePrimaryClusGammaESDPtMCPt[iCut]); + + if(fDoMesonAnalysis){ + fHistoTrueMotherInvMassPt[iCut] = new TH2F("ESD_TrueMother_InvMass_Pt","ESD_TrueMother_InvMass_Pt",800,0,0.8,250,0,25); + fTrueList[iCut]->Add(fHistoTrueMotherInvMassPt[iCut]); + fHistoTruePrimaryMotherInvMassPt[iCut] = new TH2F("ESD_TruePrimaryMother_InvMass_Pt", "ESD_TruePrimaryMother_InvMass_Pt", 800,0,0.8,250,0,25); + fHistoTruePrimaryMotherInvMassPt[iCut]->Sumw2(); + fTrueList[iCut]->Add(fHistoTruePrimaryMotherInvMassPt[iCut]); + fHistoTruePrimaryMotherW0WeightingInvMassPt[iCut] = new TH2F("ESD_TruePrimaryMotherW0Weights_InvMass_Pt", "ESD_TruePrimaryMotherW0Weights_InvMass_Pt", 800,0,0.8,250,0,25); + fHistoTruePrimaryMotherW0WeightingInvMassPt[iCut]->Sumw2(); + fTrueList[iCut]->Add(fHistoTruePrimaryMotherW0WeightingInvMassPt[iCut]); + fProfileTruePrimaryMotherWeightsInvMassPt[iCut] = new TProfile2D("ESD_TruePrimaryMotherWeights_InvMass_Pt", "ESD_TruePrimaryMotherWeights_InvMass_Pt", 800,0,0.8,250,0,25); + fProfileTruePrimaryMotherWeightsInvMassPt[iCut]->Sumw2(); + fTrueList[iCut]->Add(fProfileTruePrimaryMotherWeightsInvMassPt[iCut]); + fHistoTrueSecondaryMotherInvMassPt[iCut] = new TH2F("ESD_TrueSecondaryMother_InvMass_Pt", "ESD_TrueSecondaryMother_InvMass_Pt", 800,0,0.8,250,0,25); + fHistoTrueSecondaryMotherInvMassPt[iCut]->Sumw2(); + fTrueList[iCut]->Add(fHistoTrueSecondaryMotherInvMassPt[iCut]); + fHistoTrueSecondaryMotherFromK0sInvMassPt[iCut] = new TH2F("ESD_TrueSecondaryMotherFromK0s_InvMass_Pt","ESD_TrueSecondaryMotherFromK0s_InvMass_Pt",800,0,0.8,250,0,25); + fHistoTrueSecondaryMotherFromK0sInvMassPt[iCut]->Sumw2(); + fTrueList[iCut]->Add(fHistoTrueSecondaryMotherFromK0sInvMassPt[iCut]); + fHistoTrueSecondaryMotherFromEtaInvMassPt[iCut] = new TH2F("ESD_TrueSecondaryMotherFromEta_InvMass_Pt","ESD_TrueSecondaryMotherFromEta_InvMass_Pt",800,0,0.8,250,0,25); + fTrueList[iCut]->Add(fHistoTrueSecondaryMotherFromEtaInvMassPt[iCut]); + fHistoTrueSecondaryMotherFromLambdaInvMassPt[iCut] = new TH2F("ESD_TrueSecondaryMotherFromLambda_InvMass_Pt","ESD_TrueSecondaryMotherFromLambda_InvMass_Pt",800,0,0.8,250,0,25); + fTrueList[iCut]->Add(fHistoTrueSecondaryMotherFromLambdaInvMassPt[iCut]); + fHistoTrueMotherDalitzInvMassPt[iCut] = new TH2F("ESD_TrueDalitz_InvMass_Pt","ESD_TrueDalitz_InvMass_Pt",800,0,0.8,250,0,25); + fTrueList[iCut]->Add(fHistoTrueMotherDalitzInvMassPt[iCut]); + if (fDoMesonQA > 0){ + fHistoTruePrimaryPi0MCPtResolPt[iCut] = new TH2F("ESD_TruePrimaryPi0_MCPt_ResolPt","ESD_TruePrimaryPi0_ResolPt_MCPt",500,0.03,25,1000,-1.,1.); + fHistoTruePrimaryPi0MCPtResolPt[iCut]->Sumw2(); + SetLogBinningXTH2(fHistoTruePrimaryPi0MCPtResolPt[iCut]); + fTrueList[iCut]->Add(fHistoTruePrimaryPi0MCPtResolPt[iCut]); + fHistoTruePrimaryEtaMCPtResolPt[iCut] = new TH2F("ESD_TruePrimaryEta_MCPt_ResolPt","ESD_TruePrimaryEta_ResolPt_MCPt",500,0.03,25,1000,-1.,1.); + fHistoTruePrimaryEtaMCPtResolPt[iCut]->Sumw2(); + SetLogBinningXTH2(fHistoTruePrimaryEtaMCPtResolPt[iCut]); + fTrueList[iCut]->Add(fHistoTruePrimaryEtaMCPtResolPt[iCut]); + fHistoTrueBckGGInvMassPt[iCut] = new TH2F("ESD_TrueBckGG_InvMass_Pt","ESD_TrueBckGG_InvMass_Pt",800,0,0.8,250,0,25); + fTrueList[iCut]->Add(fHistoTrueBckGGInvMassPt[iCut]); + fHistoTrueBckContInvMassPt[iCut] = new TH2F("ESD_TrueBckCont_InvMass_Pt","ESD_TrueBckCont_InvMass_Pt",800,0,0.8,250,0,25); + fTrueList[iCut]->Add(fHistoTrueBckContInvMassPt[iCut]); + fHistoTrueK0sWithPi0DaughterMCPt[iCut] = new TH1F("ESD_TrueK0sWithPi0Daughter_MCPt","ESD_TrueK0sWithPi0Daughter_MCPt",250,0,25); + fTrueList[iCut]->Add(fHistoTrueK0sWithPi0DaughterMCPt[iCut]); + fHistoTrueEtaWithPi0DaughterMCPt[iCut] = new TH1F("ESD_TrueEtaWithPi0Daughter_MCPt","ESD_TrueEtaWithPi0Daughter_MCPt",250,0,25); + fTrueList[iCut]->Add(fHistoTrueEtaWithPi0DaughterMCPt[iCut]); + fHistoTrueLambdaWithPi0DaughterMCPt[iCut] = new TH1F("ESD_TrueLambdaWithPi0Daughter_MCPt","ESD_TrueLambdaWithPi0Daughter_MCPt",250,0,25); + fTrueList[iCut]->Add(fHistoTrueLambdaWithPi0DaughterMCPt[iCut]); + + fHistoTruePi0PtY[iCut] = new TH2F("ESD_TruePi0_Pt_Y","ESD_TruePi0_Pt_Y",150,0.03,15.,150,-1.5,1.5); + SetLogBinningXTH2(fHistoTruePi0PtY[iCut]); + fTrueList[iCut]->Add(fHistoTruePi0PtY[iCut]); + fHistoTrueEtaPtY[iCut] = new TH2F("ESD_TrueEta_Pt_Y","ESD_TrueEta_Pt_Y",150,0.03,15.,150,-1.5,1.5); + SetLogBinningXTH2(fHistoTrueEtaPtY[iCut]); + fTrueList[iCut]->Add(fHistoTrueEtaPtY[iCut]); + fHistoTruePi0PtAlpha[iCut] = new TH2F("ESD_TruePi0_Pt_Alpha","ESD_TruePi0_Pt_Alpha",150,0.03,15.,100,0,1); + SetLogBinningXTH2(fHistoTruePi0PtAlpha[iCut]); + fTrueList[iCut]->Add(fHistoTruePi0PtAlpha[iCut]); + fHistoTrueEtaPtAlpha[iCut] = new TH2F("ESD_TrueEta_Pt_Alpha","ESD_TrueEta_Pt_Alpha",150,0.03,15.,100,0,1); + SetLogBinningXTH2(fHistoTrueEtaPtAlpha[iCut]); + fTrueList[iCut]->Add(fHistoTrueEtaPtAlpha[iCut]); + + fHistoTruePi0PtOpenAngle[iCut] = new TH2F("ESD_TruePi0_Pt_OpenAngle","ESD_TruePi0_Pt_OpenAngle",150,0.03,15.,200,0,2*TMath::Pi()); + SetLogBinningXTH2(fHistoTruePi0PtOpenAngle[iCut]); + fTrueList[iCut]->Add(fHistoTruePi0PtOpenAngle[iCut]); + fHistoTrueEtaPtOpenAngle[iCut] = new TH2F("ESD_TrueEta_Pt_OpenAngle","ESD_TrueEta_Pt_OpenAngle",150,0.03,15.,200,0,2*TMath::Pi()); + SetLogBinningXTH2(fHistoTrueEtaPtOpenAngle[iCut]); + fTrueList[iCut]->Add(fHistoTrueEtaPtOpenAngle[iCut]); + + fHistoTrueConvGammaEta[iCut] = new TH1F("ESD_TrueConvGamma_Eta","ESD_TrueConvGamma_Eta",2000,-2,2); + fTrueList[iCut]->Add(fHistoTrueConvGammaEta[iCut]); + + } + } + } + } + + fV0Reader=(AliV0ReaderV1*)AliAnalysisManager::GetAnalysisManager()->GetTask("V0ReaderV1"); + if(!fV0Reader){printf("Error: No V0 Reader");return;} // GetV0Reader + + if(fV0Reader) + if((AliConversionCuts*)fV0Reader->GetConversionCuts()) + if(((AliConversionCuts*)fV0Reader->GetConversionCuts())->GetCutHistograms()) + fOutputContainer->Add(((AliConversionCuts*)fV0Reader->GetConversionCuts())->GetCutHistograms()); + + for(Int_t iCut = 0; iCutAt(iCut))) continue; + if(((AliConversionCuts*)fCutArray->At(iCut))->GetCutHistograms()){ + fCutFolder[iCut]->Add(((AliConversionCuts*)fCutArray->At(iCut))->GetCutHistograms()); + } + if(!((AliCaloPhotonCuts*)fClusterCutArray->At(iCut))) continue; + if(((AliCaloPhotonCuts*)fClusterCutArray->At(iCut))->GetCutHistograms()){ + fCutFolder[iCut]->Add(((AliCaloPhotonCuts*)fClusterCutArray->At(iCut))->GetCutHistograms()); + } + if(fDoMesonAnalysis){ + if(!((AliConversionMesonCuts*)fMesonCutArray->At(iCut))) continue; + if(((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->GetCutHistograms()){ + fCutFolder[iCut]->Add(((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->GetCutHistograms()); + } + } + } + PostData(1, fOutputContainer); +} +//_____________________________________________________________________________ +Bool_t AliAnalysisTaskGammaConvCalo::Notify() +{ + for(Int_t iCut = 0; iCutAt(iCut))->GetDoEtaShift()){ + fProfileEtaShift[iCut]->Fill(0.,(((AliConversionCuts*)fCutArray->At(iCut))->GetEtaShift())); + continue; // No Eta Shift requested, continue + } + if(((AliConversionCuts*)fCutArray->At(iCut))->GetEtaShift() == 0.0){ // Eta Shift requested but not set, get shift automatically + ((AliConversionCuts*)fCutArray->At(iCut))->GetCorrectEtaShiftFromPeriod(fV0Reader->GetPeriodName()); + fProfileEtaShift[iCut]->Fill(0.,(((AliConversionCuts*)fCutArray->At(iCut))->GetEtaShift())); + ((AliConversionCuts*)fCutArray->At(iCut))->DoEtaShift(kFALSE); // Eta Shift Set, make sure that it is called only once + continue; + } + else{ + printf(" Gamma Conversion Task %s :: Eta Shift Manually Set to %f \n\n", + (((AliConversionCuts*)fCutArray->At(iCut))->GetCutNumber()).Data(),((AliConversionCuts*)fCutArray->At(iCut))->GetEtaShift()); + fProfileEtaShift[iCut]->Fill(0.,(((AliConversionCuts*)fCutArray->At(iCut))->GetEtaShift())); + ((AliConversionCuts*)fCutArray->At(iCut))->DoEtaShift(kFALSE); // Eta Shift Set, make sure that it is called only once + } + } + + return kTRUE; +} +//_____________________________________________________________________________ +void AliAnalysisTaskGammaConvCalo::UserExec(Option_t *) +{ + // + // Called for each event + // + Int_t eventQuality = ((AliConversionCuts*)fV0Reader->GetConversionCuts())->GetEventQuality(); + if(eventQuality == 2 || eventQuality == 3){// Event Not Accepted due to MC event missing or wrong trigger for V0ReaderV1 + for(Int_t iCut = 0; iCutFill(eventQuality); + } + return; + } + + if(fIsMC) fMCEvent = MCEvent(); + if(fMCEvent == NULL) fIsMC = kFALSE; + + fInputEvent = InputEvent(); + + if(fIsMC && fInputEvent->IsA()==AliESDEvent::Class()){ + fMCStack = fMCEvent->Stack(); + if(fMCStack == NULL) fIsMC = kFALSE; + } + + fReaderGammas = fV0Reader->GetReconstructedGammas(); // Gammas from default Cut + + // ------------------- BeginEvent ---------------------------- + + AliEventplane *EventPlane = fInputEvent->GetEventplane(); + if(fIsHeavyIon ==1)fEventPlaneAngle = EventPlane->GetEventplane("V0",fInputEvent,2); + else fEventPlaneAngle=0.0; + + if(fIsMC && fInputEvent->IsA()==AliAODEvent::Class() && !(fV0Reader->AreAODsRelabeled())){ + RelabelAODPhotonCandidates(kTRUE); // In case of AODMC relabeling MC + fV0Reader->RelabelAODs(kTRUE); + } + + + for(Int_t iCut = 0; iCutAt(iCut))->IsEventAcceptedByConversionCut(fV0Reader->GetConversionCuts(),fInputEvent,fMCEvent,fIsHeavyIon); + + if(eventNotAccepted){ + // cout << "event rejected due to wrong trigger: " <Fill(eventNotAccepted); // Check Centrality, PileUp, SDD and V0AND --> Not Accepted => eventQuality = 1 + continue; + } + + if(eventQuality != 0){// Event Not Accepted + //cout << "event rejected due to: " <Fill(eventQuality); + continue; + } + + fHistoNEvents[iCut]->Fill(eventQuality); // Should be 0 here + fHistoNGoodESDTracks[iCut]->Fill(fV0Reader->GetNumberOfPrimaryTracks()); + if(((AliConversionCuts*)fCutArray->At(iCut))->IsHeavyIon() == 2) fHistoNV0Tracks[iCut]->Fill(fInputEvent->GetVZEROData()->GetMTotV0A()); + else fHistoNV0Tracks[iCut]->Fill(fInputEvent->GetVZEROData()->GetMTotV0A()+fInputEvent->GetVZEROData()->GetMTotV0C()); + + if(fIsMC){ + // Process MC Particle + if(((AliConversionCuts*)fCutArray->At(iCut))->GetSignalRejection() != 0){ + if(fInputEvent->IsA()==AliESDEvent::Class()){ + ((AliConversionCuts*)fCutArray->At(iCut))->GetNotRejectedParticles(((AliConversionCuts*)fCutArray->At(iCut))->GetSignalRejection(), + ((AliConversionCuts*)fCutArray->At(iCut))->GetAcceptedHeader(), + fMCEvent); + } + else if(fInputEvent->IsA()==AliAODEvent::Class()){ + ((AliConversionCuts*)fCutArray->At(iCut))->GetNotRejectedParticles(((AliConversionCuts*)fCutArray->At(iCut))->GetSignalRejection(), + ((AliConversionCuts*)fCutArray->At(iCut))->GetAcceptedHeader(), + fInputEvent); + } + + if(((AliConversionCuts*)fCutArray->At(iCut))->GetAcceptedHeader()){ + for(Int_t i = 0;i<(((AliConversionCuts*)fCutArray->At(iCut))->GetAcceptedHeader())->GetEntries();i++){ + TString nameBin= fHistoMCHeaders[iCut]->GetXaxis()->GetBinLabel(i+1); + if (nameBin.CompareTo("")== 0){ + TString nameHeader = ((TObjString*)((TList*)((AliConversionCuts*)fCutArray->At(iCut)) + ->GetAcceptedHeader())->At(i))->GetString(); + fHistoMCHeaders[iCut]->GetXaxis()->SetBinLabel(i+1,nameHeader.Data()); + } + } + } + } + } + if(fIsMC){ + if(fInputEvent->IsA()==AliESDEvent::Class()) + ProcessMCParticles(); + if(fInputEvent->IsA()==AliAODEvent::Class()) + ProcessAODMCParticles(); + } + + // it is in the loop to have the same conversion cut string (used also for MC stuff that should be same for V0 and Cluster) + ProcessClusters(); // process calo clusters + ProcessPhotonCandidates(); // Process this cuts gammas + + fHistoNGammaCandidates[iCut]->Fill(fGammaCandidates->GetEntries()); + fHistoNGoodESDTracksVsNGammaCanditates[iCut]->Fill(fV0Reader->GetNumberOfPrimaryTracks(),fGammaCandidates->GetEntries()); + if(fDoMesonAnalysis){ // Meson Analysis + if(((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->UseMCPSmearing() && fIsMC){ + fUnsmearedPx = new Double_t[fGammaCandidates->GetEntries()]; // Store unsmeared Momenta + fUnsmearedPy = new Double_t[fGammaCandidates->GetEntries()]; + fUnsmearedPz = new Double_t[fGammaCandidates->GetEntries()]; + fUnsmearedE = new Double_t[fGammaCandidates->GetEntries()]; + + for(Int_t gamma=0;gammaGetEntries();gamma++){ // Smear the AODPhotons in MC + fUnsmearedPx[gamma] = ((AliAODConversionPhoton*)fGammaCandidates->At(gamma))->Px(); + fUnsmearedPy[gamma] = ((AliAODConversionPhoton*)fGammaCandidates->At(gamma))->Py(); + fUnsmearedPz[gamma] = ((AliAODConversionPhoton*)fGammaCandidates->At(gamma))->Pz(); + fUnsmearedE[gamma] = ((AliAODConversionPhoton*)fGammaCandidates->At(gamma))->E(); + ((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->SmearParticle(dynamic_cast(fGammaCandidates->At(gamma))); + } + } + + PhotonTagging(); // tag PCM photons with calorimeter + + CalculatePi0Candidates(); // Combine Gammas from conversion and from calo + + if(((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->DoBGCalculation()){ + if(((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->BackgroundHandlerType() == 0){ + CalculateBackground(); // Combinatorial Background + UpdateEventByEventData(); // Store Event for mixed Events + } + else{ + CalculateBackgroundRP(); // Combinatorial Background + fBGHandlerRP[iCut]->AddEvent(fGammaCandidates,fInputEvent); // Store Event for mixed Events + fBGClusHandlerRP[iCut]->AddEvent(fClusterCandidates,fInputEvent); // Store Event for mixed Events + } + } + if(((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->UseMCPSmearing() && fIsMC){ + for(Int_t gamma=0;gammaGetEntries();gamma++){ // Smear the AODPhotons in MC + ((AliAODConversionPhoton*)fGammaCandidates->At(gamma))->SetPx(fUnsmearedPx[gamma]); // Reset Unsmeared Momenta + ((AliAODConversionPhoton*)fGammaCandidates->At(gamma))->SetPy(fUnsmearedPy[gamma]); + ((AliAODConversionPhoton*)fGammaCandidates->At(gamma))->SetPz(fUnsmearedPz[gamma]); + ((AliAODConversionPhoton*)fGammaCandidates->At(gamma))->SetE(fUnsmearedE[gamma]); + } + delete[] fUnsmearedPx; fUnsmearedPx = 0x0; + delete[] fUnsmearedPy; fUnsmearedPy = 0x0; + delete[] fUnsmearedPz; fUnsmearedPz = 0x0; + delete[] fUnsmearedE; fUnsmearedE = 0x0; + } + } + + fGammaCandidates->Clear(); // delete this cuts good gammas + fClusterCandidates->Clear(); // delete cluster candidates + } + + if(fIsMC && fInputEvent->IsA()==AliAODEvent::Class() && !(fV0Reader->AreAODsRelabeled())){ + RelabelAODPhotonCandidates(kFALSE); // Back to ESDMC Label + fV0Reader->RelabelAODs(kFALSE); + } + + PostData(1, fOutputContainer); +} + +//________________________________________________________________________ +void AliAnalysisTaskGammaConvCalo::ProcessClusters() +{ + + Int_t nclus = 0; + nclus = fInputEvent->GetNumberOfCaloClusters(); + + cout << nclus << endl; + + if(nclus == 0) return; + + // vertex + Double_t vertex[3] = {0}; + InputEvent()->GetPrimaryVertex()->GetXYZ(vertex); + + // Loop over EMCal clusters + for(Int_t i = 0; i < nclus; i++){ + + AliVCluster* clus = fInputEvent->GetCaloCluster(i); + if (!clus) continue; + if(!((AliCaloPhotonCuts*)fClusterCutArray->At(fiCut))->ClusterIsSelected(clus,fInputEvent,fIsMC)) continue; + // TLorentzvector with cluster + TLorentzVector clusterVector; + clus->GetMomentum(clusterVector,vertex); + + TLorentzVector* tmpvec = new TLorentzVector(); + tmpvec->SetPxPyPzE(clusterVector.Px(),clusterVector.Py(),clusterVector.Pz(),clusterVector.E()); + + // convert to AODConversionPhoton + AliAODConversionPhoton *PhotonCandidate=new AliAODConversionPhoton(tmpvec); + if(!PhotonCandidate) continue; + + // get MC label + Int_t mclab[2] = {0,0}; + if(fIsMC){ + mclab[0] = clus->GetLabel(); + mclab[1] = clus->GetLabel(); + cout << "mclabels: " << mclab[0] << " " << mclab[1] << endl; + } + + PhotonCandidate->SetMCLabel(mclab); + + // fIsFromMBHeader = kTRUE; + // if(fIsMC){ + // Int_t isPosFromMBHeader + // = ((AliConversionCuts*)fCutArray->At(fiCut))->IsParticleFromBGEvent(mclab, fMCStack, fInputEvent); + // if(isPosFromMBHeader == 0 && ((AliConversionCuts*)fCutArray->At(fiCut))->GetSignalRejection() != 3) continue; + // + // Int_t isNegFromMBHeader + // = ((AliConversionCuts*)fCutArray->At(fiCut))->IsParticleFromBGEvent(mclab, fMCStack, fInputEvent); + // if(isNegFromMBHeader == 0 && ((AliConversionCuts*)fCutArray->At(fiCut))->GetSignalRejection() != 3) continue; + // + // if( (isNegFromMBHeader+isPosFromMBHeader) != 4) fIsFromMBHeader = kFALSE; + // } + + fHistoClusGammaPt[fiCut]->Fill(PhotonCandidate->Pt()); + fClusterCandidates->Add(PhotonCandidate); // if no second loop is required add to events good gammas + + // if(fIsFromMBHeader){ + // fHistoConvGammaPt[fiCut]->Fill(PhotonCandidate->Pt()); + // if (fDoPhotonQA > 0){ + // fHistoConvGammaR[fiCut]->Fill(PhotonCandidate->GetConversionRadius()); + // fHistoConvGammaEta[fiCut]->Fill(PhotonCandidate->Eta()); + // } + // } + + if(fIsMC){ + ProcessTrueClusterCandidates(PhotonCandidate); + } + + // if (fIsFromMBHeader && fDoPhotonQA == 2){ + // if (fIsHeavyIon == 1 && PhotonCandidate->Pt() > 0.399 && PhotonCandidate->Pt() < 12.){ + // fPtGamma = PhotonCandidate->Pt(); + // fDCAzPhoton = PhotonCandidate->GetDCAzToPrimVtx(); + // fRConvPhoton = PhotonCandidate->GetConversionRadius(); + // fEtaPhoton = PhotonCandidate->GetPhotonEta(); + // fCharCatPhoton = PhotonCandidate->GetPhotonQuality(); + // fTreeConvGammaPtDcazCat[fiCut]->Fill(); + // } else if ( PhotonCandidate->Pt() > 0.299 && PhotonCandidate->Pt() < 16.){ + // fPtGamma = PhotonCandidate->Pt(); + // fDCAzPhoton = PhotonCandidate->GetDCAzToPrimVtx(); + // fRConvPhoton = PhotonCandidate->GetConversionRadius(); + // fEtaPhoton = PhotonCandidate->GetPhotonEta(); + // fCharCatPhoton = PhotonCandidate->GetPhotonQuality(); + // fTreeConvGammaPtDcazCat[fiCut]->Fill(); + // } + // } + + } + +} + +//________________________________________________________________________ +void AliAnalysisTaskGammaConvCalo::ProcessTrueClusterCandidates(AliAODConversionPhoton *TruePhotonCandidate) +{ + + TParticle *Photon = TruePhotonCandidate->GetMCParticle(fMCStack); + + if(Photon == NULL){ + // cout << "no photon" << endl; + return; + } + + if(Photon->GetPdgCode() != 22){ + // cout << "pdg code: " << Photon->GetPdgCode() << endl; + return; // Particle is no Photon + } + + // True Photon + if(fIsFromMBHeader){ + fHistoTrueClusGammaPt[fiCut]->Fill(TruePhotonCandidate->Pt()); + if (fDoPhotonQA > 0) fHistoTrueConvGammaEta[fiCut]->Fill(TruePhotonCandidate->Eta()); + } + + if(Photon->GetMother(0) <= fMCStack->GetNprimary()){ + // Count just primary MC Gammas as true --> For Ratio esdtruegamma / mcconvgamma + if(fIsFromMBHeader){ + fHistoTruePrimaryClusGammaPt[fiCut]->Fill(TruePhotonCandidate->Pt()); + fHistoTruePrimaryClusGammaESDPtMCPt[fiCut]->Fill(TruePhotonCandidate->Pt(),Photon->Pt()); // Allways Filled + } + } + + // maybe later + // else{ + // if(fIsFromMBHeader){ + // fCharPhotonMCInfo = 2; + // fHistoTrueSecondaryConvGammaPt[fiCut]->Fill(TruePhotonCandidate->Pt()); + // if(fMCStack->Particle(Photon->GetMother(0))->GetMother(0) > -1 && + // fMCStack->Particle(fMCStack->Particle(Photon->GetMother(0))->GetMother(0))->GetPdgCode() == 3122){ + // fHistoTrueSecondaryConvGammaFromXFromLambdaPt[fiCut]->Fill(TruePhotonCandidate->Pt()); + // fCharPhotonMCInfo = 5; + // } + // if(fMCStack->Particle(Photon->GetMother(0))->GetMother(0) > -1 && + // fMCStack->Particle(fMCStack->Particle(Photon->GetMother(0))->GetMother(0))->GetPdgCode() == 310){ + // fHistoTrueSecondaryConvGammaFromXFromK0sPt[fiCut]->Fill(TruePhotonCandidate->Pt()); + // fCharPhotonMCInfo = 4; + // } + // if(fMCStack->Particle(Photon->GetMother(0))->GetMother(0) > -1 && + // fMCStack->Particle(fMCStack->Particle(Photon->GetMother(0))->GetMother(0))->GetPdgCode() == 221){ + // fCharPhotonMCInfo = 3; + // } + // } + // } +} + + +//________________________________________________________________________ +void AliAnalysisTaskGammaConvCalo::ProcessPhotonCandidates() +{ + Int_t nV0 = 0; + TList *GammaCandidatesStepOne = new TList(); + TList *GammaCandidatesStepTwo = new TList(); + // Loop over Photon Candidates allocated by ReaderV1 + for(Int_t i = 0; i < fReaderGammas->GetEntriesFast(); i++){ + AliAODConversionPhoton* PhotonCandidate = (AliAODConversionPhoton*) fReaderGammas->At(i); + if(!PhotonCandidate) continue; + fIsFromMBHeader = kTRUE; + if(fIsMC && ((AliConversionCuts*)fCutArray->At(fiCut))->GetSignalRejection() != 0){ + Int_t isPosFromMBHeader + = ((AliConversionCuts*)fCutArray->At(fiCut))->IsParticleFromBGEvent(PhotonCandidate->GetMCLabelPositive(), fMCStack, fInputEvent); + if(isPosFromMBHeader == 0 && ((AliConversionCuts*)fCutArray->At(fiCut))->GetSignalRejection() != 3) continue; + Int_t isNegFromMBHeader + = ((AliConversionCuts*)fCutArray->At(fiCut))->IsParticleFromBGEvent(PhotonCandidate->GetMCLabelNegative(), fMCStack, fInputEvent); + if(isNegFromMBHeader == 0 && ((AliConversionCuts*)fCutArray->At(fiCut))->GetSignalRejection() != 3) continue; + + if( (isNegFromMBHeader+isPosFromMBHeader) != 4) fIsFromMBHeader = kFALSE; + } + + if(!((AliConversionCuts*)fCutArray->At(fiCut))->PhotonIsSelected(PhotonCandidate,fInputEvent)) continue; + if(!((AliConversionCuts*)fCutArray->At(fiCut))->InPlaneOutOfPlaneCut(PhotonCandidate->GetPhotonPhi(),fEventPlaneAngle)) continue; + if(!((AliConversionCuts*)fCutArray->At(fiCut))->UseElecSharingCut() && + !((AliConversionCuts*)fCutArray->At(fiCut))->UseToCloseV0sCut()){ + fGammaCandidates->Add(PhotonCandidate); // if no second loop is required add to events good gammas + + if(fIsFromMBHeader){ + fHistoConvGammaPt[fiCut]->Fill(PhotonCandidate->Pt()); + if (fDoPhotonQA > 0){ + fHistoConvGammaR[fiCut]->Fill(PhotonCandidate->GetConversionRadius()); + fHistoConvGammaEta[fiCut]->Fill(PhotonCandidate->Eta()); + } + } + if(fIsMC){ + if(fInputEvent->IsA()==AliESDEvent::Class()) + ProcessTruePhotonCandidates(PhotonCandidate); + if(fInputEvent->IsA()==AliAODEvent::Class()) + ProcessTruePhotonCandidatesAOD(PhotonCandidate); + } + if (fIsFromMBHeader && fDoPhotonQA == 2){ + if (fIsHeavyIon == 1 && PhotonCandidate->Pt() > 0.399 && PhotonCandidate->Pt() < 12.){ + fPtGamma = PhotonCandidate->Pt(); + fDCAzPhoton = PhotonCandidate->GetDCAzToPrimVtx(); + fRConvPhoton = PhotonCandidate->GetConversionRadius(); + fEtaPhoton = PhotonCandidate->GetPhotonEta(); + fCharCatPhoton = PhotonCandidate->GetPhotonQuality(); + fTreeConvGammaPtDcazCat[fiCut]->Fill(); + } else if ( PhotonCandidate->Pt() > 0.299 && PhotonCandidate->Pt() < 16.){ + fPtGamma = PhotonCandidate->Pt(); + fDCAzPhoton = PhotonCandidate->GetDCAzToPrimVtx(); + fRConvPhoton = PhotonCandidate->GetConversionRadius(); + fEtaPhoton = PhotonCandidate->GetPhotonEta(); + fCharCatPhoton = PhotonCandidate->GetPhotonQuality(); + fTreeConvGammaPtDcazCat[fiCut]->Fill(); + } + } + } else if(((AliConversionCuts*)fCutArray->At(fiCut))->UseElecSharingCut()){ // if Shared Electron cut is enabled, Fill array, add to step one + ((AliConversionCuts*)fCutArray->At(fiCut))->FillElectonLabelArray(PhotonCandidate,nV0); + nV0++; + GammaCandidatesStepOne->Add(PhotonCandidate); + } else if(!((AliConversionCuts*)fCutArray->At(fiCut))->UseElecSharingCut() && + ((AliConversionCuts*)fCutArray->At(fiCut))->UseToCloseV0sCut()){ // shared electron is disabled, step one not needed -> step two + GammaCandidatesStepTwo->Add(PhotonCandidate); + } + } + if(((AliConversionCuts*)fCutArray->At(fiCut))->UseElecSharingCut()){ + for(Int_t i = 0;iGetEntries();i++){ + AliAODConversionPhoton *PhotonCandidate= (AliAODConversionPhoton*) GammaCandidatesStepOne->At(i); + if(!PhotonCandidate) continue; + fIsFromMBHeader = kTRUE; + if(fMCStack && ((AliConversionCuts*)fCutArray->At(fiCut))->GetSignalRejection() != 0){ + Int_t isPosFromMBHeader + = ((AliConversionCuts*)fCutArray->At(fiCut))->IsParticleFromBGEvent(PhotonCandidate->GetMCLabelPositive(), fMCStack, fInputEvent); + Int_t isNegFromMBHeader + = ((AliConversionCuts*)fCutArray->At(fiCut))->IsParticleFromBGEvent(PhotonCandidate->GetMCLabelNegative(), fMCStack, fInputEvent); + if( (isNegFromMBHeader+isPosFromMBHeader) != 4) fIsFromMBHeader = kFALSE; + } + if(!((AliConversionCuts*)fCutArray->At(fiCut))->RejectSharedElectronV0s(PhotonCandidate,i,GammaCandidatesStepOne->GetEntries())) continue; + if(!((AliConversionCuts*)fCutArray->At(fiCut))->UseToCloseV0sCut()){ // To Colse v0s cut diabled, step two not needed + fGammaCandidates->Add(PhotonCandidate); + if(fIsFromMBHeader){ + fHistoConvGammaPt[fiCut]->Fill(PhotonCandidate->Pt()); + if (fDoPhotonQA > 0){ + fHistoConvGammaR[fiCut]->Fill(PhotonCandidate->GetConversionRadius()); + fHistoConvGammaEta[fiCut]->Fill(PhotonCandidate->Eta()); + } + } + } + if(fIsMC){ + if(fInputEvent->IsA()==AliESDEvent::Class()) + ProcessTruePhotonCandidates(PhotonCandidate); + if(fInputEvent->IsA()==AliAODEvent::Class()) + ProcessTruePhotonCandidatesAOD(PhotonCandidate); + } else GammaCandidatesStepTwo->Add(PhotonCandidate); // Close v0s cut enabled -> add to list two + + if (fIsFromMBHeader && fDoPhotonQA == 2){ + if (fIsHeavyIon ==1 && PhotonCandidate->Pt() > 0.399 && PhotonCandidate->Pt() < 12.){ + fPtGamma = PhotonCandidate->Pt(); + fDCAzPhoton = PhotonCandidate->GetDCAzToPrimVtx(); + fRConvPhoton = PhotonCandidate->GetConversionRadius(); + fEtaPhoton = PhotonCandidate->GetPhotonEta(); + fCharCatPhoton = PhotonCandidate->GetPhotonQuality(); + fTreeConvGammaPtDcazCat[fiCut]->Fill(); + } else if ( PhotonCandidate->Pt() > 0.299 && PhotonCandidate->Pt() < 16.){ + fPtGamma = PhotonCandidate->Pt(); + fDCAzPhoton = PhotonCandidate->GetDCAzToPrimVtx(); + fRConvPhoton = PhotonCandidate->GetConversionRadius(); + fEtaPhoton = PhotonCandidate->GetPhotonEta(); + fCharCatPhoton = PhotonCandidate->GetPhotonQuality(); + fTreeConvGammaPtDcazCat[fiCut]->Fill(); + } + } + } + } + if(((AliConversionCuts*)fCutArray->At(fiCut))->UseToCloseV0sCut()){ + for(Int_t i = 0;iGetEntries();i++){ + AliAODConversionPhoton* PhotonCandidate = (AliAODConversionPhoton*) GammaCandidatesStepTwo->At(i); + if(!PhotonCandidate) continue; + fIsFromMBHeader = kTRUE; + if(fMCStack && ((AliConversionCuts*)fCutArray->At(fiCut))->GetSignalRejection() != 0){ + Int_t isPosFromMBHeader + = ((AliConversionCuts*)fCutArray->At(fiCut))->IsParticleFromBGEvent(PhotonCandidate->GetMCLabelPositive(), fMCStack, fInputEvent); + Int_t isNegFromMBHeader + = ((AliConversionCuts*)fCutArray->At(fiCut))->IsParticleFromBGEvent(PhotonCandidate->GetMCLabelNegative(), fMCStack, fInputEvent); + if( (isNegFromMBHeader+isPosFromMBHeader) != 4) fIsFromMBHeader = kFALSE; + } + if(!((AliConversionCuts*)fCutArray->At(fiCut))->RejectToCloseV0s(PhotonCandidate,GammaCandidatesStepTwo,i)) continue; + fGammaCandidates->Add(PhotonCandidate); // Add gamma to current cut TList + if(fIsFromMBHeader){ + fHistoConvGammaPt[fiCut]->Fill(PhotonCandidate->Pt()); + if (fDoPhotonQA > 0){ + fHistoConvGammaR[fiCut]->Fill(PhotonCandidate->GetConversionRadius()); + fHistoConvGammaEta[fiCut]->Fill(PhotonCandidate->Eta()); + } + } + if(fIsMC){ + if(fInputEvent->IsA()==AliESDEvent::Class()) + ProcessTruePhotonCandidates(PhotonCandidate); + if(fInputEvent->IsA()==AliAODEvent::Class()) + ProcessTruePhotonCandidatesAOD(PhotonCandidate); + } + if (fIsFromMBHeader){ + if (fIsHeavyIon == 1 && PhotonCandidate->Pt() > 0.399 && PhotonCandidate->Pt() < 12.){ + fPtGamma = PhotonCandidate->Pt(); + fDCAzPhoton = PhotonCandidate->GetDCAzToPrimVtx(); + fRConvPhoton = PhotonCandidate->GetConversionRadius(); + fEtaPhoton = PhotonCandidate->GetPhotonEta(); + fCharCatPhoton = PhotonCandidate->GetPhotonQuality(); + fTreeConvGammaPtDcazCat[fiCut]->Fill(); + } else if ( PhotonCandidate->Pt() > 0.299 && PhotonCandidate->Pt() < 16.){ + fPtGamma = PhotonCandidate->Pt(); + fDCAzPhoton = PhotonCandidate->GetDCAzToPrimVtx(); + fRConvPhoton = PhotonCandidate->GetConversionRadius(); + fEtaPhoton = PhotonCandidate->GetPhotonEta(); + fCharCatPhoton = PhotonCandidate->GetPhotonQuality(); + fTreeConvGammaPtDcazCat[fiCut]->Fill(); + } + } + } + } + + delete GammaCandidatesStepOne; + GammaCandidatesStepOne = 0x0; + delete GammaCandidatesStepTwo; + GammaCandidatesStepTwo = 0x0; + +} +//________________________________________________________________________ +void AliAnalysisTaskGammaConvCalo::ProcessTruePhotonCandidatesAOD(AliAODConversionPhoton *TruePhotonCandidate) +{ + + Double_t magField = fInputEvent->GetMagneticField(); + if( magField < 0.0 ){ + magField = 1.0; + } + else { + magField = -1.0; + } + + TClonesArray *AODMCTrackArray = dynamic_cast(fInputEvent->FindListObject(AliAODMCParticle::StdBranchName())); + AliAODMCParticle *posDaughter = (AliAODMCParticle*) AODMCTrackArray->At(TruePhotonCandidate->GetMCLabelPositive()); + AliAODMCParticle *negDaughter = (AliAODMCParticle*) AODMCTrackArray->At(TruePhotonCandidate->GetMCLabelNegative()); + fCharPhotonMCInfo = 0; + + if(posDaughter == NULL || negDaughter == NULL) return; // One particle does not exist + Int_t pdgCode[2] = {abs(posDaughter->GetPdgCode()),abs(negDaughter->GetPdgCode())}; + + if(posDaughter->GetMother() != negDaughter->GetMother()){ + FillPhotonCombinatorialBackgroundHist(TruePhotonCandidate, pdgCode); + fCharPhotonMCInfo = 1; + return; + } + else if(posDaughter->GetMother() == -1){ + FillPhotonCombinatorialBackgroundHist(TruePhotonCandidate, pdgCode); + fCharPhotonMCInfo = 1; + return; + } + + if(pdgCode[0]!=11 || pdgCode[1]!=11){ + fCharPhotonMCInfo = 1; + return; //One Particle is not a electron + } + if(posDaughter->GetPdgCode()==negDaughter->GetPdgCode()){ + fCharPhotonMCInfo = 1; + return; // Same Charge + } + + AliAODMCParticle *Photon = (AliAODMCParticle*) AODMCTrackArray->At(posDaughter->GetMother()); + AliVTrack * electronCandidate = ((AliConversionCuts*)fCutArray->At(fiCut))->GetTrack(fInputEvent,TruePhotonCandidate->GetTrackLabelNegative() ); + AliVTrack * positronCandidate = ((AliConversionCuts*)fCutArray->At(fiCut))->GetTrack(fInputEvent,TruePhotonCandidate->GetTrackLabelPositive() ); + Double_t deltaPhi = magField * TVector2::Phi_mpi_pi( electronCandidate->Phi()-positronCandidate->Phi()); + + if(Photon->GetPdgCode() != 22){ + fHistoTrueDalitzPsiPairDeltaPhi[fiCut]->Fill(deltaPhi,TruePhotonCandidate->GetPsiPair()); + fCharPhotonMCInfo = 1; + return; // Mother is no Photon + } + + if(((posDaughter->GetMCProcessCode())) != 5 || ((negDaughter->GetMCProcessCode())) != 5){ + fCharPhotonMCInfo = 1; + return;// check if the daughters come from a conversion + } + // STILL A BUG IN ALIROOT >>8 HAS TPO BE REMOVED AFTER FIX + + + + // True Photon + if(fIsFromMBHeader){ + fHistoTrueConvGammaPt[fiCut]->Fill(TruePhotonCandidate->Pt()); + if (fDoPhotonQA > 0) fHistoTrueConvGammaEta[fiCut]->Fill(TruePhotonCandidate->Eta()); + } + fHistoTrueGammaPsiPairDeltaPhi[fiCut]->Fill(deltaPhi,TruePhotonCandidate->GetPsiPair()); + if(Photon->IsPrimary()){ + // Count just primary MC Gammas as true --> For Ratio esdtruegamma / mcconvgamma + if(fIsFromMBHeader){ + fCharPhotonMCInfo = 6; + fHistoTruePrimaryConvGammaPt[fiCut]->Fill(TruePhotonCandidate->Pt()); + fHistoTruePrimaryConvGammaESDPtMCPt[fiCut]->Fill(TruePhotonCandidate->Pt(),Photon->Pt()); // Allways Filled + } + // (Not Filled for i6, Extra Signal Gamma (parambox) are secondary) + } else { + if(fIsFromMBHeader){ + fHistoTrueSecondaryConvGammaPt[fiCut]->Fill(TruePhotonCandidate->Pt()); + fCharPhotonMCInfo = 2; + if(((AliAODMCParticle*)AODMCTrackArray->At(Photon->GetMother()))->GetMother() > -1 && + ((AliAODMCParticle*)AODMCTrackArray->At(((AliAODMCParticle*)AODMCTrackArray->At(Photon->GetMother()))->GetMother()))->GetPdgCode() == 3122){ + fCharPhotonMCInfo = 5; + fHistoTrueSecondaryConvGammaFromXFromLambdaPt[fiCut]->Fill(TruePhotonCandidate->Pt()); + } + if(((AliAODMCParticle*)AODMCTrackArray->At(Photon->GetMother()))->GetMother() > -1 && + ((AliAODMCParticle*)AODMCTrackArray->At(((AliAODMCParticle*)AODMCTrackArray->At(Photon->GetMother()))->GetMother()))->GetPdgCode() == 310){ + fCharPhotonMCInfo = 4; + fHistoTrueSecondaryConvGammaFromXFromK0sPt[fiCut]->Fill(TruePhotonCandidate->Pt()); + } + if(((AliAODMCParticle*)AODMCTrackArray->At(Photon->GetMother()))->GetMother() > -1 && + ((AliAODMCParticle*)AODMCTrackArray->At(((AliAODMCParticle*)AODMCTrackArray->At(Photon->GetMother()))->GetMother()))->GetPdgCode() == 221){ + fCharPhotonMCInfo = 3; + } + } + } + +} +//________________________________________________________________________ +void AliAnalysisTaskGammaConvCalo::ProcessTruePhotonCandidates(AliAODConversionPhoton *TruePhotonCandidate) +{ + + Double_t magField = fInputEvent->GetMagneticField(); + if( magField < 0.0 ){ + magField = 1.0; + } + else { + magField = -1.0; + } + + // Process True Photons + TParticle *posDaughter = TruePhotonCandidate->GetPositiveMCDaughter(fMCStack); + TParticle *negDaughter = TruePhotonCandidate->GetNegativeMCDaughter(fMCStack); + fCharPhotonMCInfo = 0; + + if(posDaughter == NULL || negDaughter == NULL) return; // One particle does not exist + Int_t pdgCode[2] = {abs(posDaughter->GetPdgCode()),abs(negDaughter->GetPdgCode())}; + fCharPhotonMCInfo = 1; + if(posDaughter->GetMother(0) != negDaughter->GetMother(0)){ + FillPhotonCombinatorialBackgroundHist(TruePhotonCandidate, pdgCode); + return; + } + else if(posDaughter->GetMother(0) == -1){ + FillPhotonCombinatorialBackgroundHist(TruePhotonCandidate, pdgCode); + return; + } + + if(pdgCode[0]!=11 || pdgCode[1]!=11) return; //One Particle is not a electron + + if(posDaughter->GetPdgCode()==negDaughter->GetPdgCode()) return; // Same Charge + + TParticle *Photon = TruePhotonCandidate->GetMCParticle(fMCStack); + AliVTrack * electronCandidate = ((AliConversionCuts*)fCutArray->At(fiCut))->GetTrack(fInputEvent,TruePhotonCandidate->GetTrackLabelNegative() ); + AliVTrack * positronCandidate = ((AliConversionCuts*)fCutArray->At(fiCut))->GetTrack(fInputEvent,TruePhotonCandidate->GetTrackLabelPositive() ); + Double_t deltaPhi = magField * TVector2::Phi_mpi_pi( electronCandidate->Phi()-positronCandidate->Phi()); + + if(Photon->GetPdgCode() != 22){ + fHistoTrueDalitzPsiPairDeltaPhi[fiCut]->Fill(deltaPhi,TruePhotonCandidate->GetPsiPair()); + return; // Mother is no Photon + } + + if(posDaughter->GetUniqueID() != 5 || negDaughter->GetUniqueID() !=5) return;// check if the daughters come from a conversion + + + + // True Photon + if(fIsFromMBHeader){ + fHistoTrueConvGammaPt[fiCut]->Fill(TruePhotonCandidate->Pt()); + if (fDoPhotonQA > 0) fHistoTrueConvGammaEta[fiCut]->Fill(TruePhotonCandidate->Eta()); + } + fHistoTrueGammaPsiPairDeltaPhi[fiCut]->Fill(deltaPhi,TruePhotonCandidate->GetPsiPair()); + if(posDaughter->GetMother(0) <= fMCStack->GetNprimary()){ + // Count just primary MC Gammas as true --> For Ratio esdtruegamma / mcconvgamma + if(fIsFromMBHeader){ + fCharPhotonMCInfo = 6; + fHistoTruePrimaryConvGammaPt[fiCut]->Fill(TruePhotonCandidate->Pt()); + fHistoTruePrimaryConvGammaESDPtMCPt[fiCut]->Fill(TruePhotonCandidate->Pt(),Photon->Pt()); // Allways Filled + + } + // (Not Filled for i6, Extra Signal Gamma (parambox) are secondary) + } + else{ + if(fIsFromMBHeader){ + fCharPhotonMCInfo = 2; + fHistoTrueSecondaryConvGammaPt[fiCut]->Fill(TruePhotonCandidate->Pt()); + if(fMCStack->Particle(Photon->GetMother(0))->GetMother(0) > -1 && + fMCStack->Particle(fMCStack->Particle(Photon->GetMother(0))->GetMother(0))->GetPdgCode() == 3122){ + fHistoTrueSecondaryConvGammaFromXFromLambdaPt[fiCut]->Fill(TruePhotonCandidate->Pt()); + fCharPhotonMCInfo = 5; + } + if(fMCStack->Particle(Photon->GetMother(0))->GetMother(0) > -1 && + fMCStack->Particle(fMCStack->Particle(Photon->GetMother(0))->GetMother(0))->GetPdgCode() == 310){ + fHistoTrueSecondaryConvGammaFromXFromK0sPt[fiCut]->Fill(TruePhotonCandidate->Pt()); + fCharPhotonMCInfo = 4; + } + if(fMCStack->Particle(Photon->GetMother(0))->GetMother(0) > -1 && + fMCStack->Particle(fMCStack->Particle(Photon->GetMother(0))->GetMother(0))->GetPdgCode() == 221){ + fCharPhotonMCInfo = 3; + } + } + } + + // pi0 photon + //Bool_t bpi0 = 0; + Int_t imother = Photon->GetMother(0); + AliMCParticle *McMother = static_cast(fMCEvent->GetTrack(imother)); + if(McMother->PdgCode() == 111){ + fHistoTrueConvPi0GammaPt[fiCut]->Fill(TruePhotonCandidate->Pt()); + } + +} +//________________________________________________________________________ +void AliAnalysisTaskGammaConvCalo::ProcessAODMCParticles() +{ + + TClonesArray *AODMCTrackArray = dynamic_cast(fInputEvent->FindListObject(AliAODMCParticle::StdBranchName())); + + // Loop over all primary MC particle + for(Int_t i = 0; i < AODMCTrackArray->GetEntriesFast(); i++) { + + AliAODMCParticle* particle = static_cast(AODMCTrackArray->At(i)); + if (!particle) continue; + if (!particle->IsPrimary()) continue; + + Int_t isMCFromMBHeader = -1; + if(((AliConversionCuts*)fCutArray->At(fiCut))->GetSignalRejection() != 0){ + isMCFromMBHeader + = ((AliConversionCuts*)fCutArray->At(fiCut))->IsParticleFromBGEvent(i, fMCStack, fInputEvent); + if(isMCFromMBHeader == 0 && ((AliConversionCuts*)fCutArray->At(fiCut))->GetSignalRejection() != 3) continue; + } + + if(!((AliConversionCuts*)fCutArray->At(fiCut))->InPlaneOutOfPlaneCut(particle->Phi(),fEventPlaneAngle,kFALSE)) continue; + if(((AliConversionCuts*)fCutArray->At(fiCut))->PhotonIsSelectedAODMC(particle,AODMCTrackArray,kFALSE)){ + fHistoMCAllGammaPt[fiCut]->Fill(particle->Pt()); // All MC Gamma + if(particle->GetMother() >-1){ // Meson Decay Gamma + switch((static_cast(AODMCTrackArray->At(particle->GetMother())))->GetPdgCode()){ + case 111: // Pi0 + fHistoMCDecayGammaPi0Pt[fiCut]->Fill(particle->Pt()); + break; + case 113: // Rho0 + fHistoMCDecayGammaRhoPt[fiCut]->Fill(particle->Pt()); + break; + case 221: // Eta + fHistoMCDecayGammaEtaPt[fiCut]->Fill(particle->Pt()); + break; + case 223: // Omega + fHistoMCDecayGammaOmegaPt[fiCut]->Fill(particle->Pt()); + break; + case 331: // Eta' + fHistoMCDecayGammaEtapPt[fiCut]->Fill(particle->Pt()); + break; + case 333: // Phi + fHistoMCDecayGammaPhiPt[fiCut]->Fill(particle->Pt()); + break; + case 3212: // Sigma + fHistoMCDecayGammaSigmaPt[fiCut]->Fill(particle->Pt()); + break; + } + } + } + if(((AliConversionCuts*)fCutArray->At(fiCut))->PhotonIsSelectedAODMC(particle,AODMCTrackArray,kTRUE)){ + Double_t rConv = 0; + for(Int_t daughterIndex=particle->GetDaughter(0);daughterIndex<=particle->GetDaughter(1);daughterIndex++){ + AliAODMCParticle *tmpDaughter = static_cast(AODMCTrackArray->At(daughterIndex)); + if(!tmpDaughter) continue; + if(abs(tmpDaughter->GetPdgCode()) == 11){ + rConv = sqrt( (tmpDaughter->Xv()*tmpDaughter->Xv()) + (tmpDaughter->Yv()*tmpDaughter->Yv()) ); + } + } + fHistoMCConvGammaPt[fiCut]->Fill(particle->Pt()); + if (fDoPhotonQA > 0){ + fHistoMCConvGammaR[fiCut]->Fill(rConv); + fHistoMCConvGammaEta[fiCut]->Fill(particle->Eta()); + } + } + // Converted MC Gamma + if(fDoMesonAnalysis){ + if(particle->GetPdgCode() == 310 && fDoMesonQA > 0){ + Double_t mesonY = 10.; + if(particle->E() - particle->Pz() == 0 || particle->E() + particle->Pz() == 0){ + mesonY=10.-((AliConversionCuts*)fCutArray->At(fiCut))->GetEtaShift(); + } else { + mesonY = 0.5*(TMath::Log((particle->E()+particle->Pz()) / (particle->E()-particle->Pz())))-((AliConversionCuts*)fCutArray->At(fiCut))->GetEtaShift(); + } + Float_t weightedK0s= 1; + if(((AliConversionCuts*)fCutArray->At(fiCut))->IsParticleFromBGEvent(i, fMCStack, fInputEvent)){ + if (particle->Pt()>0.005){ + weightedK0s= ((AliConversionCuts*)fCutArray->At(fiCut))->GetWeightForMeson(fV0Reader->GetPeriodName(),i, 0x0, fInputEvent); + //cout << "MC input \t"<Pt()<<"\t"<Fill(particle->Pt(),weightedK0s); + fHistoMCK0sWOWeightPt[fiCut]->Fill(particle->Pt()); + fHistoMCK0sPtY[fiCut]->Fill(particle->Pt(),mesonY,weightedK0s); + } + if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut)) + ->MesonIsSelectedAODMC(particle,AODMCTrackArray,((AliConversionCuts*)fCutArray->At(fiCut))->GetEtaShift())){ + AliAODMCParticle* daughter0 = static_cast(AODMCTrackArray->At(particle->GetDaughter(0))); + AliAODMCParticle* daughter1 = static_cast(AODMCTrackArray->At(particle->GetDaughter(1))); + Float_t weighted= 1; + if(((AliConversionCuts*)fCutArray->At(fiCut))->IsParticleFromBGEvent(i, fMCStack, fInputEvent)){ + if (particle->Pt()>0.005){ + weighted= ((AliConversionCuts*)fCutArray->At(fiCut))->GetWeightForMeson(fV0Reader->GetPeriodName(),i, 0x0, fInputEvent); + // if(particle->GetPdgCode() == 221){ + // cout << "MC input \t"<Pt()<<"\t"<E() - particle->Pz() == 0 || particle->E() + particle->Pz() == 0){ + mesonY=10.-((AliConversionCuts*)fCutArray->At(fiCut))->GetEtaShift(); + } else{ + mesonY = 0.5*(TMath::Log((particle->E()+particle->Pz()) / (particle->E()-particle->Pz())))-((AliConversionCuts*)fCutArray->At(fiCut))->GetEtaShift(); + } + + if(particle->GetPdgCode() == 111){ + fHistoMCPi0Pt[fiCut]->Fill(particle->Pt(),weighted); // All MC Pi0 + fHistoMCPi0WOWeightPt[fiCut]->Fill(particle->Pt()); + if (fDoMesonQA > 0) fHistoMCPi0PtY[fiCut]->Fill(particle->Pt(),mesonY,weighted); // All MC Pi0 + } else if(particle->GetPdgCode() == 221){ + fHistoMCEtaPt[fiCut]->Fill(particle->Pt(),weighted); // All MC Eta + fHistoMCEtaWOWeightPt[fiCut]->Fill(particle->Pt()); + if (fDoMesonQA > 0) fHistoMCEtaPtY[fiCut]->Fill(particle->Pt(),mesonY,weighted); // All MC Pi0 + } + + // Check the acceptance for both gammas + if(((AliConversionCuts*)fCutArray->At(fiCut))->PhotonIsSelectedAODMC(daughter0,AODMCTrackArray,kFALSE) && + ((AliConversionCuts*)fCutArray->At(fiCut))->PhotonIsSelectedAODMC(daughter1,AODMCTrackArray,kFALSE) && + ((AliConversionCuts*)fCutArray->At(fiCut))->InPlaneOutOfPlaneCut(daughter0->Phi(),fEventPlaneAngle,kFALSE) && + ((AliConversionCuts*)fCutArray->At(fiCut))->InPlaneOutOfPlaneCut(daughter1->Phi(),fEventPlaneAngle,kFALSE)){ + + if(particle->GetPdgCode() == 111){ + fHistoMCPi0InAccPt[fiCut]->Fill(particle->Pt(),weighted); // MC Pi0 with gamma in acc + } else if(particle->GetPdgCode() == 221){ + fHistoMCEtaInAccPt[fiCut]->Fill(particle->Pt(),weighted); // MC Eta with gamma in acc + } + } + } + } + } + +} +//________________________________________________________________________ +void AliAnalysisTaskGammaConvCalo::ProcessMCParticles() +{ + // Loop over all primary MC particle + for(Int_t i = 0; i < fMCStack->GetNprimary(); i++) { + TParticle* particle = (TParticle *)fMCStack->Particle(i); + if (!particle) continue; + + Int_t isMCFromMBHeader = -1; + if(((AliConversionCuts*)fCutArray->At(fiCut))->GetSignalRejection() != 0){ + isMCFromMBHeader + = ((AliConversionCuts*)fCutArray->At(fiCut))->IsParticleFromBGEvent(i, fMCStack, fInputEvent); + if(isMCFromMBHeader == 0 && ((AliConversionCuts*)fCutArray->At(fiCut))->GetSignalRejection() != 3) continue; + } + + if(!((AliConversionCuts*)fCutArray->At(fiCut))->InPlaneOutOfPlaneCut(particle->Phi(),fEventPlaneAngle,kFALSE)) continue; + if(((AliConversionCuts*)fCutArray->At(fiCut))->PhotonIsSelectedMC(particle,fMCStack,kFALSE)){ + fHistoMCAllGammaPt[fiCut]->Fill(particle->Pt()); // All MC Gamma + if(particle->GetMother(0) >-1){ // Meson Decay Gamma + switch(fMCStack->Particle(particle->GetMother(0))->GetPdgCode()){ + case 111: // Pi0 + fHistoMCDecayGammaPi0Pt[fiCut]->Fill(particle->Pt()); + break; + case 113: // Rho0 + fHistoMCDecayGammaRhoPt[fiCut]->Fill(particle->Pt()); + break; + case 221: // Eta + fHistoMCDecayGammaEtaPt[fiCut]->Fill(particle->Pt()); + break; + case 223: // Omega + fHistoMCDecayGammaOmegaPt[fiCut]->Fill(particle->Pt()); + break; + case 331: // Eta' + fHistoMCDecayGammaEtapPt[fiCut]->Fill(particle->Pt()); + break; + case 333: // Phi + fHistoMCDecayGammaPhiPt[fiCut]->Fill(particle->Pt()); + break; + case 3212: // Sigma + fHistoMCDecayGammaSigmaPt[fiCut]->Fill(particle->Pt()); + break; + } + } + } + if(((AliConversionCuts*)fCutArray->At(fiCut))->PhotonIsSelectedMC(particle,fMCStack,kTRUE)){ + fHistoMCConvGammaPt[fiCut]->Fill(particle->Pt()); + if (fDoPhotonQA > 0){ + fHistoMCConvGammaR[fiCut]->Fill(((TParticle*)fMCStack->Particle(particle->GetFirstDaughter()))->R()); + fHistoMCConvGammaEta[fiCut]->Fill(particle->Eta()); + } + } // Converted MC Gamma + if(fDoMesonAnalysis){ + if(particle->GetPdgCode() == 310 && fDoMesonQA > 0){ + Double_t mesonY = 10.; + if(particle->Energy() - particle->Pz() == 0 || particle->Energy() + particle->Pz() == 0){ + mesonY=10.-((AliConversionCuts*)fCutArray->At(fiCut))->GetEtaShift(); + } else{ + mesonY = 0.5*(TMath::Log((particle->Energy()+particle->Pz()) / (particle->Energy()-particle->Pz())))-((AliConversionCuts*)fCutArray->At(fiCut))->GetEtaShift(); + } + Float_t weightedK0s= 1; + if(((AliConversionCuts*)fCutArray->At(fiCut))->IsParticleFromBGEvent(i, fMCStack, fInputEvent)){ + if (particle->Pt()>0.005){ + weightedK0s= ((AliConversionCuts*)fCutArray->At(fiCut))->GetWeightForMeson(fV0Reader->GetPeriodName(),i, fMCStack, fInputEvent); + //cout << "MC input \t"<Pt()<<"\t"<IsPhysicalPrimary(i)){ + fHistoMCK0sPt[fiCut]->Fill(particle->Pt(),weightedK0s); + fHistoMCK0sWOWeightPt[fiCut]->Fill(particle->Pt()); + fHistoMCK0sPtY[fiCut]->Fill(particle->Pt(),mesonY,weightedK0s); + } + } + if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut)) + ->MesonIsSelectedMC(particle,fMCStack,((AliConversionCuts*)fCutArray->At(fiCut))->GetEtaShift())){ + TParticle* daughter0 = (TParticle*)fMCStack->Particle(particle->GetFirstDaughter()); + TParticle* daughter1 = (TParticle*)fMCStack->Particle(particle->GetLastDaughter()); + + Float_t weighted= 1; + if(((AliConversionCuts*)fCutArray->At(fiCut))->IsParticleFromBGEvent(i, fMCStack, fInputEvent)){ + if (particle->Pt()>0.005){ + weighted= ((AliConversionCuts*)fCutArray->At(fiCut))->GetWeightForMeson(fV0Reader->GetPeriodName(),i, fMCStack, fInputEvent); + // if(particle->GetPdgCode() == 221){ + // cout << "MC input \t"<Pt()<<"\t"<Energy() - particle->Pz() == 0 || particle->Energy() + particle->Pz() == 0){ + mesonY=10.-((AliConversionCuts*)fCutArray->At(fiCut))->GetEtaShift(); + } else{ + mesonY = 0.5*(TMath::Log((particle->Energy()+particle->Pz()) / (particle->Energy()-particle->Pz())))-((AliConversionCuts*)fCutArray->At(fiCut))->GetEtaShift(); + } + + if(particle->GetPdgCode() == 111){ + fHistoMCPi0Pt[fiCut]->Fill(particle->Pt(),weighted); // All MC Pi0 + fHistoMCPi0WOWeightPt[fiCut]->Fill(particle->Pt()); + if (fDoMesonQA > 0) fHistoMCPi0PtY[fiCut]->Fill(particle->Pt(),mesonY,weighted); // All MC Pi0 + } else if(particle->GetPdgCode() == 221){ + fHistoMCEtaPt[fiCut]->Fill(particle->Pt(),weighted); // All MC Eta + fHistoMCEtaWOWeightPt[fiCut]->Fill(particle->Pt()); + if (fDoMesonQA > 0) fHistoMCEtaPtY[fiCut]->Fill(particle->Pt(),mesonY,weighted); // All MC Pi0 + } + + // Check the acceptance for both gammas + if(((AliConversionCuts*)fCutArray->At(fiCut))->PhotonIsSelectedMC(daughter0,fMCStack,kFALSE) && + ((AliConversionCuts*)fCutArray->At(fiCut))->PhotonIsSelectedMC(daughter1,fMCStack,kFALSE) && + ((AliConversionCuts*)fCutArray->At(fiCut))->InPlaneOutOfPlaneCut(daughter0->Phi(),fEventPlaneAngle,kFALSE) && + ((AliConversionCuts*)fCutArray->At(fiCut))->InPlaneOutOfPlaneCut(daughter1->Phi(),fEventPlaneAngle,kFALSE)){ + + if(particle->GetPdgCode() == 111){ + fHistoMCPi0InAccPt[fiCut]->Fill(particle->Pt(),weighted); // MC Pi0 with gamma in acc + } else if(particle->GetPdgCode() == 221){ + fHistoMCEtaInAccPt[fiCut]->Fill(particle->Pt(),weighted); // MC Eta with gamma in acc + } + } + } + } + } + + if (fDoMesonQA){ + for(Int_t i = fMCStack->GetNprimary(); i < fMCStack->GetNtrack(); i++) { + TParticle* particle = (TParticle *)fMCStack->Particle(i); + if (!particle) continue; + + Int_t isMCFromMBHeader = -1; + if(((AliConversionCuts*)fCutArray->At(fiCut))->GetSignalRejection() != 0){ + isMCFromMBHeader + = ((AliConversionCuts*)fCutArray->At(fiCut))->IsParticleFromBGEvent(i, fMCStack, fInputEvent); + if(isMCFromMBHeader == 0 && ((AliConversionCuts*)fCutArray->At(fiCut))->GetSignalRejection() != 3) continue; + } + + if(fDoMesonAnalysis){ + if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut)) + ->MesonIsSelectedMC(particle,fMCStack,((AliConversionCuts*)fCutArray->At(fiCut))->GetEtaShift())){ + Float_t weighted= 1; + if(((AliConversionCuts*)fCutArray->At(fiCut))->IsParticleFromBGEvent(i, fMCStack, fInputEvent)){ + if (particle->Pt()>0.005){ + weighted= ((AliConversionCuts*)fCutArray->At(fiCut))->GetWeightForMeson(fV0Reader->GetPeriodName(),i, fMCStack, fInputEvent); + // if(particle->GetPdgCode() == 221){ + // cout << "MC input \t"<Pt()<<"\t"<GetPdgCode() == 111){ + Int_t pdgCode = ((TParticle*)fMCStack->Particle( particle->GetFirstMother() ))->GetPdgCode(); + Int_t source = GetSourceClassification(111,pdgCode); + fHistoMCSecPi0PtvsSource[fiCut]->Fill(particle->Pt(),source,weighted); // All MC Pi0 + fHistoMCSecPi0Source[fiCut]->Fill(pdgCode); + } else if(particle->GetPdgCode() == 221){ + Int_t pdgCode = ((TParticle*)fMCStack->Particle( particle->GetFirstMother() ))->GetPdgCode(); + fHistoMCSecEtaPt[fiCut]->Fill(particle->Pt(),weighted); // All MC Pi0 + fHistoMCSecEtaSource[fiCut]->Fill(pdgCode); + } + } + } + } + } +} + +//________________________________________________________________________ +void AliAnalysisTaskGammaConvCalo::PhotonTagging(){ + + // Conversion Gammas + if(fGammaCandidates->GetEntries()>0){ + + for(Int_t firstGammaIndex=0;firstGammaIndexGetEntries()-1;firstGammaIndex++){ + + // get conversion photon + AliAODConversionPhoton *gamma0=dynamic_cast(fGammaCandidates->At(firstGammaIndex)); + if (gamma0==NULL) continue; + + TLorentzVector photonVector; + photonVector.SetPxPyPzE(gamma0->GetPx(),gamma0->GetPy(),gamma0->GetPz(),gamma0->GetPhotonP()); + + Bool_t btagpi0 = 0; + Bool_t btageta = 0; + + // loop over clusters + for(Int_t secondGammaIndex = 0; secondGammaIndexGetEntries(); ++secondGammaIndex) { + + AliAODConversionPhoton *gamma1=dynamic_cast(fClusterCandidates->At(secondGammaIndex)); + if (gamma1==NULL) continue; + + TLorentzVector clusterVector; + clusterVector.SetPxPyPzE(gamma1->GetPx(),gamma1->GetPy(),gamma1->GetPz(),gamma1->GetPhotonP()); + + // do the tagging + TLorentzVector pairVector = photonVector+clusterVector; + + // see if pi0? + if((pairVector.M() > 0.11 && pairVector.M() < 0.15)){ + btagpi0 = 1; + } + // or eta + if((pairVector.M() > 0.50 && pairVector.M() < 0.6)){ + btageta = 1; + } + }// end loop over clusters + + if(btagpi0 && btageta) + fHistoConvGammaTagged[fiCut]->Fill(photonVector.Pt()); + else if(btagpi0 && !btageta) + fHistoConvGammaPi0Tagged[fiCut]->Fill(photonVector.Pt()); + else if(btageta && !btagpi0) + fHistoConvGammaEtaTagged[fiCut]->Fill(photonVector.Pt()); + else + fHistoConvGammaUntagged[fiCut]->Fill(photonVector.Pt()); + + }// end loop over gammas + }// end if + return; +} + +//________________________________________________________________________ +void AliAnalysisTaskGammaConvCalo::CalculatePi0Candidates(){ + + // Conversion Gammas + if(fGammaCandidates->GetEntries()>0){ + + // vertex + Double_t vertex[3] = {0}; + InputEvent()->GetPrimaryVertex()->GetXYZ(vertex); + + for(Int_t firstGammaIndex=0;firstGammaIndexGetEntries();firstGammaIndex++){ + AliAODConversionPhoton *gamma0=dynamic_cast(fGammaCandidates->At(firstGammaIndex)); + if (gamma0==NULL) continue; + + for(Int_t secondGammaIndex=0;secondGammaIndexGetEntries();secondGammaIndex++){ + + AliAODConversionPhoton *gamma1=dynamic_cast(fClusterCandidates->At(secondGammaIndex)); + if (gamma1==NULL) continue; + + AliAODConversionMother *pi0cand = new AliAODConversionMother(gamma0,gamma1); + pi0cand->SetLabels(firstGammaIndex,secondGammaIndex); + + if((((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->MesonIsSelected(pi0cand,kTRUE,((AliConversionCuts*)fCutArray->At(fiCut))->GetEtaShift()))){ + fHistoMotherInvMassPt[fiCut]->Fill(pi0cand->M(),pi0cand->Pt()); + + // fill new histograms + fHistoPhotonPairAll[fiCut]->Fill(pi0cand->M(),pi0cand->Pt()); + fHistoPhotonPairAllGam[fiCut]->Fill(pi0cand->M(),gamma0->Pt()); + + if(pi0cand->GetAlpha()<0.1) + fHistoMotherInvMassEalpha[fiCut]->Fill(pi0cand->M(),pi0cand->E()); + + if (fDoMesonQA > 0){ + if ( pi0cand->M() > 0.05 && pi0cand->M() < 0.17){ + fHistoMotherPi0PtY[fiCut]->Fill(pi0cand->Pt(),pi0cand->Rapidity()-((AliConversionCuts*)fCutArray->At(fiCut))->GetEtaShift()); + fHistoMotherPi0PtAlpha[fiCut]->Fill(pi0cand->Pt(),pi0cand->GetAlpha()); + fHistoMotherPi0PtOpenAngle[fiCut]->Fill(pi0cand->Pt(),pi0cand->GetOpeningAngle()); + } + if ( pi0cand->M() > 0.45 && pi0cand->M() < 0.65){ + fHistoMotherEtaPtY[fiCut]->Fill(pi0cand->Pt(),pi0cand->Rapidity()-((AliConversionCuts*)fCutArray->At(fiCut))->GetEtaShift()); + fHistoMotherEtaPtAlpha[fiCut]->Fill(pi0cand->Pt(),pi0cand->GetAlpha()); + fHistoMotherEtaPtOpenAngle[fiCut]->Fill(pi0cand->Pt(),pi0cand->GetOpeningAngle()); + } + } + if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->DoBGCalculation()){ + Int_t zbin = 0; + Int_t mbin = 0; + + if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->BackgroundHandlerType() == 0){ + zbin = fBGHandler[fiCut]->GetZBinIndex(fInputEvent->GetPrimaryVertex()->GetZ()); + if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->UseTrackMultiplicity()){ + mbin = fBGHandler[fiCut]->GetMultiplicityBinIndex(fV0Reader->GetNumberOfPrimaryTracks()); + } else { + mbin = fBGHandler[fiCut]->GetMultiplicityBinIndex(fGammaCandidates->GetEntries()); + } + } else{ + zbin = fBGHandlerRP[fiCut]->GetZBinIndex(fInputEvent->GetPrimaryVertex()->GetZ()); + if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->UseTrackMultiplicity()){ + mbin = fBGHandlerRP[fiCut]->GetMultiplicityBinIndex(fV0Reader->GetNumberOfPrimaryTracks()); + } else { + mbin = fBGHandlerRP[fiCut]->GetMultiplicityBinIndex(fGammaCandidates->GetEntries()); + } + } + Double_t sparesFill[4] = {pi0cand->M(),pi0cand->Pt(),(Double_t)zbin,(Double_t)mbin}; + fSparseMotherInvMassPtZM[fiCut]->Fill(sparesFill,1); + } + + + if(fIsMC){ + if(fInputEvent->IsA()==AliESDEvent::Class()) + ProcessTrueMesonCandidates(pi0cand,gamma0,gamma1); + if(fInputEvent->IsA()==AliAODEvent::Class()) + ProcessTrueMesonCandidatesAOD(pi0cand,gamma0,gamma1); + } + if (fDoMesonQA == 2){ + fInvMass = pi0cand->M(); + fPt = pi0cand->Pt(); + if (abs(gamma0->GetDCAzToPrimVtx()) < abs(gamma1->GetDCAzToPrimVtx())){ + fDCAzGammaMin = gamma0->GetDCAzToPrimVtx(); + fDCAzGammaMax = gamma1->GetDCAzToPrimVtx(); + } else { + fDCAzGammaMin = gamma1->GetDCAzToPrimVtx(); + fDCAzGammaMax = gamma0->GetDCAzToPrimVtx(); + } + fCharFlag = pi0cand->GetMesonQuality(); + // cout << "gamma 0: " << gamma0->GetV0Index()<< "\t" << gamma0->GetPx() << "\t" << gamma0->GetPy() << "\t" << gamma0->GetPz() << "\t" << endl; + // cout << "gamma 1: " << gamma1->GetV0Index()<< "\t"<< gamma1->GetPx() << "\t" << gamma1->GetPy() << "\t" << gamma1->GetPz() << "\t" << endl; + // cout << "pi0: "< 0.399 && fPt < 20. ) { + if (fInvMass > 0.08 && fInvMass < 0.2) fTreeMesonsInvMassPtDcazMinDcazMaxFlag[fiCut]->Fill(); + if ((fInvMass > 0.45 && fInvMass < 0.6) && (fPt > 0.999 && fPt < 20.) )fTreeMesonsInvMassPtDcazMinDcazMaxFlag[fiCut]->Fill(); + } else if (fPt > 0.299 && fPt < 20. ) { + if ( (fInvMass > 0.08 && fInvMass < 0.2) || (fInvMass > 0.45 && fInvMass < 0.6)) fTreeMesonsInvMassPtDcazMinDcazMaxFlag[fiCut]->Fill(); + } + } + } + delete pi0cand; + pi0cand=0x0; + } + } + } +} +//______________________________________________________________________ +void AliAnalysisTaskGammaConvCalo::ProcessTrueMesonCandidates(AliAODConversionMother *Pi0Candidate, AliAODConversionPhoton *TrueGammaCandidate0, AliAODConversionPhoton *TrueGammaCandidate1) +{ + // Process True Mesons + AliStack *MCStack = fMCEvent->Stack(); + fCharMesonMCInfo = 0; + if(TrueGammaCandidate0->GetV0Index()GetNumberOfV0s()){ + Bool_t isTruePi0 = kFALSE; + Bool_t isTrueEta = kFALSE; + Bool_t isTruePi0Dalitz = kFALSE; + Bool_t isTrueEtaDalitz = kFALSE; + Bool_t gamma0DalitzCand = kFALSE; + Bool_t gamma1DalitzCand = kFALSE; + Int_t gamma0MCLabel = TrueGammaCandidate0->GetMCParticleLabel(MCStack); + Int_t gamma0MotherLabel = -1; + if(gamma0MCLabel != -1){ // Gamma is Combinatorial; MC Particles don't belong to the same Mother + // Daughters Gamma 0 + TParticle * negativeMC = (TParticle*)TrueGammaCandidate0->GetNegativeMCDaughter(MCStack); + TParticle * positiveMC = (TParticle*)TrueGammaCandidate0->GetPositiveMCDaughter(MCStack); + TParticle * gammaMC0 = (TParticle*)MCStack->Particle(gamma0MCLabel); + if(abs(negativeMC->GetPdgCode())==11 && abs(positiveMC->GetPdgCode())==11){ // Electrons ... + if(negativeMC->GetUniqueID() == 5 && positiveMC->GetUniqueID() ==5){ // ... From Conversion ... + if(gammaMC0->GetPdgCode() == 22){ // ... with Gamma Mother + gamma0MotherLabel=gammaMC0->GetFirstMother(); + } + } + if(gammaMC0->GetPdgCode() ==111){ // Dalitz candidate + gamma0DalitzCand = kTRUE; + gamma0MotherLabel=-111; + } + if(gammaMC0->GetPdgCode() ==221){ // Dalitz candidate + gamma0DalitzCand = kTRUE; + gamma0MotherLabel=-221; + } + } + } + if(TrueGammaCandidate1->GetV0Index()GetNumberOfV0s()){ + Int_t gamma1MCLabel = TrueGammaCandidate1->GetMCParticleLabel(MCStack); + Int_t gamma1MotherLabel = -1; + if(gamma1MCLabel != -1){ // Gamma is Combinatorial; MC Particles don't belong to the same Mother + // Daughters Gamma 1 + TParticle * negativeMC = (TParticle*)TrueGammaCandidate1->GetNegativeMCDaughter(MCStack); + TParticle * positiveMC = (TParticle*)TrueGammaCandidate1->GetPositiveMCDaughter(MCStack); + TParticle * gammaMC1 = (TParticle*)MCStack->Particle(gamma1MCLabel); + if(abs(negativeMC->GetPdgCode())==11 && abs(positiveMC->GetPdgCode())==11){ // Electrons ... + if(negativeMC->GetUniqueID() == 5 && positiveMC->GetUniqueID() ==5){ // ... From Conversion ... + if(gammaMC1->GetPdgCode() == 22){ // ... with Gamma Mother + gamma1MotherLabel=gammaMC1->GetFirstMother(); + } + } + if(gammaMC1->GetPdgCode() ==111 ){ // Dalitz candidate + gamma1DalitzCand = kTRUE; + gamma1MotherLabel=-111; + } + if(gammaMC1->GetPdgCode() ==221){ // Dalitz candidate + gamma1DalitzCand = kTRUE; + gamma1MotherLabel=-221; + } + } + } + if(gamma0MotherLabel>=0 && gamma0MotherLabel==gamma1MotherLabel){ + if(((TParticle*)MCStack->Particle(gamma1MotherLabel))->GetPdgCode() == 111){ + isTruePi0=kTRUE; + } + if(((TParticle*)MCStack->Particle(gamma1MotherLabel))->GetPdgCode() == 221){ + isTrueEta=kTRUE; + } + } + + //Identify Dalitz candidate + if (gamma1DalitzCand || gamma0DalitzCand){ + if (gamma0DalitzCand && gamma0MCLabel >=0 && gamma0MCLabel==gamma1MotherLabel){ + if (gamma0MotherLabel == -111) isTruePi0Dalitz = kTRUE; + if (gamma0MotherLabel == -221) isTrueEtaDalitz = kTRUE; + } + if (gamma1DalitzCand && gamma1MCLabel >=0 && gamma1MCLabel==gamma0MotherLabel){ + if (gamma1MotherLabel == -111) isTruePi0Dalitz = kTRUE; + if (gamma1MotherLabel == -221) isTrueEtaDalitz = kTRUE; + } + } + + + if(isTruePi0 || isTrueEta){// True Pion or Eta + fHistoTrueMotherInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt()); + if (fDoMesonQA > 0){ + if (isTruePi0){ + if ( Pi0Candidate->M() > 0.05 && Pi0Candidate->M() < 0.17){ + fHistoTruePi0PtY[fiCut]->Fill(Pi0Candidate->Pt(),Pi0Candidate->Rapidity()-((AliConversionCuts*)fCutArray->At(fiCut))->GetEtaShift()); + fHistoTruePi0PtAlpha[fiCut]->Fill(Pi0Candidate->Pt(),Pi0Candidate->GetAlpha()); + fHistoTruePi0PtOpenAngle[fiCut]->Fill(Pi0Candidate->Pt(),Pi0Candidate->GetOpeningAngle()); + } + } else if (isTrueEta){ + if ( Pi0Candidate->M() > 0.45 && Pi0Candidate->M() < 0.65){ + fHistoTrueEtaPtY[fiCut]->Fill(Pi0Candidate->Pt(),Pi0Candidate->Rapidity()-((AliConversionCuts*)fCutArray->At(fiCut))->GetEtaShift()); + fHistoTrueEtaPtAlpha[fiCut]->Fill(Pi0Candidate->Pt(),Pi0Candidate->GetAlpha()); + fHistoTrueEtaPtOpenAngle[fiCut]->Fill(Pi0Candidate->Pt(),Pi0Candidate->GetOpeningAngle()); + } + } + } + if(gamma0MotherLabel >= MCStack->GetNprimary()){ // Secondary Meson + Int_t secMotherLabel = ((TParticle*)MCStack->Particle(gamma1MotherLabel))->GetMother(0); + Float_t weightedSec= 1; + if(((AliConversionCuts*)fCutArray->At(fiCut))->IsParticleFromBGEvent(secMotherLabel, fMCStack, fInputEvent) && MCStack->Particle(secMotherLabel)->GetPdgCode()==310){ + weightedSec= ((AliConversionCuts*)fCutArray->At(fiCut))->GetWeightForMeson(fV0Reader->GetPeriodName(),secMotherLabel, fMCStack, fInputEvent)/2.; //invariant mass is additive thus the weight for the daughters has to be devide by two for the K0s at a certain pt + //cout << "MC input \t"<Pt()<<"\t"<Fill(Pi0Candidate->M(),Pi0Candidate->Pt(),weightedSec); + fCharMesonMCInfo = 2; + if (secMotherLabel >-1){ + if(MCStack->Particle(secMotherLabel)->GetPdgCode()==310){ + fCharMesonMCInfo = 4; + fHistoTrueSecondaryMotherFromK0sInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt(),weightedSec); + if (fDoMesonQA > 0)fHistoTrueK0sWithPi0DaughterMCPt[fiCut]->Fill(MCStack->Particle(secMotherLabel)->Pt()); + } + if(MCStack->Particle(secMotherLabel)->GetPdgCode()==221){ + fCharMesonMCInfo = 3; + fHistoTrueSecondaryMotherFromEtaInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt(),weightedSec); + if (fDoMesonQA > 0)fHistoTrueEtaWithPi0DaughterMCPt[fiCut]->Fill(MCStack->Particle(secMotherLabel)->Pt()); + } + if(MCStack->Particle(secMotherLabel)->GetPdgCode()==3122){ + fCharMesonMCInfo = 7; + fHistoTrueSecondaryMotherFromLambdaInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt(),weightedSec); + if (fDoMesonQA > 0)fHistoTrueLambdaWithPi0DaughterMCPt[fiCut]->Fill(MCStack->Particle(secMotherLabel)->Pt()); + } + } + } else { // Only primary pi0 for efficiency calculation + fCharMesonMCInfo = 6; + Float_t weighted= 1; + if(((AliConversionCuts*)fCutArray->At(fiCut))->IsParticleFromBGEvent(gamma1MotherLabel, fMCStack, fInputEvent)){ + if (((TParticle*)MCStack->Particle(gamma1MotherLabel))->Pt()>0.005){ + weighted= ((AliConversionCuts*)fCutArray->At(fiCut))->GetWeightForMeson(fV0Reader->GetPeriodName(),gamma1MotherLabel, fMCStack, fInputEvent); + // cout << "rec \t " <Fill(Pi0Candidate->M(),Pi0Candidate->Pt(),weighted); + fHistoTruePrimaryMotherW0WeightingInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt()); + fProfileTruePrimaryMotherWeightsInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt(),weighted); + + + if (fDoMesonQA > 0){ + if(isTruePi0){ // Only primary pi0 for resolution + fHistoTruePrimaryPi0MCPtResolPt[fiCut]->Fill(((TParticle*)MCStack->Particle(gamma1MotherLabel))->Pt(),(Pi0Candidate->Pt()-((TParticle*)MCStack->Particle(gamma1MotherLabel))->Pt())/((TParticle*)MCStack->Particle(gamma1MotherLabel))->Pt(),weighted); + } + if (isTrueEta){ // Only primary eta for resolution + fHistoTruePrimaryEtaMCPtResolPt[fiCut]->Fill(((TParticle*)MCStack->Particle(gamma1MotherLabel))->Pt(),(Pi0Candidate->Pt()-((TParticle*)MCStack->Particle(gamma1MotherLabel))->Pt())/((TParticle*)MCStack->Particle(gamma1MotherLabel))->Pt(),weighted); + } + } + } + } else if(!isTruePi0 && !isTrueEta){ // Background + if (fDoMesonQA > 0){ + if(gamma0MotherLabel>-1 && gamma1MotherLabel>-1){ // Both Tracks are Photons and have a mother but not Pi0 or Eta + fHistoTrueBckGGInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt()); + fCharMesonMCInfo = 1; + } else { // No photon or without mother + fHistoTrueBckContInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt()); + } + } + if( isTruePi0Dalitz || isTrueEtaDalitz ){ + // Dalitz + fCharMesonMCInfo = 5; + fHistoTrueMotherDalitzInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt()); + } else if (gamma0DalitzCand || gamma1DalitzCand){ + if (fDoMesonQA > 0)fHistoTrueBckContInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt()); + } + } + } + } +} +//______________________________________________________________________ +void AliAnalysisTaskGammaConvCalo::ProcessTrueMesonCandidatesAOD(AliAODConversionMother *Pi0Candidate, AliAODConversionPhoton *TrueGammaCandidate0, AliAODConversionPhoton *TrueGammaCandidate1) +{ + + // Process True Mesons + TClonesArray *AODMCTrackArray = dynamic_cast(fInputEvent->FindListObject(AliAODMCParticle::StdBranchName())); + Bool_t isTruePi0 = kFALSE; + Bool_t isTrueEta = kFALSE; + Bool_t isTruePi0Dalitz = kFALSE; + Bool_t isTrueEtaDalitz = kFALSE; + Bool_t gamma0DalitzCand = kFALSE; + Bool_t gamma1DalitzCand = kFALSE; + + AliAODMCParticle *positiveMC = static_cast(AODMCTrackArray->At(TrueGammaCandidate0->GetMCLabelPositive())); + AliAODMCParticle *negativeMC = static_cast(AODMCTrackArray->At(TrueGammaCandidate0->GetMCLabelNegative())); + + fCharMesonMCInfo = 0; + Int_t gamma0MCLabel = -1; + Int_t gamma0MotherLabel = -1; + if(!positiveMC||!negativeMC) + return; + + if(positiveMC->GetMother()>-1&&(negativeMC->GetMother() == positiveMC->GetMother())){ + gamma0MCLabel = positiveMC->GetMother(); + } + + if(gamma0MCLabel != -1){ // Gamma is Combinatorial; MC Particles don't belong to the same Mother + // Daughters Gamma 0 + AliAODMCParticle * gammaMC0 = static_cast(AODMCTrackArray->At(gamma0MCLabel)); + if(abs(negativeMC->GetPdgCode())==11 && abs(positiveMC->GetPdgCode())==11){ // Electrons ... + if(((positiveMC->GetMCProcessCode())) == 5 && ((negativeMC->GetMCProcessCode())) == 5){ // ... From Conversion ... + if(gammaMC0->GetPdgCode() == 22){ // ... with Gamma Mother + gamma0MotherLabel=gammaMC0->GetMother(); + } + } + if(gammaMC0->GetPdgCode() ==111){ // Dalitz candidate + gamma0DalitzCand = kTRUE; + gamma0MotherLabel=-111; + } + if(gammaMC0->GetPdgCode() ==221){ // Dalitz candidate + gamma0DalitzCand = kTRUE; + gamma0MotherLabel=-221; + } + } + } + positiveMC = static_cast(AODMCTrackArray->At(TrueGammaCandidate1->GetMCLabelPositive())); + negativeMC = static_cast(AODMCTrackArray->At(TrueGammaCandidate1->GetMCLabelNegative())); + + Int_t gamma1MCLabel = -1; + Int_t gamma1MotherLabel = -1; + if(!positiveMC||!negativeMC) + return; + + if(positiveMC->GetMother()>-1&&(negativeMC->GetMother() == positiveMC->GetMother())){ + gamma1MCLabel = positiveMC->GetMother(); + } + if(gamma1MCLabel != -1){ // Gamma is Combinatorial; MC Particles don't belong to the same Mother + // Daughters Gamma 1 + AliAODMCParticle * gammaMC1 = static_cast(AODMCTrackArray->At(gamma1MCLabel)); + if(abs(negativeMC->GetPdgCode())==11 && abs(positiveMC->GetPdgCode())==11){ // Electrons ... + if(((positiveMC->GetMCProcessCode())) == 5 && ((negativeMC->GetMCProcessCode())) == 5){ // ... From Conversion ... + if(gammaMC1->GetPdgCode() == 22){ // ... with Gamma Mother + gamma1MotherLabel=gammaMC1->GetMother(); + } + } + if(gammaMC1->GetPdgCode() ==111 ){ // Dalitz candidate + gamma1DalitzCand = kTRUE; + gamma1MotherLabel=-111; + } + if(gammaMC1->GetPdgCode() ==221){ // Dalitz candidate + gamma1DalitzCand = kTRUE; + gamma1MotherLabel=-221; + } + } + } + if(gamma0MotherLabel>=0 && gamma0MotherLabel==gamma1MotherLabel){ + if(static_cast(AODMCTrackArray->At(gamma1MotherLabel))->GetPdgCode() == 111){ + isTruePi0=kTRUE; + } + if(static_cast(AODMCTrackArray->At(gamma1MotherLabel))->GetPdgCode() == 221){ + isTrueEta=kTRUE; + } + } + + //Identify Dalitz candidate + if (gamma1DalitzCand || gamma0DalitzCand){ + if (gamma0DalitzCand && gamma0MCLabel >=0 && gamma0MCLabel==gamma1MotherLabel){ + if (gamma0MotherLabel == -111) isTruePi0Dalitz = kTRUE; + if (gamma0MotherLabel == -221) isTrueEtaDalitz = kTRUE; + } + if (gamma1DalitzCand && gamma1MCLabel >=0 && gamma1MCLabel==gamma0MotherLabel){ + if (gamma1MotherLabel == -111) isTruePi0Dalitz = kTRUE; + if (gamma1MotherLabel == -221) isTrueEtaDalitz = kTRUE; + } + } + + if(isTruePi0 || isTrueEta){// True Pion or Eta + fHistoTrueMotherInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt()); + if (fDoMesonQA > 0){ + if (isTruePi0){ + if ( Pi0Candidate->M() > 0.05 && Pi0Candidate->M() < 0.17){ + fHistoTruePi0PtY[fiCut]->Fill(Pi0Candidate->Pt(),Pi0Candidate->Rapidity()-((AliConversionCuts*)fCutArray->At(fiCut))->GetEtaShift()); + fHistoTruePi0PtAlpha[fiCut]->Fill(Pi0Candidate->Pt(),Pi0Candidate->GetAlpha()); + fHistoTruePi0PtOpenAngle[fiCut]->Fill(Pi0Candidate->Pt(),Pi0Candidate->GetOpeningAngle()); + } + } else if (isTrueEta){ + if ( Pi0Candidate->M() > 0.45 && Pi0Candidate->M() < 0.65){ + fHistoTrueEtaPtY[fiCut]->Fill(Pi0Candidate->Pt(),Pi0Candidate->Rapidity()-((AliConversionCuts*)fCutArray->At(fiCut))->GetEtaShift()); + fHistoTrueEtaPtAlpha[fiCut]->Fill(Pi0Candidate->Pt(),Pi0Candidate->GetAlpha()); + fHistoTrueEtaPtOpenAngle[fiCut]->Fill(Pi0Candidate->Pt(),Pi0Candidate->GetOpeningAngle()); + } + } + } + if(!(static_cast(AODMCTrackArray->At(gamma0MotherLabel))->IsPrimary())){ // Secondary Meson + Int_t secMotherLabel = static_cast(AODMCTrackArray->At(gamma1MotherLabel))->GetMother(); + Float_t weightedSec= 1; + if(((AliConversionCuts*)fCutArray->At(fiCut))->IsParticleFromBGEvent(secMotherLabel, 0x0, fInputEvent) && static_cast(AODMCTrackArray->At(secMotherLabel))->GetPdgCode()==310){ + weightedSec= ((AliConversionCuts*)fCutArray->At(fiCut))->GetWeightForMeson(fV0Reader->GetPeriodName(),secMotherLabel, 0x0, fInputEvent)/2.; //invariant mass is additive thus the weight for the daughters has to be devide by two for the K0s at a certain pt + //cout << "MC input \t"<Pt()<<"\t"<Fill(Pi0Candidate->M(),Pi0Candidate->Pt(),weightedSec); + fCharMesonMCInfo = 2; + if (secMotherLabel >-1){ + if(static_cast(AODMCTrackArray->At(secMotherLabel))->GetPdgCode()==310){ + fCharMesonMCInfo = 4; + fHistoTrueSecondaryMotherFromK0sInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt(),weightedSec); + if (fDoMesonQA > 0)fHistoTrueK0sWithPi0DaughterMCPt[fiCut]->Fill(static_cast(AODMCTrackArray->At(secMotherLabel))->Pt()); + } + if(static_cast(AODMCTrackArray->At(secMotherLabel))->GetPdgCode()==221){ + fCharMesonMCInfo = 3; + fHistoTrueSecondaryMotherFromEtaInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt(),weightedSec); + if (fDoMesonQA > 0)fHistoTrueEtaWithPi0DaughterMCPt[fiCut]->Fill(static_cast(AODMCTrackArray->At(secMotherLabel))->Pt()); + } + if(static_cast(AODMCTrackArray->At(secMotherLabel))->GetPdgCode()==3122){ + fCharMesonMCInfo = 7; + fHistoTrueSecondaryMotherFromLambdaInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt(),weightedSec); + if (fDoMesonQA > 0)fHistoTrueLambdaWithPi0DaughterMCPt[fiCut]->Fill(static_cast(AODMCTrackArray->At(secMotherLabel))->Pt()); + } + } + }else{ // Only primary pi0 for efficiency calculation + Float_t weighted= 1; + fCharMesonMCInfo = 6; + if(((AliConversionCuts*)fCutArray->At(fiCut))->IsParticleFromBGEvent(gamma1MotherLabel, 0x0, fInputEvent)){ + if (static_cast(AODMCTrackArray->At(gamma1MotherLabel))->Pt()>0.005){ + weighted= ((AliConversionCuts*)fCutArray->At(fiCut))->GetWeightForMeson(fV0Reader->GetPeriodName(),gamma1MotherLabel, 0x0, fInputEvent); + // cout << "rec \t " <Fill(Pi0Candidate->M(),Pi0Candidate->Pt(),weighted); + fHistoTruePrimaryMotherW0WeightingInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt()); + fProfileTruePrimaryMotherWeightsInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt(),weighted); + + if (fDoMesonQA > 0){ + if(isTruePi0){ // Only primary pi0 for resolution + fHistoTruePrimaryPi0MCPtResolPt[fiCut]->Fill(static_cast(AODMCTrackArray->At(gamma1MotherLabel))->Pt(), + (Pi0Candidate->Pt()-static_cast(AODMCTrackArray->At(gamma1MotherLabel))->Pt())/static_cast(AODMCTrackArray->At(gamma1MotherLabel))->Pt(),weighted); + + } + if (isTrueEta){ // Only primary eta for resolution + fHistoTruePrimaryEtaMCPtResolPt[fiCut]->Fill(static_cast(AODMCTrackArray->At(gamma1MotherLabel))->Pt(), + (Pi0Candidate->Pt()-static_cast(AODMCTrackArray->At(gamma1MotherLabel))->Pt())/static_cast(AODMCTrackArray->At(gamma1MotherLabel))->Pt(),weighted); + } + } + } + } else if(!isTruePi0 && !isTrueEta) { // Background + if (fDoMesonQA > 0){ + if(gamma0MotherLabel>-1 && gamma1MotherLabel>-1){ // Both Tracks are Photons and have a mother but not Pi0 or Eta + fHistoTrueBckGGInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt()); + fCharMesonMCInfo = 1; + } else { // No photon or without mother + fHistoTrueBckContInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt()); + } + } + if( isTruePi0Dalitz || isTrueEtaDalitz ){ + // Dalitz + fCharMesonMCInfo = 5; + fHistoTrueMotherDalitzInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt()); + } else if (gamma0DalitzCand || gamma1DalitzCand){ + if (fDoMesonQA > 0)fHistoTrueBckContInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt()); + } + } +} +//________________________________________________________________________ +void AliAnalysisTaskGammaConvCalo::CalculateBackground(){ + + Int_t zbin= fBGClusHandler[fiCut]->GetZBinIndex(fInputEvent->GetPrimaryVertex()->GetZ()); + Int_t mbin = 0; + + if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->UseTrackMultiplicity()){ + mbin = fBGClusHandler[fiCut]->GetMultiplicityBinIndex(fV0Reader->GetNumberOfPrimaryTracks()); + } else { + mbin = fBGClusHandler[fiCut]->GetMultiplicityBinIndex(fGammaCandidates->GetEntries()); + } + + AliGammaConversionAODBGHandler::GammaConversionVertex *bgEventVertex = NULL; + if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->UseTrackMultiplicity()){ + for(Int_t nEventsInBG=0;nEventsInBGGetNBGEvents();nEventsInBG++){ + AliGammaConversionAODVector *previousEventV0s = fBGClusHandler[fiCut]->GetBGGoodV0s(zbin,mbin,nEventsInBG); + if(fMoveParticleAccordingToVertex == kTRUE){ + bgEventVertex = fBGClusHandler[fiCut]->GetBGEventVertex(zbin,mbin,nEventsInBG); + } + + for(Int_t iCurrent=0;iCurrentGetEntries();iCurrent++){ + AliAODConversionPhoton currentEventGoodV0 = *(AliAODConversionPhoton*)(fGammaCandidates->At(iCurrent)); + for(UInt_t iPrevious=0;iPrevioussize();iPrevious++){ + AliAODConversionPhoton previousGoodV0 = (AliAODConversionPhoton)(*(previousEventV0s->at(iPrevious))); + if(fMoveParticleAccordingToVertex == kTRUE){ + MoveParticleAccordingToVertex(&previousGoodV0,bgEventVertex); + } + if(((AliConversionCuts*)fCutArray->At(fiCut))->GetInPlaneOutOfPlaneCut() != 0){ + RotateParticleAccordingToEP(&previousGoodV0,bgEventVertex->fEP,fEventPlaneAngle); + } + + AliAODConversionMother *backgroundCandidate = new AliAODConversionMother(¤tEventGoodV0,&previousGoodV0); + backgroundCandidate->CalculateDistanceOfClossetApproachToPrimVtx(fInputEvent->GetPrimaryVertex()); + if((((AliConversionMesonCuts*)fMesonCutArray->At(fiCut)) + ->MesonIsSelected(backgroundCandidate,kFALSE,((AliConversionCuts*)fCutArray->At(fiCut))->GetEtaShift()))){ + fHistoMotherBackInvMassPt[fiCut]->Fill(backgroundCandidate->M(),backgroundCandidate->Pt()); + Double_t sparesFill[4] = {backgroundCandidate->M(),backgroundCandidate->Pt(),(Double_t)zbin,(Double_t)mbin}; + fSparseMotherBackInvMassPtZM[fiCut]->Fill(sparesFill,1); + } + delete backgroundCandidate; + backgroundCandidate = 0x0; + } + } + } + } else { + for(Int_t nEventsInBG=0;nEventsInBG GetNBGEvents();nEventsInBG++){ + AliGammaConversionAODVector *previousEventV0s = fBGClusHandler[fiCut]->GetBGGoodV0s(zbin,mbin,nEventsInBG); + if(previousEventV0s){ + if(fMoveParticleAccordingToVertex == kTRUE){ + bgEventVertex = fBGClusHandler[fiCut]->GetBGEventVertex(zbin,mbin,nEventsInBG); + } + for(Int_t iCurrent=0;iCurrentGetEntries();iCurrent++){ + AliAODConversionPhoton currentEventGoodV0 = *(AliAODConversionPhoton*)(fGammaCandidates->At(iCurrent)); + for(UInt_t iPrevious=0;iPrevioussize();iPrevious++){ + + AliAODConversionPhoton previousGoodV0 = (AliAODConversionPhoton)(*(previousEventV0s->at(iPrevious))); + + if(fMoveParticleAccordingToVertex == kTRUE){ + MoveParticleAccordingToVertex(&previousGoodV0,bgEventVertex); + } + if(((AliConversionCuts*)fCutArray->At(fiCut))->GetInPlaneOutOfPlaneCut() != 0){ + RotateParticleAccordingToEP(&previousGoodV0,bgEventVertex->fEP,fEventPlaneAngle); + } + + + AliAODConversionMother *backgroundCandidate = new AliAODConversionMother(¤tEventGoodV0,&previousGoodV0); + backgroundCandidate->CalculateDistanceOfClossetApproachToPrimVtx(fInputEvent->GetPrimaryVertex()); + if((((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->MesonIsSelected(backgroundCandidate,kFALSE,((AliConversionCuts*)fCutArray->At(fiCut))->GetEtaShift()))){ + fHistoMotherBackInvMassPt[fiCut]->Fill(backgroundCandidate->M(),backgroundCandidate->Pt()); + Double_t sparesFill[4] = {backgroundCandidate->M(),backgroundCandidate->Pt(),(Double_t)zbin,(Double_t)mbin}; + fSparseMotherBackInvMassPtZM[fiCut]->Fill(sparesFill,1); + } + delete backgroundCandidate; + backgroundCandidate = 0x0; + } + } + } + } + } +} + +//________________________________________________________________________ +void AliAnalysisTaskGammaConvCalo::CalculateBackgroundRP(){ + + Int_t zbin= fBGHandlerRP[fiCut]->GetZBinIndex(fInputEvent->GetPrimaryVertex()->GetZ()); + Int_t mbin = 0; + if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->UseTrackMultiplicity()){ + mbin = fBGHandlerRP[fiCut]->GetMultiplicityBinIndex(fV0Reader->GetNumberOfPrimaryTracks()); + } else { + mbin = fBGHandlerRP[fiCut]->GetMultiplicityBinIndex(fGammaCandidates->GetEntries()); + } + + + //Rotation Method + if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->UseRotationMethod()){ + // Correct for the number of rotations + // BG is for rotation the same, except for factor NRotations + Double_t weight=1./Double_t(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->GetNumberOfBGEvents()); + + for(Int_t firstGammaIndex=0;firstGammaIndexGetEntries();firstGammaIndex++){ + + AliAODConversionPhoton *gamma0=dynamic_cast(fGammaCandidates->At(firstGammaIndex)); + if (gamma0==NULL) continue; + for(Int_t secondGammaIndex=firstGammaIndex+1;secondGammaIndexGetEntries();secondGammaIndex++){ + AliAODConversionPhoton *gamma1=dynamic_cast(fGammaCandidates->At(secondGammaIndex)); + if (gamma1 == NULL) continue; + if(!((AliConversionCuts*)fCutArray->At(fiCut))->PhotonIsSelected(gamma1,fInputEvent))continue; + for(Int_t nRandom=0;nRandom<((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->GetNumberOfBGEvents();nRandom++){ + + RotateParticle(gamma1); + AliAODConversionMother backgroundCandidate(gamma0,gamma1); + backgroundCandidate.CalculateDistanceOfClossetApproachToPrimVtx(fInputEvent->GetPrimaryVertex()); + if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut)) + ->MesonIsSelected(&backgroundCandidate,kFALSE,((AliConversionCuts*)fCutArray->At(fiCut))->GetEtaShift())){ + fHistoMotherBackInvMassPt[fiCut]->Fill(backgroundCandidate.M(),backgroundCandidate.Pt()); + Double_t sparesFill[4] = {backgroundCandidate.M(),backgroundCandidate.Pt(),(Double_t)zbin,(Double_t)mbin}; + fSparseMotherBackInvMassPtZM[fiCut]->Fill(sparesFill,weight); + } + } + } + } + } + else{ + // Do Event Mixing + for(Int_t nEventsInBG=0;nEventsInBG GetNBGEvents(fGammaCandidates,fInputEvent);nEventsInBG++){ + + AliGammaConversionPhotonVector *previousEventGammas = fBGHandlerRP[fiCut]->GetBGGoodGammas(fGammaCandidates,fInputEvent,nEventsInBG); + + if(previousEventGammas){ + // test weighted background + Double_t weight=1.0; + // Correct for the number of eventmixing: + // N gammas -> (N-1) + (N-2) +(N-3) ...+ (N-(N-1)) using sum formula sum(i)=N*(N-1)/2 -> N*(N-1)/2 + // real combinations (since you cannot combine a photon with its own) + // but BG leads to N_{a}*N_{b} combinations + weight*=0.5*(Double_t(fGammaCandidates->GetEntries()-1))/Double_t(previousEventGammas->size()); + + for(Int_t iCurrent=0;iCurrentGetEntries();iCurrent++){ + AliAODConversionPhoton *gamma0 = (AliAODConversionPhoton*)(fGammaCandidates->At(iCurrent)); + for(UInt_t iPrevious=0;iPrevioussize();iPrevious++){ + + AliAODConversionPhoton *gamma1 = (AliAODConversionPhoton*)(previousEventGammas->at(iPrevious)); + + AliAODConversionMother backgroundCandidate(gamma0,gamma1); + backgroundCandidate.CalculateDistanceOfClossetApproachToPrimVtx(fInputEvent->GetPrimaryVertex()); + if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut)) + ->MesonIsSelected(&backgroundCandidate,kFALSE,((AliConversionCuts*)fCutArray->At(fiCut))->GetEtaShift())){ + fHistoMotherBackInvMassPt[fiCut]->Fill(backgroundCandidate.M(),backgroundCandidate.Pt()); + Double_t sparesFill[4] = {backgroundCandidate.M(),backgroundCandidate.Pt(),(Double_t)zbin,(Double_t)mbin}; + fSparseMotherBackInvMassPtZM[fiCut]->Fill(sparesFill,weight); + } + } + } + } + } + } +} +//________________________________________________________________________ +void AliAnalysisTaskGammaConvCalo::RotateParticle(AliAODConversionPhoton *gamma){ + Int_t fNDegreesPMBackground= ((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->NDegreesRotation(); + Double_t nRadiansPM = fNDegreesPMBackground*TMath::Pi()/180; + Double_t rotationValue = fRandom.Rndm()*2*nRadiansPM + TMath::Pi()-nRadiansPM; + gamma->RotateZ(rotationValue); +} + +//________________________________________________________________________ +void AliAnalysisTaskGammaConvCalo::RotateParticleAccordingToEP(AliAODConversionPhoton *gamma, Double_t previousEventEP, Double_t thisEventEP){ + + previousEventEP=previousEventEP+TMath::Pi(); + thisEventEP=thisEventEP+TMath::Pi(); + Double_t rotationValue= thisEventEP-previousEventEP; + gamma->RotateZ(rotationValue); +} + +//________________________________________________________________________ +void AliAnalysisTaskGammaConvCalo::MoveParticleAccordingToVertex(AliAODConversionPhoton* particle,const AliGammaConversionAODBGHandler::GammaConversionVertex *vertex){ + //see header file for documentation + + Double_t dx = vertex->fX - fInputEvent->GetPrimaryVertex()->GetX(); + Double_t dy = vertex->fY - fInputEvent->GetPrimaryVertex()->GetY(); + Double_t dz = vertex->fZ - fInputEvent->GetPrimaryVertex()->GetZ(); + + Double_t movedPlace[3] = {particle->GetConversionX() - dx,particle->GetConversionY() - dy,particle->GetConversionZ() - dz}; + particle->SetConversionPoint(movedPlace); +} +//________________________________________________________________________ +void AliAnalysisTaskGammaConvCalo::UpdateEventByEventData(){ + //see header file for documentation + if(fGammaCandidates->GetEntries() >0 ){ + if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->UseTrackMultiplicity()){ + fBGHandler[fiCut]->AddEvent(fGammaCandidates,fInputEvent->GetPrimaryVertex()->GetX(),fInputEvent->GetPrimaryVertex()->GetY(),fInputEvent->GetPrimaryVertex()->GetZ(),fV0Reader->GetNumberOfPrimaryTracks(),fEventPlaneAngle); + fBGClusHandler[fiCut]->AddEvent(fClusterCandidates,fInputEvent->GetPrimaryVertex()->GetX(),fInputEvent->GetPrimaryVertex()->GetY(),fInputEvent->GetPrimaryVertex()->GetZ(),fV0Reader->GetNumberOfPrimaryTracks(),fEventPlaneAngle); + } else { // means we use #V0s for multiplicity + fBGHandler[fiCut]->AddEvent(fGammaCandidates,fInputEvent->GetPrimaryVertex()->GetX(),fInputEvent->GetPrimaryVertex()->GetY(),fInputEvent->GetPrimaryVertex()->GetZ(),fGammaCandidates->GetEntries(),fEventPlaneAngle); + fBGClusHandler[fiCut]->AddEvent(fClusterCandidates,fInputEvent->GetPrimaryVertex()->GetX(),fInputEvent->GetPrimaryVertex()->GetY(),fInputEvent->GetPrimaryVertex()->GetZ(),fGammaCandidates->GetEntries(),fEventPlaneAngle); + } + } +} + + +//________________________________________________________________________ +void AliAnalysisTaskGammaConvCalo::FillPhotonCombinatorialBackgroundHist(AliAODConversionPhoton *TruePhotonCandidate, Int_t pdgCode[]) +{ + // Combinatorial Bck = 0 ee, 1 ep,i 2 ek, 3 ep, 4 emu, 5 pipi, 6 pik, 7 pip, 8 pimu, 9 kk, 10 kp, 11 kmu, 12 pp, 13 pmu, 14 mumu, 15 Rest + if(pdgCode[0]==11 && pdgCode[1]==11){if(fIsFromMBHeader)fHistoCombinatorialPt[fiCut]->Fill(TruePhotonCandidate->Pt(),0);} + else if( (pdgCode[0]==11 && pdgCode[1]==211) || (pdgCode[0]==211 && pdgCode[1]==11) ) + {if(fIsFromMBHeader)fHistoCombinatorialPt[fiCut]->Fill(TruePhotonCandidate->Pt(),1);} + else if( (pdgCode[0]==11 && pdgCode[1]==321) || (pdgCode[0]==321 && pdgCode[1]==11) ) + {if(fIsFromMBHeader)fHistoCombinatorialPt[fiCut]->Fill(TruePhotonCandidate->Pt(),2);} + else if( (pdgCode[0]==11 && pdgCode[1]==2212) || (pdgCode[0]==2212 && pdgCode[1]==11) ) + {if(fIsFromMBHeader)fHistoCombinatorialPt[fiCut]->Fill(TruePhotonCandidate->Pt(),3);} + else if( (pdgCode[0]==11 && pdgCode[1]==13) || (pdgCode[0]==13 && pdgCode[1]==11) ) + {if(fIsFromMBHeader)fHistoCombinatorialPt[fiCut]->Fill(TruePhotonCandidate->Pt(),4);} + else if( pdgCode[0]==211 && pdgCode[1]==211 ){if(fIsFromMBHeader)fHistoCombinatorialPt[fiCut]->Fill(TruePhotonCandidate->Pt(),5);} + else if( (pdgCode[0]==211 && pdgCode[1]==321) || (pdgCode[0]==321 && pdgCode[1]==211) ) + {if(fIsFromMBHeader)fHistoCombinatorialPt[fiCut]->Fill(TruePhotonCandidate->Pt(),6);} + else if( (pdgCode[0]==211 && pdgCode[1]==2212) || (pdgCode[0]==2212 && pdgCode[1]==211) ) + {if(fIsFromMBHeader)fHistoCombinatorialPt[fiCut]->Fill(TruePhotonCandidate->Pt(),7);} + else if( (pdgCode[0]==211 && pdgCode[1]==13) || (pdgCode[0]==13 && pdgCode[1]==211) ) + {if(fIsFromMBHeader)fHistoCombinatorialPt[fiCut]->Fill(TruePhotonCandidate->Pt(),8);} + else if( pdgCode[0]==321 && pdgCode[1]==321 ){if(fIsFromMBHeader)fHistoCombinatorialPt[fiCut]->Fill(TruePhotonCandidate->Pt(),9);} + else if( (pdgCode[0]==321 && pdgCode[1]==2212) || (pdgCode[0]==2212 && pdgCode[1]==321) ) + {if(fIsFromMBHeader)fHistoCombinatorialPt[fiCut]->Fill(TruePhotonCandidate->Pt(),10);} + else if( (pdgCode[0]==321 && pdgCode[1]==13) || (pdgCode[0]==13 && pdgCode[1]==321) ) + {if(fIsFromMBHeader)fHistoCombinatorialPt[fiCut]->Fill(TruePhotonCandidate->Pt(),11);} + else if( pdgCode[0]==2212 && pdgCode[1]==2212 ){if(fIsFromMBHeader)fHistoCombinatorialPt[fiCut]->Fill(TruePhotonCandidate->Pt(),12);} + else if( (pdgCode[0]==2212 && pdgCode[1]==13) || (pdgCode[0]==13 && pdgCode[1]==2212) ) + {if(fIsFromMBHeader)fHistoCombinatorialPt[fiCut]->Fill(TruePhotonCandidate->Pt(),13);} + else if( pdgCode[0]==13 && pdgCode[1]==13 ){if(fIsFromMBHeader)fHistoCombinatorialPt[fiCut]->Fill(TruePhotonCandidate->Pt(),14);} + else {if(fIsFromMBHeader)fHistoCombinatorialPt[fiCut]->Fill(TruePhotonCandidate->Pt(),15);} +} +//________________________________________________________________________ +void AliAnalysisTaskGammaConvCalo::RelabelAODPhotonCandidates(Bool_t mode){ + + // Relabeling For AOD Event + // ESDiD -> AODiD + // MCLabel -> AODMCLabel + + if(mode){ + fMCStackPos = new Int_t[fReaderGammas->GetEntries()]; + fMCStackNeg = new Int_t[fReaderGammas->GetEntries()]; + fESDArrayPos = new Int_t[fReaderGammas->GetEntries()]; + fESDArrayNeg = new Int_t[fReaderGammas->GetEntries()]; + } + + for(Int_t iGamma = 0;iGammaGetEntries();iGamma++){ + AliAODConversionPhoton* PhotonCandidate = (AliAODConversionPhoton*) fReaderGammas->At(iGamma); + if(!PhotonCandidate) continue; + if(!mode){// Back to ESD Labels + PhotonCandidate->SetMCLabelPositive(fMCStackPos[iGamma]); + PhotonCandidate->SetMCLabelNegative(fMCStackNeg[iGamma]); + PhotonCandidate->SetLabelPositive(fESDArrayPos[iGamma]); + PhotonCandidate->SetLabelNegative(fESDArrayNeg[iGamma]); + continue; + } + fMCStackPos[iGamma] = PhotonCandidate->GetMCLabelPositive(); + fMCStackNeg[iGamma] = PhotonCandidate->GetMCLabelNegative(); + fESDArrayPos[iGamma] = PhotonCandidate->GetTrackLabelPositive(); + fESDArrayNeg[iGamma] = PhotonCandidate->GetTrackLabelNegative(); + + Bool_t AODLabelPos = kFALSE; + Bool_t AODLabelNeg = kFALSE; + + for(Int_t i = 0; iGetNumberOfTracks();i++){ + AliAODTrack *tempDaughter = static_cast(fInputEvent->GetTrack(i)); + if(!AODLabelPos){ + if( tempDaughter->GetID() == PhotonCandidate->GetTrackLabelPositive() ){ + PhotonCandidate->SetMCLabelPositive(abs(tempDaughter->GetLabel())); + PhotonCandidate->SetLabelPositive(i); + AODLabelPos = kTRUE; + } + } + if(!AODLabelNeg){ + if( tempDaughter->GetID() == PhotonCandidate->GetTrackLabelNegative()){ + PhotonCandidate->SetMCLabelNegative(abs(tempDaughter->GetLabel())); + PhotonCandidate->SetLabelNegative(i); + AODLabelNeg = kTRUE; + } + } + if(AODLabelNeg && AODLabelPos){ + break; + } + } + if(!AODLabelPos || !AODLabelNeg){ + cout<<"WARNING!!! AOD TRACKS NOT FOUND FOR"<GetXaxis(); + Int_t bins = axisafter->GetNbins(); + Double_t from = axisafter->GetXmin(); + Double_t to = axisafter->GetXmax(); + Double_t *newbins = new Double_t[bins+1]; + newbins[0] = from; + Double_t factor = TMath::Power(to/from, 1./bins); + for(Int_t i=1; i<=bins; ++i) newbins[i] = factor * newbins[i-1]; + axisafter->Set(bins, newbins); + delete [] newbins; +} + +//________________________________________________________________________ +void AliAnalysisTaskGammaConvCalo::Terminate(const Option_t *) +{ + + //fOutputContainer->Print(); // Will crash on GRID +} + +//________________________________________________________________________ +Int_t AliAnalysisTaskGammaConvCalo::GetSourceClassification(Int_t daughter, Int_t pdgCode){ + + if (daughter == 111) { + if (abs(pdgCode) == 310) return 1; // k0s + else if (abs(pdgCode) == 3122) return 2; // Lambda + else if (abs(pdgCode) == 130) return 3; // K0L + else if (abs(pdgCode) == 2212) return 4; // proton + else if (abs(pdgCode) == 2112) return 5; // neutron + else if (abs(pdgCode) == 211) return 6; // pion + else if (abs(pdgCode) == 321) return 7; // kaon + else if (abs(pdgCode) == 113 || abs(pdgCode) == 213 ) return 8; // rho 0,+,- + else if (abs(pdgCode) == 3222 || abs(pdgCode) == 3212 || abs(pdgCode) == 3112 ) return 9; // Sigma + else if (abs(pdgCode) == 2224 || abs(pdgCode) == 2214 || abs(pdgCode) == 2114 || abs(pdgCode) == 1114 ) return 10; // Delta + else if (abs(pdgCode) == 313 || abs(pdgCode) == 323 ) return 11; // K* + else return 15; + } + return 15; + +} + +// //________________________________________________________________________ +// Double_t AliAnalysisTaskGammaConvCalo::GetMaxCellEnergy(const AliVCluster *cluster, Short_t &id) const +// { +// // Get maximum energy of attached cell. +// +// id = -1; +// Double_t maxe = 0; +// Int_t ncells = cluster->GetNCells(); +// if (fEsdCells) { +// for (Int_t i=0; iGetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i))); +// if (e>maxe) { +// maxe = e; +// id = cluster->GetCellAbsId(i); +// } +// } +// } else { +// for (Int_t i=0; iGetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i))); +// if (e>maxe) +// maxe = e; +// id = cluster->GetCellAbsId(i); +// } +// } +// return maxe; +// } diff --git a/PWGGA/GammaConv/AliAnalysisTaskGammaConvCalo.h b/PWGGA/GammaConv/AliAnalysisTaskGammaConvCalo.h new file mode 100644 index 00000000000..5e4e01786cf --- /dev/null +++ b/PWGGA/GammaConv/AliAnalysisTaskGammaConvCalo.h @@ -0,0 +1,300 @@ +#ifndef ALIANLYSISTASKGAMMACONVCALO_cxx +#define ALIANLYSISTASKGAMMACONVCALO_cxx + +#include "AliAnalysisTaskSE.h" +#include "AliESDtrack.h" +#include "AliV0ReaderV1.h" +#include "AliKFConversionPhoton.h" +#include "AliGammaConversionAODBGHandler.h" +#include "AliConversionAODBGHandlerRP.h" +#include "AliCaloPhotonCuts.h" +#include "AliConversionMesonCuts.h" +#include "AliAnalysisManager.h" +#include "TProfile2D.h" +#include "TH3.h" +#include "TH3F.h" + +class AliAnalysisTaskGammaConvCalo : public AliAnalysisTaskSE { + public: + + AliAnalysisTaskGammaConvCalo(); + AliAnalysisTaskGammaConvCalo(const char *name); + virtual ~AliAnalysisTaskGammaConvCalo(); + + virtual void UserCreateOutputObjects(); + virtual Bool_t Notify(); + virtual void UserExec(Option_t *); + virtual void Terminate(const Option_t*); + void InitBack(); + + void SetIsHeavyIon(Int_t flag){ + fIsHeavyIon = flag; + } + + // base functions for selecting photon and meson candidates in reconstructed data + void ProcessClusters(); + void ProcessPhotonCandidates(); + void PhotonTagging(); + void CalculatePi0Candidates(); + + // MC functions + void SetIsMC(Bool_t isMC){fIsMC=isMC;} + void ProcessMCParticles(); + void ProcessAODMCParticles(); + void RelabelAODPhotonCandidates(Bool_t mode); + void ProcessTruePhotonCandidates( AliAODConversionPhoton* TruePhotonCandidate); + void ProcessTrueClusterCandidates( AliAODConversionPhoton* TruePhotonCandidate); + void ProcessTruePhotonCandidatesAOD( AliAODConversionPhoton* TruePhotonCandidate); + void ProcessTrueMesonCandidates( AliAODConversionMother *Pi0Candidate, AliAODConversionPhoton *TrueGammaCandidate0, AliAODConversionPhoton *TrueGammaCandidate1); + void ProcessTrueMesonCandidatesAOD(AliAODConversionMother *Pi0Candidate, AliAODConversionPhoton *TrueGammaCandidate0, AliAODConversionPhoton *TrueGammaCandidate1); + + // switches for additional analysis streams or outputs + void SetDoMesonAnalysis(Bool_t flag){fDoMesonAnalysis = flag;} + void SetDoMesonQA(Int_t flag){fDoMesonQA = flag;} + void SetDoPhotonQA(Int_t flag){fDoPhotonQA = flag;} + + // Setting the cut lists for the conversion photons + void SetConversionCutList(Int_t nCuts, TList *CutArray){ + fnCuts = nCuts; + fCutArray = CutArray; + } + + // Setting the cut lists for the calo photons + void SetCaloCutList(Int_t nCuts, TList *CutArray){ + fnCuts = nCuts; + fClusterCutArray = CutArray; + } + + // Setting the cut lists for the meson + void SetMesonCutList(Int_t nCuts, TList *CutArray){ + fnCuts = nCuts; + fMesonCutArray = CutArray; + } + + // emcal functions + Double_t GetMaxCellEnergy(const AliVCluster *c) const { Short_t id=-1; return GetMaxCellEnergy(c,id); } + Double_t GetMaxCellEnergy(const AliVCluster *c, Short_t &id) const; + + // BG HandlerSettings + void CalculateBackground(); + void CalculateBackgroundRP(); + void RotateParticle(AliAODConversionPhoton *gamma); + void RotateParticleAccordingToEP(AliAODConversionPhoton *gamma, Double_t previousEventEP, Double_t thisEventEP); + void SetMoveParticleAccordingToVertex(Bool_t flag){fMoveParticleAccordingToVertex = flag;} + void FillPhotonCombinatorialBackgroundHist(AliAODConversionPhoton *TruePhotonCandidate, Int_t pdgCode[]); + void MoveParticleAccordingToVertex(AliAODConversionPhoton* particle,const AliGammaConversionAODBGHandler::GammaConversionVertex *vertex); + void UpdateEventByEventData(); + + // Additional functions for convenience + void SetLogBinningXTH2(TH2* histoRebin); + Int_t GetSourceClassification(Int_t daughter, Int_t pdgCode); + + protected: + AliV0ReaderV1 *fV0Reader; // basic photon Selection Task + AliGammaConversionAODBGHandler **fBGHandler; // BG handler for Conversion + AliConversionAODBGHandlerRP **fBGHandlerRP; // BG handler for Conversion (possibility to mix with respect to RP) + AliGammaConversionAODBGHandler **fBGClusHandler; // BG handler for Cluster + AliConversionAODBGHandlerRP **fBGClusHandlerRP; // BG handler for Cluster (possibility to mix with respect to RP) + AliVEvent *fInputEvent; // current event + AliMCEvent *fMCEvent; // corresponding MC event + AliStack *fMCStack; // stack belonging to MC event +// AliESDEvent *fEsdEv; //!pointer to input esd event +// AliAODEvent *fAodEv; //!pointer to input aod event +// TObjArray *fEsdClusters; //!pointer to esd clusters +// AliESDCaloCells *fEsdCells; //!pointer to esd cells +// TObjArray *fAodClusters; //!pointer to aod clusters +// AliAODCaloCells *fAodCells; //!pointer to aod cells + TList **fCutFolder; // Array of lists for containers belonging to cut + TList **fESDList; // Array of lists with histograms with reconstructed properties + TList **fBackList; // Array of lists with BG THnSparseF + TList **fMotherList; // Array of lists with Signal THnSparseF + TList **fPhotonDCAList; // Array of lists with photon dca trees + TList **fMesonDCAList; // Array of lists with meson dca trees + TList **fTrueList; // Array of lists with histograms with MC validated reconstructed properties + TList **fMCList; // Array of lists with histograms with pure MC information + TList **fHeaderNameList; // Array of lists with header names for MC header selection + TList **fTagOutputList; //!Array of lists of output histograms for tagged photons + TList *fOutputContainer; // Output container + TClonesArray *fReaderGammas; // Array with conversion photons selected by V0Reader Cut + TList *fGammaCandidates; // current list of photon candidates + TList *fClusterCandidates; //! current list of cluster candidates + TList *fCutArray; // List with Conversion Cuts + AliConversionCuts *fConversionCuts; // ConversionCutObject + TList *fClusterCutArray; // List with Cluster Cuts + AliCaloPhotonCuts *fCaloPhotonCuts; // CaloPhotonCutObject + TList *fMesonCutArray; // List with Meson Cuts + AliConversionMesonCuts *fMesonCuts; // MesonCutObject + + //histograms for Conversions reconstructed quantities + TH1F **fHistoConvGammaPt; //! histogram conversion photon pT + TH1F **fHistoConvGammaR; //! histogram conversion photon R + TH1F **fHistoConvGammaEta; //! histogram conversion photon Eta + TTree **fTreeConvGammaPtDcazCat; //! tree with dca for conversions + Float_t fPtGamma; //! pt of conversion for tree + Float_t fDCAzPhoton; //! dcaz of conversion for tree + Float_t fRConvPhoton; //! R of conversion for tree + Float_t fEtaPhoton; //! eta of conversion for tree + UChar_t fCharCatPhoton; //! category of conversion for tree + UChar_t fCharPhotonMCInfo; //! MC info of conversion for tree + // 0: garbage, + // 1: background + // 2: secondary photon not from eta or k0s, + // 3: secondary photon from eta, + // 4: secondary photon from k0s, + // 5: dalitz + // 6: primary gamma + //histograms for mesons reconstructed quantities + TH2F **fHistoMotherInvMassPt; //! array of histogram with signal + BG for same event photon pairs, inv Mass, pt + THnSparseF **fSparseMotherInvMassPtZM; //! array of THnSparseF with signal + BG for same event photon pairs, inv Mass, pt + TH2F **fHistoMotherBackInvMassPt; //! array of histogram with BG for mixed event photon pairs, inv Mass, pt + THnSparseF **fSparseMotherBackInvMassPtZM; //! array of THnSparseF with BG for same event photon pairs, inv Mass, pt + TH2F **fHistoMotherInvMassEalpha; //! array of histograms with alpha cut of 0.1 for inv mass vs pt + TH2F **fHistoMotherPi0PtY; //! array of histograms with invariant mass cut of 0.05 && pi0cand->M() < 0.17, pt, Y + TH2F **fHistoMotherEtaPtY; //! array of histograms with invariant mass cut of 0.45 && pi0cand->M() < 0.65, pt, Y + TH2F **fHistoMotherPi0PtAlpha; //! array of histograms with invariant mass cut of 0.05 && pi0cand->M() < 0.17, pt, alpha + TH2F **fHistoMotherEtaPtAlpha; //! array of histograms with invariant mass cut of 0.45 && pi0cand->M() < 0.65, pt, alpha + TH2F **fHistoMotherPi0PtOpenAngle; //! array of histograms with invariant mass cut of 0.05 && pi0cand->M() < 0.17, pt, openAngle + TH2F **fHistoMotherEtaPtOpenAngle; //! array of histograms with invariant mass cut of 0.45 && pi0cand->M() < 0.65, pt, openAngle + TTree **fTreeMesonsInvMassPtDcazMinDcazMaxFlag; //! array of trees with dca information for mesons + Float_t fInvMass; // inv mass for meson tree + Float_t fPt; // pt for meson tree + Float_t fDCAzGammaMin; // dcaz for meson tree gamma 1 + Float_t fDCAzGammaMax; // dcaz for meson tree gamma 2 + UChar_t fCharFlag; // category of meson for tree + UChar_t fCharMesonMCInfo; // MC information meson for tree + // 0: garbage, + // 1: background + // 2: secondary meson not from eta or k0s, + // 3: secondary meson from eta, + // 4: secondary meson from k0s, + // 5: dalitz + // 6: primary meson gamma-gamma-channel + + // histograms for rec photons tagged by Calo + TH1F **fHistoConvGammaUntagged; //! array of histo for untagged photon candidates vs pt + TH1F **fHistoConvGammaTagged; //! array of histo for tagged photon candidates vs pt + TH1F **fHistoConvGammaPi0Tagged; //! array of histo for tagged photon candidates vs pt + TH1F **fHistoConvGammaEtaTagged; //! array of histo for tagged photon candidates vs pt + TH2F **fHistoPhotonPairAll; //! array of histo for pairs + TH2F **fHistoPhotonPairAllGam; //! array of histo for pairs vs. pt of converted photon + // histograms for rec photon clusters + TH1F ** fHistoClusGammaPt; //! array of histos with cluster, pt + + //histograms for pure MC quantities + TH1I **fHistoMCHeaders; //! array of histos for header names + TH1F **fHistoMCAllGammaPt; //! array of histos with all gamma, pT + TH1F **fHistoMCDecayGammaPi0Pt; //! array of histos with decay gamma from pi0, pT + TH1F **fHistoMCDecayGammaRhoPt; //! array of histos with decay gamma from rho, pT + TH1F **fHistoMCDecayGammaEtaPt; //! array of histos with decay gamma from eta, pT + TH1F **fHistoMCDecayGammaOmegaPt; //! array of histos with decay gamma from omega, pT + TH1F **fHistoMCDecayGammaEtapPt; //! array of histos with decay gamma from eta', pT + TH1F **fHistoMCDecayGammaPhiPt; //! array of histos with decay gamma from phi, pT + TH1F **fHistoMCDecayGammaSigmaPt; //! array of histos with decay gamma from Sigma0, pT + TH1F **fHistoMCConvGammaPt; //! array of histos with converted gamma, pT + TH1F **fHistoMCConvGammaR; //! array of histos with converted gamma, R + TH1F **fHistoMCConvGammaEta; //! array of histos with converted gamma, Eta + TH1F **fHistoMCPi0Pt; //! array of histos with weighted pi0, pT + TH1F **fHistoMCPi0WOWeightPt; //! array of histos with unweighted pi0, pT + TH1F **fHistoMCEtaPt; //! array of histos with weighted eta, pT + TH1F **fHistoMCEtaWOWeightPt; //! array of histos with unweighted eta, pT + TH1F **fHistoMCPi0InAccPt; //! array of histos with weighted pi0 in acceptance, pT + TH1F **fHistoMCEtaInAccPt; //! array of histos with weighted eta in acceptance, pT + TH2F **fHistoMCPi0PtY; //! array of histos with weighted pi0, pT, Y + TH2F **fHistoMCEtaPtY; //! array of histos with weighted eta, pT, Y + TH1F **fHistoMCK0sPt; //! array of histos with weighted K0s, pT + TH1F **fHistoMCK0sWOWeightPt; //! array of histos with unweighted K0s, pT + TH2F **fHistoMCK0sPtY; //! array of histos with weighted K0s, pT, Y + TH2F **fHistoMCSecPi0PtvsSource; //! array of histos with secondary pi0, pT, source + TH1F **fHistoMCSecPi0Source; //! array of histos with secondary pi0, source + TH1F **fHistoMCSecEtaPt; //! array of histos with secondary eta, pT + TH1F **fHistoMCSecEtaSource; //! array of histos with secondary eta, source + // MC validated reconstructed quantities mesons + TH2F **fHistoTrueMotherInvMassPt; //! array of histos with validated mothers, invMass, pt + TH2F **fHistoTruePrimaryMotherInvMassPt; //! array of histos with validated weighted primary mothers, invMass, pt + TH2F **fHistoTruePrimaryMotherW0WeightingInvMassPt; //! array of histos with validated unweighted primary mothers, invMass, pt + TProfile2D **fProfileTruePrimaryMotherWeightsInvMassPt; //! array of profiles with weights for validated primary mothers, invMass, pt + TH2F **fHistoTruePrimaryPi0MCPtResolPt; //! array of histos with validated weighted primary pi0, MCpt, resol pt + TH2F **fHistoTruePrimaryEtaMCPtResolPt; //! array of histos with validated weighted primary eta, MCpt, resol pt + TH2F **fHistoTrueSecondaryMotherInvMassPt; //! array of histos with validated secondary mothers, invMass, pt + TH2F **fHistoTrueSecondaryMotherFromK0sInvMassPt; //! array of histos with validated secondary mothers from K0s, invMass, pt + TH1F **fHistoTrueK0sWithPi0DaughterMCPt; //! array of histos with K0s with reconstructed pi0 as daughter, pt + TH2F **fHistoTrueSecondaryMotherFromEtaInvMassPt; //! array of histos with validated secondary mothers from eta, invMass, pt + TH1F **fHistoTrueEtaWithPi0DaughterMCPt; //! array of histos with eta with reconstructed pi0 as daughter, pt + TH2F **fHistoTrueSecondaryMotherFromLambdaInvMassPt; //! array of histos with validated secondary mothers from Lambda, invMass, pt + TH1F **fHistoTrueLambdaWithPi0DaughterMCPt; //! array of histos with lambda with reconstructed pi0 as daughter, pt + TH2F **fHistoTrueBckGGInvMassPt; //! array of histos with pure gamma gamma combinatorial BG, invMass, pt + TH2F **fHistoTrueBckContInvMassPt; //! array of histos with contamination BG, invMass, pt + TH2F **fHistoTruePi0PtY; //! array of histos with validated pi0, pt, Y + TH2F **fHistoTrueEtaPtY; //! array of histos with validated eta, pt, Y + TH2F **fHistoTruePi0PtAlpha; //! array of histos with validated pi0, pt, alpha + TH2F **fHistoTrueEtaPtAlpha; //! array of histos with validated eta, pt, alpha + TH2F **fHistoTruePi0PtOpenAngle; //! array of histos with validated pi0, pt, openAngle + TH2F **fHistoTrueEtaPtOpenAngle; //! array of histos with validated eta, pt, openAngle + TH2F **fHistoTrueMotherDalitzInvMassPt; //! array of histos with validated mother, but Dalitz decay, invMass, pt + // MC validated reconstructed quantities photons + TH1F **fHistoTrueConvGammaPt; //! array of histos with validated conversion photon, pt + TH1F **fHistoTrueConvPi0GammaPt; //! array of histos with validated conversion photon from pi0, pt + TH1F **fHistoTrueConvGammaEta; //! array of histos with validated conversion photon, eta + TH2F **fHistoCombinatorialPt; //! array of histos with combinatorial BG, pt, source + TH1F **fHistoTruePrimaryConvGammaPt; //! array of histos with validated primary conversion photon, pt + TH2F **fHistoTruePrimaryConvGammaESDPtMCPt; //! array of histos with validated primary conversion photon, rec pt, mc pt + TH1F **fHistoTrueSecondaryConvGammaPt; //! array of histos with validated secondary conversion photon, pt + TH1F **fHistoTrueSecondaryConvGammaFromXFromK0sPt; //! array of histos with validated secondary conversion photon from K0s, pt + TH1F **fHistoTrueSecondaryConvGammaFromXFromLambdaPt;//! array of histos with validated secondary conversion photon from Lambda, pt + TH2F **fHistoTrueDalitzPsiPairDeltaPhi; //! array of histos with validated dalitz virtual photon, delta phi, psi pair + TH2F **fHistoTrueGammaPsiPairDeltaPhi; //! array of histos with validated conversion photon, delta phi, psi pair + TH1F ** fHistoTrueClusGammaPt; //! array of histos with validated cluster, pt + TH1F ** fHistoTruePrimaryClusGammaPt; //! array of histos with validated primary cluster, pt + TH2F ** fHistoTruePrimaryClusGammaESDPtMCPt; //! array of histos with validated primary cluster, rec Pt, MC pt + + // event histograms + TH1I **fHistoNEvents; //! array of histos with event information + TH1I **fHistoNGoodESDTracks; //! array of histos with number of good tracks (2010 Standard track cuts) + TH1I **fHistoNGammaCandidates; //! array of histos with number of gamma candidates per event + TH2F **fHistoNGoodESDTracksVsNGammaCanditates; //! array of histos with number of good tracks vs gamma candidates + TH1I **fHistoNV0Tracks; //! array of histos with V0 counts + TProfile **fProfileEtaShift; //! array of profiles with eta shift + + // additional variables + Double_t fEventPlaneAngle; // EventPlaneAngle + TRandom3 fRandom; // random + Int_t fNGammaCandidates; // number of gamma candidates in event + Double_t *fUnsmearedPx; //[fNGammaCandidates] + Double_t *fUnsmearedPy; //[fNGammaCandidates] + Double_t *fUnsmearedPz; //[fNGammaCandidates] + Double_t *fUnsmearedE; //[fNGammaCandidates] + Int_t *fMCStackPos; //[fNGammaCandidates] + Int_t *fMCStackNeg; //[fNGammaCandidates] + Int_t *fESDArrayPos; //[fNGammaCandidates] + Int_t *fESDArrayNeg; //[fNGammaCandidates] + Int_t fnCuts; // number of cuts to be analysed in parallel + Int_t fiCut; // current cut + Bool_t fMoveParticleAccordingToVertex; // boolean for BG calculation + Int_t fIsHeavyIon; // switch for pp = 0, PbPb = 1, pPb = 2 + Bool_t fDoMesonAnalysis; // flag for meson analysis + Int_t fDoMesonQA; // flag for meson QA + Int_t fDoPhotonQA; // flag for photon QA + Bool_t fIsFromMBHeader; // flag for MC headers + Bool_t fIsMC; // flag for MC information + + + // cluster cut variables + Double_t fMinE; + Int_t fNminCells; + Double_t fEMCm02cut; + //double fMinErat = 0; + //double fMinEcc = 0; + + + // TString fClusName; // cluster branch name (def="") + // const TObjArray *fRecPoints; // pointer to rec points (AliAnalysisTaskEMCALClusterizeFast) + // const TClonesArray *fDigits; // pointer to digits (AliAnalysisTaskEMCALClusterizeFast) + + private: + AliAnalysisTaskGammaConvCalo(const AliAnalysisTaskGammaConvCalo&); // Prevent copy-construction + AliAnalysisTaskGammaConvCalo &operator=(const AliAnalysisTaskGammaConvCalo&); // Prevent assignment + + ClassDef(AliAnalysisTaskGammaConvCalo, 1); +}; + +#endif diff --git a/PWGGA/GammaConv/AliCaloPhotonCuts.cxx b/PWGGA/GammaConv/AliCaloPhotonCuts.cxx new file mode 100644 index 00000000000..0ec31f28f59 --- /dev/null +++ b/PWGGA/GammaConv/AliCaloPhotonCuts.cxx @@ -0,0 +1,1268 @@ +/************************************************************************** + * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * * + * Authors: Friederike Bock, Baldo Sahlmueller * + * Version 1.0 * + * * + * 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. * + **************************************************************************/ + +//////////////////////////////////////////////// +//--------------------------------------------- +// Class handling all kinds of selection cuts for +// Photon from EMCAL clusters +//--------------------------------------------- +//////////////////////////////////////////////// + +#include "AliCaloPhotonCuts.h" +#include "AliAnalysisManager.h" +#include "AliInputEventHandler.h" +#include "AliMCEventHandler.h" +#include "AliAODHandler.h" +#include "TH1.h" +#include "TH2.h" +#include "TF1.h" +#include "AliStack.h" +#include "AliAODConversionMother.h" +#include "TObjString.h" +#include "AliAODEvent.h" +#include "AliESDEvent.h" +#include "AliCentrality.h" +#include "TList.h" +#include "TFile.h" +#include "AliLog.h" +#include "AliV0ReaderV1.h" +#include "AliAODMCParticle.h" +#include "AliAODMCHeader.h" + +class iostream; + +using namespace std; + +ClassImp(AliCaloPhotonCuts) + + +const char* AliCaloPhotonCuts::fgkCutNames[AliCaloPhotonCuts::kNCuts] = { + "ClusterType", //0 + "EtaMin", //1 + "EtaMax", //2 + "PhiMin", //3 + "PhiMax", //4 + "DistanceToBadChannel", //5 + "Timing", //6 + "TrackMatching", //7 + "ExoticCell", //8 + "MinEnergy", //9 + "MinNCells", //10 + "MinM02", //11 + "MaxM02", //12 + "MinM20", //13 + "MaxM20", //14 + "MaximumDispersion", //15 + "NLM" //16 +}; + + +//________________________________________________________________________ +AliCaloPhotonCuts::AliCaloPhotonCuts(const char *name,const char *title) : + AliAnalysisCuts(name,title), + fHistograms(NULL), + fClusterType(0), + fMinEtaCut(-10), + fMaxEtaCut(10), + fMinPhiCut(-10000), + fMaxPhiCut(-10000), + fMinDistanceToBadChannel(0), + fMaxTimeDiff(10e10), + fMinDistTrackToCluster(0), + fExoticCell(0), + fMinEnergy(0), + fMinNCells(0), + fMaxM02(1000), + fMinM02(0), + fMaxM20(1000), + fMinM20(0), + fMaxDispersion(1000), + fMinNLM(0), + fMaxNLM(1000), + fCutString(NULL), + fHistCutIndex(NULL), + fHistAcceptanceCuts(NULL), + fHistClusterIdentificationCuts(NULL), + fHistClusterEtavsPhiBeforeAcc(NULL), + fHistClusterEtavsPhiAfterAcc(NULL), + fHistClusterEtavsPhiAfterQA(NULL), + fHistDistanceToBadChannelBeforeAcc(NULL), + fHistDistanceToBadChannelAfterAcc(NULL), + fHistClusterTimevsEBeforeQA(NULL), + fHistClusterTimevsEAfterQA(NULL), + fHistExoticCellBeforeQA(NULL), + fHistExoticCellAfterQA(NULL), + fHistDistanceTrackToClusterBeforeQA(NULL), + fHistDistanceTrackToClusterAfterQA(NULL), + fHistEnergyOfClusterBeforeQA(NULL), + fHistEnergyOfClusterAfterQA(NULL), + fHistNCellsBeforeQA(NULL), + fHistNCellsAfterQA(NULL), + fHistM02BeforeQA(NULL), + fHistM02AfterQA(NULL), + fHistM20BeforeQA(NULL), + fHistM20AfterQA(NULL), + fHistDispersionBeforeQA(NULL), + fHistDispersionAfterQA(NULL), + fHistNLMBeforeQA(NULL), + fHistNLMAfterQA(NULL) +{ + for(Int_t jj=0;jjSetOwner(kTRUE); + if(name=="")fHistograms->SetName(Form("CaloCuts_%s",GetCutNumber().Data())); + else fHistograms->SetName(Form("%s_%s",name.Data(),GetCutNumber().Data())); + } + + // IsPhotonSelected + fHistCutIndex=new TH1F(Form("IsPhotonSelected %s",GetCutNumber().Data()),"IsPhotonSelected",5,-0.5,4.5); + fHistCutIndex->GetXaxis()->SetBinLabel(kPhotonIn+1,"in"); + fHistCutIndex->GetXaxis()->SetBinLabel(kDetector+1,"detector"); + fHistCutIndex->GetXaxis()->SetBinLabel(kAcceptance+1,"acceptance"); + fHistCutIndex->GetXaxis()->SetBinLabel(kClusterQuality+1,"cluster QA"); + fHistCutIndex->GetXaxis()->SetBinLabel(kPhotonOut+1,"out"); + fHistograms->Add(fHistCutIndex); + + // Acceptance Cuts + fHistAcceptanceCuts=new TH1F(Form("AcceptanceCuts %s",GetCutNumber().Data()),"AcceptanceCuts",5,-0.5,4.5); + fHistAcceptanceCuts->GetXaxis()->SetBinLabel(1,"in"); + fHistAcceptanceCuts->GetXaxis()->SetBinLabel(2,"eta"); + fHistAcceptanceCuts->GetXaxis()->SetBinLabel(3,"phi"); + fHistAcceptanceCuts->GetXaxis()->SetBinLabel(4,"distance to bad channel"); + fHistAcceptanceCuts->GetXaxis()->SetBinLabel(5,"out"); + fHistograms->Add(fHistAcceptanceCuts); + + // Cluster Cuts + fHistClusterIdentificationCuts =new TH1F(Form("ClusterQualityCuts %s",GetCutNumber().Data()),"ClusterQualityCuts",10,-0.5,9.5); + fHistClusterIdentificationCuts->GetXaxis()->SetBinLabel(1,"in"); + fHistClusterIdentificationCuts->GetXaxis()->SetBinLabel(2,"timing"); + fHistClusterIdentificationCuts->GetXaxis()->SetBinLabel(3,"track matching"); + fHistClusterIdentificationCuts->GetXaxis()->SetBinLabel(4,"Exotics"); + fHistClusterIdentificationCuts->GetXaxis()->SetBinLabel(5,"minimum energy"); + fHistClusterIdentificationCuts->GetXaxis()->SetBinLabel(6,"M02"); + fHistClusterIdentificationCuts->GetXaxis()->SetBinLabel(7,"M20"); + fHistClusterIdentificationCuts->GetXaxis()->SetBinLabel(8,"dispersion"); + fHistClusterIdentificationCuts->GetXaxis()->SetBinLabel(9,"NLM"); + fHistClusterIdentificationCuts->GetXaxis()->SetBinLabel(10,"out"); + fHistograms->Add(fHistClusterIdentificationCuts); + + // Acceptance related histogramms + fHistClusterEtavsPhiBeforeAcc=new TH2F(Form("EtaPhi_beforeAcceptance %s",GetCutNumber().Data()),"EtaPhi_beforeAcceptance",462,-TMath::Pi(),TMath::Pi(),110,-0.7,0.7); + fHistograms->Add(fHistClusterEtavsPhiBeforeAcc); + fHistClusterEtavsPhiAfterAcc=new TH2F(Form("EtaPhi_afterAcceptance %s",GetCutNumber().Data()),"EtaPhi_afterAcceptance",462,-TMath::Pi(),TMath::Pi(),110,-0.7,0.7); + fHistograms->Add(fHistClusterEtavsPhiAfterAcc); + fHistClusterEtavsPhiAfterQA=new TH2F(Form("EtaPhi_afterClusterQA %s",GetCutNumber().Data()),"EtaPhi_afterClusterQA",462,-TMath::Pi(),TMath::Pi(),110,-0.7,0.7); + fHistograms->Add(fHistClusterEtavsPhiAfterQA); + fHistDistanceToBadChannelBeforeAcc = new TH1F(Form("DistanceToBadChannel_beforeAcceptance %s",GetCutNumber().Data()),"DistanceToBadChannel_beforeAcceptance",200,0,40); + fHistograms->Add(fHistDistanceToBadChannelBeforeAcc); + fHistDistanceToBadChannelAfterAcc = new TH1F(Form("DistanceToBadChannel_afterAcceptance %s",GetCutNumber().Data()),"DistanceToBadChannel_afterAcceptance",200,0,40); + fHistograms->Add(fHistDistanceToBadChannelAfterAcc); + + // Cluster quality related histograms + fHistClusterTimevsEBeforeQA=new TH2F(Form("ClusterTimeVsE_beforeClusterQA %s",GetCutNumber().Data()),"ClusterTimeVsE_beforeClusterQA",400,-10e-6,10e-6,100,0.,40); + fHistograms->Add(fHistClusterTimevsEBeforeQA); + fHistClusterTimevsEAfterQA=new TH2F(Form("ClusterTimeVsE_afterClusterQA %s",GetCutNumber().Data()),"ClusterTimeVsE_afterClusterQA",400,-10e-6,10e-6,100,0.,40); + fHistograms->Add(fHistClusterTimevsEAfterQA); + fHistExoticCellBeforeQA=new TH2F(Form("ExoticCell_beforeClusterQA %s",GetCutNumber().Data()),"ExoticCell_beforeClusterQA",400,0,40,50,0.75,1); + fHistograms->Add(fHistExoticCellBeforeQA); + fHistExoticCellAfterQA=new TH2F(Form("ExoticCell_afterClusterQA %s",GetCutNumber().Data()),"ExoticCell_afterClusterQA",400,0,40,50,0.75,1); + fHistograms->Add(fHistExoticCellAfterQA); + fHistDistanceTrackToClusterBeforeQA = new TH1F(Form("DistanceToTrack_beforeClusterQA %s",GetCutNumber().Data()),"DistanceToTrack_beforeClusterQA",200,0,40); + fHistograms->Add(fHistDistanceTrackToClusterBeforeQA); + fHistDistanceTrackToClusterAfterQA = new TH1F(Form("DistanceToTrack_afterClusterQA %s",GetCutNumber().Data()),"DistanceToTrack_afterClusterQA",200,0,40); + fHistograms->Add(fHistDistanceTrackToClusterAfterQA); + fHistEnergyOfClusterBeforeQA = new TH1F(Form("EnergyOfCluster_beforeClusterQA %s",GetCutNumber().Data()),"EnergyOfCluster_beforeClusterQA",300,0,30); + fHistograms->Add(fHistEnergyOfClusterBeforeQA); + fHistEnergyOfClusterAfterQA = new TH1F(Form("EnergyOfCluster_afterClusterQA %s",GetCutNumber().Data()),"EnergyOfCluster_afterClusterQA",300,0,30); + fHistograms->Add(fHistEnergyOfClusterAfterQA); + fHistNCellsBeforeQA = new TH1F(Form("NCellPerCluster_beforeClusterQA %s",GetCutNumber().Data()),"NCellPerCluster_beforeClusterQA",50,0,50); + fHistograms->Add(fHistNCellsBeforeQA); + fHistNCellsAfterQA = new TH1F(Form("NCellPerCluster_afterClusterQA %s",GetCutNumber().Data()),"NCellPerCluster_afterClusterQA",50,0,50); + fHistograms->Add(fHistNCellsAfterQA); + fHistM02BeforeQA = new TH1F(Form("M02_beforeClusterQA %s",GetCutNumber().Data()),"M02_beforeClusterQA",100,0,5); + fHistograms->Add(fHistM02BeforeQA); + fHistM02AfterQA = new TH1F(Form("M02_afterClusterQA %s",GetCutNumber().Data()),"M02_afterClusterQA",100,0,5); + fHistograms->Add(fHistM02AfterQA); + fHistM20BeforeQA = new TH1F(Form("M20_beforeClusterQA %s",GetCutNumber().Data()),"M20_beforeClusterQA",100,0,2.5); + fHistograms->Add(fHistM20BeforeQA); + fHistM20AfterQA = new TH1F(Form("M20_afterClusterQA %s",GetCutNumber().Data()),"M20_afterClusterQA",100,0,2.5); + fHistograms->Add(fHistM20AfterQA); + fHistDispersionBeforeQA = new TH1F(Form("Dispersion_beforeClusterQA %s",GetCutNumber().Data()),"Dispersion_beforeClusterQA",100,0,4); + fHistograms->Add(fHistDispersionBeforeQA); + fHistDispersionAfterQA = new TH1F(Form("Dispersion_afterClusterQA %s",GetCutNumber().Data()),"Dispersion_afterClusterQA",100,0,4); + fHistograms->Add(fHistDispersionAfterQA); + fHistNLMBeforeQA = new TH1F(Form("NLM_beforeClusterQA %s",GetCutNumber().Data()),"NLM_beforeClusterQA",10,0,10); + fHistograms->Add(fHistNLMBeforeQA); + fHistNLMAfterQA = new TH1F(Form("NLM_afterClusterQA %s",GetCutNumber().Data()),"NLM_afterClusterQA",10,0,10); + fHistograms->Add(fHistNLMAfterQA); + + TH1::AddDirectory(kTRUE); +} + +/* +///________________________________________________________________________ +Bool_t AliCaloPhotonCuts::ClusterIsSelectedMC(TParticle *particle,AliStack *fMCStack,Bool_t checkForConvertedGamma){ + // MonteCarlo Photon Selection + + if(!fMCStack)return kFALSE; + + if (particle->GetPdgCode() == 22){ + + + if( particle->Eta() > (fEtaCut) || particle->Eta() < (-fEtaCut) ) + return kFALSE; + if(fEtaCutMin>-0.1){ + if( particle->Eta() < (fEtaCutMin) && particle->Eta() > (-fEtaCutMin) ) + return kFALSE; + } + + if(particle->GetMother(0) >-1 && fMCStack->Particle(particle->GetMother(0))->GetPdgCode() == 22){ + return kFALSE; // no photon as mothers! + } + + if(particle->GetMother(0) >= fMCStack->GetNprimary()){ + return kFALSE; // the gamma has a mother, and it is not a primary particle + } + + if(!checkForConvertedGamma) return kTRUE; // return in case of accepted gamma + + // looking for conversion gammas (electron + positron from pairbuilding (= 5) ) + TParticle* ePos = NULL; + TParticle* eNeg = NULL; + + if(particle->GetNDaughters() >= 2){ + for(Int_t daughterIndex=particle->GetFirstDaughter();daughterIndex<=particle->GetLastDaughter();daughterIndex++){ + TParticle *tmpDaughter = fMCStack->Particle(daughterIndex); + if(tmpDaughter->GetUniqueID() == 5){ + if(tmpDaughter->GetPdgCode() == 11){ + eNeg = tmpDaughter; + } else if(tmpDaughter->GetPdgCode() == -11){ + ePos = tmpDaughter; + } + } + } + } + + if(ePos == NULL || eNeg == NULL){ // means we do not have two daughters from pair production + return kFALSE; + } + + if(ePos->Pt()Pt()Eta() > (fEtaCut) || ePos->Eta() < (-fEtaCut) || + eNeg->Eta() > (fEtaCut) || eNeg->Eta() < (-fEtaCut) ) + return kFALSE; + + if(fEtaCutMin > -0.1){ + if( (ePos->Eta() < (fEtaCutMin) && ePos->Eta() > (-fEtaCutMin)) || + (eNeg->Eta() < (fEtaCutMin) && eNeg->Eta() > (-fEtaCutMin)) ) + return kFALSE; + } + + if(ePos->R()>fMaxR){ + return kFALSE; // cuts on distance from collision point + } + + if(abs(ePos->Vz()) > fMaxZ){ + return kFALSE; // outside material + } + if(abs(eNeg->Vz()) > fMaxZ){ + return kFALSE; // outside material + } + + if( ePos->R() <= ((abs(ePos->Vz()) * fLineCutZRSlope) - fLineCutZValue)){ + return kFALSE; // line cut to exclude regions where we do not reconstruct + } else if ( fEtaCutMin != -0.1 && ePos->R() >= ((abs(ePos->Vz()) * fLineCutZRSlopeMin) - fLineCutZValueMin)){ + return kFALSE; + } + + if( eNeg->R() <= ((abs(eNeg->Vz()) * fLineCutZRSlope) - fLineCutZValue)){ + return kFALSE; // line cut to exclude regions where we do not reconstruct + } else if ( fEtaCutMin != -0.1 && eNeg->R() >= ((abs(eNeg->Vz()) * fLineCutZRSlopeMin) - fLineCutZValueMin)){ + return kFALSE; + } + + return kTRUE; + //if(AcceptanceCut(particle,ePos,eNeg))return kTRUE; + } + return kFALSE; +} +///________________________________________________________________________ +Bool_t AliCaloPhotonCuts::ClusterIsSelectedAODMC(AliAODMCParticle *particle,TClonesArray *aodmcArray,Bool_t checkForConvertedGamma){ + // MonteCarlo Photon Selection + + if(!aodmcArray)return kFALSE; + + if (particle->GetPdgCode() == 22){ + if( particle->Eta() > (fEtaCut) || particle->Eta() < (-fEtaCut) ) + return kFALSE; + if(fEtaCutMin>-0.1){ + if( particle->Eta() < (fEtaCutMin) && particle->Eta() > (-fEtaCutMin) ) + return kFALSE; + } + + if(particle->GetMother() > -1){ + if((static_cast(aodmcArray->At(particle->GetMother())))->GetPdgCode() == 22){ + return kFALSE; // no photon as mothers! + } + if(!(static_cast(aodmcArray->At(particle->GetMother()))->IsPrimary())){ + return kFALSE; // the gamma has a mother, and it is not a primary particle + } + } + + if(!checkForConvertedGamma) return kTRUE; // return in case of accepted gamma + + // looking for conversion gammas (electron + positron from pairbuilding (= 5) ) + AliAODMCParticle* ePos = NULL; + AliAODMCParticle* eNeg = NULL; + + if(particle->GetNDaughters() >= 2){ + for(Int_t daughterIndex=particle->GetDaughter(0);daughterIndex<=particle->GetDaughter(1);daughterIndex++){ + AliAODMCParticle *tmpDaughter = static_cast(aodmcArray->At(daughterIndex)); + if(!tmpDaughter) continue; + if(((tmpDaughter->GetMCProcessCode())) == 5){ // STILL A BUG IN ALIROOT >>8 HAS TPO BE REMOVED AFTER FIX + if(tmpDaughter->GetPdgCode() == 11){ + eNeg = tmpDaughter; + } else if(tmpDaughter->GetPdgCode() == -11){ + ePos = tmpDaughter; + } + } + } + } + + if(ePos == NULL || eNeg == NULL){ // means we do not have two daughters from pair production + return kFALSE; + } + + if(ePos->Pt()Pt()Eta() > (fEtaCut) || ePos->Eta() < (-fEtaCut) || + eNeg->Eta() > (fEtaCut) || eNeg->Eta() < (-fEtaCut) ) + return kFALSE; + + if(fEtaCutMin > -0.1){ + if( (ePos->Eta() < (fEtaCutMin) && ePos->Eta() > (-fEtaCutMin)) || + (eNeg->Eta() < (fEtaCutMin) && eNeg->Eta() > (-fEtaCutMin)) ) + return kFALSE; + } + + Double_t rPos = sqrt( (ePos->Xv()*ePos->Xv()) + (ePos->Yv()*ePos->Yv()) ); + Double_t rNeg = sqrt( (eNeg->Xv()*eNeg->Xv()) + (eNeg->Yv()*eNeg->Yv()) ); + + if(rPos>fMaxR){ + return kFALSE; // cuts on distance from collision point + } + if(abs(ePos->Zv()) > fMaxZ){ + return kFALSE; // outside material + } + if(abs(eNeg->Zv()) > fMaxZ){ + return kFALSE; // outside material + } + + if( rPos <= ((abs(ePos->Zv()) * fLineCutZRSlope) - fLineCutZValue)){ + return kFALSE; // line cut to exclude regions where we do not reconstruct + } else if ( fEtaCutMin != -0.1 && rPos >= ((abs(ePos->Zv()) * fLineCutZRSlopeMin) - fLineCutZValueMin)){ + return kFALSE; + } + + if( rNeg <= ((abs(eNeg->Zv()) * fLineCutZRSlope) - fLineCutZValue)){ + return kFALSE; // line cut to exclude regions where we do not reconstruct + } else if ( fEtaCutMin != -0.1 && rNeg >= ((abs(eNeg->Zv()) * fLineCutZRSlopeMin) - fLineCutZValueMin)){ + return kFALSE; + } + + return kTRUE; + //if(AcceptanceCut(particle,ePos,eNeg))return kTRUE; + } + return kFALSE; +}*/ + + + +///________________________________________________________________________ +// This function selects the clusters based on their quality criteria +///________________________________________________________________________ +Bool_t AliCaloPhotonCuts::ClusterQualityCuts(AliVCluster* cluster, AliVEvent *event, Bool_t isMC) +{ // Specific Photon Cuts + + Int_t cutIndex = 0; + if(fHistClusterIdentificationCuts)fHistClusterIdentificationCuts->Fill(cutIndex); + cutIndex++; + + + // Fill Histos before Cuts + if(fHistClusterTimevsEBeforeQA) fHistClusterTimevsEBeforeQA->Fill(cluster->GetTOF(), cluster->E()); +// if(fHistExoticCellBeforeQA) fHistExoticCellBeforeQA->Fill(cluster->E(), ); + if(fHistDistanceTrackToClusterBeforeQA) fHistDistanceTrackToClusterBeforeQA->Fill(cluster->GetEmcCpvDistance()); + if(fHistEnergyOfClusterBeforeQA) fHistEnergyOfClusterBeforeQA->Fill(cluster->E()); + if(fHistNCellsBeforeQA) fHistNCellsBeforeQA->Fill(cluster->GetNCells()); + if(fHistM02BeforeQA) fHistM02BeforeQA->Fill(cluster->GetM02()); + if(fHistM20BeforeQA) fHistM20BeforeQA->Fill(cluster->GetM20()); + if(fHistDispersionBeforeQA) fHistDispersionBeforeQA->Fill(cluster->GetDispersion()); +// if(fHistNLMBeforeQA) fHistNLMBeforeQA->Fill(cluster->GetNExMax()); + + // Check wether timing is ok + if(abs(cluster->GetTOF()) > fMaxTimeDiff && !isMC){ + if(fHistClusterIdentificationCuts)fHistClusterIdentificationCuts->Fill(cutIndex); //1 + return kFALSE; + } + cutIndex++; //2, next cut + + // Minimum distance to track + if(cluster->GetEmcCpvDistance() < fMinDistTrackToCluster){ + if(fHistClusterIdentificationCuts)fHistClusterIdentificationCuts->Fill(cutIndex); //2 + return kFALSE; + } + cutIndex++;//3, next cut + + // exotic cell cut --IMPLEMENT LATER--- +// if(!AcceptanceCuts(photon)){ +// if(fHistClusterIdentificationCuts)fHistClusterIdentificationCuts->Fill(cutIndex); //3 +// return kFALSE; +// } + cutIndex++; //4, next cut + + // minimum cell energy cut + if(cluster->E() < fMinEnergy){ + if(fHistClusterIdentificationCuts)fHistClusterIdentificationCuts->Fill(cutIndex); //4 + return kFALSE; + } + cutIndex++; //5, next cut + + // minimum number of cells + if(cluster->GetNCells() < fMinNCells) { + if(fHistClusterIdentificationCuts)fHistClusterIdentificationCuts->Fill(cutIndex); //5 + return kFALSE; + } + cutIndex++; //6, next cut + + // M02 cut + if( cluster->GetM02()< fMinM02 || cluster->GetM02() > fMaxM02 ) { + if(fHistClusterIdentificationCuts)fHistClusterIdentificationCuts->Fill(cutIndex); //6 + return kFALSE; + } + cutIndex++; //7, next cut + + // M20 cut + if( cluster->GetM20()< fMinM20 || cluster->GetM20() > fMaxM20 ) { + if(fHistClusterIdentificationCuts)fHistClusterIdentificationCuts->Fill(cutIndex); //7 + return kFALSE; + } + cutIndex++; //8, next cut + + // dispersion cut + if( cluster->GetDispersion()> fMaxDispersion) { + if(fHistClusterIdentificationCuts)fHistClusterIdentificationCuts->Fill(cutIndex); //8 + return kFALSE; + } + cutIndex++; //9, next cut + + // NLM cut --IMPLEMENT LATER--- + if( cluster->GetDispersion()> fMaxDispersion) { + if(fHistClusterIdentificationCuts)fHistClusterIdentificationCuts->Fill(cutIndex); //9 + return kFALSE; + } + cutIndex++; //9, next cut + + // DONE with selecting photons + if(fHistClusterIdentificationCuts)fHistClusterIdentificationCuts->Fill(cutIndex); //10 + + // Histos after Cuts + Double_t vertex[3] = {0}; + event->GetPrimaryVertex()->GetXYZ(vertex); + // TLorentzvector with cluster + TLorentzVector clusterVector; + cluster->GetMomentum(clusterVector,vertex); + Double_t etaCluster = clusterVector.Eta(); + Double_t phiCluster = clusterVector.Phi(); + + if(fHistClusterEtavsPhiAfterQA) fHistClusterEtavsPhiAfterQA->Fill(phiCluster,etaCluster); + if(fHistClusterTimevsEAfterQA) fHistClusterTimevsEAfterQA->Fill(cluster->GetTOF(), cluster->E()); +// if(fHistExoticCellAfterQA) fHistExoticCellAfterQA->Fill(cluster->E(), ); + if(fHistDistanceTrackToClusterAfterQA) fHistDistanceTrackToClusterAfterQA->Fill(cluster->GetEmcCpvDistance()); + if(fHistEnergyOfClusterAfterQA) fHistEnergyOfClusterAfterQA->Fill(cluster->E()); + if(fHistNCellsAfterQA) fHistNCellsAfterQA->Fill(cluster->GetNCells()); + if(fHistM02AfterQA) fHistM02AfterQA->Fill(cluster->GetM02()); + if(fHistM20AfterQA) fHistM20AfterQA->Fill(cluster->GetM20()); + if(fHistDispersionAfterQA) fHistDispersionAfterQA->Fill(cluster->GetDispersion()); +// if(fHistNLMBeforeQA) fHistNLMAfterQA->Fill(cluster->GetNExMax()); + + return kTRUE; + +} + + +///________________________________________________________________________ +Bool_t AliCaloPhotonCuts::ClusterIsSelected(AliVCluster *cluster, AliVEvent * event, Bool_t isMC) +{ + //Selection of Reconstructed photon clusters with Calorimeters + + FillClusterCutIndex(kPhotonIn); + + Double_t vertex[3] = {0}; + event->GetPrimaryVertex()->GetXYZ(vertex); + // TLorentzvector with cluster + TLorentzVector clusterVector; + cluster->GetMomentum(clusterVector,vertex); + Double_t etaCluster = clusterVector.Eta(); + Double_t phiCluster = clusterVector.Phi(); + + // Histos before cuts + if(fHistClusterEtavsPhiBeforeAcc) fHistClusterEtavsPhiBeforeAcc->Fill(phiCluster,etaCluster); + + // Cluster Selection - 0= accept any calo cluster + if (fClusterType > 0){ + //Select EMCAL cluster + if (fClusterType == 1 && !cluster->IsEMCAL()){ + FillClusterCutIndex(kDetector); + } + //Select PHOS cluster + if (fClusterType == 2 && !cluster->IsPHOS()){ + FillClusterCutIndex(kDetector); + } + } + + // Acceptance Cuts + if(!AcceptanceCuts(cluster,event)){ + FillClusterCutIndex(kAcceptance); + return kFALSE; + } + // Cluster Quality Cuts + if(!ClusterQualityCuts(cluster,event,isMC)){ + FillClusterCutIndex(kClusterQuality); + return kFALSE; + } + + // Photon passed cuts + FillClusterCutIndex(kPhotonOut); + return kTRUE; +} + + +///________________________________________________________________________ +Bool_t AliCaloPhotonCuts::AcceptanceCuts(AliVCluster *cluster, AliVEvent* event) +{ + // Exclude certain areas for photon reconstruction + + Int_t cutIndex=0; + if(fHistAcceptanceCuts)fHistAcceptanceCuts->Fill(cutIndex); + cutIndex++; + + + Double_t vertex[3] = {0}; + event->GetPrimaryVertex()->GetXYZ(vertex); + // TLorentzvector with cluster + TLorentzVector clusterVector; + cluster->GetMomentum(clusterVector,vertex); + Double_t etaCluster = clusterVector.Eta(); + Double_t phiCluster = clusterVector.Phi(); + + // check eta range + if (fMinEtaCut !=-10 && fMaxEtaCut !=10 ){ + if (etaCluster < fMinEtaCut || etaCluster > fMaxEtaCut){ + if(fHistAcceptanceCuts)fHistAcceptanceCuts->Fill(cutIndex); + return kFALSE; + } + } + cutIndex++; + + // check phi range + if (fMinPhiCut !=-10000 && fMaxPhiCut !=-10000 ){ + if (phiCluster < fMinPhiCut || phiCluster > fMaxEtaCut){ + if(fHistAcceptanceCuts)fHistAcceptanceCuts->Fill(cutIndex); + return kFALSE; + } + } + cutIndex++; + + // check distance to bad channel + if (cluster->GetDistanceToBadChannel() < fMinDistanceToBadChannel){ + if(fHistAcceptanceCuts)fHistAcceptanceCuts->Fill(cutIndex); + return kFALSE; + } + + if(fHistAcceptanceCuts)fHistAcceptanceCuts->Fill(cutIndex); + + // Histos after cuts + if(fHistClusterEtavsPhiAfterAcc) fHistClusterEtavsPhiAfterAcc->Fill(phiCluster,etaCluster); + + return kTRUE; +} + +///________________________________________________________________________ +Bool_t AliCaloPhotonCuts::UpdateCutString() { + ///Update the cut string (if it has been created yet) + + if(fCutString && fCutString->GetString().Length() == kNCuts) { + fCutString->SetString(GetCutNumber()); + } else { + return kFALSE; + } + return kTRUE; +} + +///________________________________________________________________________ +Bool_t AliCaloPhotonCuts::InitializeCutsFromCutString(const TString analysisCutSelection ) { + // Initialize Cuts from a given Cut string + AliInfo(Form("Set CaloCut Number: %s",analysisCutSelection.Data())); + if(analysisCutSelection.Length()!=kNCuts) { + AliError(Form("Cut selection has the wrong length! size is %d, number of cuts is %d", analysisCutSelection.Length(), kNCuts)); + return kFALSE; + } + if(!analysisCutSelection.IsDigit()){ + AliError("Cut selection contains characters"); + return kFALSE; + } + + const char *cutSelection = analysisCutSelection.Data(); + #define ASSIGNARRAY(i) fCuts[i] = cutSelection[i] - '0' + for(Int_t ii=0;ii %3.2f\n", fMinDistanceToBadChannel ); + + printf("Cluster Quality cuts: \n"); + printf("\t time difference < %3.2f\n", fMaxTimeDiff ); + printf("\t min distance to track > %3.2f\n", fMinDistTrackToCluster ); + printf("\t E_{cluster} > %3.2f\n", fMinEnergy ); + printf("\t %3.2f < M02 < %3.2f\n", fMinM02, fMaxM02 ); + +} + +// EMCAL acceptance 2011 +// 1.39626, 3.125 (phi) +// -0.66687,,0.66465 + + +///________________________________________________________________________ +Bool_t AliCaloPhotonCuts::SetClusterTypeCut(Int_t clusterType) +{ // Set Cut + switch(clusterType){ + case 0: // all clusters + fClusterType=0; + break; + case 1: // EMCAL clusters + fClusterType=1; + break; + case 2: // PHOS clusters + fClusterType=2; + break; + default: + AliError(Form("ClusterTypeCut not defined %d",clusterType)); + return kFALSE; + } + return kTRUE; +} + +//___________________________________________________________________ +Bool_t AliCaloPhotonCuts::SetMinEtaCut(Int_t minEta) +{ + switch(minEta){ + case 0: + fMinEtaCut=-0.6687; + break; + case 1: + fClusterType=-0.5; + break; + case 2: + fClusterType=-2; + break; + default: + AliError(Form("MinEta Cut not defined %d",minEta)); + return kFALSE; + } + return kTRUE; +} + + +//___________________________________________________________________ +Bool_t AliCaloPhotonCuts::SetMaxEtaCut(Int_t maxEta) +{ + switch(maxEta){ + case 0: + fMaxEtaCut=0.66465; + break; + case 1: + fMaxEtaCut=0.5; + break; + case 2: + fMaxEtaCut=2; + break; + default: + AliError(Form("MaxEta Cut not defined %d",maxEta)); + return kFALSE; + } + return kTRUE; +} + +//___________________________________________________________________ +Bool_t AliCaloPhotonCuts::SetMinPhiCut(Int_t minPhi) +{ + switch(minPhi){ + case 0: + fMinPhiCut=-10000; + break; + case 1: + fMinPhiCut=1.39626; + break; + default: + AliError(Form("MinPhi Cut not defined %d",minPhi)); + return kFALSE; + } + return kTRUE; +} + +//___________________________________________________________________ +Bool_t AliCaloPhotonCuts::SetMaxPhiCut(Int_t maxPhi) +{ + switch(maxPhi){ + case 0: + fMaxPhiCut=-10000; + break; + case 1: + fMaxPhiCut=3.125; + break; + default: + AliError(Form("Max Phi Cut not defined %d",maxPhi)); + return kFALSE; + } + return kTRUE; +} + +//___________________________________________________________________ +Bool_t AliCaloPhotonCuts::SetDistanceToBadChannelCut(Int_t distanceToBadChannel) +{ + switch(distanceToBadChannel){ + case 0: + fMinDistanceToBadChannel=0; + break; + case 1: + fMinDistanceToBadChannel=5; + break; + default: + AliError(Form("minimum distance to bad channel Cut not defined %d",distanceToBadChannel)); + return kFALSE; + } + return kTRUE; +} + +//___________________________________________________________________ +Bool_t AliCaloPhotonCuts::SetTimingCut(Int_t timing) +{ + switch(timing){ + case 0: + fMaxTimeDiff=500; + break; + case 1: + fMaxTimeDiff=10e-7; //1000ns + break; + case 2: + fMaxTimeDiff=50e-8; //500ns + break; + case 3: + fMaxTimeDiff=20e-8; //200ns + break; + case 4: + fMaxTimeDiff=10e-8; //100ns + break; + case 5: + fMaxTimeDiff=50e-9; //50ns + break; + + default: + AliError(Form("Timing Cut not defined %d",timing)); + return kFALSE; + } + return kTRUE; +} + +//___________________________________________________________________ +Bool_t AliCaloPhotonCuts::SetTrackMatchingCut(Int_t trackMatching) +{ + switch(trackMatching){ + case 0: + fMinDistTrackToCluster=0; + break; + case 1: + fMinDistTrackToCluster=5; + break; + default: + AliError(Form("Track Matching Cut not defined %d",trackMatching)); + return kFALSE; + } + return kTRUE; +} + +//___________________________________________________________________ +Bool_t AliCaloPhotonCuts::SetExoticCellCut(Int_t exoticCell) +{ + switch(exoticCell){ + case 0: + fExoticCell=0; + break; + case 1: + fExoticCell=5; + break; + default: + AliError(Form("Exotic cell Cut not defined %d",exoticCell)); + return kFALSE; + } + return kTRUE; +} + +//___________________________________________________________________ +Bool_t AliCaloPhotonCuts::SetMinEnergyCut(Int_t minEnergy) +{ + switch(minEnergy){ + case 0: + fMinEnergy=0; + break; + case 1: + fMinEnergy=0.05; + break; + case 2: + fMinEnergy=0.1; + break; + case 3: + fMinEnergy=0.15; + break; + case 4: + fMinEnergy=0.2; + break; + case 5: + fMinEnergy=0.3; + break; + case 6: + fMinEnergy=0.5; + break; + case 7: + fMinEnergy=0.75; + break; + case 8: + fMinEnergy=1.; + break; + case 9: + fMinEnergy=1.25; + break; + default: + AliError(Form("Minimum Energy Cut not defined %d",minEnergy)); + return kFALSE; + } + return kTRUE; +} + +//___________________________________________________________________ +Bool_t AliCaloPhotonCuts::SetMinNCellsCut(Int_t minNCells) +{ + switch(minNCells){ + case 0: + fMinNCells=0; + break; + case 1: + fMinNCells=1; + break; + case 2: + fMinNCells=2; + break; + case 3: + fMinNCells=3; + break; + case 4: + fMinNCells=4; + break; + case 5: + fMinNCells=5; + break; + case 6: + fMinNCells=6; + break; + + default: + AliError(Form("Min N cells Cut not defined %d",minNCells)); + return kFALSE; + } + return kTRUE; +} + +//___________________________________________________________________ +Bool_t AliCaloPhotonCuts::SetMaxM02(Int_t maxM02) +{ + switch(maxM02){ + case 0: + fMaxM02=100; + break; + case 1: + fMaxM02=1.; + break; + case 2: + fMaxM02=0.7; + break; + case 3: + fMaxM02=0.5; + break; + case 4: + fMaxM02=0.4; + break; + + default: + AliError(Form("Max M02 Cut not defined %d",maxM02)); + return kFALSE; + } + return kTRUE; +} + +//___________________________________________________________________ +Bool_t AliCaloPhotonCuts::SetMinM02(Int_t minM02) +{ + switch(minM02){ + case 0: + fMinM02=0; + break; + case 1: + fMinM02=0.002; + break; + default: + AliError(Form("Min M02 not defined %d",minM02)); + return kFALSE; + } + return kTRUE; +} + +//___________________________________________________________________ +Bool_t AliCaloPhotonCuts::SetMaxM20(Int_t maxM20) +{ + switch(maxM20){ + case 0: + fMaxM20=100; + break; + case 1: + fMaxM20=0.5; + break; + default: + AliError(Form("Max M20 Cut not defined %d",maxM20)); + return kFALSE; + } + return kTRUE; +} + +//___________________________________________________________________ +Bool_t AliCaloPhotonCuts::SetMinM20(Int_t minM20) +{ + switch(minM20){ + case 0: + fMinM20=0; + break; + case 1: + fMinM20=0.002; + break; + default: + AliError(Form("Min M20 Cut not defined %d",minM20)); + return kFALSE; + } + return kTRUE; +} + +//___________________________________________________________________ +Bool_t AliCaloPhotonCuts::SetDispersion(Int_t dispersion) +{ + switch(dispersion){ + case 0: + fMaxDispersion =100; + break; + case 1: + fMaxDispersion=2.; + break; + default: + AliError(Form("Maximum Dispersion Cut not defined %d",dispersion)); + return kFALSE; + } + return kTRUE; +} + +//___________________________________________________________________ +Bool_t AliCaloPhotonCuts::SetNLM(Int_t nlm) +{ + switch(nlm){ + case 0: + fMinNLM =0; + fMaxNLM =100; + break; + case 1: + fMinNLM =0; + fMaxNLM =1; + break; + default: + AliError(Form("NLM Cut not defined %d",nlm)); + return kFALSE; + } + return kTRUE; +} + +///________________________________________________________________________ +TString AliCaloPhotonCuts::GetCutNumber(){ + // returns TString with current cut number + TString a(kNCuts); + for(Int_t ii=0;iiFill(photoncut);} + + ///Cut functions + Bool_t AcceptanceCuts(AliVCluster* cluster, AliVEvent *event); + Bool_t ClusterQualityCuts(AliVCluster* cluster,AliVEvent *event, Bool_t isMC); + + // Set Individual Cuts + Bool_t SetClusterTypeCut(Int_t); + Bool_t SetMinEtaCut(Int_t); + Bool_t SetMaxEtaCut(Int_t); + Bool_t SetMinPhiCut(Int_t); + Bool_t SetMaxPhiCut(Int_t); + Bool_t SetDistanceToBadChannelCut(Int_t); + Bool_t SetTimingCut(Int_t); + Bool_t SetTrackMatchingCut(Int_t); + Bool_t SetExoticCellCut(Int_t); + Bool_t SetMinEnergyCut(Int_t); + Bool_t SetMinNCellsCut(Int_t); + Bool_t SetMaxM02(Int_t); + Bool_t SetMinM02(Int_t); + Bool_t SetMaxM20(Int_t); + Bool_t SetMinM20(Int_t); + Bool_t SetDispersion(Int_t); + Bool_t SetNLM(Int_t); + + protected: + TList *fHistograms; + + //cuts + Int_t fClusterType; // which cluster do we have + Double_t fMinEtaCut; // min eta cut + Double_t fMaxEtaCut; // max eta cut + Double_t fMinPhiCut; // phi cut + Double_t fMaxPhiCut; // phi cut + Double_t fMinDistanceToBadChannel; // minimum distance to bad channel + Double_t fMaxTimeDiff; // maximum time difference to triggered collision + Double_t fMinDistTrackToCluster; // minimum distance between track and cluster + Double_t fExoticCell; // exotic cell cut + Double_t fMinEnergy; // minium energy per cluster + Int_t fMinNCells; // minimum number of cells + Double_t fMaxM02; // maximum M02 + Double_t fMinM02; // minimum M02 + Double_t fMaxM20; // maximum M20 + Double_t fMinM20; // minimum M20 + Double_t fMaxDispersion; // maximum dispersion + Int_t fMinNLM; // minimum number of local maxima in cluster + Int_t fMaxNLM; // maximum number of local maxima in cluster + + // CutString + TObjString *fCutString; // cut number used for analysis + + // Histograms + TH1F *fHistCutIndex; // bookkeeping for cuts + TH1F *fHistAcceptanceCuts; // bookkeeping for acceptance cuts + TH1F *fHistClusterIdentificationCuts; // bookkeeping for cluster identification cuts + + TH2F* fHistClusterEtavsPhiBeforeAcc; // eta-phi-distribution before acceptance cuts + TH2F* fHistClusterEtavsPhiAfterAcc; // eta-phi-distribution of all after acceptance cuts + TH2F* fHistClusterEtavsPhiAfterQA; // eta-phi-distribution of all after cluster quality cuts + TH1F* fHistDistanceToBadChannelBeforeAcc; // distance to bad channel before acceptance cuts + TH1F* fHistDistanceToBadChannelAfterAcc; // distance to bad channel after acceptance cuts + TH2F* fHistClusterTimevsEBeforeQA; // Cluster time vs E before cluster quality cuts + TH2F* fHistClusterTimevsEAfterQA; // Cluster time vs E after cluster quality cuts + TH2F* fHistExoticCellBeforeQA; // Exotic cell: 1-Ecross/E cell vs Ecluster before acceptance cuts + TH2F* fHistExoticCellAfterQA; // Exotic cell: 1-Ecross/E cell vs Ecluster after cluster quality cuts + TH1F* fHistDistanceTrackToClusterBeforeQA; // distance cluster to track before acceptance cuts + TH1F* fHistDistanceTrackToClusterAfterQA; // distance cluster to track after cluster quality cuts + TH1F* fHistEnergyOfClusterBeforeQA; // enery per cluster before acceptance cuts + TH1F* fHistEnergyOfClusterAfterQA; // enery per cluster after cluster quality cuts + TH1F* fHistNCellsBeforeQA; // number of cells per cluster before acceptance cuts + TH1F* fHistNCellsAfterQA; // number of cells per cluster after cluster quality cuts + TH1F* fHistM02BeforeQA; // M02 before acceptance cuts + TH1F* fHistM02AfterQA; // M02 after cluster quality cuts + TH1F* fHistM20BeforeQA; // M20 before acceptance cuts + TH1F* fHistM20AfterQA; // M20 after cluster quality cuts + TH1F* fHistDispersionBeforeQA; // dispersion before acceptance cuts + TH1F* fHistDispersionAfterQA; // dispersion after cluster quality cuts + TH1F* fHistNLMBeforeQA; // number of local maxima in cluster before acceptance cuts + TH1F* fHistNLMAfterQA; // number of local maxima in cluster after cluster quality cuts + + private: + + ClassDef(AliCaloPhotonCuts,1) +}; + +#endif diff --git a/PWGGA/GammaConv/macros/AddTask_GammaConvCalo_PbPb.C b/PWGGA/GammaConv/macros/AddTask_GammaConvCalo_PbPb.C new file mode 100644 index 00000000000..6fd0adae669 --- /dev/null +++ b/PWGGA/GammaConv/macros/AddTask_GammaConvCalo_PbPb.C @@ -0,0 +1,205 @@ +void AddTask_GammaConvCalo_PbPb( Int_t trainConfig = 1, //change different set of cuts + Bool_t isMC = kFALSE, //run MC + Int_t enableQAMesonTask = 0, //enable QA in AliAnalysisTaskGammaConvV1 + Int_t enableQAPhotonTask = 0, // enable additional QA task + TString fileNameInputForWeighting = "MCSpectraInput.root", // path to file for weigting input + Int_t headerSelectionInt = 0, // 1 pi0 header, 2 eta header, 3 both (only for "named" boxes) + TString cutnumberAODBranch = "1000000060084000001500000", + TString periodName = "LHC13d2", //name of the period for added signals and weighting + Bool_t doWeighting = kFALSE //enable Weighting + ) { + + // ================= Load Librariers ================================= + gSystem->Load("libCore.so"); + gSystem->Load("libTree.so"); + gSystem->Load("libGeom.so"); + gSystem->Load("libVMC.so"); + gSystem->Load("libPhysics.so"); + gSystem->Load("libMinuit"); + gSystem->Load("libSTEERBase"); + gSystem->Load("libESD"); + gSystem->Load("libAOD"); + gSystem->Load("libANALYSIS"); + gSystem->Load("libANALYSISalice"); + gSystem->Load("libPWGGAGammaConv.so"); + gSystem->Load("libCDB.so"); + gSystem->Load("libSTEER.so"); + gSystem->Load("libSTEERBase.so"); + gSystem->Load("libTENDER.so"); + gSystem->Load("libTENDERSupplies.so"); + + // ================== GetAnalysisManager =============================== + AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager(); + if (!mgr) { + Error(Form("AddTask_GammaConvV1_%i",trainConfig), "No analysis manager found."); + return ; + } + + // ================== GetInputEventHandler ============================= + AliVEventHandler *inputHandler=mgr->GetInputEventHandler(); + + //========= Add PID Reponse to ANALYSIS manager ==== + if(!(AliPIDResponse*)mgr->GetTask("PIDResponseTask")){ + gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskPIDResponse.C"); + AddTaskPIDResponse(isMC,1,0,4,0,"",1,1,4); + } + + //========= Set Cutnumber for V0Reader ================================ + TString cutnumber = "1000000000084001001500000000"; + AliAnalysisDataContainer *cinput = mgr->GetCommonInputContainer(); + + //========= Add V0 Reader to ANALYSIS manager if not yet existent ===== + if( !(AliV0ReaderV1*)mgr->GetTask("V0ReaderV1") ){ + AliV0ReaderV1 *fV0ReaderV1 = new AliV0ReaderV1("V0ReaderV1"); + + fV0ReaderV1->SetUseOwnXYZCalculation(kTRUE); + fV0ReaderV1->SetCreateAODs(kFALSE);// AOD Output + fV0ReaderV1->SetUseAODConversionPhoton(kTRUE); + + if (!mgr) { + Error("AddTask_V0ReaderV1", "No analysis manager found."); + return; + } + + // Set AnalysisCut Number + AliConversionCuts *fCuts=NULL; + if(cutnumber!=""){ + fCuts= new AliConversionCuts(cutnumber.Data(),cutnumber.Data()); + fCuts->SetPreSelectionCutFlag(kTRUE); + if(fCuts->InitializeCutsFromCutString(cutnumber.Data())){ + fV0ReaderV1->SetConversionCuts(fCuts); + fCuts->SetFillCutHistograms("",kTRUE); + } + } + if(inputHandler->IsA()==AliAODInputHandler::Class()){ + // AOD mode + fV0ReaderV1->SetDeltaAODBranchName(Form("GammaConv_%s_gamma",cutnumberAODBranch.Data())); + } + fV0ReaderV1->Init(); + + AliLog::SetGlobalLogLevel(AliLog::kFatal); + + //connect input V0Reader + mgr->AddTask(fV0ReaderV1); + mgr->ConnectInput(fV0ReaderV1,0,cinput); + + } + + //================================================ + //========= Add task to the ANALYSIS manager ===== + //================================================ + AliAnalysisTaskGammaConvCalo *task=NULL; + task= new AliAnalysisTaskGammaConvCalo(Form("GammaConvCalo_%i",trainConfig)); + task->SetIsHeavyIon(1); + task->SetIsMC(isMC); + // Cut Numbers to use in Analysis + Int_t numberOfCuts = 5; + + TString *cutarray = new TString[numberOfCuts]; + TString *clustercutarray = new TString[numberOfCuts]; + TString *mesonCutArray = new TString[numberOfCuts]; + + // meson cuts + // meson type (Dalitz or not), BG scheme, pool depth, rotation degrees, rapidity cut, radius cut, alpha, chi2, shared electrons, reject to close v0, MC smearing, dca, dca, dca + + if (trainConfig == 1){ + cutarray[ 0] = "6010001002092970028250400000"; clustercutarray[0] = "10000040022030000"; mesonCutArray[ 0] = "01525065000000"; // 0-5% + cutarray[ 1] = "6120001002092970028250400000"; clustercutarray[1] = "10000040022030000"; mesonCutArray[ 1] = "01525065000000"; // 5-10% + cutarray[ 2] = "5010001002092970028250400000"; clustercutarray[2] = "10000040022030000"; mesonCutArray[ 2] = "01525065000000"; // 0-10% + cutarray[ 3] = "5240001002092970028250400000"; clustercutarray[3] = "10000040022030000"; mesonCutArray[ 3] = "01525065000000"; // 20-40% + cutarray[ 4] = "5250001002092970028250400000"; clustercutarray[4] = "10000040022030000"; mesonCutArray[ 4] = "01525065000000"; // 20-50% + } else { + Error(Form("GammaConvCalo_%i",trainConfig), "wrong trainConfig variable no cuts have been specified for the configuration"); + return; + } + + TList *ConvCutList = new TList(); + TList *MesonCutList = new TList(); + + TList *HeaderList = new TList(); + if (periodName.CompareTo("LHC13d2")==0){ + TObjString *Header1 = new TObjString("pi0_1"); + HeaderList->Add(Header1); + // TObjString *Header3 = new TObjString("eta_2"); + // HeaderList->Add(Header3); + + } else if (periodName.CompareTo("LHC12a17x_fix")==0){ + TObjString *Header1 = new TObjString("PARAM"); + HeaderList->Add(Header1); + } else if (periodName.CompareTo("LHC14a1a")==0){ + if (headerSelectionInt == 1){ + TObjString *Header1 = new TObjString("pi0_1"); + HeaderList->Add(Header1); + } else if (headerSelectionInt == 2){ + TObjString *Header1 = new TObjString("eta_2"); + HeaderList->Add(Header1); + } else { + TObjString *Header1 = new TObjString("pi0_1"); + HeaderList->Add(Header1); + TObjString *Header2 = new TObjString("eta_2"); + HeaderList->Add(Header2); + } + } else if (periodName.CompareTo("LHC14a1b")==0 || periodName.CompareTo("LHC14a1c")==0){ + TObjString *Header1 = new TObjString("BOX"); + HeaderList->Add(Header1); + } + + ConvCutList->SetOwner(kTRUE); + AliConversionCuts **analysisCuts = new AliConversionCuts*[numberOfCuts]; + ClusterCutList->SetOwner(kTRUE); + AliCaloPhotonCuts **analysisClusterCuts = new AliCaloPhotonCuts*[numberOfCuts]; + MesonCutList->SetOwner(kTRUE); + AliConversionMesonCuts **analysisMesonCuts = new AliConversionMesonCuts*[numberOfCuts]; + + for(Int_t i = 0; iSetUseReweightingWithHistogramFromFile(kTRUE, kTRUE, kFALSE,fileNameInputForWeighting, Form("Pi0_Hijing_%s_PbPb_2760GeV_0005TPC",periodName.Data()), Form("Eta_Hijing_%s_PbPb_2760GeV_0005TPC",periodName.Data()), "","Pi0_Fit_Data_PbPb_2760GeV_0005V0M","Eta_Fit_Data_PbPb_2760GeV_0005V0M"); + if ( i == 1 && doWeighting) analysisCuts[i]->SetUseReweightingWithHistogramFromFile(kTRUE, kTRUE, kFALSE,fileNameInputForWeighting, Form("Pi0_Hijing_%s_PbPb_2760GeV_0510TPC",periodName.Data()), Form("Eta_Hijing_%s_PbPb_2760GeV_0510TPC",periodName.Data()), "","Pi0_Fit_Data_PbPb_2760GeV_0510V0M","Eta_Fit_Data_PbPb_2760GeV_0510V0M"); + if ( i == 2 && doWeighting) analysisCuts[i]->SetUseReweightingWithHistogramFromFile(kTRUE, kTRUE, kFALSE,fileNameInputForWeighting, Form("Pi0_Hijing_%s_PbPb_2760GeV_0010TPC",periodName.Data()), Form("Eta_Hijing_%s_PbPb_2760GeV_0010TPC",periodName.Data()), "","Pi0_Fit_Data_PbPb_2760GeV_0010V0M","Eta_Fit_Data_PbPb_2760GeV_0010V0M"); + if ( i == 3 && doWeighting) analysisCuts[i]->SetUseReweightingWithHistogramFromFile(kTRUE, kTRUE, kFALSE,fileNameInputForWeighting, Form("Pi0_Hijing_%s_PbPb_2760GeV_2040TPC",periodName.Data()), Form("Eta_Hijing_%s_PbPb_2760GeV_2040TPC",periodName.Data()), "","Pi0_Fit_Data_PbPb_2760GeV_2040V0M","Eta_Fit_Data_PbPb_2760GeV_2040V0M"); + if ( i == 4 && doWeighting) analysisCuts[i]->SetUseReweightingWithHistogramFromFile(kTRUE, kTRUE, kFALSE,fileNameInputForWeighting, Form("Pi0_Hijing_%s_PbPb_2760GeV_2050TPC",periodName.Data()), Form("Eta_Hijing_%s_PbPb_2760GeV_2050TPC",periodName.Data()), "","Pi0_Fit_Data_PbPb_2760GeV_2050V0M","Eta_Fit_Data_PbPb_2760GeV_2050V0M"); + } + } + analysisCuts[i]->InitializeCutsFromCutString(cutarray[i].Data()); + if (periodName.CompareTo("LHC14a1b") ==0 || periodName.CompareTo("LHC14a1c") ==0 ){ + if (headerSelectionInt == 1) analysisCuts[i]->SetAddedSignalPDGCode(111); + if (headerSelectionInt == 2) analysisCuts[i]->SetAddedSignalPDGCode(221); + } + ConvCutList->Add(analysisCuts[i]); + + analysisClusterCuts[i] = new AliCaloPhotonCuts(); + analysisClusterCuts[i]->InitializeCutsFromCutString(clustercutarray[i].Data()); + ClusterCutList->Add(analysisClusterCuts[i]); + analysisClusterCuts[i]->SetFillCutHistograms(""); + + analysisCuts[i]->SetFillCutHistograms("",kFALSE); + analysisMesonCuts[i] = new AliConversionMesonCuts(); + analysisMesonCuts[i]->InitializeCutsFromCutString(mesonCutArray[i].Data()); + MesonCutList->Add(analysisMesonCuts[i]); + analysisMesonCuts[i]->SetFillCutHistograms(""); + analysisCuts[i]->SetAcceptedHeader(HeaderList); + } + + task->SetConversionCutList(numberOfCuts,ConvCutList); + task->SetCaloCutList(numberOfCuts,ClusterCutList); + task->SetMesonCutList(numberOfCuts,MesonCutList); + task->SetMoveParticleAccordingToVertex(kTRUE); + task->SetDoMesonAnalysis(kTRUE); + task->SetDoMesonQA(enableQAMesonTask); //Attention new switch for Pi0 QA + task->SetDoPhotonQA(enableQAPhotonTask); //Attention new switch small for Photon QA + + //connect containers + AliAnalysisDataContainer *coutput = + mgr->CreateContainer(Form("GammaConvCalo_%i",trainConfig), TList::Class(), + AliAnalysisManager::kOutputContainer,Form("GammaConvCalo_%i.root",trainConfig)); + + mgr->AddTask(task); + mgr->ConnectInput(task,0,cinput); + mgr->ConnectOutput(task,1,coutput); + + return; + +} diff --git a/PWGGA/GammaConv/macros/AddTask_GammaConvCalo_pPb.C b/PWGGA/GammaConv/macros/AddTask_GammaConvCalo_pPb.C new file mode 100644 index 00000000000..94de2918b69 --- /dev/null +++ b/PWGGA/GammaConv/macros/AddTask_GammaConvCalo_pPb.C @@ -0,0 +1,181 @@ +void AddTask_GammaConvCalo_pPb( Int_t trainConfig = 1, //change different set of cuts + Bool_t isMC = kFALSE, //run MC + Int_t enableQAMesonTask = 0, //enable QA in AliAnalysisTaskGammaConvV1 + Int_t enableQAPhotonTask = 0, // enable additional QA task + TString fileNameInputForWeighting = "MCSpectraInput.root", // path to file for weigting input + Int_t doWeightingPart = 0, //enable Weighting + TString generatorName = "DPMJET", + TString cutnumberAODBranch = "8000000060084000001500000" // cutnumber for AOD branch + ) { + + // ================= Load Librariers ================================= + gSystem->Load("libCore.so"); + gSystem->Load("libTree.so"); + gSystem->Load("libGeom.so"); + gSystem->Load("libVMC.so"); + gSystem->Load("libPhysics.so"); + gSystem->Load("libMinuit"); + gSystem->Load("libSTEERBase"); + gSystem->Load("libESD"); + gSystem->Load("libAOD"); + gSystem->Load("libANALYSIS"); + gSystem->Load("libANALYSISalice"); + gSystem->Load("libPWGGAGammaConv.so"); + gSystem->Load("libCDB.so"); + gSystem->Load("libSTEER.so"); + gSystem->Load("libSTEERBase.so"); + gSystem->Load("libTENDER.so"); + gSystem->Load("libTENDERSupplies.so"); + + // ================== GetAnalysisManager =============================== + AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager(); + if (!mgr) { + Error(Form("AddTask_GammaConvV1_%i",trainConfig), "No analysis manager found."); + return ; + } + + // ================== GetInputEventHandler ============================= + AliVEventHandler *inputHandler=mgr->GetInputEventHandler(); + + //========= Add PID Reponse to ANALYSIS manager ==== + if(!(AliPIDResponse*)mgr->GetTask("PIDResponseTask")){ + gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskPIDResponse.C"); + AddTaskPIDResponse(isMC,1,0,4,0,"",1,1,4); + } + + Printf("here \n"); + + //========= Set Cutnumber for V0Reader ================================ + TString cutnumber = "8000000060084001001500000000"; + AliAnalysisDataContainer *cinput = mgr->GetCommonInputContainer(); + + //========= Add V0 Reader to ANALYSIS manager if not yet existent ===== + if( !(AliV0ReaderV1*)mgr->GetTask("V0ReaderV1") ){ + AliV0ReaderV1 *fV0ReaderV1 = new AliV0ReaderV1("V0ReaderV1"); + + fV0ReaderV1->SetUseOwnXYZCalculation(kTRUE); + fV0ReaderV1->SetCreateAODs(kFALSE);// AOD Output + fV0ReaderV1->SetUseAODConversionPhoton(kTRUE); + + if (!mgr) { + Error("AddTask_V0ReaderV1", "No analysis manager found."); + return; + } + + // Set AnalysisCut Number + AliConversionCuts *fCuts=NULL; + if(cutnumber!=""){ + fCuts= new AliConversionCuts(cutnumber.Data(),cutnumber.Data()); + fCuts->SetPreSelectionCutFlag(kTRUE); + if(fCuts->InitializeCutsFromCutString(cutnumber.Data())){ + fV0ReaderV1->SetConversionCuts(fCuts); + fCuts->SetFillCutHistograms("",kTRUE); + } + } + if(inputHandler->IsA()==AliAODInputHandler::Class()){ + // AOD mode + fV0ReaderV1->SetDeltaAODBranchName(Form("GammaConv_%s_gamma",cutnumberAODBranch.Data())); + } + fV0ReaderV1->Init(); + + AliLog::SetGlobalLogLevel(AliLog::kFatal); + + //connect input V0Reader + mgr->AddTask(fV0ReaderV1); + mgr->ConnectInput(fV0ReaderV1,0,cinput); + + } + + //================================================ + //========= Add task to the ANALYSIS manager ===== + //================================================ + AliAnalysisTaskGammaConvCalo *task=NULL; + task= new AliAnalysisTaskGammaConvCalo(Form("GammaConvCalo_%i",trainConfig)); + task->SetIsHeavyIon(2); + task->SetIsMC(isMC); + // Cut Numbers to use in Analysis + Int_t numberOfCuts = 2; + + TString *cutarray = new TString[numberOfCuts]; + TString *clustercutarray = new TString[numberOfCuts]; + TString *mesonCutArray = new TString[numberOfCuts]; + + // cluster cuts + // 0 "ClusterType", 1 "EtaMin", 2 "EtaMax", 3 "PhiMin", 4 "PhiMax", 5 "DistanceToBadChannel", 6 "Timing", 7 "TrackMatching", 8 "ExoticCell", + // 9 "MinEnergy", 10 "MinNCells", 11 "MinM02", 12 "MaxM02", 13 "MinM20", 14 "MaxM20", 15 "MaximumDispersion", 16 "NLM" + + if (trainConfig == 1){ + cutarray[ 0] = "8000001002092970028250400000"; clustercutarray[0] = "10000040022030000"; mesonCutArray[0] = "01525065000000"; //standart cut, kINT7 + cutarray[ 1] = "8005201002092970028250400000"; clustercutarray[1] = "10000040022030000"; mesonCutArray[1] = "01525065000000"; //standard cut, kEMC7 + } else { + Error(Form("GammaConvCalo_%i",trainConfig), "wrong trainConfig variable no cuts have been specified for the configuration"); + return; + } + + TList *ConvCutList = new TList(); + TList *ClusterCutList = new TList(); + TList *MesonCutList = new TList(); + + TList *HeaderList = new TList(); + if (doWeightingPart==1) { + TObjString *Header1 = new TObjString("pi0_1"); + HeaderList->Add(Header1); + } + if (doWeightingPart==2){ + TObjString *Header3 = new TObjString("eta_2"); + HeaderList->Add(Header3); + } + if (doWeightingPart==3) { + TObjString *Header1 = new TObjString("pi0_1"); + HeaderList->Add(Header1); + TObjString *Header3 = new TObjString("eta_2"); + HeaderList->Add(Header3); + } + + ConvCutList->SetOwner(kTRUE); + AliConversionCuts **analysisCuts = new AliConversionCuts*[numberOfCuts]; + ClusterCutList->SetOwner(kTRUE); + AliCaloPhotonCuts **analysisClusterCuts = new AliCaloPhotonCuts*[numberOfCuts]; + MesonCutList->SetOwner(kTRUE); + AliConversionMesonCuts **analysisMesonCuts = new AliConversionMesonCuts*[numberOfCuts]; + + for(Int_t i = 0; iInitializeCutsFromCutString(cutarray[i].Data()); + ConvCutList->Add(analysisCuts[i]); + analysisCuts[i]->SetFillCutHistograms("",kFALSE); + + Printf("here %d %s\n", i, clustercutarray[i].Data() ); + + analysisClusterCuts[i] = new AliCaloPhotonCuts(); + analysisClusterCuts[i]->InitializeCutsFromCutString(clustercutarray[i].Data()); + ClusterCutList->Add(analysisClusterCuts[i]); + analysisClusterCuts[i]->SetFillCutHistograms(""); + + analysisMesonCuts[i] = new AliConversionMesonCuts(); + analysisMesonCuts[i]->InitializeCutsFromCutString(mesonCutArray[i].Data()); + MesonCutList->Add(analysisMesonCuts[i]); + analysisMesonCuts[i]->SetFillCutHistograms(""); + analysisCuts[i]->SetAcceptedHeader(HeaderList); + } + + task->SetConversionCutList(numberOfCuts,ConvCutList); + task->SetCaloCutList(numberOfCuts,ClusterCutList); + task->SetMesonCutList(numberOfCuts,MesonCutList); + task->SetMoveParticleAccordingToVertex(kTRUE); + task->SetDoMesonAnalysis(kTRUE); + task->SetDoMesonQA(enableQAMesonTask); //Attention new switch for Pi0 QA + task->SetDoPhotonQA(enableQAPhotonTask); //Attention new switch small for Photon QA + + //connect containers + AliAnalysisDataContainer *coutput = + mgr->CreateContainer(Form("GammaConvCalo_%i",trainConfig), TList::Class(), + AliAnalysisManager::kOutputContainer,Form("GammaConvCalo_%i.root",trainConfig)); + + mgr->AddTask(task); + mgr->ConnectInput(task,0,cinput); + mgr->ConnectOutput(task,1,coutput); + + return; + +} diff --git a/PWGGA/GammaConv/macros/AddTask_GammaConvCalo_pp.C b/PWGGA/GammaConv/macros/AddTask_GammaConvCalo_pp.C new file mode 100644 index 00000000000..915b2bc6bae --- /dev/null +++ b/PWGGA/GammaConv/macros/AddTask_GammaConvCalo_pp.C @@ -0,0 +1,162 @@ +void AddTask_GammaConvCalo_pp( Int_t trainConfig = 1, //change different set of cuts + Bool_t isMC = kFALSE, //run MC + Int_t enableQAMesonTask = 1, //enable QA in AliAnalysisTaskGammaConvV1 + Int_t enableQAPhotonTask = 1, // enable additional QA task + TString fileNameInputForWeighting = "MCSpectraInput.root", // path to file for weigting input + TString cutnumberAODBranch = "0000000060084001001500000" + ) { + + // ================= Load Librariers ================================= + gSystem->Load("libCore.so"); + gSystem->Load("libTree.so"); + gSystem->Load("libGeom.so"); + gSystem->Load("libVMC.so"); + gSystem->Load("libPhysics.so"); + gSystem->Load("libMinuit"); + gSystem->Load("libSTEERBase"); + gSystem->Load("libESD"); + gSystem->Load("libAOD"); + gSystem->Load("libANALYSIS"); + gSystem->Load("libANALYSISalice"); + gSystem->Load("libPWGGAGammaConv.so"); + gSystem->Load("libCDB.so"); + gSystem->Load("libSTEER.so"); + gSystem->Load("libSTEERBase.so"); + gSystem->Load("libTENDER.so"); + gSystem->Load("libTENDERSupplies.so"); + + // ================== GetAnalysisManager =============================== + AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager(); + if (!mgr) { + Error(Form("AddTask_GammaConvV1_%i",trainConfig), "No analysis manager found."); + return ; + } + + // ================== GetInputEventHandler ============================= + AliVEventHandler *inputHandler=mgr->GetInputEventHandler(); + + //========= Add PID Reponse to ANALYSIS manager ==== + if(!(AliPIDResponse*)mgr->GetTask("PIDResponseTask")){ + gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskPIDResponse.C"); + AddTaskPIDResponse(isMC,1,0,4,0,"",1,1,4); + } + + //========= Set Cutnumber for V0Reader ================================ + TString cutnumber = "0000000002084000002200000000"; + AliAnalysisDataContainer *cinput = mgr->GetCommonInputContainer(); + + //========= Add V0 Reader to ANALYSIS manager if not yet existent ===== + if( !(AliV0ReaderV1*)mgr->GetTask("V0ReaderV1") ){ + AliV0ReaderV1 *fV0ReaderV1 = new AliV0ReaderV1("V0ReaderV1"); + + fV0ReaderV1->SetUseOwnXYZCalculation(kTRUE); + fV0ReaderV1->SetCreateAODs(kFALSE);// AOD Output + fV0ReaderV1->SetUseAODConversionPhoton(kTRUE); + + if (!mgr) { + Error("AddTask_V0ReaderV1", "No analysis manager found."); + return; + } + + // Set AnalysisCut Number + AliConversionCuts *fCuts=NULL; + if(cutnumber!=""){ + fCuts= new AliConversionCuts(cutnumber.Data(),cutnumber.Data()); + fCuts->SetPreSelectionCutFlag(kTRUE); + if(fCuts->InitializeCutsFromCutString(cutnumber.Data())){ + fV0ReaderV1->SetConversionCuts(fCuts); + fCuts->SetFillCutHistograms("",kTRUE); + } + } + if(inputHandler->IsA()==AliAODInputHandler::Class()){ + // AOD mode + fV0ReaderV1->SetDeltaAODBranchName(Form("GammaConv_%s_gamma",cutnumberAODBranch.Data())); + } + fV0ReaderV1->Init(); + + AliLog::SetGlobalLogLevel(AliLog::kFatal); + + //connect input V0Reader + mgr->AddTask(fV0ReaderV1); + mgr->ConnectInput(fV0ReaderV1,0,cinput); + + } + + //================================================ + //========= Add task to the ANALYSIS manager ===== + //================================================ + AliAnalysisTaskGammaConvCalo *task=NULL; + task= new AliAnalysisTaskGammaConvCalo(Form("GammaConvCalo_%i",trainConfig)); + task->SetIsHeavyIon(0); + task->SetIsMC(isMC); + // Cut Numbers to use in Analysis + Int_t numberOfCuts = 3; + + TString *cutarray = new TString[numberOfCuts]; + TString *clustercutarray = new TString[numberOfCuts]; + TString *mesonCutArray = new TString[numberOfCuts]; + + // meson cuts + // meson type (Dalitz or not), BG scheme, pool depth, rotation degrees, rapidity cut, radius cut, alpha, chi2, shared electrons, reject to close v0, MC smearing, dca, dca, dca + + if (trainConfig == 1){ + cutarray[ 0] = "0000001002092970028250400000"; clustercutarray[0] = "10000040022030000"; mesonCutArray[0] = "01525065000000"; //standard cut LHC11h pp 2.76TeV, kMB + cutarray[ 1] = "0005101002092970028250400000"; clustercutarray[1] = "10000040022030000"; mesonCutArray[1] = "01525065000000"; //standard cut LHC11h pp 2.76TeV, kEMC1 + cutarray[ 2] = "0002001002092970028250400000"; clustercutarray[2] = "10000040022030000"; mesonCutArray[2] = "01525065000000"; //standard cut LHC11h pp 2.76TeV, SDD V0OR + } else { + Error(Form("GammaConvCalo_%i",trainConfig), "wrong trainConfig variable no cuts have been specified for the configuration"); + return; + } + + TList *ConvCutList = new TList(); + TList *MesonCutList = new TList(); + + TList *HeaderList = new TList(); + TObjString *Header1 = new TObjString("BOX"); + HeaderList->Add(Header1); + + ConvCutList->SetOwner(kTRUE); + AliConversionCuts **analysisCuts = new AliConversionCuts*[numberOfCuts]; + ClusterCutList->SetOwner(kTRUE); + AliCaloPhotonCuts **analysisClusterCuts = new AliCaloPhotonCuts*[numberOfCuts]; + MesonCutList->SetOwner(kTRUE); + AliConversionMesonCuts **analysisMesonCuts = new AliConversionMesonCuts*[numberOfCuts]; + + for(Int_t i = 0; iInitializeCutsFromCutString(cutarray[i].Data()); + ConvCutList->Add(analysisCuts[i]); + analysisCuts[i]->SetFillCutHistograms("",kFALSE); + + analysisClusterCuts[i] = new AliCaloPhotonCuts(); + analysisClusterCuts[i]->InitializeCutsFromCutString(clustercutarray[i].Data()); + ClusterCutList->Add(analysisClusterCuts[i]); + analysisClusterCuts[i]->SetFillCutHistograms(""); + + analysisMesonCuts[i] = new AliConversionMesonCuts(); + analysisMesonCuts[i]->InitializeCutsFromCutString(mesonCutArray[i].Data()); + MesonCutList->Add(analysisMesonCuts[i]); + analysisMesonCuts[i]->SetFillCutHistograms(""); + analysisCuts[i]->SetAcceptedHeader(HeaderList); + } + + task->SetConversionCutList(numberOfCuts,ConvCutList); + task->SetCaloCutList(numberOfCuts,ClusterCutList); + task->SetMesonCutList(numberOfCuts,MesonCutList); + task->SetMoveParticleAccordingToVertex(kTRUE); + task->SetDoMesonAnalysis(kTRUE); + task->SetDoMesonQA(enableQAMesonTask); //Attention new switch for Pi0 QA + task->SetDoPhotonQA(enableQAPhotonTask); //Attention new switch small for Photon QA + + //connect containers + AliAnalysisDataContainer *coutput = + mgr->CreateContainer(Form("GammaConvCalo_%i",trainConfig), TList::Class(), + AliAnalysisManager::kOutputContainer,Form("GammaConvCalo_%i.root",trainConfig)); + + mgr->AddTask(task); + mgr->ConnectInput(task,0,cinput); + mgr->ConnectOutput(task,1,coutput); + + return; + +} diff --git a/PWGGA/GammaConv/macros/AddTask_GammaConvV1_pPb.C b/PWGGA/GammaConv/macros/AddTask_GammaConvV1_pPb.C index b964f498945..010619992b5 100644 --- a/PWGGA/GammaConv/macros/AddTask_GammaConvV1_pPb.C +++ b/PWGGA/GammaConv/macros/AddTask_GammaConvV1_pPb.C @@ -992,15 +992,15 @@ void AddTask_GammaConvV1_pPb( Int_t trainConfig = 1, //change different set of cutarray[ 2] = "8000012007092170008260430000"; mesonCutArray[ 2] = "01621035009000"; // min R =35 cm & only photon quality 2 cutarray[ 3] = "8000012007092170008260440000"; mesonCutArray[ 3] = "01621035009000"; // min R =35 cm & only photon quality 3 } else if (trainConfig == 179) { - cutarray[ 0] = "8000011002092170008260400000"; mesonCutArray[ 0] = "01628035009000"; //new standard rapidity 0-0.5 in cms - cutarray[ 1] = "8000011005092170008260400000"; mesonCutArray[ 1] = "01621035009000"; //new standard RCut=10cm - cutarray[ 2] = "8000011008092170008260400000"; mesonCutArray[ 2] = "01621035009000"; //new standard RCut=12.5cm - cutarray[ 3] = "8000011006092170008260400000"; mesonCutArray[ 3] = "01621035009000"; //new standard RCut=20cm + cutarray[ 0] = "8000011002092170008260400000"; mesonCutArray[ 0] = "01628035009000"; //new standard rapidity 0-0.5 in cms + cutarray[ 1] = "8000011005092170008260400000"; mesonCutArray[ 1] = "01621035009000"; //new standard RCut=10cm + cutarray[ 2] = "8000011008092170008260400000"; mesonCutArray[ 2] = "01621035009000"; //new standard RCut=12.5cm + cutarray[ 3] = "8000011006092170008260400000"; mesonCutArray[ 3] = "01621035009000"; //new standard RCut=20cm } else if (trainConfig == 180) { - cutarray[ 0] = "8000012002092170008260400000"; mesonCutArray[ 0] = "01628035009000"; //new standard rapidity 0-0.5 in cms - cutarray[ 1] = "8000012005092170008260400000"; mesonCutArray[ 1] = "01621035009000"; //new standard RCut=10cm - cutarray[ 2] = "8000012008092170008260400000"; mesonCutArray[ 2] = "01621035009000"; //new standard RCut=12.5cm - cutarray[ 3] = "8000012006092170008260400000"; mesonCutArray[ 3] = "01621035009000"; //new standard RCut=20cm + cutarray[ 0] = "8000012002092170008260400000"; mesonCutArray[ 0] = "01628035009000"; //new standard rapidity 0-0.5 in cms + cutarray[ 1] = "8000012005092170008260400000"; mesonCutArray[ 1] = "01621035009000"; //new standard RCut=10cm + cutarray[ 2] = "8000012008092170008260400000"; mesonCutArray[ 2] = "01621035009000"; //new standard RCut=12.5cm + cutarray[ 3] = "8000012006092170008260400000"; mesonCutArray[ 3] = "01621035009000"; //new standard RCut=20cm } else { Error(Form("GammaConvV1_%i",trainConfig), "wrong trainConfig variable no cuts have been specified for the configuration"); return; diff --git a/PWGGA/PWGGAGammaConvLinkDef.h b/PWGGA/PWGGAGammaConvLinkDef.h index 9f52c2e31f7..9e42b410cc3 100644 --- a/PWGGA/PWGGAGammaConvLinkDef.h +++ b/PWGGA/PWGGAGammaConvLinkDef.h @@ -5,17 +5,18 @@ #pragma link off all functions; // Base classes -#pragma link C++ class AliConversionPhotonBase++; +#pragma link C++ class AliConversionPhotonBase+; #pragma link C++ class AliAODConversionParticle+; -#pragma link C++ class AliAODConversionMother++; -#pragma link C++ class AliAODConversionPhoton++; -#pragma link C++ class AliKFConversionPhoton++; -#pragma link C++ class AliKFConversionMother++; -#pragma link C++ class AliConversionCuts++; -#pragma link C++ class AliConversionSelection++; +#pragma link C++ class AliAODConversionMother+; +#pragma link C++ class AliAODConversionPhoton+; +#pragma link C++ class AliKFConversionPhoton+; +#pragma link C++ class AliKFConversionMother+; +#pragma link C++ class AliCaloPhotonCuts+; +#pragma link C++ class AliConversionCuts+; +#pragma link C++ class AliConversionSelection+; #pragma link C++ class AliV0ReaderV1+; #pragma link C++ class AliConversionAODBGHandlerRP+; -#pragma link C++ class AliConversionTrackCuts++; +#pragma link C++ class AliConversionTrackCuts+; #pragma link C++ class AliConversionMesonCuts+; #pragma link C++ class AliDalitzElectronCuts+; #pragma link C++ class AliDalitzElectronSelector+; @@ -30,9 +31,9 @@ #pragma link C++ class AliAnalysisTaskResolution+; #pragma link C++ class AliAnaConvIsolation+; -#pragma link C++ class AliAnaConvCorrBase++; -#pragma link C++ class AliAnaConvCorrPion++; -#pragma link C++ class AliAnaConvCorrPhoton++; +#pragma link C++ class AliAnaConvCorrBase+; +#pragma link C++ class AliAnaConvCorrPion+; +#pragma link C++ class AliAnaConvCorrPhoton+; #pragma link C++ class AliAnalysisTaskdPhi+; // Old tasks @@ -41,6 +42,6 @@ #pragma link C++ class AliPrimaryPionSelector+; #pragma link C++ class AliPrimaryPionCuts+; #pragma link C++ class AliAnalysisTaskEtaToPiPlPiMiGamma+; - +#pragma link C++ class AliAnalysisTaskGammaConvCalo+; #endif -- 2.39.3