+++ /dev/null
-#include "AliV0ReaderV1.h"
-#include "AliKFParticle.h"
-#include "AliAODv0.h"
-#include "AliESDv0.h"
-#include "AliAODEvent.h"
-#include "AliESDEvent.h"
-#include "AliKFParticle.h"
-#include "AliKFConversionPhoton.h"
-#include "AliAODConversionPhoton.h"
-#include "AliConversionPhotonBase.h"
-#include "TVector.h"
-#include "AliKFVertex.h"
-#include "AliAODTrack.h"
-#include "AliESDtrack.h"
-#include "AliAnalysisManager.h"
-#include "AliInputEventHandler.h"
-#include "AliAODHandler.h"
-#include "AliPIDResponse.h"
-#include "TH1.h"
-#include "TH2.h"
-#include "TChain.h"
-#include "AliStack.h"
-
-class iostream;
-
-
-using namespace std;
-
-ClassImp(AliV0ReaderV1)
-
-//________________________________________________________________________
-AliV0ReaderV1::AliV0ReaderV1(const char *name) : AliAnalysisTaskSE(name),
-fConversionGammas(NULL),
-fESDEvent(NULL),
-fAODEvent(NULL),
-fMCStack(NULL),
-fOutputList(NULL),
-fCurrentMotherKFCandidate(NULL),
-fCurrentPositiveKFParticle(NULL),
-fCurrentNegativeKFParticle(NULL),
-fCurrentTrackLabels(NULL),
-fCurrentV0Index(-1),
-fNCentralityBins(5),
- fCentralityBin(0),
- fBGHandler(NULL),
- fVertexZ(-999),
- fCentrality(-1),
- fEPAngle(-1),
-fMaxVertexZ(10),
-fMaxR(180),// 100 meter(outside of ALICE)
-fMinR(5),// 100 meter(outside of ALICE)
-fEtaCut(0.9),
-fEtaCutMin(-0.1),
-fPtCut(0.),
-fSinglePtCut(0.),
-fMaxZ(240),
-fMinClsTPC(0.),
-fMinClsTPCToF(0.),
-fLineCutZRSlope(0.),
-fLineCutZValue(7),
-fLineCutZRSlopeMin(0.),
-fLineCutZValueMin(-2),
-fChi2CutConversion(30),
-fPIDProbabilityCutNegativeParticle(0),
-fPIDProbabilityCutPositiveParticle(0),
-fDodEdxSigmaCut(kTRUE),
-fDoTOFsigmaCut(kFALSE), // RRnewTOF
-fPIDTRDEfficiency(0.95),
-fDoTRDPID(kFALSE),
-fPIDnSigmaAboveElectronLine(5),
-fPIDnSigmaBelowElectronLine(-3),
-fTofPIDnSigmaAboveElectronLine(100), // RRnewTOF
-fTofPIDnSigmaBelowElectronLine(-100), // RRnewTOF
-fPIDnSigmaAbovePionLine(0),
-fPIDnSigmaAbovePionLineHighPt(-100),
-fPIDMinPnSigmaAbovePionLine(1),
-fPIDMaxPnSigmaAbovePionLine(3),
-fDoKaonRejectionLowP(kTRUE),
- fDoProtonRejectionLowP(kTRUE),
- fDoPionRejectionLowP(kTRUE),
- fPIDnSigmaAtLowPAroundKaonLine(0),
- fPIDnSigmaAtLowPAroundProtonLine(0),
- fPIDnSigmaAtLowPAroundPionLine(0),
- fPIDMinPKaonRejectionLowP(1.5),
- fPIDMinPProtonRejectionLowP(2),
- fPIDMinPPionRejectionLowP(0.5),
- fDoQtGammaSelection(kTRUE),
- fDoHighPtQtGammaSelection(kTRUE), // RRnew
- fQtMax(0.05),
- fHighPtQtMax(0.06), // RRnew
- fPtBorderForQt(2.5), // RRnew
- fNSigmaMass(0.),
- fUseImprovedVertex(kTRUE),
- fUseOwnXYZCalculation(kTRUE),
- fUseConstructGamma(kFALSE),
- fUseEtaMinCut(kFALSE),
- fUseOnFlyV0Finder(kTRUE),
- fDoPhotonAsymmetryCut(kTRUE),
- fMinPPhotonAsymmetryCut(100.),
-fMinPhotonAsymmetry(0.),
- kUseAODConversionPhoton(kFALSE),
- fIsHeavyIon(kTRUE),
- fCreateAOD(kFALSE),
- fDeltaAODFilename("AliAODGammaConversion.root")
-{
-
- fLineCutZRSlope=tan(2*atan(exp(-fEtaCut)));
- // Input slot #0 works with a TChain
- DefineInput(0, TChain::Class());
- // Output slot #0 id reserved by the base class for AOD
- // Output slot #1 writes into a TH1 container
- DefineOutput(1, TList::Class());
-
- fCurrentTrackLabels=new Int_t[2];
-}
-
-
-//________________________________________________________________________
-AliV0ReaderV1::~AliV0ReaderV1()
-{
- if(fConversionGammas){
- fConversionGammas->Delete();// Clear Objects
- delete fConversionGammas;
- fConversionGammas=0x0;
- }
-
- if(fCurrentTrackLabels){
- delete[] fCurrentTrackLabels;
- fCurrentTrackLabels=NULL;}
-
-}
-
-//________________________________________________________________________
-void AliV0ReaderV1::UserCreateOutputObjects()
-{
-
-if(fOutputList != NULL){
- delete fOutputList;
- fOutputList = NULL;
- }
- if(fOutputList == NULL){
- fOutputList = new TList();
- fOutputList->SetOwner(kTRUE);
- }
-
- TList *fCutList=new TList();
- fCutList->SetName("GammaReconstruction");
- fCutList->SetOwner(kTRUE);
- fOutputList->Add(fCutList);
-
-//GammaMass-plots
-Int_t kGCnXBinsGammaMass = 4000;
-Double_t kGCfirstXBinGammaMass= 0.;
-Double_t kGClastXBinGammaMass = 1.;
-
- Int_t kGCnYBinsSpectra = 250;
- Double_t kGCfirstYBinSpectra = 0.;
- Double_t kGClastYBinSpectra = 25.;
-
-// Process Gammas Histograms
-
-hV0CurrentFinder=new TH1F("ESD_V0sCurrentFinder_InvMass","V0sCurrentFinder",kGCnXBinsGammaMass,kGCfirstXBinGammaMass,kGClastXBinGammaMass);
-fCutList->Add(hV0CurrentFinder);
-hV0AllArmenteros=new TH2F("ESD_V0sCurrentFinder_Armenteros","Armenteros Alpha Qt",200,-1,1,250,0,0.25);
-fCutList->Add(hV0AllArmenteros);
-hV0Good=new TH1F("ESD_GoodV0s_InvMass","GoodV0s",kGCnXBinsGammaMass,kGCfirstXBinGammaMass,kGClastXBinGammaMass);
-fCutList->Add(hV0Good);
-hV0GoodArmenteros=new TH2F("ESD_GoodV0s_Armenteros","Armenteros Alpha Qt",200,-1,1,250,0,0.25);
-fCutList->Add(hV0GoodArmenteros);
-
-
-// Track Cuts
-hV0CutLikeSign=new TH1F("ESD_CutLikeSign_InvMass","LikeSign",kGCnXBinsGammaMass,kGCfirstXBinGammaMass,kGClastXBinGammaMass);
-fCutList->Add(hV0CutLikeSign);
-hV0CutRefit=new TH1F("ESD_CutRefit_InvMass","No TPC refit",kGCnXBinsGammaMass,kGCfirstXBinGammaMass,kGClastXBinGammaMass);
-fCutList->Add(hV0CutRefit);
-hV0CutKinks=new TH1F("ESD_CutKink_InvMass","Kinks",kGCnXBinsGammaMass,kGCfirstXBinGammaMass,kGClastXBinGammaMass);
-fCutList->Add(hV0CutKinks);
-hV0CutMinNclsTPCToF=new TH1F("ESD_CutMinNClsTPCToF_InvMass","Min Ncls TPC ToF",kGCnXBinsGammaMass,kGCfirstXBinGammaMass,kGClastXBinGammaMass);
-fCutList->Add(hV0CutMinNclsTPCToF);
-
-// Event Cuts
-/*hV0CutNContributors=new TH1F("ESD_CutNContributors_InvMass","NContributors<=0",kGCnXBinsGammaMass,kGCfirstXBinGammaMass,kGClastXBinGammaMass);
-fCutList->Add(hV0CutNContributors);
-hV0CutVertexZ=new TH1F("ESD_CutVertexZ_InvMass","VertexZ",kGCnXBinsGammaMass,kGCfirstXBinGammaMass,kGClastXBinGammaMass);
-fCutList->Add(hV0CutVertexZ);
-*/
-
-
-// dEdx Cuts
-
-hV0CutdEdxElectron=new TH1F("ESD_CutdEdxSigmaElectronLine_InvMass" ,"dedx ElectronLine" , kGCnXBinsGammaMass, kGCfirstXBinGammaMass, kGClastXBinGammaMass);
-fCutList->Add(hV0CutdEdxElectron);
-hV0CutdEdxPion=new TH1F("ESD_CutdEdxSigmaPionLine_InvMass" ,"dedx PionLine" , kGCnXBinsGammaMass, kGCfirstXBinGammaMass, kGClastXBinGammaMass);
-fCutList->Add(hV0CutdEdxPion);
-hV0CutdEdxKaonLowP=new TH1F("ESD_CutKaonRejectionLowP_InvMass" ,"dedx KaonRejection LowP" , kGCnXBinsGammaMass, kGCfirstXBinGammaMass, kGClastXBinGammaMass);
-fCutList->Add(hV0CutdEdxKaonLowP);
-hV0CutdEdxProtonLowP=new TH1F("ESD_CutProtonRejectionLowP_InvMass" ,"dedx ProtonRejection LowP" , kGCnXBinsGammaMass, kGCfirstXBinGammaMass, kGClastXBinGammaMass);
-fCutList->Add(hV0CutdEdxProtonLowP);
-hV0CutdEdxPionLowP=new TH1F("ESD_CutPionRejectionLowP_InvMass" ,"dedx PionRejection LowP" , kGCnXBinsGammaMass, kGCfirstXBinGammaMass, kGClastXBinGammaMass);
-fCutList->Add(hV0CutdEdxPionLowP);
-hV0CutdEdxTOFElectron=new TH1F("ESD_CutTOFsigmaElec_InvMass", "ESD_CutTOFsigmaElec_InvMass",kGCnXBinsGammaMass, kGCfirstXBinGammaMass, kGClastXBinGammaMass);
-fCutList->Add(hV0CutdEdxTOFElectron);
-hV0CutdEdxTRD=new TH1F("ESD_CutTRD_InvMass", "ESD_CutTRD_InvMass",kGCnXBinsGammaMass, kGCfirstXBinGammaMass, kGClastXBinGammaMass);
-fCutList->Add(hV0CutdEdxTRD);
-hGammadEdxbefore=new TH2F("Gamma_dEdx_before","dEdx Gamma before" ,kGCnYBinsSpectra, kGCfirstYBinSpectra, kGClastYBinSpectra,400, 0,200);
-fCutList->Add(hGammadEdxbefore);
-hGammadEdxafter=new TH2F("Gamma_dEdx_after","dEdx Gamma after" ,kGCnYBinsSpectra, kGCfirstYBinSpectra, kGClastYBinSpectra,400, 0,200);
-fCutList->Add(hGammadEdxafter);
-
-
-// Armenteros
-
-hV0CutQt=new TH1F("ESD_CutQt_InvMass","ESD_CutQt_InvMass",kGCnXBinsGammaMass, kGCfirstXBinGammaMass, kGClastXBinGammaMass);
-fCutList->Add(hV0CutQt);
-
-// Kinematic Cuts
-
-hV0CutR=new TH1F("ESD_CutR_InvMass" ,"Above RMax" , kGCnXBinsGammaMass, kGCfirstXBinGammaMass, kGClastXBinGammaMass);
-fCutList->Add(hV0CutR);
-hV0CutMinR=new TH1F("ESD_CutMinR_InvMass" ,"Above RMax" , kGCnXBinsGammaMass, kGCfirstXBinGammaMass, kGClastXBinGammaMass);
-fCutList->Add(hV0CutMinR);
-hV0CutLine=new TH1F("ESD_CutLine_InvMass" ,"Out of reconstruction area" , kGCnXBinsGammaMass, kGCfirstXBinGammaMass, kGClastXBinGammaMass);
-fCutList->Add(hV0CutLine);
-hV0CutZ=new TH1F("ESD_CutZ_InvMass" ,"Out of reconstruction area" , kGCnXBinsGammaMass, kGCfirstXBinGammaMass, kGClastXBinGammaMass);
-fCutList->Add(hV0CutZ);
-hV0CutEta=new TH1F("ESD_CutEta_InvMass" ,"Above #eta max" , kGCnXBinsGammaMass, kGCfirstXBinGammaMass, kGClastXBinGammaMass);
-fCutList->Add(hV0CutEta);
-
-
-hV0CutSinglePt=new TH1F("ESD_CutSinglePt_InvMass" ,"Below p_{t} min" , kGCnXBinsGammaMass, kGCfirstXBinGammaMass, kGClastXBinGammaMass);
-fCutList->Add(hV0CutSinglePt);
-hV0CutNDF=new TH1F("ESD_CutNDF_InvMass" ,"#chi^{2} > Max" , kGCnXBinsGammaMass, kGCfirstXBinGammaMass, kGClastXBinGammaMass);
-fCutList->Add(hV0CutNDF);
-hV0CutChi2=new TH1F("ESD_CutChi2_InvMass" ,"#chi^{2} > Max" , kGCnXBinsGammaMass, kGCfirstXBinGammaMass, kGClastXBinGammaMass);
-fCutList->Add(hV0CutChi2);
-hV0CutPt= new TH1F("ESD_CutPt_InvMass" ,"Below p_{t} min" , kGCnXBinsGammaMass, kGCfirstXBinGammaMass, kGClastXBinGammaMass);
-fCutList->Add(hV0CutPt);
-
-// Asymmetry Cut
-
-hV0CutAsymmetry=new TH1F("ESD_CutPhotonAsymmetry_InvMass" ,"Out of reconstruction area" , kGCnXBinsGammaMass, kGCfirstXBinGammaMass, kGClastXBinGammaMass);
-fCutList->Add(hV0CutAsymmetry);
-
-// PID Prob
-
-hV0CutPIDProb=new TH1F("ESD_CutPIDProb_InvMass" ,"wrong TPC PID" , kGCnXBinsGammaMass, kGCfirstXBinGammaMass, kGClastXBinGammaMass);
-fCutList->Add(hV0CutPIDProb);
-
-// Event Info
-
- // Other
- TList *fOtherList=new TList();
- fOtherList->SetName("EventInfo");
- fOtherList->SetOwner(kTRUE);
- fOutputList->Add(fOtherList);
-
-
- hV0EventCuts=new TH1F("ESD_EventCuts","Event Cuts",10,-0.5,9.5);
- fOtherList->Add(hV0EventCuts);
- hNEvents=new TH1F("NEvents_vs_Centrality","NEvents vs Centrality",fNCentralityBins,-0.5,fNCentralityBins-0.5);
- fOtherList->Add(hNEvents);
- hCentrality=new TH1F("Centrality","Centrality",100,0,100);
- fOtherList->Add(hCentrality);
- hVertexZ=new TH1F("VertexZ","VertexZ",1000,-50,50);
- fOtherList->Add(hVertexZ);
-
-// QA
-
- TList *fGammaList=new TList();
- fGammaList->SetName("GammaInfo");
- fGammaList->SetOwner(kTRUE);
- fOutputList->Add(fGammaList);
-
- hGammaPt=new TH1F*[fNCentralityBins];
-
- for(int i=0;i<fNCentralityBins;i++){
- hGammaPt[i]=new TH1F(Form("GammaSpectrum_RECO_Pt_%d",i),"Reco Gamma Pt",kGCnYBinsSpectra, kGCfirstYBinSpectra, kGClastYBinSpectra);
- fGammaList->Add(hGammaPt[i]);
- }
-
- hGammaPhi=new TH1F("Phi_Gamma","Phi Gamma" ,36, 0, 2*TMath::Pi());
- fGammaList->Add(hGammaPhi);
- hGammaConversionMapXY=new TH2F("Gamma_ConversionMap_XY","Conversion Point xy",400,-200,200,400,-200,200);
- fGammaList->Add(hGammaConversionMapXY);
- hGammaConversionMapZR=new TH2F("Gamma_ConversionMap_ZR","Conversion Point zr",500,-250,250,180,0,180);
- fGammaList->Add(hGammaConversionMapZR);
-
-
- // FILL MC PART only if MC is available
- if(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()){
-
- hMCPtTRUE=new TH1F*[fNCentralityBins];
- hMCPtRECOTRUE=new TH1F*[fNCentralityBins];
-
- for(int i=0;i<fNCentralityBins;i++){
- hMCPtTRUE[i]=new TH1F(Form("GammaSpectrum_MC_Pt_%d",i),"TRUE Gamma Pt",kGCnYBinsSpectra, kGCfirstYBinSpectra, kGClastYBinSpectra);
- fGammaList->Add(hMCPtTRUE[i]);
- hMCPtRECOTRUE[i]=new TH1F(Form("GammaSpectrum_RECOTRUE_Pt_%d",i),"True Reco Gamma Pt",kGCnYBinsSpectra, kGCfirstYBinSpectra, kGClastYBinSpectra);
- fGammaList->Add(hMCPtRECOTRUE[i]);
- }
-
- // Call Sumw2 Option
- for (Int_t i=0; i<fGammaList->GetEntries(); i++) {
- TH1 *h1 = dynamic_cast<TH1*>(fGammaList->At(i));
- if (h1){h1->Sumw2();}
- }
-
- hMCPtResolution=new TH2F("Resolution_Gamma_dPt_Pt","dPt vs Pt", kGCnYBinsSpectra, kGCfirstYBinSpectra, kGClastYBinSpectra,200,-10,10);
- fGammaList->Add(hMCPtResolution);
- hMCPtResolutionPhi=new TH2F("Resolution_Gamma_dPt_Phi","dPt vs Phi",180,0,2*TMath::Pi(),200,-10,10);
- fGammaList->Add(hMCPtResolutionPhi);
- hMCRResolutionvsR=new TH2F("Resolution_dRAbs_VS_R","dR vs R", 720,0,360,100,-5,5);
- fGammaList->Add(hMCRResolutionvsR);
- hMCZResolutionvsZ=new TH2F("Resolution_dZAbs_VS_Z","dZ vs Z", 200,-50,50,100,-5,5);
- fGammaList->Add(hMCZResolutionvsZ);
-
-}
-// Gamma Output
-
-if(fCreateAOD){kUseAODConversionPhoton=kTRUE;}
-
-if(fConversionGammas == NULL){
- if(kUseAODConversionPhoton){
- fConversionGammas = new TClonesArray("AliAODConversionPhoton",100);}
- else{
- fConversionGammas = new TClonesArray("AliKFConversionPhoton",100);}
-}
-fConversionGammas->Delete();//Reset the TClonesArray
-
-// Create AODs
-
-if(fCreateAOD){
- fConversionGammas->SetName(Form("GammaConv_gamma"));
-
- AddAODBranch("TClonesArray", &fConversionGammas, fDeltaAODFilename.Data());
- AliAnalysisManager::GetAnalysisManager()->RegisterExtraFile(fDeltaAODFilename.Data());
-}
-
- PostData(1, fOutputList);
-
-}
-//________________________________________________________________________
-void AliV0ReaderV1::UserExec(Option_t *){
-
- fAODEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
- fESDEvent = dynamic_cast<AliESDEvent*>(fInputEvent);
- if(!fAODEvent&&!fESDEvent) {
- AliError("No Input event");
- return;
- }
-
- fMCStack=NULL;
- if(fMCEvent){
- fMCStack = fMCEvent->Stack();}
-
-
- if(fESDEvent)AliKFParticle::SetField(fESDEvent->GetMagneticField());
- if(fAODEvent)AliKFParticle::SetField(fAODEvent->GetMagneticField());
-
- fConversionGammas->Delete();//Reset the TClonesArray
-
- // Event Cuts
-
- EventCuts();
-
- if(EventIsSelected()){
-
- // Process V0s
- for(fCurrentV0Index=0;fCurrentV0Index<fInputEvent->GetNumberOfV0s();fCurrentV0Index++){
- if(CheckV0Status()){
-
- ProcessV0();
- }
- }
- ProcessMCGammasForEfficiency();
-
- // Set AOD Output
-
- ///Make sure delta aod is filled if standard aod is filled (for synchronization when reading aod with standard aod)
- if(fCreateAOD) {
- AliAODHandler * aodhandler = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
- if (aodhandler && aodhandler->GetFillAOD()) {
- AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillExtension(kTRUE);
- }}
- }
-
- PostData(1, fOutputList);
-}
- ///________________________________________________________________________
-Bool_t AliV0ReaderV1::IsGammaCandidate(AliConversionPhotonBase *fPhotonCandidate)
-{
-
- // Fill Histos before Cuts
- hV0CurrentFinder->Fill(fPhotonCandidate->GetPhotonMass());
- hV0AllArmenteros->Fill(fPhotonCandidate->GetArmenterosAlpha(),fPhotonCandidate->GetArmenterosQt());
-
- Bool_t passcuts=kTRUE;
-
- // Gamma selection based on QT from Armenteros
- if(fDoQtGammaSelection == kTRUE){
- if(!ArmenterosQtCut(fPhotonCandidate))return kFALSE;//passcuts=kFALSE;
- }
-
- // Chi Cut
-
- if(fPhotonCandidate->GetChi2perNDF() > fChi2CutConversion || fPhotonCandidate->GetChi2perNDF() <=0){
- hV0CutChi2->Fill(fPhotonCandidate->GetPhotonMass());
- return kFALSE;
- }
-
- // Reconstruction Acceptance Cuts
- if(!AcceptanceCuts(fPhotonCandidate))return kFALSE;//passcuts=kFALSE;
-
-
- // Track Cuts
- if(!TrackCuts(fPhotonCandidate))return kFALSE;//passcuts=kFALSE;
-
-
- // PID Cuts
- if(!dEdxCuts(fPhotonCandidate))return kFALSE;//passcuts=kFALSE;
-
- // Asymmetry Cut
- if(fDoPhotonAsymmetryCut == kTRUE){
- if(!AsymmetryCut(fPhotonCandidate))return kFALSE;//passcuts=kFALSE;
- }
-
- //Check the pid probability
-
- if(!PIDProbabilityCut(fPhotonCandidate))return kFALSE;//passcuts=kFALSE;
-
-return passcuts;
-
-}
-
-///________________________________________________________________________
-const AliExternalTrackParam *AliV0ReaderV1::GetExternalTrackParam(Int_t charge){
-
-
- if(!(charge==1||charge==-1)){AliError("Charge not defined");return 0x0;}
-
- Int_t label;
- if(charge>0)label=0;
- else label=1;
- // Check for sign flip
-
- if(fESDEvent){
- AliESDv0 *fCurrentV0=dynamic_cast<AliESDv0*>(fESDEvent->GetV0(fCurrentV0Index));
- if(fCurrentV0){
- if(!fCurrentV0->GetParamN()||!fCurrentV0->GetParamP())return 0x0;
- if(!fESDEvent->GetTrack(fCurrentV0->GetNindex())||!fESDEvent->GetTrack(fCurrentV0->GetPindex()))return 0x0;
- if((fESDEvent->GetTrack(fCurrentV0->GetPindex()))->Charge()==charge){
- fCurrentTrackLabels[label]=fCurrentV0->GetPindex();
- return fCurrentV0->GetParamP();}
- if((fESDEvent->GetTrack(fCurrentV0->GetNindex()))->Charge()==charge){
- fCurrentTrackLabels[label]=fCurrentV0->GetNindex();
- return fCurrentV0->GetParamN();}
- }
-
- }
-
- if(fAODEvent){
-/* AliAODv0 *fCurrentV0=dynamic_cast<AliAODv0*>(fAODEvent->GetV0(fCurrentV0Index));
- if(fCurrentV0){
- if(!fCurrentV0->GetParamN()||!fCurrentV0->GetParamP())return 0x0;
- if(!fAODEvent->GetTrack(fCurrentV0->GetNegID())||!fAODEvent->GetTrack(fCurrentV0->GetPosID()))return 0x0;
- if((fAODEvent->GetTrack(fCurrentV0->GetPosID()))->Charge()==1){
- fCurrentTrackLabels[label]=fCurrentV0->GetPosID();
- return fCurrentV0->GetParamP();}
- if((fAODEvent->GetTrack(fCurrentV0->GetNegID()))->Charge()==1){
- fCurrentTrackLabels[label]=fCurrentV0->GetNegID();
- return fCurrentV0->GetParamN();}
- }
- */ }
- return 0x0;
-}
-
-///________________________________________________________________________
-Bool_t AliV0ReaderV1::CheckV0Status()
-{
- if(fESDEvent){
- AliESDv0 *fCurrentV0=(AliESDv0*)(fESDEvent->GetV0(fCurrentV0Index));
- if(!fCurrentV0){
- printf("Requested V0 does not exist");
- return kFALSE;}
- //checks if on the fly mode is set
- if(fCurrentV0->GetOnFlyStatus()==fUseOnFlyV0Finder)return kTRUE;
- }
-
- if(fAODEvent){
- AliAODv0 *fCurrentV0=dynamic_cast<AliAODv0*>(fAODEvent->GetV0(fCurrentV0Index));
- if(!fCurrentV0){
- AliWarning("Requested V0 does not exist");
- return kFALSE;}
-
- //checks if on the fly mode is set
- if(fCurrentV0->GetOnFlyStatus()==fUseOnFlyV0Finder)return kTRUE;
- }
- return kFALSE;
-}
-
-
-
-
-///________________________________________________________________________
-void AliV0ReaderV1::ProcessV0(){
-
- // Reset TrackLabels
- fCurrentTrackLabels[0]=-1;
- fCurrentTrackLabels[1]=-1;
-
-
- // cout<<"V0ReaderV1 ProcessV0 "<<fCurrentV0Index<<endl;
-
- // Get Daughter KF Particles
-
- const AliExternalTrackParam *fCurrentExternalTrackParamPositive=GetExternalTrackParamP();
- const AliExternalTrackParam *fCurrentExternalTrackParamNegative=GetExternalTrackParamN();
-
- if(fCurrentExternalTrackParamPositive&&fCurrentExternalTrackParamNegative){
-
- fCurrentNegativeKFParticle=new AliKFParticle(*(fCurrentExternalTrackParamNegative),11);
- fCurrentPositiveKFParticle=new AliKFParticle(*(fCurrentExternalTrackParamPositive),-11);
-
- }
- //else{hV0CutLikeSign->Fill(fCurrentMotherKFCandidate->M());}// Like Sign error is already catched here
-
-
- // Reconstruct Gamma
-
- if(fCurrentNegativeKFParticle&&fCurrentPositiveKFParticle){
-
- if(fUseConstructGamma==kTRUE){
-
- fCurrentMotherKFCandidate = new AliKFConversionPhoton();
- fCurrentMotherKFCandidate->ConstructGamma(*fCurrentNegativeKFParticle,*fCurrentPositiveKFParticle);
- }else{
- fCurrentMotherKFCandidate = new AliKFConversionPhoton(*fCurrentNegativeKFParticle,*fCurrentPositiveKFParticle);
- fCurrentMotherKFCandidate->SetMassConstraint(0,fNSigmaMass);
- }
-
- if(fCurrentNegativeKFParticle){delete fCurrentNegativeKFParticle;
- fCurrentNegativeKFParticle=0x0;}
- if(fCurrentPositiveKFParticle){ delete fCurrentPositiveKFParticle;
- fCurrentPositiveKFParticle=0x0;}
-
-
- // Update Vertex
- if(fUseImprovedVertex == kTRUE){
- AliKFVertex primaryVertexImproved(*GetPrimaryVertex());
- primaryVertexImproved+=*fCurrentMotherKFCandidate;
- fCurrentMotherKFCandidate->SetProductionVertex(primaryVertexImproved);
- }
-
- // Set Track Labels
-
- fCurrentMotherKFCandidate->SetV0Index(fCurrentV0Index);
- fCurrentMotherKFCandidate->SetTrackLabels(fCurrentTrackLabels[0],fCurrentTrackLabels[1]);
-
- //Set MC Label
-
- if(fMCStack){
- Int_t labeln=TMath::Abs(GetTrack(fCurrentMotherKFCandidate->GetTrackLabelPositive())->GetLabel());
- Int_t labelp=TMath::Abs(GetTrack(fCurrentMotherKFCandidate->GetTrackLabelNegative())->GetLabel());
-
- TParticle *fNegativeMCParticle = fMCStack->Particle(labeln);
- TParticle *fPositiveMCParticle = fMCStack->Particle(labelp);
-
- if(fPositiveMCParticle&&fNegativeMCParticle){
- fCurrentMotherKFCandidate->SetMCLabelPositive(labelp);
- fCurrentMotherKFCandidate->SetMCLabelNegative(labeln);
- }
- }
-
-
-
- //Add PID information with ESD tender (AOD implementation is not complete)
-
-
- AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
- AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
- AliPIDResponse *fPIDResponse = (AliPIDResponse*)inputHandler->GetPIDResponse();
-
- if(fESDEvent){
- Int_t labelp=((AliESDv0*)fESDEvent->GetV0(fCurrentMotherKFCandidate->GetV0Index()))->GetNindex();
- Int_t labeln=((AliESDv0*)fESDEvent->GetV0(fCurrentMotherKFCandidate->GetV0Index()))->GetNindex();
-
- AliESDtrack *trackpos=fESDEvent->GetTrack(labelp);
- AliESDtrack *trackneg=fESDEvent->GetTrack(labeln);
-
- if(trackpos&&trackneg){
-
- Float_t fNSigmadEdxPositive[5];
- Float_t fNSigmadEdxNegative[5];
-
- fNSigmadEdxPositive[0]=fPIDResponse->NumberOfSigmasTPC(trackpos,AliPID::kElectron);
- fNSigmadEdxPositive[1]=fPIDResponse->NumberOfSigmasTPC(trackpos,AliPID::kMuon);
- fNSigmadEdxPositive[2]=fPIDResponse->NumberOfSigmasTPC(trackpos,AliPID::kPion);
- fNSigmadEdxPositive[3]=fPIDResponse->NumberOfSigmasTPC(trackpos,AliPID::kKaon);
- fNSigmadEdxPositive[4]=fPIDResponse->NumberOfSigmasTPC(trackpos,AliPID::kProton);
-
- fNSigmadEdxNegative[0]=fPIDResponse->NumberOfSigmasTPC(trackneg,AliPID::kElectron);
- fNSigmadEdxNegative[1]=fPIDResponse->NumberOfSigmasTPC(trackneg,AliPID::kMuon);
- fNSigmadEdxNegative[2]=fPIDResponse->NumberOfSigmasTPC(trackneg,AliPID::kPion);
- fNSigmadEdxNegative[3]=fPIDResponse->NumberOfSigmasTPC(trackneg,AliPID::kKaon);
- fNSigmadEdxNegative[4]=fPIDResponse->NumberOfSigmasTPC(trackneg,AliPID::kProton);
-
- fCurrentMotherKFCandidate->SetNSigmadEdx(fNSigmadEdxPositive,fNSigmadEdxNegative);
- }
- }
-
- // Calculate ConversionPoint
-
- if(fUseOwnXYZCalculation){
-
- Double_t convpos[3]={0,0,0};
- GetConversionPoint(fCurrentExternalTrackParamPositive,fCurrentExternalTrackParamNegative,convpos);
- fCurrentMotherKFCandidate->SetConversionPoint(convpos);
-
- }
-
- // if(kTRUE){
- if(IsGammaCandidate(fCurrentMotherKFCandidate)){
- // Fill Histos after Cuts
-
- // Process MC
-
- ProcessMC(fCurrentMotherKFCandidate);
-
-
- hV0Good->Fill(fCurrentMotherKFCandidate->M());
- hV0GoodArmenteros->Fill(fCurrentMotherKFCandidate->GetArmenterosAlpha(),fCurrentMotherKFCandidate->GetArmenterosQt());
-
-
- // Set Mass Zero for Gammas
- SetGammaMassZero();
-
- // Add Gamma to the TClonesArray
-
- if(kUseAODConversionPhoton){
- new((*fConversionGammas)[fConversionGammas->GetEntriesFast()]) AliAODConversionPhoton(fCurrentMotherKFCandidate);
- }
- else{
- new((*fConversionGammas)[fConversionGammas->GetEntriesFast()]) AliKFConversionPhoton(*fCurrentMotherKFCandidate);
- }
- // Fill QA Histos
-
- hGammaPhi->Fill(fCurrentMotherKFCandidate->Phi());
- hGammaPt[fCentralityBin]->Fill(fCurrentMotherKFCandidate->Pt());
- hGammaConversionMapXY->Fill(fCurrentMotherKFCandidate->GetConversionX(),fCurrentMotherKFCandidate->GetConversionY());
- hGammaConversionMapZR->Fill(fCurrentMotherKFCandidate->GetConversionZ(),fCurrentMotherKFCandidate->GetConversionRadius());
-
- }
- delete fCurrentMotherKFCandidate;
- fCurrentMotherKFCandidate=NULL;
- }
-
-}
-
-///________________________________________________________________________
-Bool_t AliV0ReaderV1::ArmenterosQtCut(AliConversionPhotonBase *fPhotonCandidate)
-{
- if(fDoHighPtQtGammaSelection){
- if(fPhotonCandidate->GetPhotonPt() < fPtBorderForQt){
- if(fPhotonCandidate->GetArmenterosQt()>fQtMax){
- hV0CutQt->Fill(fPhotonCandidate->GetPhotonMass());
- return kFALSE;
- }
- } else {
- if(fPhotonCandidate->GetArmenterosQt()>fHighPtQtMax){
- hV0CutQt->Fill(fPhotonCandidate->GetPhotonMass());
- return kFALSE;
- }
- }
- } else {
-
- if(fPhotonCandidate->GetArmenterosQt()>fQtMax){
- hV0CutQt->Fill(fPhotonCandidate->GetPhotonMass());
- return kFALSE;
- }
- }
- return kTRUE;
-}
-
-///________________________________________________________________________
-Bool_t AliV0ReaderV1::AcceptanceCuts(AliConversionPhotonBase *fPhotonCandidate)
-{
- AliVTrack *fCurrentNegativeTrack=GetTrack(fPhotonCandidate->GetTrackLabelNegative());
- AliVTrack *fCurrentPositiveTrack=GetTrack(fPhotonCandidate->GetTrackLabelPositive());
-
- if(fPhotonCandidate->GetConversionRadius()>fMaxR){ // cuts on distance from collision point
- hV0CutR->Fill(fPhotonCandidate->GetPhotonMass());
- return kFALSE;
- }
-
- if(fPhotonCandidate->GetConversionRadius()<fMinR){ // cuts on distance from collision point
- hV0CutMinR->Fill(fPhotonCandidate->GetPhotonMass());
- return kFALSE;
- }
-
- if(fPhotonCandidate->GetConversionRadius() <= ((TMath::Abs(fPhotonCandidate->GetConversionZ())*fLineCutZRSlope)-fLineCutZValue)){
- hV0CutLine->Fill(fPhotonCandidate->GetPhotonMass());
- return kFALSE;
- }
- /*else if (fUseEtaMinCut && fPhotonCandidate->GetConversionRadius() >= ((TMath::Abs(fPhotonCandidate->GetConversionZ())*fLineCutZRSlopeMin)-fLineCutZValueMin )){
- hV0CutLine->Fill(fPhotonCandidate->GetPhotonMass());
- return kFALSE;
- }*/
-
-
- if(TMath::Abs(fPhotonCandidate->GetConversionZ()) > fMaxZ ){ // cuts out regions where we do not reconstruct
- hV0CutZ->Fill(fPhotonCandidate->GetPhotonMass());
- return kFALSE;
- }
-
- if(TMath::Abs(fPhotonCandidate->GetPhotonEta())> fEtaCut || TMath::Abs(fPhotonCandidate->GetPhotonEta())< fEtaCutMin){
- hV0CutEta->Fill(fPhotonCandidate->GetPhotonMass());
- return kFALSE;
- }
-
- if(TMath::Abs(fCurrentNegativeTrack->Eta())> fEtaCut || TMath::Abs(fCurrentNegativeTrack->Eta())< fEtaCutMin){
- hV0CutEta->Fill(fPhotonCandidate->GetPhotonMass());
- return kFALSE;
- }
-
- if(TMath::Abs(fCurrentPositiveTrack->Eta())> fEtaCut || TMath::Abs(fCurrentPositiveTrack->Eta())< fEtaCutMin){
- hV0CutEta->Fill(fPhotonCandidate->GetPhotonMass());
- return kFALSE;
- }
-
- if( fCurrentNegativeTrack->Pt()< fSinglePtCut || fCurrentNegativeTrack->Pt()< fSinglePtCut){
- hV0CutSinglePt->Fill(fPhotonCandidate->GetPhotonMass());
- return kFALSE;
- }
-
-
- if(fPhotonCandidate->GetPhotonPt()<fPtCut){
- hV0CutPt->Fill(fPhotonCandidate->GetPhotonMass());
- return kFALSE;
- }
-return kTRUE;
-}
-
-
-
-///________________________________________________________________________
-Bool_t AliV0ReaderV1::TrackCuts(AliConversionPhotonBase *fPhotonCandidate){
-
- Bool_t passtrackcuts=kTRUE;
-
-
-
- if(fESDEvent){
-
- AliESDtrack *fCurrentNegativeESDTrack=(AliESDtrack*)fESDEvent->GetTrack(fPhotonCandidate->GetTrackLabelNegative());
- AliESDtrack *fCurrentPositiveESDTrack=(AliESDtrack*)fESDEvent->GetTrack(fPhotonCandidate->GetTrackLabelPositive());
-
- if(!fCurrentNegativeESDTrack||!fCurrentPositiveESDTrack)return kFALSE;
-
- // avoid like sign
- if(fCurrentNegativeESDTrack->Charge() == fCurrentPositiveESDTrack->Charge()){
-
- hV0CutLikeSign->Fill(fPhotonCandidate->GetPhotonMass());
- passtrackcuts=kFALSE;
- }
-
- if( (!(fCurrentNegativeESDTrack->IsOn(AliESDtrack::kTPCrefit))||(!(fCurrentPositiveESDTrack->IsOn(AliESDtrack::kTPCrefit))))){
- hV0CutRefit->Fill(fPhotonCandidate->GetPhotonMass());
- passtrackcuts=kFALSE;
-
- }
-
- if( fCurrentNegativeESDTrack->GetKinkIndex(0) > 0 ||
- fCurrentPositiveESDTrack->GetKinkIndex(0) > 0) {
- passtrackcuts=kFALSE;
- }
-
- if(fCurrentNegativeESDTrack->GetNcls(1) < fMinClsTPC || fCurrentPositiveESDTrack->GetNcls(1) < fMinClsTPC ){
- passtrackcuts=kFALSE;
- hV0CutKinks->Fill(fPhotonCandidate->GetPhotonMass());
- }
-
- /* Double_t negclsToF = 0.;
- if (!fUseCorrectedTPCClsInfo ){
- if(fCurrentNegativeESDTrack->GetTPCNclsF()!=0 ){
- negclsToF = (Double_t)fCurrentNegativeESDTrack->GetNcls(1)/(Double_t)fCurrentNegativeESDTrack->GetTPCNclsF();
- }
- } else {
- negclsToF = fCurrentNegativeESDTrack->GetTPCClusterInfo(2,0,GetFirstTPCRow(GetXYRadius()));
- }
-
- Double_t posclsToF = 0.;
- if (!fUseCorrectedTPCClsInfo ){
- if(fCurrentTrack->GetTPCNclsF()!=0 ){
- posclsToF = (Double_t)fCurrentTrack->GetNcls(1)/(Double_t)fCurrentTrack->GetTPCNclsF();
- }
- }else{
- posclsToF = fCurrentTrack->GetTPCClusterInfo(2,0,GetFirstTPCRow(GetXYRadius()));
- }
-
- if( negclsToF < fMinClsTPCToF || posclsToF < fMinClsTPCToF ){
- hV0CutMinNclsTPCToF->Fill(fPhotonCandidate->GetPhotonMass());
- passtrackcuts=kFALSE; }
- */
-
-
-
-
- }
-
- if(fAODEvent){
-
- AliAODTrack *fCurrentNegativeESDTrack=(AliAODTrack*)fAODEvent->GetTrack(fPhotonCandidate->GetTrackLabelNegative());
- AliAODTrack *fCurrentPositiveESDTrack=(AliAODTrack*)fAODEvent->GetTrack(fPhotonCandidate->GetTrackLabelPositive());
-
- if(!fCurrentNegativeESDTrack||!fCurrentPositiveESDTrack)return kFALSE;
-
- // avoid like sign
- if(fCurrentNegativeESDTrack->Charge() == fCurrentPositiveESDTrack->Charge()){
-
- hV0CutLikeSign->Fill(fPhotonCandidate->GetPhotonMass());
- passtrackcuts=kFALSE;
- }
-
- if( !(fCurrentNegativeESDTrack->IsOn(AliESDtrack::kTPCrefit))){
- hV0CutRefit->Fill(fPhotonCandidate->GetPhotonMass());
- passtrackcuts=kFALSE;
- }
-
- if( !(fCurrentPositiveESDTrack->IsOn(AliESDtrack::kTPCrefit))){
- hV0CutRefit->Fill(fPhotonCandidate->GetPhotonMass());
- passtrackcuts=kFALSE;
- }
-
- // to be implemented
- /*
- if( fCurrentNegativeESDTrack->GetKinkIndex(0) > 0 ||
- fCurrentPositiveESDTrack->GetKinkIndex(0) > 0) {
- }*/
-
- if(fCurrentNegativeESDTrack->GetNcls(1) < fMinClsTPC || fCurrentPositiveESDTrack->GetNcls(1) < fMinClsTPC ){
- passtrackcuts=kFALSE;}
-
-
- }
-
- return passtrackcuts;
-}
-
-///________________________________________________________________________
-Bool_t AliV0ReaderV1::dEdxCuts(AliConversionPhotonBase *fPhotonCandidate){
-
- AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
- AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
- AliPIDResponse *fPIDResponse = (AliPIDResponse*)inputHandler->GetPIDResponse();
-
- AliVTrack *fCurrentTrack=0x0;
-
- for(Int_t ilabel=0;ilabel<2;ilabel++){
-
- fCurrentTrack=GetTrack(fPhotonCandidate->GetTrackLabel(ilabel));
-
- if(!fCurrentTrack)return kFALSE;
-
- // Fill dEdx before cuts
-
- hGammadEdxbefore->Fill(fCurrentTrack->P(),fCurrentTrack->GetTPCsignal());
-
- if(fDodEdxSigmaCut == kTRUE){
- if( fPIDResponse->NumberOfSigmasTPC(fCurrentTrack,AliPID::kElectron)<fPIDnSigmaBelowElectronLine ||
- fPIDResponse->NumberOfSigmasTPC(fCurrentTrack,AliPID::kElectron)>fPIDnSigmaAboveElectronLine){
-
- hV0CutdEdxElectron->Fill(fPhotonCandidate->GetPhotonMass());
- return kFALSE;
- }
-
- if( fCurrentTrack->P()>fPIDMinPnSigmaAbovePionLine && fCurrentTrack->P()<fPIDMaxPnSigmaAbovePionLine ){
- if(fPIDResponse->NumberOfSigmasTPC(fCurrentTrack,AliPID::kElectron)>fPIDnSigmaBelowElectronLine &&
- fPIDResponse->NumberOfSigmasTPC(fCurrentTrack,AliPID::kElectron)<fPIDnSigmaAboveElectronLine&&
- fPIDResponse->NumberOfSigmasTPC(fCurrentTrack,AliPID::kPion)<fPIDnSigmaAbovePionLine){
-
- hV0CutdEdxPion->Fill(fPhotonCandidate->GetPhotonMass());
- return kFALSE;
- }
- }
-
- // High Pt
- if( fCurrentTrack->P()>fPIDMaxPnSigmaAbovePionLine ){
- if(fPIDResponse->NumberOfSigmasTPC(fCurrentTrack,AliPID::kElectron)>fPIDnSigmaBelowElectronLine &&
- fPIDResponse->NumberOfSigmasTPC(fCurrentTrack,AliPID::kElectron)<fPIDnSigmaAboveElectronLine&&
- fPIDResponse->NumberOfSigmasTPC(fCurrentTrack,AliPID::kPion)<fPIDnSigmaAbovePionLineHighPt){
-
- hV0CutdEdxPion->Fill(fPhotonCandidate->GetPhotonMass());
- return kFALSE;
- }
- }
- }
-
- if(fDoKaonRejectionLowP == kTRUE){
- if(fCurrentTrack->P()<fPIDMinPKaonRejectionLowP ){
- if( TMath::Abs(fPIDResponse->NumberOfSigmasTPC(fCurrentTrack,AliPID::kKaon))<fPIDnSigmaAtLowPAroundKaonLine){
- hV0CutdEdxKaonLowP->Fill(fPhotonCandidate->GetPhotonMass());
- return kFALSE;
- }
- }
- }
-
- if(fDoProtonRejectionLowP == kTRUE){
- if( fCurrentTrack->P()<fPIDMinPProtonRejectionLowP ){
- if( TMath::Abs(fPIDResponse->NumberOfSigmasTPC(fCurrentTrack,AliPID::kProton))<fPIDnSigmaAtLowPAroundProtonLine){
- hV0CutdEdxProtonLowP->Fill(fPhotonCandidate->GetPhotonMass());
- return kFALSE;
- }
- }
- }
-
- if(fDoPionRejectionLowP == kTRUE){
- if( fCurrentTrack->P()<fPIDMinPPionRejectionLowP ){
- if( TMath::Abs(fPIDResponse->NumberOfSigmasTPC(fCurrentTrack,AliPID::kPion))<fPIDnSigmaAtLowPAroundPionLine){
- hV0CutdEdxPionLowP->Fill(fPhotonCandidate->GetPhotonMass());
- return kFALSE;
- }
- }
- }
-
-
- if( fDoTOFsigmaCut == kTRUE ){ // RRnewTOF start /////////////////////////////////////////////////////////////////////////////
-
- if((fPIDResponse->NumberOfSigmasTOF(fCurrentTrack,AliPID::kElectron)>fTofPIDnSigmaAboveElectronLine) || (fPIDResponse->NumberOfSigmasTOF(fCurrentTrack,AliPID::kElectron)<fTofPIDnSigmaBelowElectronLine)){
- hV0CutdEdxTOFElectron->Fill(fPhotonCandidate->GetPhotonMass());
- return kFALSE;
- }
- } /////////////////////////////// RRnewTOF end ///////////////////////////////////////////////////////////////////////////////
-
- // Apply TRD PID
- if(fDoTRDPID){
- if(!fPIDResponse->IdentifiedAsElectronTRD(fCurrentTrack,fPIDTRDEfficiency)){
- hV0CutdEdxTRD->Fill(fPhotonCandidate->GetPhotonMass());
- return kFALSE;
- }
- }
-
- // Fill dEdx Histogram after Cuts
-
- hGammadEdxafter->Fill(fCurrentTrack->P(),fCurrentTrack->GetTPCsignal());
-
- }
-
- return kTRUE;
-}
-
-///________________________________________________________________________
-Bool_t AliV0ReaderV1::AsymmetryCut(AliConversionPhotonBase *fPhotonCandidate)
-{
-
- for(Int_t ii=0;ii<2;ii++){
- AliVTrack *fCurrentTrack=GetTrack(fPhotonCandidate->GetTrackLabel(ii));
-
- if( fCurrentTrack->P()>fMinPPhotonAsymmetryCut ){
- Double_t trackNegAsy=0;
- if (fPhotonCandidate->GetPhotonP()!=0.){
- trackNegAsy= fCurrentTrack->P()/fPhotonCandidate->GetPhotonP();
- }
- if( trackNegAsy<fMinPhotonAsymmetry ||trackNegAsy>(1.- fMinPhotonAsymmetry)){
- hV0CutAsymmetry->Fill(fPhotonCandidate->GetPhotonMass());
- return kFALSE;
- }
- }
- }
-
- return kTRUE;
-}
-
-
-///________________________________________________________________________
-
-Int_t AliV0ReaderV1::GetNumberOfContributorsVtx(){
- // returns number of contributors to the vertex
- if(fESDEvent){
- if(fESDEvent->GetPrimaryVertexTracks()->GetNContributors()>0) {
- return fESDEvent->GetPrimaryVertexTracks()->GetNContributors();
- }
-
- if(fESDEvent->GetPrimaryVertexTracks()->GetNContributors()<1) {
- // return 0;
- //-AM test pi0s without SPD only vertex
- if(fESDEvent->GetPrimaryVertexSPD()->GetNContributors()>0) {
- return fESDEvent->GetPrimaryVertexSPD()->GetNContributors();
-
- }
- if(fESDEvent->GetPrimaryVertexSPD()->GetNContributors()<1) {
- // cout<<"number of contributors from bad vertex type::"<< fESDEvent->GetPrimaryVertex()->GetName() << endl;
- return 0;
- }
- }
- }
- if(fAODEvent){
- if(fAODEvent->GetPrimaryVertex()->GetNContributors()>0) {
- return fESDEvent->GetPrimaryVertex()->GetNContributors();
- }
- if(fAODEvent->GetPrimaryVertex()->GetNContributors()<1) {
- if(fAODEvent->GetPrimaryVertexSPD()->GetNContributors()>0) {
- return fAODEvent->GetPrimaryVertexSPD()->GetNContributors();
-
- }
- if(fAODEvent->GetPrimaryVertexSPD()->GetNContributors()<1) {
- AliWarning(Form("Number of contributors from bad vertex type:: %s",fAODEvent->GetPrimaryVertex()->GetName()));
- return 0;
- }
- }
- }
-
- return 0;
-
-}
-
-///________________________________________________________________________
-Double_t AliV0ReaderV1::GetCentrality(){
- if(fAODEvent){
- if(fAODEvent->GetHeader()){return fAODEvent->GetHeader()->GetCentrality();}
- }
-
- if(fESDEvent){
- return ((AliCentrality *)fESDEvent->GetCentrality())->GetCentralityPercentile("V0M");
- }
-
- return -1;
-}
-
-///________________________________________________________________________
-
-AliEventplane *AliV0ReaderV1::GetEventPlane(){
- if(fESDEvent){
- return fESDEvent->GetEventplane();
- }
- if(fAODEvent){
- return fAODEvent->GetEventplane();
- }
-return 0x0;
-}
-
-//________________________________________________________________________
-Bool_t AliV0ReaderV1::SetEventPlane()
-{
- if(fMCStack||!fIsHeavyIon){
- // NO EP in MC mode and pp mode
- fEPAngle=0;
- return kTRUE;
- }
-
- if(fIsHeavyIon){
- AliEventplane *fEP=GetEventPlane();
-
- if(!fEP)return kFALSE;
-
- fEPAngle=fEP->GetEventplane("Q");
-
- return kTRUE;
- }
- return kFALSE;
-}
-
-///________________________________________________________________________
-Bool_t AliV0ReaderV1::EventCuts(){
-
- fEventIsSelected=kTRUE;
-
- // Z Vertex Position Cut
-
- if(!VertexZCut()){
- hV0EventCuts->Fill(0);
- fEventIsSelected=kFALSE;;
- }
- // Number of Contributors Cut
-
- if(GetNumberOfContributorsVtx()<=0) {
- hV0EventCuts->Fill(1);
- fEventIsSelected=kFALSE;;
- }
- // Centrality Selection
-
- if(!CentralitySelection()){
- hV0EventCuts->Fill(2);
- fEventIsSelected=kFALSE;;
- }
-
- // Event Plane
- if(!SetEventPlane()){hV0EventCuts->Fill(4);}
-
- // Fill Event Histograms
-
- if(fEventIsSelected){
- // Fill Event Histograms
- hV0EventCuts->Fill(9);
- hVertexZ->Fill(fVertexZ);
- hCentrality->Fill(fCentrality);
- hNEvents->Fill(fCentralityBin);
- }
- else{hV0EventCuts->Fill(8);}
-
- return fEventIsSelected;
-}
-
-///________________________________________________________________________
-Bool_t AliV0ReaderV1::VertexZCut(){
-
- fVertexZ=GetPrimaryVertex()->GetZ();
-
- if(fBGHandler){
- if(fBGHandler->GetZBinIndex(fVertexZ)<0)return kFALSE;
- }
- else{
- if(fVertexZ>fMaxVertexZ)return kFALSE;
- }
- return kTRUE;
-}
-
-///________________________________________________________________________
-AliVTrack *AliV0ReaderV1::GetTrack(Int_t label){
- if(fESDEvent){
- return (AliESDtrack*)fESDEvent->GetTrack(label);
- }
- if(fAODEvent)return (AliAODTrack*)fAODEvent->GetTrack(label);
- return 0x0;
-}
-
-///________________________________________________________________________
-Bool_t AliV0ReaderV1::PIDProbabilityCut(AliConversionPhotonBase *fPhotonCandidate){
-
- if(fESDEvent){
-
- Bool_t iResult=kFALSE;
-
- Double_t *posProbArray = new Double_t[AliPID::kSPECIES];
- Double_t *negProbArray = new Double_t[AliPID::kSPECIES];
-
- AliESDtrack* negTrack = (AliESDtrack*)fESDEvent->GetTrack(fPhotonCandidate->GetTrackLabelNegative());
- AliESDtrack* posTrack = (AliESDtrack*)fESDEvent->GetTrack(fPhotonCandidate->GetTrackLabelPositive());
-
- if(negProbArray && posProbArray){
-
- negTrack->GetTPCpid(negProbArray);
- posTrack->GetTPCpid(posProbArray);
-
- if(negProbArray[AliPID::kElectron]>=fPIDProbabilityCutNegativeParticle && posProbArray[AliPID::kElectron]>=fPIDProbabilityCutPositiveParticle){
- iResult=kTRUE;
- }
- else{hV0CutPIDProb->Fill(fPhotonCandidate->GetPhotonMass());}
- }
-
- delete [] posProbArray;
- delete [] negProbArray;
- return iResult;
-
- }
- if(fAODEvent){
- // not possible to implement
- return kTRUE;}
- return kFALSE;
-}
-
-///________________________________________________________________________
-const AliVertex *AliV0ReaderV1::GetPrimaryVertex() {
- if(fESDEvent)return fESDEvent->GetPrimaryVertex();
- if(fAODEvent)return const_cast<const AliVertex*>(dynamic_cast<AliVertex*>(fAODEvent->GetPrimaryVertex()));
-return 0x0;
-}
-
-///________________________________________________________________________
-Bool_t AliV0ReaderV1::GetHelixCenter(const AliExternalTrackParam *track, Double_t b,Int_t charge, Double_t center[2]){
- // see header file for documentation
-
- Double_t helix[6];
- track->GetHelixParameters(helix,b);
-
- Double_t xpos = helix[5];
- Double_t ypos = helix[0];
- Double_t radius = TMath::Abs(1./helix[4]);
- Double_t phi = helix[2];
-
- if(phi < 0){
- phi = phi + 2*TMath::Pi();
- }
-
- phi -= TMath::Pi()/2.;
- Double_t xpoint = radius * TMath::Cos(phi);
- Double_t ypoint = radius * TMath::Sin(phi);
-
- if(b<0){
- if(charge > 0){
- xpoint = - xpoint;
- ypoint = - ypoint;
- }
-
- if(charge < 0){
- xpoint = xpoint;
- ypoint = ypoint;
- }
- }
- if(b>0){
- if(charge > 0){
- xpoint = xpoint;
- ypoint = ypoint;
- }
-
- if(charge < 0){
- xpoint = - xpoint;
- ypoint = - ypoint;
- }
- }
- center[0] = xpos + xpoint;
- center[1] = ypos + ypoint;
-
- return 1;
-}
-///________________________________________________________________________
-Bool_t AliV0ReaderV1::GetConversionPoint(const AliExternalTrackParam *pparam,const AliExternalTrackParam *nparam,Double_t convpos[3]){
-
- if(!pparam||!nparam)return kFALSE;
-
- Double_t helixcenterpos[2];
- GetHelixCenter(pparam,GetMagneticField(),pparam->Charge(),helixcenterpos);
-
- Double_t helixcenterneg[2];
- GetHelixCenter(nparam,GetMagneticField(),nparam->Charge(),helixcenterneg);
-
- Double_t helixpos[6];
- pparam->GetHelixParameters(helixpos,GetMagneticField());
- Double_t posradius = TMath::Abs(1./helixpos[4]);
-
- Double_t helixneg[6];
- nparam->GetHelixParameters(helixneg,GetMagneticField());
- Double_t negradius = TMath::Abs(1./helixneg[4]);
-
- // Calculate xy-position
-
- Double_t xpos = helixcenterpos[0];
- Double_t ypos = helixcenterpos[1];
- Double_t xneg = helixcenterneg[0];
- Double_t yneg = helixcenterneg[1];
-
- convpos[0] = (xpos*negradius + xneg*posradius)/(negradius+posradius);
- convpos[1] = (ypos*negradius+ yneg*posradius)/(negradius+posradius);
-
-
- // Calculate z-position
-
- Double_t deltaXPos = convpos[0] - xpos;
- Double_t deltaYPos = convpos[1] - ypos;
-
- Double_t deltaXNeg = convpos[0] - xneg;
- Double_t deltaYNeg = convpos[1] - yneg;
-
- Double_t alphaPos = TMath::Pi() + TMath::ATan2(-deltaYPos,-deltaXPos);
- Double_t alphaNeg = TMath::Pi() + TMath::ATan2(-deltaYNeg,-deltaXNeg);
-
- Double_t vertexXNeg = xneg + TMath::Abs(negradius)*
- TMath::Cos(alphaNeg);
- Double_t vertexYNeg = yneg + TMath::Abs(negradius)*
- TMath::Sin(alphaNeg);
-
- Double_t vertexXPos = xpos + TMath::Abs(posradius)*
- TMath::Cos(alphaPos);
- Double_t vertexYPos = ypos + TMath::Abs(posradius)*
- TMath::Sin(alphaPos);
-
- Double_t x0neg = helixneg[5];
- Double_t y0neg = helixneg[0];
-
- Double_t x0pos = helixpos[5];
- Double_t y0pos = helixpos[0];
-
- Double_t dNeg = TMath::Sqrt((vertexXNeg - x0neg)*(vertexXNeg - x0neg)
- +(vertexYNeg - y0neg)*(vertexYNeg - y0neg));
-
- Double_t dPos = TMath::Sqrt((vertexXPos - x0pos)*(vertexXPos - x0pos)
- +(vertexYPos - y0pos)*(vertexYPos - y0pos));
-
- Double_t rNeg = TMath::Sqrt(negradius*negradius -
- dNeg*dNeg/4.);
-
- Double_t rPos = TMath::Sqrt(posradius*posradius -
- dPos*dPos/4.);
-
- Double_t deltabetaNeg = 2*(TMath::Pi() + TMath::ATan2(-dNeg/2.,-rNeg));
- Double_t deltabetaPos = 2*(TMath::Pi() + TMath::ATan2(-dPos/2.,-rPos));
-
- Double_t deltaUNeg = negradius*deltabetaNeg;
- Double_t deltaUPos = posradius*deltabetaPos;
-
- Double_t zphaseNeg = nparam->GetZ() + deltaUNeg * nparam->GetTgl();
- Double_t zphasePos = pparam->GetZ() + deltaUPos * pparam->GetTgl();
-
- convpos[2] = (zphasePos*negradius+zphaseNeg*posradius)/(negradius+posradius);
-
-return kTRUE;
-}
-///________________________________________________________________________
-void AliV0ReaderV1::ProcessMC(AliKFConversionPhoton *fCurrentReconstructedGamma){
-
- if(!fMCStack)return;
-
- TParticle *fMotherMCParticle=NULL;
- TParticle *fNegativeMCParticle=NULL;
- TParticle *fPositiveMCParticle=NULL;
-
- // Get MC Particles
- fMotherMCParticle = fCurrentReconstructedGamma->GetMCParticle(fMCStack);
- fNegativeMCParticle = fCurrentReconstructedGamma->GetNegativeMCDaughter(fMCStack);
- fPositiveMCParticle = fCurrentReconstructedGamma->GetPositiveMCDaughter(fMCStack);
-
- if(fPositiveMCParticle&&fNegativeMCParticle&&fMotherMCParticle){
-
- // Check if it is a true photon
-
- if(fMotherMCParticle->GetPdgCode()==22){
-
- hMCPtRECOTRUE[fCentralityBin]->Fill(fCurrentReconstructedGamma->GetPt());
-
- // Pt Resolution
-
- Double_t mcpt = fMotherMCParticle->Pt();
- Double_t esdpt = fCurrentReconstructedGamma->GetPt();
- Double_t resdPt = 0.;
- if(mcpt > 0){
- resdPt = ((esdpt - mcpt)/mcpt)*100.;
- } else if(mcpt < 0){
- AliWarning("Pt of MC particle is negative, this will cause wrong calculation of resPt");
- }
-
-
- hMCPtResolution->Fill(mcpt,resdPt);
- hMCPtResolutionPhi->Fill(fMotherMCParticle->Phi(),resdPt);
-
- // Conversion Point Resolution
-
- Double_t resdR = 0.;
- if(fNegativeMCParticle->R() != 0){
- resdR = ((fCurrentReconstructedGamma->GetConversionRadius() - fNegativeMCParticle->R())/fNegativeMCParticle->R())*100.;
- }
- Double_t resdRAbs = 0.;
- resdRAbs = (fCurrentReconstructedGamma->GetConversionRadius() - fNegativeMCParticle->R());
-
- hMCRResolutionvsR->Fill(fNegativeMCParticle->R(),resdRAbs);
-
- // fHistograms->FillHistogram("Resolution_dR", fV0Reader->GetNegativeMCParticle()->R(), resdR);
- // fHistograms->FillHistogram("Resolution_MC_R", fV0Reader->GetNegativeMCParticle()->R());
- // fHistograms->FillHistogram("Resolution_ESD_R", fV0Reader->GetXYRadius());
- // fHistograms->FillHistogram("Resolution_R_dPt", fV0Reader->GetNegativeMCParticle()->R(), resdPt);
-
- Double_t resdZ = 0.;
- if(fNegativeMCParticle->Vz() != 0){
- resdZ = ((fCurrentReconstructedGamma->GetZ() -fNegativeMCParticle->Vz())/fNegativeMCParticle->Vz())*100.;
- }
- Double_t resdZAbs = 0.;
- resdZAbs = fCurrentReconstructedGamma->GetZ() -fNegativeMCParticle->Vz();
-
- hMCZResolutionvsZ->Fill( fNegativeMCParticle->Vz(), resdZAbs);
- }
- }
-
-}
-
-///________________________________________________________________________
-void AliV0ReaderV1::ProcessMCGammasForEfficiency(){
-
- if(!fMCStack)return;
-
- for (Int_t iTracks = 0; iTracks < fMCStack->GetNprimary(); iTracks++) {
- TParticle* particle = (TParticle *)fMCStack->Particle(iTracks);
-
- //process the gammas
-
- if(IsMCConversionGammaInAcceptance(particle)){
-
- hMCPtTRUE[fCentralityBin]->Fill(particle->Pt());
-
- }
- }
-}
-
-
-///________________________________________________________________________
-Bool_t AliV0ReaderV1::IsMCConversionGammaInAcceptance(TParticle *particle){
- if(!fMCStack)return kFALSE;
-
- if (particle->GetPdgCode() == 22){
- if(TMath::Abs(particle->Eta())> fEtaCut || TMath::Abs(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
- }
-
- // looking for conversion (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(AcceptanceCut(particle,ePos,eNeg))return kTRUE;
- }
- return kFALSE;
-}
-
-///________________________________________________________________________
-Bool_t AliV0ReaderV1::AcceptanceCut(TParticle *particle, TParticle * ePos,TParticle* eNeg){
-
- // cuts on distance from collision point
-
- if(particle->R()>fMaxR){
- return kFALSE;}
-
- if(ePos->R()>fMaxR){
- return kFALSE;
- }
-
- if(ePos->R()<fMinR){
- return kFALSE;
- }
-
- if( ePos->R() <= ((TMath::Abs(ePos->Vz())*fLineCutZRSlope)-fLineCutZValue)){
- return kFALSE;
- }
- /*else if (fUseEtaMinCut && ePos->R() >= ((TMath::Abs(ePos->Vz())*fLineCutZRSlopeMin)-fLineCutZValueMin )){
- return kFALSE;
- } */
-
- if(TMath::Abs(eNeg->Vz()) > fMaxZ ){ // cuts out regions where we do not reconstruct
- return kFALSE;
- }
-
- if(eNeg->Vz()!=ePos->Vz()||eNeg->R()!=ePos->R()){
- return kFALSE;
- }
-
- if(TMath::Abs(ePos->Vz()) > fMaxZ ){ // cuts out regions where we do not reconstruct
- return kFALSE;
- }
-
- if(TMath::Abs(particle->Eta())> fEtaCut || TMath::Abs(particle->Eta())< fEtaCutMin){
- return kFALSE;
- }
-
- if(TMath::Abs(ePos->Eta())> fEtaCut || TMath::Abs(ePos->Eta())< fEtaCutMin){
- return kFALSE;
- }
-
- if(TMath::Abs(eNeg->Eta())> fEtaCut || TMath::Abs(eNeg->Eta())< fEtaCutMin){
- return kFALSE;
- }
-
- if( ePos->Pt()< fSinglePtCut || eNeg->Pt()< fSinglePtCut){
- return kFALSE;
- }
-
- if(particle->Pt()<fPtCut){
- return kFALSE;
- }
-
- return kTRUE;
-}
-///________________________________________________________________________
-void AliV0ReaderV1::PrintCuts(){
-
- cout<<"V0 Reader initialized with following settings"<<endl;
-
-
- cout<<"Acceptance Eta:"<<endl;
- cout<<fEtaCutMin<<" < eta < "<<fEtaCut<<endl;
- cout<<"Conversion Point"<<endl;
- cout<<"Z <"<<fMaxZ<<endl;
- cout<<fMinR<<" < R < "<<fMaxR<<endl;
- cout<<"Line Cut Slope"<<fLineCutZRSlope<<" ZValue "<<fLineCutZValue<<endl;
-
- cout<<"Pt Gamma > "<<fPtCut<<endl;
- cout<<"Pt Daughters > "<<fSinglePtCut<<endl;
-
- cout<<"Armenteros Qt Cut"<<endl;
-
- if(fDoHighPtQtGammaSelection){
- cout<<" qt < "<<fQtMax<<" for pt < "<<fPtBorderForQt<<endl;
- cout<<" qt < "<<fHighPtQtMax<<" for pt > "<<fPtBorderForQt<<endl;
- }
- else{
- cout<<" qt < "<<fQtMax<<endl;
- }
- cout<<"Chi2perNDF > "<<fChi2CutConversion<<endl;
-
-
-
-}
-
-//_______________________________________________________________________
-
-Bool_t AliV0ReaderV1::CentralitySelection(){
-
- fCentralityBin=-1;
-
- if(!fIsHeavyIon){fCentralityBin=0;fCentrality=0;return kTRUE;}
-
- fCentrality=GetCentrality();
-
- if(fIsHeavyIon){
-
- if(fBGHandler){
- fCentralityBin=fBGHandler->GetCentralityBinIndex(Int_t(fCentrality));
- }
- else{
- Double_t fCentralityBins[fNCentralityBins+1];
- for(int i=0;i<fNCentralityBins;i++){
- fCentralityBins[i]=i*100/Double_t(fNCentralityBins);
- }
-
- for(int i=0;i<fNCentralityBins;i++){
- if(fCentrality>fCentralityBins[i]&&fCentrality<fCentralityBins[i+1]){
- fCentralityBin=i;
- return kTRUE;}
- }
- }
-
- }
- if(fCentralityBin>=0&&fCentrality>=0){
-
- return kTRUE;
- }
-
- AliWarning("Centrality not defined");
- return kFALSE;
-
-}
-
-//________________________________________________________________________
-void AliV0ReaderV1::Terminate(Option_t *)
-{
-
- printf(Form("V0ReaderV1: V0s processed: %4.0f reconstructed photons: %4.0f \n",hV0CurrentFinder->GetEntries(),hV0Good->GetEntries()));
-
-}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff