-/**************************************************************************\r
- * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *\r
- * *\r
- * Author: The ALICE Off-line Project. *\r
- * Contributors are mentioned in the code where appropriate. *\r
- * *\r
- * Permission to use, copy, modify and distribute this software and its *\r
- * documentation strictly for non-commercial purposes is hereby granted *\r
- * without fee, provided that the above copyright notice appears in all *\r
- * copies and that both the copyright notice and this permission notice *\r
- * appear in the supporting documentation. The authors make no claims *\r
- * about the suitability of this software for any purpose. It is *\r
- * provided "as is" without express or implied warranty. *\r
- **************************************************************************/\r
-\r
-#include "AliTwoParticlePIDCorr.h"\r
-#include "AliVParticle.h"\r
-#include "TFormula.h"\r
-#include "TAxis.h"\r
-#include "TChain.h"\r
-#include "TTree.h"\r
-#include "TH1F.h"\r
-#include "TH2F.h"\r
-#include "TH3F.h"\r
-#include "TList.h"\r
-#include "TFile.h"\r
-\r
-#include "AliCentrality.h"\r
-#include "Riostream.h"\r
-\r
-#include <TSpline.h>\r
-#include <AliPID.h>\r
-#include "AliESDpid.h"\r
-#include "AliAODpidUtil.h"\r
-#include <AliPIDResponse.h>\r
-\r
-#include <AliAnalysisManager.h>\r
-#include <AliInputEventHandler.h>\r
-#include "AliAODInputHandler.h"\r
-\r
-#include "AliAnalysisTaskSE.h"\r
-#include "AliAnalysisManager.h"\r
-#include "AliCentrality.h"\r
-\r
-#include "AliVEvent.h"\r
-#include "AliAODEvent.h"\r
-#include "AliAODTrack.h"\r
-#include "AliVTrack.h"\r
-\r
-#include "THnSparse.h"\r
-\r
-#include "AliAODMCHeader.h"\r
-#include "AliAODMCParticle.h"\r
-#include "AliMCEventHandler.h"\r
-#include "AliMCEvent.h"\r
-#include "AliMCParticle.h"\r
-#include "TParticle.h"\r
-#include <TDatabasePDG.h>\r
-#include <TParticlePDG.h>\r
-\r
-#include "AliGenCocktailEventHeader.h"\r
-#include "AliGenEventHeader.h"\r
-\r
-#include "AliEventPoolManager.h"\r
-//#include "AliAnalysisUtils.h"\r
-using namespace AliPIDNameSpace;\r
-using namespace std;\r
-\r
-ClassImp(AliTwoParticlePIDCorr)\r
-ClassImp(LRCParticlePID)\r
-//________________________________________________________________________\r
-AliTwoParticlePIDCorr::AliTwoParticlePIDCorr() // All data members should be initialised here\r
-:AliAnalysisTaskSE(),\r
- fOutput(0),\r
- fCentralityMethod("V0A"),\r
- fSampleType("pPb"),\r
-fnTracksVertex(1), // QA tracks pointing to principal vertex (= 3 default)\r
- trkVtx(0),\r
- zvtx(0),\r
- fFilterBit(768),\r
- fzvtxcut(10.0),\r
- ffilltrigassoUNID(kFALSE),\r
- ffilltrigUNIDassoID(kFALSE),\r
- ffilltrigIDassoUNID(kTRUE),\r
- ffilltrigIDassoID(kFALSE),\r
- ffilltrigIDassoIDMCTRUTH(kFALSE),\r
- fPtOrderMCTruth(kFALSE),\r
- fTriggerSpeciesSelection(kFALSE),\r
- fAssociatedSpeciesSelection(kFALSE),\r
- fTriggerSpecies(SpPion),\r
- fAssociatedSpecies(SpPion),\r
- fCustomBinning(""),\r
- fBinningString(""),\r
- fcontainPIDtrig(kTRUE),\r
- fcontainPIDasso(kFALSE),\r
-//frejectPileUp(kFALSE),\r
- fminPt(0.2),\r
- fmaxPt(10.0),\r
- fmineta(-0.8),\r
- fmaxeta(0.8),\r
- fminprotonsigmacut(-6.0),\r
- fmaxprotonsigmacut(-3.0),\r
- fminpionsigmacut(0.0),\r
- fmaxpionsigmacut(4.0),\r
- fselectprimaryTruth(kTRUE),\r
- fonlyprimarydatareco(kFALSE),\r
- fdcacut(kFALSE),\r
- fdcacutvalue(3.0),\r
- ffillhistQAReco(kFALSE),\r
- ffillhistQATruth(kFALSE),\r
- kTrackVariablesPair(0),\r
- fminPtTrig(0),\r
- fmaxPtTrig(0),\r
- fminPtComboeff(2.0),\r
- fmaxPtComboeff(4.0), \r
- fminPtAsso(0),\r
- fmaxPtAsso(0), \r
- fhistcentrality(0),\r
- fEventCounter(0),\r
- fEtaSpectrasso(0),\r
- fphiSpectraasso(0),\r
- MCtruthpt(0),\r
- MCtrutheta(0),\r
- MCtruthphi(0),\r
- MCtruthpionpt(0),\r
- MCtruthpioneta(0),\r
- MCtruthpionphi(0),\r
- MCtruthkaonpt(0),\r
- MCtruthkaoneta(0),\r
- MCtruthkaonphi(0),\r
- MCtruthprotonpt(0),\r
- MCtruthprotoneta(0),\r
- MCtruthprotonphi(0),\r
- fPioncont(0),\r
- fKaoncont(0),\r
- fProtoncont(0),\r
- fCentralityCorrelation(0x0),\r
- fHistoTPCdEdx(0x0),\r
- fHistoTOFbeta(0x0),\r
- fTPCTOFPion3d(0),\r
- fTPCTOFKaon3d(0),\r
- fTPCTOFProton3d(0),\r
- fPionPt(0),\r
- fPionEta(0),\r
- fPionPhi(0),\r
- fKaonPt(0),\r
- fKaonEta(0),\r
- fKaonPhi(0),\r
- fProtonPt(0),\r
- fProtonEta(0),\r
- fProtonPhi(0),\r
- fCorrelatonTruthPrimary(0),\r
- fCorrelatonTruthPrimarymix(0),\r
- fTHnCorrUNID(0),\r
- fTHnCorrUNIDmix(0),\r
- fTHnCorrID(0),\r
- fTHnCorrIDmix(0),\r
- fTHnCorrIDUNID(0),\r
- fTHnCorrIDUNIDmix(0),\r
- fTHnTrigcount(0),\r
- fTHnTrigcountMCTruthPrim(0),\r
- fPoolMgr(0x0),\r
- fArrayMC(0),\r
- fAnalysisType("MCAOD"), \r
- fefffilename(""),\r
- twoTrackEfficiencyCutValue(0.02),\r
-//fControlConvResoncances(0),\r
- fPID(NULL),\r
- eventno(0),\r
- fPtTOFPIDmin(0.6),\r
- fPtTOFPIDmax(4.0),\r
- fRequestTOFPID(kTRUE),\r
- fPIDType(NSigmaTPCTOF),\r
- fNSigmaPID(3),\r
- fUseExclusiveNSigma(kFALSE),\r
- fRemoveTracksT0Fill(kFALSE),\r
-fSelectCharge(0),\r
-fTriggerSelectCharge(0),\r
-fAssociatedSelectCharge(0),\r
-fTriggerRestrictEta(-1),\r
-fEtaOrdering(kFALSE),\r
-fCutConversions(kFALSE),\r
-fCutResonances(kFALSE),\r
-fRejectResonanceDaughters(-1),\r
- fOnlyOneEtaSide(0),\r
-fInjectedSignals(kFALSE),\r
- fRemoveWeakDecays(kFALSE),\r
-fRemoveDuplicates(kFALSE),\r
- fapplyTrigefficiency(kFALSE),\r
- fapplyAssoefficiency(kFALSE),\r
- ffillefficiency(kFALSE),\r
- fmesoneffrequired(kFALSE),\r
- fkaonprotoneffrequired(kFALSE),\r
-//fAnalysisUtils(0x0),\r
- fDCAXYCut(0) \r
-\r
-{ \r
- for ( Int_t i = 0; i < 16; i++) { \r
- fHistQA[i] = NULL;\r
- }\r
-\r
- for ( Int_t i = 0; i < 6; i++ ){\r
- fTHnrecomatchedallPid[i] = NULL;\r
- fTHngenprimPidTruth[i] = NULL;\r
- effcorection[i]=NULL;\r
- //effmap[i]=NULL;\r
-\r
- }\r
-\r
- for(Int_t ipart=0;ipart<NSpecies;ipart++)\r
- for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++)\r
- fnsigmas[ipart][ipid]=999.;\r
-\r
- for(Int_t ipart=0;ipart<NSpecies;ipart++) {fHasDoubleCounting[ipart]=kFALSE;}\r
-\r
- }\r
-//________________________________________________________________________\r
-AliTwoParticlePIDCorr::AliTwoParticlePIDCorr(const char *name) // All data members should be initialised here\r
- :AliAnalysisTaskSE(name),\r
- fOutput(0),\r
- fCentralityMethod("V0A"),\r
- fSampleType("pPb"),\r
- fnTracksVertex(1), // QA tracks pointing to principal vertex (= 3 default)\r
- trkVtx(0),\r
- zvtx(0),\r
- fFilterBit(768),\r
- fzvtxcut(10.0),\r
- ffilltrigassoUNID(kFALSE),\r
- ffilltrigUNIDassoID(kFALSE),\r
- ffilltrigIDassoUNID(kTRUE),\r
- ffilltrigIDassoID(kFALSE),\r
- ffilltrigIDassoIDMCTRUTH(kFALSE),\r
- fPtOrderMCTruth(kFALSE),\r
- fTriggerSpeciesSelection(kFALSE),\r
- fAssociatedSpeciesSelection(kFALSE),\r
- fTriggerSpecies(SpPion),\r
- fAssociatedSpecies(SpPion),\r
- fCustomBinning(""),\r
- fBinningString(""),\r
- fcontainPIDtrig(kTRUE),\r
- fcontainPIDasso(kFALSE),\r
- // frejectPileUp(kFALSE),\r
- fminPt(0.2),\r
- fmaxPt(10.0),\r
- fmineta(-0.8),\r
- fmaxeta(0.8),\r
- fminprotonsigmacut(-6.0),\r
- fmaxprotonsigmacut(-3.0),\r
- fminpionsigmacut(0.0),\r
- fmaxpionsigmacut(4.0),\r
- fselectprimaryTruth(kTRUE),\r
- fonlyprimarydatareco(kFALSE),\r
- fdcacut(kFALSE),\r
- fdcacutvalue(3.0),\r
- ffillhistQAReco(kFALSE),\r
- ffillhistQATruth(kFALSE),\r
- kTrackVariablesPair(0),\r
- fminPtTrig(0),\r
- fmaxPtTrig(0),\r
- fminPtComboeff(2.0),\r
- fmaxPtComboeff(4.0), \r
- fminPtAsso(0),\r
- fmaxPtAsso(0), \r
- fhistcentrality(0),\r
- fEventCounter(0),\r
- fEtaSpectrasso(0),\r
- fphiSpectraasso(0),\r
- MCtruthpt(0),\r
- MCtrutheta(0),\r
- MCtruthphi(0),\r
- MCtruthpionpt(0),\r
- MCtruthpioneta(0),\r
- MCtruthpionphi(0),\r
- MCtruthkaonpt(0),\r
- MCtruthkaoneta(0),\r
- MCtruthkaonphi(0),\r
- MCtruthprotonpt(0),\r
- MCtruthprotoneta(0),\r
- MCtruthprotonphi(0),\r
- fPioncont(0),\r
- fKaoncont(0),\r
- fProtoncont(0),\r
- fCentralityCorrelation(0x0),\r
- fHistoTPCdEdx(0x0),\r
- fHistoTOFbeta(0x0),\r
- fTPCTOFPion3d(0),\r
- fTPCTOFKaon3d(0),\r
- fTPCTOFProton3d(0),\r
- fPionPt(0),\r
- fPionEta(0),\r
- fPionPhi(0),\r
- fKaonPt(0),\r
- fKaonEta(0),\r
- fKaonPhi(0),\r
- fProtonPt(0),\r
- fProtonEta(0),\r
- fProtonPhi(0),\r
- fCorrelatonTruthPrimary(0),\r
- fCorrelatonTruthPrimarymix(0),\r
- fTHnCorrUNID(0),\r
- fTHnCorrUNIDmix(0),\r
- fTHnCorrID(0),\r
- fTHnCorrIDmix(0),\r
- fTHnCorrIDUNID(0),\r
- fTHnCorrIDUNIDmix(0),\r
- fTHnTrigcount(0),\r
- fTHnTrigcountMCTruthPrim(0),\r
- fPoolMgr(0x0),\r
- fArrayMC(0),\r
- fAnalysisType("MCAOD"),\r
- fefffilename(""), \r
- twoTrackEfficiencyCutValue(0.02),\r
-//fControlConvResoncances(0),\r
- fPID(NULL),\r
- eventno(0),\r
- fPtTOFPIDmin(0.6),\r
- fPtTOFPIDmax(4.0),\r
- fRequestTOFPID(kTRUE),\r
- fPIDType(NSigmaTPCTOF),\r
- fNSigmaPID(3),\r
- fUseExclusiveNSigma(kFALSE),\r
- fRemoveTracksT0Fill(kFALSE),\r
-fSelectCharge(0),\r
-fTriggerSelectCharge(0),\r
-fAssociatedSelectCharge(0),\r
-fTriggerRestrictEta(-1),\r
-fEtaOrdering(kFALSE),\r
-fCutConversions(kFALSE),\r
-fCutResonances(kFALSE),\r
-fRejectResonanceDaughters(-1),\r
- fOnlyOneEtaSide(0),\r
-fInjectedSignals(kFALSE),\r
- fRemoveWeakDecays(kFALSE),\r
-fRemoveDuplicates(kFALSE),\r
- fapplyTrigefficiency(kFALSE),\r
- fapplyAssoefficiency(kFALSE),\r
- ffillefficiency(kFALSE),\r
- fmesoneffrequired(kFALSE),\r
- fkaonprotoneffrequired(kFALSE),\r
- //fAnalysisUtils(0x0),\r
- fDCAXYCut(0) \r
-{\r
- \r
- for ( Int_t i = 0; i < 16; i++) { \r
- fHistQA[i] = NULL;\r
- }\r
- \r
-for ( Int_t i = 0; i < 6; i++ ){\r
- fTHnrecomatchedallPid[i] = NULL;\r
- fTHngenprimPidTruth[i] = NULL;\r
- effcorection[i]=NULL;\r
- //effmap[i]=NULL;\r
-\r
- }\r
-\r
- for(Int_t ipart=0;ipart<NSpecies;ipart++)\r
- for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++)\r
- fnsigmas[ipart][ipid]=999.;\r
-\r
- for(Int_t ipart=0;ipart<NSpecies;ipart++) {fHasDoubleCounting[ipart]=kFALSE;}\r
-\r
- // The last in the above list should not have a comma after it\r
- \r
- // Constructor\r
- // Define input and output slots here (never in the dummy constructor)\r
- // Input slot #0 works with a TChain - it is connected to the default input container\r
- // Output slot #1 writes into a TH1 container\r
- \r
- DefineOutput(1, TList::Class()); // for output list\r
-\r
-}\r
-\r
-//________________________________________________________________________\r
-AliTwoParticlePIDCorr::~AliTwoParticlePIDCorr()\r
-{\r
- // Destructor. Clean-up the output list, but not the histograms that are put inside\r
- // (the list is owner and will clean-up these histograms). Protect in PROOF case.\r
- if (fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {\r
- delete fOutput;\r
-\r
- }\r
- if (fPID) delete fPID;\r
- \r
- }\r
-//________________________________________________________________________\r
-Float_t AliTwoParticlePIDCorr::PhiRange(Float_t DPhi)\r
-\r
-{\r
- //\r
- // Puts the argument in the range [-pi/2,3 pi/2].\r
- //\r
- \r
- if (DPhi < -TMath::Pi()/2) DPhi += 2*TMath::Pi();\r
- if (DPhi > 3*TMath::Pi()/2) DPhi -= 2*TMath::Pi(); \r
-\r
- return DPhi;\r
- \r
-}\r
-//________________________________________________________________________\r
-void AliTwoParticlePIDCorr::UserCreateOutputObjects()\r
-{\r
- // Create histograms\r
- // Called once (on the worker node)\r
- AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();\r
- AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());\r
- fPID = inputHandler->GetPIDResponse();\r
-\r
- //AliAnalysisUtils *fUtils = new AliAnalysisUtils();\r
-\r
-//get the efficiency correction map\r
-\r
-\r
- fOutput = new TList();\r
- fOutput->SetOwner(); // IMPORTANT! \r
-\r
- Int_t centmultbins=10;\r
- Double_t centmultmin=0.0;\r
- Double_t centmultmax=100.0;\r
- if(fSampleType=="pPb" || fSampleType=="PbPb")\r
- {\r
- centmultbins=10;\r
- centmultmin=0.0;\r
- centmultmax=100.0;\r
- }\r
- if(fSampleType=="pp")\r
- {\r
- centmultbins=10;\r
- centmultmin=0.0;\r
- centmultmax=200.0;\r
- }\r
-\r
-fhistcentrality=new TH1F("fhistcentrality","fhistcentrality",centmultbins*4,centmultmin,centmultmax);\r
-fOutput->Add(fhistcentrality);\r
-\r
- fEventCounter = new TH1F("fEventCounter","EventCounter", 10, 0.5,10.5);\r
- fEventCounter->GetXaxis()->SetBinLabel(1,"Event Accesed");\r
- fEventCounter->GetXaxis()->SetBinLabel(2,"Have a vertex");\r
- fEventCounter->GetXaxis()->SetBinLabel(5,"After vertex Cut");\r
- fEventCounter->GetXaxis()->SetBinLabel(6,"Within 0-100% centrality");\r
- fEventCounter->GetXaxis()->SetBinLabel(7,"Event Analyzed");\r
- //fEventCounter->GetXaxis()->SetBinLabel(8,"Event Analysis finished");\r
- fOutput->Add(fEventCounter);\r
- \r
-fEtaSpectrasso=new TH2F("fEtaSpectraasso","fEtaSpectraasso",180,-0.9,0.9,100,0.,20. );\r
-fOutput->Add(fEtaSpectrasso);\r
-\r
-fphiSpectraasso=new TH2F("fphiSpectraasso","fphiSpectraasso",72,0,2*TMath::Pi(),100,0.,20.);\r
-fOutput->Add(fphiSpectraasso);\r
-\r
-\r
- if(fSampleType=="pPb" || fSampleType=="PbPb"){ fCentralityCorrelation = new TH2D("fCentralityCorrelation", ";centrality;multiplicity", 101, 0, 101, 20000, 0,40000);\r
- fOutput->Add(fCentralityCorrelation);\r
- }\r
-\r
-fHistoTPCdEdx = new TH2F("hHistoTPCdEdx", ";p_{T} (GeV/c);dE/dx (au.)",200,0.0,10.0,500, 0., 500.);\r
-fOutput->Add(fHistoTPCdEdx);\r
-fHistoTOFbeta = new TH2F(Form("hHistoTOFbeta"), ";p_{T} (GeV/c);v/c",70, 0., 7., 500, 0.1, 1.1);\r
- fOutput->Add(fHistoTOFbeta);\r
- \r
- fTPCTOFPion3d=new TH3F ("fTPCTOFpion3d", "fTPCTOFpion3d",100,0., 10., 120,-60.,60.,120,-60.,60);\r
- fOutput->Add(fTPCTOFPion3d);\r
- \r
- fTPCTOFKaon3d=new TH3F ("fTPCTOFKaon3d", "fTPCTOFKaon3d",100,0., 10., 120,-60.,60.,120,-60.,60);\r
- fOutput->Add(fTPCTOFKaon3d);\r
-\r
- fTPCTOFProton3d=new TH3F ("fTPCTOFProton3d", "fTPCTOFProton3d",100,0., 10., 120,-60.,60.,120,-60.,60);\r
- fOutput->Add(fTPCTOFProton3d);\r
-\r
-if(ffillhistQAReco)\r
- {\r
- fPionPt = new TH1F("fHistQAPionPt","p_{T} distribution",200,0.,10.);\r
- fOutput->Add(fPionPt);\r
- fPionEta= new TH1F("fHistQAPionEta","#eta distribution",360,-1.8,1.8);\r
- fOutput->Add(fPionEta);\r
- fPionPhi = new TH1F("fHistQAPionPhi","#phi distribution",340,0,6.8);\r
- fOutput->Add(fPionPhi);\r
- \r
- fKaonPt = new TH1F("fHistQAKaonPt","p_{T} distribution",200,0.,10.);\r
- fOutput->Add(fKaonPt);\r
- fKaonEta= new TH1F("fHistQAKaonEta","#eta distribution",360,-1.8,1.8);\r
- fOutput->Add(fKaonEta);\r
- fKaonPhi = new TH1F("fHistQAKaonPhi","#phi distribution",340,0,6.8);\r
- fOutput->Add(fKaonPhi);\r
- \r
- fProtonPt = new TH1F("fHistQAProtonPt","p_{T} distribution",200,0.,10.);\r
- fOutput->Add(fProtonPt);\r
- fProtonEta= new TH1F("fHistQAProtonEta","#eta distribution",360,-1.8,1.8);\r
- fOutput->Add(fProtonEta);\r
- fProtonPhi= new TH1F("fHistQAProtonPhi","#phi distribution",340,0,6.8);\r
- fOutput->Add(fProtonPhi);\r
- }\r
-\r
- fHistQA[0] = new TH1F("fHistQAvx", "Histo Vx All ", 50, -5., 5.);\r
- fHistQA[1] = new TH1F("fHistQAvy", "Histo Vy All", 50, -5., 5.);\r
- fHistQA[2] = new TH1F("fHistQAvz", "Histo Vz All", 50, -25., 25.); \r
- fHistQA[3] = new TH1F("fHistQAvxA", "Histo Vx After Cut ", 50, -5., 5.);\r
- fHistQA[4] = new TH1F("fHistQAvyA", "Histo Vy After Cut", 50, -5., 5.);\r
- fHistQA[5] = new TH1F("fHistQAvzA", "Histo Vz After Cut", 50, -25., 25.);\r
- fHistQA[6] = new TH1F("fHistQADcaXyC", "Histo DCAxy after cut", 50, -5., 5.);\r
- fHistQA[7] = new TH1F("fHistQADcaZC", "Histo DCAz after cut", 50, -5., 5.); \r
- fHistQA[8] = new TH1F("fHistQAPt","p_{T} distribution",200,0.,10.);\r
- fHistQA[9] = new TH1F("fHistQAEta","#eta distribution",360,-1.8,1.8);\r
- fHistQA[10] = new TH1F("fHistQAPhi","#phi distribution",340,0,6.8);\r
- fHistQA[11] = new TH1F("fHistQANCls","Number of TPC cluster",200,0,200);\r
- fHistQA[13] = new TH1F("fHistQAChi2","Chi2 per NDF",100,0,10);\r
- fHistQA[12] = new TH1F("fHistQANCls1","Number of TPC cluster1",200,0,200);\r
- fHistQA[14] = new TH1F("nCrossedRowsTPC","Number of TPC ccrossed rows",200,0,200);\r
- fHistQA[15] = new TH1F("ratioCrossedRowsOverFindableClustersTPC","Number of TPC ccrossed rows find clusters",200,0,2);\r
-\r
-for(Int_t i = 0; i < 16; i++)\r
- {\r
- fOutput->Add(fHistQA[i]);\r
- }\r
-\r
- kTrackVariablesPair=6 ;\r
-\r
- if(fcontainPIDtrig && !fcontainPIDasso) kTrackVariablesPair=7;\r
- \r
- if(!fcontainPIDtrig && fcontainPIDasso) kTrackVariablesPair=7;\r
- \r
- if(fcontainPIDtrig && fcontainPIDasso) kTrackVariablesPair=8;\r
- \r
- \r
-// two particle histograms\r
- Int_t iBinPair[kTrackVariablesPair]; // binning for track variables\r
- Double_t* dBinsPair[kTrackVariablesPair]; // bins for track variables \r
- TString* axisTitlePair; // axis titles for track variables\r
- axisTitlePair=new TString[kTrackVariablesPair];\r
-\r
- TString defaultBinningStr;\r
- defaultBinningStr = "eta: -1.0, -0.9, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0\n"\r
- "p_t_assoc: 0.5, 0.75, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 8.0,10.0\n"\r
- "p_t_leading_course: 0.5, 1.0, 2.0, 3.0, 4.0, 6.0, 8.0,10.0\n"\r
- "p_t_eff:0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.25, 1.5, 1.75, 2.0, 2.25, 2.5, 2.75, 3.0, 3.25, 3.5, 3.75, 4.0, 4.5, 5.0,5.5, 6.0, 7.0, 8.0,9.0,10.0\n"\r
- "vertex: -10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10\n"\r
- "delta_phi: -1.570796, -1.483530, -1.396263, -1.308997, -1.221730, -1.134464, -1.047198, -0.959931, -0.872665, -0.785398, -0.698132, -0.610865, -0.523599, -0.436332, -0.349066, -0.261799, -0.174533, -0.087266, 0.0, 0.087266, 0.174533, 0.261799, 0.349066, 0.436332, 0.523599, 0.610865, 0.698132, 0.785398, 0.872665, 0.959931, 1.047198, 1.134464, 1.221730, 1.308997, 1.396263, 1.483530, 1.570796, 1.658063, 1.745329, 1.832596, 1.919862, 2.007129, 2.094395, 2.181662, 2.268928, 2.356194, 2.443461, 2.530727, 2.617994, 2.705260, 2.792527, 2.879793, 2.967060, 3.054326, 3.141593, 3.228859, 3.316126, 3.403392, 3.490659, 3.577925, 3.665191, 3.752458, 3.839724, 3.926991, 4.014257, 4.101524, 4.188790, 4.276057, 4.363323, 4.450590, 4.537856, 4.625123, 4.712389\n" // this binning starts at -pi/2 and is modulo 3 \r
- "delta_eta: -2.4, -2.3, -2.2, -2.1, -2.0, -1.9, -1.8, -1.7, -1.6, -1.5, -1.4, -1.3, -1.2, -1.1, -1.0, -0.9, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0,2.1, 2.2, 2.3, 2.4\n"\r
- "multiplicity: 0, 1, 2, 3, 4, 5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100.1\n";\r
-\r
- if(fcontainPIDtrig){\r
- defaultBinningStr += "PIDTrig: -0.5,0.5,1.5,2.5,3.5\n"; // course\r
- }\r
- if(fcontainPIDasso){\r
- defaultBinningStr += "PIDAsso: -0.5,0.5,1.5,2.5,3.5\n"; // course\r
- }\r
- // =========================================================\r
- // Customization (adopted from AliUEHistograms)\r
- // =========================================================\r
-\r
- TObjArray* lines = defaultBinningStr.Tokenize("\n");\r
- for (Int_t i=0; i<lines->GetEntriesFast(); i++)\r
- {\r
- TString line(lines->At(i)->GetName());\r
- TString tag = line(0, line.Index(":")+1);\r
- if (!fCustomBinning.BeginsWith(tag) && !fCustomBinning.Contains(TString("\n") + tag))\r
- fBinningString += line + "\n";\r
- else\r
- AliInfo(Form("Using custom binning for %s", tag.Data()));\r
- }\r
- delete lines;\r
- fBinningString += fCustomBinning;\r
- \r
- AliInfo(Form("Used THn Binning:\n%s",fBinningString.Data()));\r
-\r
- // =========================================================\r
- // Now set the bins\r
- // =========================================================\r
-\r
- dBinsPair[0] = GetBinning(fBinningString, "multiplicity", iBinPair[0]);\r
- axisTitlePair[0] = " multiplicity";\r
-\r
- dBinsPair[1] = GetBinning(fBinningString, "vertex", iBinPair[1]);\r
- axisTitlePair[1] = "v_{Z} (cm)"; \r
-\r
- dBinsPair[2] = GetBinning(fBinningString, "p_t_leading_course", iBinPair[2]);\r
- axisTitlePair[2] = "p_{T,trig.} (GeV/c)"; \r
-\r
- dBinsPair[3] = GetBinning(fBinningString, "p_t_assoc", iBinPair[3]);\r
- axisTitlePair[3] = "p_{T,assoc.} (GeV/c)";\r
-\r
- dBinsPair[4] = GetBinning(fBinningString, "delta_eta", iBinPair[4]);\r
- axisTitlePair[4] = "#Delta#eta"; \r
-\r
- dBinsPair[5] = GetBinning(fBinningString, "delta_phi", iBinPair[5]);\r
- axisTitlePair[5] = "#Delta#varphi (rad)"; \r
-\r
- if(fcontainPIDtrig && !fcontainPIDasso){\r
- dBinsPair[6] = GetBinning(fBinningString, "PIDTrig", iBinPair[6]);\r
- axisTitlePair[6] = "PIDTrig"; \r
- }\r
-\r
- if(!fcontainPIDtrig && fcontainPIDasso){\r
- dBinsPair[6] = GetBinning(fBinningString, "PIDAsso", iBinPair[6]);\r
- axisTitlePair[6] = "PIDAsso"; \r
- }\r
-\r
-if(fcontainPIDtrig && fcontainPIDasso){\r
-\r
- dBinsPair[6] = GetBinning(fBinningString, "PIDTrig", iBinPair[6]);\r
- axisTitlePair[6] = "PIDTrig";\r
-\r
- dBinsPair[7] = GetBinning(fBinningString, "PIDAsso", iBinPair[7]);\r
- axisTitlePair[7] = "PIDAsso"; \r
- }\r
- \r
- Int_t nEtaBin = -1;\r
- Double_t* EtaBin = GetBinning(fBinningString, "eta", nEtaBin);\r
- \r
- Int_t nPteffbin = -1;\r
- Double_t* Pteff = GetBinning(fBinningString, "p_t_eff", nPteffbin);\r
-\r
-\r
- fminPtTrig=dBinsPair[2][0];\r
- fmaxPtTrig=dBinsPair[2][iBinPair[2]];\r
- fminPtAsso=dBinsPair[3][0];\r
- fmaxPtAsso=dBinsPair[3][iBinPair[3]];\r
-\r
- //fminPtComboeff=fminPtTrig;***then this value will be fixed ,even Setter can't change it's value\r
- //fmaxPtComboeff=fmaxPtTrig;\r
-//THnSparses for calculation of efficiency\r
-\r
- if((fAnalysisType =="MCAOD") && ffillefficiency) {\r
- const Int_t nDim = 4;// cent zvtx pt eta\r
- Int_t fBinsCh[nDim] = {iBinPair[0], iBinPair[1], nPteffbin ,nEtaBin};//**********************change it\r
- Double_t fMinCh[nDim] = { dBinsPair[0][0],dBinsPair[1][0], Pteff[0], EtaBin[0] };\r
- Double_t fMaxCh[nDim] = { dBinsPair[0][iBinPair[0]], dBinsPair[1][iBinPair[1]], Pteff[nPteffbin], EtaBin[nEtaBin]};\r
-\r
-TString Histrename;\r
-for(Int_t jj=0;jj<6;jj++)//PID type binning\r
- {\r
- Histrename="fTHnrecomatchedallPid";Histrename+=jj;\r
- fTHnrecomatchedallPid[jj] = new THnSparseF(Histrename.Data(),"cent:zvtx::Pt:eta", nDim, fBinsCh, fMinCh, fMaxCh); \r
- fTHnrecomatchedallPid[jj]->Sumw2(); \r
- fTHnrecomatchedallPid[jj]->GetAxis(0)->Set(iBinPair[0], dBinsPair[0]);\r
- fTHnrecomatchedallPid[jj]->GetAxis(0)->SetTitle("Centrality");\r
- fTHnrecomatchedallPid[jj]->GetAxis(1)->Set(iBinPair[1],dBinsPair[1]);\r
- fTHnrecomatchedallPid[jj]->GetAxis(1)->SetTitle("zvtx"); \r
- fTHnrecomatchedallPid[jj]->GetAxis(2)->Set(nPteffbin, Pteff); \r
- fTHnrecomatchedallPid[jj]->GetAxis(2)->SetTitle("p_{T} (GeV/c)");\r
- fTHnrecomatchedallPid[jj]->GetAxis(3)->Set(nEtaBin,EtaBin);\r
- fTHnrecomatchedallPid[jj]->GetAxis(3)->SetTitle("#eta");\r
- fOutput->Add(fTHnrecomatchedallPid[jj]);\r
-\r
-Histrename="fTHngenprimPidTruth";Histrename+=jj;\r
- fTHngenprimPidTruth[jj] = new THnSparseF(Histrename.Data(),"cent:zvtx::Pt:eta", nDim, fBinsCh, fMinCh, fMaxCh);\r
- fTHngenprimPidTruth[jj]->Sumw2(); \r
- fTHngenprimPidTruth[jj]->GetAxis(0)->Set(iBinPair[0],dBinsPair[0]);\r
- fTHngenprimPidTruth[jj]->GetAxis(0)->SetTitle("Centrality");\r
- fTHngenprimPidTruth[jj]->GetAxis(1)->Set(iBinPair[1], dBinsPair[1]);\r
- fTHngenprimPidTruth[jj]->GetAxis(1)->SetTitle("zvtx"); \r
- fTHngenprimPidTruth[jj]->GetAxis(2)->Set(nPteffbin, Pteff); \r
- fTHngenprimPidTruth[jj]->GetAxis(2)->SetTitle("p_{T} (GeV/c)");\r
- fTHngenprimPidTruth[jj]->GetAxis(3)->Set(nEtaBin, EtaBin);\r
- fTHngenprimPidTruth[jj]->GetAxis(3)->SetTitle("#eta");\r
- fOutput->Add(fTHngenprimPidTruth[jj]);\r
- }\r
- }\r
-\r
- Int_t fBins[kTrackVariablesPair];\r
- Double_t fMin[kTrackVariablesPair];\r
- Double_t fMax[kTrackVariablesPair];\r
-\r
-//ThnSparses for Correlation plots(data & MC)\r
- for(Int_t i=0;i<kTrackVariablesPair;i++)\r
- {\r
- fBins[i] =iBinPair[i];\r
- fMin[i]= dBinsPair[i][0];\r
- fMax[i] = dBinsPair[i][iBinPair[i]];\r
- }\r
- if(ffilltrigassoUNID)\r
- {\r
- fTHnCorrUNID = new THnSparseF("fTHnCorrUNID","cent:zvtx:pttrig:ptasso:deltaeta:deltaphi", kTrackVariablesPair, fBins, fMin, fMax); \r
- fTHnCorrUNID->Sumw2();\r
- for(Int_t i=0; i<kTrackVariablesPair;i++){\r
- SetAsymmetricBin(fTHnCorrUNID,i,dBinsPair[i],iBinPair[i],axisTitlePair[i]);\r
- }\r
- fOutput->Add(fTHnCorrUNID);\r
- \r
-fTHnCorrUNIDmix = new THnSparseF("fTHnCorrUNIDmix","cent:zvtx:pttrig:ptasso:deltaeta:deltaphi", kTrackVariablesPair, fBins, fMin, fMax); \r
- fTHnCorrUNIDmix->Sumw2();\r
-for(Int_t i=0; i<kTrackVariablesPair;i++){\r
- SetAsymmetricBin(fTHnCorrUNIDmix,i,dBinsPair[i],iBinPair[i],axisTitlePair[i]);\r
- }\r
- fOutput->Add(fTHnCorrUNIDmix);\r
- }\r
-\r
- if(ffilltrigIDassoID)\r
- {\r
- fTHnCorrID = new THnSparseF("fTHnCorrID","cent:zvtx:pttrig:ptasso:deltaeta:deltaphi", kTrackVariablesPair, fBins, fMin, fMax); \r
- fTHnCorrID->Sumw2();\r
-for(Int_t i=0; i<kTrackVariablesPair;i++){\r
- SetAsymmetricBin(fTHnCorrID,i,dBinsPair[i],iBinPair[i],axisTitlePair[i]);\r
- }\r
- fOutput->Add(fTHnCorrID);\r
-\r
-fTHnCorrIDmix = new THnSparseF("fTHnCorrIDmix","cent:zvtx:pttrig:ptasso:deltaeta:deltaphi", kTrackVariablesPair, fBins, fMin, fMax); \r
- fTHnCorrIDmix->Sumw2();\r
-for(Int_t i=0; i<kTrackVariablesPair;i++){\r
- SetAsymmetricBin(fTHnCorrIDmix,i,dBinsPair[i],iBinPair[i],axisTitlePair[i]);\r
- }\r
- fOutput->Add(fTHnCorrIDmix);\r
- }\r
-\r
- if(ffilltrigUNIDassoID || ffilltrigIDassoUNID)//***********a bit tricky, be careful\r
- {\r
- fTHnCorrIDUNID = new THnSparseF("fTHnCorrIDUNID","cent:zvtx:pttrig:ptasso:deltaeta:deltaphi", kTrackVariablesPair, fBins, fMin, fMax); \r
- fTHnCorrIDUNID->Sumw2();\r
-for(Int_t i=0; i<kTrackVariablesPair;i++){\r
- SetAsymmetricBin(fTHnCorrIDUNID,i,dBinsPair[i],iBinPair[i],axisTitlePair[i]);\r
- }\r
- fOutput->Add(fTHnCorrIDUNID);\r
-\r
- fTHnCorrIDUNIDmix = new THnSparseF("fTHnCorrIDUNIDmix","cent:zvtx:pttrig:ptasso:deltaeta:deltaphi", kTrackVariablesPair, fBins, fMin, fMax); \r
- fTHnCorrIDUNIDmix->Sumw2();\r
-for(Int_t i=0; i<kTrackVariablesPair;i++){\r
- SetAsymmetricBin(fTHnCorrIDUNIDmix,i,dBinsPair[i],iBinPair[i],axisTitlePair[i]);\r
- }\r
- fOutput->Add(fTHnCorrIDUNIDmix);\r
- }\r
-\r
-\r
-\r
- //ThnSparse for Correlation plots(truth MC)\r
- if((fAnalysisType == "MCAOD") && ffilltrigIDassoIDMCTRUTH) {//remember that in this case uidentified means other than pions, kaons, protons\r
- fCorrelatonTruthPrimary = new THnSparseF("fCorrelatonTruthPrimary","cent:zvtx:pttrig:ptasso:deltaeta:deltaphi", kTrackVariablesPair, fBins, fMin, fMax); \r
- fCorrelatonTruthPrimary->Sumw2();\r
-for(Int_t i=0; i<kTrackVariablesPair;i++){\r
- SetAsymmetricBin(fCorrelatonTruthPrimary,i,dBinsPair[i],iBinPair[i],axisTitlePair[i]);\r
- }\r
- fOutput->Add(fCorrelatonTruthPrimary);\r
-\r
- fCorrelatonTruthPrimarymix = new THnSparseF("fCorrelatonTruthPrimarymix","cent:zvtx:pttrig:ptasso:deltaeta:deltaphi", kTrackVariablesPair, fBins, fMin, fMax); \r
- fCorrelatonTruthPrimarymix->Sumw2();\r
-for(Int_t i=0; i<kTrackVariablesPair;i++){\r
- SetAsymmetricBin(fCorrelatonTruthPrimarymix,i,dBinsPair[i],iBinPair[i],axisTitlePair[i]);\r
- }\r
- fOutput->Add(fCorrelatonTruthPrimarymix);\r
- }\r
-\r
-\r
- //binning for trigger no. counting\r
-\r
- Double_t* fMint;\r
- Double_t* fMaxt;\r
- Int_t* fBinst;\r
- Int_t dims=3;\r
- if(fcontainPIDtrig) dims=4;\r
- fMint= new Double_t[dims];\r
- fMaxt= new Double_t[dims];\r
- fBinst= new Int_t[dims];\r
- for(Int_t i=0; i<3;i++)\r
- {\r
- fBinst[i]=iBinPair[i];\r
- fMint[i]=dBinsPair[i][0];\r
- fMaxt[i]=dBinsPair[i][iBinPair[i]];\r
- }\r
- if(fcontainPIDtrig){\r
- fBinst[3]=iBinPair[6];\r
- fMint[3]=dBinsPair[6][0];\r
- fMaxt[3]=dBinsPair[6][iBinPair[6]];\r
- }\r
-\r
- //ThSparse for trigger counting(data & reco MC)\r
- if(ffilltrigassoUNID || ffilltrigUNIDassoID || ffilltrigIDassoUNID || ffilltrigIDassoID)\r
- {\r
- fTHnTrigcount = new THnSparseF("fTHnTrigcount","cent:zvtx:pt", dims, fBinst, fMint, fMaxt); \r
- fTHnTrigcount->Sumw2();\r
-for(Int_t i=0; i<3;i++){\r
- SetAsymmetricBin(fTHnTrigcount,i,dBinsPair[i],iBinPair[i],axisTitlePair[i]);\r
- }\r
- if(fcontainPIDtrig) SetAsymmetricBin(fTHnTrigcount,3,dBinsPair[6],iBinPair[6],axisTitlePair[6]);\r
- fOutput->Add(fTHnTrigcount);\r
- }\r
- \r
- if((fAnalysisType =="MCAOD") && ffilltrigIDassoIDMCTRUTH) {\r
- //ThSparse for trigger counting(truth MC)\r
-fTHnTrigcountMCTruthPrim = new THnSparseF("fTHnTrigcountMCTruthPrim","cent:zvtx:pt", dims, fBinst, fMint, fMaxt); \r
- fTHnTrigcountMCTruthPrim->Sumw2();\r
-for(Int_t i=0; i<3;i++){\r
- SetAsymmetricBin(fTHnTrigcountMCTruthPrim,i,dBinsPair[i],iBinPair[i],axisTitlePair[i]);\r
- }\r
-if(fcontainPIDtrig) SetAsymmetricBin(fTHnTrigcountMCTruthPrim,3,dBinsPair[6],iBinPair[6],axisTitlePair[6]);\r
- fOutput->Add(fTHnTrigcountMCTruthPrim);\r
- }\r
-\r
-if(fAnalysisType=="MCAOD"){\r
- if(ffillhistQATruth)\r
- {\r
- MCtruthpt=new TH1F ("MCtruthpt","ptdistributiontruthprim",100,0.,10.);\r
- fOutput->Add(MCtruthpt);\r
-\r
- MCtrutheta=new TH1F ("MCtrutheta","etadistributiontruthprim",360,-1.8,1.8);\r
- fOutput->Add(MCtrutheta);\r
-\r
- MCtruthphi=new TH1F ("MCtruthphi","phidisttruthprim",340,0,6.8);\r
- fOutput->Add(MCtruthphi);\r
-\r
- MCtruthpionpt=new TH1F ("MCtruthpionpt","MCtruthpionpt",100,0.,10.);\r
- fOutput->Add(MCtruthpionpt);\r
-\r
- MCtruthpioneta=new TH1F ("MCtruthpioneta","MCtruthpioneta",360,-1.8,1.8);\r
- fOutput->Add(MCtruthpioneta);\r
-\r
- MCtruthpionphi=new TH1F ("MCtruthpionphi","MCtruthpionphi",340,0,6.8);\r
- fOutput->Add(MCtruthpionphi);\r
-\r
- MCtruthkaonpt=new TH1F ("MCtruthkaonpt","MCtruthkaonpt",100,0.,10.);\r
- fOutput->Add(MCtruthkaonpt);\r
-\r
- MCtruthkaoneta=new TH1F ("MCtruthkaoneta","MCtruthkaoneta",360,-1.8,1.8);\r
- fOutput->Add(MCtruthkaoneta);\r
-\r
- MCtruthkaonphi=new TH1F ("MCtruthkaonphi","MCtruthkaonphi",340,0,6.8);\r
- fOutput->Add(MCtruthkaonphi);\r
-\r
- MCtruthprotonpt=new TH1F ("MCtruthprotonpt","MCtruthprotonpt",100,0.,10.);\r
- fOutput->Add(MCtruthprotonpt);\r
-\r
- MCtruthprotoneta=new TH1F ("MCtruthprotoneta","MCtruthprotoneta",360,-1.8,1.8);\r
- fOutput->Add(MCtruthprotoneta);\r
-\r
- MCtruthprotonphi=new TH1F ("MCtruthprotonphi","MCtruthprotonphi",340,0,6.8);\r
- fOutput->Add(MCtruthprotonphi);\r
- }\r
- fPioncont=new TH2F("fPioncont", "fPioncont",10,-0.5,9.5,100,0.,10.);\r
- fOutput->Add(fPioncont);\r
-\r
- fKaoncont=new TH2F("fKaoncont","fKaoncont",10,-0.5,9.5,100,0.,10.);\r
- fOutput->Add(fKaoncont);\r
-\r
- fProtoncont=new TH2F("fProtoncont","fProtoncont",10,-0.5,9.5,100,0.,10.);\r
- fOutput->Add(fProtoncont);\r
- }\r
-\r
-//Mixing\r
-DefineEventPool();\r
-\r
- if(fapplyTrigefficiency || fapplyAssoefficiency)\r
- {\r
- const Int_t nDimt = 4;// cent zvtx pt eta\r
- Int_t fBinsCht[nDimt] = {iBinPair[0], iBinPair[1], nPteffbin ,nEtaBin};//*************change it\r
- Double_t fMinCht[nDimt] = { dBinsPair[0][0],dBinsPair[1][0], Pteff[0], EtaBin[0] };\r
- Double_t fMaxCht[nDimt] = {dBinsPair[0][iBinPair[0]], dBinsPair[1][iBinPair[1]], Pteff[nPteffbin], EtaBin[nEtaBin]};\r
-\r
- TString Histrexname;\r
-for(Int_t jj=0;jj<6;jj++)// PID type binning\r
- {\r
- Histrexname="effcorection";Histrexname+=jj;\r
- effcorection[jj] = new THnSparseF(Histrexname.Data(),"cent:zvtx::Pt:eta", nDimt, fBinsCht, fMinCht, fMaxCht);\r
- effcorection[jj]->Sumw2(); \r
- effcorection[jj]->GetAxis(0)->Set(iBinPair[0], dBinsPair[0]);\r
- effcorection[jj]->GetAxis(0)->SetTitle("Centrality");\r
- effcorection[jj]->GetAxis(1)->Set( iBinPair[1],dBinsPair[1]);\r
- effcorection[jj]->GetAxis(1)->SetTitle("zvtx"); \r
- effcorection[jj]->GetAxis(2)->Set(nPteffbin, Pteff); \r
- effcorection[jj]->GetAxis(2)->SetTitle("p_{T} (GeV/c)");\r
- effcorection[jj]->GetAxis(3)->Set( nEtaBin,EtaBin);\r
- effcorection[jj]->GetAxis(3)->SetTitle("#eta");\r
- fOutput->Add(effcorection[jj]);\r
- }\r
-// TFile *fsifile = new TFile(fefffilename,"READ");\r
- TFile *fileT=TFile::Open(fefffilename);\r
- TString Nameg;\r
-for(Int_t jj=0;jj<6;jj++)//type binning\r
- {\r
-Nameg="effmap";Nameg+=jj;\r
-//effcorection[jj] = (THnSparseF*)fsifile->Get(Nameg.Data());\r
-effcorection[jj] = (THnSparseF*)fileT->Get(Nameg.Data());\r
-\r
-//effcorection[jj]->SetDirectory(0);//****************************not present in case oh THnF\r
- }\r
-//fsifile->Close();\r
-fileT->Close();\r
-\r
- }\r
- \r
-//fControlConvResoncances = new TH2F("fControlConvResoncances", ";id;delta mass", 3, -0.5, 2.5, 100, -0.1, 0.1);\r
-// fOutput->Add(fControlConvResoncances);\r
-\r
- \r
- PostData(1, fOutput); // Post data for ALL output slots >0 here, to get at least an empty histogram\r
-}\r
-//-------------------------------------------------------------------------------\r
-void AliTwoParticlePIDCorr::UserExec( Option_t * ){\r
-\r
- \r
- if(fAnalysisType == "AOD") {\r
-\r
- doAODevent();\r
- \r
- }//AOD--analysis-----\r
-\r
- else if(fAnalysisType == "MCAOD") {\r
- \r
- doMCAODevent();\r
- \r
- }\r
- \r
- else return;\r
- \r
-}\r
-//-------------------------------------------------------------------------\r
-void AliTwoParticlePIDCorr::doMCAODevent() \r
-{\r
- AliVEvent *event = InputEvent();\r
- if (!event) { Printf("ERROR: Could not retrieve event"); return; }\r
- AliAODEvent* aod = dynamic_cast<AliAODEvent*>(event);\r
- if (!aod) {\r
- AliError("Cannot get the AOD event");\r
- return;\r
- }\r
- \r
-// count all events(physics triggered) \r
- fEventCounter->Fill(1);\r
-\r
- // get centrality object and check quality(valid for p-Pb and Pb-Pb)\r
- Double_t cent_v0=0.0;\r
-\r
- if(fSampleType=="pPb" || fSampleType=="PbPb")\r
- {\r
- AliCentrality *centrality=0;\r
- if(aod) \r
- centrality = aod->GetHeader()->GetCentralityP();\r
- // if (centrality->GetQuality() != 0) return ;\r
-\r
- if(centrality)\r
- { \r
- cent_v0 = centrality->GetCentralityPercentile(fCentralityMethod);\r
- }\r
- else\r
- {\r
- cent_v0= -1;\r
- }\r
- }\r
-\r
-//check the PIDResponse handler\r
- if (!fPID) return;\r
-\r
-// get mag. field required for twotrack efficiency cut\r
- Float_t bSign = 0;\r
- bSign = (aod->GetMagneticField() > 0) ? 1 : -1;\r
-\r
- //check for TClonesArray(truth track MC information)\r
-fArrayMC = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));\r
- if (!fArrayMC) {\r
- AliFatal("Error: MC particles branch not found!\n");\r
- return;\r
- }\r
- \r
- //check for AliAODMCHeader(truth event MC information)\r
- AliAODMCHeader *header=NULL;\r
- header=(AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName()); \r
- if(!header) {\r
- printf("MC header branch not found!\n");\r
- return;\r
- }\r
- \r
-//Only consider MC events within the vtx-z region used also as cut on the reconstructed vertex\r
-Float_t zVtxmc =header->GetVtxZ();\r
- if(TMath::Abs(zVtxmc)>fzvtxcut) return;\r
-\r
- // For productions with injected signals, figure out above which label to skip particles/tracks\r
- Int_t skipParticlesAbove = 0;\r
-\r
- if (fInjectedSignals)\r
- {\r
- AliGenEventHeader* eventHeader = 0;\r
- Int_t headers = 0;\r
-\r
-// AOD\r
- if (!header)\r
- AliFatal("fInjectedSignals set but no MC header found");\r
- \r
- headers = header->GetNCocktailHeaders();\r
- eventHeader = header->GetCocktailHeader(0);\r
-\r
- if (!eventHeader)\r
- {\r
- // We avoid AliFatal here, because the AOD productions sometimes have events where the MC header is missing \r
- // (due to unreadable Kinematics) and we don't want to loose the whole job because of a few events\r
- AliError("First event header not found. Skipping this event.");\r
- //fHistos->FillEvent(centrality, AliUEHist::kCFStepAnaTopology);\r
- return;\r
- }\r
-skipParticlesAbove = eventHeader->NProduced();\r
- AliInfo(Form("Injected signals in this event (%d headers). Keeping events of %s. Will skip particles/tracks above %d.", headers, eventHeader->ClassName(), skipParticlesAbove));\r
- }\r
-\r
- //vertex selection(is it fine for PP?)\r
- if ( fSampleType=="pPb"){\r
- trkVtx = aod->GetPrimaryVertex();\r
- if (!trkVtx || trkVtx->GetNContributors()<=0) return;\r
- TString vtxTtl = trkVtx->GetTitle();\r
- if (!vtxTtl.Contains("VertexerTracks")) return;\r
- zvtx = trkVtx->GetZ();\r
- const AliAODVertex* spdVtx = aod->GetPrimaryVertexSPD();\r
- if (!spdVtx || spdVtx->GetNContributors()<=0) return;\r
- TString vtxTyp = spdVtx->GetTitle();\r
- Double_t cov[6]={0};\r
- spdVtx->GetCovarianceMatrix(cov);\r
- Double_t zRes = TMath::Sqrt(cov[5]);\r
- if (vtxTyp.Contains("vertexer:Z") && (zRes>0.25)) return;\r
- if (TMath::Abs(spdVtx->GetZ() - trkVtx->GetZ())>0.5) return;\r
- }\r
- else if(fSampleType=="PbPb" || fSampleType=="pp") {//for pp and pb-pb case\r
- Int_t nVertex = aod->GetNumberOfVertices();\r
- if( nVertex > 0 ) { \r
- trkVtx = aod->GetPrimaryVertex();\r
- Int_t nTracksPrim = trkVtx->GetNContributors();\r
- zvtx = trkVtx->GetZ();\r
- //if (fDebug > 1)AliInfo(Form(" Vertex in = %f with %d particles by %s data ...",zVertex,nTracksPrim,vertex->GetName()));\r
- // Reject TPC only vertex\r
- TString name(trkVtx->GetName());\r
- if (name.CompareTo("PrimaryVertex") && name.CompareTo("SPDVertex"))return;\r
-\r
- // Select a quality vertex by number of tracks?\r
- if( nTracksPrim < fnTracksVertex ) {\r
- //if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");\r
- return ;\r
- }\r
- // TODO remove vertexer Z events with dispersion > 0.02: Doesn't work for AOD at present\r
- //if (strcmp(vertex->GetTitle(), "AliVertexerZ") == 0 && vertex->GetDispersion() > 0.02)\r
- // return kFALSE;\r
- // if (fDebug > 1) AliInfo(" Primary-vertex Selection: event ACCEPTED...");\r
- }\r
- else return;\r
-\r
- }\r
- else return;//as there is no proper sample type\r
-\r
-\r
- fHistQA[0]->Fill((trkVtx->GetX()));fHistQA[1]->Fill((trkVtx->GetY()));fHistQA[2]->Fill((trkVtx->GetZ())); //for trkVtx only before vertex cut |zvtx|<10 cm\r
-\r
- //count events having a proper vertex\r
- fEventCounter->Fill(2);\r
-\r
- if (TMath::Abs(zvtx) > fzvtxcut) return;\r
-\r
-fHistQA[3]->Fill((trkVtx->GetX()));fHistQA[4]->Fill((trkVtx->GetY()));fHistQA[5]->Fill((trkVtx->GetZ()));//after vertex cut for trkVtx only\r
-\r
-//now we have events passed physics trigger, centrality,eventzvtx cut \r
-\r
-//count events after vertex cut\r
- fEventCounter->Fill(5);\r
- \r
-if(!aod) return; //for safety\r
-\r
-if (fSampleType=="pPb" || fSampleType=="PbPb") if (cent_v0 < 0) return;//for pp case it is the multiplicity\r
-\r
- Double_t nooftrackstruth=0.0;//in case of pp this will give the multiplicity(for truth case) after the track loop(only for unidentified particles that pass kinematic cuts)\r
-\r
- Double_t cent_v0_truth=0.0;\r
-\r
- //TObjArray* tracksMCtruth=0;\r
-TObjArray* tracksMCtruth=new TObjArray;//for truth MC particles with PID,here unidentified means any particle other than pion, kaon or proton(Basicaly Spundefined of AliHelperPID)******WARNING::different from data and reco MC\r
- tracksMCtruth->SetOwner(kTRUE); //***********************************IMPORTANT!\r
-\r
- eventno++;\r
-\r
- //There is a small difference between zvtx and zVtxmc?????? \r
- //cout<<"***********************************************zVtxmc="<<zVtxmc<<endl;\r
- //cout<<"***********************************************zvtx="<<zvtx<<endl;\r
- \r
-//now process the truth particles(for both efficiency & correlation function)\r
-Int_t nMCTrack = fArrayMC->GetEntriesFast();\r
- \r
-for (Int_t iMC = 0; iMC < nMCTrack; iMC++) \r
-{ //MC truth track loop starts\r
- \r
-AliAODMCParticle *partMC = (AliAODMCParticle*) fArrayMC->At(iMC);\r
- \r
- if(!partMC){\r
- AliError(Form("ERROR: Could not retrieve AODMCtrack %d",iMC));\r
- continue;\r
- }\r
-\r
-//consider only charged particles\r
- if(partMC->Charge() == 0) continue;\r
-\r
-//consider only primary particles; neglect all secondary particles including from weak decays\r
- if(fselectprimaryTruth && !partMC->IsPhysicalPrimary()) continue;\r
-\r
-\r
-//remove injected signals(primaries above <maxLabel>)\r
- if (fInjectedSignals && partMC->GetLabel() >= skipParticlesAbove) continue;\r
-\r
-//remove duplicates\r
- Bool_t isduplicate=kFALSE;\r
- if (fRemoveDuplicates)\r
- { \r
- for (Int_t j=iMC+1; j<nMCTrack; ++j) \r
- {//2nd trutuh loop starts\r
-AliAODMCParticle *partMC2 = (AliAODMCParticle*) fArrayMC->At(j);\r
- if(!partMC2){\r
- AliError(Form("ERROR: Could not retrieve AODMCtrack %d",j));\r
- continue;\r
- } \r
- if (partMC->GetLabel() == partMC2->GetLabel())\r
- {\r
-isduplicate=kTRUE;\r
- break; \r
- } \r
- }//2nd truth loop ends\r
- }\r
- if(fRemoveDuplicates && isduplicate) continue;//remove duplicates\r
-\r
-//give only kinematic cuts at the truth level \r
- if (partMC->Eta() < fmineta || partMC->Eta() > fmaxeta) continue;\r
- if (partMC->Pt() < fminPt || partMC->Pt() > fmaxPt) continue;\r
-\r
- if(!partMC) continue;//for safety\r
-\r
- //To determine multiplicity in case of PP\r
- nooftrackstruth++;\r
- //cout<<"**************************************"<<TMath::Abs(partMC->GetLabel())<<endl;\r
-//only physical primary(all/unidentified) \r
-if(ffillhistQATruth)\r
- {\r
- MCtruthpt->Fill(partMC->Pt());\r
- MCtrutheta->Fill(partMC->Eta());\r
- MCtruthphi->Fill(partMC->Phi());\r
- }\r
- //get particle ID\r
-Int_t pdgtruth=((AliAODMCParticle*)partMC)->GetPdgCode();\r
-Int_t particletypeTruth=-999;\r
- if (TMath::Abs(pdgtruth)==211)\r
- {\r
- particletypeTruth=SpPion;\r
-if(ffillhistQATruth)\r
- {\r
- MCtruthpionpt->Fill(partMC->Pt());\r
- MCtruthpioneta->Fill(partMC->Eta());\r
- MCtruthpionphi->Fill(partMC->Phi());\r
- }\r
- }\r
- if (TMath::Abs(pdgtruth)==321)\r
- {\r
- particletypeTruth=SpKaon;\r
-if(ffillhistQATruth)\r
- {\r
- MCtruthkaonpt->Fill(partMC->Pt());\r
- MCtruthkaoneta->Fill(partMC->Eta());\r
- MCtruthkaonphi->Fill(partMC->Phi());\r
- }\r
- }\r
-if(TMath::Abs(pdgtruth)==2212)\r
- {\r
- particletypeTruth=SpProton;\r
-if(ffillhistQATruth)\r
- {\r
- MCtruthprotonpt->Fill(partMC->Pt());\r
- MCtruthprotoneta->Fill(partMC->Eta());\r
- MCtruthprotonphi->Fill(partMC->Phi());\r
- }\r
- }\r
- if(TMath::Abs(pdgtruth)!=211 && TMath::Abs(pdgtruth)!=321 && TMath::Abs(pdgtruth)!=2212) particletypeTruth=unidentified;//*********************WARNING:: situation is different from reco MC and data case(here we don't have SpUndefined particles,because here unidentified=SpUndefined)\r
-\r
- // -- Fill THnSparse for efficiency and contamination calculation\r
- if (fSampleType=="pp") cent_v0=15.0;//integrated over multiplicity(so put any fixed value for each track so that practically means there is only one bin in multiplicity i.e. multiplicity intregated out )**************Important\r
-\r
- Double_t primmctruth[4] = {cent_v0, zVtxmc,partMC->Pt(), partMC->Eta()};\r
- if(ffillefficiency)\r
- {\r
- fTHngenprimPidTruth[5]->Fill(primmctruth);//for all primary truth particles(4)\r
-if (TMath::Abs(pdgtruth)==211 || TMath::Abs(pdgtruth)==321) fTHngenprimPidTruth[3]->Fill(primmctruth);//for primary truth mesons(3)\r
-if (TMath::Abs(pdgtruth)==2212 || TMath::Abs(pdgtruth)==321) fTHngenprimPidTruth[4]->Fill(primmctruth);//for primary truth kaons+protons(4)\r
- if (TMath::Abs(pdgtruth)==211) fTHngenprimPidTruth[0]->Fill(primmctruth);//for pions\r
- if (TMath::Abs(pdgtruth)==321) fTHngenprimPidTruth[1]->Fill(primmctruth);//for kaons\r
- if(TMath::Abs(pdgtruth)==2212) fTHngenprimPidTruth[2]->Fill(primmctruth);//for protons\r
- }\r
- \r
- Float_t effmatrixtruth=1.0;//In Truth MC, no case of efficiency correction so it should be always 1.0\r
-if(partMC->Pt()>=fminPtAsso || partMC->Pt()<=fmaxPtTrig)//to reduce memory consumption in pool\r
- {\r
-LRCParticlePID* copy6 = new LRCParticlePID(particletypeTruth,partMC->Charge(),partMC->Pt(),partMC->Eta(), partMC->Phi(),effmatrixtruth);\r
-//copy6->SetUniqueID(eventno * 100000 + TMath::Abs(partMC->GetLabel()));\r
- copy6->SetUniqueID(eventno * 100000 + (Int_t)nooftrackstruth);\r
- tracksMCtruth->Add(copy6);//************** TObjArray used for truth correlation function calculation\r
- }\r
- }//MC truth track loop ends\r
-\r
-//*********************still in event loop\r
-\r
- if (fSampleType=="pp") cent_v0_truth=nooftrackstruth;\r
- else cent_v0_truth=cent_v0;//the notation cent_v0 is reserved for reco/data case\r
-\r
- //now cent_v0_truth should be used for all correlation function calculation\r
-\r
-if(nooftrackstruth>0.0 && ffilltrigIDassoIDMCTRUTH)\r
- {\r
- //Fill Correlations for MC truth particles(same event)\r
-if(tracksMCtruth && tracksMCtruth->GetEntriesFast()>0)//hadron triggered correlation\r
- Fillcorrelation(tracksMCtruth,0,cent_v0_truth,zVtxmc,bSign,fPtOrderMCTruth,kFALSE,kFALSE,"trigIDassoIDMCTRUTH");//mixcase=kFALSE for same event case\r
-\r
-//start mixing\r
-AliEventPool* pool2 = fPoolMgr->GetEventPool(cent_v0_truth, zVtxmc+200);\r
-if (pool2 && pool2->IsReady())\r
- {//start mixing only when pool->IsReady\r
-if(tracksMCtruth && tracksMCtruth->GetEntriesFast()>0)\r
- {//proceed only when no. of trigger particles >0 in current event\r
-for (Int_t jMix=0; jMix<pool2->GetCurrentNEvents(); jMix++) \r
- { //pool event loop start\r
- TObjArray* bgTracks6 = pool2->GetEvent(jMix);\r
- if(!bgTracks6) continue;\r
- Fillcorrelation(tracksMCtruth,bgTracks6,cent_v0_truth,zVtxmc,bSign,fPtOrderMCTruth,kFALSE,kTRUE,"trigIDassoIDMCTRUTH");//mixcase=kTRUE for mixing case\r
- \r
- }// pool event loop ends mixing case\r
- }//if(trackstrig && trackstrig->GetEntriesFast()>0) condition ends mixing case\r
-} //if pool->IsReady() condition ends mixing case\r
-\r
- //still in main event loop\r
-\r
- if(tracksMCtruth){\r
-if(pool2) pool2->UpdatePool(CloneAndReduceTrackList(tracksMCtruth));//ownership of tracksasso is with pool now, don't delete it\r
- }\r
- }\r
-\r
- //still in main event loop\r
-\r
-if(tracksMCtruth) delete tracksMCtruth;\r
-\r
-//now deal with reco tracks\r
- //TObjArray* tracksUNID=0;\r
- TObjArray* tracksUNID = new TObjArray;\r
- tracksUNID->SetOwner(kTRUE);\r
-\r
- //TObjArray* tracksID=0;\r
- TObjArray* tracksID = new TObjArray;\r
- tracksID->SetOwner(kTRUE);\r
-\r
-\r
- Float_t bSign1=aod->GetHeader()->GetMagneticField() ;//used for reconstructed track dca cut\r
-\r
- Double_t trackscount=0.0;\r
-\r
-// loop over reconstructed tracks \r
- for (Int_t itrk = 0; itrk < aod->GetNumberOfTracks(); itrk++) \r
-{ // reconstructed track loop starts\r
- AliAODTrack* track = dynamic_cast<AliAODTrack*>(aod->GetTrack(itrk));\r
- if (!track) continue;\r
- //get the corresponding MC track at the truth level (doing reco matching)\r
- AliAODMCParticle* recomatched = static_cast<AliAODMCParticle*>(fArrayMC->At(TMath::Abs(track->GetLabel())));\r
- if(!recomatched) continue;//if a reco track doesn't have corresponding truth track at generated level is a fake track(label==0), ignore it\r
-\r
-//remove injected signals \r
- if(fInjectedSignals)\r
- {\r
- AliAODMCParticle* mother = recomatched;\r
-\r
- while (!mother->IsPhysicalPrimary())\r
- {// find the primary mother;the first stable mother is searched and checked if it is <= <maxLabel>\r
- if (mother->GetMother() < 0)\r
- {\r
- mother = 0;\r
- break;\r
- }\r
- \r
- mother =(AliAODMCParticle*) fArrayMC->At(((AliAODMCParticle*)mother)->GetMother());\r
- if (!mother)\r
- break;\r
- }\r
- if (!mother)\r
- {\r
- Printf("WARNING: No mother found for particle %d:", recomatched->GetLabel());\r
- continue;\r
- }\r
- if (mother->GetLabel() >= skipParticlesAbove) continue;//remove injected signals(primaries above <maxLabel>)\r
- }//remove injected signals\r
-\r
- if (fRemoveWeakDecays && ((AliAODMCParticle*) recomatched)->IsSecondaryFromWeakDecay()) continue;//remove weak decays\r
- \r
- Bool_t isduplicate2=kFALSE;\r
-if (fRemoveDuplicates)\r
- {\r
- for (Int_t j =itrk+1; j < aod->GetNumberOfTracks(); j++) \r
- {//2nd loop starts\r
- AliAODTrack* track2 = dynamic_cast<AliAODTrack*>(aod->GetTrack(j));\r
- if (!track2) continue;\r
- AliAODMCParticle* recomatched2 = static_cast<AliAODMCParticle*>(fArrayMC->At(TMath::Abs(track2->GetLabel())));\r
-if(!recomatched2) continue;\r
-\r
-if (track->GetLabel() == track2->GetLabel())\r
- {\r
-isduplicate2=kTRUE;\r
- break; \r
- }\r
- }//2nd loop ends\r
- }\r
- if(fRemoveDuplicates && isduplicate2) continue;//remove duplicates\r
- \r
- fHistQA[11]->Fill(track->GetTPCNcls());\r
- Int_t tracktype=ClassifyTrack(track,trkVtx,bSign1);//dcacut=kFALSE,onlyprimary=kFALSE\r
-\r
- if(tracktype==0) continue; \r
- if(tracktype==1)//tracks "not" passed AliAODTrack::kPrimary at reconstructed level & have proper TPC PID response(?)\r
-{\r
- if(!track) continue;//for safety\r
- //accepted all(primaries+secondary) reconstructed tracks(pt 0.2 to 10.0,,eta -0.8 to 0.8)\r
- trackscount++;\r
-\r
-//check for eta , phi holes\r
- fEtaSpectrasso->Fill(track->Eta(),track->Pt());\r
- fphiSpectraasso->Fill(track->Phi(),track->Pt());\r
-\r
- Int_t particletypeMC=-9999;\r
-\r
-//tag all particles as unidentified\r
- particletypeMC=unidentified;\r
-\r
- Float_t effmatrix=1.;\r
-\r
-// -- Fill THnSparse for efficiency calculation\r
- if (fSampleType=="pp") cent_v0=15.0;//integrated over multiplicity(so put any fixed value for each track so that practically means there is only one bin in multiplicityi.e multiplicity intregated out )**************Important\r
- //NOTE:: this will be used for fillinfg THnSparse of efficiency & also to get the the track by track eff. factor on the fly(only in case of pp)\r
-\r
- //Clone & Reduce track list(TObjArray) for unidentified particles\r
-if(track->Pt()>=fminPtAsso || track->Pt()<=fmaxPtTrig)//to reduce memory consumption in pool\r
- {\r
- if (fapplyTrigefficiency || fapplyAssoefficiency)//get the trackingefficiency x contamination factor for unidentified particles\r
- effmatrix=GetTrackbyTrackeffvalue(track,cent_v0,zvtx,particletypeMC);\r
- LRCParticlePID* copy = new LRCParticlePID(particletypeMC,track->Charge(),track->Pt(),track->Eta(), track->Phi(),effmatrix);\r
- copy->SetUniqueID(eventno * 100000 +(Int_t)trackscount);\r
- tracksUNID->Add(copy);//track information Storage for UNID correlation function(tracks that pass the filterbit & kinematic cuts only)\r
- }\r
- Double_t allrecomatchedpid[4] = {cent_v0, zVtxmc,recomatched->Pt(), recomatched->Eta()};\r
-if(ffillefficiency) fTHnrecomatchedallPid[5]->Fill(allrecomatchedpid);//for all\r
-\r
- //now start the particle identification process:)\r
-\r
-//get the pdg code of the corresponding truth particle\r
- Int_t pdgCode = ((AliAODMCParticle*)recomatched)->GetPdgCode();\r
-\r
-Float_t dEdx = track->GetTPCsignal();\r
- fHistoTPCdEdx->Fill(track->Pt(), dEdx);\r
-\r
- if(HasTOFPID(track))\r
-{\r
-Float_t beta = GetBeta(track);\r
-fHistoTOFbeta->Fill(track->Pt(), beta);\r
- }\r
-\r
- //do track identification(nsigma method)\r
- particletypeMC=GetParticle(track);//******************************problem is here\r
-\r
-//2-d TPCTOF map(for each Pt interval)\r
- if(HasTOFPID(track)){\r
- fTPCTOFPion3d->Fill(track->Pt(),fnsigmas[SpPion][NSigmaTOF],fnsigmas[SpPion][NSigmaTPC]);\r
- fTPCTOFKaon3d->Fill(track->Pt(),fnsigmas[SpKaon][NSigmaTOF],fnsigmas[SpKaon][NSigmaTPC]);\r
- fTPCTOFProton3d->Fill(track->Pt(),fnsigmas[SpProton][NSigmaTOF],fnsigmas[SpProton][NSigmaTPC]); \r
- }\r
- if(particletypeMC==SpUndefined) continue;\r
-\r
- //Pt, Eta , Phi distribution of the reconstructed identified particles\r
-if(ffillhistQAReco)\r
- {\r
-if (particletypeMC==SpPion)\r
- {\r
- fPionPt->Fill(track->Pt());\r
- fPionEta->Fill(track->Eta());\r
- fPionPhi->Fill(track->Phi());\r
- }\r
-if (particletypeMC==SpKaon)\r
- {\r
- fKaonPt->Fill(track->Pt());\r
- fKaonEta->Fill(track->Eta());\r
- fKaonPhi->Fill(track->Phi());\r
- }\r
-if (particletypeMC==SpProton)\r
- {\r
- fProtonPt->Fill(track->Pt());\r
- fProtonEta->Fill(track->Eta());\r
- fProtonPhi->Fill(track->Phi());\r
- }\r
- }\r
-\r
- //for misidentification fraction calculation(do it with tuneonPID)\r
- if(particletypeMC==SpPion )\r
- {\r
- if(TMath::Abs(pdgCode)==211) fPioncont->Fill(1.,track->Pt());\r
- if(TMath::Abs(pdgCode)==321) fPioncont->Fill(3.,track->Pt());\r
- if(TMath::Abs(pdgCode)==2212) fPioncont->Fill(5.,track->Pt());\r
- if(TMath::Abs(pdgCode)!=211 && TMath::Abs(pdgCode)!=321 && TMath::Abs(pdgCode)!=2212) fPioncont->Fill(7.,track->Pt());\r
- }\r
-if(particletypeMC==SpKaon )\r
- {\r
- if(TMath::Abs(pdgCode)==211) fKaoncont->Fill(1.,track->Pt());\r
- if(TMath::Abs(pdgCode)==321) fKaoncont->Fill(3.,track->Pt());\r
- if(TMath::Abs(pdgCode)==2212) fKaoncont->Fill(5.,track->Pt());\r
- if(TMath::Abs(pdgCode)!=211 && TMath::Abs(pdgCode)!=321 && TMath::Abs(pdgCode)!=2212) fKaoncont->Fill(7.,track->Pt());\r
- }\r
- if(particletypeMC==SpProton )\r
- {\r
- if(TMath::Abs(pdgCode)==211) fProtoncont->Fill(1.,track->Pt());\r
- if(TMath::Abs(pdgCode)==321) fProtoncont->Fill(3.,track->Pt());\r
- if(TMath::Abs(pdgCode)==2212) fProtoncont->Fill(5.,track->Pt());\r
- if(TMath::Abs(pdgCode)!=211 && TMath::Abs(pdgCode)!=321 && TMath::Abs(pdgCode)!=2212) fProtoncont->Fill(7.,track->Pt());\r
- }\r
-\r
- //fill tracking efficiency\r
- if(ffillefficiency)\r
- {\r
- if(particletypeMC==SpPion || particletypeMC==SpKaon)\r
- {\r
-if(TMath::Abs(pdgCode)==211 || TMath::Abs(pdgCode)==321) fTHnrecomatchedallPid[3]->Fill(allrecomatchedpid);//for mesons\r
- }\r
- if(particletypeMC==SpKaon || particletypeMC==SpProton)\r
- {\r
-if(TMath::Abs(pdgCode)==321 || TMath::Abs(pdgCode)==2212) fTHnrecomatchedallPid[4]->Fill(allrecomatchedpid);//for kaons+protons\r
- }\r
- if(particletypeMC==SpPion && TMath::Abs(pdgCode)==211) fTHnrecomatchedallPid[0]->Fill(allrecomatchedpid);//for pions \r
- if(particletypeMC==SpKaon && TMath::Abs(pdgCode)==321) fTHnrecomatchedallPid[1]->Fill(allrecomatchedpid);//for kaons\r
- if(particletypeMC==SpProton && TMath::Abs(pdgCode)==2212) fTHnrecomatchedallPid[2]->Fill(allrecomatchedpid);//for protons\r
- }\r
-\r
-if(track->Pt()>=fminPtAsso || track->Pt()<=fmaxPtTrig)//to reduce memory consumption in pool\r
- {\r
-if (fapplyTrigefficiency || fapplyAssoefficiency)\r
- effmatrix=GetTrackbyTrackeffvalue(track,cent_v0,zvtx,particletypeMC);//get the tracking eff x TOF matching eff x PID eff x contamination factor for identified particles \r
- LRCParticlePID* copy1 = new LRCParticlePID(particletypeMC,track->Charge(),track->Pt(),track->Eta(), track->Phi(),effmatrix);\r
- copy1->SetUniqueID(eventno * 100000 + (Int_t)trackscount);\r
- tracksID->Add(copy1);\r
- }\r
- }// if(tracktype==1) condition structure ands\r
-\r
-}//reco track loop ends\r
-\r
- //*************************************************************still in event loop\r
- \r
-//same event delta-eta-deltaphi plot \r
-if(fSampleType=="pp") cent_v0=trackscount;//multiplicity\r
-\r
-if(trackscount>0.0)\r
- { \r
-//fill the centrality/multiplicity distribution of the selected events\r
- fhistcentrality->Fill(cent_v0);//*********************************WARNING::binning of cent_v0 is different for pp and pPb/PbPb case\r
-\r
- if (fSampleType=="pPb" || fSampleType=="PbPb") fCentralityCorrelation->Fill(cent_v0, trackscount);//only with unidentified tracks(i.e before PID selection);;;;;can be used to remove centrality outliers??????\r
-\r
-//count selected events having centrality betn 0-100%\r
- fEventCounter->Fill(6);\r
-\r
- //same event delte-eta, delta-phi plot\r
-if(tracksUNID && tracksUNID->GetEntriesFast()>0)//hadron triggered correlation\r
- {//same event calculation starts\r
- if(ffilltrigassoUNID) Fillcorrelation(tracksUNID,0,cent_v0,zvtx,bSign,kTRUE,kTRUE,kFALSE,"trigassoUNID");//mixcase=kFALSE (hadron-hadron correlation)\r
- if(tracksID && tracksID->GetEntriesFast()>0 && ffilltrigUNIDassoID) Fillcorrelation(tracksUNID,tracksID,cent_v0,zvtx,bSign,kFALSE,kTRUE,kFALSE,"trigUNIDassoID");//mixcase=kFALSE (hadron-ID correlation)\r
- }\r
-\r
-if(tracksID && tracksID->GetEntriesFast()>0)//ID triggered correlation\r
- {//same event calculation starts\r
- if(tracksUNID && tracksUNID->GetEntriesFast()>0 && ffilltrigIDassoUNID) Fillcorrelation(tracksID,tracksUNID,cent_v0,zvtx,bSign,kFALSE,kTRUE,kFALSE,"trigIDassoUNID");//mixcase=kFALSE (ID-hadron correlation)\r
- if(ffilltrigIDassoID) Fillcorrelation(tracksID,0,cent_v0,zvtx,bSign,kFALSE,kTRUE,kFALSE,"trigIDassoID");//mixcase=kFALSE (ID-ID correlation)\r
- }\r
-\r
-//still in main event loop\r
-//start mixing\r
- if(ffilltrigassoUNID || ffilltrigIDassoUNID){//mixing with unidentified particles\r
-AliEventPool* pool = fPoolMgr->GetEventPool(cent_v0, zvtx);//In the pool there is tracksUNID(i.e associateds are unidentified)\r
-if (pool && pool->IsReady())\r
- {//start mixing only when pool->IsReady\r
- for (Int_t jMix=0; jMix<pool->GetCurrentNEvents(); jMix++) \r
- { //pool event loop start\r
- TObjArray* bgTracks = pool->GetEvent(jMix);\r
- if(!bgTracks) continue;\r
- if(ffilltrigassoUNID && tracksUNID && tracksUNID->GetEntriesFast()>0)//*******************************hadron trggered mixing\r
- Fillcorrelation(tracksUNID,bgTracks,cent_v0,zvtx,bSign,kTRUE,kTRUE,kTRUE,"trigassoUNID");//mixcase=kTRUE\r
- if(ffilltrigIDassoUNID && tracksID && tracksID->GetEntriesFast()>0)//***********************************ID trggered mixing\r
- Fillcorrelation(tracksID,bgTracks,cent_v0,zvtx,bSign,kFALSE,kTRUE,kTRUE,"trigIDassoUNID");//mixcase=kTRUE \r
- }// pool event loop ends mixing case\r
-\r
-} //if pool->IsReady() condition ends mixing case\r
- if(tracksUNID) {\r
-if(pool)\r
- pool->UpdatePool(CloneAndReduceTrackList(tracksUNID));\r
- }\r
- }//mixing with unidentified particles\r
-\r
- if(ffilltrigUNIDassoID || ffilltrigIDassoID){//mixing with identified particles\r
-AliEventPool* pool1 = fPoolMgr->GetEventPool(cent_v0, zvtx+100);//In the pool1 there is tracksID(i.e associateds are identified)\r
-if (pool1 && pool1->IsReady())\r
- {//start mixing only when pool->IsReady\r
-for (Int_t jMix=0; jMix<pool1->GetCurrentNEvents(); jMix++) \r
- { //pool event loop start\r
- TObjArray* bgTracks2 = pool1->GetEvent(jMix);\r
- if(!bgTracks2) continue;\r
-if(ffilltrigUNIDassoID && tracksUNID && tracksUNID->GetEntriesFast()>0)\r
- Fillcorrelation(tracksUNID,bgTracks2,cent_v0,zvtx,bSign,kFALSE,kTRUE,kTRUE,"trigUNIDassoID");//mixcase=kTRUE \r
-if(ffilltrigIDassoID && tracksID && tracksID->GetEntriesFast()>0)\r
- Fillcorrelation(tracksID,bgTracks2,cent_v0,zvtx,bSign,kFALSE,kTRUE,kTRUE,"trigIDassoID");//mixcase=kTRUE\r
-\r
- }// pool event loop ends mixing case\r
-} //if pool1->IsReady() condition ends mixing case\r
-\r
-if(tracksID) {\r
-if(pool1) \r
- pool1->UpdatePool(CloneAndReduceTrackList(tracksID));//ownership of tracksasso is with pool now, don't delete it(tracksUNID is with pool)\r
- }\r
- }//mixing with identified particles\r
-\r
- //no. of events analyzed\r
-fEventCounter->Fill(7);\r
- }\r
-\r
-if(tracksUNID) delete tracksUNID;\r
-\r
-if(tracksID) delete tracksID;\r
-\r
-\r
-PostData(1, fOutput);\r
-\r
-}\r
-//________________________________________________________________________\r
-void AliTwoParticlePIDCorr::doAODevent() \r
-{\r
-\r
- //get AOD\r
- AliVEvent *event = InputEvent();\r
- if (!event) { Printf("ERROR: Could not retrieve event"); return; }\r
- AliAODEvent* aod = dynamic_cast<AliAODEvent*>(event);\r
- if (!aod) {\r
- AliError("Cannot get the AOD event");\r
- return;\r
- }\r
-\r
-// count all events \r
- fEventCounter->Fill(1);\r
-\r
- // get centrality object and check quality\r
- Double_t cent_v0=0;\r
-\r
- if(fSampleType=="pPb" || fSampleType=="PbPb")\r
- {\r
- AliCentrality *centrality=0;\r
- if(aod) \r
- centrality = aod->GetHeader()->GetCentralityP();\r
- // if (centrality->GetQuality() != 0) return ;\r
-\r
- if(centrality)\r
- { \r
- cent_v0 = centrality->GetCentralityPercentile(fCentralityMethod);\r
- }\r
- else\r
- {\r
- cent_v0= -1;\r
- }\r
- }\r
-\r
- Float_t bSign = (aod->GetMagneticField() > 0) ? 1 : -1;//for two track efficiency cut in correlation function calculation\r
- Float_t bSign1=aod->GetHeader()->GetMagneticField() ;//for dca cut in ClassifyTrack(), i.e in track loop\r
-\r
-// Pileup selection ************************************************\r
-// if(frejectPileUp) //applicable only for TPC only tracks,not for hybrid tracks?.\r
- // {\r
- //if (fAnalysisUtils && fAnalysisUtils->IsPileUpEvent(aod)) return;\r
- // }\r
-\r
- if (!fPID) return;//this should be available with each event even if we don't do PID selection\r
- \r
- //vertex selection(is it fine for PP?)\r
- if ( fSampleType=="pPb"){\r
- trkVtx = aod->GetPrimaryVertex();\r
- if (!trkVtx || trkVtx->GetNContributors()<=0) return;\r
- TString vtxTtl = trkVtx->GetTitle();\r
- if (!vtxTtl.Contains("VertexerTracks")) return;\r
- zvtx = trkVtx->GetZ();\r
- const AliAODVertex* spdVtx = aod->GetPrimaryVertexSPD();\r
- if (!spdVtx || spdVtx->GetNContributors()<=0) return;\r
- TString vtxTyp = spdVtx->GetTitle();\r
- Double_t cov[6]={0};\r
- spdVtx->GetCovarianceMatrix(cov);\r
- Double_t zRes = TMath::Sqrt(cov[5]);\r
- if (vtxTyp.Contains("vertexer:Z") && (zRes>0.25)) return;\r
- if (TMath::Abs(spdVtx->GetZ() - trkVtx->GetZ())>0.5) return;\r
- }\r
- else if(fSampleType=="PbPb" || fSampleType=="pp") {//for pp and pb-pb case\r
- Int_t nVertex = aod->GetNumberOfVertices();\r
- if( nVertex > 0 ) { \r
- trkVtx = aod->GetPrimaryVertex();\r
- Int_t nTracksPrim = trkVtx->GetNContributors();\r
- zvtx = trkVtx->GetZ();\r
- //if (fDebug > 1)AliInfo(Form(" Vertex in = %f with %d particles by %s data ...",zVertex,nTracksPrim,vertex->GetName()));\r
- // Reject TPC only vertex\r
- TString name(trkVtx->GetName());\r
- if (name.CompareTo("PrimaryVertex") && name.CompareTo("SPDVertex"))return;\r
-\r
- // Select a quality vertex by number of tracks?\r
- if( nTracksPrim < fnTracksVertex ) {\r
- //if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");\r
- return ;\r
- }\r
- // TODO remove vertexer Z events with dispersion > 0.02: Doesn't work for AOD at present\r
- //if (strcmp(vertex->GetTitle(), "AliVertexerZ") == 0 && vertex->GetDispersion() > 0.02)\r
- // return kFALSE;\r
- // if (fDebug > 1) AliInfo(" Primary-vertex Selection: event ACCEPTED...");\r
- }\r
- else return;\r
-\r
- }\r
- else return;//as there is no proper sample type\r
-\r
-\r
- fHistQA[0]->Fill((trkVtx->GetX()));fHistQA[1]->Fill((trkVtx->GetY()));fHistQA[2]->Fill((trkVtx->GetZ())); //for trkVtx only before vertex cut |zvtx|<10 cm\r
-\r
-//count events having a proper vertex\r
- fEventCounter->Fill(2);\r
-\r
- if (TMath::Abs(zvtx) > fzvtxcut) return;\r
-\r
-//count events after vertex cut\r
- fEventCounter->Fill(5);\r
-\r
-\r
- //if(!fAnalysisUtils->IsVertexSelected2013pA(aod)) return;\r
- \r
- fHistQA[3]->Fill((trkVtx->GetX()));fHistQA[4]->Fill((trkVtx->GetY()));fHistQA[5]->Fill((trkVtx->GetZ()));//after vertex cut,for trkVtx only\r
-\r
-\r
- if(!aod) return; //for safety\r
- \r
-if(fSampleType=="pPb" || fSampleType=="PbPb") if (cent_v0 < 0) return;//for pp case it is the multiplicity\r
-\r
- TObjArray* tracksUNID= new TObjArray;//track info before doing PID\r
- tracksUNID->SetOwner(kTRUE); // IMPORTANT!\r
-\r
- TObjArray* tracksID= new TObjArray;//only pions, kaons,protonsi.e. after doing the PID selection\r
- tracksID->SetOwner(kTRUE); // IMPORTANT!\r
- \r
- Double_t trackscount=0.0;//in case of pp this will give the multiplicity after the track loop(only for unidentified particles that pass the filterbit and kinematic cuts)\r
-\r
- eventno++;\r
-\r
- for (Int_t itrk = 0; itrk < aod->GetNumberOfTracks(); itrk++) \r
-{ //track loop starts for TObjArray(containing track and event information) filling; used for correlation function calculation \r
- AliAODTrack* track = dynamic_cast<AliAODTrack*>(aod->GetTrack(itrk));\r
- if (!track) continue;\r
- fHistQA[11]->Fill(track->GetTPCNcls());\r
- Int_t particletype=-9999;//required for PID filling\r
- Int_t tracktype=ClassifyTrack(track,trkVtx,bSign1);//dcacut=kFALSE,onlyprimary=kFALSE\r
- if(tracktype!=1) continue; \r
-\r
- if(!track) continue;//for safety\r
-\r
-//check for eta , phi holes\r
- fEtaSpectrasso->Fill(track->Eta(),track->Pt());\r
- fphiSpectraasso->Fill(track->Phi(),track->Pt());\r
-\r
- trackscount++;\r
-\r
- //if no applyefficiency , set the eff factor=1.0\r
- Float_t effmatrix=1.0;\r
-\r
- //tag all particles as unidentified that passed filterbit & kinematic cuts \r
- particletype=unidentified;\r
-\r
-\r
- if (fSampleType=="pp") cent_v0=15.0;//integrated over multiplicity [i.e each track has multiplicity 15.0](so put any fixed value for each track so that practically means there is only one bin in multiplicityi.e multiplicity intregated out )**************Important for efficiency related issues\r
-\r
-\r
- //to reduce memory consumption in pool\r
- if(track->Pt()>=fminPtAsso || track->Pt()<=fmaxPtTrig) \r
- {\r
- //Clone & Reduce track list(TObjArray) for unidentified particles\r
- if (fapplyTrigefficiency || fapplyAssoefficiency)//get the trackingefficiency x contamination factor for unidentified particles\r
- effmatrix=GetTrackbyTrackeffvalue(track,cent_v0,zvtx,particletype);\r
- LRCParticlePID* copy = new LRCParticlePID(particletype,track->Charge(),track->Pt(),track->Eta(), track->Phi(),effmatrix);\r
- copy->SetUniqueID(eventno * 100000 + (Int_t)trackscount);\r
- tracksUNID->Add(copy);//track information Storage for UNID correlation function(tracks that pass the filterbit & kinematic cuts only)\r
- }\r
-\r
-//now start the particle identificaion process:) \r
-\r
-//track passing filterbit 768 have proper TPC response,or need to be checked explicitly before doing PID????\r
-\r
- Float_t dEdx = track->GetTPCsignal();\r
- fHistoTPCdEdx->Fill(track->Pt(), dEdx);\r
-\r
- //fill beta vs Pt plots only for tracks having proper TOF response(much less tracks compared to the no. that pass the filterbit & kinematic cuts)\r
- if(HasTOFPID(track))\r
-{\r
- Float_t beta = GetBeta(track);\r
- fHistoTOFbeta->Fill(track->Pt(), beta);\r
- }\r
-\r
-\r
-//track identification(using nsigma method)\r
- particletype=GetParticle(track);//*******************************change may be required(It should return only pion,kaon, proton and Spundefined; NOT unidentifed***************be careful)\r
-\r
-\r
-//2-d TPCTOF map(for each Pt interval)\r
- if(HasTOFPID(track)){\r
- fTPCTOFPion3d->Fill(track->Pt(),fnsigmas[SpPion][NSigmaTOF],fnsigmas[SpPion][NSigmaTPC]);\r
- fTPCTOFKaon3d->Fill(track->Pt(),fnsigmas[SpKaon][NSigmaTOF],fnsigmas[SpKaon][NSigmaTPC]);\r
- fTPCTOFProton3d->Fill(track->Pt(),fnsigmas[SpProton][NSigmaTOF],fnsigmas[SpProton][NSigmaTPC]); \r
- }\r
-\r
-//ignore the Spundefined particles as they also contain pion, kaon, proton outside the nsigma cut(also if tracks don't have proper TOF PID in a certain Pt interval) & these tracks are actually counted when we do the the efficiency correction, so considering them as unidentified particles & doing the efficiency correction(i.e defining unidentified=pion+Kaon+proton+SpUndefined is right only without efficiency correction) for them will be two times wrong!!!!! \r
- if (particletype==SpUndefined) continue;//this condition creating a modulated structure in delphi projection in mixed event case(only when we are dealing with identified particles i.e. tracksID)!!!!!!!!!!!\r
-\r
- //Pt, Eta , Phi distribution of the reconstructed identified particles\r
-if(ffillhistQAReco)\r
- {\r
-if (particletype==SpPion)\r
- {\r
- fPionPt->Fill(track->Pt());\r
- fPionEta->Fill(track->Eta());\r
- fPionPhi->Fill(track->Phi());\r
- }\r
-if (particletype==SpKaon)\r
- {\r
- fKaonPt->Fill(track->Pt());\r
- fKaonEta->Fill(track->Eta());\r
- fKaonPhi->Fill(track->Phi());\r
- }\r
-if (particletype==SpProton)\r
- {\r
- fProtonPt->Fill(track->Pt());\r
- fProtonEta->Fill(track->Eta());\r
- fProtonPhi->Fill(track->Phi());\r
- }\r
- }\r
- \r
-if(track->Pt()>=fminPtAsso || track->Pt()<=fmaxPtTrig) \r
- {\r
-if (fapplyTrigefficiency || fapplyAssoefficiency)\r
- effmatrix=GetTrackbyTrackeffvalue(track,cent_v0,zvtx,particletype);//get the tracking eff x TOF matching eff x PID eff x contamination factor for identified particles; Bool_t mesoneffrequired=kFALSE\r
- LRCParticlePID* copy1 = new LRCParticlePID(particletype,track->Charge(),track->Pt(),track->Eta(), track->Phi(),effmatrix);\r
- copy1->SetUniqueID(eventno * 100000 + (Int_t)trackscount);\r
- tracksID->Add(copy1);\r
- }\r
-} //track loop ends but still in event loop\r
-\r
-if(trackscount<1.0){\r
- if(tracksUNID) delete tracksUNID;\r
- if(tracksID) delete tracksID;\r
- return;\r
- }\r
-\r
-// cout<<fminPtAsso<<"***"<<fmaxPtAsso<<"*********************"<<fminPtTrig<<"***"<<fmaxPtTrig<<"*****************"<<kTrackVariablesPair<<endl;\r
-if(fSampleType=="pp") cent_v0=trackscount;//multiplicity\r
- \r
-//fill the centrality/multiplicity distribution of the selected events\r
- fhistcentrality->Fill(cent_v0);//*********************************WARNING::binning of cent_v0 is different for pp and pPb/PbPb case\r
-\r
-if(fSampleType=="pPb" || fSampleType=="PbPb") fCentralityCorrelation->Fill(cent_v0, trackscount);//only with unidentified tracks(i.e before PID selection);;;;;can be used to remove centrality outliers??????\r
-\r
-//count selected events having centrality betn 0-100%\r
- fEventCounter->Fill(6);\r
-\r
-//same event delta-eta-deltaphi plot \r
-\r
-if(tracksUNID && tracksUNID->GetEntriesFast()>0)//hadron triggered correlation\r
- {//same event calculation starts\r
- if(ffilltrigassoUNID) Fillcorrelation(tracksUNID,0,cent_v0,zvtx,bSign,kTRUE,kTRUE,kFALSE,"trigassoUNID");//mixcase=kFALSE (hadron-hadron correlation)\r
- if(tracksID && tracksID->GetEntriesFast()>0 && ffilltrigUNIDassoID) Fillcorrelation(tracksUNID,tracksID,cent_v0,zvtx,bSign,kFALSE,kTRUE,kFALSE,"trigUNIDassoID");//mixcase=kFALSE (hadron-ID correlation)\r
- }\r
-\r
-if(tracksID && tracksID->GetEntriesFast()>0)//ID triggered correlation\r
- {//same event calculation starts\r
- if(tracksUNID && tracksUNID->GetEntriesFast()>0 && ffilltrigIDassoUNID) Fillcorrelation(tracksID,tracksUNID,cent_v0,zvtx,bSign,kFALSE,kTRUE,kFALSE,"trigIDassoUNID");//mixcase=kFALSE (ID-hadron correlation)\r
- if(ffilltrigIDassoID) Fillcorrelation(tracksID,0,cent_v0,zvtx,bSign,kFALSE,kTRUE,kFALSE,"trigIDassoID");//mixcase=kFALSE (ID-ID correlation)\r
- }\r
-\r
-//still in main event loop\r
-\r
-\r
-//start mixing\r
- if(ffilltrigassoUNID || ffilltrigIDassoUNID){//mixing with unidentified particles\r
-AliEventPool* pool = fPoolMgr->GetEventPool(cent_v0, zvtx);//In the pool there is tracksUNID(i.e associateds are unidentified)\r
-if (pool && pool->IsReady())\r
- {//start mixing only when pool->IsReady\r
- for (Int_t jMix=0; jMix<pool->GetCurrentNEvents(); jMix++) \r
- { //pool event loop start\r
- TObjArray* bgTracks = pool->GetEvent(jMix);\r
- if(!bgTracks) continue;\r
- if(ffilltrigassoUNID && tracksUNID && tracksUNID->GetEntriesFast()>0)//*******************************hadron trggered mixing\r
- Fillcorrelation(tracksUNID,bgTracks,cent_v0,zvtx,bSign,kTRUE,kTRUE,kTRUE,"trigassoUNID");//mixcase=kTRUE\r
- if(ffilltrigIDassoUNID && tracksID && tracksID->GetEntriesFast()>0)//***********************************ID trggered mixing\r
- Fillcorrelation(tracksID,bgTracks,cent_v0,zvtx,bSign,kFALSE,kTRUE,kTRUE,"trigIDassoUNID");//mixcase=kTRUE \r
- }// pool event loop ends mixing case\r
-\r
-} //if pool->IsReady() condition ends mixing case\r
- if(tracksUNID) {\r
-if(pool)\r
- pool->UpdatePool(CloneAndReduceTrackList(tracksUNID));\r
- }\r
- }//mixing with unidentified particles\r
-\r
-\r
- if(ffilltrigUNIDassoID || ffilltrigIDassoID){//mixing with identified particles\r
-AliEventPool* pool1 = fPoolMgr->GetEventPool(cent_v0, zvtx+100);//In the pool1 there is tracksID(i.e associateds are identified)\r
-if (pool1 && pool1->IsReady())\r
- {//start mixing only when pool->IsReady\r
-for (Int_t jMix=0; jMix<pool1->GetCurrentNEvents(); jMix++) \r
- { //pool event loop start\r
- TObjArray* bgTracks2 = pool1->GetEvent(jMix);\r
- if(!bgTracks2) continue;\r
-if(ffilltrigUNIDassoID && tracksUNID && tracksUNID->GetEntriesFast()>0)\r
- Fillcorrelation(tracksUNID,bgTracks2,cent_v0,zvtx,bSign,kFALSE,kTRUE,kTRUE,"trigUNIDassoID");//mixcase=kTRUE \r
-if(ffilltrigIDassoID && tracksID && tracksID->GetEntriesFast()>0)\r
- Fillcorrelation(tracksID,bgTracks2,cent_v0,zvtx,bSign,kFALSE,kTRUE,kTRUE,"trigIDassoID");//mixcase=kTRUE\r
-\r
- }// pool event loop ends mixing case\r
-} //if pool1->IsReady() condition ends mixing case\r
-\r
-if(tracksID) {\r
-if(pool1) \r
- pool1->UpdatePool(CloneAndReduceTrackList(tracksID));//ownership of tracksasso is with pool now, don't delete it(tracksUNID is with pool)\r
- }\r
- }//mixing with identified particles\r
-\r
-\r
- //no. of events analyzed\r
-fEventCounter->Fill(7);\r
- \r
-//still in main event loop\r
-\r
-\r
-if(tracksUNID) delete tracksUNID;\r
-\r
-if(tracksID) delete tracksID;\r
-\r
-\r
-PostData(1, fOutput);\r
-\r
-} // *************************event loop ends******************************************//_______________________________________________________________________\r
-TObjArray* AliTwoParticlePIDCorr::CloneAndReduceTrackList(TObjArray* tracks)\r
-{\r
- // clones a track list by using AliDPhiBasicParticle which uses much less memory (used for event mixing)\r
- \r
- TObjArray* tracksClone = new TObjArray;\r
- tracksClone->SetOwner(kTRUE);\r
- \r
- for (Int_t i=0; i<tracks->GetEntriesFast(); i++)\r
- {\r
- LRCParticlePID* particle = (LRCParticlePID*) tracks->UncheckedAt(i);\r
- LRCParticlePID* copy100 = new LRCParticlePID(particle->getparticle(),particle->Charge(), particle->Pt(),particle->Eta(), particle->Phi(), particle->geteffcorrectionval());\r
- copy100->SetUniqueID(particle->GetUniqueID());\r
- tracksClone->Add(copy100);\r
- }\r
- \r
- return tracksClone;\r
-}\r
-\r
-//--------------------------------------------------------------------------------\r
-void AliTwoParticlePIDCorr::Fillcorrelation(TObjArray *trackstrig,TObjArray *tracksasso,Double_t cent,Float_t vtx,Float_t bSign,Bool_t fPtOrder,Bool_t twoTrackEfficiencyCut,Bool_t mixcase,TString fillup)\r
-{\r
-\r
- //before calling this function check that either trackstrig & tracksasso are available \r
-\r
- // Eta() is extremely time consuming, therefore cache it for the inner loop here:\r
- TObjArray* input = (tracksasso) ? tracksasso : trackstrig;\r
- TArrayF eta(input->GetEntriesFast());\r
- for (Int_t i=0; i<input->GetEntriesFast(); i++)\r
- eta[i] = ((LRCParticlePID*) input->UncheckedAt(i))->Eta();\r
-\r
- //if(trackstrig)\r
- Int_t jmax=trackstrig->GetEntriesFast();\r
- if(tracksasso)\r
- jmax=tracksasso->GetEntriesFast();\r
-\r
-// identify K, Lambda candidates and flag those particles\r
- // a TObject bit is used for this\r
-const UInt_t kResonanceDaughterFlag = 1 << 14;\r
- if (fRejectResonanceDaughters > 0)\r
- {\r
- Double_t resonanceMass = -1;\r
- Double_t massDaughter1 = -1;\r
- Double_t massDaughter2 = -1;\r
- const Double_t interval = 0.02;\r
- switch (fRejectResonanceDaughters)\r
- {\r
- case 1: resonanceMass = 0.9; massDaughter1 = 0.1396; massDaughter2 = 0.9383; break; // method test\r
- case 2: resonanceMass = 0.4976; massDaughter1 = 0.1396; massDaughter2 = massDaughter1; break; // k0\r
- case 3: resonanceMass = 1.115; massDaughter1 = 0.1396; massDaughter2 = 0.9383; break; // lambda\r
- default: AliFatal(Form("Invalid setting %d", fRejectResonanceDaughters));\r
- } \r
-\r
-for (Int_t i=0; i<trackstrig->GetEntriesFast(); i++)\r
- trackstrig->UncheckedAt(i)->ResetBit(kResonanceDaughterFlag);\r
- for (Int_t i=0; tracksasso->GetEntriesFast(); i++)\r
- tracksasso->UncheckedAt(i)->ResetBit(kResonanceDaughterFlag);\r
-\r
- for (Int_t i=0; i<trackstrig->GetEntriesFast(); i++)\r
- {\r
- LRCParticlePID *trig=(LRCParticlePID*)(trackstrig->UncheckedAt(i));\r
-for (Int_t j=0; tracksasso->GetEntriesFast(); j++)\r
- {\r
- LRCParticlePID *asso=(LRCParticlePID*)(tracksasso->UncheckedAt(j));\r
-\r
- // check if both particles point to the same element (does not occur for mixed events, but if subsets are mixed within the same event)\r
-if (trig->IsEqual(asso)) continue;\r
-\r
-if (trig->Charge() * asso->Charge() > 0) continue;\r
-\r
- Float_t mass = GetInvMassSquaredCheap(trig->Pt(), trig->Eta(), trig->Phi(), asso->Pt(), asso->Eta(), asso->Phi(), massDaughter1, massDaughter2);\r
- \r
-if (TMath::Abs(mass - resonanceMass*resonanceMass) < interval*5)\r
- {\r
- mass = GetInvMassSquared(trig->Pt(), trig->Eta(), trig->Phi(), asso->Pt(), asso->Eta(), asso->Phi(), massDaughter1, massDaughter2);\r
-\r
- if (mass > (resonanceMass-interval)*(resonanceMass-interval) && mass < (resonanceMass+interval)*(resonanceMass+interval))\r
- {\r
- trig->SetBit(kResonanceDaughterFlag);\r
- asso->SetBit(kResonanceDaughterFlag);\r
- \r
-// Printf("Flagged %d %d %f", i, j, TMath::Sqrt(mass));\r
- }\r
- }\r
- }\r
- }\r
- }\r
-\r
- //two particle correlation filling\r
-\r
-for(Int_t i=0;i<trackstrig->GetEntriesFast();i++)\r
- { //trigger loop starts\r
- LRCParticlePID *trig=(LRCParticlePID*)(trackstrig->UncheckedAt(i));\r
- if(!trig) continue;\r
- Float_t trigpt=trig->Pt();\r
- //to avoid overflow qnd underflow\r
- if(trigpt<fminPtTrig || trigpt>fmaxPtTrig) continue;\r
- Int_t particlepidtrig=trig->getparticle();\r
- if(fTriggerSpeciesSelection){ if (particlepidtrig!=fTriggerSpecies) continue;}\r
-\r
- Float_t trigeta=trig->Eta();\r
-\r
- // some optimization\r
- if (fTriggerRestrictEta > 0 && TMath::Abs(trigeta) > fTriggerRestrictEta)\r
- continue;\r
-\r
-if (fOnlyOneEtaSide != 0)\r
- {\r
- if (fOnlyOneEtaSide * trigeta < 0)\r
- continue;\r
- }\r
- if (fTriggerSelectCharge != 0)\r
- if (trig->Charge() * fTriggerSelectCharge < 0)\r
- continue;\r
- \r
- if (fRejectResonanceDaughters > 0)\r
- if (trig->TestBit(kResonanceDaughterFlag)) continue;\r
-\r
- Float_t trigphi=trig->Phi();\r
- Float_t trackefftrig=1.0;\r
- if(fapplyTrigefficiency) trackefftrig=trig->geteffcorrectionval();\r
- // cout<<"*******************trackefftrig="<<trackefftrig<<endl;\r
- Double_t* trigval;\r
- Int_t dim=3;\r
- if(fcontainPIDtrig) dim=4;\r
- trigval= new Double_t[dim];\r
- trigval[0] = cent;\r
- trigval[1] = vtx;\r
- trigval[2] = trigpt;\r
-if(fcontainPIDtrig) trigval[3] = particlepidtrig;\r
-\r
- //filling only for same event case(mixcase=kFALSE)\r
-\r
- //trigger particle counting only when mixcase=kFALSE i.e for same event correlation function calculation\r
-if(mixcase==kFALSE) \r
- {\r
- if(ffilltrigassoUNID==kTRUE && ffilltrigUNIDassoID==kTRUE){\r
- if(fillup=="trigassoUNID" ) fTHnTrigcount->Fill(trigval,1.0/trackefftrig); \r
- }\r
- if(ffilltrigassoUNID==kTRUE && ffilltrigUNIDassoID==kFALSE){\r
- if(fillup=="trigassoUNID" ) fTHnTrigcount->Fill(trigval,1.0/trackefftrig); \r
- }\r
-if(ffilltrigassoUNID==kFALSE && ffilltrigUNIDassoID==kTRUE){\r
- if(fillup=="trigUNIDassoID") fTHnTrigcount->Fill(trigval,1.0/trackefftrig); \r
- }\r
- //ensure that trigIDassoID , trigassoUNID, trigIDassoUNID & trigUNIDassoID case FillCorrelation called only once in the event loop for same event correlation function calculation, otherwise there will be multiple counting of pion, kaon,proton,unidentified\r
-if(ffilltrigIDassoUNID==kTRUE && ffilltrigIDassoID==kTRUE){\r
- if(fillup=="trigIDassoID") fTHnTrigcount->Fill(trigval,1.0/trackefftrig); \r
- }\r
- if(ffilltrigIDassoUNID==kTRUE && ffilltrigIDassoID==kFALSE){\r
- if(fillup=="trigIDassoUNID" ) fTHnTrigcount->Fill(trigval,1.0/trackefftrig); \r
- }\r
-if(ffilltrigIDassoUNID==kFALSE && ffilltrigIDassoID==kTRUE){\r
- if(fillup=="trigIDassoID") fTHnTrigcount->Fill(trigval,1.0/trackefftrig); \r
- }\r
-\r
- if(fillup=="trigIDassoIDMCTRUTH") fTHnTrigcountMCTruthPrim->Fill(trigval,1.0/trackefftrig); //In truth MC case "Unidentified" means any particle other than pion,kaon or proton and no efficiency correction(default value 1.0)************************be careful!!!! \r
- }\r
-\r
- //asso loop starts within trigger loop\r
- for(Int_t j=0;j<jmax;j++)\r
- {\r
- LRCParticlePID *asso=0;\r
- if(!tracksasso)\r
- asso=(LRCParticlePID*)(trackstrig->UncheckedAt(j));\r
- else\r
- asso=(LRCParticlePID*)(tracksasso->UncheckedAt(j));\r
-\r
- if(!asso) continue;\r
-\r
- //to avoid overflow qnd underflow\r
- if(asso->Pt()<fminPtAsso || asso->Pt()>fmaxPtAsso) continue;//***********************Important\r
-\r
- if(fmaxPtAsso==fminPtTrig) {if(asso->Pt()==fminPtTrig) continue;}//******************Think about it!\r
-\r
- if(!tracksasso && i==j) continue;\r
-\r
- // check if both particles point to the same element (does not occur for mixed events, but if subsets are mixed within the same event)\r
- // if (tracksasso && trig->IsEqual(asso)) continue;\r
-\r
- if (tracksasso && (trig->GetUniqueID()==asso->GetUniqueID())) continue;\r
- \r
- if (fPtOrder)\r
- if (asso->Pt() >= trig->Pt()) continue;\r
-\r
- Int_t particlepidasso=asso->getparticle(); \r
- if(fAssociatedSpeciesSelection){ if (particlepidasso!=fAssociatedSpecies) continue;}\r
- \r
-\r
-if (fAssociatedSelectCharge != 0)\r
-if (asso->Charge() * fAssociatedSelectCharge < 0) continue;\r
- \r
- if (fSelectCharge > 0)\r
- {\r
- // skip like sign\r
- if (fSelectCharge == 1 && asso->Charge() * trig->Charge() > 0)\r
- continue;\r
- \r
- // skip unlike sign\r
- if (fSelectCharge == 2 && asso->Charge() * trig->Charge() < 0)\r
- continue;\r
- }\r
-\r
-if (fEtaOrdering)\r
- {\r
- if (trigeta < 0 && asso->Eta() < trigeta)\r
- continue;\r
- if (trigeta > 0 && asso->Eta() > trigeta)\r
- continue;\r
- }\r
-\r
-if (fRejectResonanceDaughters > 0)\r
- if (asso->TestBit(kResonanceDaughterFlag))\r
- {\r
-// Printf("Skipped j=%d", j);\r
- continue;\r
- }\r
-\r
- // conversions\r
- if (fCutConversions && asso->Charge() * trig->Charge() < 0)\r
- {\r
- Float_t mass = GetInvMassSquaredCheap(trig->Pt(), trigeta, trig->Phi(), asso->Pt(),eta[j], asso->Phi(), 0.510e-3, 0.510e-3);\r
- \r
- if (mass < 0.1)\r
- {\r
- mass = GetInvMassSquared(trig->Pt(), trigeta, trig->Phi(), asso->Pt(), eta[j], asso->Phi(), 0.510e-3, 0.510e-3);\r
- \r
- //fControlConvResoncances->Fill(0.0, mass);\r
-\r
- if (mass < 0.04*0.04) \r
- continue;\r
- }\r
- }\r
-\r
- // K0s\r
- if (fCutResonances && asso->Charge() * trig->Charge() < 0)\r
- {\r
- Float_t mass = GetInvMassSquaredCheap(trig->Pt(), trigeta, trig->Phi(), asso->Pt(), eta[j], asso->Phi(), 0.1396, 0.1396);\r
- \r
- const Float_t kK0smass = 0.4976;\r
- \r
- if (TMath::Abs(mass - kK0smass*kK0smass) < 0.1)\r
- {\r
- mass = GetInvMassSquared(trig->Pt(), trigeta, trig->Phi(), asso->Pt(),eta[j], asso->Phi(), 0.1396, 0.1396);\r
- \r
- //fControlConvResoncances->Fill(1, mass - kK0smass*kK0smass);\r
-\r
- if (mass > (kK0smass-0.02)*(kK0smass-0.02) && mass < (kK0smass+0.02)*(kK0smass+0.02))\r
- continue;\r
- }\r
- }\r
- \r
- // Lambda\r
- if (fCutResonances && asso->Charge() * trig->Charge() < 0)\r
- {\r
- Float_t mass1 = GetInvMassSquaredCheap(trig->Pt(), trigeta, trig->Phi(), asso->Pt(), eta[j], asso->Phi(), 0.1396, 0.9383);\r
- Float_t mass2 = GetInvMassSquaredCheap(trig->Pt(), trigeta, trig->Phi(), asso->Pt(),eta[j] , asso->Phi(), 0.9383, 0.1396);\r
- \r
- const Float_t kLambdaMass = 1.115;\r
-\r
- if (TMath::Abs(mass1 - kLambdaMass*kLambdaMass) < 0.1)\r
- {\r
- mass1 = GetInvMassSquared(trig->Pt(), trigeta, trig->Phi(), asso->Pt(),eta[j], asso->Phi(), 0.1396, 0.9383);\r
-\r
- //fControlConvResoncances->Fill(2, mass1 - kLambdaMass*kLambdaMass);\r
- \r
- if (mass1 > (kLambdaMass-0.02)*(kLambdaMass-0.02) && mass1 < (kLambdaMass+0.02)*(kLambdaMass+0.02))\r
- continue;\r
- }\r
- if (TMath::Abs(mass2 - kLambdaMass*kLambdaMass) < 0.1)\r
- {\r
- mass2 = GetInvMassSquared(trig->Pt(), trigeta, trig->Phi(), asso->Pt(),eta[j] , asso->Phi(), 0.9383, 0.1396);\r
-\r
- //fControlConvResoncances->Fill(2, mass2 - kLambdaMass*kLambdaMass);\r
-\r
- if (mass2 > (kLambdaMass-0.02)*(kLambdaMass-0.02) && mass2 < (kLambdaMass+0.02)*(kLambdaMass+0.02))\r
- continue;\r
- }\r
- }\r
-\r
- if (twoTrackEfficiencyCut)\r
- {\r
- // the variables & cuthave been developed by the HBT group \r
- // see e.g. https://indico.cern.ch/materialDisplay.py?contribId=36&sessionId=6&materialId=slides&confId=142700\r
- Float_t phi1 = trig->Phi();\r
- Float_t pt1 = trig->Pt();\r
- Float_t charge1 = trig->Charge();\r
- Float_t phi2 = asso->Phi();\r
- Float_t pt2 = asso->Pt();\r
- Float_t charge2 = asso->Charge();\r
-\r
- Float_t deta= trigeta - eta[j]; \r
- \r
- // optimization\r
- if (TMath::Abs(deta) < twoTrackEfficiencyCutValue * 2.5 * 3)\r
- {\r
-\r
- // check first boundaries to see if is worth to loop and find the minimum\r
- Float_t dphistar1 = GetDPhiStar(phi1, pt1, charge1, phi2, pt2, charge2, 0.8, bSign);\r
- Float_t dphistar2 = GetDPhiStar(phi1, pt1, charge1, phi2, pt2, charge2, 2.5, bSign);\r
-\r
- const Float_t kLimit = twoTrackEfficiencyCutValue * 3;\r
-\r
- Float_t dphistarminabs = 1e5;\r
- Float_t dphistarmin = 1e5;\r
-\r
- if (TMath::Abs(dphistar1) < kLimit || TMath::Abs(dphistar2) < kLimit || dphistar1 * dphistar2 < 0)\r
- {\r
- for (Double_t rad=0.8; rad<2.51; rad+=0.01) \r
- {\r
- Float_t dphistar = GetDPhiStar(phi1, pt1, charge1, phi2, pt2, charge2, rad, bSign);\r
-\r
- Float_t dphistarabs = TMath::Abs(dphistar);\r
-\r
- if (dphistarabs < dphistarminabs)\r
- {\r
- dphistarmin = dphistar;\r
- dphistarminabs = dphistarabs;\r
- }\r
- }\r
-\r
-if (dphistarminabs < twoTrackEfficiencyCutValue && TMath::Abs(deta) < twoTrackEfficiencyCutValue)\r
- {\r
-// Printf("Removed track pair %d %d with %f %f %f %f %f %f %f %f %f", i, j, deta, dphistarminabs, phi1, pt1, charge1, phi2, pt2, charge2, bSign);\r
- continue;\r
- }\r
-//fTwoTrackDistancePt[1]->Fill(deta, dphistarmin, TMath::Abs(pt1 - pt2));\r
-\r
- }\r
- }\r
- }\r
-\r
- Float_t trackeffasso=1.0;\r
- if(fapplyAssoefficiency) trackeffasso=asso->geteffcorrectionval();\r
- //cout<<"*******************trackeffasso="<<trackeffasso<<endl;\r
- Float_t deleta=trigeta-eta[j];\r
- Float_t delphi=PhiRange(trigphi-asso->Phi()); \r
-\r
- //here get the two particle efficiency correction factor\r
- Float_t effweight=trackefftrig*trackeffasso;\r
- //cout<<"*******************effweight="<<effweight<<endl;\r
- Double_t* vars;\r
- vars= new Double_t[kTrackVariablesPair];\r
- vars[0]=cent;\r
- vars[1]=vtx;\r
- vars[2]=trigpt;\r
- vars[3]=asso->Pt();\r
- vars[4]=deleta;\r
- vars[5]=delphi;\r
-if(fcontainPIDtrig && !fcontainPIDasso) vars[6]=particlepidtrig;\r
-if(!fcontainPIDtrig && fcontainPIDasso) vars[6]=particlepidasso;\r
- if(fcontainPIDtrig && fcontainPIDasso){\r
- vars[6]=particlepidtrig;\r
- vars[7]=particlepidasso;\r
- }\r
-\r
- //Fill Correlation ThnSparses\r
- if(fillup=="trigassoUNID")\r
- {\r
- if(mixcase==kFALSE) fTHnCorrUNID->Fill(vars,1.0/effweight);\r
- if(mixcase==kTRUE) fTHnCorrUNIDmix->Fill(vars,1.0/effweight);\r
- }\r
- if(fillup=="trigIDassoID")\r
- {\r
- if(mixcase==kFALSE) fTHnCorrID->Fill(vars,1.0/effweight);\r
- if(mixcase==kTRUE) fTHnCorrIDmix->Fill(vars,1.0/effweight);\r
- }\r
- if(fillup=="trigIDassoIDMCTRUTH")//******************************************for TRUTH MC particles\r
- {//in this case effweight should be 1 always \r
- if(mixcase==kFALSE) fCorrelatonTruthPrimary->Fill(vars,1.0/effweight); \r
- if(mixcase==kTRUE) fCorrelatonTruthPrimarymix->Fill(vars,1.0/effweight);\r
- } \r
- if(fillup=="trigIDassoUNID" || fillup=="trigUNIDassoID")//****************************be careful\r
- {\r
- if(mixcase==kFALSE) fTHnCorrIDUNID->Fill(vars,1.0/effweight);\r
- if(mixcase==kTRUE) fTHnCorrIDUNIDmix->Fill(vars,1.0/effweight);\r
- }\r
- \r
-delete[] vars;\r
- }//asso loop ends \r
-delete[] trigval;\r
- }//trigger loop ends \r
-\r
-}\r
-\r
-//--------------------------------------------------------------------------------\r
-Float_t AliTwoParticlePIDCorr::GetTrackbyTrackeffvalue(AliAODTrack* track,Double_t cent,Float_t evzvtx, Int_t parpid)\r
-{\r
- //This function is called only when applyefficiency=kTRUE; also ensure that "track" is present before calling that function\r
- Int_t effVars[4];\r
- Float_t effvalue=1.; \r
-\r
- if(parpid==unidentified)\r
- {\r
- effVars[0] = effcorection[5]->GetAxis(0)->FindBin(cent);\r
- effVars[1] = effcorection[5]->GetAxis(1)->FindBin(evzvtx); \r
- effVars[2] = effcorection[5]->GetAxis(2)->FindBin(track->Pt()); \r
- effVars[3] = effcorection[5]->GetAxis(3)->FindBin(track->Eta()); \r
- effvalue=effcorection[5]->GetBinContent(effVars);\r
- }\r
-if(parpid==SpPion || parpid==SpKaon)\r
- {\r
- if(fmesoneffrequired && !fkaonprotoneffrequired && track->Pt()>=fminPtComboeff && track->Pt()<=fmaxPtComboeff)\r
- {\r
- effVars[0] = effcorection[3]->GetAxis(0)->FindBin(cent);\r
- effVars[1] = effcorection[3]->GetAxis(1)->FindBin(evzvtx); \r
- effVars[2] = effcorection[3]->GetAxis(2)->FindBin(track->Pt()); \r
- effVars[3] = effcorection[3]->GetAxis(3)->FindBin(track->Eta());\r
- effvalue=effcorection[3]->GetBinContent(effVars);\r
- }\r
- else{\r
- if(parpid==SpPion)\r
- {\r
- effVars[0] = effcorection[0]->GetAxis(0)->FindBin(cent);\r
- effVars[1] = effcorection[0]->GetAxis(1)->FindBin(evzvtx); \r
- effVars[2] = effcorection[0]->GetAxis(2)->FindBin(track->Pt()); \r
- effVars[3] = effcorection[0]->GetAxis(3)->FindBin(track->Eta()); \r
- effvalue=effcorection[0]->GetBinContent(effVars);\r
- }\r
- \r
- if(parpid==SpKaon)\r
- {\r
- effVars[0] = effcorection[1]->GetAxis(0)->FindBin(cent);\r
- effVars[1] = effcorection[1]->GetAxis(1)->FindBin(evzvtx); \r
- effVars[2] = effcorection[1]->GetAxis(2)->FindBin(track->Pt()); \r
- effVars[3] = effcorection[1]->GetAxis(3)->FindBin(track->Eta()); \r
- effvalue=effcorection[1]->GetBinContent(effVars);\r
- }\r
- }\r
- } \r
- \r
- if(parpid==SpProton)\r
- {\r
- effVars[0] = effcorection[2]->GetAxis(0)->FindBin(cent);\r
- effVars[1] = effcorection[2]->GetAxis(1)->FindBin(evzvtx); \r
- effVars[2] = effcorection[2]->GetAxis(2)->FindBin(track->Pt()); \r
- effVars[3] = effcorection[2]->GetAxis(3)->FindBin(track->Eta()); \r
- effvalue=effcorection[2]->GetBinContent(effVars);\r
- }\r
-\r
- if(fkaonprotoneffrequired && !fmesoneffrequired && track->Pt()>=fminPtComboeff && track->Pt()<=fmaxPtComboeff)\r
- {\r
- if(parpid==SpProton || parpid==SpKaon)\r
- {\r
- effVars[0] = effcorection[4]->GetAxis(0)->FindBin(cent);\r
- effVars[1] = effcorection[4]->GetAxis(1)->FindBin(evzvtx); \r
- effVars[2] = effcorection[4]->GetAxis(2)->FindBin(track->Pt()); \r
- effVars[3] = effcorection[4]->GetAxis(3)->FindBin(track->Eta());\r
- effvalue=effcorection[4]->GetBinContent(effVars);\r
- }\r
- } \r
- // Printf("%d %d %d %d %f", effVars[0], effVars[1], effVars[2], effVars[3], fEfficiencyCorrectionAssociated->GetBinContent(effVars));\r
- if(effvalue==0.) effvalue=1.;\r
-\r
- return effvalue; \r
-\r
-}\r
-//-----------------------------------------------------------------------\r
-\r
-Int_t AliTwoParticlePIDCorr::ClassifyTrack(AliAODTrack* track,AliAODVertex* vertex,Float_t magfield)\r
-{ \r
- \r
- if(!track) return 0;\r
- Bool_t trackOK = track->TestFilterBit(fFilterBit);\r
- if(!trackOK) return 0;\r
- //select only primary traks(for data & reco MC tracks) \r
- if(fonlyprimarydatareco && track->GetType()!=AliAODTrack::kPrimary) return 0;\r
- if(track->Charge()==0) return 0;\r
- fHistQA[12]->Fill(track->GetTPCNcls()); \r
- Float_t dxy, dz; \r
- dxy = track->DCA();\r
- dz = track->ZAtDCA();\r
- fHistQA[6]->Fill(dxy);\r
- fHistQA[7]->Fill(dz);\r
- Float_t chi2ndf = track->Chi2perNDF();\r
- fHistQA[13]->Fill(chi2ndf); \r
- Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);\r
- fHistQA[14]->Fill(nCrossedRowsTPC); \r
- //Float_t ratioCrossedRowsOverFindableClustersTPC = 1.0;\r
- if (track->GetTPCNclsF()>0) {\r
- Float_t ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();\r
- fHistQA[15]->Fill(ratioCrossedRowsOverFindableClustersTPC);\r
- }\r
-//accepted tracks \r
- Float_t pt=track->Pt();\r
- if(pt< fminPt || pt> fmaxPt) return 0;\r
- if(TMath::Abs(track->Eta())> fmaxeta) return 0;\r
- if(track->Phi()<0. || track->Phi()>2*TMath::Pi()) return 0;\r
- //if (!HasTPCPID(track)) return 0;//trigger & associated particles must have TPC PID,Is it required\r
-// DCA XY\r
- if (fdcacut && fDCAXYCut)\r
- {\r
- if (!vertex)\r
- return 0;\r
- \r
- Double_t pos[2];\r
- Double_t covar[3];\r
- AliAODTrack* clone =(AliAODTrack*) track->Clone();\r
- Bool_t success = clone->PropagateToDCA(vertex, magfield, fdcacutvalue, pos, covar);\r
- delete clone;\r
- if (!success)\r
- return 0;\r
-\r
-// Printf("%f", ((AliAODTrack*)part)->DCA());\r
-// Printf("%f", pos[0]);\r
- if (TMath::Abs(pos[0]) > fDCAXYCut->Eval(track->Pt()))\r
- return 0;\r
- }\r
-\r
- fHistQA[8]->Fill(pt);\r
- fHistQA[9]->Fill(track->Eta());\r
- fHistQA[10]->Fill(track->Phi());\r
- return 1;\r
- }\r
- //________________________________________________________________________________\r
-void AliTwoParticlePIDCorr::CalculateNSigmas(AliAODTrack *track) \r
-{\r
-//This function is called within the func GetParticle() for accepted tracks only i.e.after call of Classifytrack() & for those tracks which have proper TPC PID response . combined nsigma(circular) cut only for particles having pt upto 4.0 Gev/c and beyond that use the asymmetric nsigma cut around pion's mean position in TPC ( while filling the TObjArray for trig & asso )\r
-Float_t pt=track->Pt();\r
-\r
-//it is assumed that every track that passed the filterbit have proper TPC response(!!)\r
-Float_t nsigmaTPCkPion =fPID->NumberOfSigmasTPC(track, AliPID::kPion);\r
-Float_t nsigmaTPCkKaon =fPID->NumberOfSigmasTPC(track, AliPID::kKaon);\r
-Float_t nsigmaTPCkProton =fPID->NumberOfSigmasTPC(track, AliPID::kProton);\r
-\r
-Float_t nsigmaTOFkProton=999.,nsigmaTOFkKaon=999.,nsigmaTOFkPion=999.;\r
-Float_t nsigmaTPCTOFkProton=999.,nsigmaTPCTOFkKaon=999.,nsigmaTPCTOFkPion=999.;\r
-\r
- if(HasTOFPID(track) && pt>fPtTOFPIDmin)\r
- {\r
-\r
-nsigmaTOFkPion =fPID->NumberOfSigmasTOF(track, AliPID::kPion);\r
-nsigmaTOFkKaon =fPID->NumberOfSigmasTOF(track, AliPID::kKaon);\r
-nsigmaTOFkProton =fPID->NumberOfSigmasTOF(track, AliPID::kProton);\r
-//---combined\r
-nsigmaTPCTOFkPion = TMath::Sqrt(nsigmaTPCkPion*nsigmaTPCkPion+nsigmaTOFkPion*nsigmaTOFkPion);\r
-nsigmaTPCTOFkKaon = TMath::Sqrt(nsigmaTPCkKaon*nsigmaTPCkKaon+nsigmaTOFkKaon*nsigmaTOFkKaon);\r
-nsigmaTPCTOFkProton = TMath::Sqrt(nsigmaTPCkProton*nsigmaTPCkProton+nsigmaTOFkProton*nsigmaTOFkProton);\r
-\r
-\r
- }\r
-else{\r
- // --- combined\r
- // if TOF is missing and below fPtTOFPID only the TPC information is used\r
- nsigmaTPCTOFkProton = TMath::Abs(nsigmaTPCkProton);\r
- nsigmaTPCTOFkKaon = TMath::Abs(nsigmaTPCkKaon);\r
- nsigmaTPCTOFkPion = TMath::Abs(nsigmaTPCkPion);\r
-\r
- }\r
-\r
-//set data member fnsigmas\r
- fnsigmas[SpPion][NSigmaTPC]=nsigmaTPCkPion;\r
- fnsigmas[SpKaon][NSigmaTPC]=nsigmaTPCkKaon;\r
- fnsigmas[SpProton][NSigmaTPC]=nsigmaTPCkProton;\r
-\r
- //for all tracks below fPtTOFPIDmin and also for tracks above fPtTOFPIDmin without proper TOF response these TOF nsigma values will be 999.\r
- fnsigmas[SpPion][NSigmaTOF]=nsigmaTOFkPion;\r
- fnsigmas[SpKaon][NSigmaTOF]=nsigmaTOFkKaon;\r
- fnsigmas[SpProton][NSigmaTOF]=nsigmaTOFkProton;\r
-\r
- //for all tracks below fPtTOFPIDmin and also for tracks above fPtTOFPIDmin without proper TOF response these TPCTOF nsigma values will be TMath::Abs(TPC only nsigma)\r
- fnsigmas[SpPion][NSigmaTPCTOF]=nsigmaTPCTOFkPion;\r
- fnsigmas[SpKaon][NSigmaTPCTOF]=nsigmaTPCTOFkKaon;\r
- fnsigmas[SpProton][NSigmaTPCTOF]=nsigmaTPCTOFkProton;\r
-\r
-\r
-}\r
-//----------------------------------------------------------------------------\r
-Int_t AliTwoParticlePIDCorr::FindMinNSigma(AliAODTrack *track) \r
-{\r
- //this function is always called after calling the function CalculateNSigmas(AliAODTrack *track)\r
-if(fRequestTOFPID && track->Pt()>fPtTOFPIDmin && track->Pt()<=fPtTOFPIDmax && (!HasTOFPID(track)) )return SpUndefined;//so any track having Pt>0.6 && Pt<=4.0 Gev withot having proper TOF response will be defined as SpUndefined\r
-//get the identity of the particle with the minimum Nsigma\r
- Float_t nsigmaPion=999., nsigmaKaon=999., nsigmaProton=999.;\r
- switch (fPIDType){\r
- case NSigmaTPC:\r
- nsigmaProton = TMath::Abs(fnsigmas[SpProton][NSigmaTPC]);\r
- nsigmaKaon = TMath::Abs(fnsigmas[SpKaon][NSigmaTPC]) ;\r
- nsigmaPion = TMath::Abs(fnsigmas[SpPion][NSigmaTPC]) ;\r
- break;\r
- case NSigmaTOF:\r
- nsigmaProton = TMath::Abs(fnsigmas[SpProton][NSigmaTOF]);\r
- nsigmaKaon = TMath::Abs(fnsigmas[SpKaon][NSigmaTOF]) ;\r
- nsigmaPion = TMath::Abs(fnsigmas[SpPion][NSigmaTOF]) ;\r
- break;\r
- case NSigmaTPCTOF://In case of no TOF matching the combined nsigma is the TPC one\r
- nsigmaProton = TMath::Abs(fnsigmas[SpProton][NSigmaTPCTOF]);\r
- nsigmaKaon = TMath::Abs(fnsigmas[SpKaon][NSigmaTPCTOF]) ;\r
- nsigmaPion = TMath::Abs(fnsigmas[SpPion][NSigmaTPCTOF]) ;\r
- break;\r
- }\r
-\r
- if(track->Pt()<=fPtTOFPIDmax) {\r
- // guess the particle based on the smaller nsigma (within fNSigmaPID)\r
- if( ( nsigmaKaon==nsigmaPion ) && ( nsigmaKaon==nsigmaProton )) return SpUndefined;//it is the default value for the three\r
-if( ( nsigmaKaon < nsigmaPion ) && ( nsigmaKaon < nsigmaProton ) && (nsigmaKaon < fNSigmaPID)) return SpKaon;\r
-if( ( nsigmaPion < nsigmaKaon ) && ( nsigmaPion < nsigmaProton ) && (nsigmaPion < fNSigmaPID)) return SpPion;\r
-if( ( nsigmaProton < nsigmaKaon ) && ( nsigmaProton < nsigmaPion ) && (nsigmaProton < fNSigmaPID)) return SpProton;\r
-\r
-// else, return undefined\r
- return SpUndefined;\r
- }\r
- else {//asymmetric nsigma cut around pion's mean position for tracks having Pt>4.0 Gev\r
- if(fminprotonsigmacut<fnsigmas[SpPion][NSigmaTPC] && fnsigmas[SpPion][NSigmaTPC]<fmaxprotonsigmacut) return SpProton;\r
- if(fminpionsigmacut<fnsigmas[SpPion][NSigmaTPC] && fnsigmas[SpPion][NSigmaTPC]<fmaxpionsigmacut) return SpPion;\r
-// else, return undefined(here we don't detect kaons)\r
- return SpUndefined;\r
- }\r
-\r
-}\r
-\r
-//------------------------------------------------------------------------------------------\r
-Bool_t* AliTwoParticlePIDCorr::GetDoubleCounting(AliAODTrack * trk){ \r
- //this function is always called after calling the function CalculateNSigmas(AliAODTrack *track)\r
-\r
- //if a particle has double counting set fHasDoubleCounting[ipart]=kTRUE\r
- //fill DC histos\r
- for(Int_t ipart=0;ipart<NSpecies;ipart++)fHasDoubleCounting[ipart]=kFALSE;//array with kTRUE for second (or third) identity of the track\r
- \r
- Int_t MinNSigma=FindMinNSigma(trk);//not filling the NSigmaRec histos\r
- \r
- \r
- if(MinNSigma==SpUndefined)return fHasDoubleCounting;//in case of undefined no Double counting\r
- \r
- Float_t nsigmaPion=999., nsigmaKaon=999., nsigmaProton=999.;\r
- switch (fPIDType) {\r
- case NSigmaTPC:\r
- nsigmaProton = TMath::Abs(fnsigmas[SpProton][NSigmaTPC]);\r
- nsigmaKaon = TMath::Abs(fnsigmas[SpKaon][NSigmaTPC]) ;\r
- nsigmaPion = TMath::Abs(fnsigmas[SpPion][NSigmaTPC]) ;\r
- break;\r
- case NSigmaTOF:\r
- nsigmaProton = TMath::Abs(fnsigmas[SpProton][NSigmaTOF]);\r
- nsigmaKaon = TMath::Abs(fnsigmas[SpKaon][NSigmaTOF]) ;\r
- nsigmaPion = TMath::Abs(fnsigmas[SpPion][NSigmaTOF]) ;\r
- break;\r
- case NSigmaTPCTOF://In case of no TOF matching the combined nsigma is the TPC one\r
- nsigmaProton = TMath::Abs(fnsigmas[SpProton][NSigmaTPCTOF]);\r
- nsigmaKaon = TMath::Abs(fnsigmas[SpKaon][NSigmaTPCTOF]) ;\r
- nsigmaPion = TMath::Abs(fnsigmas[SpPion][NSigmaTPCTOF]) ;\r
- break;\r
- }\r
-\r
- //there is chance of overlapping only for particles having pt below 4.0 GEv\r
- if(trk->Pt()<=4.0){\r
- if(nsigmaPion<fNSigmaPID && MinNSigma!=SpPion)fHasDoubleCounting[SpPion]=kTRUE;\r
- if(nsigmaKaon<fNSigmaPID && MinNSigma!=SpKaon)fHasDoubleCounting[SpKaon]=kTRUE;\r
- if(nsigmaProton<fNSigmaPID && MinNSigma!=SpProton)fHasDoubleCounting[SpProton]=kTRUE;\r
- \r
- }\r
-\r
- \r
- return fHasDoubleCounting;\r
-}\r
-\r
-//////////////////////////////////////////////////////////////////////////////////////////////////\r
-Int_t AliTwoParticlePIDCorr::GetParticle(AliAODTrack * trk){ \r
- //return the specie according to the minimum nsigma value\r
- //no double counting, this has to be evaluated using CheckDoubleCounting()\r
- //Printf("fPtTOFPID %.1f, fRequestTOFPID %d, fNSigmaPID %.1f, fPIDType %d",fPtTOFPID,fRequestTOFPID,fNSigmaPID,fPIDType);\r
- \r
- CalculateNSigmas(trk);//fill the data member fnsigmas with the nsigmas value [ipart][iPID]\r
- \r
- if(fUseExclusiveNSigma){\r
- Bool_t *HasDC;\r
- HasDC=GetDoubleCounting(trk);\r
- for(Int_t ipart=0;ipart<NSpecies;ipart++){\r
- if(HasDC[ipart]==kTRUE) return SpUndefined;\r
- }\r
- return FindMinNSigma(trk);//NSigmaRec distr filled here\r
- }\r
- else return FindMinNSigma(trk);//NSigmaRec distr filled here\r
- \r
-}\r
-\r
-//-------------------------------------------------------------------------------------\r
-Bool_t\r
-AliTwoParticlePIDCorr::HasTPCPID(AliAODTrack *track) const\r
-{ \r
- // check PID signal \r
- AliPIDResponse::EDetPidStatus statustpc = fPID->CheckPIDStatus(AliPIDResponse::kTPC,track);\r
- if(statustpc!=AliPIDResponse::kDetPidOk) return kFALSE;\r
- //ULong_t status=track->GetStatus();\r
- //if (!( (status & AliAODTrack::kTPCpid ) == AliAODTrack::kTPCpid )) return kFALSE;//remove light nuclei\r
- //if (track->GetTPCsignal() <= 0.) return kFALSE;\r
- // if(track->GetTPCsignalN() < 60) return kFALSE;//tracks with TPCsignalN< 60 have questionable dEdx,cutting on TPCsignalN > 70 or > 60 shouldn't make too much difference in statistics,also it is IMO safe to use TPC also for MIPs.\r
- \r
- return kTRUE; \r
-}\r
-//___________________________________________________________\r
-\r
-Bool_t\r
-AliTwoParticlePIDCorr::HasTOFPID(AliAODTrack *track) const\r
-{\r
- // check TOF matched track \r
- //ULong_t status=track->GetStatus();\r
- //if (!( (status & AliAODTrack::kITSin ) == AliAODTrack::kITSin )) return kFALSE;\r
- AliPIDResponse::EDetPidStatus statustof = fPID->CheckPIDStatus(AliPIDResponse::kTOF,track);\r
- if(statustof!= AliPIDResponse::kDetPidOk) return kFALSE;\r
- if(track->Pt()<=fPtTOFPIDmin) return kFALSE;\r
- //if(!((status & AliAODTrack::kTOFpid ) == AliAODTrack::kTOFpid )) return kFALSE;\r
- //Float_t probMis = fPIDresponse->GetTOFMismatchProbability(track);\r
- // if (probMis > 0.01) return kFALSE;\r
-if(fRemoveTracksT0Fill)\r
- {\r
-Int_t startTimeMask = fPID->GetTOFResponse().GetStartTimeMask(track->P());\r
- if (startTimeMask < 0)return kFALSE; \r
- }\r
- return kTRUE;\r
-}\r
-\r
-//________________________________________________________________________\r
-Float_t AliTwoParticlePIDCorr :: GetBeta(AliAODTrack *track)\r
-{\r
- //it is called only when TOF PID is available\r
- Double_t p = track->P();\r
- Double_t time=track->GetTOFsignal()-fPID->GetTOFResponse().GetStartTime(p);\r
- Double_t timei[5];\r
- track->GetIntegratedTimes(timei);\r
- return timei[0]/time;\r
-}\r
-//------------------------------------------------------------------------------------------------------\r
-\r
-Float_t AliTwoParticlePIDCorr::GetInvMassSquared(Float_t pt1, Float_t eta1, Float_t phi1, Float_t pt2, Float_t eta2, Float_t phi2, Float_t m0_1, Float_t m0_2)\r
-{\r
- // calculate inv mass squared\r
- // same can be achieved, but with more computing time with\r
- /*TLorentzVector photon, p1, p2;\r
- p1.SetPtEtaPhiM(triggerParticle->Pt(), triggerEta, triggerParticle->Phi(), 0.510e-3);\r
- p2.SetPtEtaPhiM(particle->Pt(), eta[j], particle->Phi(), 0.510e-3);\r
- photon = p1+p2;\r
- photon.M()*/\r
- \r
- Float_t tantheta1 = 1e10;\r
- \r
- if (eta1 < -1e-10 || eta1 > 1e-10)\r
- tantheta1 = 2 * TMath::Exp(-eta1) / ( 1 - TMath::Exp(-2*eta1));\r
- \r
- Float_t tantheta2 = 1e10;\r
- if (eta2 < -1e-10 || eta2 > 1e-10)\r
- tantheta2 = 2 * TMath::Exp(-eta2) / ( 1 - TMath::Exp(-2*eta2));\r
- \r
- Float_t e1squ = m0_1 * m0_1 + pt1 * pt1 * (1.0 + 1.0 / tantheta1 / tantheta1);\r
- Float_t e2squ = m0_2 * m0_2 + pt2 * pt2 * (1.0 + 1.0 / tantheta2 / tantheta2);\r
- \r
- Float_t mass2 = m0_1 * m0_1 + m0_2 * m0_2 + 2 * ( TMath::Sqrt(e1squ * e2squ) - ( pt1 * pt2 * ( TMath::Cos(phi1 - phi2) + 1.0 / tantheta1 / tantheta2 ) ) );\r
- \r
- return mass2;\r
-}\r
-//---------------------------------------------------------------------------------\r
-\r
-Float_t AliTwoParticlePIDCorr::GetInvMassSquaredCheap(Float_t pt1, Float_t eta1, Float_t phi1, Float_t pt2, Float_t eta2, Float_t phi2, Float_t m0_1, Float_t m0_2)\r
-{\r
- // calculate inv mass squared approximately\r
- \r
- Float_t tantheta1 = 1e10;\r
- \r
- if (eta1 < -1e-10 || eta1 > 1e-10)\r
- {\r
- Float_t expTmp = 1.0-eta1+eta1*eta1/2-eta1*eta1*eta1/6+eta1*eta1*eta1*eta1/24;\r
- tantheta1 = 2.0 * expTmp / ( 1.0 - expTmp*expTmp);\r
- }\r
- \r
- Float_t tantheta2 = 1e10;\r
- if (eta2 < -1e-10 || eta2 > 1e-10)\r
- {\r
- Float_t expTmp = 1.0-eta2+eta2*eta2/2-eta2*eta2*eta2/6+eta2*eta2*eta2*eta2/24;\r
- tantheta2 = 2.0 * expTmp / ( 1.0 - expTmp*expTmp);\r
- }\r
- \r
- Float_t e1squ = m0_1 * m0_1 + pt1 * pt1 * (1.0 + 1.0 / tantheta1 / tantheta1);\r
- Float_t e2squ = m0_2 * m0_2 + pt2 * pt2 * (1.0 + 1.0 / tantheta2 / tantheta2);\r
- \r
- // fold onto 0...pi\r
- Float_t deltaPhi = TMath::Abs(phi1 - phi2);\r
- while (deltaPhi > TMath::TwoPi())\r
- deltaPhi -= TMath::TwoPi();\r
- if (deltaPhi > TMath::Pi())\r
- deltaPhi = TMath::TwoPi() - deltaPhi;\r
- \r
- Float_t cosDeltaPhi = 0;\r
- if (deltaPhi < TMath::Pi()/3)\r
- cosDeltaPhi = 1.0 - deltaPhi*deltaPhi/2 + deltaPhi*deltaPhi*deltaPhi*deltaPhi/24;\r
- else if (deltaPhi < 2*TMath::Pi()/3)\r
- cosDeltaPhi = -(deltaPhi - TMath::Pi()/2) + 1.0/6 * TMath::Power((deltaPhi - TMath::Pi()/2), 3);\r
- else\r
- cosDeltaPhi = -1.0 + 1.0/2.0*(deltaPhi - TMath::Pi())*(deltaPhi - TMath::Pi()) - 1.0/24.0 * TMath::Power(deltaPhi - TMath::Pi(), 4);\r
- \r
- Float_t mass2 = m0_1 * m0_1 + m0_2 * m0_2 + 2 * ( TMath::Sqrt(e1squ * e2squ) - ( pt1 * pt2 * ( cosDeltaPhi + 1.0 / tantheta1 / tantheta2 ) ) );\r
- \r
-// Printf(Form("%f %f %f %f %f %f %f %f %f", pt1, eta1, phi1, pt2, eta2, phi2, m0_1, m0_2, mass2));\r
- \r
- return mass2;\r
-}\r
-//--------------------------------------------------------------------------------\r
-Float_t AliTwoParticlePIDCorr::GetDPhiStar(Float_t phi1, Float_t pt1, Float_t charge1, Float_t phi2, Float_t pt2, Float_t charge2, Float_t radius, Float_t bSign)\r
-{ \r
- //\r
- // calculates dphistar\r
- //\r
- \r
- Float_t dphistar = phi1 - phi2 - charge1 * bSign * TMath::ASin(0.075 * radius / pt1) + charge2 * bSign * TMath::ASin(0.075 * radius / pt2);\r
- \r
- static const Double_t kPi = TMath::Pi();\r
- \r
- // circularity\r
-// if (dphistar > 2 * kPi)\r
-// dphistar -= 2 * kPi;\r
-// if (dphistar < -2 * kPi)\r
-// dphistar += 2 * kPi;\r
- \r
- if (dphistar > kPi)\r
- dphistar = kPi * 2 - dphistar;\r
- if (dphistar < -kPi)\r
- dphistar = -kPi * 2 - dphistar;\r
- if (dphistar > kPi) // might look funny but is needed\r
- dphistar = kPi * 2 - dphistar;\r
- \r
- return dphistar;\r
-}\r
-//_________________________________________________________________________\r
-void AliTwoParticlePIDCorr ::DefineEventPool()\r
-{\r
-const Int_t MaxNofEvents=1000;\r
-const Int_t MaxNofTracks=50000;\r
-const Int_t NofVrtxBins=10+(1+10)*2;\r
-Double_t ZvrtxBins[NofVrtxBins+1]={ -10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10, \r
- 90, 92, 94, 96, 98, 100, 102, 104, 106, 108, 110, \r
- 190, 192, 194, 196, 198, 200, 202, 204, 206, 208, 210 \r
- };\r
- if (fSampleType=="pp"){\r
-const Int_t NofCentBins=5;\r
-Double_t CentralityBins[NofCentBins+1]={0.,10., 20., 40.,80.,200.1};\r
-fPoolMgr = new AliEventPoolManager(MaxNofEvents,MaxNofTracks,NofCentBins,CentralityBins,NofVrtxBins,ZvrtxBins);\r
- }\r
-if (fSampleType=="pPb" || fSampleType=="PbPb")\r
- {\r
-const Int_t NofCentBins=15;\r
-Double_t CentralityBins[NofCentBins+1]={0., 1., 2., 3., 4., 5., 10., 20., 30., 40., 50., 60., 70., 80., 90., 100.1 };\r
-fPoolMgr = new AliEventPoolManager(MaxNofEvents,MaxNofTracks,NofCentBins,CentralityBins,NofVrtxBins,ZvrtxBins);\r
- }\r
-fPoolMgr->SetTargetValues(MaxNofTracks, 0.1, 5);\r
-\r
-//if(!fPoolMgr) return kFALSE;\r
-//return kTRUE;\r
-\r
-}\r
-//------------------------------------------------------------------------\r
-Double_t* AliTwoParticlePIDCorr::GetBinning(const char* configuration, const char* tag, Int_t& nBins)\r
-{\r
- // This method is a copy from AliUEHist::GetBinning\r
- // takes the binning from <configuration> identified by <tag>\r
- // configuration syntax example:\r
- // eta: 2.4, -2.3, -2.2, -2.1, -2.0, -1.9, -1.8, -1.7, -1.6, -1.5, -1.4, -1.3, -1.2, -1.1, -1.0, -0.9, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, 2.1, 2.2, 2.3, 2.4\r
- // phi: .....\r
- //\r
- // returns bin edges which have to be deleted by the caller\r
- \r
- TString config(configuration);\r
- TObjArray* lines = config.Tokenize("\n");\r
- for (Int_t i=0; i<lines->GetEntriesFast(); i++)\r
- {\r
- TString line(lines->At(i)->GetName());\r
- if (line.BeginsWith(TString(tag) + ":"))\r
- {\r
- line.Remove(0, strlen(tag) + 1);\r
- line.ReplaceAll(" ", "");\r
- TObjArray* binning = line.Tokenize(",");\r
- Double_t* bins = new Double_t[binning->GetEntriesFast()];\r
- for (Int_t j=0; j<binning->GetEntriesFast(); j++)\r
- bins[j] = TString(binning->At(j)->GetName()).Atof();\r
- \r
- nBins = binning->GetEntriesFast() - 1;\r
-\r
- delete binning;\r
- delete lines;\r
- return bins;\r
- }\r
- }\r
- \r
- delete lines;\r
- AliFatal(Form("Tag %s not found in %s", tag, configuration));\r
- return 0;\r
-}\r
-//_______________________________________________________________________________\r
-void AliTwoParticlePIDCorr::SetAsymmetricBin(THnSparse *h,Int_t dim,Double_t *arraybin,Int_t arraybinsize,TString axisTitle) \r
-{\r
- TAxis *axis = 0x0;\r
- axis = h->GetAxis(dim);\r
- axis->Set(arraybinsize,arraybin);\r
- axis->SetTitle(axisTitle);\r
-}\r
-\r
-//________________________________________________________________________\r
-void AliTwoParticlePIDCorr::Terminate(Option_t *) \r
-{\r
- // Draw result to screen, or perform fitting, normalizations\r
- // Called once at the end of the query\r
- fOutput = dynamic_cast<TList*> (GetOutputData(1));\r
- if(!fOutput) { Printf("ERROR: could not retrieve TList fOutput"); return; }\r
- \r
- \r
-}\r
-//------------------------------------------------------------------ \r
- \r
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: The ALICE Off-line Project. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+
+#include "AliTwoParticlePIDCorr.h"
+#include "AliVParticle.h"
+#include "TFormula.h"
+#include "TAxis.h"
+#include "TChain.h"
+#include "TTree.h"
+#include "TH1F.h"
+#include "TH2F.h"
+#include "TH3F.h"
+#include "TList.h"
+#include "TFile.h"
+
+#include "AliCentrality.h"
+#include "Riostream.h"
+
+#include <TSpline.h>
+#include <AliPID.h>
+#include "AliESDpid.h"
+#include "AliAODpidUtil.h"
+#include <AliPIDResponse.h>
+
+#include <AliAnalysisManager.h>
+#include <AliInputEventHandler.h>
+#include "AliAODInputHandler.h"
+
+#include "AliAnalysisTaskSE.h"
+#include "AliAnalysisManager.h"
+#include "AliCentrality.h"
+
+#include "AliVEvent.h"
+#include "AliAODEvent.h"
+#include "AliAODTrack.h"
+#include "AliVTrack.h"
+
+#include "THnSparse.h"
+
+#include "AliAODMCHeader.h"
+#include "AliAODMCParticle.h"
+#include "AliMCEventHandler.h"
+#include "AliMCEvent.h"
+#include "AliMCParticle.h"
+#include "TParticle.h"
+#include <TDatabasePDG.h>
+#include <TParticlePDG.h>
+
+#include "AliGenCocktailEventHeader.h"
+#include "AliGenEventHeader.h"
+
+#include "AliEventPoolManager.h"
+//#include "AliAnalysisUtils.h"
+using namespace AliPIDNameSpace;
+using namespace std;
+
+ClassImp(AliTwoParticlePIDCorr)
+ClassImp(LRCParticlePID)
+//________________________________________________________________________
+AliTwoParticlePIDCorr::AliTwoParticlePIDCorr() // All data members should be initialised here
+:AliAnalysisTaskSE(),
+ fOutput(0),
+ fCentralityMethod("V0A"),
+ fSampleType("pPb"),
+fnTracksVertex(1), // QA tracks pointing to principal vertex (= 3 default)
+ trkVtx(0),
+ zvtx(0),
+ fFilterBit(768),
+ fzvtxcut(10.0),
+ ffilltrigassoUNID(kFALSE),
+ ffilltrigUNIDassoID(kFALSE),
+ ffilltrigIDassoUNID(kTRUE),
+ ffilltrigIDassoID(kFALSE),
+ ffilltrigIDassoIDMCTRUTH(kFALSE),
+ fPtOrderMCTruth(kFALSE),
+ fTriggerSpeciesSelection(kFALSE),
+ fAssociatedSpeciesSelection(kFALSE),
+ fTriggerSpecies(SpPion),
+ fAssociatedSpecies(SpPion),
+ fCustomBinning(""),
+ fBinningString(""),
+ fcontainPIDtrig(kTRUE),
+ fcontainPIDasso(kFALSE),
+//frejectPileUp(kFALSE),
+ fminPt(0.2),
+ fmaxPt(10.0),
+ fmineta(-0.8),
+ fmaxeta(0.8),
+ fminprotonsigmacut(-6.0),
+ fmaxprotonsigmacut(-3.0),
+ fminpionsigmacut(0.0),
+ fmaxpionsigmacut(4.0),
+ fselectprimaryTruth(kTRUE),
+ fonlyprimarydatareco(kFALSE),
+ fdcacut(kFALSE),
+ fdcacutvalue(3.0),
+ ffillhistQAReco(kFALSE),
+ ffillhistQATruth(kFALSE),
+ kTrackVariablesPair(0),
+ fminPtTrig(0),
+ fmaxPtTrig(0),
+ fminPtComboeff(2.0),
+ fmaxPtComboeff(4.0),
+ fminPtAsso(0),
+ fmaxPtAsso(0),
+ fhistcentrality(0),
+ fEventCounter(0),
+ fEtaSpectrasso(0),
+ fphiSpectraasso(0),
+ MCtruthpt(0),
+ MCtrutheta(0),
+ MCtruthphi(0),
+ MCtruthpionpt(0),
+ MCtruthpioneta(0),
+ MCtruthpionphi(0),
+ MCtruthkaonpt(0),
+ MCtruthkaoneta(0),
+ MCtruthkaonphi(0),
+ MCtruthprotonpt(0),
+ MCtruthprotoneta(0),
+ MCtruthprotonphi(0),
+ fPioncont(0),
+ fKaoncont(0),
+ fProtoncont(0),
+ fCentralityCorrelation(0x0),
+ fHistoTPCdEdx(0x0),
+ fHistoTOFbeta(0x0),
+ fTPCTOFPion3d(0),
+ fTPCTOFKaon3d(0),
+ fTPCTOFProton3d(0),
+ fPionPt(0),
+ fPionEta(0),
+ fPionPhi(0),
+ fKaonPt(0),
+ fKaonEta(0),
+ fKaonPhi(0),
+ fProtonPt(0),
+ fProtonEta(0),
+ fProtonPhi(0),
+ fCorrelatonTruthPrimary(0),
+ fCorrelatonTruthPrimarymix(0),
+ fTHnCorrUNID(0),
+ fTHnCorrUNIDmix(0),
+ fTHnCorrID(0),
+ fTHnCorrIDmix(0),
+ fTHnCorrIDUNID(0),
+ fTHnCorrIDUNIDmix(0),
+ fTHnTrigcount(0),
+ fTHnTrigcountMCTruthPrim(0),
+ fPoolMgr(0x0),
+ fArrayMC(0),
+ fAnalysisType("MCAOD"),
+ fefffilename(""),
+ twoTrackEfficiencyCutValue(0.02),
+//fControlConvResoncances(0),
+ fPID(NULL),
+ eventno(0),
+ fPtTOFPIDmin(0.6),
+ fPtTOFPIDmax(4.0),
+ fRequestTOFPID(kTRUE),
+ fPIDType(NSigmaTPCTOF),
+ fNSigmaPID(3),
+ fUseExclusiveNSigma(kFALSE),
+ fRemoveTracksT0Fill(kFALSE),
+fSelectCharge(0),
+fTriggerSelectCharge(0),
+fAssociatedSelectCharge(0),
+fTriggerRestrictEta(-1),
+fEtaOrdering(kFALSE),
+fCutConversions(kFALSE),
+fCutResonances(kFALSE),
+fRejectResonanceDaughters(-1),
+ fOnlyOneEtaSide(0),
+fInjectedSignals(kFALSE),
+ fRemoveWeakDecays(kFALSE),
+fRemoveDuplicates(kFALSE),
+ fapplyTrigefficiency(kFALSE),
+ fapplyAssoefficiency(kFALSE),
+ ffillefficiency(kFALSE),
+ fmesoneffrequired(kFALSE),
+ fkaonprotoneffrequired(kFALSE),
+//fAnalysisUtils(0x0),
+ fDCAXYCut(0)
+
+{
+ for ( Int_t i = 0; i < 16; i++) {
+ fHistQA[i] = NULL;
+ }
+
+ for ( Int_t i = 0; i < 6; i++ ){
+ fTHnrecomatchedallPid[i] = NULL;
+ fTHngenprimPidTruth[i] = NULL;
+ effcorection[i]=NULL;
+ //effmap[i]=NULL;
+
+ }
+
+ for(Int_t ipart=0;ipart<NSpecies;ipart++)
+ for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++)
+ fnsigmas[ipart][ipid]=999.;
+
+ for(Int_t ipart=0;ipart<NSpecies;ipart++) {fHasDoubleCounting[ipart]=kFALSE;}
+
+ }
+//________________________________________________________________________
+AliTwoParticlePIDCorr::AliTwoParticlePIDCorr(const char *name) // All data members should be initialised here
+ :AliAnalysisTaskSE(name),
+ fOutput(0),
+ fCentralityMethod("V0A"),
+ fSampleType("pPb"),
+ fnTracksVertex(1), // QA tracks pointing to principal vertex (= 3 default)
+ trkVtx(0),
+ zvtx(0),
+ fFilterBit(768),
+ fzvtxcut(10.0),
+ ffilltrigassoUNID(kFALSE),
+ ffilltrigUNIDassoID(kFALSE),
+ ffilltrigIDassoUNID(kTRUE),
+ ffilltrigIDassoID(kFALSE),
+ ffilltrigIDassoIDMCTRUTH(kFALSE),
+ fPtOrderMCTruth(kFALSE),
+ fTriggerSpeciesSelection(kFALSE),
+ fAssociatedSpeciesSelection(kFALSE),
+ fTriggerSpecies(SpPion),
+ fAssociatedSpecies(SpPion),
+ fCustomBinning(""),
+ fBinningString(""),
+ fcontainPIDtrig(kTRUE),
+ fcontainPIDasso(kFALSE),
+ // frejectPileUp(kFALSE),
+ fminPt(0.2),
+ fmaxPt(10.0),
+ fmineta(-0.8),
+ fmaxeta(0.8),
+ fminprotonsigmacut(-6.0),
+ fmaxprotonsigmacut(-3.0),
+ fminpionsigmacut(0.0),
+ fmaxpionsigmacut(4.0),
+ fselectprimaryTruth(kTRUE),
+ fonlyprimarydatareco(kFALSE),
+ fdcacut(kFALSE),
+ fdcacutvalue(3.0),
+ ffillhistQAReco(kFALSE),
+ ffillhistQATruth(kFALSE),
+ kTrackVariablesPair(0),
+ fminPtTrig(0),
+ fmaxPtTrig(0),
+ fminPtComboeff(2.0),
+ fmaxPtComboeff(4.0),
+ fminPtAsso(0),
+ fmaxPtAsso(0),
+ fhistcentrality(0),
+ fEventCounter(0),
+ fEtaSpectrasso(0),
+ fphiSpectraasso(0),
+ MCtruthpt(0),
+ MCtrutheta(0),
+ MCtruthphi(0),
+ MCtruthpionpt(0),
+ MCtruthpioneta(0),
+ MCtruthpionphi(0),
+ MCtruthkaonpt(0),
+ MCtruthkaoneta(0),
+ MCtruthkaonphi(0),
+ MCtruthprotonpt(0),
+ MCtruthprotoneta(0),
+ MCtruthprotonphi(0),
+ fPioncont(0),
+ fKaoncont(0),
+ fProtoncont(0),
+ fCentralityCorrelation(0x0),
+ fHistoTPCdEdx(0x0),
+ fHistoTOFbeta(0x0),
+ fTPCTOFPion3d(0),
+ fTPCTOFKaon3d(0),
+ fTPCTOFProton3d(0),
+ fPionPt(0),
+ fPionEta(0),
+ fPionPhi(0),
+ fKaonPt(0),
+ fKaonEta(0),
+ fKaonPhi(0),
+ fProtonPt(0),
+ fProtonEta(0),
+ fProtonPhi(0),
+ fCorrelatonTruthPrimary(0),
+ fCorrelatonTruthPrimarymix(0),
+ fTHnCorrUNID(0),
+ fTHnCorrUNIDmix(0),
+ fTHnCorrID(0),
+ fTHnCorrIDmix(0),
+ fTHnCorrIDUNID(0),
+ fTHnCorrIDUNIDmix(0),
+ fTHnTrigcount(0),
+ fTHnTrigcountMCTruthPrim(0),
+ fPoolMgr(0x0),
+ fArrayMC(0),
+ fAnalysisType("MCAOD"),
+ fefffilename(""),
+ twoTrackEfficiencyCutValue(0.02),
+//fControlConvResoncances(0),
+ fPID(NULL),
+ eventno(0),
+ fPtTOFPIDmin(0.6),
+ fPtTOFPIDmax(4.0),
+ fRequestTOFPID(kTRUE),
+ fPIDType(NSigmaTPCTOF),
+ fNSigmaPID(3),
+ fUseExclusiveNSigma(kFALSE),
+ fRemoveTracksT0Fill(kFALSE),
+fSelectCharge(0),
+fTriggerSelectCharge(0),
+fAssociatedSelectCharge(0),
+fTriggerRestrictEta(-1),
+fEtaOrdering(kFALSE),
+fCutConversions(kFALSE),
+fCutResonances(kFALSE),
+fRejectResonanceDaughters(-1),
+ fOnlyOneEtaSide(0),
+fInjectedSignals(kFALSE),
+ fRemoveWeakDecays(kFALSE),
+fRemoveDuplicates(kFALSE),
+ fapplyTrigefficiency(kFALSE),
+ fapplyAssoefficiency(kFALSE),
+ ffillefficiency(kFALSE),
+ fmesoneffrequired(kFALSE),
+ fkaonprotoneffrequired(kFALSE),
+ //fAnalysisUtils(0x0),
+ fDCAXYCut(0)
+{
+
+ for ( Int_t i = 0; i < 16; i++) {
+ fHistQA[i] = NULL;
+ }
+
+for ( Int_t i = 0; i < 6; i++ ){
+ fTHnrecomatchedallPid[i] = NULL;
+ fTHngenprimPidTruth[i] = NULL;
+ effcorection[i]=NULL;
+ //effmap[i]=NULL;
+
+ }
+
+ for(Int_t ipart=0;ipart<NSpecies;ipart++)
+ for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++)
+ fnsigmas[ipart][ipid]=999.;
+
+ for(Int_t ipart=0;ipart<NSpecies;ipart++) {fHasDoubleCounting[ipart]=kFALSE;}
+
+ // The last in the above list should not have a comma after it
+
+ // Constructor
+ // Define input and output slots here (never in the dummy constructor)
+ // Input slot #0 works with a TChain - it is connected to the default input container
+ // Output slot #1 writes into a TH1 container
+
+ DefineOutput(1, TList::Class()); // for output list
+
+}
+
+//________________________________________________________________________
+AliTwoParticlePIDCorr::~AliTwoParticlePIDCorr()
+{
+ // Destructor. Clean-up the output list, but not the histograms that are put inside
+ // (the list is owner and will clean-up these histograms). Protect in PROOF case.
+ if (fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
+ delete fOutput;
+
+ }
+ if (fPID) delete fPID;
+
+ }
+//________________________________________________________________________
+Float_t AliTwoParticlePIDCorr::PhiRange(Float_t DPhi)
+
+{
+ //
+ // Puts the argument in the range [-pi/2,3 pi/2].
+ //
+
+ if (DPhi < -TMath::Pi()/2) DPhi += 2*TMath::Pi();
+ if (DPhi > 3*TMath::Pi()/2) DPhi -= 2*TMath::Pi();
+
+ return DPhi;
+
+}
+//________________________________________________________________________
+void AliTwoParticlePIDCorr::UserCreateOutputObjects()
+{
+ // Create histograms
+ // Called once (on the worker node)
+ AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
+ AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
+ fPID = inputHandler->GetPIDResponse();
+
+ //AliAnalysisUtils *fUtils = new AliAnalysisUtils();
+
+//get the efficiency correction map
+
+
+ fOutput = new TList();
+ fOutput->SetOwner(); // IMPORTANT!
+
+ Int_t centmultbins=10;
+ Double_t centmultmin=0.0;
+ Double_t centmultmax=100.0;
+ if(fSampleType=="pPb" || fSampleType=="PbPb")
+ {
+ centmultbins=10;
+ centmultmin=0.0;
+ centmultmax=100.0;
+ }
+ if(fSampleType=="pp")
+ {
+ centmultbins=10;
+ centmultmin=0.0;
+ centmultmax=200.0;
+ }
+
+fhistcentrality=new TH1F("fhistcentrality","fhistcentrality",centmultbins*4,centmultmin,centmultmax);
+fOutput->Add(fhistcentrality);
+
+ fEventCounter = new TH1F("fEventCounter","EventCounter", 10, 0.5,10.5);
+ fEventCounter->GetXaxis()->SetBinLabel(1,"Event Accesed");
+ fEventCounter->GetXaxis()->SetBinLabel(2,"Have a vertex");
+ fEventCounter->GetXaxis()->SetBinLabel(5,"After vertex Cut");
+ fEventCounter->GetXaxis()->SetBinLabel(6,"Within 0-100% centrality");
+ fEventCounter->GetXaxis()->SetBinLabel(7,"Event Analyzed");
+ //fEventCounter->GetXaxis()->SetBinLabel(8,"Event Analysis finished");
+ fOutput->Add(fEventCounter);
+
+fEtaSpectrasso=new TH2F("fEtaSpectraasso","fEtaSpectraasso",180,-0.9,0.9,100,0.,20. );
+fOutput->Add(fEtaSpectrasso);
+
+fphiSpectraasso=new TH2F("fphiSpectraasso","fphiSpectraasso",72,0,2*TMath::Pi(),100,0.,20.);
+fOutput->Add(fphiSpectraasso);
+
+
+ if(fSampleType=="pPb" || fSampleType=="PbPb"){ fCentralityCorrelation = new TH2D("fCentralityCorrelation", ";centrality;multiplicity", 101, 0, 101, 20000, 0,40000);
+ fOutput->Add(fCentralityCorrelation);
+ }
+
+fHistoTPCdEdx = new TH2F("hHistoTPCdEdx", ";p_{T} (GeV/c);dE/dx (au.)",200,0.0,10.0,500, 0., 500.);
+fOutput->Add(fHistoTPCdEdx);
+fHistoTOFbeta = new TH2F(Form("hHistoTOFbeta"), ";p_{T} (GeV/c);v/c",70, 0., 7., 500, 0.1, 1.1);
+ fOutput->Add(fHistoTOFbeta);
+
+ fTPCTOFPion3d=new TH3F ("fTPCTOFpion3d", "fTPCTOFpion3d",100,0., 10., 120,-60.,60.,120,-60.,60);
+ fOutput->Add(fTPCTOFPion3d);
+
+ fTPCTOFKaon3d=new TH3F ("fTPCTOFKaon3d", "fTPCTOFKaon3d",100,0., 10., 120,-60.,60.,120,-60.,60);
+ fOutput->Add(fTPCTOFKaon3d);
+
+ fTPCTOFProton3d=new TH3F ("fTPCTOFProton3d", "fTPCTOFProton3d",100,0., 10., 120,-60.,60.,120,-60.,60);
+ fOutput->Add(fTPCTOFProton3d);
+
+if(ffillhistQAReco)
+ {
+ fPionPt = new TH1F("fHistQAPionPt","p_{T} distribution",200,0.,10.);
+ fOutput->Add(fPionPt);
+ fPionEta= new TH1F("fHistQAPionEta","#eta distribution",360,-1.8,1.8);
+ fOutput->Add(fPionEta);
+ fPionPhi = new TH1F("fHistQAPionPhi","#phi distribution",340,0,6.8);
+ fOutput->Add(fPionPhi);
+
+ fKaonPt = new TH1F("fHistQAKaonPt","p_{T} distribution",200,0.,10.);
+ fOutput->Add(fKaonPt);
+ fKaonEta= new TH1F("fHistQAKaonEta","#eta distribution",360,-1.8,1.8);
+ fOutput->Add(fKaonEta);
+ fKaonPhi = new TH1F("fHistQAKaonPhi","#phi distribution",340,0,6.8);
+ fOutput->Add(fKaonPhi);
+
+ fProtonPt = new TH1F("fHistQAProtonPt","p_{T} distribution",200,0.,10.);
+ fOutput->Add(fProtonPt);
+ fProtonEta= new TH1F("fHistQAProtonEta","#eta distribution",360,-1.8,1.8);
+ fOutput->Add(fProtonEta);
+ fProtonPhi= new TH1F("fHistQAProtonPhi","#phi distribution",340,0,6.8);
+ fOutput->Add(fProtonPhi);
+ }
+
+ fHistQA[0] = new TH1F("fHistQAvx", "Histo Vx All ", 50, -5., 5.);
+ fHistQA[1] = new TH1F("fHistQAvy", "Histo Vy All", 50, -5., 5.);
+ fHistQA[2] = new TH1F("fHistQAvz", "Histo Vz All", 50, -25., 25.);
+ fHistQA[3] = new TH1F("fHistQAvxA", "Histo Vx After Cut ", 50, -5., 5.);
+ fHistQA[4] = new TH1F("fHistQAvyA", "Histo Vy After Cut", 50, -5., 5.);
+ fHistQA[5] = new TH1F("fHistQAvzA", "Histo Vz After Cut", 50, -25., 25.);
+ fHistQA[6] = new TH1F("fHistQADcaXyC", "Histo DCAxy after cut", 50, -5., 5.);
+ fHistQA[7] = new TH1F("fHistQADcaZC", "Histo DCAz after cut", 50, -5., 5.);
+ fHistQA[8] = new TH1F("fHistQAPt","p_{T} distribution",200,0.,10.);
+ fHistQA[9] = new TH1F("fHistQAEta","#eta distribution",360,-1.8,1.8);
+ fHistQA[10] = new TH1F("fHistQAPhi","#phi distribution",340,0,6.8);
+ fHistQA[11] = new TH1F("fHistQANCls","Number of TPC cluster",200,0,200);
+ fHistQA[13] = new TH1F("fHistQAChi2","Chi2 per NDF",100,0,10);
+ fHistQA[12] = new TH1F("fHistQANCls1","Number of TPC cluster1",200,0,200);
+ fHistQA[14] = new TH1F("nCrossedRowsTPC","Number of TPC ccrossed rows",200,0,200);
+ fHistQA[15] = new TH1F("ratioCrossedRowsOverFindableClustersTPC","Number of TPC ccrossed rows find clusters",200,0,2);
+
+for(Int_t i = 0; i < 16; i++)
+ {
+ fOutput->Add(fHistQA[i]);
+ }
+
+ kTrackVariablesPair=6 ;
+
+ if(fcontainPIDtrig && !fcontainPIDasso) kTrackVariablesPair=7;
+
+ if(!fcontainPIDtrig && fcontainPIDasso) kTrackVariablesPair=7;
+
+ if(fcontainPIDtrig && fcontainPIDasso) kTrackVariablesPair=8;
+
+
+// two particle histograms
+ Int_t iBinPair[kTrackVariablesPair]; // binning for track variables
+ Double_t* dBinsPair[kTrackVariablesPair]; // bins for track variables
+ TString* axisTitlePair; // axis titles for track variables
+ axisTitlePair=new TString[kTrackVariablesPair];
+
+ TString defaultBinningStr;
+ defaultBinningStr = "eta: -1.0, -0.9, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0\n"
+ "p_t_assoc: 0.5, 0.75, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 8.0,10.0\n"
+ "p_t_leading_course: 0.5, 1.0, 2.0, 3.0, 4.0, 6.0, 8.0,10.0\n"
+ "p_t_eff:0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.25, 1.5, 1.75, 2.0, 2.25, 2.5, 2.75, 3.0, 3.25, 3.5, 3.75, 4.0, 4.5, 5.0,5.5, 6.0, 7.0, 8.0,9.0,10.0\n"
+ "vertex: -10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10\n"
+ "delta_phi: -1.570796, -1.483530, -1.396263, -1.308997, -1.221730, -1.134464, -1.047198, -0.959931, -0.872665, -0.785398, -0.698132, -0.610865, -0.523599, -0.436332, -0.349066, -0.261799, -0.174533, -0.087266, 0.0, 0.087266, 0.174533, 0.261799, 0.349066, 0.436332, 0.523599, 0.610865, 0.698132, 0.785398, 0.872665, 0.959931, 1.047198, 1.134464, 1.221730, 1.308997, 1.396263, 1.483530, 1.570796, 1.658063, 1.745329, 1.832596, 1.919862, 2.007129, 2.094395, 2.181662, 2.268928, 2.356194, 2.443461, 2.530727, 2.617994, 2.705260, 2.792527, 2.879793, 2.967060, 3.054326, 3.141593, 3.228859, 3.316126, 3.403392, 3.490659, 3.577925, 3.665191, 3.752458, 3.839724, 3.926991, 4.014257, 4.101524, 4.188790, 4.276057, 4.363323, 4.450590, 4.537856, 4.625123, 4.712389\n" // this binning starts at -pi/2 and is modulo 3
+ "delta_eta: -2.4, -2.3, -2.2, -2.1, -2.0, -1.9, -1.8, -1.7, -1.6, -1.5, -1.4, -1.3, -1.2, -1.1, -1.0, -0.9, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0,2.1, 2.2, 2.3, 2.4\n"
+ "multiplicity: 0, 1, 2, 3, 4, 5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100.1\n";
+
+ if(fcontainPIDtrig){
+ defaultBinningStr += "PIDTrig: -0.5,0.5,1.5,2.5,3.5\n"; // course
+ }
+ if(fcontainPIDasso){
+ defaultBinningStr += "PIDAsso: -0.5,0.5,1.5,2.5,3.5\n"; // course
+ }
+ // =========================================================
+ // Customization (adopted from AliUEHistograms)
+ // =========================================================
+
+ TObjArray* lines = defaultBinningStr.Tokenize("\n");
+ for (Int_t i=0; i<lines->GetEntriesFast(); i++)
+ {
+ TString line(lines->At(i)->GetName());
+ TString tag = line(0, line.Index(":")+1);
+ if (!fCustomBinning.BeginsWith(tag) && !fCustomBinning.Contains(TString("\n") + tag))
+ fBinningString += line + "\n";
+ else
+ AliInfo(Form("Using custom binning for %s", tag.Data()));
+ }
+ delete lines;
+ fBinningString += fCustomBinning;
+
+ AliInfo(Form("Used THn Binning:\n%s",fBinningString.Data()));
+
+ // =========================================================
+ // Now set the bins
+ // =========================================================
+
+ dBinsPair[0] = GetBinning(fBinningString, "multiplicity", iBinPair[0]);
+ axisTitlePair[0] = " multiplicity";
+
+ dBinsPair[1] = GetBinning(fBinningString, "vertex", iBinPair[1]);
+ axisTitlePair[1] = "v_{Z} (cm)";
+
+ dBinsPair[2] = GetBinning(fBinningString, "p_t_leading_course", iBinPair[2]);
+ axisTitlePair[2] = "p_{T,trig.} (GeV/c)";
+
+ dBinsPair[3] = GetBinning(fBinningString, "p_t_assoc", iBinPair[3]);
+ axisTitlePair[3] = "p_{T,assoc.} (GeV/c)";
+
+ dBinsPair[4] = GetBinning(fBinningString, "delta_eta", iBinPair[4]);
+ axisTitlePair[4] = "#Delta#eta";
+
+ dBinsPair[5] = GetBinning(fBinningString, "delta_phi", iBinPair[5]);
+ axisTitlePair[5] = "#Delta#varphi (rad)";
+
+ if(fcontainPIDtrig && !fcontainPIDasso){
+ dBinsPair[6] = GetBinning(fBinningString, "PIDTrig", iBinPair[6]);
+ axisTitlePair[6] = "PIDTrig";
+ }
+
+ if(!fcontainPIDtrig && fcontainPIDasso){
+ dBinsPair[6] = GetBinning(fBinningString, "PIDAsso", iBinPair[6]);
+ axisTitlePair[6] = "PIDAsso";
+ }
+
+if(fcontainPIDtrig && fcontainPIDasso){
+
+ dBinsPair[6] = GetBinning(fBinningString, "PIDTrig", iBinPair[6]);
+ axisTitlePair[6] = "PIDTrig";
+
+ dBinsPair[7] = GetBinning(fBinningString, "PIDAsso", iBinPair[7]);
+ axisTitlePair[7] = "PIDAsso";
+ }
+
+ Int_t nEtaBin = -1;
+ Double_t* EtaBin = GetBinning(fBinningString, "eta", nEtaBin);
+
+ Int_t nPteffbin = -1;
+ Double_t* Pteff = GetBinning(fBinningString, "p_t_eff", nPteffbin);
+
+
+ fminPtTrig=dBinsPair[2][0];
+ fmaxPtTrig=dBinsPair[2][iBinPair[2]];
+ fminPtAsso=dBinsPair[3][0];
+ fmaxPtAsso=dBinsPair[3][iBinPair[3]];
+
+ //fminPtComboeff=fminPtTrig;***then this value will be fixed ,even Setter can't change it's value
+ //fmaxPtComboeff=fmaxPtTrig;
+//THnSparses for calculation of efficiency
+
+ if((fAnalysisType =="MCAOD") && ffillefficiency) {
+ const Int_t nDim = 4;// cent zvtx pt eta
+ Int_t fBinsCh[nDim] = {iBinPair[0], iBinPair[1], nPteffbin ,nEtaBin};//**********************change it
+ Double_t fMinCh[nDim] = { dBinsPair[0][0],dBinsPair[1][0], Pteff[0], EtaBin[0] };
+ Double_t fMaxCh[nDim] = { dBinsPair[0][iBinPair[0]], dBinsPair[1][iBinPair[1]], Pteff[nPteffbin], EtaBin[nEtaBin]};
+
+TString Histrename;
+for(Int_t jj=0;jj<6;jj++)//PID type binning
+ {
+ Histrename="fTHnrecomatchedallPid";Histrename+=jj;
+ fTHnrecomatchedallPid[jj] = new THnSparseF(Histrename.Data(),"cent:zvtx::Pt:eta", nDim, fBinsCh, fMinCh, fMaxCh);
+ fTHnrecomatchedallPid[jj]->Sumw2();
+ fTHnrecomatchedallPid[jj]->GetAxis(0)->Set(iBinPair[0], dBinsPair[0]);
+ fTHnrecomatchedallPid[jj]->GetAxis(0)->SetTitle("Centrality");
+ fTHnrecomatchedallPid[jj]->GetAxis(1)->Set(iBinPair[1],dBinsPair[1]);
+ fTHnrecomatchedallPid[jj]->GetAxis(1)->SetTitle("zvtx");
+ fTHnrecomatchedallPid[jj]->GetAxis(2)->Set(nPteffbin, Pteff);
+ fTHnrecomatchedallPid[jj]->GetAxis(2)->SetTitle("p_{T} (GeV/c)");
+ fTHnrecomatchedallPid[jj]->GetAxis(3)->Set(nEtaBin,EtaBin);
+ fTHnrecomatchedallPid[jj]->GetAxis(3)->SetTitle("#eta");
+ fOutput->Add(fTHnrecomatchedallPid[jj]);
+
+Histrename="fTHngenprimPidTruth";Histrename+=jj;
+ fTHngenprimPidTruth[jj] = new THnSparseF(Histrename.Data(),"cent:zvtx::Pt:eta", nDim, fBinsCh, fMinCh, fMaxCh);
+ fTHngenprimPidTruth[jj]->Sumw2();
+ fTHngenprimPidTruth[jj]->GetAxis(0)->Set(iBinPair[0],dBinsPair[0]);
+ fTHngenprimPidTruth[jj]->GetAxis(0)->SetTitle("Centrality");
+ fTHngenprimPidTruth[jj]->GetAxis(1)->Set(iBinPair[1], dBinsPair[1]);
+ fTHngenprimPidTruth[jj]->GetAxis(1)->SetTitle("zvtx");
+ fTHngenprimPidTruth[jj]->GetAxis(2)->Set(nPteffbin, Pteff);
+ fTHngenprimPidTruth[jj]->GetAxis(2)->SetTitle("p_{T} (GeV/c)");
+ fTHngenprimPidTruth[jj]->GetAxis(3)->Set(nEtaBin, EtaBin);
+ fTHngenprimPidTruth[jj]->GetAxis(3)->SetTitle("#eta");
+ fOutput->Add(fTHngenprimPidTruth[jj]);
+ }
+ }
+
+ Int_t fBins[kTrackVariablesPair];
+ Double_t fMin[kTrackVariablesPair];
+ Double_t fMax[kTrackVariablesPair];
+
+//ThnSparses for Correlation plots(data & MC)
+ for(Int_t i=0;i<kTrackVariablesPair;i++)
+ {
+ fBins[i] =iBinPair[i];
+ fMin[i]= dBinsPair[i][0];
+ fMax[i] = dBinsPair[i][iBinPair[i]];
+ }
+ if(ffilltrigassoUNID)
+ {
+ fTHnCorrUNID = new THnSparseF("fTHnCorrUNID","cent:zvtx:pttrig:ptasso:deltaeta:deltaphi", kTrackVariablesPair, fBins, fMin, fMax);
+ fTHnCorrUNID->Sumw2();
+ for(Int_t i=0; i<kTrackVariablesPair;i++){
+ SetAsymmetricBin(fTHnCorrUNID,i,dBinsPair[i],iBinPair[i],axisTitlePair[i]);
+ }
+ fOutput->Add(fTHnCorrUNID);
+
+fTHnCorrUNIDmix = new THnSparseF("fTHnCorrUNIDmix","cent:zvtx:pttrig:ptasso:deltaeta:deltaphi", kTrackVariablesPair, fBins, fMin, fMax);
+ fTHnCorrUNIDmix->Sumw2();
+for(Int_t i=0; i<kTrackVariablesPair;i++){
+ SetAsymmetricBin(fTHnCorrUNIDmix,i,dBinsPair[i],iBinPair[i],axisTitlePair[i]);
+ }
+ fOutput->Add(fTHnCorrUNIDmix);
+ }
+
+ if(ffilltrigIDassoID)
+ {
+ fTHnCorrID = new THnSparseF("fTHnCorrID","cent:zvtx:pttrig:ptasso:deltaeta:deltaphi", kTrackVariablesPair, fBins, fMin, fMax);
+ fTHnCorrID->Sumw2();
+for(Int_t i=0; i<kTrackVariablesPair;i++){
+ SetAsymmetricBin(fTHnCorrID,i,dBinsPair[i],iBinPair[i],axisTitlePair[i]);
+ }
+ fOutput->Add(fTHnCorrID);
+
+fTHnCorrIDmix = new THnSparseF("fTHnCorrIDmix","cent:zvtx:pttrig:ptasso:deltaeta:deltaphi", kTrackVariablesPair, fBins, fMin, fMax);
+ fTHnCorrIDmix->Sumw2();
+for(Int_t i=0; i<kTrackVariablesPair;i++){
+ SetAsymmetricBin(fTHnCorrIDmix,i,dBinsPair[i],iBinPair[i],axisTitlePair[i]);
+ }
+ fOutput->Add(fTHnCorrIDmix);
+ }
+
+ if(ffilltrigUNIDassoID || ffilltrigIDassoUNID)//***********a bit tricky, be careful
+ {
+ fTHnCorrIDUNID = new THnSparseF("fTHnCorrIDUNID","cent:zvtx:pttrig:ptasso:deltaeta:deltaphi", kTrackVariablesPair, fBins, fMin, fMax);
+ fTHnCorrIDUNID->Sumw2();
+for(Int_t i=0; i<kTrackVariablesPair;i++){
+ SetAsymmetricBin(fTHnCorrIDUNID,i,dBinsPair[i],iBinPair[i],axisTitlePair[i]);
+ }
+ fOutput->Add(fTHnCorrIDUNID);
+
+ fTHnCorrIDUNIDmix = new THnSparseF("fTHnCorrIDUNIDmix","cent:zvtx:pttrig:ptasso:deltaeta:deltaphi", kTrackVariablesPair, fBins, fMin, fMax);
+ fTHnCorrIDUNIDmix->Sumw2();
+for(Int_t i=0; i<kTrackVariablesPair;i++){
+ SetAsymmetricBin(fTHnCorrIDUNIDmix,i,dBinsPair[i],iBinPair[i],axisTitlePair[i]);
+ }
+ fOutput->Add(fTHnCorrIDUNIDmix);
+ }
+
+
+
+ //ThnSparse for Correlation plots(truth MC)
+ if((fAnalysisType == "MCAOD") && ffilltrigIDassoIDMCTRUTH) {//remember that in this case uidentified means other than pions, kaons, protons
+ fCorrelatonTruthPrimary = new THnSparseF("fCorrelatonTruthPrimary","cent:zvtx:pttrig:ptasso:deltaeta:deltaphi", kTrackVariablesPair, fBins, fMin, fMax);
+ fCorrelatonTruthPrimary->Sumw2();
+for(Int_t i=0; i<kTrackVariablesPair;i++){
+ SetAsymmetricBin(fCorrelatonTruthPrimary,i,dBinsPair[i],iBinPair[i],axisTitlePair[i]);
+ }
+ fOutput->Add(fCorrelatonTruthPrimary);
+
+ fCorrelatonTruthPrimarymix = new THnSparseF("fCorrelatonTruthPrimarymix","cent:zvtx:pttrig:ptasso:deltaeta:deltaphi", kTrackVariablesPair, fBins, fMin, fMax);
+ fCorrelatonTruthPrimarymix->Sumw2();
+for(Int_t i=0; i<kTrackVariablesPair;i++){
+ SetAsymmetricBin(fCorrelatonTruthPrimarymix,i,dBinsPair[i],iBinPair[i],axisTitlePair[i]);
+ }
+ fOutput->Add(fCorrelatonTruthPrimarymix);
+ }
+
+
+ //binning for trigger no. counting
+
+ Double_t* fMint;
+ Double_t* fMaxt;
+ Int_t* fBinst;
+ Int_t dims=3;
+ if(fcontainPIDtrig) dims=4;
+ fMint= new Double_t[dims];
+ fMaxt= new Double_t[dims];
+ fBinst= new Int_t[dims];
+ for(Int_t i=0; i<3;i++)
+ {
+ fBinst[i]=iBinPair[i];
+ fMint[i]=dBinsPair[i][0];
+ fMaxt[i]=dBinsPair[i][iBinPair[i]];
+ }
+ if(fcontainPIDtrig){
+ fBinst[3]=iBinPair[6];
+ fMint[3]=dBinsPair[6][0];
+ fMaxt[3]=dBinsPair[6][iBinPair[6]];
+ }
+
+ //ThSparse for trigger counting(data & reco MC)
+ if(ffilltrigassoUNID || ffilltrigUNIDassoID || ffilltrigIDassoUNID || ffilltrigIDassoID)
+ {
+ fTHnTrigcount = new THnSparseF("fTHnTrigcount","cent:zvtx:pt", dims, fBinst, fMint, fMaxt);
+ fTHnTrigcount->Sumw2();
+for(Int_t i=0; i<3;i++){
+ SetAsymmetricBin(fTHnTrigcount,i,dBinsPair[i],iBinPair[i],axisTitlePair[i]);
+ }
+ if(fcontainPIDtrig) SetAsymmetricBin(fTHnTrigcount,3,dBinsPair[6],iBinPair[6],axisTitlePair[6]);
+ fOutput->Add(fTHnTrigcount);
+ }
+
+ if((fAnalysisType =="MCAOD") && ffilltrigIDassoIDMCTRUTH) {
+ //ThSparse for trigger counting(truth MC)
+fTHnTrigcountMCTruthPrim = new THnSparseF("fTHnTrigcountMCTruthPrim","cent:zvtx:pt", dims, fBinst, fMint, fMaxt);
+ fTHnTrigcountMCTruthPrim->Sumw2();
+for(Int_t i=0; i<3;i++){
+ SetAsymmetricBin(fTHnTrigcountMCTruthPrim,i,dBinsPair[i],iBinPair[i],axisTitlePair[i]);
+ }
+if(fcontainPIDtrig) SetAsymmetricBin(fTHnTrigcountMCTruthPrim,3,dBinsPair[6],iBinPair[6],axisTitlePair[6]);
+ fOutput->Add(fTHnTrigcountMCTruthPrim);
+ }
+
+if(fAnalysisType=="MCAOD"){
+ if(ffillhistQATruth)
+ {
+ MCtruthpt=new TH1F ("MCtruthpt","ptdistributiontruthprim",100,0.,10.);
+ fOutput->Add(MCtruthpt);
+
+ MCtrutheta=new TH1F ("MCtrutheta","etadistributiontruthprim",360,-1.8,1.8);
+ fOutput->Add(MCtrutheta);
+
+ MCtruthphi=new TH1F ("MCtruthphi","phidisttruthprim",340,0,6.8);
+ fOutput->Add(MCtruthphi);
+
+ MCtruthpionpt=new TH1F ("MCtruthpionpt","MCtruthpionpt",100,0.,10.);
+ fOutput->Add(MCtruthpionpt);
+
+ MCtruthpioneta=new TH1F ("MCtruthpioneta","MCtruthpioneta",360,-1.8,1.8);
+ fOutput->Add(MCtruthpioneta);
+
+ MCtruthpionphi=new TH1F ("MCtruthpionphi","MCtruthpionphi",340,0,6.8);
+ fOutput->Add(MCtruthpionphi);
+
+ MCtruthkaonpt=new TH1F ("MCtruthkaonpt","MCtruthkaonpt",100,0.,10.);
+ fOutput->Add(MCtruthkaonpt);
+
+ MCtruthkaoneta=new TH1F ("MCtruthkaoneta","MCtruthkaoneta",360,-1.8,1.8);
+ fOutput->Add(MCtruthkaoneta);
+
+ MCtruthkaonphi=new TH1F ("MCtruthkaonphi","MCtruthkaonphi",340,0,6.8);
+ fOutput->Add(MCtruthkaonphi);
+
+ MCtruthprotonpt=new TH1F ("MCtruthprotonpt","MCtruthprotonpt",100,0.,10.);
+ fOutput->Add(MCtruthprotonpt);
+
+ MCtruthprotoneta=new TH1F ("MCtruthprotoneta","MCtruthprotoneta",360,-1.8,1.8);
+ fOutput->Add(MCtruthprotoneta);
+
+ MCtruthprotonphi=new TH1F ("MCtruthprotonphi","MCtruthprotonphi",340,0,6.8);
+ fOutput->Add(MCtruthprotonphi);
+ }
+ fPioncont=new TH2F("fPioncont", "fPioncont",10,-0.5,9.5,100,0.,10.);
+ fOutput->Add(fPioncont);
+
+ fKaoncont=new TH2F("fKaoncont","fKaoncont",10,-0.5,9.5,100,0.,10.);
+ fOutput->Add(fKaoncont);
+
+ fProtoncont=new TH2F("fProtoncont","fProtoncont",10,-0.5,9.5,100,0.,10.);
+ fOutput->Add(fProtoncont);
+ }
+
+//Mixing
+DefineEventPool();
+
+ if(fapplyTrigefficiency || fapplyAssoefficiency)
+ {
+ const Int_t nDimt = 4;// cent zvtx pt eta
+ Int_t fBinsCht[nDimt] = {iBinPair[0], iBinPair[1], nPteffbin ,nEtaBin};//*************change it
+ Double_t fMinCht[nDimt] = { dBinsPair[0][0],dBinsPair[1][0], Pteff[0], EtaBin[0] };
+ Double_t fMaxCht[nDimt] = {dBinsPair[0][iBinPair[0]], dBinsPair[1][iBinPair[1]], Pteff[nPteffbin], EtaBin[nEtaBin]};
+
+ TString Histrexname;
+for(Int_t jj=0;jj<6;jj++)// PID type binning
+ {
+ Histrexname="effcorection";Histrexname+=jj;
+ effcorection[jj] = new THnSparseF(Histrexname.Data(),"cent:zvtx::Pt:eta", nDimt, fBinsCht, fMinCht, fMaxCht);
+ effcorection[jj]->Sumw2();
+ effcorection[jj]->GetAxis(0)->Set(iBinPair[0], dBinsPair[0]);
+ effcorection[jj]->GetAxis(0)->SetTitle("Centrality");
+ effcorection[jj]->GetAxis(1)->Set( iBinPair[1],dBinsPair[1]);
+ effcorection[jj]->GetAxis(1)->SetTitle("zvtx");
+ effcorection[jj]->GetAxis(2)->Set(nPteffbin, Pteff);
+ effcorection[jj]->GetAxis(2)->SetTitle("p_{T} (GeV/c)");
+ effcorection[jj]->GetAxis(3)->Set( nEtaBin,EtaBin);
+ effcorection[jj]->GetAxis(3)->SetTitle("#eta");
+ fOutput->Add(effcorection[jj]);
+ }
+// TFile *fsifile = new TFile(fefffilename,"READ");
+ TFile *fileT=TFile::Open(fefffilename);
+ TString Nameg;
+for(Int_t jj=0;jj<6;jj++)//type binning
+ {
+Nameg="effmap";Nameg+=jj;
+//effcorection[jj] = (THnSparseF*)fsifile->Get(Nameg.Data());
+effcorection[jj] = (THnSparseF*)fileT->Get(Nameg.Data());
+
+//effcorection[jj]->SetDirectory(0);//****************************not present in case oh THnF
+ }
+//fsifile->Close();
+fileT->Close();
+
+ }
+
+//fControlConvResoncances = new TH2F("fControlConvResoncances", ";id;delta mass", 3, -0.5, 2.5, 100, -0.1, 0.1);
+// fOutput->Add(fControlConvResoncances);
+
+
+ PostData(1, fOutput); // Post data for ALL output slots >0 here, to get at least an empty histogram
+}
+//-------------------------------------------------------------------------------
+void AliTwoParticlePIDCorr::UserExec( Option_t * ){
+
+
+ if(fAnalysisType == "AOD") {
+
+ doAODevent();
+
+ }//AOD--analysis-----
+
+ else if(fAnalysisType == "MCAOD") {
+
+ doMCAODevent();
+
+ }
+
+ else return;
+
+}
+//-------------------------------------------------------------------------
+void AliTwoParticlePIDCorr::doMCAODevent()
+{
+ AliVEvent *event = InputEvent();
+ if (!event) { Printf("ERROR: Could not retrieve event"); return; }
+ AliAODEvent* aod = dynamic_cast<AliAODEvent*>(event);
+ if (!aod) {
+ AliError("Cannot get the AOD event");
+ return;
+ }
+
+// count all events(physics triggered)
+ fEventCounter->Fill(1);
+
+ // get centrality object and check quality(valid for p-Pb and Pb-Pb)
+ Double_t cent_v0=0.0;
+
+ if(fSampleType=="pPb" || fSampleType=="PbPb")
+ {
+ AliCentrality *centrality=0;
+ if(aod)
+ centrality = aod->GetHeader()->GetCentralityP();
+ // if (centrality->GetQuality() != 0) return ;
+
+ if(centrality)
+ {
+ cent_v0 = centrality->GetCentralityPercentile(fCentralityMethod);
+ }
+ else
+ {
+ cent_v0= -1;
+ }
+ }
+
+//check the PIDResponse handler
+ if (!fPID) return;
+
+// get mag. field required for twotrack efficiency cut
+ Float_t bSign = 0;
+ bSign = (aod->GetMagneticField() > 0) ? 1 : -1;
+
+ //check for TClonesArray(truth track MC information)
+fArrayMC = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
+ if (!fArrayMC) {
+ AliFatal("Error: MC particles branch not found!\n");
+ return;
+ }
+
+ //check for AliAODMCHeader(truth event MC information)
+ AliAODMCHeader *header=NULL;
+ header=(AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
+ if(!header) {
+ printf("MC header branch not found!\n");
+ return;
+ }
+
+//Only consider MC events within the vtx-z region used also as cut on the reconstructed vertex
+Float_t zVtxmc =header->GetVtxZ();
+ if(TMath::Abs(zVtxmc)>fzvtxcut) return;
+
+ // For productions with injected signals, figure out above which label to skip particles/tracks
+ Int_t skipParticlesAbove = 0;
+
+ if (fInjectedSignals)
+ {
+ AliGenEventHeader* eventHeader = 0;
+ Int_t headers = 0;
+
+// AOD
+ if (!header)
+ AliFatal("fInjectedSignals set but no MC header found");
+
+ headers = header->GetNCocktailHeaders();
+ eventHeader = header->GetCocktailHeader(0);
+
+ if (!eventHeader)
+ {
+ // We avoid AliFatal here, because the AOD productions sometimes have events where the MC header is missing
+ // (due to unreadable Kinematics) and we don't want to loose the whole job because of a few events
+ AliError("First event header not found. Skipping this event.");
+ //fHistos->FillEvent(centrality, AliUEHist::kCFStepAnaTopology);
+ return;
+ }
+skipParticlesAbove = eventHeader->NProduced();
+ AliInfo(Form("Injected signals in this event (%d headers). Keeping events of %s. Will skip particles/tracks above %d.", headers, eventHeader->ClassName(), skipParticlesAbove));
+ }
+
+ //vertex selection(is it fine for PP?)
+ if ( fSampleType=="pPb"){
+ trkVtx = aod->GetPrimaryVertex();
+ if (!trkVtx || trkVtx->GetNContributors()<=0) return;
+ TString vtxTtl = trkVtx->GetTitle();
+ if (!vtxTtl.Contains("VertexerTracks")) return;
+ zvtx = trkVtx->GetZ();
+ const AliAODVertex* spdVtx = aod->GetPrimaryVertexSPD();
+ if (!spdVtx || spdVtx->GetNContributors()<=0) return;
+ TString vtxTyp = spdVtx->GetTitle();
+ Double_t cov[6]={0};
+ spdVtx->GetCovarianceMatrix(cov);
+ Double_t zRes = TMath::Sqrt(cov[5]);
+ if (vtxTyp.Contains("vertexer:Z") && (zRes>0.25)) return;
+ if (TMath::Abs(spdVtx->GetZ() - trkVtx->GetZ())>0.5) return;
+ }
+ else if(fSampleType=="PbPb" || fSampleType=="pp") {//for pp and pb-pb case
+ Int_t nVertex = aod->GetNumberOfVertices();
+ if( nVertex > 0 ) {
+ trkVtx = aod->GetPrimaryVertex();
+ Int_t nTracksPrim = trkVtx->GetNContributors();
+ zvtx = trkVtx->GetZ();
+ //if (fDebug > 1)AliInfo(Form(" Vertex in = %f with %d particles by %s data ...",zVertex,nTracksPrim,vertex->GetName()));
+ // Reject TPC only vertex
+ TString name(trkVtx->GetName());
+ if (name.CompareTo("PrimaryVertex") && name.CompareTo("SPDVertex"))return;
+
+ // Select a quality vertex by number of tracks?
+ if( nTracksPrim < fnTracksVertex ) {
+ //if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
+ return ;
+ }
+ // TODO remove vertexer Z events with dispersion > 0.02: Doesn't work for AOD at present
+ //if (strcmp(vertex->GetTitle(), "AliVertexerZ") == 0 && vertex->GetDispersion() > 0.02)
+ // return kFALSE;
+ // if (fDebug > 1) AliInfo(" Primary-vertex Selection: event ACCEPTED...");
+ }
+ else return;
+
+ }
+ else return;//as there is no proper sample type
+
+
+ fHistQA[0]->Fill((trkVtx->GetX()));fHistQA[1]->Fill((trkVtx->GetY()));fHistQA[2]->Fill((trkVtx->GetZ())); //for trkVtx only before vertex cut |zvtx|<10 cm
+
+ //count events having a proper vertex
+ fEventCounter->Fill(2);
+
+ if (TMath::Abs(zvtx) > fzvtxcut) return;
+
+fHistQA[3]->Fill((trkVtx->GetX()));fHistQA[4]->Fill((trkVtx->GetY()));fHistQA[5]->Fill((trkVtx->GetZ()));//after vertex cut for trkVtx only
+
+//now we have events passed physics trigger, centrality,eventzvtx cut
+
+//count events after vertex cut
+ fEventCounter->Fill(5);
+
+if(!aod) return; //for safety
+
+if (fSampleType=="pPb" || fSampleType=="PbPb") if (cent_v0 < 0) return;//for pp case it is the multiplicity
+
+ Double_t nooftrackstruth=0.0;//in case of pp this will give the multiplicity(for truth case) after the track loop(only for unidentified particles that pass kinematic cuts)
+
+ Double_t cent_v0_truth=0.0;
+
+ //TObjArray* tracksMCtruth=0;
+TObjArray* tracksMCtruth=new TObjArray;//for truth MC particles with PID,here unidentified means any particle other than pion, kaon or proton(Basicaly Spundefined of AliHelperPID)******WARNING::different from data and reco MC
+ tracksMCtruth->SetOwner(kTRUE); //***********************************IMPORTANT!
+
+ eventno++;
+
+ //There is a small difference between zvtx and zVtxmc??????
+ //cout<<"***********************************************zVtxmc="<<zVtxmc<<endl;
+ //cout<<"***********************************************zvtx="<<zvtx<<endl;
+
+//now process the truth particles(for both efficiency & correlation function)
+Int_t nMCTrack = fArrayMC->GetEntriesFast();
+
+for (Int_t iMC = 0; iMC < nMCTrack; iMC++)
+{ //MC truth track loop starts
+
+AliAODMCParticle *partMC = (AliAODMCParticle*) fArrayMC->At(iMC);
+
+ if(!partMC){
+ AliError(Form("ERROR: Could not retrieve AODMCtrack %d",iMC));
+ continue;
+ }
+
+//consider only charged particles
+ if(partMC->Charge() == 0) continue;
+
+//consider only primary particles; neglect all secondary particles including from weak decays
+ if(fselectprimaryTruth && !partMC->IsPhysicalPrimary()) continue;
+
+
+//remove injected signals(primaries above <maxLabel>)
+ if (fInjectedSignals && partMC->GetLabel() >= skipParticlesAbove) continue;
+
+//remove duplicates
+ Bool_t isduplicate=kFALSE;
+ if (fRemoveDuplicates)
+ {
+ for (Int_t j=iMC+1; j<nMCTrack; ++j)
+ {//2nd trutuh loop starts
+AliAODMCParticle *partMC2 = (AliAODMCParticle*) fArrayMC->At(j);
+ if(!partMC2){
+ AliError(Form("ERROR: Could not retrieve AODMCtrack %d",j));
+ continue;
+ }
+ if (partMC->GetLabel() == partMC2->GetLabel())
+ {
+isduplicate=kTRUE;
+ break;
+ }
+ }//2nd truth loop ends
+ }
+ if(fRemoveDuplicates && isduplicate) continue;//remove duplicates
+
+//give only kinematic cuts at the truth level
+ if (partMC->Eta() < fmineta || partMC->Eta() > fmaxeta) continue;
+ if (partMC->Pt() < fminPt || partMC->Pt() > fmaxPt) continue;
+
+ if(!partMC) continue;//for safety
+
+ //To determine multiplicity in case of PP
+ nooftrackstruth++;
+ //cout<<"**************************************"<<TMath::Abs(partMC->GetLabel())<<endl;
+//only physical primary(all/unidentified)
+if(ffillhistQATruth)
+ {
+ MCtruthpt->Fill(partMC->Pt());
+ MCtrutheta->Fill(partMC->Eta());
+ MCtruthphi->Fill(partMC->Phi());
+ }
+ //get particle ID
+Int_t pdgtruth=((AliAODMCParticle*)partMC)->GetPdgCode();
+Int_t particletypeTruth=-999;
+ if (TMath::Abs(pdgtruth)==211)
+ {
+ particletypeTruth=SpPion;
+if(ffillhistQATruth)
+ {
+ MCtruthpionpt->Fill(partMC->Pt());
+ MCtruthpioneta->Fill(partMC->Eta());
+ MCtruthpionphi->Fill(partMC->Phi());
+ }
+ }
+ if (TMath::Abs(pdgtruth)==321)
+ {
+ particletypeTruth=SpKaon;
+if(ffillhistQATruth)
+ {
+ MCtruthkaonpt->Fill(partMC->Pt());
+ MCtruthkaoneta->Fill(partMC->Eta());
+ MCtruthkaonphi->Fill(partMC->Phi());
+ }
+ }
+if(TMath::Abs(pdgtruth)==2212)
+ {
+ particletypeTruth=SpProton;
+if(ffillhistQATruth)
+ {
+ MCtruthprotonpt->Fill(partMC->Pt());
+ MCtruthprotoneta->Fill(partMC->Eta());
+ MCtruthprotonphi->Fill(partMC->Phi());
+ }
+ }
+ if(TMath::Abs(pdgtruth)!=211 && TMath::Abs(pdgtruth)!=321 && TMath::Abs(pdgtruth)!=2212) particletypeTruth=unidentified;//*********************WARNING:: situation is different from reco MC and data case(here we don't have SpUndefined particles,because here unidentified=SpUndefined)
+
+ // -- Fill THnSparse for efficiency and contamination calculation
+ if (fSampleType=="pp") cent_v0=15.0;//integrated over multiplicity(so put any fixed value for each track so that practically means there is only one bin in multiplicity i.e. multiplicity intregated out )**************Important
+
+ Double_t primmctruth[4] = {cent_v0, zVtxmc,partMC->Pt(), partMC->Eta()};
+ if(ffillefficiency)
+ {
+ fTHngenprimPidTruth[5]->Fill(primmctruth);//for all primary truth particles(4)
+if (TMath::Abs(pdgtruth)==211 || TMath::Abs(pdgtruth)==321) fTHngenprimPidTruth[3]->Fill(primmctruth);//for primary truth mesons(3)
+if (TMath::Abs(pdgtruth)==2212 || TMath::Abs(pdgtruth)==321) fTHngenprimPidTruth[4]->Fill(primmctruth);//for primary truth kaons+protons(4)
+ if (TMath::Abs(pdgtruth)==211) fTHngenprimPidTruth[0]->Fill(primmctruth);//for pions
+ if (TMath::Abs(pdgtruth)==321) fTHngenprimPidTruth[1]->Fill(primmctruth);//for kaons
+ if(TMath::Abs(pdgtruth)==2212) fTHngenprimPidTruth[2]->Fill(primmctruth);//for protons
+ }
+
+ Float_t effmatrixtruth=1.0;//In Truth MC, no case of efficiency correction so it should be always 1.0
+if(partMC->Pt()>=fminPtAsso || partMC->Pt()<=fmaxPtTrig)//to reduce memory consumption in pool
+ {
+LRCParticlePID* copy6 = new LRCParticlePID(particletypeTruth,partMC->Charge(),partMC->Pt(),partMC->Eta(), partMC->Phi(),effmatrixtruth);
+//copy6->SetUniqueID(eventno * 100000 + TMath::Abs(partMC->GetLabel()));
+ copy6->SetUniqueID(eventno * 100000 + (Int_t)nooftrackstruth);
+ tracksMCtruth->Add(copy6);//************** TObjArray used for truth correlation function calculation
+ }
+ }//MC truth track loop ends
+
+//*********************still in event loop
+
+ if (fSampleType=="pp") cent_v0_truth=nooftrackstruth;
+ else cent_v0_truth=cent_v0;//the notation cent_v0 is reserved for reco/data case
+
+ //now cent_v0_truth should be used for all correlation function calculation
+
+if(nooftrackstruth>0.0 && ffilltrigIDassoIDMCTRUTH)
+ {
+ //Fill Correlations for MC truth particles(same event)
+if(tracksMCtruth && tracksMCtruth->GetEntriesFast()>0)//hadron triggered correlation
+ Fillcorrelation(tracksMCtruth,0,cent_v0_truth,zVtxmc,bSign,fPtOrderMCTruth,kFALSE,kFALSE,"trigIDassoIDMCTRUTH");//mixcase=kFALSE for same event case
+
+//start mixing
+AliEventPool* pool2 = fPoolMgr->GetEventPool(cent_v0_truth, zVtxmc+200);
+if (pool2 && pool2->IsReady())
+ {//start mixing only when pool->IsReady
+if(tracksMCtruth && tracksMCtruth->GetEntriesFast()>0)
+ {//proceed only when no. of trigger particles >0 in current event
+for (Int_t jMix=0; jMix<pool2->GetCurrentNEvents(); jMix++)
+ { //pool event loop start
+ TObjArray* bgTracks6 = pool2->GetEvent(jMix);
+ if(!bgTracks6) continue;
+ Fillcorrelation(tracksMCtruth,bgTracks6,cent_v0_truth,zVtxmc,bSign,fPtOrderMCTruth,kFALSE,kTRUE,"trigIDassoIDMCTRUTH");//mixcase=kTRUE for mixing case
+
+ }// pool event loop ends mixing case
+ }//if(trackstrig && trackstrig->GetEntriesFast()>0) condition ends mixing case
+} //if pool->IsReady() condition ends mixing case
+
+ //still in main event loop
+
+ if(tracksMCtruth){
+if(pool2) pool2->UpdatePool(CloneAndReduceTrackList(tracksMCtruth));//ownership of tracksasso is with pool now, don't delete it
+ }
+ }
+
+ //still in main event loop
+
+if(tracksMCtruth) delete tracksMCtruth;
+
+//now deal with reco tracks
+ //TObjArray* tracksUNID=0;
+ TObjArray* tracksUNID = new TObjArray;
+ tracksUNID->SetOwner(kTRUE);
+
+ //TObjArray* tracksID=0;
+ TObjArray* tracksID = new TObjArray;
+ tracksID->SetOwner(kTRUE);
+
+
+ Float_t bSign1=aod->GetHeader()->GetMagneticField() ;//used for reconstructed track dca cut
+
+ Double_t trackscount=0.0;
+
+// loop over reconstructed tracks
+ for (Int_t itrk = 0; itrk < aod->GetNumberOfTracks(); itrk++)
+{ // reconstructed track loop starts
+ AliAODTrack* track = dynamic_cast<AliAODTrack*>(aod->GetTrack(itrk));
+ if (!track) continue;
+ //get the corresponding MC track at the truth level (doing reco matching)
+ AliAODMCParticle* recomatched = static_cast<AliAODMCParticle*>(fArrayMC->At(TMath::Abs(track->GetLabel())));
+ if(!recomatched) continue;//if a reco track doesn't have corresponding truth track at generated level is a fake track(label==0), ignore it
+
+//remove injected signals
+ if(fInjectedSignals)
+ {
+ AliAODMCParticle* mother = recomatched;
+
+ while (!mother->IsPhysicalPrimary())
+ {// find the primary mother;the first stable mother is searched and checked if it is <= <maxLabel>
+ if (mother->GetMother() < 0)
+ {
+ mother = 0;
+ break;
+ }
+
+ mother =(AliAODMCParticle*) fArrayMC->At(((AliAODMCParticle*)mother)->GetMother());
+ if (!mother)
+ break;
+ }
+ if (!mother)
+ {
+ Printf("WARNING: No mother found for particle %d:", recomatched->GetLabel());
+ continue;
+ }
+ if (mother->GetLabel() >= skipParticlesAbove) continue;//remove injected signals(primaries above <maxLabel>)
+ }//remove injected signals
+
+ if (fRemoveWeakDecays && ((AliAODMCParticle*) recomatched)->IsSecondaryFromWeakDecay()) continue;//remove weak decays
+
+ Bool_t isduplicate2=kFALSE;
+if (fRemoveDuplicates)
+ {
+ for (Int_t j =itrk+1; j < aod->GetNumberOfTracks(); j++)
+ {//2nd loop starts
+ AliAODTrack* track2 = dynamic_cast<AliAODTrack*>(aod->GetTrack(j));
+ if (!track2) continue;
+ AliAODMCParticle* recomatched2 = static_cast<AliAODMCParticle*>(fArrayMC->At(TMath::Abs(track2->GetLabel())));
+if(!recomatched2) continue;
+
+if (track->GetLabel() == track2->GetLabel())
+ {
+isduplicate2=kTRUE;
+ break;
+ }
+ }//2nd loop ends
+ }
+ if(fRemoveDuplicates && isduplicate2) continue;//remove duplicates
+
+ fHistQA[11]->Fill(track->GetTPCNcls());
+ Int_t tracktype=ClassifyTrack(track,trkVtx,bSign1);//dcacut=kFALSE,onlyprimary=kFALSE
+
+ if(tracktype==0) continue;
+ if(tracktype==1)//tracks "not" passed AliAODTrack::kPrimary at reconstructed level & have proper TPC PID response(?)
+{
+ if(!track) continue;//for safety
+ //accepted all(primaries+secondary) reconstructed tracks(pt 0.2 to 10.0,,eta -0.8 to 0.8)
+ trackscount++;
+
+//check for eta , phi holes
+ fEtaSpectrasso->Fill(track->Eta(),track->Pt());
+ fphiSpectraasso->Fill(track->Phi(),track->Pt());
+
+ Int_t particletypeMC=-9999;
+
+//tag all particles as unidentified
+ particletypeMC=unidentified;
+
+ Float_t effmatrix=1.;
+
+// -- Fill THnSparse for efficiency calculation
+ if (fSampleType=="pp") cent_v0=15.0;//integrated over multiplicity(so put any fixed value for each track so that practically means there is only one bin in multiplicityi.e multiplicity intregated out )**************Important
+ //NOTE:: this will be used for fillinfg THnSparse of efficiency & also to get the the track by track eff. factor on the fly(only in case of pp)
+
+ //Clone & Reduce track list(TObjArray) for unidentified particles
+if(track->Pt()>=fminPtAsso || track->Pt()<=fmaxPtTrig)//to reduce memory consumption in pool
+ {
+ if (fapplyTrigefficiency || fapplyAssoefficiency)//get the trackingefficiency x contamination factor for unidentified particles
+ effmatrix=GetTrackbyTrackeffvalue(track,cent_v0,zvtx,particletypeMC);
+ LRCParticlePID* copy = new LRCParticlePID(particletypeMC,track->Charge(),track->Pt(),track->Eta(), track->Phi(),effmatrix);
+ copy->SetUniqueID(eventno * 100000 +(Int_t)trackscount);
+ tracksUNID->Add(copy);//track information Storage for UNID correlation function(tracks that pass the filterbit & kinematic cuts only)
+ }
+ Double_t allrecomatchedpid[4] = {cent_v0, zVtxmc,recomatched->Pt(), recomatched->Eta()};
+if(ffillefficiency) fTHnrecomatchedallPid[5]->Fill(allrecomatchedpid);//for all
+
+ //now start the particle identification process:)
+
+//get the pdg code of the corresponding truth particle
+ Int_t pdgCode = ((AliAODMCParticle*)recomatched)->GetPdgCode();
+
+Float_t dEdx = track->GetTPCsignal();
+ fHistoTPCdEdx->Fill(track->Pt(), dEdx);
+
+ if(HasTOFPID(track))
+{
+Float_t beta = GetBeta(track);
+fHistoTOFbeta->Fill(track->Pt(), beta);
+ }
+
+ //do track identification(nsigma method)
+ particletypeMC=GetParticle(track);//******************************problem is here
+
+//2-d TPCTOF map(for each Pt interval)
+ if(HasTOFPID(track)){
+ fTPCTOFPion3d->Fill(track->Pt(),fnsigmas[SpPion][NSigmaTOF],fnsigmas[SpPion][NSigmaTPC]);
+ fTPCTOFKaon3d->Fill(track->Pt(),fnsigmas[SpKaon][NSigmaTOF],fnsigmas[SpKaon][NSigmaTPC]);
+ fTPCTOFProton3d->Fill(track->Pt(),fnsigmas[SpProton][NSigmaTOF],fnsigmas[SpProton][NSigmaTPC]);
+ }
+ if(particletypeMC==SpUndefined) continue;
+
+ //Pt, Eta , Phi distribution of the reconstructed identified particles
+if(ffillhistQAReco)
+ {
+if (particletypeMC==SpPion)
+ {
+ fPionPt->Fill(track->Pt());
+ fPionEta->Fill(track->Eta());
+ fPionPhi->Fill(track->Phi());
+ }
+if (particletypeMC==SpKaon)
+ {
+ fKaonPt->Fill(track->Pt());
+ fKaonEta->Fill(track->Eta());
+ fKaonPhi->Fill(track->Phi());
+ }
+if (particletypeMC==SpProton)
+ {
+ fProtonPt->Fill(track->Pt());
+ fProtonEta->Fill(track->Eta());
+ fProtonPhi->Fill(track->Phi());
+ }
+ }
+
+ //for misidentification fraction calculation(do it with tuneonPID)
+ if(particletypeMC==SpPion )
+ {
+ if(TMath::Abs(pdgCode)==211) fPioncont->Fill(1.,track->Pt());
+ if(TMath::Abs(pdgCode)==321) fPioncont->Fill(3.,track->Pt());
+ if(TMath::Abs(pdgCode)==2212) fPioncont->Fill(5.,track->Pt());
+ if(TMath::Abs(pdgCode)!=211 && TMath::Abs(pdgCode)!=321 && TMath::Abs(pdgCode)!=2212) fPioncont->Fill(7.,track->Pt());
+ }
+if(particletypeMC==SpKaon )
+ {
+ if(TMath::Abs(pdgCode)==211) fKaoncont->Fill(1.,track->Pt());
+ if(TMath::Abs(pdgCode)==321) fKaoncont->Fill(3.,track->Pt());
+ if(TMath::Abs(pdgCode)==2212) fKaoncont->Fill(5.,track->Pt());
+ if(TMath::Abs(pdgCode)!=211 && TMath::Abs(pdgCode)!=321 && TMath::Abs(pdgCode)!=2212) fKaoncont->Fill(7.,track->Pt());
+ }
+ if(particletypeMC==SpProton )
+ {
+ if(TMath::Abs(pdgCode)==211) fProtoncont->Fill(1.,track->Pt());
+ if(TMath::Abs(pdgCode)==321) fProtoncont->Fill(3.,track->Pt());
+ if(TMath::Abs(pdgCode)==2212) fProtoncont->Fill(5.,track->Pt());
+ if(TMath::Abs(pdgCode)!=211 && TMath::Abs(pdgCode)!=321 && TMath::Abs(pdgCode)!=2212) fProtoncont->Fill(7.,track->Pt());
+ }
+
+ //fill tracking efficiency
+ if(ffillefficiency)
+ {
+ if(particletypeMC==SpPion || particletypeMC==SpKaon)
+ {
+if(TMath::Abs(pdgCode)==211 || TMath::Abs(pdgCode)==321) fTHnrecomatchedallPid[3]->Fill(allrecomatchedpid);//for mesons
+ }
+ if(particletypeMC==SpKaon || particletypeMC==SpProton)
+ {
+if(TMath::Abs(pdgCode)==321 || TMath::Abs(pdgCode)==2212) fTHnrecomatchedallPid[4]->Fill(allrecomatchedpid);//for kaons+protons
+ }
+ if(particletypeMC==SpPion && TMath::Abs(pdgCode)==211) fTHnrecomatchedallPid[0]->Fill(allrecomatchedpid);//for pions
+ if(particletypeMC==SpKaon && TMath::Abs(pdgCode)==321) fTHnrecomatchedallPid[1]->Fill(allrecomatchedpid);//for kaons
+ if(particletypeMC==SpProton && TMath::Abs(pdgCode)==2212) fTHnrecomatchedallPid[2]->Fill(allrecomatchedpid);//for protons
+ }
+
+if(track->Pt()>=fminPtAsso || track->Pt()<=fmaxPtTrig)//to reduce memory consumption in pool
+ {
+if (fapplyTrigefficiency || fapplyAssoefficiency)
+ effmatrix=GetTrackbyTrackeffvalue(track,cent_v0,zvtx,particletypeMC);//get the tracking eff x TOF matching eff x PID eff x contamination factor for identified particles
+ LRCParticlePID* copy1 = new LRCParticlePID(particletypeMC,track->Charge(),track->Pt(),track->Eta(), track->Phi(),effmatrix);
+ copy1->SetUniqueID(eventno * 100000 + (Int_t)trackscount);
+ tracksID->Add(copy1);
+ }
+ }// if(tracktype==1) condition structure ands
+
+}//reco track loop ends
+
+ //*************************************************************still in event loop
+
+//same event delta-eta-deltaphi plot
+if(fSampleType=="pp") cent_v0=trackscount;//multiplicity
+
+if(trackscount>0.0)
+ {
+//fill the centrality/multiplicity distribution of the selected events
+ fhistcentrality->Fill(cent_v0);//*********************************WARNING::binning of cent_v0 is different for pp and pPb/PbPb case
+
+ if (fSampleType=="pPb" || fSampleType=="PbPb") fCentralityCorrelation->Fill(cent_v0, trackscount);//only with unidentified tracks(i.e before PID selection);;;;;can be used to remove centrality outliers??????
+
+//count selected events having centrality betn 0-100%
+ fEventCounter->Fill(6);
+
+ //same event delte-eta, delta-phi plot
+if(tracksUNID && tracksUNID->GetEntriesFast()>0)//hadron triggered correlation
+ {//same event calculation starts
+ if(ffilltrigassoUNID) Fillcorrelation(tracksUNID,0,cent_v0,zvtx,bSign,kTRUE,kTRUE,kFALSE,"trigassoUNID");//mixcase=kFALSE (hadron-hadron correlation)
+ if(tracksID && tracksID->GetEntriesFast()>0 && ffilltrigUNIDassoID) Fillcorrelation(tracksUNID,tracksID,cent_v0,zvtx,bSign,kFALSE,kTRUE,kFALSE,"trigUNIDassoID");//mixcase=kFALSE (hadron-ID correlation)
+ }
+
+if(tracksID && tracksID->GetEntriesFast()>0)//ID triggered correlation
+ {//same event calculation starts
+ if(tracksUNID && tracksUNID->GetEntriesFast()>0 && ffilltrigIDassoUNID) Fillcorrelation(tracksID,tracksUNID,cent_v0,zvtx,bSign,kFALSE,kTRUE,kFALSE,"trigIDassoUNID");//mixcase=kFALSE (ID-hadron correlation)
+ if(ffilltrigIDassoID) Fillcorrelation(tracksID,0,cent_v0,zvtx,bSign,kFALSE,kTRUE,kFALSE,"trigIDassoID");//mixcase=kFALSE (ID-ID correlation)
+ }
+
+//still in main event loop
+//start mixing
+ if(ffilltrigassoUNID || ffilltrigIDassoUNID){//mixing with unidentified particles
+AliEventPool* pool = fPoolMgr->GetEventPool(cent_v0, zvtx);//In the pool there is tracksUNID(i.e associateds are unidentified)
+if (pool && pool->IsReady())
+ {//start mixing only when pool->IsReady
+ for (Int_t jMix=0; jMix<pool->GetCurrentNEvents(); jMix++)
+ { //pool event loop start
+ TObjArray* bgTracks = pool->GetEvent(jMix);
+ if(!bgTracks) continue;
+ if(ffilltrigassoUNID && tracksUNID && tracksUNID->GetEntriesFast()>0)//*******************************hadron trggered mixing
+ Fillcorrelation(tracksUNID,bgTracks,cent_v0,zvtx,bSign,kTRUE,kTRUE,kTRUE,"trigassoUNID");//mixcase=kTRUE
+ if(ffilltrigIDassoUNID && tracksID && tracksID->GetEntriesFast()>0)//***********************************ID trggered mixing
+ Fillcorrelation(tracksID,bgTracks,cent_v0,zvtx,bSign,kFALSE,kTRUE,kTRUE,"trigIDassoUNID");//mixcase=kTRUE
+ }// pool event loop ends mixing case
+
+} //if pool->IsReady() condition ends mixing case
+ if(tracksUNID) {
+if(pool)
+ pool->UpdatePool(CloneAndReduceTrackList(tracksUNID));
+ }
+ }//mixing with unidentified particles
+
+ if(ffilltrigUNIDassoID || ffilltrigIDassoID){//mixing with identified particles
+AliEventPool* pool1 = fPoolMgr->GetEventPool(cent_v0, zvtx+100);//In the pool1 there is tracksID(i.e associateds are identified)
+if (pool1 && pool1->IsReady())
+ {//start mixing only when pool->IsReady
+for (Int_t jMix=0; jMix<pool1->GetCurrentNEvents(); jMix++)
+ { //pool event loop start
+ TObjArray* bgTracks2 = pool1->GetEvent(jMix);
+ if(!bgTracks2) continue;
+if(ffilltrigUNIDassoID && tracksUNID && tracksUNID->GetEntriesFast()>0)
+ Fillcorrelation(tracksUNID,bgTracks2,cent_v0,zvtx,bSign,kFALSE,kTRUE,kTRUE,"trigUNIDassoID");//mixcase=kTRUE
+if(ffilltrigIDassoID && tracksID && tracksID->GetEntriesFast()>0)
+ Fillcorrelation(tracksID,bgTracks2,cent_v0,zvtx,bSign,kFALSE,kTRUE,kTRUE,"trigIDassoID");//mixcase=kTRUE
+
+ }// pool event loop ends mixing case
+} //if pool1->IsReady() condition ends mixing case
+
+if(tracksID) {
+if(pool1)
+ pool1->UpdatePool(CloneAndReduceTrackList(tracksID));//ownership of tracksasso is with pool now, don't delete it(tracksUNID is with pool)
+ }
+ }//mixing with identified particles
+
+ //no. of events analyzed
+fEventCounter->Fill(7);
+ }
+
+if(tracksUNID) delete tracksUNID;
+
+if(tracksID) delete tracksID;
+
+
+PostData(1, fOutput);
+
+}
+//________________________________________________________________________
+void AliTwoParticlePIDCorr::doAODevent()
+{
+
+ //get AOD
+ AliVEvent *event = InputEvent();
+ if (!event) { Printf("ERROR: Could not retrieve event"); return; }
+ AliAODEvent* aod = dynamic_cast<AliAODEvent*>(event);
+ if (!aod) {
+ AliError("Cannot get the AOD event");
+ return;
+ }
+
+// count all events
+ fEventCounter->Fill(1);
+
+ // get centrality object and check quality
+ Double_t cent_v0=0;
+
+ if(fSampleType=="pPb" || fSampleType=="PbPb")
+ {
+ AliCentrality *centrality=0;
+ if(aod)
+ centrality = aod->GetHeader()->GetCentralityP();
+ // if (centrality->GetQuality() != 0) return ;
+
+ if(centrality)
+ {
+ cent_v0 = centrality->GetCentralityPercentile(fCentralityMethod);
+ }
+ else
+ {
+ cent_v0= -1;
+ }
+ }
+
+ Float_t bSign = (aod->GetMagneticField() > 0) ? 1 : -1;//for two track efficiency cut in correlation function calculation
+ Float_t bSign1=aod->GetHeader()->GetMagneticField() ;//for dca cut in ClassifyTrack(), i.e in track loop
+
+// Pileup selection ************************************************
+// if(frejectPileUp) //applicable only for TPC only tracks,not for hybrid tracks?.
+ // {
+ //if (fAnalysisUtils && fAnalysisUtils->IsPileUpEvent(aod)) return;
+ // }
+
+ if (!fPID) return;//this should be available with each event even if we don't do PID selection
+
+ //vertex selection(is it fine for PP?)
+ if ( fSampleType=="pPb"){
+ trkVtx = aod->GetPrimaryVertex();
+ if (!trkVtx || trkVtx->GetNContributors()<=0) return;
+ TString vtxTtl = trkVtx->GetTitle();
+ if (!vtxTtl.Contains("VertexerTracks")) return;
+ zvtx = trkVtx->GetZ();
+ const AliAODVertex* spdVtx = aod->GetPrimaryVertexSPD();
+ if (!spdVtx || spdVtx->GetNContributors()<=0) return;
+ TString vtxTyp = spdVtx->GetTitle();
+ Double_t cov[6]={0};
+ spdVtx->GetCovarianceMatrix(cov);
+ Double_t zRes = TMath::Sqrt(cov[5]);
+ if (vtxTyp.Contains("vertexer:Z") && (zRes>0.25)) return;
+ if (TMath::Abs(spdVtx->GetZ() - trkVtx->GetZ())>0.5) return;
+ }
+ else if(fSampleType=="PbPb" || fSampleType=="pp") {//for pp and pb-pb case
+ Int_t nVertex = aod->GetNumberOfVertices();
+ if( nVertex > 0 ) {
+ trkVtx = aod->GetPrimaryVertex();
+ Int_t nTracksPrim = trkVtx->GetNContributors();
+ zvtx = trkVtx->GetZ();
+ //if (fDebug > 1)AliInfo(Form(" Vertex in = %f with %d particles by %s data ...",zVertex,nTracksPrim,vertex->GetName()));
+ // Reject TPC only vertex
+ TString name(trkVtx->GetName());
+ if (name.CompareTo("PrimaryVertex") && name.CompareTo("SPDVertex"))return;
+
+ // Select a quality vertex by number of tracks?
+ if( nTracksPrim < fnTracksVertex ) {
+ //if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
+ return ;
+ }
+ // TODO remove vertexer Z events with dispersion > 0.02: Doesn't work for AOD at present
+ //if (strcmp(vertex->GetTitle(), "AliVertexerZ") == 0 && vertex->GetDispersion() > 0.02)
+ // return kFALSE;
+ // if (fDebug > 1) AliInfo(" Primary-vertex Selection: event ACCEPTED...");
+ }
+ else return;
+
+ }
+ else return;//as there is no proper sample type
+
+
+ fHistQA[0]->Fill((trkVtx->GetX()));fHistQA[1]->Fill((trkVtx->GetY()));fHistQA[2]->Fill((trkVtx->GetZ())); //for trkVtx only before vertex cut |zvtx|<10 cm
+
+//count events having a proper vertex
+ fEventCounter->Fill(2);
+
+ if (TMath::Abs(zvtx) > fzvtxcut) return;
+
+//count events after vertex cut
+ fEventCounter->Fill(5);
+
+
+ //if(!fAnalysisUtils->IsVertexSelected2013pA(aod)) return;
+
+ fHistQA[3]->Fill((trkVtx->GetX()));fHistQA[4]->Fill((trkVtx->GetY()));fHistQA[5]->Fill((trkVtx->GetZ()));//after vertex cut,for trkVtx only
+
+
+ if(!aod) return; //for safety
+
+if(fSampleType=="pPb" || fSampleType=="PbPb") if (cent_v0 < 0) return;//for pp case it is the multiplicity
+
+ TObjArray* tracksUNID= new TObjArray;//track info before doing PID
+ tracksUNID->SetOwner(kTRUE); // IMPORTANT!
+
+ TObjArray* tracksID= new TObjArray;//only pions, kaons,protonsi.e. after doing the PID selection
+ tracksID->SetOwner(kTRUE); // IMPORTANT!
+
+ Double_t trackscount=0.0;//in case of pp this will give the multiplicity after the track loop(only for unidentified particles that pass the filterbit and kinematic cuts)
+
+ eventno++;
+
+ for (Int_t itrk = 0; itrk < aod->GetNumberOfTracks(); itrk++)
+{ //track loop starts for TObjArray(containing track and event information) filling; used for correlation function calculation
+ AliAODTrack* track = dynamic_cast<AliAODTrack*>(aod->GetTrack(itrk));
+ if (!track) continue;
+ fHistQA[11]->Fill(track->GetTPCNcls());
+ Int_t particletype=-9999;//required for PID filling
+ Int_t tracktype=ClassifyTrack(track,trkVtx,bSign1);//dcacut=kFALSE,onlyprimary=kFALSE
+ if(tracktype!=1) continue;
+
+ if(!track) continue;//for safety
+
+//check for eta , phi holes
+ fEtaSpectrasso->Fill(track->Eta(),track->Pt());
+ fphiSpectraasso->Fill(track->Phi(),track->Pt());
+
+ trackscount++;
+
+ //if no applyefficiency , set the eff factor=1.0
+ Float_t effmatrix=1.0;
+
+ //tag all particles as unidentified that passed filterbit & kinematic cuts
+ particletype=unidentified;
+
+
+ if (fSampleType=="pp") cent_v0=15.0;//integrated over multiplicity [i.e each track has multiplicity 15.0](so put any fixed value for each track so that practically means there is only one bin in multiplicityi.e multiplicity intregated out )**************Important for efficiency related issues
+
+
+ //to reduce memory consumption in pool
+ if(track->Pt()>=fminPtAsso || track->Pt()<=fmaxPtTrig)
+ {
+ //Clone & Reduce track list(TObjArray) for unidentified particles
+ if (fapplyTrigefficiency || fapplyAssoefficiency)//get the trackingefficiency x contamination factor for unidentified particles
+ effmatrix=GetTrackbyTrackeffvalue(track,cent_v0,zvtx,particletype);
+ LRCParticlePID* copy = new LRCParticlePID(particletype,track->Charge(),track->Pt(),track->Eta(), track->Phi(),effmatrix);
+ copy->SetUniqueID(eventno * 100000 + (Int_t)trackscount);
+ tracksUNID->Add(copy);//track information Storage for UNID correlation function(tracks that pass the filterbit & kinematic cuts only)
+ }
+
+//now start the particle identificaion process:)
+
+//track passing filterbit 768 have proper TPC response,or need to be checked explicitly before doing PID????
+
+ Float_t dEdx = track->GetTPCsignal();
+ fHistoTPCdEdx->Fill(track->Pt(), dEdx);
+
+ //fill beta vs Pt plots only for tracks having proper TOF response(much less tracks compared to the no. that pass the filterbit & kinematic cuts)
+ if(HasTOFPID(track))
+{
+ Float_t beta = GetBeta(track);
+ fHistoTOFbeta->Fill(track->Pt(), beta);
+ }
+
+
+//track identification(using nsigma method)
+ particletype=GetParticle(track);//*******************************change may be required(It should return only pion,kaon, proton and Spundefined; NOT unidentifed***************be careful)
+
+
+//2-d TPCTOF map(for each Pt interval)
+ if(HasTOFPID(track)){
+ fTPCTOFPion3d->Fill(track->Pt(),fnsigmas[SpPion][NSigmaTOF],fnsigmas[SpPion][NSigmaTPC]);
+ fTPCTOFKaon3d->Fill(track->Pt(),fnsigmas[SpKaon][NSigmaTOF],fnsigmas[SpKaon][NSigmaTPC]);
+ fTPCTOFProton3d->Fill(track->Pt(),fnsigmas[SpProton][NSigmaTOF],fnsigmas[SpProton][NSigmaTPC]);
+ }
+
+//ignore the Spundefined particles as they also contain pion, kaon, proton outside the nsigma cut(also if tracks don't have proper TOF PID in a certain Pt interval) & these tracks are actually counted when we do the the efficiency correction, so considering them as unidentified particles & doing the efficiency correction(i.e defining unidentified=pion+Kaon+proton+SpUndefined is right only without efficiency correction) for them will be two times wrong!!!!!
+ if (particletype==SpUndefined) continue;//this condition creating a modulated structure in delphi projection in mixed event case(only when we are dealing with identified particles i.e. tracksID)!!!!!!!!!!!
+
+ //Pt, Eta , Phi distribution of the reconstructed identified particles
+if(ffillhistQAReco)
+ {
+if (particletype==SpPion)
+ {
+ fPionPt->Fill(track->Pt());
+ fPionEta->Fill(track->Eta());
+ fPionPhi->Fill(track->Phi());
+ }
+if (particletype==SpKaon)
+ {
+ fKaonPt->Fill(track->Pt());
+ fKaonEta->Fill(track->Eta());
+ fKaonPhi->Fill(track->Phi());
+ }
+if (particletype==SpProton)
+ {
+ fProtonPt->Fill(track->Pt());
+ fProtonEta->Fill(track->Eta());
+ fProtonPhi->Fill(track->Phi());
+ }
+ }
+
+if(track->Pt()>=fminPtAsso || track->Pt()<=fmaxPtTrig)
+ {
+if (fapplyTrigefficiency || fapplyAssoefficiency)
+ effmatrix=GetTrackbyTrackeffvalue(track,cent_v0,zvtx,particletype);//get the tracking eff x TOF matching eff x PID eff x contamination factor for identified particles; Bool_t mesoneffrequired=kFALSE
+ LRCParticlePID* copy1 = new LRCParticlePID(particletype,track->Charge(),track->Pt(),track->Eta(), track->Phi(),effmatrix);
+ copy1->SetUniqueID(eventno * 100000 + (Int_t)trackscount);
+ tracksID->Add(copy1);
+ }
+} //track loop ends but still in event loop
+
+if(trackscount<1.0){
+ if(tracksUNID) delete tracksUNID;
+ if(tracksID) delete tracksID;
+ return;
+ }
+
+// cout<<fminPtAsso<<"***"<<fmaxPtAsso<<"*********************"<<fminPtTrig<<"***"<<fmaxPtTrig<<"*****************"<<kTrackVariablesPair<<endl;
+if(fSampleType=="pp") cent_v0=trackscount;//multiplicity
+
+//fill the centrality/multiplicity distribution of the selected events
+ fhistcentrality->Fill(cent_v0);//*********************************WARNING::binning of cent_v0 is different for pp and pPb/PbPb case
+
+if(fSampleType=="pPb" || fSampleType=="PbPb") fCentralityCorrelation->Fill(cent_v0, trackscount);//only with unidentified tracks(i.e before PID selection);;;;;can be used to remove centrality outliers??????
+
+//count selected events having centrality betn 0-100%
+ fEventCounter->Fill(6);
+
+//same event delta-eta-deltaphi plot
+
+if(tracksUNID && tracksUNID->GetEntriesFast()>0)//hadron triggered correlation
+ {//same event calculation starts
+ if(ffilltrigassoUNID) Fillcorrelation(tracksUNID,0,cent_v0,zvtx,bSign,kTRUE,kTRUE,kFALSE,"trigassoUNID");//mixcase=kFALSE (hadron-hadron correlation)
+ if(tracksID && tracksID->GetEntriesFast()>0 && ffilltrigUNIDassoID) Fillcorrelation(tracksUNID,tracksID,cent_v0,zvtx,bSign,kFALSE,kTRUE,kFALSE,"trigUNIDassoID");//mixcase=kFALSE (hadron-ID correlation)
+ }
+
+if(tracksID && tracksID->GetEntriesFast()>0)//ID triggered correlation
+ {//same event calculation starts
+ if(tracksUNID && tracksUNID->GetEntriesFast()>0 && ffilltrigIDassoUNID) Fillcorrelation(tracksID,tracksUNID,cent_v0,zvtx,bSign,kFALSE,kTRUE,kFALSE,"trigIDassoUNID");//mixcase=kFALSE (ID-hadron correlation)
+ if(ffilltrigIDassoID) Fillcorrelation(tracksID,0,cent_v0,zvtx,bSign,kFALSE,kTRUE,kFALSE,"trigIDassoID");//mixcase=kFALSE (ID-ID correlation)
+ }
+
+//still in main event loop
+
+
+//start mixing
+ if(ffilltrigassoUNID || ffilltrigIDassoUNID){//mixing with unidentified particles
+AliEventPool* pool = fPoolMgr->GetEventPool(cent_v0, zvtx);//In the pool there is tracksUNID(i.e associateds are unidentified)
+if (pool && pool->IsReady())
+ {//start mixing only when pool->IsReady
+ for (Int_t jMix=0; jMix<pool->GetCurrentNEvents(); jMix++)
+ { //pool event loop start
+ TObjArray* bgTracks = pool->GetEvent(jMix);
+ if(!bgTracks) continue;
+ if(ffilltrigassoUNID && tracksUNID && tracksUNID->GetEntriesFast()>0)//*******************************hadron trggered mixing
+ Fillcorrelation(tracksUNID,bgTracks,cent_v0,zvtx,bSign,kTRUE,kTRUE,kTRUE,"trigassoUNID");//mixcase=kTRUE
+ if(ffilltrigIDassoUNID && tracksID && tracksID->GetEntriesFast()>0)//***********************************ID trggered mixing
+ Fillcorrelation(tracksID,bgTracks,cent_v0,zvtx,bSign,kFALSE,kTRUE,kTRUE,"trigIDassoUNID");//mixcase=kTRUE
+ }// pool event loop ends mixing case
+
+} //if pool->IsReady() condition ends mixing case
+ if(tracksUNID) {
+if(pool)
+ pool->UpdatePool(CloneAndReduceTrackList(tracksUNID));
+ }
+ }//mixing with unidentified particles
+
+
+ if(ffilltrigUNIDassoID || ffilltrigIDassoID){//mixing with identified particles
+AliEventPool* pool1 = fPoolMgr->GetEventPool(cent_v0, zvtx+100);//In the pool1 there is tracksID(i.e associateds are identified)
+if (pool1 && pool1->IsReady())
+ {//start mixing only when pool->IsReady
+for (Int_t jMix=0; jMix<pool1->GetCurrentNEvents(); jMix++)
+ { //pool event loop start
+ TObjArray* bgTracks2 = pool1->GetEvent(jMix);
+ if(!bgTracks2) continue;
+if(ffilltrigUNIDassoID && tracksUNID && tracksUNID->GetEntriesFast()>0)
+ Fillcorrelation(tracksUNID,bgTracks2,cent_v0,zvtx,bSign,kFALSE,kTRUE,kTRUE,"trigUNIDassoID");//mixcase=kTRUE
+if(ffilltrigIDassoID && tracksID && tracksID->GetEntriesFast()>0)
+ Fillcorrelation(tracksID,bgTracks2,cent_v0,zvtx,bSign,kFALSE,kTRUE,kTRUE,"trigIDassoID");//mixcase=kTRUE
+
+ }// pool event loop ends mixing case
+} //if pool1->IsReady() condition ends mixing case
+
+if(tracksID) {
+if(pool1)
+ pool1->UpdatePool(CloneAndReduceTrackList(tracksID));//ownership of tracksasso is with pool now, don't delete it(tracksUNID is with pool)
+ }
+ }//mixing with identified particles
+
+
+ //no. of events analyzed
+fEventCounter->Fill(7);
+
+//still in main event loop
+
+
+if(tracksUNID) delete tracksUNID;
+
+if(tracksID) delete tracksID;
+
+
+PostData(1, fOutput);
+
+} // *************************event loop ends******************************************//_______________________________________________________________________
+TObjArray* AliTwoParticlePIDCorr::CloneAndReduceTrackList(TObjArray* tracks)
+{
+ // clones a track list by using AliDPhiBasicParticle which uses much less memory (used for event mixing)
+
+ TObjArray* tracksClone = new TObjArray;
+ tracksClone->SetOwner(kTRUE);
+
+ for (Int_t i=0; i<tracks->GetEntriesFast(); i++)
+ {
+ LRCParticlePID* particle = (LRCParticlePID*) tracks->UncheckedAt(i);
+ LRCParticlePID* copy100 = new LRCParticlePID(particle->getparticle(),particle->Charge(), particle->Pt(),particle->Eta(), particle->Phi(), particle->geteffcorrectionval());
+ copy100->SetUniqueID(particle->GetUniqueID());
+ tracksClone->Add(copy100);
+ }
+
+ return tracksClone;
+}
+
+//--------------------------------------------------------------------------------
+void AliTwoParticlePIDCorr::Fillcorrelation(TObjArray *trackstrig,TObjArray *tracksasso,Double_t cent,Float_t vtx,Float_t bSign,Bool_t fPtOrder,Bool_t twoTrackEfficiencyCut,Bool_t mixcase,TString fillup)
+{
+
+ //before calling this function check that either trackstrig & tracksasso are available
+
+ // Eta() is extremely time consuming, therefore cache it for the inner loop here:
+ TObjArray* input = (tracksasso) ? tracksasso : trackstrig;
+ TArrayF eta(input->GetEntriesFast());
+ for (Int_t i=0; i<input->GetEntriesFast(); i++)
+ eta[i] = ((LRCParticlePID*) input->UncheckedAt(i))->Eta();
+
+ //if(trackstrig)
+ Int_t jmax=trackstrig->GetEntriesFast();
+ if(tracksasso)
+ jmax=tracksasso->GetEntriesFast();
+
+// identify K, Lambda candidates and flag those particles
+ // a TObject bit is used for this
+const UInt_t kResonanceDaughterFlag = 1 << 14;
+ if (fRejectResonanceDaughters > 0)
+ {
+ Double_t resonanceMass = -1;
+ Double_t massDaughter1 = -1;
+ Double_t massDaughter2 = -1;
+ const Double_t interval = 0.02;
+ switch (fRejectResonanceDaughters)
+ {
+ case 1: resonanceMass = 0.9; massDaughter1 = 0.1396; massDaughter2 = 0.9383; break; // method test
+ case 2: resonanceMass = 0.4976; massDaughter1 = 0.1396; massDaughter2 = massDaughter1; break; // k0
+ case 3: resonanceMass = 1.115; massDaughter1 = 0.1396; massDaughter2 = 0.9383; break; // lambda
+ default: AliFatal(Form("Invalid setting %d", fRejectResonanceDaughters));
+ }
+
+for (Int_t i=0; i<trackstrig->GetEntriesFast(); i++)
+ trackstrig->UncheckedAt(i)->ResetBit(kResonanceDaughterFlag);
+ for (Int_t i=0; tracksasso->GetEntriesFast(); i++)
+ tracksasso->UncheckedAt(i)->ResetBit(kResonanceDaughterFlag);
+
+ for (Int_t i=0; i<trackstrig->GetEntriesFast(); i++)
+ {
+ LRCParticlePID *trig=(LRCParticlePID*)(trackstrig->UncheckedAt(i));
+for (Int_t j=0; tracksasso->GetEntriesFast(); j++)
+ {
+ LRCParticlePID *asso=(LRCParticlePID*)(tracksasso->UncheckedAt(j));
+
+ // check if both particles point to the same element (does not occur for mixed events, but if subsets are mixed within the same event)
+if (trig->IsEqual(asso)) continue;
+
+if (trig->Charge() * asso->Charge() > 0) continue;
+
+ Float_t mass = GetInvMassSquaredCheap(trig->Pt(), trig->Eta(), trig->Phi(), asso->Pt(), asso->Eta(), asso->Phi(), massDaughter1, massDaughter2);
+
+if (TMath::Abs(mass - resonanceMass*resonanceMass) < interval*5)
+ {
+ mass = GetInvMassSquared(trig->Pt(), trig->Eta(), trig->Phi(), asso->Pt(), asso->Eta(), asso->Phi(), massDaughter1, massDaughter2);
+
+ if (mass > (resonanceMass-interval)*(resonanceMass-interval) && mass < (resonanceMass+interval)*(resonanceMass+interval))
+ {
+ trig->SetBit(kResonanceDaughterFlag);
+ asso->SetBit(kResonanceDaughterFlag);
+
+// Printf("Flagged %d %d %f", i, j, TMath::Sqrt(mass));
+ }
+ }
+ }
+ }
+ }
+
+ //two particle correlation filling
+
+for(Int_t i=0;i<trackstrig->GetEntriesFast();i++)
+ { //trigger loop starts
+ LRCParticlePID *trig=(LRCParticlePID*)(trackstrig->UncheckedAt(i));
+ if(!trig) continue;
+ Float_t trigpt=trig->Pt();
+ //to avoid overflow qnd underflow
+ if(trigpt<fminPtTrig || trigpt>fmaxPtTrig) continue;
+ Int_t particlepidtrig=trig->getparticle();
+ if(fTriggerSpeciesSelection){ if (particlepidtrig!=fTriggerSpecies) continue;}
+
+ Float_t trigeta=trig->Eta();
+
+ // some optimization
+ if (fTriggerRestrictEta > 0 && TMath::Abs(trigeta) > fTriggerRestrictEta)
+ continue;
+
+if (fOnlyOneEtaSide != 0)
+ {
+ if (fOnlyOneEtaSide * trigeta < 0)
+ continue;
+ }
+ if (fTriggerSelectCharge != 0)
+ if (trig->Charge() * fTriggerSelectCharge < 0)
+ continue;
+
+ if (fRejectResonanceDaughters > 0)
+ if (trig->TestBit(kResonanceDaughterFlag)) continue;
+
+ Float_t trigphi=trig->Phi();
+ Float_t trackefftrig=1.0;
+ if(fapplyTrigefficiency) trackefftrig=trig->geteffcorrectionval();
+ // cout<<"*******************trackefftrig="<<trackefftrig<<endl;
+ Double_t* trigval;
+ Int_t dim=3;
+ if(fcontainPIDtrig) dim=4;
+ trigval= new Double_t[dim];
+ trigval[0] = cent;
+ trigval[1] = vtx;
+ trigval[2] = trigpt;
+if(fcontainPIDtrig) trigval[3] = particlepidtrig;
+
+ //filling only for same event case(mixcase=kFALSE)
+
+ //trigger particle counting only when mixcase=kFALSE i.e for same event correlation function calculation
+if(mixcase==kFALSE)
+ {
+ if(ffilltrigassoUNID==kTRUE && ffilltrigUNIDassoID==kTRUE){
+ if(fillup=="trigassoUNID" ) fTHnTrigcount->Fill(trigval,1.0/trackefftrig);
+ }
+ if(ffilltrigassoUNID==kTRUE && ffilltrigUNIDassoID==kFALSE){
+ if(fillup=="trigassoUNID" ) fTHnTrigcount->Fill(trigval,1.0/trackefftrig);
+ }
+if(ffilltrigassoUNID==kFALSE && ffilltrigUNIDassoID==kTRUE){
+ if(fillup=="trigUNIDassoID") fTHnTrigcount->Fill(trigval,1.0/trackefftrig);
+ }
+ //ensure that trigIDassoID , trigassoUNID, trigIDassoUNID & trigUNIDassoID case FillCorrelation called only once in the event loop for same event correlation function calculation, otherwise there will be multiple counting of pion, kaon,proton,unidentified
+if(ffilltrigIDassoUNID==kTRUE && ffilltrigIDassoID==kTRUE){
+ if(fillup=="trigIDassoID") fTHnTrigcount->Fill(trigval,1.0/trackefftrig);
+ }
+ if(ffilltrigIDassoUNID==kTRUE && ffilltrigIDassoID==kFALSE){
+ if(fillup=="trigIDassoUNID" ) fTHnTrigcount->Fill(trigval,1.0/trackefftrig);
+ }
+if(ffilltrigIDassoUNID==kFALSE && ffilltrigIDassoID==kTRUE){
+ if(fillup=="trigIDassoID") fTHnTrigcount->Fill(trigval,1.0/trackefftrig);
+ }
+
+ if(fillup=="trigIDassoIDMCTRUTH") fTHnTrigcountMCTruthPrim->Fill(trigval,1.0/trackefftrig); //In truth MC case "Unidentified" means any particle other than pion,kaon or proton and no efficiency correction(default value 1.0)************************be careful!!!!
+ }
+
+ //asso loop starts within trigger loop
+ for(Int_t j=0;j<jmax;j++)
+ {
+ LRCParticlePID *asso=0;
+ if(!tracksasso)
+ asso=(LRCParticlePID*)(trackstrig->UncheckedAt(j));
+ else
+ asso=(LRCParticlePID*)(tracksasso->UncheckedAt(j));
+
+ if(!asso) continue;
+
+ //to avoid overflow qnd underflow
+ if(asso->Pt()<fminPtAsso || asso->Pt()>fmaxPtAsso) continue;//***********************Important
+
+ if(fmaxPtAsso==fminPtTrig) {if(asso->Pt()==fminPtTrig) continue;}//******************Think about it!
+
+ if(!tracksasso && i==j) continue;
+
+ // check if both particles point to the same element (does not occur for mixed events, but if subsets are mixed within the same event)
+ // if (tracksasso && trig->IsEqual(asso)) continue;
+
+ if (tracksasso && (trig->GetUniqueID()==asso->GetUniqueID())) continue;
+
+ if (fPtOrder)
+ if (asso->Pt() >= trig->Pt()) continue;
+
+ Int_t particlepidasso=asso->getparticle();
+ if(fAssociatedSpeciesSelection){ if (particlepidasso!=fAssociatedSpecies) continue;}
+
+
+if (fAssociatedSelectCharge != 0)
+if (asso->Charge() * fAssociatedSelectCharge < 0) continue;
+
+ if (fSelectCharge > 0)
+ {
+ // skip like sign
+ if (fSelectCharge == 1 && asso->Charge() * trig->Charge() > 0)
+ continue;
+
+ // skip unlike sign
+ if (fSelectCharge == 2 && asso->Charge() * trig->Charge() < 0)
+ continue;
+ }
+
+if (fEtaOrdering)
+ {
+ if (trigeta < 0 && asso->Eta() < trigeta)
+ continue;
+ if (trigeta > 0 && asso->Eta() > trigeta)
+ continue;
+ }
+
+if (fRejectResonanceDaughters > 0)
+ if (asso->TestBit(kResonanceDaughterFlag))
+ {
+// Printf("Skipped j=%d", j);
+ continue;
+ }
+
+ // conversions
+ if (fCutConversions && asso->Charge() * trig->Charge() < 0)
+ {
+ Float_t mass = GetInvMassSquaredCheap(trig->Pt(), trigeta, trig->Phi(), asso->Pt(),eta[j], asso->Phi(), 0.510e-3, 0.510e-3);
+
+ if (mass < 0.1)
+ {
+ mass = GetInvMassSquared(trig->Pt(), trigeta, trig->Phi(), asso->Pt(), eta[j], asso->Phi(), 0.510e-3, 0.510e-3);
+
+ //fControlConvResoncances->Fill(0.0, mass);
+
+ if (mass < 0.04*0.04)
+ continue;
+ }
+ }
+
+ // K0s
+ if (fCutResonances && asso->Charge() * trig->Charge() < 0)
+ {
+ Float_t mass = GetInvMassSquaredCheap(trig->Pt(), trigeta, trig->Phi(), asso->Pt(), eta[j], asso->Phi(), 0.1396, 0.1396);
+
+ const Float_t kK0smass = 0.4976;
+
+ if (TMath::Abs(mass - kK0smass*kK0smass) < 0.1)
+ {
+ mass = GetInvMassSquared(trig->Pt(), trigeta, trig->Phi(), asso->Pt(),eta[j], asso->Phi(), 0.1396, 0.1396);
+
+ //fControlConvResoncances->Fill(1, mass - kK0smass*kK0smass);
+
+ if (mass > (kK0smass-0.02)*(kK0smass-0.02) && mass < (kK0smass+0.02)*(kK0smass+0.02))
+ continue;
+ }
+ }
+
+ // Lambda
+ if (fCutResonances && asso->Charge() * trig->Charge() < 0)
+ {
+ Float_t mass1 = GetInvMassSquaredCheap(trig->Pt(), trigeta, trig->Phi(), asso->Pt(), eta[j], asso->Phi(), 0.1396, 0.9383);
+ Float_t mass2 = GetInvMassSquaredCheap(trig->Pt(), trigeta, trig->Phi(), asso->Pt(),eta[j] , asso->Phi(), 0.9383, 0.1396);
+
+ const Float_t kLambdaMass = 1.115;
+
+ if (TMath::Abs(mass1 - kLambdaMass*kLambdaMass) < 0.1)
+ {
+ mass1 = GetInvMassSquared(trig->Pt(), trigeta, trig->Phi(), asso->Pt(),eta[j], asso->Phi(), 0.1396, 0.9383);
+
+ //fControlConvResoncances->Fill(2, mass1 - kLambdaMass*kLambdaMass);
+
+ if (mass1 > (kLambdaMass-0.02)*(kLambdaMass-0.02) && mass1 < (kLambdaMass+0.02)*(kLambdaMass+0.02))
+ continue;
+ }
+ if (TMath::Abs(mass2 - kLambdaMass*kLambdaMass) < 0.1)
+ {
+ mass2 = GetInvMassSquared(trig->Pt(), trigeta, trig->Phi(), asso->Pt(),eta[j] , asso->Phi(), 0.9383, 0.1396);
+
+ //fControlConvResoncances->Fill(2, mass2 - kLambdaMass*kLambdaMass);
+
+ if (mass2 > (kLambdaMass-0.02)*(kLambdaMass-0.02) && mass2 < (kLambdaMass+0.02)*(kLambdaMass+0.02))
+ continue;
+ }
+ }
+
+ if (twoTrackEfficiencyCut)
+ {
+ // the variables & cuthave been developed by the HBT group
+ // see e.g. https://indico.cern.ch/materialDisplay.py?contribId=36&sessionId=6&materialId=slides&confId=142700
+ Float_t phi1 = trig->Phi();
+ Float_t pt1 = trig->Pt();
+ Float_t charge1 = trig->Charge();
+ Float_t phi2 = asso->Phi();
+ Float_t pt2 = asso->Pt();
+ Float_t charge2 = asso->Charge();
+
+ Float_t deta= trigeta - eta[j];
+
+ // optimization
+ if (TMath::Abs(deta) < twoTrackEfficiencyCutValue * 2.5 * 3)
+ {
+
+ // check first boundaries to see if is worth to loop and find the minimum
+ Float_t dphistar1 = GetDPhiStar(phi1, pt1, charge1, phi2, pt2, charge2, 0.8, bSign);
+ Float_t dphistar2 = GetDPhiStar(phi1, pt1, charge1, phi2, pt2, charge2, 2.5, bSign);
+
+ const Float_t kLimit = twoTrackEfficiencyCutValue * 3;
+
+ Float_t dphistarminabs = 1e5;
+ Float_t dphistarmin = 1e5;
+
+ if (TMath::Abs(dphistar1) < kLimit || TMath::Abs(dphistar2) < kLimit || dphistar1 * dphistar2 < 0)
+ {
+ for (Double_t rad=0.8; rad<2.51; rad+=0.01)
+ {
+ Float_t dphistar = GetDPhiStar(phi1, pt1, charge1, phi2, pt2, charge2, rad, bSign);
+
+ Float_t dphistarabs = TMath::Abs(dphistar);
+
+ if (dphistarabs < dphistarminabs)
+ {
+ dphistarmin = dphistar;
+ dphistarminabs = dphistarabs;
+ }
+ }
+
+if (dphistarminabs < twoTrackEfficiencyCutValue && TMath::Abs(deta) < twoTrackEfficiencyCutValue)
+ {
+// Printf("Removed track pair %d %d with %f %f %f %f %f %f %f %f %f", i, j, deta, dphistarminabs, phi1, pt1, charge1, phi2, pt2, charge2, bSign);
+ continue;
+ }
+//fTwoTrackDistancePt[1]->Fill(deta, dphistarmin, TMath::Abs(pt1 - pt2));
+
+ }
+ }
+ }
+
+ Float_t trackeffasso=1.0;
+ if(fapplyAssoefficiency) trackeffasso=asso->geteffcorrectionval();
+ //cout<<"*******************trackeffasso="<<trackeffasso<<endl;
+ Float_t deleta=trigeta-eta[j];
+ Float_t delphi=PhiRange(trigphi-asso->Phi());
+
+ //here get the two particle efficiency correction factor
+ Float_t effweight=trackefftrig*trackeffasso;
+ //cout<<"*******************effweight="<<effweight<<endl;
+ Double_t* vars;
+ vars= new Double_t[kTrackVariablesPair];
+ vars[0]=cent;
+ vars[1]=vtx;
+ vars[2]=trigpt;
+ vars[3]=asso->Pt();
+ vars[4]=deleta;
+ vars[5]=delphi;
+if(fcontainPIDtrig && !fcontainPIDasso) vars[6]=particlepidtrig;
+if(!fcontainPIDtrig && fcontainPIDasso) vars[6]=particlepidasso;
+ if(fcontainPIDtrig && fcontainPIDasso){
+ vars[6]=particlepidtrig;
+ vars[7]=particlepidasso;
+ }
+
+ //Fill Correlation ThnSparses
+ if(fillup=="trigassoUNID")
+ {
+ if(mixcase==kFALSE) fTHnCorrUNID->Fill(vars,1.0/effweight);
+ if(mixcase==kTRUE) fTHnCorrUNIDmix->Fill(vars,1.0/effweight);
+ }
+ if(fillup=="trigIDassoID")
+ {
+ if(mixcase==kFALSE) fTHnCorrID->Fill(vars,1.0/effweight);
+ if(mixcase==kTRUE) fTHnCorrIDmix->Fill(vars,1.0/effweight);
+ }
+ if(fillup=="trigIDassoIDMCTRUTH")//******************************************for TRUTH MC particles
+ {//in this case effweight should be 1 always
+ if(mixcase==kFALSE) fCorrelatonTruthPrimary->Fill(vars,1.0/effweight);
+ if(mixcase==kTRUE) fCorrelatonTruthPrimarymix->Fill(vars,1.0/effweight);
+ }
+ if(fillup=="trigIDassoUNID" || fillup=="trigUNIDassoID")//****************************be careful
+ {
+ if(mixcase==kFALSE) fTHnCorrIDUNID->Fill(vars,1.0/effweight);
+ if(mixcase==kTRUE) fTHnCorrIDUNIDmix->Fill(vars,1.0/effweight);
+ }
+
+delete[] vars;
+ }//asso loop ends
+delete[] trigval;
+ }//trigger loop ends
+
+}
+
+//--------------------------------------------------------------------------------
+Float_t AliTwoParticlePIDCorr::GetTrackbyTrackeffvalue(AliAODTrack* track,Double_t cent,Float_t evzvtx, Int_t parpid)
+{
+ //This function is called only when applyefficiency=kTRUE; also ensure that "track" is present before calling that function
+ Int_t effVars[4];
+ Float_t effvalue=1.;
+
+ if(parpid==unidentified)
+ {
+ effVars[0] = effcorection[5]->GetAxis(0)->FindBin(cent);
+ effVars[1] = effcorection[5]->GetAxis(1)->FindBin(evzvtx);
+ effVars[2] = effcorection[5]->GetAxis(2)->FindBin(track->Pt());
+ effVars[3] = effcorection[5]->GetAxis(3)->FindBin(track->Eta());
+ effvalue=effcorection[5]->GetBinContent(effVars);
+ }
+if(parpid==SpPion || parpid==SpKaon)
+ {
+ if(fmesoneffrequired && !fkaonprotoneffrequired && track->Pt()>=fminPtComboeff && track->Pt()<=fmaxPtComboeff)
+ {
+ effVars[0] = effcorection[3]->GetAxis(0)->FindBin(cent);
+ effVars[1] = effcorection[3]->GetAxis(1)->FindBin(evzvtx);
+ effVars[2] = effcorection[3]->GetAxis(2)->FindBin(track->Pt());
+ effVars[3] = effcorection[3]->GetAxis(3)->FindBin(track->Eta());
+ effvalue=effcorection[3]->GetBinContent(effVars);
+ }
+ else{
+ if(parpid==SpPion)
+ {
+ effVars[0] = effcorection[0]->GetAxis(0)->FindBin(cent);
+ effVars[1] = effcorection[0]->GetAxis(1)->FindBin(evzvtx);
+ effVars[2] = effcorection[0]->GetAxis(2)->FindBin(track->Pt());
+ effVars[3] = effcorection[0]->GetAxis(3)->FindBin(track->Eta());
+ effvalue=effcorection[0]->GetBinContent(effVars);
+ }
+
+ if(parpid==SpKaon)
+ {
+ effVars[0] = effcorection[1]->GetAxis(0)->FindBin(cent);
+ effVars[1] = effcorection[1]->GetAxis(1)->FindBin(evzvtx);
+ effVars[2] = effcorection[1]->GetAxis(2)->FindBin(track->Pt());
+ effVars[3] = effcorection[1]->GetAxis(3)->FindBin(track->Eta());
+ effvalue=effcorection[1]->GetBinContent(effVars);
+ }
+ }
+ }
+
+ if(parpid==SpProton)
+ {
+ effVars[0] = effcorection[2]->GetAxis(0)->FindBin(cent);
+ effVars[1] = effcorection[2]->GetAxis(1)->FindBin(evzvtx);
+ effVars[2] = effcorection[2]->GetAxis(2)->FindBin(track->Pt());
+ effVars[3] = effcorection[2]->GetAxis(3)->FindBin(track->Eta());
+ effvalue=effcorection[2]->GetBinContent(effVars);
+ }
+
+ if(fkaonprotoneffrequired && !fmesoneffrequired && track->Pt()>=fminPtComboeff && track->Pt()<=fmaxPtComboeff)
+ {
+ if(parpid==SpProton || parpid==SpKaon)
+ {
+ effVars[0] = effcorection[4]->GetAxis(0)->FindBin(cent);
+ effVars[1] = effcorection[4]->GetAxis(1)->FindBin(evzvtx);
+ effVars[2] = effcorection[4]->GetAxis(2)->FindBin(track->Pt());
+ effVars[3] = effcorection[4]->GetAxis(3)->FindBin(track->Eta());
+ effvalue=effcorection[4]->GetBinContent(effVars);
+ }
+ }
+ // Printf("%d %d %d %d %f", effVars[0], effVars[1], effVars[2], effVars[3], fEfficiencyCorrectionAssociated->GetBinContent(effVars));
+ if(effvalue==0.) effvalue=1.;
+
+ return effvalue;
+
+}
+//-----------------------------------------------------------------------
+
+Int_t AliTwoParticlePIDCorr::ClassifyTrack(AliAODTrack* track,AliAODVertex* vertex,Float_t magfield)
+{
+
+ if(!track) return 0;
+ Bool_t trackOK = track->TestFilterBit(fFilterBit);
+ if(!trackOK) return 0;
+ //select only primary traks(for data & reco MC tracks)
+ if(fonlyprimarydatareco && track->GetType()!=AliAODTrack::kPrimary) return 0;
+ if(track->Charge()==0) return 0;
+ fHistQA[12]->Fill(track->GetTPCNcls());
+ Float_t dxy, dz;
+ dxy = track->DCA();
+ dz = track->ZAtDCA();
+ fHistQA[6]->Fill(dxy);
+ fHistQA[7]->Fill(dz);
+ Float_t chi2ndf = track->Chi2perNDF();
+ fHistQA[13]->Fill(chi2ndf);
+ Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
+ fHistQA[14]->Fill(nCrossedRowsTPC);
+ //Float_t ratioCrossedRowsOverFindableClustersTPC = 1.0;
+ if (track->GetTPCNclsF()>0) {
+ Float_t ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
+ fHistQA[15]->Fill(ratioCrossedRowsOverFindableClustersTPC);
+ }
+//accepted tracks
+ Float_t pt=track->Pt();
+ if(pt< fminPt || pt> fmaxPt) return 0;
+ if(TMath::Abs(track->Eta())> fmaxeta) return 0;
+ if(track->Phi()<0. || track->Phi()>2*TMath::Pi()) return 0;
+ //if (!HasTPCPID(track)) return 0;//trigger & associated particles must have TPC PID,Is it required
+// DCA XY
+ if (fdcacut && fDCAXYCut)
+ {
+ if (!vertex)
+ return 0;
+
+ Double_t pos[2];
+ Double_t covar[3];
+ AliAODTrack* clone =(AliAODTrack*) track->Clone();
+ Bool_t success = clone->PropagateToDCA(vertex, magfield, fdcacutvalue, pos, covar);
+ delete clone;
+ if (!success)
+ return 0;
+
+// Printf("%f", ((AliAODTrack*)part)->DCA());
+// Printf("%f", pos[0]);
+ if (TMath::Abs(pos[0]) > fDCAXYCut->Eval(track->Pt()))
+ return 0;
+ }
+
+ fHistQA[8]->Fill(pt);
+ fHistQA[9]->Fill(track->Eta());
+ fHistQA[10]->Fill(track->Phi());
+ return 1;
+ }
+ //________________________________________________________________________________
+void AliTwoParticlePIDCorr::CalculateNSigmas(AliAODTrack *track)
+{
+//This function is called within the func GetParticle() for accepted tracks only i.e.after call of Classifytrack() & for those tracks which have proper TPC PID response . combined nsigma(circular) cut only for particles having pt upto 4.0 Gev/c and beyond that use the asymmetric nsigma cut around pion's mean position in TPC ( while filling the TObjArray for trig & asso )
+Float_t pt=track->Pt();
+
+//it is assumed that every track that passed the filterbit have proper TPC response(!!)
+Float_t nsigmaTPCkPion =fPID->NumberOfSigmasTPC(track, AliPID::kPion);
+Float_t nsigmaTPCkKaon =fPID->NumberOfSigmasTPC(track, AliPID::kKaon);
+Float_t nsigmaTPCkProton =fPID->NumberOfSigmasTPC(track, AliPID::kProton);
+
+Float_t nsigmaTOFkProton=999.,nsigmaTOFkKaon=999.,nsigmaTOFkPion=999.;
+Float_t nsigmaTPCTOFkProton=999.,nsigmaTPCTOFkKaon=999.,nsigmaTPCTOFkPion=999.;
+
+ if(HasTOFPID(track) && pt>fPtTOFPIDmin)
+ {
+
+nsigmaTOFkPion =fPID->NumberOfSigmasTOF(track, AliPID::kPion);
+nsigmaTOFkKaon =fPID->NumberOfSigmasTOF(track, AliPID::kKaon);
+nsigmaTOFkProton =fPID->NumberOfSigmasTOF(track, AliPID::kProton);
+//---combined
+nsigmaTPCTOFkPion = TMath::Sqrt(nsigmaTPCkPion*nsigmaTPCkPion+nsigmaTOFkPion*nsigmaTOFkPion);
+nsigmaTPCTOFkKaon = TMath::Sqrt(nsigmaTPCkKaon*nsigmaTPCkKaon+nsigmaTOFkKaon*nsigmaTOFkKaon);
+nsigmaTPCTOFkProton = TMath::Sqrt(nsigmaTPCkProton*nsigmaTPCkProton+nsigmaTOFkProton*nsigmaTOFkProton);
+
+
+ }
+else{
+ // --- combined
+ // if TOF is missing and below fPtTOFPID only the TPC information is used
+ nsigmaTPCTOFkProton = TMath::Abs(nsigmaTPCkProton);
+ nsigmaTPCTOFkKaon = TMath::Abs(nsigmaTPCkKaon);
+ nsigmaTPCTOFkPion = TMath::Abs(nsigmaTPCkPion);
+
+ }
+
+//set data member fnsigmas
+ fnsigmas[SpPion][NSigmaTPC]=nsigmaTPCkPion;
+ fnsigmas[SpKaon][NSigmaTPC]=nsigmaTPCkKaon;
+ fnsigmas[SpProton][NSigmaTPC]=nsigmaTPCkProton;
+
+ //for all tracks below fPtTOFPIDmin and also for tracks above fPtTOFPIDmin without proper TOF response these TOF nsigma values will be 999.
+ fnsigmas[SpPion][NSigmaTOF]=nsigmaTOFkPion;
+ fnsigmas[SpKaon][NSigmaTOF]=nsigmaTOFkKaon;
+ fnsigmas[SpProton][NSigmaTOF]=nsigmaTOFkProton;
+
+ //for all tracks below fPtTOFPIDmin and also for tracks above fPtTOFPIDmin without proper TOF response these TPCTOF nsigma values will be TMath::Abs(TPC only nsigma)
+ fnsigmas[SpPion][NSigmaTPCTOF]=nsigmaTPCTOFkPion;
+ fnsigmas[SpKaon][NSigmaTPCTOF]=nsigmaTPCTOFkKaon;
+ fnsigmas[SpProton][NSigmaTPCTOF]=nsigmaTPCTOFkProton;
+
+
+}
+//----------------------------------------------------------------------------
+Int_t AliTwoParticlePIDCorr::FindMinNSigma(AliAODTrack *track)
+{
+ //this function is always called after calling the function CalculateNSigmas(AliAODTrack *track)
+if(fRequestTOFPID && track->Pt()>fPtTOFPIDmin && track->Pt()<=fPtTOFPIDmax && (!HasTOFPID(track)) )return SpUndefined;//so any track having Pt>0.6 && Pt<=4.0 Gev withot having proper TOF response will be defined as SpUndefined
+//get the identity of the particle with the minimum Nsigma
+ Float_t nsigmaPion=999., nsigmaKaon=999., nsigmaProton=999.;
+ switch (fPIDType){
+ case NSigmaTPC:
+ nsigmaProton = TMath::Abs(fnsigmas[SpProton][NSigmaTPC]);
+ nsigmaKaon = TMath::Abs(fnsigmas[SpKaon][NSigmaTPC]) ;
+ nsigmaPion = TMath::Abs(fnsigmas[SpPion][NSigmaTPC]) ;
+ break;
+ case NSigmaTOF:
+ nsigmaProton = TMath::Abs(fnsigmas[SpProton][NSigmaTOF]);
+ nsigmaKaon = TMath::Abs(fnsigmas[SpKaon][NSigmaTOF]) ;
+ nsigmaPion = TMath::Abs(fnsigmas[SpPion][NSigmaTOF]) ;
+ break;
+ case NSigmaTPCTOF://In case of no TOF matching the combined nsigma is the TPC one
+ nsigmaProton = TMath::Abs(fnsigmas[SpProton][NSigmaTPCTOF]);
+ nsigmaKaon = TMath::Abs(fnsigmas[SpKaon][NSigmaTPCTOF]) ;
+ nsigmaPion = TMath::Abs(fnsigmas[SpPion][NSigmaTPCTOF]) ;
+ break;
+ }
+
+ if(track->Pt()<=fPtTOFPIDmax) {
+ // guess the particle based on the smaller nsigma (within fNSigmaPID)
+ if( ( nsigmaKaon==nsigmaPion ) && ( nsigmaKaon==nsigmaProton )) return SpUndefined;//it is the default value for the three
+if( ( nsigmaKaon < nsigmaPion ) && ( nsigmaKaon < nsigmaProton ) && (nsigmaKaon < fNSigmaPID)) return SpKaon;
+if( ( nsigmaPion < nsigmaKaon ) && ( nsigmaPion < nsigmaProton ) && (nsigmaPion < fNSigmaPID)) return SpPion;
+if( ( nsigmaProton < nsigmaKaon ) && ( nsigmaProton < nsigmaPion ) && (nsigmaProton < fNSigmaPID)) return SpProton;
+
+// else, return undefined
+ return SpUndefined;
+ }
+ else {//asymmetric nsigma cut around pion's mean position for tracks having Pt>4.0 Gev
+ if(fminprotonsigmacut<fnsigmas[SpPion][NSigmaTPC] && fnsigmas[SpPion][NSigmaTPC]<fmaxprotonsigmacut) return SpProton;
+ if(fminpionsigmacut<fnsigmas[SpPion][NSigmaTPC] && fnsigmas[SpPion][NSigmaTPC]<fmaxpionsigmacut) return SpPion;
+// else, return undefined(here we don't detect kaons)
+ return SpUndefined;
+ }
+
+}
+
+//------------------------------------------------------------------------------------------
+Bool_t* AliTwoParticlePIDCorr::GetDoubleCounting(AliAODTrack * trk){
+ //this function is always called after calling the function CalculateNSigmas(AliAODTrack *track)
+
+ //if a particle has double counting set fHasDoubleCounting[ipart]=kTRUE
+ //fill DC histos
+ for(Int_t ipart=0;ipart<NSpecies;ipart++)fHasDoubleCounting[ipart]=kFALSE;//array with kTRUE for second (or third) identity of the track
+
+ Int_t MinNSigma=FindMinNSigma(trk);//not filling the NSigmaRec histos
+
+
+ if(MinNSigma==SpUndefined)return fHasDoubleCounting;//in case of undefined no Double counting
+
+ Float_t nsigmaPion=999., nsigmaKaon=999., nsigmaProton=999.;
+ switch (fPIDType) {
+ case NSigmaTPC:
+ nsigmaProton = TMath::Abs(fnsigmas[SpProton][NSigmaTPC]);
+ nsigmaKaon = TMath::Abs(fnsigmas[SpKaon][NSigmaTPC]) ;
+ nsigmaPion = TMath::Abs(fnsigmas[SpPion][NSigmaTPC]) ;
+ break;
+ case NSigmaTOF:
+ nsigmaProton = TMath::Abs(fnsigmas[SpProton][NSigmaTOF]);
+ nsigmaKaon = TMath::Abs(fnsigmas[SpKaon][NSigmaTOF]) ;
+ nsigmaPion = TMath::Abs(fnsigmas[SpPion][NSigmaTOF]) ;
+ break;
+ case NSigmaTPCTOF://In case of no TOF matching the combined nsigma is the TPC one
+ nsigmaProton = TMath::Abs(fnsigmas[SpProton][NSigmaTPCTOF]);
+ nsigmaKaon = TMath::Abs(fnsigmas[SpKaon][NSigmaTPCTOF]) ;
+ nsigmaPion = TMath::Abs(fnsigmas[SpPion][NSigmaTPCTOF]) ;
+ break;
+ }
+
+ //there is chance of overlapping only for particles having pt below 4.0 GEv
+ if(trk->Pt()<=4.0){
+ if(nsigmaPion<fNSigmaPID && MinNSigma!=SpPion)fHasDoubleCounting[SpPion]=kTRUE;
+ if(nsigmaKaon<fNSigmaPID && MinNSigma!=SpKaon)fHasDoubleCounting[SpKaon]=kTRUE;
+ if(nsigmaProton<fNSigmaPID && MinNSigma!=SpProton)fHasDoubleCounting[SpProton]=kTRUE;
+
+ }
+
+
+ return fHasDoubleCounting;
+}
+
+//////////////////////////////////////////////////////////////////////////////////////////////////
+Int_t AliTwoParticlePIDCorr::GetParticle(AliAODTrack * trk){
+ //return the specie according to the minimum nsigma value
+ //no double counting, this has to be evaluated using CheckDoubleCounting()
+ //Printf("fPtTOFPID %.1f, fRequestTOFPID %d, fNSigmaPID %.1f, fPIDType %d",fPtTOFPID,fRequestTOFPID,fNSigmaPID,fPIDType);
+
+ CalculateNSigmas(trk);//fill the data member fnsigmas with the nsigmas value [ipart][iPID]
+
+ if(fUseExclusiveNSigma){
+ Bool_t *HasDC;
+ HasDC=GetDoubleCounting(trk);
+ for(Int_t ipart=0;ipart<NSpecies;ipart++){
+ if(HasDC[ipart]==kTRUE) return SpUndefined;
+ }
+ return FindMinNSigma(trk);//NSigmaRec distr filled here
+ }
+ else return FindMinNSigma(trk);//NSigmaRec distr filled here
+
+}
+
+//-------------------------------------------------------------------------------------
+Bool_t
+AliTwoParticlePIDCorr::HasTPCPID(AliAODTrack *track) const
+{
+ // check PID signal
+ AliPIDResponse::EDetPidStatus statustpc = fPID->CheckPIDStatus(AliPIDResponse::kTPC,track);
+ if(statustpc!=AliPIDResponse::kDetPidOk) return kFALSE;
+ //ULong_t status=track->GetStatus();
+ //if (!( (status & AliAODTrack::kTPCpid ) == AliAODTrack::kTPCpid )) return kFALSE;//remove light nuclei
+ //if (track->GetTPCsignal() <= 0.) return kFALSE;
+ // if(track->GetTPCsignalN() < 60) return kFALSE;//tracks with TPCsignalN< 60 have questionable dEdx,cutting on TPCsignalN > 70 or > 60 shouldn't make too much difference in statistics,also it is IMO safe to use TPC also for MIPs.
+
+ return kTRUE;
+}
+//___________________________________________________________
+
+Bool_t
+AliTwoParticlePIDCorr::HasTOFPID(AliAODTrack *track) const
+{
+ // check TOF matched track
+ //ULong_t status=track->GetStatus();
+ //if (!( (status & AliAODTrack::kITSin ) == AliAODTrack::kITSin )) return kFALSE;
+ AliPIDResponse::EDetPidStatus statustof = fPID->CheckPIDStatus(AliPIDResponse::kTOF,track);
+ if(statustof!= AliPIDResponse::kDetPidOk) return kFALSE;
+ if(track->Pt()<=fPtTOFPIDmin) return kFALSE;
+ //if(!((status & AliAODTrack::kTOFpid ) == AliAODTrack::kTOFpid )) return kFALSE;
+ //Float_t probMis = fPIDresponse->GetTOFMismatchProbability(track);
+ // if (probMis > 0.01) return kFALSE;
+if(fRemoveTracksT0Fill)
+ {
+Int_t startTimeMask = fPID->GetTOFResponse().GetStartTimeMask(track->P());
+ if (startTimeMask < 0)return kFALSE;
+ }
+ return kTRUE;
+}
+
+//________________________________________________________________________
+Float_t AliTwoParticlePIDCorr :: GetBeta(AliAODTrack *track)
+{
+ //it is called only when TOF PID is available
+ Double_t p = track->P();
+ Double_t time=track->GetTOFsignal()-fPID->GetTOFResponse().GetStartTime(p);
+ Double_t timei[5];
+ track->GetIntegratedTimes(timei);
+ return timei[0]/time;
+}
+//------------------------------------------------------------------------------------------------------
+
+Float_t AliTwoParticlePIDCorr::GetInvMassSquared(Float_t pt1, Float_t eta1, Float_t phi1, Float_t pt2, Float_t eta2, Float_t phi2, Float_t m0_1, Float_t m0_2)
+{
+ // calculate inv mass squared
+ // same can be achieved, but with more computing time with
+ /*TLorentzVector photon, p1, p2;
+ p1.SetPtEtaPhiM(triggerParticle->Pt(), triggerEta, triggerParticle->Phi(), 0.510e-3);
+ p2.SetPtEtaPhiM(particle->Pt(), eta[j], particle->Phi(), 0.510e-3);
+ photon = p1+p2;
+ photon.M()*/
+
+ Float_t tantheta1 = 1e10;
+
+ if (eta1 < -1e-10 || eta1 > 1e-10)
+ tantheta1 = 2 * TMath::Exp(-eta1) / ( 1 - TMath::Exp(-2*eta1));
+
+ Float_t tantheta2 = 1e10;
+ if (eta2 < -1e-10 || eta2 > 1e-10)
+ tantheta2 = 2 * TMath::Exp(-eta2) / ( 1 - TMath::Exp(-2*eta2));
+
+ Float_t e1squ = m0_1 * m0_1 + pt1 * pt1 * (1.0 + 1.0 / tantheta1 / tantheta1);
+ Float_t e2squ = m0_2 * m0_2 + pt2 * pt2 * (1.0 + 1.0 / tantheta2 / tantheta2);
+
+ Float_t mass2 = m0_1 * m0_1 + m0_2 * m0_2 + 2 * ( TMath::Sqrt(e1squ * e2squ) - ( pt1 * pt2 * ( TMath::Cos(phi1 - phi2) + 1.0 / tantheta1 / tantheta2 ) ) );
+
+ return mass2;
+}
+//---------------------------------------------------------------------------------
+
+Float_t AliTwoParticlePIDCorr::GetInvMassSquaredCheap(Float_t pt1, Float_t eta1, Float_t phi1, Float_t pt2, Float_t eta2, Float_t phi2, Float_t m0_1, Float_t m0_2)
+{
+ // calculate inv mass squared approximately
+
+ Float_t tantheta1 = 1e10;
+
+ if (eta1 < -1e-10 || eta1 > 1e-10)
+ {
+ Float_t expTmp = 1.0-eta1+eta1*eta1/2-eta1*eta1*eta1/6+eta1*eta1*eta1*eta1/24;
+ tantheta1 = 2.0 * expTmp / ( 1.0 - expTmp*expTmp);
+ }
+
+ Float_t tantheta2 = 1e10;
+ if (eta2 < -1e-10 || eta2 > 1e-10)
+ {
+ Float_t expTmp = 1.0-eta2+eta2*eta2/2-eta2*eta2*eta2/6+eta2*eta2*eta2*eta2/24;
+ tantheta2 = 2.0 * expTmp / ( 1.0 - expTmp*expTmp);
+ }
+
+ Float_t e1squ = m0_1 * m0_1 + pt1 * pt1 * (1.0 + 1.0 / tantheta1 / tantheta1);
+ Float_t e2squ = m0_2 * m0_2 + pt2 * pt2 * (1.0 + 1.0 / tantheta2 / tantheta2);
+
+ // fold onto 0...pi
+ Float_t deltaPhi = TMath::Abs(phi1 - phi2);
+ while (deltaPhi > TMath::TwoPi())
+ deltaPhi -= TMath::TwoPi();
+ if (deltaPhi > TMath::Pi())
+ deltaPhi = TMath::TwoPi() - deltaPhi;
+
+ Float_t cosDeltaPhi = 0;
+ if (deltaPhi < TMath::Pi()/3)
+ cosDeltaPhi = 1.0 - deltaPhi*deltaPhi/2 + deltaPhi*deltaPhi*deltaPhi*deltaPhi/24;
+ else if (deltaPhi < 2*TMath::Pi()/3)
+ cosDeltaPhi = -(deltaPhi - TMath::Pi()/2) + 1.0/6 * TMath::Power((deltaPhi - TMath::Pi()/2), 3);
+ else
+ cosDeltaPhi = -1.0 + 1.0/2.0*(deltaPhi - TMath::Pi())*(deltaPhi - TMath::Pi()) - 1.0/24.0 * TMath::Power(deltaPhi - TMath::Pi(), 4);
+
+ Float_t mass2 = m0_1 * m0_1 + m0_2 * m0_2 + 2 * ( TMath::Sqrt(e1squ * e2squ) - ( pt1 * pt2 * ( cosDeltaPhi + 1.0 / tantheta1 / tantheta2 ) ) );
+
+// Printf(Form("%f %f %f %f %f %f %f %f %f", pt1, eta1, phi1, pt2, eta2, phi2, m0_1, m0_2, mass2));
+
+ return mass2;
+}
+//--------------------------------------------------------------------------------
+Float_t AliTwoParticlePIDCorr::GetDPhiStar(Float_t phi1, Float_t pt1, Float_t charge1, Float_t phi2, Float_t pt2, Float_t charge2, Float_t radius, Float_t bSign)
+{
+ //
+ // calculates dphistar
+ //
+
+ Float_t dphistar = phi1 - phi2 - charge1 * bSign * TMath::ASin(0.075 * radius / pt1) + charge2 * bSign * TMath::ASin(0.075 * radius / pt2);
+
+ static const Double_t kPi = TMath::Pi();
+
+ // circularity
+// if (dphistar > 2 * kPi)
+// dphistar -= 2 * kPi;
+// if (dphistar < -2 * kPi)
+// dphistar += 2 * kPi;
+
+ if (dphistar > kPi)
+ dphistar = kPi * 2 - dphistar;
+ if (dphistar < -kPi)
+ dphistar = -kPi * 2 - dphistar;
+ if (dphistar > kPi) // might look funny but is needed
+ dphistar = kPi * 2 - dphistar;
+
+ return dphistar;
+}
+//_________________________________________________________________________
+void AliTwoParticlePIDCorr ::DefineEventPool()
+{
+const Int_t MaxNofEvents=1000;
+const Int_t MaxNofTracks=50000;
+const Int_t NofVrtxBins=10+(1+10)*2;
+Double_t ZvrtxBins[NofVrtxBins+1]={ -10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10,
+ 90, 92, 94, 96, 98, 100, 102, 104, 106, 108, 110,
+ 190, 192, 194, 196, 198, 200, 202, 204, 206, 208, 210
+ };
+ if (fSampleType=="pp"){
+const Int_t NofCentBins=5;
+Double_t CentralityBins[NofCentBins+1]={0.,10., 20., 40.,80.,200.1};
+fPoolMgr = new AliEventPoolManager(MaxNofEvents,MaxNofTracks,NofCentBins,CentralityBins,NofVrtxBins,ZvrtxBins);
+ }
+if (fSampleType=="pPb" || fSampleType=="PbPb")
+ {
+const Int_t NofCentBins=15;
+Double_t CentralityBins[NofCentBins+1]={0., 1., 2., 3., 4., 5., 10., 20., 30., 40., 50., 60., 70., 80., 90., 100.1 };
+fPoolMgr = new AliEventPoolManager(MaxNofEvents,MaxNofTracks,NofCentBins,CentralityBins,NofVrtxBins,ZvrtxBins);
+ }
+fPoolMgr->SetTargetValues(MaxNofTracks, 0.1, 5);
+
+//if(!fPoolMgr) return kFALSE;
+//return kTRUE;
+
+}
+//------------------------------------------------------------------------
+Double_t* AliTwoParticlePIDCorr::GetBinning(const char* configuration, const char* tag, Int_t& nBins)
+{
+ // This method is a copy from AliUEHist::GetBinning
+ // takes the binning from <configuration> identified by <tag>
+ // configuration syntax example:
+ // eta: 2.4, -2.3, -2.2, -2.1, -2.0, -1.9, -1.8, -1.7, -1.6, -1.5, -1.4, -1.3, -1.2, -1.1, -1.0, -0.9, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, 2.1, 2.2, 2.3, 2.4
+ // phi: .....
+ //
+ // returns bin edges which have to be deleted by the caller
+
+ TString config(configuration);
+ TObjArray* lines = config.Tokenize("\n");
+ for (Int_t i=0; i<lines->GetEntriesFast(); i++)
+ {
+ TString line(lines->At(i)->GetName());
+ if (line.BeginsWith(TString(tag) + ":"))
+ {
+ line.Remove(0, strlen(tag) + 1);
+ line.ReplaceAll(" ", "");
+ TObjArray* binning = line.Tokenize(",");
+ Double_t* bins = new Double_t[binning->GetEntriesFast()];
+ for (Int_t j=0; j<binning->GetEntriesFast(); j++)
+ bins[j] = TString(binning->At(j)->GetName()).Atof();
+
+ nBins = binning->GetEntriesFast() - 1;
+
+ delete binning;
+ delete lines;
+ return bins;
+ }
+ }
+
+ delete lines;
+ AliFatal(Form("Tag %s not found in %s", tag, configuration));
+ return 0;
+}
+//_______________________________________________________________________________
+void AliTwoParticlePIDCorr::SetAsymmetricBin(THnSparse *h,Int_t dim,Double_t *arraybin,Int_t arraybinsize,TString axisTitle)
+{
+ TAxis *axis = 0x0;
+ axis = h->GetAxis(dim);
+ axis->Set(arraybinsize,arraybin);
+ axis->SetTitle(axisTitle);
+}
+
+//________________________________________________________________________
+void AliTwoParticlePIDCorr::Terminate(Option_t *)
+{
+ // Draw result to screen, or perform fitting, normalizations
+ // Called once at the end of the query
+ fOutput = dynamic_cast<TList*> (GetOutputData(1));
+ if(!fOutput) { Printf("ERROR: could not retrieve TList fOutput"); return; }
+
+
+}
+//------------------------------------------------------------------
+
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff