* about the suitability of this software for any purpose. It is *
* provided "as is" without express or implied warranty. *
**************************************************************************/
-//Source code::dphicorrelations, TaskBFpsi, AliHelperPID
#include "AliTwoParticlePIDCorr.h"
#include "AliVParticle.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TH3F.h"
+#include "TProfile.h"
#include "TList.h"
#include "TFile.h"
#include "TGrid.h"
-
+#include "TExMap.h"
#include "AliCentrality.h"
#include "Riostream.h"
+#include "AliAnalysisDataSlot.h"
+ #include "AliAnalysisDataContainer.h"
+
#include "AliTHn.h"
#include "AliCFContainer.h"
#include "THn.h"
#include "THnSparse.h"
-
+#include "TBits.h"
#include <TSpline.h>
#include <AliPID.h>
#include "AliESDpid.h"
#include <AliPIDResponse.h>
#include "AliPIDCombined.h"
-#include <AliAnalysisManager.h>
#include <AliInputEventHandler.h>
#include "AliAODInputHandler.h"
#include "AliAODEvent.h"
#include "AliAODTrack.h"
#include "AliVTrack.h"
+#include "AliAODv0.h"
+#include "AliAODcascade.h"
#include "THnSparse.h"
#include "AliGenCocktailEventHeader.h"
#include "AliGenEventHeader.h"
#include "AliCollisionGeometry.h"
+#include "AliOADBContainer.h"
#include "AliEventPoolManager.h"
#include "AliAnalysisUtils.h"
const char * kPIDTypeName[]={"TPC","TOF","TPC-TOF"} ;
const char * kDetectorName[]={"ITS","TPC","TOF"} ;
const char * kParticleSpeciesName[]={"Pions","Kaons","Protons","Undefined"} ;
+//Source code::dphicorrelations,VnV0, TaskBFpsi, AliHelperPID,
//________________________________________________________________________
AliTwoParticlePIDCorr::AliTwoParticlePIDCorr() // All data members should be initialised here
:AliAnalysisTaskSE(),
fOutput(0),
fOutputList(0),
+ fList(0),
fCentralityMethod("V0A"),
+ fPPVsMultUtils(kFALSE),
fSampleType("pPb"),
fRequestEventPlane(kFALSE),
+ fRequestEventPlanemixing(kFALSE),
fnTracksVertex(1), // QA tracks pointing to principal vertex (= 3 default)
trkVtx(0),
zvtx(0),
fFilterBit(768),
fTrackStatus(0),
fSharedClusterCut(-1),
+ fSharedTPCmapCut(-1),
+ fSharedfraction_Pair_cut(-1),
fVertextype(1),
skipParticlesAbove(0),
fzvtxcut(10.0),
+ fVxMax_MC(0.3),
+ fVyMax_MC(0.3),
+ fVzMax_MC(10.),
ffilltrigassoUNID(kFALSE),
ffilltrigUNIDassoID(kFALSE),
ffilltrigIDassoUNID(kTRUE),
fmaxPtAsso(0),
fmincentmult(0),
fmaxcentmult(0),
+ fPriHistShare(0),
fhistcentrality(0),
+ fhistImpactParm(0),
+ fhistImpactParmvsMult(0x0),
+ fNchNpartCorr(0x0),
fEventCounter(0),
fEtaSpectrasso(0),
fphiSpectraasso(0),
fEventnobaryon(0),
fEventnomeson(0),
fhistJetTrigestimate(0),
+fTwoTrackDistancePtdip(0x0),
+fTwoTrackDistancePtdipmix(0x0),
fCentralityCorrelation(0x0),
fHistVZEROAGainEqualizationMap(0),
fHistVZEROCGainEqualizationMap(0),
fHistVZEROCvsVZEROAmultiplicity(0x0),
fHistEQVZEROCvsEQVZEROAmultiplicity(0x0),
fHistVZEROSignal(0x0),
-fHistEventPlaneReco(0x0),
fHistEventPlaneTruth(0x0),
fHistPsiMinusPhi(0x0),
+fEventPlanePID(0x0),
+evplaneMC(999.),
+ fgPsi2v0a(999.),
+ fgPsi2v0c(999.),
+ fgPsi2tpc(999.),
+ fgPsi3v0a(999.),
+ fgPsi3v0c(999.),
+ fgPsi3tpc(999.),
+ fgPsi2v0aMC(999.),
+ fgPsi2v0cMC(999.),
+ fgPsi2tpcMC(999.),
+ fgPsi3v0aMC(999.),
+ fgPsi3v0cMC(999.),
+ fgPsi3tpcMC(999.),
+ gReactionPlane(999.),
+ fV2(kTRUE),
+ fV3(kFALSE),
+ fIsAfter2011(kTRUE),
+ fRun(-1),
+ fNcluster(70),
+ fEPdet("V0A"),
+ fMultV0(NULL),
+ fV0Cpol(100),
+ fV0Apol(100),
+ fHResTPCv0A2(NULL),
+fHResTPCv0C2(NULL),
+fHResv0Cv0A2(NULL),
+fHResTPCv0A3(NULL),
+fHResTPCv0C3(NULL),
+fHResv0Cv0A3(NULL),
+ fHResMA2(NULL),
+fHResMC2(NULL),
+fHResAC2(NULL),
+fHResMA3(NULL),
+fHResMC3(NULL),
+fHResAC3(NULL),
+fPhiRPTPC(NULL),
+fPhiRPTPCv3(NULL),
+fPhiRPv0A(NULL),
+fPhiRPv0C(NULL),
+fPhiRPv0Av3(NULL),
+fPhiRPv0Cv3(NULL),
fControlConvResoncances(0),
fHistoTPCdEdx(0x0),
fHistoTOFbeta(0x0),
fAnalysisType("AOD"),
fefffilename(""),
ftwoTrackEfficiencyCutDataReco(kTRUE),
+fTwoTrackCutMinRadius(0.8),
+fTwoTrackCutMaxRadius(2.5),
twoTrackEfficiencyCutValue(0.02),
fPID(NULL),
fPIDCombined(NULL),
fmesoneffrequired(kFALSE),
fkaonprotoneffrequired(kFALSE),
fAnalysisUtils(0x0),
- fDCAXYCut(0)
+ fDCAXYCut(0),
+ fV0TrigCorr(kFALSE),
+ fUsev0DaughterPID(kFALSE),
+ fMinPtDaughter(1.0),// v0 related cut starts here
+ fMaxPtDaughter(4.0),
+ fDCAToPrimVtx(0.1),
+ fMaxDCADaughter(1.0),
+ fMinCPA(0.998),
+ lMax(100),
+ fHistRawPtCentInvK0s(0x0),
+ fHistRawPtCentInvLambda(0x0),
+ fHistRawPtCentInvAntiLambda(0x0),
+ fHistFinalPtCentInvK0s(0x0),
+ fHistFinalPtCentInvLambda(0x0),
+ fHistFinalPtCentInvAntiLambda(0x0),
+ NCtau(3.0),
+fCutctauK0s(2.68),
+ fCutctauLambda(7.8),
+ fCutctauAntiLambda(7.8),
+ fRapCutK0s(0.7),
+ fRapCutLambda(0.7),
+fDaugNClsTPC(70),
+fFracTPCcls(0)
{
for ( Int_t i = 0; i < 16; i++) {
fTrackHistEfficiency[i] = NULL;
effcorection[i]=NULL;
//effmap[i]=NULL;
-
}
+ for ( Int_t i = 0; i < 2; i++ ){
+ fTwoTrackDistancePt[i]=NULL;
+ fTwoTrackDistancePtmix[i]=NULL;
+}
for(Int_t ipart=0;ipart<NSpecies;ipart++)
for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++)
for(Int_t ipart=0;ipart<NSpecies;ipart++) {fHasDoubleCounting[ipart]=kFALSE;}
+ for(Int_t i = 0; i != 2; ++i)
+ for(Int_t j = 0; j != 2; ++j)
+ for(Int_t iC = 0; iC < 9; iC++){
+ fMeanQ[iC][i][j] = 0.;
+ fWidthQ[iC][i][j] = 1.;
+ fMeanQv3[iC][i][j] = 0.;
+ fWidthQv3[iC][i][j] = 1.;
+ }
+
}
//________________________________________________________________________
AliTwoParticlePIDCorr::AliTwoParticlePIDCorr(const char *name) // All data members should be initialised here
:AliAnalysisTaskSE(name),
fOutput(0),
fOutputList(0),
+ fList(0),
fCentralityMethod("V0A"),
+ fPPVsMultUtils(kFALSE),
fSampleType("pPb"),
fRequestEventPlane(kFALSE),
+ fRequestEventPlanemixing(kFALSE),
fnTracksVertex(1), // QA tracks pointing to principal vertex (= 3 default)
trkVtx(0),
zvtx(0),
fFilterBit(768),
fTrackStatus(0),
fSharedClusterCut(-1),
+ fSharedTPCmapCut(-1),
+ fSharedfraction_Pair_cut(-1),
fVertextype(1),
skipParticlesAbove(0),
fzvtxcut(10.0),
+ fVxMax_MC(0.3),
+ fVyMax_MC(0.3),
+ fVzMax_MC(10.),
ffilltrigassoUNID(kFALSE),
ffilltrigUNIDassoID(kFALSE),
ffilltrigIDassoUNID(kTRUE),
fminPtAsso(0),
fmaxPtAsso(0),
fmincentmult(0),
- fmaxcentmult(0),
+ fmaxcentmult(0),
+ fPriHistShare(0),
fhistcentrality(0),
+ fhistImpactParm(0),
+ fhistImpactParmvsMult(0x0),
+ fNchNpartCorr(0x0),
fEventCounter(0),
fEtaSpectrasso(0),
fphiSpectraasso(0),
fEventnobaryon(0),
fEventnomeson(0),
fhistJetTrigestimate(0),
+fTwoTrackDistancePtdip(0x0),
+fTwoTrackDistancePtdipmix(0x0),
fCentralityCorrelation(0x0),
fHistVZEROAGainEqualizationMap(0),
fHistVZEROCGainEqualizationMap(0),
fHistVZEROCvsVZEROAmultiplicity(0x0),
fHistEQVZEROCvsEQVZEROAmultiplicity(0x0),
fHistVZEROSignal(0x0),
-fHistEventPlaneReco(0x0),
fHistEventPlaneTruth(0x0),
-fHistPsiMinusPhi(0x0),
+ fHistPsiMinusPhi(0x0),
+fEventPlanePID(0x0),
+evplaneMC(999.),
+ fgPsi2v0a(999.),
+ fgPsi2v0c(999.),
+ fgPsi2tpc(999.),
+ fgPsi3v0a(999.),
+ fgPsi3v0c(999.),
+ fgPsi3tpc(999.),
+ fgPsi2v0aMC(999.),
+ fgPsi2v0cMC(999.),
+ fgPsi2tpcMC(999.),
+ fgPsi3v0aMC(999.),
+ fgPsi3v0cMC(999.),
+ fgPsi3tpcMC(999.),
+ gReactionPlane(999.),
+ fV2(kTRUE),
+ fV3(kFALSE),
+ fIsAfter2011(kTRUE),
+ fRun(-1),
+ fNcluster(70),
+ fEPdet("V0A"),
+ fMultV0(NULL),
+ fV0Cpol(100),
+ fV0Apol(100),
+ fHResTPCv0A2(NULL),
+fHResTPCv0C2(NULL),
+fHResv0Cv0A2(NULL),
+fHResTPCv0A3(NULL),
+fHResTPCv0C3(NULL),
+fHResv0Cv0A3(NULL),
+ fHResMA2(NULL),
+fHResMC2(NULL),
+fHResAC2(NULL),
+fHResMA3(NULL),
+fHResMC3(NULL),
+fHResAC3(NULL),
+fPhiRPTPC(NULL),
+fPhiRPTPCv3(NULL),
+fPhiRPv0A(NULL),
+fPhiRPv0C(NULL),
+fPhiRPv0Av3(NULL),
+fPhiRPv0Cv3(NULL),
fControlConvResoncances(0),
fHistoTPCdEdx(0x0),
fHistoTOFbeta(0x0),
fAnalysisType("AOD"),
fefffilename(""),
ftwoTrackEfficiencyCutDataReco(kTRUE),
+fTwoTrackCutMinRadius(0.8),
+fTwoTrackCutMaxRadius(2.5),
twoTrackEfficiencyCutValue(0.02),
fPID(NULL),
fPIDCombined(NULL),
fmesoneffrequired(kFALSE),
fkaonprotoneffrequired(kFALSE),
fAnalysisUtils(0x0),
- fDCAXYCut(0)
+ fDCAXYCut(0),
+ fV0TrigCorr(kFALSE),
+ fUsev0DaughterPID(kFALSE),
+ fMinPtDaughter(1.0),// v0 related cut starts here
+ fMaxPtDaughter(4.0),
+ fDCAToPrimVtx(0.1),
+ fMaxDCADaughter(1.0),
+ fMinCPA(0.998),
+ lMax(100),
+ fHistRawPtCentInvK0s(0x0),
+ fHistRawPtCentInvLambda(0x0),
+ fHistRawPtCentInvAntiLambda(0x0),
+ fHistFinalPtCentInvK0s(0x0),
+ fHistFinalPtCentInvLambda(0x0),
+ fHistFinalPtCentInvAntiLambda(0x0),
+ NCtau(3.0),
+fCutctauK0s(2.68),
+ fCutctauLambda(7.8),
+ fCutctauAntiLambda(7.8),
+ fRapCutK0s(0.7),
+ fRapCutLambda(0.7),
+fDaugNClsTPC(70),
+fFracTPCcls(0)
+
+ // The last in the above list should not have a comma after it
+
{
for ( Int_t i = 0; i < 16; i++) {
//effmap[i]=NULL;
}
+for ( Int_t i = 0; i < 2; i++ ){
+ fTwoTrackDistancePt[i]=NULL;
+ fTwoTrackDistancePtmix[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
+ for(Int_t i = 0; i != 2; ++i)
+ for(Int_t j = 0; j != 2; ++j)
+ for(Int_t iC = 0; iC < 9; iC++){
+ fMeanQ[iC][i][j] = 0.;
+ fWidthQ[iC][i][j] = 1.;
+ fMeanQv3[iC][i][j] = 0.;
+ fWidthQv3[iC][i][j] = 1.;
+ }
// 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
-
+ DefineInput(0, TChain::Class());
+
DefineOutput(1, TList::Class()); // for output list
DefineOutput(2, TList::Class());
+ DefineOutput(3, TList::Class());
}
}
+if(fRequestEventPlane){
+if (fList && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
+ delete fList;
+ }
+ }
+
if (fPID) delete fPID;
if (fPIDCombined) delete fPIDCombined;
{
// 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
// global switch disabling the reference
// (to avoid "Replacing existing TH1" if several wagons are created in train)
Bool_t oldStatus = TH1::AddDirectoryStatus();
TH1::AddDirectory(kFALSE);
+ const Int_t nPsiTOF = 10;
+ const Int_t nCentrBin = 9;
+
+
fOutput = new TList();
fOutput->SetOwner(); // IMPORTANT!
fOutputList = new TList;
fOutputList->SetOwner();
fOutputList->SetName("PIDQAList");
-
+
+ if(fRequestEventPlane){
+ fList = new TList;
+ fList->SetOwner();
+ fList->SetName("EPQAList");
+ }
fEventCounter = new TH1F("fEventCounter","EventCounter", 19, 0.5,19.5);
fEventCounter->GetXaxis()->SetBinLabel(1,"Event Accesed");
fEventCounter->GetXaxis()->SetBinLabel(3,"After PileUP Cut");//only for Data
fEventCounter->GetXaxis()->SetBinLabel(5,"Have A Vertex");
fEventCounter->GetXaxis()->SetBinLabel(7,"After vertex Cut");
- fEventCounter->GetXaxis()->SetBinLabel(9,"After centrality flattening");
- fEventCounter->GetXaxis()->SetBinLabel(11,"Within 0-100% centrality");
- fEventCounter->GetXaxis()->SetBinLabel(13,"Event Analyzed");
+ fEventCounter->GetXaxis()->SetBinLabel(9,"Getting centrality");
+ fEventCounter->GetXaxis()->SetBinLabel(11,"After centrality flattening");
+ fEventCounter->GetXaxis()->SetBinLabel(13,"Within 0-100% centrality");
+ fEventCounter->GetXaxis()->SetBinLabel(15,"Event Analyzed");
//fEventCounter->GetXaxis()->SetBinLabel(8,"Event Analysis finished");
fOutput->Add(fEventCounter);
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);
+ if(fSampleType=="pPb" || fSampleType=="PbPb" || fPPVsMultUtils==kTRUE || fCentralityMethod == "MC_b"){ fCentralityCorrelation = new TH2D("fCentralityCorrelation", ";centrality_ImpactParam;multiplicity", 101, 0, 101, 20000, 0,40000);
fOutput->Add(fCentralityCorrelation);
}
fOutput->Add(fHistCentStats);
}
-if(fCentralityMethod.EndsWith("_MANUAL"))
+if(fSampleType=="pp_2_76" || fCentralityMethod.EndsWith("_MANUAL") || (fSampleType=="pp_7" && fPPVsMultUtils==kFALSE))
{
fhistcentrality=new TH1F("fhistcentrality","referencemultiplicity",30001,-0.5,30000.5);
fOutput->Add(fhistcentrality);
fhistcentrality=new TH1F("fhistcentrality","centrality",220,-5,105);
fOutput->Add(fhistcentrality);
}
+ if(fCentralityMethod=="MC_b"){
+fhistImpactParm=new TH1F("fhistImpactParm","Impact_Parameter",300,0,300);
+fOutput->Add(fhistImpactParm);
+fhistImpactParmvsMult=new TH2F("fhistImpactParmvsMult","fhistImpactParmvsMult",300,0,300,5001,-0.5,50000.5);
+fOutput->Add(fhistImpactParmvsMult);
+ }
-if(fCentralityMethod.EndsWith("_MANUAL"))
- {
-TString gmultName[4] = {"V0A_MANUAL","V0C_MANUAL","V0M_MANUAL","TRACKS_MANUAL"};
+ if(fAnalysisType =="MCAOD" || fAnalysisType =="MC"){
+fNchNpartCorr=new TH2F("fNchNpartCorr","fNchNpartCorr",500,0.0,500.0,5001,-0.5,50000.5);
+fOutput->Add(fNchNpartCorr);
+ }
+
+ TString gmultName[4] = {"V0A_MANUAL","V0C_MANUAL","V0M_MANUAL","TRACKS_MANUAL"};
fHistRefmult = new TH2F("fHistRefmult",
"Reference multiplicity",
4,-0.5,3.5,10000,0,20000);
for(Int_t i = 1; i <= 4; i++){
fHistRefmult->GetXaxis()->SetBinLabel(i,gmultName[i-1].Data());
- //fHistCentStatsUsed->GetXaxis()->SetBinLabel(i,gCentName[i-1].Data());
}
fOutput->Add(fHistRefmult);
-
+ if(fCentralityMethod == "V0A_MANUAL" || fCentralityMethod == "V0M_MANUAL" || fCentralityMethod == "V0C_MANUAL" ){
//TPC vs EQVZERO multiplicity
fHistEQVZEROvsTPCmultiplicity = new TH2F("fHistEQVZEROvsTPCmultiplicity","EqVZERO vs TPC multiplicity",10001,-0.5,10000.5,4001,-0.5,4000.5);
fHistEQVZEROvsTPCmultiplicity->GetXaxis()->SetTitle("EqVZERO multiplicity (a.u.)");
fHistVZEROSignal = new TH2F("fHistVZEROSignal","VZERO signal vs VZERO channel;VZERO channel; Signal (a.u.)",64,0.5,64.5,3001,-0.5,30000.5);
fOutput->Add(fHistVZEROSignal);
-}
+ }
- if(fSampleType=="PbPb" && fRequestEventPlane){
-//Event plane
- fHistEventPlaneReco = new TH2F("fHistEventPlaneReco",";#Psi_{2} [deg.];Centrality percentile;Counts",100,0,360.,220,-5,105);
- fOutput->Add(fHistEventPlaneReco);
+ if(fRequestEventPlane){
//Event plane
- fHistEventPlaneTruth = new TH2F("fHistEventPlaneTruth",";#Psi_{2} [deg.];Centrality percentile;Counts",100,0,360.,220,-5,105);
- fOutput->Add(fHistEventPlaneTruth);
-
+
fHistPsiMinusPhi = new TH2D("fHistPsiMinusPhi","",4,-0.5,3.5,100,0,2.*TMath::Pi());
- fOutput->Add(fHistPsiMinusPhi);
+ fList->Add(fHistPsiMinusPhi);
+
+ fEventPlanePID = new TH3F("fEventPlanePID",";centrality;eventplane;PID",20,0.0,100.0,4,-0.5,3.5,4,-0.5,3.5);
+ fList->Add(fEventPlanePID);
}
if(ffillhistQAReco)
{
- fPionPt = new TH1F("fHistQAPionPt","p_{T} distribution",200,0.,10.);
+ fPionPt = new TH1F("fPionPt","p_{T} distribution",200,0.,10.);
fOutputList->Add(fPionPt);
- fPionEta= new TH1F("fHistQAPionEta","#eta distribution",360,-1.8,1.8);
+ fPionEta= new TH1F("fPionEta","#eta distribution",360,-1.8,1.8);
fOutputList->Add(fPionEta);
- fPionPhi = new TH1F("fHistQAPionPhi","#phi distribution",340,0,6.8);
+ fPionPhi = new TH1F("fPionPhi","#phi distribution",340,0,6.8);
fOutputList->Add(fPionPhi);
- fKaonPt = new TH1F("fHistQAKaonPt","p_{T} distribution",200,0.,10.);
+ fKaonPt = new TH1F("fKaonPt","p_{T} distribution",200,0.,10.);
fOutputList->Add(fKaonPt);
- fKaonEta= new TH1F("fHistQAKaonEta","#eta distribution",360,-1.8,1.8);
+ fKaonEta= new TH1F("fKaonEta","#eta distribution",360,-1.8,1.8);
fOutputList->Add(fKaonEta);
- fKaonPhi = new TH1F("fHistQAKaonPhi","#phi distribution",340,0,6.8);
+ fKaonPhi = new TH1F("fKaonPhi","#phi distribution",340,0,6.8);
fOutputList->Add(fKaonPhi);
- fProtonPt = new TH1F("fHistQAProtonPt","p_{T} distribution",200,0.,10.);
+ fProtonPt = new TH1F("fProtonPt","p_{T} distribution",200,0.,10.);
fOutputList->Add(fProtonPt);
- fProtonEta= new TH1F("fHistQAProtonEta","#eta distribution",360,-1.8,1.8);
+ fProtonEta= new TH1F("fProtonEta","#eta distribution",360,-1.8,1.8);
fOutputList->Add(fProtonEta);
- fProtonPhi= new TH1F("fHistQAProtonPhi","#phi distribution",340,0,6.8);
+ fProtonPhi= new TH1F("fProtonPhi","#phi distribution",340,0,6.8);
fOutputList->Add(fProtonPhi);
}
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]);
}
+ fPriHistShare = new TH1F ("fPriHistShare","Shared clusters, primaries;#shared clusters;counts",160,0,160);
+ fOutput->Add(fPriHistShare);
+
Int_t eventplaneaxis=0;
if (fRequestEventPlane) eventplaneaxis=1;
"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";
+ "multiplicity: 0, 1, 2, 3, 4, 5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100.1\n"
+ "multiplicity_mixing: 0., 1., 2., 3., 4., 5., 10., 20., 30., 40., 50., 60., 70., 80., 90., 100.1\n";
+
+ if(fV0TrigCorr){
+defaultBinningStr += "InvariantMass:0.200,0.300,0.398,0.399,0.4,0.401,0.402,0.403,0.404,0.405,0.406,0.407,0.408,0.409,0.41,0.411,0.412,0.413,0.414,0.415,0.416,0.417,0.418,0.419,0.42,0.421,0.422,0.423,0.424,0.425,0.426,0.427,0.428,0.429,0.43,0.431,0.432,0.433,0.434,0.435,0.436,0.437,0.438,0.439,0.44,0.441,0.442,0.443,0.444,0.445,0.446,0.447,0.448,0.449,0.45,0.451,0.452,0.453,0.454,0.455,0.456,0.457,0.458,0.459,0.46,0.461,0.462,0.463,0.464,0.465,0.466,0.467,0.468,0.469,0.47,0.471,0.472,0.473,0.474,0.475,0.476,0.477,0.478,0.479,0.48,0.481,0.482,0.483,0.484,0.485,0.486,0.487,0.488,0.489,0.49,0.491,0.492,0.493,0.494,0.495,0.496,0.497,0.498,0.499,0.5,0.501,0.502,0.503,0.504,0.505,0.506,0.507,0.508,0.509,0.51,0.511,0.512,0.513,0.514,0.515,0.516,0.517,0.518,0.519,0.52,0.521,0.522,0.523,0.524,0.525,0.526,0.527,0.528,0.529,0.53,0.531,0.532,0.533,0.534,0.535,0.536,0.537,0.538,0.539,0.54,0.541,0.542,0.543,0.544,0.545,0.546,0.547,0.548,0.549,0.55,0.551,0.552,0.553,0.554,0.555,0.556,0.557,0.558,0.559,0.56,0.561,0.562,0.563,0.564,0.565,0.566,0.567,0.568,0.569,0.57,0.571,0.572,0.573,0.574,0.575,0.576,0.577,0.578,0.579,0.58,0.581,0.582,0.583,0.584,0.585,0.586,0.587,0.588,0.589,0.59,0.591,0.592,0.593,0.594,0.595,0.596,0.597,0.598,0.599,0.600,0.700,0.800,0.900,1.000,1.065,1.066,1.067,1.068,1.069,1.07,1.071,1.072,1.073,1.074,1.075,1.076,1.077,1.078,1.079,1.08,1.081,1.082,1.083,1.084,1.085,1.086,1.087,1.088,1.089,1.09,1.091,1.092,1.093,1.094,1.095,1.096,1.097,1.098,1.099,1.1,1.101,1.102,1.103,1.104,1.105,1.106,1.107,1.108,1.109,1.11,1.111,1.112,1.113,1.114,1.115,1.116,1.117,1.118,1.119,1.12,1.121,1.122,1.123,1.124,1.125,1.126,1.127,1.128,1.129,1.13,1.131,1.132,1.133,1.134,1.135,1.136,1.137,1.138,1.139,1.14,1.141,1.142,1.143,1.144,1.145,1.146,1.147,1.148,1.149,1.15,1.151,1.152,1.153,1.154,1.155,1.156,1.157,1.158,1.159,1.16,1.161,1.162,1.163,1.164,1.165\n";
+ }
if(fRequestEventPlane){
defaultBinningStr += "eventPlane: -0.5,0.5,1.5,2.5,3.5\n"; // Event Plane Bins (Psi: -0.5->0.5 (in plane), 0.5->1.5 (intermediate), 1.5->2.5 (out of plane), 2.5->3.5 (rest))
}
-
+ if(fRequestEventPlanemixing){
+defaultBinningStr += "eventPlanemixing: 0.0*TMath::DegToRad(), 30.0*TMath::DegToRad(), 60.0*TMath::DegToRad(), 90.0*TMath::DegToRad(), 120.0*TMath::DegToRad(),150.0*TMath::DegToRad(),180.1*TMath::DegToRad()\n";
+ }
if(fcontainPIDtrig){
- defaultBinningStr += "PIDTrig: -0.5,0.5,1.5,2.5,3.5\n"; // course
+ defaultBinningStr += "PIDTrig: -0.5,0.5,1.5,2.5,3.5,4.5,5.5,6.5,7.5,8.5,9.5,10.5\n"; // course
}
if(fcontainPIDasso){
defaultBinningStr += "PIDAsso: -0.5,0.5,1.5,2.5,3.5\n"; // course
}
if(fcontainPIDtrig && !fcontainPIDasso){
+ if(fV0TrigCorr){
+ dBinsPair[dim_val] = GetBinning(fBinningString, "InvariantMass", iBinPair[dim_val]);
+ axisTitlePair[dim_val] = "InvariantMass";
+ }
+ else{
dBinsPair[dim_val] = GetBinning(fBinningString, "PIDTrig", iBinPair[dim_val]);
axisTitlePair[dim_val] = "PIDTrig";
+ }
if(SetChargeAxis==2){
dBinsPair[dim_val+1] = GetBinning(fBinningString, "TrigCharge", iBinPair[dim_val+1]);
axisTitlePair[dim_val+1] = "TrigCharge";
if(fcontainPIDtrig && fcontainPIDasso){
+ if(fV0TrigCorr){
+ dBinsPair[dim_val] = GetBinning(fBinningString, "InvariantMass", iBinPair[dim_val]);
+ axisTitlePair[dim_val] = "InvariantMass";
+ }
+ else{
dBinsPair[dim_val] = GetBinning(fBinningString, "PIDTrig", iBinPair[dim_val]);
- axisTitlePair[dim_val] = "PIDTrig";
+ axisTitlePair[dim_val] = "PIDTrig";
+ }
dBinsPair[dim_val+1] = GetBinning(fBinningString, "PIDAsso", iBinPair[dim_val+1]);
axisTitlePair[dim_val+1] = "PIDAsso";
Int_t nPteffbin = -1;
Double_t* Pteff = GetBinning(fBinningString, "p_t_eff", nPteffbin);
+ Int_t multmixbin = -1;
+ Double_t* multmix = GetBinning(fBinningString, "multiplicity_mixing", multmixbin);
+
+ //Set the limits from custom binning
fminPtTrig=dBinsPair[2][0];
fmaxPtTrig=dBinsPair[2][iBinPair[2]];
fminPtAsso=dBinsPair[3][0];
fmincentmult=dBinsPair[0][0];
fmaxcentmult=dBinsPair[0][iBinPair[0]];
+ //event pool manager
+Int_t MaxNofEvents=1000;
+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(fRequestEventPlanemixing){
+ // Event plane angle (Psi) bins for event mixing
+
+ Int_t nPsiBins=-1;;
+ Double_t* psibins = GetBinning(fBinningString, "eventPlanemixing", nPsiBins);
+ fPoolMgr = new AliEventPoolManager(MaxNofEvents,fMaxNofMixingTracks,multmixbin,multmix,NofVrtxBins,ZvrtxBins, nPsiBins, psibins);
+ if(psibins) delete [] psibins;
+ }
+
+ else{
+ const Int_t nPsiBinsd=1;
+ Double_t psibinsd[nPsiBinsd+1]={0.0, 2000.0};
+fPoolMgr = new AliEventPoolManager(MaxNofEvents,fMaxNofMixingTracks,multmixbin,multmix,NofVrtxBins,ZvrtxBins, nPsiBinsd, psibinsd);
+
+ }
+fPoolMgr->SetTargetValues(fMaxNofMixingTracks, 0.1, 5);
+
+ if(!fPoolMgr){
+ AliError("Event Mixing required, but Pool Manager not initialized...");
+ return;
+ }
+
//fminPtComboeff=fminPtTrig;***then this value will be fixed ,even Setter can't change it's value
//fmaxPtComboeff=fmaxPtTrig;
//THnSparses for calculation of efficiency
//ThnSparse for Correlation plots(truth MC)
- if((fAnalysisType == "MCAOD") && ffilltrigIDassoIDMCTRUTH) {//remember that in this case uidentified means other than pions, kaons, protons
+ if(ffilltrigIDassoIDMCTRUTH) {//remember that in this case uidentified means other than pions, kaons, protons
fCorrelatonTruthPrimary = new AliTHn("fCorrelatonTruthPrimary", title, anaSteps, kTrackVariablesPair, iBinPair);
for (Int_t j=0; j<kTrackVariablesPair; j++) {
fOutput->Add(fTHnTrigcount);
}
- if((fAnalysisType =="MCAOD") && ffilltrigIDassoIDMCTRUTH) {
+ if(ffilltrigIDassoIDMCTRUTH) {
//AliTHns for trigger counting(truth MC)
fTHnTrigcountMCTruthPrim = new AliTHn("fTHnTrigcountMCTruthPrim", "fTHnTrigcountMCTruthPrim", 2, dims, fBinst); //2 steps;;;;0->same event;;;;;1->mixed event
for(Int_t i=0; i<dims;i++){
- fTHnTrigcount->SetBinLimits(i, dBinsTrig[i]);
- fTHnTrigcount->SetVarTitle(i, axisTitleTrig[i]);
+ fTHnTrigcountMCTruthPrim->SetBinLimits(i, dBinsTrig[i]);
+ fTHnTrigcountMCTruthPrim->SetVarTitle(i, axisTitleTrig[i]);
}
fOutput->Add(fTHnTrigcountMCTruthPrim);
}
-if(fAnalysisType=="MCAOD"){
+if(fAnalysisType=="MCAOD" || fAnalysisType=="MC"){
if(ffillhistQATruth)
{
MCtruthpt=new TH1F ("MCtruthpt","ptdistributiontruthprim",100,0.,10.);
fhistJetTrigestimate=new TH2F("fhistJetTrigestimate","fhistJetTrigestimate",iBinPair[0],dBinsPair[0],6,-0.5,5.5);
fOutput->Add(fhistJetTrigestimate);
+ fTwoTrackDistancePtdip = new TH3F("fTwoTrackDistancePtdip", ";#Delta#eta;#Delta#varphi;#Delta p_{T}", 36, -1.8, 1.8, 72,-TMath::Pi()/2, 3*TMath::Pi()/2, 40, 0, 10);
+ fOutput->Add(fTwoTrackDistancePtdip);
+
+fTwoTrackDistancePtdipmix = new TH3F("fTwoTrackDistancePtdipmix", ";#Delta#eta;#Delta#varphi;#Delta p_{T}", 36, -1.8, 1.8, 72,-TMath::Pi()/2, 3*TMath::Pi()/2, 40, 0, 10);
+ fOutput->Add(fTwoTrackDistancePtdipmix);
+ TString Histttrname;
+for(Int_t jj=0;jj<2;jj++)// PID type binning
+ {
+ Histttrname="fTwoTrackDistancePt";Histttrname+=jj;
+ fTwoTrackDistancePt[jj] = new TH3F(Histttrname.Data(), ";#Delta#eta;#Delta#varphi^{*}_{min};#Delta p_{T}", 100, -0.15, 0.15, 100, -0.05, 0.05, 20, 0, 10);
+ fOutput->Add(fTwoTrackDistancePt[jj]);
+
+ Histttrname="fTwoTrackDistancePtmix";Histttrname+=jj;
+ fTwoTrackDistancePtmix[jj] = new TH3F(Histttrname.Data(), ";#Delta#eta;#Delta#varphi^{*}_{min};#Delta p_{T}", 100, -0.15, 0.15, 100, -0.05, 0.05, 20, 0, 10);
+ fOutput->Add(fTwoTrackDistancePtmix[jj]);
+ }
//Mixing
-DefineEventPool();
+//DefineEventPool();
if(fapplyTrigefficiency || fapplyAssoefficiency)
{
fileT->Close();
}
+
+ delete [] EtaBin;
+ delete [] Pteff;
+ delete [] multmix;
+
+ //******************************************************************V0 plots*********************************************//
+ if(fV0TrigCorr){
+ //histos for v0
+ //Double_t BinsV0[]={1.0,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0,2.2,2.4,2.6,2.8,3.0,3.2,3.4,3.6,3.8,4.0,4.5,5.0,5.5,6.5,8.0};
+ //const Int_t nbinsV0 =sizeof(BinsV0)/sizeof(Double_t)-1;
+
+
+
+ fHistRawPtCentInvK0s= new TH3F("fHistRawPtCentInvK0s", "K^{0}_{s}: mass vs #it{p}_{T};Mass (GeV/#it{c}^2);#it{p}_{T} (GeV/#it{c});centrality",100,0.398,0.598,100,0.0,10.0,100,0.0,100.);
+fOutput->Add(fHistRawPtCentInvK0s);
+
+
+ fHistRawPtCentInvLambda= new TH3F("fHistRawPtCentInvLambda", "#Lambda: mass vs #it{p}_{T};Mass (GeV/#it{c}^2);#it{p}_{T} (GeV/#it{c});centrality",100,1.065,1.165,100,0.0,10.0,100,0.0,100.);
+fOutput->Add(fHistRawPtCentInvLambda);
+
+
+ fHistRawPtCentInvAntiLambda= new TH3F("fHistRawPtCentInvAntiLambda", "#bar{#Lambda} : mass vs #it{p}_{T};Mass (GeV/#it{c}^2);#it{p}_{T} (GeV/#it{c});centrality",100,1.065,1.165,100,0.0,10.0,100,0.0,100.);
+fOutput->Add(fHistRawPtCentInvAntiLambda);
+
+
+ fHistFinalPtCentInvK0s= new TH3F("fHistFinalPtCentInvK0s", "K^{0}_{s}: mass vs #it{p}_{T};Mass (GeV/#it{c}^2);#it{p}_{T} (GeV/#it{c});centrality",100,0.398,0.598,100,0.0,10.0,100,0.0,100.);
+fOutput->Add(fHistFinalPtCentInvK0s);
+
+
+ fHistFinalPtCentInvLambda= new TH3F("fHistFinalPtCentInvLambda", "#Lambda: mass vs #it{p}_{T};Mass (GeV/#it{c}^2);#it{p}_{T} (GeV/#it{c});centrality",100,1.065,1.165,100,0.0,10.0,100,0.0,100.);
+fOutput->Add(fHistFinalPtCentInvLambda);
+
+
+ fHistFinalPtCentInvAntiLambda= new TH3F("fHistFinalPtCentInvAntiLambda", "#bar{#Lambda} : mass vs #it{p}_{T};Mass (GeV/#it{c}^2);#it{p}_{T} (GeV/#it{c});centrality",100,1.065,1.165,100,0.0,10.0,100,0.0,100.);
+fOutput->Add(fHistFinalPtCentInvAntiLambda);
+
+ }
+
+ //*************************************************************EP plots***********************************************//
+ if(fRequestEventPlane){
+ // TProfile for resolutions 3 subevents (V0A, V0C, TPC)
+ // v2
+ fHResTPCv0A2 = new TProfile("hResTPCv0A2","",nCentrBin,0,nCentrBin);
+ fHResTPCv0C2 = new TProfile("hResTPCv0C2","",nCentrBin,0,nCentrBin);
+ fHResv0Cv0A2 = new TProfile("hResv0Cv0A2","",nCentrBin,0,nCentrBin);
+
+ fList->Add(fHResTPCv0A2);
+ fList->Add(fHResTPCv0C2);
+ fList->Add(fHResv0Cv0A2);
+
+ // v3
+ fHResTPCv0A3 = new TProfile("hResTPCv0A3","",nCentrBin,0,nCentrBin);
+ fHResTPCv0C3 = new TProfile("hResTPCv0C3","",nCentrBin,0,nCentrBin);
+ fHResv0Cv0A3 = new TProfile("hResv0Cv0A3","",nCentrBin,0,nCentrBin);
+
+ fList->Add(fHResTPCv0A3);
+ fList->Add(fHResTPCv0C3);
+ fList->Add(fHResv0Cv0A3);
+
+ // MC as in the dataEP resolution (but using MC tracks)
+ if(fAnalysisType == "MCAOD" && fV2){
+ fHResMA2 = new TProfile("hResMA2","",nCentrBin,0,nCentrBin);
+ fHResMC2 = new TProfile("hResMC2","",nCentrBin,0,nCentrBin);
+ fHResAC2 = new TProfile("hResAC2","",nCentrBin,0,nCentrBin);
+ fList->Add(fHResMA2);
+ fList->Add(fHResMC2);
+ fList->Add(fHResAC2);
+ }
+ if(fAnalysisType == "MCAOD" && fV3){
+ fHResMA3 = new TProfile("hResMA3","",nCentrBin,0,nCentrBin);
+ fHResMC3 = new TProfile("hResMC3","",nCentrBin,0,nCentrBin);
+ fHResAC3 = new TProfile("hResAC3","",nCentrBin,0,nCentrBin);
+ fList->Add(fHResMA3);
+ fList->Add(fHResMC3);
+ fList->Add(fHResAC3);
+ }
+
+
+ // V0A and V0C event plane distributions
+ //v2
+ fPhiRPTPC = new TH2F("fPhiRPTPCv2","#phi distribution of EP TPC;centrality;#phi (rad)",nCentrBin,0,nCentrBin,nPsiTOF,-TMath::Pi()/2,TMath::Pi()/2);
+ fPhiRPTPCv3 = new TH2F("fPhiRPTPCv3","#phi distribution of EP TPC;centrality;#phi (rad)",nCentrBin,0,nCentrBin,nPsiTOF,-TMath::Pi()/3,TMath::Pi()/3);
+ fList->Add(fPhiRPTPC);
+ fList->Add(fPhiRPTPCv3);
+
+ fPhiRPv0A = new TH2F("fPhiRPv0Av2","#phi distribution of EP VZERO-A;centrality;#phi (rad)",nCentrBin,0,nCentrBin,nPsiTOF,-TMath::Pi()/2,TMath::Pi()/2);
+ fPhiRPv0C = new TH2F("fPhiRPv0Cv2","#phi distribution of EP VZERO-C;centrality;#phi (rad)",nCentrBin,0,nCentrBin,nPsiTOF,-TMath::Pi()/2,TMath::Pi()/2);
+ fList->Add(fPhiRPv0A);
+ fList->Add(fPhiRPv0C);
+
+ //v3
+ fPhiRPv0Av3 = new TH2F("fPhiRPv0Av3","#phi distribution of EP VZERO-A;centrality;#phi (rad)",nCentrBin,0,nCentrBin,nPsiTOF,-TMath::Pi()/3,TMath::Pi()/3);
+ fPhiRPv0Cv3 = new TH2F("fPhiRPv0Cv3","#phi distribution of EP VZERO-C;centrality;#phi (rad)",nCentrBin,0,nCentrBin,nPsiTOF,-TMath::Pi()/3,TMath::Pi()/3);
+ fList->Add(fPhiRPv0Av3);
+ fList->Add(fPhiRPv0Cv3);
+
+ fHistEventPlaneTruth = new TH2F("fHistEventPlaneTruth","#phi distribution of EP MCTRUTHheader;centrality;#phi (rad)",nCentrBin,0,nCentrBin,nPsiTOF,-TMath::Pi()/2,TMath::Pi()/2);
+ fList->Add(fHistEventPlaneTruth);
+
+ }
//*****************************************************PIDQA histos*****************************************************//
}
//nsigmaDC plot
- if(fUseExclusiveNSigma) {
for(Int_t ipart=0;ipart<NSpecies;ipart++){
for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
Double_t miny=-10;
fOutputList->Add(fHistoNSigma);
}
}
- }
+
//nsigmaMC plot
if (fAnalysisType == "MCAOD"){
for(Int_t ipart=0;ipart<NSpecies;ipart++){
PostData(1, fOutput); // Post data for ALL output slots >0 here, to get at least an empty histogram
PostData(2, fOutputList);
-
+ if(fRequestEventPlane) PostData(3, fList);
AliInfo("Finished setting up the Output");
- TH1::AddDirectory(oldStatus);
-
-
-
+ TH1::AddDirectory(oldStatus);
}
//-------------------------------------------------------------------------------
void AliTwoParticlePIDCorr::UserExec( Option_t * ){
-
if(fAnalysisType == "AOD") {
}//AOD--analysis-----
- else if(fAnalysisType == "MCAOD") {
+ else if(fAnalysisType == "MCAOD" || fAnalysisType == "MC") {
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");
+
+ // get the event (for generator level: MCEvent())
+ AliVEvent* event = NULL;
+ if(fAnalysisType == "MC") {
+ event = dynamic_cast<AliVEvent*>(MCEvent());
+ }
+ else{
+ event = dynamic_cast<AliVEvent*>(InputEvent());
+ }
+ if(!event) {
+ AliError("eventMain not available");
return;
}
-
-// count all events(physics triggered)
- fEventCounter->Fill(1);
+
+ Double_t Inv_mass=0.0;//has no meaning for pions, kaons and protons(just set 0.0) to fill the LRCParticlePID position
+ evplaneMC=999.;
+ fgPsi2v0aMC=999.;
+ fgPsi2v0cMC=999.;
+ fgPsi2tpcMC=999.;
+ fgPsi3v0aMC=999.;
+ fgPsi3v0cMC=999.;
+ fgPsi3tpcMC=999.;
+ gReactionPlane = 999.;
// get centrality object and check quality(valid for p-Pb and Pb-Pb; coming soon for pp 7 TeV)
Double_t cent_v0=-1.0;
Double_t effcent=1.0;
Double_t refmultReco =0.0;
- Double_t gReactionPlane = -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
+ 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)
- 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(fAnalysisType == "MC"){
+
+ AliMCEvent *gMCEvent = dynamic_cast<AliMCEvent*>(event);
- 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;
+ if(!gMCEvent) {
+ AliError("mcEvent not available");
+ 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));
- }
+// count all events(physics triggered)
+ fEventCounter->Fill(1);
- if (fSampleType=="pp" && fCentralityMethod.EndsWith("_MANUAL"))
- {
- //make the event selection with reco vertex cut and centrality cut and return the value of the centrality
- Double_t refmultTruth = GetAcceptedEventMultiplicity(aod,kTRUE); //incase of ref multiplicity it will return the truth MC ref mullt value; need to determine the ref mult value separately for reco Mc case; in case of centrality this is final and fine
- refmultReco = GetAcceptedEventMultiplicity(aod,kFALSE);
- if(refmultTruth<=0 || refmultReco<=0) return;
- cent_v0=refmultTruth;
- }
- else {
- cent_v0=GetAcceptedEventMultiplicity(aod,kTRUE); //centrality value; 2nd argument has no meaning
- if(cent_v0<0.) return;
- }
+ AliGenEventHeader *header = dynamic_cast<AliGenEventHeader*>(gMCEvent->GenEventHeader());
+ if(!header) return;
+ TArrayF gVertexArray;
+ header->PrimaryVertex(gVertexArray);
+ Float_t zVtxmc =gVertexArray.At(2);
+ //cout<<"*****************************************************************************************************hi I am here"<<endl;
- effcent=cent_v0;// This will be required for efficiency THn filling(specially in case of pp)
+
+ cent_v0=GetAcceptedEventMultiplicity((AliVEvent*)gMCEvent,kFALSE); //b value; 2nd argument has no meaning
+
+ if(cent_v0<0.) return;//mainly returns impact parameter
- //get the event plane in case of PbPb
- if(fSampleType=="PbPb"){
+ //get the event plane in case of PbPb
if(fRequestEventPlane){
- gReactionPlane = GetEventPlane(aod,kTRUE);//get the truth event plane
- fHistEventPlaneTruth->Fill(gReactionPlane,cent_v0);
- if(gReactionPlane < 0) return;
+ gReactionPlane=GetEventPlane((AliVEvent*)gMCEvent,kTRUE,cent_v0);//get the truth event plane,middle argument has no meaning in this case
+ if(gReactionPlane==999.) return;
}
- }
- 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)
- //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;
- }
+for (Int_t iTracks = 0; iTracks < gMCEvent->GetNumberOfPrimaries(); iTracks++) {
+ AliMCParticle* partMC = dynamic_cast<AliMCParticle *>(gMCEvent->GetTrack(iTracks));
+ if (!partMC) {
+ AliError(Form("Could not receive particle %d", iTracks));
+ continue;
+ }
+//exclude non stable particles
+ if(fselectprimaryTruth && !(gMCEvent->IsPhysicalPrimary(iTracks))) 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
+//give only kinematic cuts at the generator level
if (partMC->Eta() < fmineta || partMC->Eta() > fmaxeta) continue;
if (partMC->Pt() < fminPt || partMC->Pt() > fmaxPt) continue;
if(!partMC) continue;//for safety
+ TParticle *particle = partMC->Particle();
+ if(!particle) continue;
+ Int_t particletypeTruth=-999;
+
+ Int_t pdgtruth = particle->GetPdgCode();
+
//To determine multiplicity in case of PP
nooftrackstruth++;
//cout<<"**************************************"<<TMath::Abs(partMC->GetLabel())<<endl;
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(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" && fCentralityMethod.EndsWith("_MANUAL")) effcent=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
+ if(fRequestEventPlane){
+ FillPIDEventPlane(cent_v0,particletypeTruth,partMC->Phi(),gReactionPlane);
+ }
+ /*
+//Exclude resonances
+ if(fExcludeResonancesInMC) {
+ TParticle *particle = track->Particle();
+ if(!particle) continue;
+
+ Bool_t kExcludeParticle = kFALSE;
+ Int_t gMotherIndex = particle->GetFirstMother();
+ if(gMotherIndex != -1) {
+ AliMCParticle* motherTrack = dynamic_cast<AliMCParticle *>(event->GetTrack(gMotherIndex));
+ if(motherTrack) {
+ TParticle *motherParticle = motherTrack->Particle();
+ if(motherParticle) {
+ Int_t pdgCodeOfMother = motherParticle->GetPdgCode();
+ //if((pdgCodeOfMother == 113)||(pdgCodeOfMother == 213)||(pdgCodeOfMother == 221)||(pdgCodeOfMother == 223)||(pdgCodeOfMother == 331)||(pdgCodeOfMother == 333)) {
+ }
+ if(pdgCodeOfMother == 113 // rho0
+ || pdgCodeOfMother == 213 || pdgCodeOfMother == -213 // rho+
+ // || pdgCodeOfMother == 221 // eta
+ // || pdgCodeOfMother == 331 // eta'
+ // || pdgCodeOfMother == 223 // omega
+ // || pdgCodeOfMother == 333 // phi
+ || pdgCodeOfMother == 311 || pdgCodeOfMother == -311 // K0
+ // || pdgCodeOfMother == 313 || pdgCodeOfMother == -313 // K0*
+ // || pdgCodeOfMother == 323 || pdgCodeOfMother == -323 // K+*
+ || pdgCodeOfMother == 3122 || pdgCodeOfMother == -3122 // Lambda
+ || pdgCodeOfMother == 111 // pi0 Dalitz
+ ) {
+ kExcludeParticle = kTRUE;
+ }
+ }
+ }
+ }
+
+ //Exclude from the analysis decay products of rho0, rho+, eta, eta' and phi
+ if(kExcludeParticle) continue;
+ }
+
+ //Exclude electrons with PDG
+ if(fExcludeElectronsInMC) {
+
+ TParticle *particle = track->Particle();
+
+ if (particle){
+ if(TMath::Abs(particle->GetPdgCode()) == 11) continue;
+ }
+ }
+ */
- Double_t primmctruth[4] = {effcent, zVtxmc,partMC->Pt(), partMC->Eta()};
- if(ffillefficiency)
- {
- fTrackHistEfficiency[5]->Fill(primmctruth,0);//for all primary truth particles(4)
- if (TMath::Abs(pdgtruth)==211 || TMath::Abs(pdgtruth)==321) fTrackHistEfficiency[3]->Fill(primmctruth,0);//for primary truth mesons(3)
- if (TMath::Abs(pdgtruth)==2212 || TMath::Abs(pdgtruth)==321) fTrackHistEfficiency[4]->Fill(primmctruth,0);//for primary truth kaons+protons(4)
- if (TMath::Abs(pdgtruth)==211) fTrackHistEfficiency[0]->Fill(primmctruth,0);//for pions
- if (TMath::Abs(pdgtruth)==321) fTrackHistEfficiency[1]->Fill(primmctruth,0);//for kaons
- if(TMath::Abs(pdgtruth)==2212) fTrackHistEfficiency[2]->Fill(primmctruth,0);//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()<=fmaxPtAsso) || (partMC->Pt()>=fminPtTrig && partMC->Pt()<=fmaxPtTrig))//to reduce memory consumption in pool
{
if(partMC->Charge()>0) chargeval=1;
if(partMC->Charge()<0) chargeval=-1;
if(chargeval==0) continue;
-LRCParticlePID* copy6 = new LRCParticlePID(particletypeTruth,chargeval,partMC->Pt(),partMC->Eta(), partMC->Phi(),effmatrixtruth);
+ const TBits *clustermap=0;
+ const TBits *sharemap=0;
+ LRCParticlePID* copy6 = new LRCParticlePID(particletypeTruth,Inv_mass,chargeval,partMC->Pt(),partMC->Eta(), partMC->Phi(),effmatrixtruth,clustermap,sharemap);
//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
+ }//track loop ends
-//*********************still in event loop
+ if(nooftrackstruth>0.0){
+if (fSampleType=="pPb" || fSampleType=="PbPb" || fPPVsMultUtils==kTRUE || fCentralityMethod == "MC_b") fCentralityCorrelation->Fill(cent_v0, nooftrackstruth);//only with unidentified tracks(i.e before PID selection);;;;;can be used to remove centrality outliers??????
- if (fSampleType=="PbPb"){
- if (fRandomizeReactionPlane)//only for TRuth MC??
+AliCollisionGeometry* collGeometry = dynamic_cast<AliCollisionGeometry*>(gMCEvent->GenEventHeader());
+ if(collGeometry){
+Float_t NpartProj= collGeometry-> ProjectileParticipants();
+Float_t NpartTarg = collGeometry->TargetParticipants();
+ Float_t Npart= (NpartProj + NpartTarg);
+fNchNpartCorr->Fill(Npart,nooftrackstruth);
+ }
+ }
+
+ if (fRandomizeReactionPlane)//only for TRuth MC??
{
Double_t centralityDigits = cent_v0*1000. - (Int_t)(cent_v0*1000.);
Double_t angle = TMath::TwoPi() * centralityDigits;
AliInfo(Form("Shifting phi of all tracks by %f (digits %f)", angle, centralityDigits));
ShiftTracks(tracksMCtruth, angle);
}
- }
+
Float_t weghtval=1.0;
+ Float_t bSign = 0;
if(nooftrackstruth>0.0 && ffilltrigIDassoIDMCTRUTH)
{
Fillcorrelation(gReactionPlane,tracksMCtruth,0,cent_v0,zVtxmc,weghtval,kFALSE,bSign,fPtOrderMCTruth,kFALSE,kFALSE,"trigIDassoIDMCTRUTH");//mixcase=kFALSE for same event case
//start mixing
-AliEventPool* pool2 = fPoolMgr->GetEventPool(cent_v0, zVtxmc+200);
+ AliEventPool* pool2 = fPoolMgr->GetEventPool(cent_v0, zVtxmc+200, gReactionPlane);
if (pool2 && pool2->IsReady())
{//start mixing only when pool->IsReady
if(tracksMCtruth && tracksMCtruth->GetEntriesFast()>0)
if(tracksMCtruth) delete tracksMCtruth;
-//now deal with reco tracks
-//detrmine the ref mult in case of Reco(not required if we get centrality info from AliCentrality)
- if (fSampleType=="pp" && fCentralityMethod.EndsWith("_MANUAL")) cent_v0=refmultReco;
- effcent=cent_v0;// This will be required for efficiency THn filling(specially in case of pp)
+ }//MC type
- if(fSampleType=="PbPb"){
- if(fRequestEventPlane){
- gReactionPlane = GetEventPlane(aod,kFALSE);//get the reconstructed event plane
- fHistEventPlaneReco->Fill(gReactionPlane,cent_v0);
- if(gReactionPlane < 0) return;
- }
+//"MC" type analysis is finished but still in event loop
+
+ else{//if(fAnalysisType=="MCAOD")
+
+ AliAODEvent* aod = dynamic_cast<AliAODEvent*>(event);
+ if (!aod) {
+ AliError("Cannot get the AOD event");
+ return;
}
- //TObjArray* tracksUNID=0;
- TObjArray* tracksUNID = new TObjArray;
- tracksUNID->SetOwner(kTRUE);
+
+// count all events(physics triggered)
+ fEventCounter->Fill(1);
- //TObjArray* tracksID=0;
- TObjArray* tracksID = new TObjArray;
- tracksID->SetOwner(kTRUE);
+
+//check the PIDResponse handler
+ fPID = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetPIDResponse();
+ if (!fPID) AliFatal("This Task needs the PID response attached to the inputHandler");
+ //if (!fPID) return;
+// get mag. field required for twotrack efficiency cut
+ Float_t bSign = 0;
+ bSign = (aod->GetMagneticField() > 0) ? 1 : -1;
- Float_t bSign1=aod->GetHeader()->GetMagneticField() ;//used for reconstructed track dca cut
+ //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
- 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
+ if (fInjectedSignals)
+ {
+ AliGenEventHeader* eventHeader = 0;
+ Int_t headers = 0;
-//remove injected signals
- if(fInjectedSignals)
- {
- AliAODMCParticle* mother = recomatched;
+// AOD
+ if (!header)
+ AliFatal("fInjectedSignals set but no MC header found");
+
+ headers = header->GetNCocktailHeaders();
+ eventHeader = header->GetCocktailHeader(0);
- 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)
+ if (!eventHeader)
{
- Printf("WARNING: No mother found for particle %d:", recomatched->GetLabel());
- continue;
+ // 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;
}
- 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;
+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));
+ }
-if (track->GetLabel() == track2->GetLabel())
+ if (fSampleType=="pp_2_76" || fCentralityMethod.EndsWith("_MANUAL") || (fSampleType=="pp_7" && fPPVsMultUtils==kFALSE))
{
-isduplicate2=kTRUE;
- break;
- }
- }//2nd loop ends
+ //make the event selection with reco vertex cut and centrality cut and return the value of the centrality
+ Double_t refmultTruth = GetAcceptedEventMultiplicity((AliVEvent*)aod,kTRUE); //incase of ref multiplicity it will return the truth MC ref mullt value; need to determine the ref mult value separately for reco Mc case; in case of centrality this is final and fine
+ refmultReco = GetAcceptedEventMultiplicity((AliVEvent*)aod,kFALSE);
+ if(refmultTruth<=0 || refmultReco<=0) return;
+ cent_v0=refmultTruth;
}
- if(fRemoveDuplicates && isduplicate2) continue;//remove duplicates
-
- fHistQA[11]->Fill(track->GetTPCNcls());
- Int_t tracktype=ClassifyTrack(track,trkVtx,bSign1,kTRUE);//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;
+ else {
+ cent_v0=GetAcceptedEventMultiplicity((AliVEvent*)aod,kFALSE); //centrality value; 2nd argument has no meaning
+ }
-//tag all particles as unidentified
- particletypeMC=unidentified;
+ if(cent_v0<0.) return;
+ effcent=cent_v0;// This will be required for efficiency THn filling(specially in case of pp)
- Float_t effmatrix=1.;
+ //get the event plane in case of PbPb
+ if(fRequestEventPlane){
+ gReactionPlane=GetEventPlane((AliVEvent*)aod,kTRUE,cent_v0);//get the truth event plane
+ if(gReactionPlane==999.) return;
+ }
+
+ //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!
-// -- Fill THnSparse for efficiency calculation
- if (fSampleType=="pp" && fCentralityMethod.EndsWith("_MANUAL")) effcent=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
- //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)
+ eventno++;
- //Clone & Reduce track list(TObjArray) for unidentified particles
-if((track->Pt()>=fminPtAsso && track->Pt()<=fmaxPtAsso) || (track->Pt()>=fminPtTrig && track->Pt()<=fmaxPtTrig))//to reduce memory consumption in pool
- {
- Short_t chargeval=0;
- if(track->Charge()>0) chargeval=1;
- if(track->Charge()<0) chargeval=-1;
- if(chargeval==0) continue;
- if (fapplyTrigefficiency || fapplyAssoefficiency)//get the trackingefficiency x contamination factor for unidentified particles
- effmatrix=GetTrackbyTrackeffvalue(track,effcent,zvtx,particletypeMC);
- LRCParticlePID* copy = new LRCParticlePID(particletypeMC,chargeval,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)
- }
+ //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;
+ }
-//get the pdg code of the corresponding truth particle
- Int_t pdgCode = ((AliAODMCParticle*)recomatched)->GetPdgCode();
+//consider only charged particles
+ if(partMC->Charge() == 0) continue;
- Double_t allrecomatchedpid[4] = {effcent, zVtxmc,recomatched->Pt(), recomatched->Eta()};
- if(ffillefficiency) {
-fTrackHistEfficiency[5]->Fill(allrecomatchedpid,2);//for allreco matched
- if(TMath::Abs(pdgCode)==211 || TMath::Abs(pdgCode)==321) fTrackHistEfficiency[3]->Fill(allrecomatchedpid,2);//for mesons
- if(TMath::Abs(pdgCode)==321 || TMath::Abs(pdgCode)==2212) fTrackHistEfficiency[4]->Fill(allrecomatchedpid,2);//for kaons+protons
- if(TMath::Abs(pdgCode)==211) fTrackHistEfficiency[0]->Fill(allrecomatchedpid,2);//for pions
- if(TMath::Abs(pdgCode)==321) fTrackHistEfficiency[1]->Fill(allrecomatchedpid,2);//for kaons
- if(TMath::Abs(pdgCode)==2212) fTrackHistEfficiency[2]->Fill(allrecomatchedpid,2);//for protons
+//consider only primary particles; neglect all secondary particles including from weak decays
+ if(fselectprimaryTruth && !partMC->IsPhysicalPrimary()) continue;
- if (((AliAODMCParticle*)recomatched)->IsPhysicalPrimary()) {
- fTrackHistEfficiency[5]->Fill(allrecomatchedpid,1);//for primreco matched
- if(TMath::Abs(pdgCode)==211 || TMath::Abs(pdgCode)==321) fTrackHistEfficiency[3]->Fill(allrecomatchedpid,1);//for mesons
- if(TMath::Abs(pdgCode)==321 || TMath::Abs(pdgCode)==2212) fTrackHistEfficiency[4]->Fill(allrecomatchedpid,1);//for kaons+protons
- if( TMath::Abs(pdgCode)==211) fTrackHistEfficiency[0]->Fill(allrecomatchedpid,1);//for pions
- if( TMath::Abs(pdgCode)==321) fTrackHistEfficiency[1]->Fill(allrecomatchedpid,1);//for kaons
- if( TMath::Abs(pdgCode)==2212) fTrackHistEfficiency[2]->Fill(allrecomatchedpid,1);//for protons
- }
- }
- //now start the particle identification process:)
+//remove injected signals(primaries above <maxLabel>)
+ if (fInjectedSignals && partMC->GetLabel() >= skipParticlesAbove) continue;
-Float_t dEdx = track->GetTPCsignal();
- fHistoTPCdEdx->Fill(track->Pt(), dEdx);
+//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
- if(HasTOFPID(track))
-{
-Double_t beta = GetBeta(track);
-fHistoTOFbeta->Fill(track->Pt(), beta);
- }
+//give only kinematic cuts at the truth level
+ if (partMC->Eta() < fmineta || partMC->Eta() > fmaxeta) continue;
+ if (partMC->Pt() < fminPt || partMC->Pt() > fmaxPt) continue;
-//do track identification(nsigma method)
- particletypeMC=GetParticle(track,fFIllPIDQAHistos);//******************************problem is here
+ if(!partMC) continue;//for safety
-switch(TMath::Abs(pdgCode)){
- case 2212:
- if(fFIllPIDQAHistos){
- for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
- if((ipid!=NSigmaTPC) && (!HasTOFPID(track)))continue;//not filling TOF and combined if no TOF PID
- TH2F *h=GetHistogram2D(Form("NSigmaMC_%d_%d",SpProton,ipid));
- h->Fill(track->Pt(),fnsigmas[SpProton][ipid]);
+ //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)
+
+ if(fRequestEventPlane){
+ FillPIDEventPlane(cent_v0,particletypeTruth,partMC->Phi(),gReactionPlane);
+ }
+
+ // -- Fill THnSparse for efficiency and contamination calculation
+ if (fSampleType=="pp_2_76" || fCentralityMethod.EndsWith("_MANUAL") || (fSampleType=="pp_7" && fPPVsMultUtils==kFALSE)) effcent=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] = {effcent, zVtxmc,partMC->Pt(), partMC->Eta()};
+ if(ffillefficiency)
+ {
+ fTrackHistEfficiency[5]->Fill(primmctruth,0);//for all primary truth particles(4)
+ if (TMath::Abs(pdgtruth)==211 || TMath::Abs(pdgtruth)==321) fTrackHistEfficiency[3]->Fill(primmctruth,0);//for primary truth mesons(3)
+ if (TMath::Abs(pdgtruth)==2212 || TMath::Abs(pdgtruth)==321) fTrackHistEfficiency[4]->Fill(primmctruth,0);//for primary truth kaons+protons(4)
+ if (TMath::Abs(pdgtruth)==211) fTrackHistEfficiency[0]->Fill(primmctruth,0);//for pions
+ if (TMath::Abs(pdgtruth)==321) fTrackHistEfficiency[1]->Fill(primmctruth,0);//for kaons
+ if (TMath::Abs(pdgtruth)==2212) fTrackHistEfficiency[2]->Fill(primmctruth,0);//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()<=fmaxPtAsso) || (partMC->Pt()>=fminPtTrig && partMC->Pt()<=fmaxPtTrig))//to reduce memory consumption in pool
+ {
+ Short_t chargeval=0;
+ if(partMC->Charge()>0) chargeval=1;
+ if(partMC->Charge()<0) chargeval=-1;
+ if(chargeval==0) continue;
+ const TBits *clustermap=0;
+ const TBits *sharemap=0;
+ LRCParticlePID* copy6 = new LRCParticlePID(particletypeTruth,Inv_mass,chargeval,partMC->Pt(),partMC->Eta(), partMC->Phi(),effmatrixtruth,clustermap,sharemap);
+//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 (fRandomizeReactionPlane)//only for TRuth MC??
+ {
+ Double_t centralityDigits = cent_v0*1000. - (Int_t)(cent_v0*1000.);
+ Double_t angle = TMath::TwoPi() * centralityDigits;
+ AliInfo(Form("Shifting phi of all tracks by %f (digits %f)", angle, centralityDigits));
+ ShiftTracks(tracksMCtruth, angle);
+ }
+
+if(nooftrackstruth>0.0){
+if (fSampleType=="pPb" || fSampleType=="PbPb" || fPPVsMultUtils==kTRUE || fCentralityMethod == "MC_b") fCentralityCorrelation->Fill(cent_v0, nooftrackstruth);//only with unidentified tracks(i.e before PID selection);;;;;can be used to remove centrality outliers??????
+
+ AliGenEventHeader* eventHeader = header->GetCocktailHeader(0); // get first MC header from either ESD/AOD (including cocktail header if available)
+ if (eventHeader)
+ {
+AliCollisionGeometry* collGeometry = dynamic_cast <AliCollisionGeometry*> (eventHeader);
+ if(collGeometry){
+Float_t NpartProj= collGeometry-> ProjectileParticipants();
+Float_t NpartTarg = collGeometry->TargetParticipants();
+Float_t Npart= (NpartProj + NpartTarg);
+fNchNpartCorr->Fill(Npart,nooftrackstruth);
+ }
+ }
+ }
+
+ Float_t weghtval=1.0;
+
+if(nooftrackstruth>0.0 && ffilltrigIDassoIDMCTRUTH)
+ {
+ //Fill Correlations for MC truth particles(same event)
+if(tracksMCtruth && tracksMCtruth->GetEntriesFast()>0)//hadron triggered correlation
+ Fillcorrelation(gReactionPlane,tracksMCtruth,0,cent_v0,zVtxmc,weghtval,kFALSE,bSign,fPtOrderMCTruth,kFALSE,kFALSE,"trigIDassoIDMCTRUTH");//mixcase=kFALSE for same event case
+
+//start mixing
+ AliEventPool* pool2 = fPoolMgr->GetEventPool(cent_v0, zVtxmc+200, gReactionPlane);
+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
+ Float_t nmix=(Float_t)pool2->GetCurrentNEvents();
+for (Int_t jMix=0; jMix<pool2->GetCurrentNEvents(); jMix++)
+ { //pool event loop start
+ TObjArray* bgTracks6 = pool2->GetEvent(jMix);
+ if(!bgTracks6) continue;
+ Fillcorrelation(gReactionPlane,tracksMCtruth,bgTracks6,cent_v0,zVtxmc,nmix,(jMix == 0),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
+
+
+//Float_t bSign1=aod->GetHeader()->GetMagneticField() ;//used for reconstructed track dca cut
+ Float_t bSign1=aod->GetMagneticField() ;//used for reconstructed track dca cut
+
+
+//detrmine the ref mult in case of Reco(not required if we get centrality info from AliCentrality)
+ if (fSampleType=="pp_2_76" || fCentralityMethod.EndsWith("_MANUAL") || (fSampleType=="pp_7" && fPPVsMultUtils==kFALSE)) cent_v0=refmultReco;
+ effcent=cent_v0;// This will be required for efficiency THn filling(specially in case of pp)
+
+ if(fRequestEventPlane){
+ gReactionPlane = GetEventPlane((AliVEvent*)aod,kFALSE,cent_v0);//get the reconstructed event plane
+ if(gReactionPlane==999.) return;
+ }
+
+
+ TExMap *trackMap = new TExMap();
+
+// --- track loop for mapping matrix
+ if(fFilterBit==128)
+ {
+ 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;
+ Int_t tracktype=ClassifyTrack(track,trkVtx,bSign1,kFALSE);//don't fill the histos here
+ if(tracktype!=1) continue;
+
+ if(!track) continue;//for safety
+
+ Int_t gid = track->GetID();
+ trackMap->Add(gid,itrk);
+ }//track looop ends
+ }
+
+
+
+ //TObjArray* tracksUNID=0;
+ TObjArray* tracksUNID = new TObjArray;
+ tracksUNID->SetOwner(kTRUE);
+
+ //TObjArray* tracksID=0;
+ TObjArray* tracksID = new TObjArray;
+ tracksID->SetOwner(kTRUE);
+
+
+//get the selected v0 particle TObjArray
+ TObjArray* tracksIDV0 = new TObjArray;
+ tracksIDV0->SetOwner(kTRUE);
+
+ 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,kTRUE);//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)
+
+ if(fV0TrigCorr){
+if(IsTrackFromV0(aod,track)) continue;// remove auto correlation
+ }
+
+ AliAODTrack *PIDtrack=track;//for PID purpose, mainly important for TPC only tracks
+
+ if(fFilterBit==128){
+Int_t gid1 = track->GetID();
+//if(gid1>=0) PIDtrack = track;
+ PIDtrack =(AliAODTrack*) aod->GetTrack(trackMap->GetValue(-1-gid1));
+if(!PIDtrack) continue;//for safety; so that each of the TPC only tracks have corresponding global track along with it
+ }
+
+ 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_2_76" || fCentralityMethod.EndsWith("_MANUAL") || (fSampleType=="pp_7" && fPPVsMultUtils==kFALSE)) effcent=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
+ //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()<=fmaxPtAsso) || (track->Pt()>=fminPtTrig && track->Pt()<=fmaxPtTrig))//to reduce memory consumption in pool
+ {
+ Short_t chargeval=0;
+ if(track->Charge()>0) chargeval=1;
+ if(track->Charge()<0) chargeval=-1;
+ if(chargeval==0) continue;
+ if (fapplyTrigefficiency || fapplyAssoefficiency)//get the trackingefficiency x contamination factor for unidentified particles
+ effmatrix=GetTrackbyTrackeffvalue(track,effcent,zvtx,particletypeMC);
+ LRCParticlePID* copy = new LRCParticlePID(particletypeMC,Inv_mass,chargeval,track->Pt(),track->Eta(), track->Phi(),effmatrix,track->GetTPCClusterMapPtr(),track->GetTPCSharedMapPtr());
+ 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)
+ }
+
+//get the pdg code of the corresponding truth particle
+ Int_t pdgCode = ((AliAODMCParticle*)recomatched)->GetPdgCode();
+
+ Double_t allrecomatchedpid[4] = {effcent, zVtxmc,recomatched->Pt(), recomatched->Eta()};
+ if(ffillefficiency) {
+fTrackHistEfficiency[5]->Fill(allrecomatchedpid,2);//for allreco matched
+ if(TMath::Abs(pdgCode)==211 || TMath::Abs(pdgCode)==321) fTrackHistEfficiency[3]->Fill(allrecomatchedpid,2);//for mesons
+ if(TMath::Abs(pdgCode)==321 || TMath::Abs(pdgCode)==2212) fTrackHistEfficiency[4]->Fill(allrecomatchedpid,2);//for kaons+protons
+ if(TMath::Abs(pdgCode)==211) fTrackHistEfficiency[0]->Fill(allrecomatchedpid,2);//for pions
+ if(TMath::Abs(pdgCode)==321) fTrackHistEfficiency[1]->Fill(allrecomatchedpid,2);//for kaons
+ if(TMath::Abs(pdgCode)==2212) fTrackHistEfficiency[2]->Fill(allrecomatchedpid,2);//for protons
+
+ if (((AliAODMCParticle*)recomatched)->IsPhysicalPrimary()) {
+ fTrackHistEfficiency[5]->Fill(allrecomatchedpid,1);//for primreco matched
+ if(TMath::Abs(pdgCode)==211 || TMath::Abs(pdgCode)==321) fTrackHistEfficiency[3]->Fill(allrecomatchedpid,1);//for mesons
+ if(TMath::Abs(pdgCode)==321 || TMath::Abs(pdgCode)==2212) fTrackHistEfficiency[4]->Fill(allrecomatchedpid,1);//for kaons+protons
+ if( TMath::Abs(pdgCode)==211) fTrackHistEfficiency[0]->Fill(allrecomatchedpid,1);//for pions
+ if( TMath::Abs(pdgCode)==321) fTrackHistEfficiency[1]->Fill(allrecomatchedpid,1);//for kaons
+ if( TMath::Abs(pdgCode)==2212) fTrackHistEfficiency[2]->Fill(allrecomatchedpid,1);//for protons
+ }
+ }
+
+switch(TMath::Abs(pdgCode)){
+ case 2212:
+ if(fFIllPIDQAHistos){
+ for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
+ if((ipid!=NSigmaTPC) && (!HasTOFPID(PIDtrack)))continue;//not filling TOF and combined if no TOF PID
+ TH2F *h=GetHistogram2D(Form("NSigmaMC_%d_%d",SpProton,ipid));
+ h->Fill(track->Pt(),fnsigmas[SpProton][ipid]);
}
}
break;
case 321:
if(fFIllPIDQAHistos){
for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
- if((ipid!=NSigmaTPC) && (!HasTOFPID(track)))continue;//not filling TOF and combined if no TOF PID
+ if((ipid!=NSigmaTPC) && (!HasTOFPID(PIDtrack)))continue;//not filling TOF and combined if no TOF PID
TH2F *h=GetHistogram2D(Form("NSigmaMC_%d_%d",SpKaon,ipid));
h->Fill(track->Pt(),fnsigmas[SpKaon][ipid]);
}
case 211:
if(fFIllPIDQAHistos){
for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
- if((ipid!=NSigmaTPC) && (!HasTOFPID(track)))continue;//not filling TOF and combined if no TOF PID
+ if((ipid!=NSigmaTPC) && (!HasTOFPID(PIDtrack)))continue;//not filling TOF and combined if no TOF PID
TH2F *h=GetHistogram2D(Form("NSigmaMC_%d_%d",SpPion,ipid));
h->Fill(track->Pt(),fnsigmas[SpPion][ipid]);
}
}
+ //now start the particle identification process:)
+
+Float_t dEdx = PIDtrack->GetTPCsignal();
+ fHistoTPCdEdx->Fill(track->Pt(), dEdx);
+
+ if(HasTOFPID(PIDtrack))
+{
+Double_t beta = GetBeta(PIDtrack);
+fHistoTOFbeta->Fill(track->Pt(), beta);
+ }
+
+ //remove the tracks which don't have proper TOF response-otherwise the misIDentification rate values will be wrong
+if(fRequestTOFPID && track->Pt()>fPtTOFPIDmin && (!HasTOFPID(PIDtrack)) ) continue;
+
+
+//do track identification(nsigma method)
+ particletypeMC=GetParticle(PIDtrack,fFIllPIDQAHistos);//******************************problem is here
+
//2-d TPCTOF map(for each Pt interval)
- if(HasTOFPID(track)){
+ if(HasTOFPID(PIDtrack)){
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]);
}
- //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());
- }
- }
+ //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());
+ }
+ }
+
+ //fill tracking efficiency
+ if(ffillefficiency)
+ {
+ if(particletypeMC==SpPion || particletypeMC==SpKaon)
+ {
+ if(TMath::Abs(pdgCode)==211 || TMath::Abs(pdgCode)==321) {
+ fTrackHistEfficiency[3]->Fill(allrecomatchedpid,4);//for mesons
+ if (((AliAODMCParticle*)recomatched)->IsPhysicalPrimary()) fTrackHistEfficiency[3]->Fill(allrecomatchedpid,3);//for mesons
+ }
+ }
+ if(particletypeMC==SpKaon || particletypeMC==SpProton)
+ {
+ if(TMath::Abs(pdgCode)==321 || TMath::Abs(pdgCode)==2212) {
+ fTrackHistEfficiency[4]->Fill(allrecomatchedpid,4);//for kaons+protons
+ if (((AliAODMCParticle*)recomatched)->IsPhysicalPrimary()) fTrackHistEfficiency[4]->Fill(allrecomatchedpid,3);
+ }
+ }
+ if(particletypeMC==SpPion && TMath::Abs(pdgCode)==211) {
+ fTrackHistEfficiency[0]->Fill(allrecomatchedpid,4);//for pions
+ if (((AliAODMCParticle*)recomatched)->IsPhysicalPrimary()) fTrackHistEfficiency[0]->Fill(allrecomatchedpid,3);
+ }
+ if(particletypeMC==SpKaon && TMath::Abs(pdgCode)==321) {
+ fTrackHistEfficiency[1]->Fill(allrecomatchedpid,4);//for kaons
+if (((AliAODMCParticle*)recomatched)->IsPhysicalPrimary()) fTrackHistEfficiency[1]->Fill(allrecomatchedpid,3);
+ }
+ if(particletypeMC==SpProton && TMath::Abs(pdgCode)==2212){
+ fTrackHistEfficiency[2]->Fill(allrecomatchedpid,4);//for protons
+if (((AliAODMCParticle*)recomatched)->IsPhysicalPrimary()) fTrackHistEfficiency[2]->Fill(allrecomatchedpid,3);
+ }
+ }
+
//for misidentification fraction calculation(do it with tuneonPID)
if(particletypeMC==SpPion )
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());
}
- if(particletypeMC==SpUndefined )
+ if(particletypeMC==SpUndefined )//these undefined are not due to absence of proper TOF response, rather due to the PID method only
{
if(TMath::Abs(pdgCode)==211) fUNIDcont->Fill(1.,track->Pt());
if(TMath::Abs(pdgCode)==321) fUNIDcont->Fill(3.,track->Pt());
if(particletypeMC==SpUndefined) continue;
-
- //fill tracking efficiency
- if(ffillefficiency)
- {
- if(particletypeMC==SpPion || particletypeMC==SpKaon)
- {
- if(TMath::Abs(pdgCode)==211 || TMath::Abs(pdgCode)==321) {
- fTrackHistEfficiency[3]->Fill(allrecomatchedpid,4);//for mesons
- if (((AliAODMCParticle*)recomatched)->IsPhysicalPrimary()) fTrackHistEfficiency[3]->Fill(allrecomatchedpid,3);//for mesons
- }
- }
- if(particletypeMC==SpKaon || particletypeMC==SpProton)
- {
- if(TMath::Abs(pdgCode)==321 || TMath::Abs(pdgCode)==2212) {
- fTrackHistEfficiency[4]->Fill(allrecomatchedpid,4);//for kaons+protons
- if (((AliAODMCParticle*)recomatched)->IsPhysicalPrimary()) fTrackHistEfficiency[4]->Fill(allrecomatchedpid,3);
- }
- }
- if(particletypeMC==SpPion && TMath::Abs(pdgCode)==211) {
- fTrackHistEfficiency[0]->Fill(allrecomatchedpid,4);//for pions
- if (((AliAODMCParticle*)recomatched)->IsPhysicalPrimary()) fTrackHistEfficiency[0]->Fill(allrecomatchedpid,3);
- }
- if(particletypeMC==SpKaon && TMath::Abs(pdgCode)==321) {
- fTrackHistEfficiency[1]->Fill(allrecomatchedpid,4);//for kaons
-if (((AliAODMCParticle*)recomatched)->IsPhysicalPrimary()) fTrackHistEfficiency[1]->Fill(allrecomatchedpid,3);
- }
- if(particletypeMC==SpProton && TMath::Abs(pdgCode)==2212){
- fTrackHistEfficiency[2]->Fill(allrecomatchedpid,4);//for protons
-if (((AliAODMCParticle*)recomatched)->IsPhysicalPrimary()) fTrackHistEfficiency[2]->Fill(allrecomatchedpid,3);
- }
- }
+ if(fRequestEventPlane){
+ FillPIDEventPlane(cent_v0,particletypeMC,track->Phi(),gReactionPlane);
+ }
if((track->Pt()>=fminPtAsso && track->Pt()<=fmaxPtAsso) || (track->Pt()>=fminPtTrig && track->Pt()<=fmaxPtTrig))//to reduce memory consumption in pool
{
if(chargeval==0) continue;
if (fapplyTrigefficiency || fapplyAssoefficiency)
effmatrix=GetTrackbyTrackeffvalue(track,effcent,zvtx,particletypeMC);//get the tracking eff x TOF matching eff x PID eff x contamination factor for identified particles
- LRCParticlePID* copy1 = new LRCParticlePID(particletypeMC,chargeval,track->Pt(),track->Eta(), track->Phi(),effmatrix);
+ LRCParticlePID* copy1 = new LRCParticlePID(particletypeMC,Inv_mass,chargeval,track->Pt(),track->Eta(), track->Phi(),effmatrix,track->GetTPCClusterMapPtr(),track->GetTPCSharedMapPtr());
copy1->SetUniqueID(eventno * 100000 + (Int_t)trackscount);
tracksID->Add(copy1);
}
//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??????
+ if (fSampleType=="pPb" || fSampleType=="PbPb" || fPPVsMultUtils==kTRUE || fCentralityMethod == "MC_b") 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(11);
+ fEventCounter->Fill(13);
//***************************************event no. counting
Bool_t isbaryontrig=kFALSE;
if (ismesontrig) fEventnomeson->Fill(cent_v0,zvtx);
}
+
+ if(fV0TrigCorr){
+ tracksIDV0=GetV0Particles((AliVEvent*) aod, cent_v0);
+ if(tracksIDV0->GetEntriesFast()<=0) return;
+ }
//same event delte-eta, delta-phi plot
if(tracksUNID && tracksUNID->GetEntriesFast()>0)//hadron triggered correlation
{//same event calculation starts
if(tracksID && tracksID->GetEntriesFast()>0)//ID triggered correlation
{//same event calculation starts
- if(tracksUNID && tracksUNID->GetEntriesFast()>0 && ffilltrigIDassoUNID) Fillcorrelation(gReactionPlane,tracksID,tracksUNID,cent_v0,zvtx,weghtval,kFALSE,bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kFALSE,"trigIDassoUNID");//mixcase=kFALSE (ID-hadron correlation)
+ if(tracksUNID && tracksUNID->GetEntriesFast()>0 && ffilltrigIDassoUNID) {
+ if(fV0TrigCorr) Fillcorrelation(gReactionPlane,tracksIDV0,tracksUNID,cent_v0,zvtx,weghtval,kFALSE,bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kFALSE,"trigIDassoUNID");//mixcase=kFALSE (ID-hadron correlation)
+
+ else Fillcorrelation(gReactionPlane,tracksID,tracksUNID,cent_v0,zvtx,weghtval,kFALSE,bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kFALSE,"trigIDassoUNID");//mixcase=kFALSE (ID-hadron correlation)
+ }
if(ffilltrigIDassoID) Fillcorrelation(gReactionPlane,tracksID,0,cent_v0,zvtx,weghtval,kFALSE,bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,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)
+ AliEventPool* pool = fPoolMgr->GetEventPool(cent_v0, zvtx,gReactionPlane);//In the pool there is tracksUNID(i.e associateds are unidentified)
if (pool && pool->IsReady())
{//start mixing only when pool->IsReady
Float_t nmix1=(Float_t)pool->GetCurrentNEvents();
if(ffilltrigassoUNID && tracksUNID && tracksUNID->GetEntriesFast()>0)//*******************************hadron trggered mixing
Fillcorrelation(gReactionPlane,tracksUNID,bgTracks,cent_v0,zvtx,nmix1,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigassoUNID");//mixcase=kTRUE
if(ffilltrigIDassoUNID && tracksID && tracksID->GetEntriesFast()>0)//***********************************ID trggered mixing
- Fillcorrelation(gReactionPlane,tracksID,bgTracks,cent_v0,zvtx,nmix1,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigIDassoUNID");//mixcase=kTRUE
+ {
+ if(fV0TrigCorr) Fillcorrelation(gReactionPlane,tracksIDV0,bgTracks,cent_v0,zvtx,nmix1,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigIDassoUNID");//mixcase=kTRUE
+
+ else Fillcorrelation(gReactionPlane,tracksID,bgTracks,cent_v0,zvtx,nmix1,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigIDassoUNID");//mixcase=kTRUE
+ }
}// pool event loop ends mixing case
} //if pool->IsReady() condition ends mixing case
}//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)
+ AliEventPool* pool1 = fPoolMgr->GetEventPool(cent_v0, zvtx+100,gReactionPlane);//In the pool1 there is tracksID(i.e associateds are identified)
if (pool1 && pool1->IsReady())
{//start mixing only when pool->IsReady
Float_t nmix2=(Float_t)pool1->GetCurrentNEvents();
}//mixing with identified particles
//no. of events analyzed
-fEventCounter->Fill(13);
+fEventCounter->Fill(15);
}
if(tracksUNID) delete tracksUNID;
if(tracksID) delete tracksID;
+if(fV0TrigCorr) {
+ if(tracksIDV0) delete tracksIDV0;
+ }
-PostData(1, fOutput);
+}//AOD || MCAOD condition
+
+//still in the main event loop
}
//________________________________________________________________________
AliError("Cannot get the AOD event");
return;
}
+ Double_t Inv_mass=0.0;
// count all events
fEventCounter->Fill(1);
-if (!fPID) return;//this should be available with each event even if we don't do PID selection
-
+ fPID = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetPIDResponse();
+ if (!fPID) AliFatal("This Task needs the PID response attached to the inputHandler");
+ //if (!fPID) return;//this should be available with each event even if we don't do PID selection
+
+ fgPsi2v0a=999.;
+ fgPsi2v0c=999.;
+ fgPsi2tpc=999.;
+ fgPsi3v0a=999.;
+ fgPsi3v0c=999.;
+ fgPsi3tpc=999.;
+ gReactionPlane = 999.;
+
Double_t cent_v0= -999.;
Double_t effcent=1.0;
- Double_t gReactionPlane = -1.;
Float_t bSign = 0.;
Double_t trackscount=0;//counts particles passed filterbit cuts and kinematic cuts used in this analysis
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
+ Float_t bSign1=aod->GetMagneticField() ;//for dca cut in ClassifyTrack(), i.e in track loop
// check event cuts and fill event histograms and return the centrality or reference multiplicity value
- if((cent_v0 = GetAcceptedEventMultiplicity(aod,kFALSE)) < 0){
+ if((cent_v0 = GetAcceptedEventMultiplicity((AliVEvent*)aod,kFALSE)) < 0){
return;
}
-
+ effcent=cent_v0;//required for efficiency correction case********Extremely Important
//get the event plane in case of PbPb
- if(fSampleType=="PbPb"){
if(fRequestEventPlane){
- gReactionPlane = GetEventPlane(aod,kFALSE);
- fHistEventPlaneReco->Fill(gReactionPlane,cent_v0);
- if(gReactionPlane < 0) return;
+ gReactionPlane = GetEventPlane((AliVEvent*)aod,kFALSE,cent_v0);
+ if(gReactionPlane==999.) return;
}
- }
+
+
+TExMap *trackMap = new TExMap();
+// --- track loop for mapping matrix
+ if(fFilterBit==128)
+ {
+ 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;
+ Int_t tracktype=ClassifyTrack(track,trkVtx,bSign1,kFALSE);//don't fill the histos here
+ if(tracktype!=1) continue;
+
+ if(!track) continue;//for safety
+
+ Int_t gid = track->GetID();
+ trackMap->Add(gid,itrk);
+ }//track looop ends
+ }
TObjArray* tracksUNID= new TObjArray;//track info before doing PID
tracksUNID->SetOwner(kTRUE); // IMPORTANT!
TObjArray* tracksID= new TObjArray;//only pions, kaons,protons i.e. after doing the PID selection
tracksID->SetOwner(kTRUE); // IMPORTANT!
+ //get the selected v0 particle TObjArray
+ TObjArray* tracksIDV0= new TObjArray;
+ tracksIDV0->SetOwner(kTRUE); // IMPORTANT!
eventno++;
if(!track) continue;//for safety
+ if(fV0TrigCorr){
+if(IsTrackFromV0(aod,track)) continue;// remove auto correlation
+ }
+AliAODTrack *PIDtrack=track;//for PID purpose, mainly important for TPC only tracks
+
+ if(fFilterBit==128){
+Int_t gid1 = track->GetID();
+//if(gid1>=0) PIDtrack = track;
+ PIDtrack =(AliAODTrack*) aod->GetTrack(trackMap->GetValue(-1-gid1));
+if(!PIDtrack) continue;//for safety; so that each of the TPC only tracks have corresponding global track along with it
+ }
+
//check for eta , phi holes
fEtaSpectrasso->Fill(track->Eta(),track->Pt());
fphiSpectraasso->Fill(track->Phi(),track->Pt());
if(track->Pt()>=fmaxPtTrig) fTrigPtJet=kTRUE;
- if (fSampleType=="pp") effcent=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
+if (fSampleType=="pp_2_76" || fCentralityMethod.EndsWith("_MANUAL") || (fSampleType=="pp_7" && fPPVsMultUtils==kFALSE)) effcent=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(chargeval==0) continue;
if (fapplyTrigefficiency || fapplyAssoefficiency)//get the trackingefficiency x contamination factor for unidentified particles
effmatrix=GetTrackbyTrackeffvalue(track,effcent,zvtx,particletype);
- LRCParticlePID* copy = new LRCParticlePID(particletype,chargeval,track->Pt(),track->Eta(), track->Phi(),effmatrix);
+ LRCParticlePID* copy = new LRCParticlePID(particletype,Inv_mass,chargeval,track->Pt(),track->Eta(), track->Phi(),effmatrix,track->GetTPCClusterMapPtr(),track->GetTPCSharedMapPtr());
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)
}
//track passing filterbit 768 have proper TPC response,or need to be checked explicitly before doing PID????
- Float_t dEdx = track->GetTPCsignal();
+ Float_t dEdx = PIDtrack->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))
+ if(HasTOFPID(PIDtrack))
{
- Double_t beta = GetBeta(track);
+ Double_t beta = GetBeta(PIDtrack);
fHistoTOFbeta->Fill(track->Pt(), beta);
}
-
+ //remove the tracks which don't have proper TOF response-otherwise the misIDentification rate values will be wrong(in MC)
+if(fRequestTOFPID && track->Pt()>fPtTOFPIDmin && (!HasTOFPID(PIDtrack)) ) continue;
+
//track identification(using nsigma method)
- particletype=GetParticle(track,fFIllPIDQAHistos);//*******************************change may be required(It should return only pion,kaon, proton and Spundefined; NOT unidentifed***************be careful)
-
+ particletype=GetParticle(PIDtrack,fFIllPIDQAHistos);//*******************************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)){
+ if(HasTOFPID(PIDtrack)){
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)!!!!!!!!!!!
+ if(fRequestEventPlane){
+ FillPIDEventPlane(cent_v0,particletype,track->Phi(),gReactionPlane);
+ }
+
+
//Pt, Eta , Phi distribution of the reconstructed identified particles
if(ffillhistQAReco)
{
if(chargeval==0) continue;
if (fapplyTrigefficiency || fapplyAssoefficiency)
effmatrix=GetTrackbyTrackeffvalue(track,effcent,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,chargeval,track->Pt(),track->Eta(), track->Phi(),effmatrix);
+ LRCParticlePID* copy1 = new LRCParticlePID(particletype,Inv_mass,chargeval,track->Pt(),track->Eta(), track->Phi(),effmatrix,track->GetTPCClusterMapPtr(),track->GetTPCSharedMapPtr());
copy1->SetUniqueID(eventno * 100000 + (Int_t)trackscount);
tracksID->Add(copy1);
}
//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??????
+if(fSampleType=="pPb" || fSampleType=="PbPb" || fPPVsMultUtils==kTRUE) 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(11);
+ fEventCounter->Fill(13);
//***************************************event no. counting
Bool_t isbaryontrig=kFALSE;
if (ismesontrig) fEventnomeson->Fill(cent_v0,zvtx);
}
+//Get the TObjArray of V0 particles
+ if(fV0TrigCorr){
+ tracksIDV0=GetV0Particles((AliVEvent*) aod,cent_v0);
+ if(tracksIDV0->GetEntriesFast()<=0) return;
+ }
//same event delta-eta-deltaphi plot
-
if(tracksUNID && tracksUNID->GetEntriesFast()>0)//hadron triggered correlation
{//same event calculation starts
if(ffilltrigassoUNID) Fillcorrelation(gReactionPlane,tracksUNID,0,cent_v0,zvtx,weightval,kFALSE,bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kFALSE,"trigassoUNID");//mixcase=kFALSE (hadron-hadron correlation)
- if(tracksID && tracksID->GetEntriesFast()>0 && ffilltrigUNIDassoID) Fillcorrelation(gReactionPlane,tracksUNID,tracksID,cent_v0,zvtx,weightval,kFALSE,bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kFALSE,"trigUNIDassoID");//mixcase=kFALSE (hadron-ID correlation)
+ if(tracksID && tracksID->GetEntriesFast()>0 && ffilltrigUNIDassoID) Fillcorrelation(gReactionPlane,tracksUNID,tracksID,cent_v0,zvtx,weightval,kFALSE,bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,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(gReactionPlane,tracksID,tracksUNID,cent_v0,zvtx,weightval,kFALSE,bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kFALSE,"trigIDassoUNID");//mixcase=kFALSE (ID-hadron correlation)
+ if(tracksUNID && tracksUNID->GetEntriesFast()>0 && ffilltrigIDassoUNID){
+if(fV0TrigCorr) Fillcorrelation(gReactionPlane,tracksIDV0,tracksUNID,cent_v0,zvtx,weightval,kFALSE,bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kFALSE,"trigIDassoUNID");//mixcase=kFALSE (ID-hadron correlation)
+else Fillcorrelation(gReactionPlane,tracksID,tracksUNID,cent_v0,zvtx,weightval,kFALSE,bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kFALSE,"trigIDassoUNID");//mixcase=kFALSE (ID-hadron correlation)
+ }
if(ffilltrigIDassoID) Fillcorrelation(gReactionPlane,tracksID,0,cent_v0,zvtx,weightval,kFALSE,bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kFALSE,"trigIDassoID");//mixcase=kFALSE (ID-ID correlation)
}
//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)
+AliEventPool* pool = fPoolMgr->GetEventPool(cent_v0, zvtx,gReactionPlane);//In the pool there is tracksUNID(i.e associateds are unidentified)
if (pool && pool->IsReady())
{//start mixing only when pool->IsReady
Float_t nmix1=(Float_t)pool->GetCurrentNEvents();
if(ffilltrigassoUNID && tracksUNID && tracksUNID->GetEntriesFast()>0)//*******************************hadron trggered mixing
Fillcorrelation(gReactionPlane,tracksUNID,bgTracks,cent_v0,zvtx,nmix1,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigassoUNID");//mixcase=kTRUE
if(ffilltrigIDassoUNID && tracksID && tracksID->GetEntriesFast()>0)//***********************************ID trggered mixing
- Fillcorrelation(gReactionPlane,tracksID,bgTracks,cent_v0,zvtx,nmix1,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigIDassoUNID");//mixcase=kTRUE
+ {
+if(fV0TrigCorr) Fillcorrelation(gReactionPlane,tracksIDV0,bgTracks,cent_v0,zvtx,nmix1,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigIDassoUNID");//mixcase=kTRUE
+else Fillcorrelation(gReactionPlane,tracksID,bgTracks,cent_v0,zvtx,nmix1,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigIDassoUNID");//mixcase=kTRUE
+ }
}// pool event loop ends mixing case
} //if pool->IsReady() condition ends mixing case
if(tracksUNID) {
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)
+ AliEventPool* pool1 = fPoolMgr->GetEventPool(cent_v0, zvtx+100,gReactionPlane);//In the pool1 there is tracksID(i.e associateds are identified)
if (pool1 && pool1->IsReady())
{//start mixing only when pool->IsReady
Float_t nmix2=(Float_t)pool1->GetCurrentNEvents();
//no. of events analyzed
-fEventCounter->Fill(13);
+fEventCounter->Fill(15);
//still in main event loop
if(tracksID) delete tracksID;
+if(fV0TrigCorr) {
+ if(tracksIDV0) delete tracksIDV0;
+ }
-PostData(1, fOutput);
-
-} // *************************event loop ends******************************************//_______________________________________________________________________
+} // *************************event loop ends******************************************
+//_______________________________________________________________________
TObjArray* AliTwoParticlePIDCorr::CloneAndReduceTrackList(TObjArray* tracks)
{
// clones a track list by using AliDPhiBasicParticle which uses much less memory (used for event mixing)
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());
+ LRCParticlePID* copy100 = new LRCParticlePID(particle->getparticle(),particle->GetInvMass(),particle->Charge(), particle->Pt(),particle->Eta(), particle->Phi(), particle->geteffcorrectionval(),particle->GetTPCPadMap(),particle->GetTPCSharedMap());
copy100->SetUniqueID(particle->GetUniqueID());
tracksClone->Add(copy100);
}
}
//--------------------------------------------------------------------------------
-void AliTwoParticlePIDCorr::Fillcorrelation(Double_t gReactionPlane,TObjArray *trackstrig,TObjArray *tracksasso,Double_t cent,Float_t vtx,Float_t weight,Bool_t firstTime,Float_t bSign,Bool_t fPtOrder,Bool_t twoTrackEfficiencyCut,Bool_t mixcase,TString fillup)
+void AliTwoParticlePIDCorr::Fillcorrelation(Float_t ReactionPlane,TObjArray *trackstrig,TObjArray *tracksasso,Double_t cent,Float_t vtx,Float_t weight,Bool_t firstTime,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
Float_t trigpt=trig->Pt();
//to avoid overflow qnd underflow
if(trigpt<fminPtTrig || trigpt>fmaxPtTrig) continue;
+ Double_t ParticlePID_InvMass=0.0;
+ if(fV0TrigCorr) ParticlePID_InvMass=trig->GetInvMass();
+ else{
Int_t particlepidtrig=trig->getparticle();
- if(fTriggerSpeciesSelection){ if (particlepidtrig!=fTriggerSpecies) continue;}
+ ParticlePID_InvMass=(Double_t) particlepidtrig;
+ if(fTriggerSpeciesSelection){ if (particlepidtrig!=fTriggerSpecies) continue;}//***********************************forks,lam.Alam their PID numbers have no meaning, only their Inv_mass will be stored
+ }
Float_t trigeta=trig->Eta();
Double_t gPsiMinusPhi = 0.;
Double_t gPsiMinusPhiBin = -10.;
if(fRequestEventPlane){
- gPsiMinusPhi = TMath::Abs(trigphi - gReactionPlane);
- //in-plane
+ gPsiMinusPhi = TMath::Abs(trigphi - ReactionPlane);
+ //in-plane(Note thet event plane angle has to be defined within 0 to 180 degree(do not use this if event ), otherwise the definition of in plane and out plane particles is wrong)
if((gPsiMinusPhi <= 7.5*TMath::DegToRad())||
+ (gPsiMinusPhi >= 352.5*TMath::DegToRad())||
+ ((172.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 187.5*TMath::DegToRad())))
+ gPsiMinusPhiBin = 0.0;
+ /*
+ if((gPsiMinusPhi <= 7.5*TMath::DegToRad())||
((172.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 187.5*TMath::DegToRad())))
gPsiMinusPhiBin = 0.0;
+ */
//intermediate
else if(((37.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 52.5*TMath::DegToRad()))||
((127.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 142.5*TMath::DegToRad()))||
if(fRequestEventPlane){
trigval[3] = gPsiMinusPhiBin;
- if(fcontainPIDtrig && SetChargeAxis==0) trigval[4] = particlepidtrig;
+ if(fcontainPIDtrig && SetChargeAxis==0) trigval[4] = ParticlePID_InvMass;
if(!fcontainPIDtrig && SetChargeAxis==2) trigval[4] = trig->Charge();
if(fcontainPIDtrig && SetChargeAxis==2) {
- trigval[4] = particlepidtrig;
+ trigval[4] = ParticlePID_InvMass;
trigval[5] = trig->Charge();
}
}
if(!fRequestEventPlane){
- if(fcontainPIDtrig && SetChargeAxis==0) trigval[3] = particlepidtrig;
+ if(fcontainPIDtrig && SetChargeAxis==0) trigval[3] = ParticlePID_InvMass;
if(!fcontainPIDtrig && SetChargeAxis==2) trigval[3] = trig->Charge();
if(fcontainPIDtrig && SetChargeAxis==2) {
- trigval[3] = particlepidtrig;
+ trigval[3] = ParticlePID_InvMass;
trigval[4] = trig->Charge();
}
}
if(!tracksasso && j==i) 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,i.e. both Trig and asso TObjArray belongs to the same Pi range but say Trig is Unidentified but asso is identified then the serial no. wise particles are not same and and j==i doesn't aplly)
+ // check if both particles point to the same element (does not occur for mixed events, but if subsets are mixed within the same event,i.e. both Trig and asso TObjArray belongs to the same Pt range but say Trig is Unidentified but asso is identified then the serial no. wise particles are not same and and j==i doesn't aplly)
// if (tracksasso && trig->IsEqual(asso)) continue;
if (tracksasso && (trig->GetUniqueID()==asso->GetUniqueID())) continue;
{
// 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);
+ Float_t dphistar1 = GetDPhiStar(phi1, pt1, charge1, phi2, pt2, charge2, fTwoTrackCutMinRadius, bSign);
+ Float_t dphistar2 = GetDPhiStar(phi1, pt1, charge1, phi2, pt2, charge2, fTwoTrackCutMaxRadius, bSign);
const Float_t kLimit = twoTrackEfficiencyCutValue * 3;
if (TMath::Abs(dphistar1) < kLimit || TMath::Abs(dphistar2) < kLimit || dphistar1 * dphistar2 < 0)
{
- for (Double_t rad=0.8; rad<2.51; rad+=0.01)
+ for (Double_t rad=fTwoTrackCutMinRadius; rad<=fTwoTrackCutMaxRadius; rad+=0.01)
{
Float_t dphistar = GetDPhiStar(phi1, pt1, charge1, phi2, pt2, charge2, rad, bSign);
dphistarminabs = dphistarabs;
}
}
+ if(mixcase==kFALSE) fTwoTrackDistancePt[0]->Fill(deta, dphistarmin, TMath::Abs(pt1 - pt2));//for same event
+ if(mixcase==kTRUE) fTwoTrackDistancePtmix[0]->Fill(deta, dphistarmin, TMath::Abs(pt1 - pt2));//for mixed event
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));
+ if(mixcase==kFALSE) fTwoTrackDistancePt[1]->Fill(deta, dphistarmin, TMath::Abs(pt1 - pt2));//for same event
+ if(mixcase==kTRUE) fTwoTrackDistancePtmix[1]->Fill(deta, dphistarmin, TMath::Abs(pt1 - pt2));//for mixed event
}
}
}
+ //pair sharedfraction cut(only between the trig and asso track)
+ if(fillup!="trigIDassoIDMCTRUTH")//******************************************NOT for TRUTH MC particles
+ {
+ if(fSharedfraction_Pair_cut>=0){
+ Bool_t passsharedfractionpaircut=CalculateSharedFraction(trig->GetTPCPadMap(),asso->GetTPCPadMap(),trig->GetTPCSharedMap(),asso->GetTPCSharedMap());
+ if(!passsharedfractionpaircut) continue;
+ }
+ }
Float_t weightperevent=weight;
Float_t trackeffasso=1.0;
if(fapplyAssoefficiency) trackeffasso=asso->geteffcorrectionval();
Float_t deleta=trigeta-eta[j];
Float_t delphi=PhiRange(trigphi-asso->Phi());
+ Float_t delpt=trigpt-asso->Pt();
+ //fill it with/without two track efficiency cut
+ if(mixcase==kFALSE) fTwoTrackDistancePtdip->Fill(deleta, delphi, TMath::Abs(delpt));//for same event
+ if(mixcase==kTRUE) fTwoTrackDistancePtdipmix->Fill(deleta, delphi, TMath::Abs(delpt));//for mixed event
+
//here get the two particle efficiency correction factor
Float_t effweight=trackefftrig*trackeffasso*weightperevent;
// if(mixcase==kFALSE) cout<<"*******************effweight="<<effweight<<endl;
vars[dimension+1]=asso->Charge();
}
if(fcontainPIDtrig && !fcontainPIDasso){
- vars[dimension]=particlepidtrig;
+ vars[dimension]=ParticlePID_InvMass;
if(SetChargeAxis==2){
vars[dimension+1]=trig->Charge();
vars[dimension+2]=asso->Charge();
}
}
if(fcontainPIDtrig && fcontainPIDasso){
- vars[dimension]=particlepidtrig;
+ vars[dimension]=ParticlePID_InvMass;
vars[dimension+1]=particlepidasso;
if(SetChargeAxis==2){
vars[dimension+2]=trig->Charge();
}
}
-//--------------------------------------------------------------------------------
+//------------------------------------------------------------------------------------------------
+Bool_t AliTwoParticlePIDCorr:: CalculateSharedFraction(const TBits *triggerPadMap,const TBits *assocPadMap,const TBits *triggerShareMap,const TBits *assocShareMap)
+{//source code-AliFemtoShareQualityPairCut.cxx
+Double_t nofhits=0;
+Double_t nofsharedhits=0;
+
+for(UInt_t imap=0;imap< (triggerPadMap->GetNbits() );imap++)
+{
+//if they are in same pad
+//cout<<triggerPadMap->TestBitNumber(imap)<<" "<< assocPadMap->TestBitNumber(imap)<<endl;
+if (triggerPadMap->TestBitNumber(imap) &&
+ assocPadMap->TestBitNumber(imap))
+{
+//if they share
+//cout<<triggerShareMap->TestBitNumber(imap)<<" "<<assocShareMap->TestBitNumber(imap)<<endl;
+if (triggerShareMap->TestBitNumber(imap) &&
+ assocShareMap->TestBitNumber(imap))
+{
+ //cout<<triggerShareMap->TestBitNumber(imap)<<" "<<assocShareMap->TestBitNumber(imap)<<endl;
+nofhits+=2;
+nofsharedhits+=2;
+}
+
+
+
+//not shared
+ else {
+
+ nofhits+=2;
+ }
+
+
+}
+//different pad
+
+//cout<< (triggerPadMap->TestBitNumber(imap) || assocPadMap->TestBitNumber(imap))<<endl;
+else if (triggerPadMap->TestBitNumber(imap) ||
+ assocPadMap->TestBitNumber(imap)) {
+ // One track has a hit, the other does not
+
+ nofhits++;
+ //cout<<"No hits :"<<nofhits<<endl;
+
+ }
+
+
+
+}
+
+Double_t SharedFraction=0.0;
+if(nofhits>0) SharedFraction=(nofsharedhits/nofhits);
+
+//cout<<"Fraction shared hits :"<<SharedFraction<<endl;
+
+if(SharedFraction>fSharedfraction_Pair_cut) return kFALSE;
+
+return kTRUE;
+
+}
+
+//________________________________________________________________________________________________
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
if (frac > fSharedClusterCut)
return 0;
}
+
+ // Rejects tracks with shared clusters after filling a control histogram
+ // This overload is used for primaries
+
+ // Get the shared maps
+ const TBits sharedMap = track->GetTPCSharedMap();
+ // Fill a control histogram
+ fPriHistShare->Fill(sharedMap.CountBits());
+
+ // Reject shared clusters
+ if (fSharedTPCmapCut >= 0)
+ {
+ if((sharedMap.CountBits()) >= 1) return 0;// Bad track, has too many shared clusters!
+ }
+
if (fill) fHistQA[8]->Fill(pt);
if (fill) fHistQA[9]->Fill(track->Eta());
if (fill) fHistQA[10]->Fill(track->Phi());
return dphistar;
}
-//_________________________________________________________________________
-void AliTwoParticlePIDCorr ::DefineEventPool()
-{
-Int_t MaxNofEvents=1000;
-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=4;
- Double_t CentralityBins[NofCentBins+1]={0.,20., 40.,60.,200.1};//Is This binning is fine for pp, or we don't require them....
-fPoolMgr = new AliEventPoolManager(MaxNofEvents,fMaxNofMixingTracks,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,fMaxNofMixingTracks,NofCentBins,CentralityBins,NofVrtxBins,ZvrtxBins);
- }
-fPoolMgr->SetTargetValues(fMaxNofMixingTracks, 0.1, 5);
-
-//if(!fPoolMgr) return kFALSE;
-//return kTRUE;
-}
//------------------------------------------------------------------------
Double_t* AliTwoParticlePIDCorr::GetBinning(const char* configuration, const char* tag, Int_t& nBins)
{
return 1.0;
}
//________________________________________________________________________
-Double_t AliTwoParticlePIDCorr::GetReferenceMultiplicityVZEROFromAOD(AliAODEvent *event){
- //Function that returns the reference multiplicity from AODs (data or reco MC)
+Double_t AliTwoParticlePIDCorr::GetReferenceMultiplicityVZEROFromAOD(AliVEvent *mainevent){
+ //Function that returns the reference multiplicity from AODs (data or reco MC, Not for Truth)
//Different ref. mult. implemented: V0M, V0A, V0C, TPC
+ if(!mainevent) return -1;
+
+ AliAODEvent* event = dynamic_cast<AliAODEvent*>(mainevent);
+
Double_t gRefMultiplicity = 0., gRefMultiplicityTPC = 0.;
Double_t gRefMultiplicityVZERO = 0., gRefMultiplicityVZEROA = 0., gRefMultiplicityVZEROC = 0.;
}//track looop ends
+ if(fCentralityMethod == "V0A_MANUAL" || fCentralityMethod == "V0M_MANUAL" || fCentralityMethod == "V0C_MANUAL" ){
//VZERO segmentation in two detectors (0-31: VZERO-C, 32-63: VZERO-A)
for(Int_t iChannel = 0; iChannel < 64; iChannel++) {
fHistVZEROSignal->Fill(iChannel,event->GetVZEROEqMultiplicity(iChannel));
//EQVZEROC vs EQVZEROA multiplicity
fHistEQVZEROCvsEQVZEROAmultiplicity->Fill(gRefMultiplicityVZEROC,gRefMultiplicityVZEROA);
+ }
+ fHistRefmult->Fill(3.,gRefMultiplicityTPC);
+ fHistRefmult->Fill(2.,gRefMultiplicityVZERO);
+ fHistRefmult->Fill(0.,gRefMultiplicityVZEROA);
+ fHistRefmult->Fill(1.,gRefMultiplicityVZEROC);
- if(fCentralityMethod == "TRACKS_MANUAL")
- gRefMultiplicity = gRefMultiplicityTPC;
- else if(fCentralityMethod == "V0M_MANUAL")
- gRefMultiplicity = gRefMultiplicityVZERO;
- else if(fCentralityMethod == "V0A_MANUAL")
- gRefMultiplicity = gRefMultiplicityVZEROA;
- else if(fCentralityMethod == "V0C_MANUAL")
- gRefMultiplicity = gRefMultiplicityVZEROC;
-
- //ref mult QA
- fHistRefmult->Fill(0.,gRefMultiplicityVZEROA);
- fHistRefmult->Fill(1.,gRefMultiplicityVZEROC);
- fHistRefmult->Fill(2.,gRefMultiplicityVZERO);
- fHistRefmult->Fill(3.,gRefMultiplicityTPC);
+ if(fCentralityMethod == "TRACKS_MANUAL") gRefMultiplicity = gRefMultiplicityTPC;
+
+ else if(fCentralityMethod == "V0M_MANUAL") gRefMultiplicity = gRefMultiplicityVZERO;
-
+ else if(fCentralityMethod == "V0A_MANUAL") gRefMultiplicity = gRefMultiplicityVZEROA;
+
+ else if(fCentralityMethod == "V0C_MANUAL") gRefMultiplicity = gRefMultiplicityVZEROC;
+
+ else gRefMultiplicity = gRefMultiplicityTPC;
+
return gRefMultiplicity;
}
//-------------------------------------------------------------------------------------------------------
-Double_t AliTwoParticlePIDCorr::GetRefMultiOrCentrality(AliAODEvent *event, Bool_t truth){
+Double_t AliTwoParticlePIDCorr::GetRefMultiOrCentrality(AliVEvent *mainevent, Bool_t truth){
- if(!event) return -1;
+ if(!mainevent) return -1;
// get centrality object and check quality
Double_t cent_v0=-1;
- Double_t nooftrackstruth=0;
+ Bool_t shift_to_TRACKS_MANUAL=kFALSE;//in case of wrong setting automatic shift to Tracks_Manual method
+
+ Double_t gRefMultiplicityTPC_Truth = 0.;
+ Double_t gRefMultiplicityVZERO_Truth = 0., gRefMultiplicityVZEROA_Truth = 0., gRefMultiplicityVZEROC_Truth = 0.;
+
+ if(fAnalysisType == "AOD"|| fAnalysisType == "MCAOD") { //centrality in AOD header //++++++++++++++
+ AliAODEvent* event = dynamic_cast<AliAODEvent*>(mainevent);
if(fCentralityMethod=="V0M" || fCentralityMethod=="V0A" || fCentralityMethod=="V0C" || fCentralityMethod=="CL1" || fCentralityMethod=="ZNA" || fCentralityMethod=="V0AEq" || fCentralityMethod=="V0CEq" || fCentralityMethod=="V0MEq")//for PbPb, pPb, pp7TeV(still to be introduced)//data or RecoMC and also for TRUTH
{
+
+if(fSampleType=="pp_7" && fPPVsMultUtils)
+ {//for pp 7 TeV case only using Alianalysisutils class
+ if(fAnalysisUtils) cent_v0 = fAnalysisUtils->GetMultiplicityPercentile((AliVEvent*)event,fCentralityMethod);
+ else cent_v0 = -1;
+ fHistCentStats->Fill(0.,fAnalysisUtils->GetMultiplicityPercentile((AliVEvent*)event,"V0A"));
+ fHistCentStats->Fill(1.,fAnalysisUtils->GetMultiplicityPercentile((AliVEvent*)event,"V0C"));
+ fHistCentStats->Fill(2.,fAnalysisUtils->GetMultiplicityPercentile((AliVEvent*)event,"V0M"));
+ fHistCentStats->Fill(3.,fAnalysisUtils->GetMultiplicityPercentile((AliVEvent*)event,"V0AEq"));//only available for LHC10d at present (Quantile info)
+ fHistCentStats->Fill(4.,fAnalysisUtils->GetMultiplicityPercentile((AliVEvent*)event,"V0CEq"));//only available for LHC10d at present (Quantile info)
+ fHistCentStats->Fill(5.,fAnalysisUtils->GetMultiplicityPercentile((AliVEvent*)event,"V0MEq"));//only available for LHC10d at present (Quantile info)
+ }
+
+else if(fSampleType=="pPb" || fSampleType=="PbPb")
+ {
AliCentrality *centralityObj=0;
AliAODHeader *header = (AliAODHeader*) event->GetHeader();
if(!header) return -1;
centralityObj = header->GetCentralityP();
// if (centrality->GetQuality() != 0) return ;
-
- if(centralityObj)
- {
+ if(centralityObj){
fHistCentStats->Fill(0.,centralityObj->GetCentralityPercentile("V0A"));
fHistCentStats->Fill(1.,centralityObj->GetCentralityPercentile("V0C"));
fHistCentStats->Fill(2.,centralityObj->GetCentralityPercentile("V0M"));
-if(fSampleType=="pp") fHistCentStats->Fill(3.,centralityObj->GetCentralityPercentile("V0AEq"));//only available for LHC10d at present (Quantile info)
-if(fSampleType=="pp") fHistCentStats->Fill(4.,centralityObj->GetCentralityPercentile("V0CEq"));//only available for LHC10d at present (Quantile info)
-if(fSampleType=="pp") fHistCentStats->Fill(5.,centralityObj->GetCentralityPercentile("V0MEq"));//only available for LHC10d at present (Quantile info)
+ fHistCentStats->Fill(3.,centralityObj->GetCentralityPercentile("V0AEq"));//only available for LHC10d at present (Quantile info)
+ fHistCentStats->Fill(4.,centralityObj->GetCentralityPercentile("V0CEq"));//only available for LHC10d at present (Quantile info)
+ fHistCentStats->Fill(5.,centralityObj->GetCentralityPercentile("V0MEq"));//only available for LHC10d at present (Quantile info)
-if(fSampleType=="pPb" || fSampleType=="PbPb") fHistCentStats->Fill(6.,centralityObj->GetCentralityPercentile("CL1"));
-if(fSampleType=="pPb" || fSampleType=="PbPb") fHistCentStats->Fill(7.,centralityObj->GetCentralityPercentile("ZNA"));
+ fHistCentStats->Fill(6.,centralityObj->GetCentralityPercentile("CL1"));
+ fHistCentStats->Fill(7.,centralityObj->GetCentralityPercentile("ZNA"));
- cent_v0 = centralityObj->GetCentralityPercentile(fCentralityMethod);
- }
+ cent_v0 = centralityObj->GetCentralityPercentile(fCentralityMethod);
+ }
else cent_v0= -1;
+ }
+ else shift_to_TRACKS_MANUAL=kTRUE;
+
}//centralitymethod condition
- else if(fCentralityMethod=="V0M_MANUAL" || fCentralityMethod=="V0A_MANUAL" || fCentralityMethod=="V0C_MANUAL" || fCentralityMethod=="TRACKS_MANUAL")//data or RecoMc and also for TRUTH
+ else if(fCentralityMethod=="V0M_MANUAL" || fCentralityMethod=="V0A_MANUAL" || fCentralityMethod=="V0C_MANUAL" || fCentralityMethod=="TRACKS_MANUAL" || shift_to_TRACKS_MANUAL)//data or RecoMc and also for TRUTH
{
if(!truth){//for data or RecoMC
- cent_v0 = GetReferenceMultiplicityVZEROFromAOD(event);
+ cent_v0 = GetReferenceMultiplicityVZEROFromAOD((AliVEvent*)event);
}//for data or RecoMC
- if(truth){//condition for TRUTH case
+ if(truth && (fAnalysisType == "MCAOD")){//condition for TRUTH case
//check for TClonesArray(truth track MC information)
fArrayMC = dynamic_cast<TClonesArray*>(event->FindListObject(AliAODMCParticle::StdBranchName()));
if (!fArrayMC) {
if(fRemoveDuplicates && isduplicate) continue;//remove duplicates
- if (fCentralityMethod=="V0M_MANUAL") {
- if(partMC->Eta() > 5.1 || partMC->Eta() < 2.8) continue;
- if (partMC->Eta() < -3.7 || partMC->Eta() > -1.7) continue;
-}
- else if (fCentralityMethod=="V0A_MANUAL") {
- if(partMC->Eta() > 5.1 || partMC->Eta() < 2.8) continue;}
- else if (fCentralityMethod=="V0C_MANUAL") {
- if(partMC->Eta() > -1.7 || partMC->Eta() < -3.7) continue;}
- else if (fCentralityMethod=="TRACKS_MANUAL") {
- if (partMC->Eta() < fmineta || partMC->Eta() > fmaxeta) continue;
- if (partMC->Pt() < fminPt || partMC->Pt() > fmaxPt) continue;
+ // if (fCentralityMethod=="V0M_MANUAL")
+ if((partMC->Eta() < 5.1 && partMC->Eta() > 2.8) || (partMC->Eta() > -3.7 && partMC->Eta() < -1.7)) gRefMultiplicityVZERO_Truth+=1;
+ // else if (fCentralityMethod=="V0A_MANUAL") {
+ if(partMC->Eta() < 5.1 && partMC->Eta() > 2.8) gRefMultiplicityVZEROA_Truth+=1;
+ // else if (fCentralityMethod=="V0C_MANUAL") {
+ if(partMC->Eta() < -1.7 && partMC->Eta() > -3.7) gRefMultiplicityVZEROC_Truth+=1;
+ //else if (fCentralityMethod=="TRACKS_MANUAL") {
+ if (partMC->Eta() > fmineta && partMC->Eta() < fmaxeta) {
+ if (partMC->Pt() > fminPt && partMC->Pt() < fmaxPt) gRefMultiplicityTPC_Truth+=1;
}
- else{//basically returns the tracks manual case
-//give only kinematic cuts at the truth level
- if (partMC->Eta() < fmineta || partMC->Eta() > fmaxeta) continue;
- if (partMC->Pt() < fminPt || partMC->Pt() > fmaxPt) continue;
- }
+
+ }//truth track loop ends
- //To determine multiplicity in case of PP
- nooftrackstruth+= 1;;
+ fHistRefmult->Fill(3.,gRefMultiplicityTPC_Truth);
+ fHistRefmult->Fill(2.,gRefMultiplicityVZERO_Truth);
+ fHistRefmult->Fill(0.,gRefMultiplicityVZEROA_Truth);
+ fHistRefmult->Fill(1.,gRefMultiplicityVZEROC_Truth);
- }//truth track loop ends
- cent_v0=nooftrackstruth;
+ if(fCentralityMethod == "TRACKS_MANUAL") cent_v0=gRefMultiplicityTPC_Truth;
+
+ else if(fCentralityMethod == "V0M_MANUAL") cent_v0=gRefMultiplicityVZERO_Truth;
+
+ else if(fCentralityMethod == "V0A_MANUAL") cent_v0=gRefMultiplicityVZEROA_Truth;
+
+ else if(fCentralityMethod == "V0C_MANUAL") cent_v0=gRefMultiplicityVZEROC_Truth;
+
+ else cent_v0=gRefMultiplicityTPC_Truth;
}//condition for TRUTH case
}//end of MANUAL method
- else if ((fAnalysisType == "MCAOD") && (fCentralityMethod == "MC_b"))//TRUTH MC
+ else if ((fAnalysisType == "MCAOD") && (fCentralityMethod == "MC_b"))//TRUTH MC in AOD production
{
AliAODMCHeader* header = (AliAODMCHeader*) event->GetList()->FindObject(AliAODMCHeader::StdBranchName());
if (!header)
AliCollisionGeometry* collGeometry = dynamic_cast<AliCollisionGeometry*> (eventHeader);
- if (collGeometry) cent_v0 = collGeometry->ImpactParameter();
+ if (collGeometry) {
+ cent_v0 = collGeometry->ImpactParameter();
+ fhistImpactParm->Fill(cent_v0);
+ }
else cent_v0=-1.;
}//end of Impact parameter method
//else return -1
else cent_v0=-1.;
+}//AOD OR MCAOD condition
+
+
+else if(fAnalysisType == "MC"){
+ Double_t gImpactParameter = -1.;
+ AliMCEvent *gMCEvent = dynamic_cast<AliMCEvent*>(mainevent);
+ if(gMCEvent){
+ AliCollisionGeometry* headerH = dynamic_cast<AliCollisionGeometry*>(gMCEvent->GenEventHeader());
+ if(headerH){
+ gImpactParameter = headerH->ImpactParameter();
+
+ for(Int_t iParticle = 0; iParticle < gMCEvent->GetNumberOfPrimaries(); iParticle++) {
+ AliMCParticle* track = dynamic_cast<AliMCParticle *>(gMCEvent->GetTrack(iParticle));
+ if (!track) {
+ AliError(Form("Could not receive particle %d", iParticle));
+ continue;
+ }
+
+ //exclude non stable particles
+ if(!(gMCEvent->IsPhysicalPrimary(iParticle))) continue;
+
+ if(track->Charge() == 0) continue;
+
+ // if (fCentralityMethod=="V0M_MANUAL")
+ if((track->Eta() < 5.1 && track->Eta() > 2.8) || (track->Eta() > -3.7 && track->Eta() < -1.7)) gRefMultiplicityVZERO_Truth+=1;
+ // else if (fCentralityMethod=="V0A_MANUAL") {
+ if(track->Eta() < 5.1 && track->Eta() > 2.8) gRefMultiplicityVZEROA_Truth+=1;
+ // else if (fCentralityMethod=="V0C_MANUAL") {
+ if(track->Eta() < -1.7 && track->Eta() > -3.7) gRefMultiplicityVZEROC_Truth+=1;
+ //else if (fCentralityMethod=="TRACKS_MANUAL") {
+ if (track->Eta() > fmineta && track->Eta() < fmaxeta) {
+ if (track->Pt() > fminPt && track->Pt() < fmaxPt) gRefMultiplicityTPC_Truth+=1;}
+
+ }//loop over primaries
+
+ fHistRefmult->Fill(3.,gRefMultiplicityTPC_Truth);
+ fHistRefmult->Fill(2.,gRefMultiplicityVZERO_Truth);
+ fHistRefmult->Fill(0.,gRefMultiplicityVZEROA_Truth);
+ fHistRefmult->Fill(1.,gRefMultiplicityVZEROC_Truth);
+ if (fCentralityMethod == "MC_b"){
+ cent_v0=gImpactParameter;
+ fhistImpactParm->Fill(gImpactParameter);
+ fhistImpactParmvsMult->Fill(gImpactParameter,gRefMultiplicityTPC_Truth);
+ }
+ else if(fCentralityMethod == "TRACKS_MANUAL") cent_v0=gRefMultiplicityTPC_Truth;
+
+ else if(fCentralityMethod == "V0M_MANUAL") cent_v0=gRefMultiplicityVZERO_Truth;
+
+ else if(fCentralityMethod == "V0A_MANUAL") cent_v0=gRefMultiplicityVZEROA_Truth;
+
+ else if(fCentralityMethod == "V0C_MANUAL") cent_v0=gRefMultiplicityVZEROC_Truth;
+
+ else cent_v0=gImpactParameter;//default value is the impact parameter
+ }//MC event header
+ }//MC event cast
+ else cent_v0 = -1.;
+ }//MC condition
+
+ else{
+ cent_v0 = -1.;
+ }
return cent_v0;
}
//-----------------------------------------------------------------------------------------
-Double_t AliTwoParticlePIDCorr::GetAcceptedEventMultiplicity(AliAODEvent *aod,Bool_t truth){
-
- if(!aod) return -1;
+Double_t AliTwoParticlePIDCorr::GetAcceptedEventMultiplicity(AliVEvent *event,Bool_t truth){
+ //do the event selection(zvtx, pileup, centrality/multiplicity cut) and then return the value of the centrality of that event
+ if(!event) return -1;
Float_t gRefMultiplicity = -1.;
- // check first event in chunk (is not needed for new reconstructions)
+ //***********************************SOURCE CODE-TASKBFPsi
+
+ // Event Vertex MC
+ if(fAnalysisType == "MC") {
+ AliMCEvent *mcevent = dynamic_cast<AliMCEvent*>(event);
+ if(!mcevent) {
+ AliError("mcEvent not available");
+ return -1.;
+ }
+
+ if(mcevent){
+ AliGenEventHeader *header = dynamic_cast<AliGenEventHeader*>(mcevent->GenEventHeader());
+ if(header){
+ TArrayF gVertexArray;
+ header->PrimaryVertex(gVertexArray);
+ //count events having a proper vertex
+ fEventCounter->Fill(5);
+
+fHistQA[0]->Fill((gVertexArray.At(0)));fHistQA[1]->Fill((gVertexArray.At(1)));fHistQA[2]->Fill((gVertexArray.At(2))); //for trkVtx only before vertex cut |zvtx|<10 cm
+
+ if(TMath::Abs(gVertexArray.At(0)) < fVxMax_MC) {
+ if(TMath::Abs(gVertexArray.At(1)) < fVyMax_MC) {
+ if(TMath::Abs(gVertexArray.At(2)) < fVzMax_MC) {
+//count events after vertex cut
+ fEventCounter->Fill(7);
+ fHistQA[3]->Fill((gVertexArray.At(0)));fHistQA[4]->Fill((gVertexArray.At(1)));fHistQA[5]->Fill((gVertexArray.At(2)));//after vertex cut,for trkVtx only
+
+ // get the reference multiplicty or centrality
+ gRefMultiplicity = GetRefMultiOrCentrality((AliVEvent*)mcevent,kFALSE);//2nd argument has no meaning
+
+ if(gRefMultiplicity<0) return -1;
+
+ // take events only within the multiplicity class mentioned in the custom binning
+ if(gRefMultiplicity < fmincentmult || gRefMultiplicity > fmaxcentmult) return -1;
+
+//count events having proper centrality/ref multiplicity
+ fEventCounter->Fill(9);
+
+ }//Vz cut
+ }//Vy cut
+ }//Vx cut
+ }//header
+ }//MC event object
+ }//MC
+
+ else if(fAnalysisType == "MCAOD" || fAnalysisType == "AOD"){// if(fAnalysisType == "MCAOD" || fAnalysisType == "AOD"
+ //vertex selection(is it fine for PP?)
+ AliAODEvent* aod = dynamic_cast<AliAODEvent*>(event);
+
+
+ // check first event in chunk (is not needed for new reconstructions)
if(fCheckFirstEventInChunk){
AliAnalysisUtils ut;
if(ut.IsFirstEventInChunk(aod))
//count events after pileup selection
fEventCounter->Fill(3);
- //vertex selection(is it fine for PP?)
- if ( fVertextype==1){//for pPb basically if(!fAnalysisUtils->IsVertexSelected2013pA(aod)) return;
- trkVtx = aod->GetPrimaryVertex();
+ if (fVertextype==1){//for pPb basically if(!fAnalysisUtils->IsVertexSelected2013pA(aod)) return;
+ trkVtx = aod->GetPrimaryVertex();
if (!trkVtx || trkVtx->GetNContributors()<=0) return -1;
TString vtxTtl = trkVtx->GetTitle();
if (!vtxTtl.Contains("VertexerTracks")) return -1;
}
else if(fVertextype==0){//default case
- trkVtx = aod->GetPrimaryVertex();
+ trkVtx =(AliAODVertex*) aod->GetPrimaryVertex();
if (!trkVtx || trkVtx->GetNContributors()<=0) return -1;//proper number of contributors
zvtx = trkVtx->GetZ();
Double32_t fCov[6];
fHistQA[3]->Fill((trkVtx->GetX()));fHistQA[4]->Fill((trkVtx->GetY()));fHistQA[5]->Fill((trkVtx->GetZ()));//after vertex cut,for trkVtx only
//get the centrality or multiplicity
- if(truth) {gRefMultiplicity = GetRefMultiOrCentrality(aod,kTRUE);}//kTRUE-->for Truth case(only meaningful in case of ref multiplicity)
+ if(truth) {gRefMultiplicity = GetRefMultiOrCentrality((AliVEvent*)aod,kTRUE);}//kTRUE-->for Truth case(only meaningful in case of ref multiplicity,in case of centrality it has no meaning)
- else {gRefMultiplicity = GetRefMultiOrCentrality(aod,kFALSE);}//kFALSE-->for data and RecoMc case(only meaningful in case of ref multiplicity)
+ else {gRefMultiplicity = GetRefMultiOrCentrality((AliVEvent*)aod,kFALSE);}//kFALSE-->for data and RecoMc case(only meaningful in case of ref multiplicity,in case of centrality it has no meaning)
if(gRefMultiplicity<0) return -1;
fEventCounter->Fill(9);
-// centrality weighting (optional for 2011 if central and semicentral triggers are used);only for data and recoMC
- if (fCentralityWeights && !AcceptEventCentralityWeight(gRefMultiplicity))//**********************
- {
- AliInfo(Form("Rejecting event because of centrality weighting: %f", gRefMultiplicity));
- return -1;
- }
+// centrality weighting (optional for 2011 if central and semicentral triggers are used);only for data and recoMC
+ if (fCentralityWeights && !AcceptEventCentralityWeight(gRefMultiplicity))//**********************
+ {
+ AliInfo(Form("Rejecting event because of centrality weighting: %f", gRefMultiplicity));
+ return -1;
+ }
+
+//count events after rejection due to centrality weighting
+ fEventCounter->Fill(11);
+
+ }
+ else gRefMultiplicity=-1;
+
+ return gRefMultiplicity;
+
+}
+//--------------------------------------------------------------------------------------------------------
+Float_t AliTwoParticlePIDCorr::GetEventPlane(AliVEvent *mainevent,Bool_t truth, Double_t v0Centr)
+{
+ Float_t eventplane=999.;
+ // Get the event plane
+ if(!mainevent) return 999.;
+
+
+ //MC: from reaction plane
+ if(fAnalysisType == "MC"){
+ if(!mainevent) {
+ AliError("mcEvent not available");
+ return 999.;
+ }
+
+ AliMCEvent *gMCEvent = dynamic_cast<AliMCEvent*>(mainevent);
+ if(gMCEvent){
+ AliCollisionGeometry* headerH = dynamic_cast<AliCollisionGeometry*>(gMCEvent->GenEventHeader());
+ if (headerH) {
+ Int_t iC = -1;
+ // Impact parameter bins(it is only for Pb-Pb)
+ if(v0Centr < 3.50) iC = 0;
+ else if(v0Centr < 4.94) iC = 1;
+ else if(v0Centr < 6.98) iC = 2;
+ else if(v0Centr < 8.55) iC = 3;
+ else if(v0Centr < 9.88) iC = 4;
+ else if(v0Centr < 11.04) iC = 5;
+ else if(v0Centr < 12.09) iC = 6;
+ else if(v0Centr < 13.05) iC = 7;
+ else iC = 8;
+
+ eventplane = headerH->ReactionPlaneAngle();
+ if(eventplane > TMath::Pi()/2 && eventplane <= TMath::Pi()*3/2) eventplane-=TMath::Pi();
+ if(eventplane > TMath::Pi()*3/2) eventplane-=2*TMath::Pi();
+ fHistEventPlaneTruth->Fill(iC,eventplane);
+ //gReactionPlane *= TMath::RadToDeg();
+ }//MC header
+ }//MC event cast
+ }//MC
+
+ else if(fAnalysisType == "MCAOD" || fAnalysisType == "AOD") {
+ //reset Q vector info
+
+ AliAODEvent* event = dynamic_cast<AliAODEvent*>(mainevent);
+
+
+ Int_t run = event->GetRunNumber();
+
+ if(run != fRun){
+ // Load the calibrations run dependent
+ if(! fIsAfter2011) OpenInfoCalbration(run);
+ fRun=run;
+ }
+
-//count events after rejection due to centrality weighting
- fEventCounter->Fill(11);
+ Int_t iC = -1;
+ if (v0Centr > 80) return 999.; // analysis only for 0-80% centrality classes
+ // centrality bins
+ if(v0Centr < 5) iC = 0;
+ else if(v0Centr < 10) iC = 1;
+ else if(v0Centr < 20) iC = 2;
+ else if(v0Centr < 30) iC = 3;
+ else if(v0Centr < 40) iC = 4;
+ else if(v0Centr < 50) iC = 5;
+ else if(v0Centr < 60) iC = 6;
+ else if(v0Centr < 70) iC = 7;
+ else iC = 8;
- return gRefMultiplicity;
-}
-//--------------------------------------------------------------------------------------------------------
-Double_t AliTwoParticlePIDCorr::GetEventPlane(AliAODEvent *event,Bool_t truth){
- // Get the event plane
+ Int_t iCcal = iC;
+ //reset Q vector info
+ Double_t Qxa2 = 0, Qya2 = 0;
+ Double_t Qxc2 = 0, Qyc2 = 0;
+ Double_t Qxa3 = 0, Qya3 = 0;
+ Double_t Qxc3 = 0, Qyc3 = 0;
- Float_t gVZEROEventPlane = -10.;
- Float_t gReactionPlane = -10.;
- Double_t qxTot = 0.0, qyTot = 0.0;
//MC: from reaction plane
if(truth)
{
AliAODMCHeader* header = (AliAODMCHeader*) event->GetList()->FindObject(AliAODMCHeader::StdBranchName());
if (header){
-
- AliGenEventHeader* eventHeader = header->GetCocktailHeader(0); // get first MC header from either ESD/AOD (including cocktail header if available)
+ evplaneMC = header->GetReactionPlaneAngle();//[0, 360]
+ //make it within [-pi/2,pi/2] to make it general
+ if(evplaneMC > TMath::Pi()/2 && evplaneMC <= TMath::Pi()*3/2) evplaneMC-=TMath::Pi();
+ if(evplaneMC > TMath::Pi()*3/2) evplaneMC-=2*TMath::Pi();
+ fHistEventPlaneTruth->Fill(iC,evplaneMC);
+ /*
+ AliGenEventHeader* eventHeader = header->GetCocktailHeader(0); // get first MC header from either ESD/AOD (including cocktail header if available)
if (eventHeader)
{
AliCollisionGeometry* collGeometry = dynamic_cast<AliCollisionGeometry*> (eventHeader);
- if (collGeometry) gReactionPlane = collGeometry->ReactionPlaneAngle();
- }
+ if (collGeometry){//get the reaction plane from MC header
+ gReactionPlane = collGeometry->ReactionPlaneAngle();//[0,180]
+ }
+ }
+ */
+ //taken from vnv0 code(get the TPC, V0A, V0C event plane using truth tracks)
+ TClonesArray *mcArray = NULL;
+ mcArray = (TClonesArray*)event->GetList()->FindObject(AliAODMCParticle::StdBranchName());
+ if(mcArray){
+ Float_t QxMCv2[3] = {0,0,0};
+ Float_t QyMCv2[3] = {0,0,0};
+ Float_t QxMCv3[3] = {0,0,0};
+ Float_t QyMCv3[3] = {0,0,0};
+ Float_t EvPlaneMCV2[3] = {0,0,0};
+ Float_t EvPlaneMCV3[3] = {0,0,0};
+ Float_t etaMin[3] = {2.8,-3.6,-0.8}; // A-side, C-side M-barrel
+ Float_t etaMax[3] = {4.88,-1.8,0.8};
+
+ // analysis on MC tracks
+ Int_t nMCtrack = mcArray->GetEntries() ;
+
+ // EP computation with MC tracks
+ for(Int_t iT=0;iT < nMCtrack;iT++){
+ AliAODMCParticle *mctr = (AliAODMCParticle*) mcArray->At(iT);
+ if(!mctr || !(mctr->IsPrimary()) || !(mctr->Charge()) || mctr->Pt() < 0.2) continue;
+
+ Float_t eta = mctr->Eta();
+ for(Int_t iD=0;iD<3;iD++){
+ if(eta > etaMin[iD] && eta < etaMax[iD]){
+ Float_t phi = mctr->Phi();
+ QxMCv2[iD] += TMath::Cos(2*phi);
+ QyMCv2[iD] += TMath::Sin(2*phi);
+ QxMCv3[iD] += TMath::Cos(3*phi);
+ QyMCv3[iD] += TMath::Sin(3*phi);
+ }
+ }
+ }
+
+ EvPlaneMCV2[0] = TMath::ATan2(QyMCv2[0],QxMCv2[0])/2.;
+ EvPlaneMCV2[1] = TMath::ATan2(QyMCv2[1],QxMCv2[1])/2.;
+ EvPlaneMCV2[2] = TMath::ATan2(QyMCv2[2],QxMCv2[2])/2.;
+ fHResMA2->Fill(Double_t(iC), TMath::Cos(2*(EvPlaneMCV2[2]-EvPlaneMCV2[0])));
+ fHResMC2->Fill(Double_t(iC), TMath::Cos(2*(EvPlaneMCV2[2]-EvPlaneMCV2[1])));
+ fHResAC2->Fill(Double_t(iC), TMath::Cos(2*(EvPlaneMCV2[0]-EvPlaneMCV2[1])));
+ fgPsi2v0aMC = EvPlaneMCV2[0];
+ fgPsi2v0cMC = EvPlaneMCV2[1];
+ fgPsi2tpcMC = EvPlaneMCV2[2];
+
+
+ EvPlaneMCV3[0] = TMath::ATan2(QyMCv3[0],QxMCv3[0])/3.;
+ EvPlaneMCV3[1] = TMath::ATan2(QyMCv3[1],QxMCv3[1])/3.;
+ EvPlaneMCV3[2] = TMath::ATan2(QyMCv3[2],QxMCv3[2])/3.;
+ fHResMA3->Fill(Double_t(iC), TMath::Cos(3*(EvPlaneMCV3[2]-EvPlaneMCV3[0])));
+ fHResMC3->Fill(Double_t(iC), TMath::Cos(3*(EvPlaneMCV3[2]-EvPlaneMCV3[1])));
+ fHResAC3->Fill(Double_t(iC), TMath::Cos(3*(EvPlaneMCV3[0]-EvPlaneMCV3[1])));
+ fgPsi3v0aMC = EvPlaneMCV3[0];
+ fgPsi3v0cMC = EvPlaneMCV3[1];
+ fgPsi3tpcMC = EvPlaneMCV3[2];
+
+ }
}
+
}
else{
-
+ Int_t nAODTracks = event->GetNumberOfTracks();
+
+// TPC EP needed for resolution studies (TPC subevent)
+ //AliEventplane * ep = (fAOD->GetHeader())->GetEventplaneP();
+ //Double_t psiTPC = ep->GetEventplane("Q", fAOD, 2); // in range of [0, pi]
+ Double_t Qx2 = 0, Qy2 = 0;
+ Double_t Qx3 = 0, Qy3 = 0;
+
+ for(Int_t iT = 0; iT < nAODTracks; iT++) {
+
+ AliAODTrack* aodTrack =(AliAODTrack*) event->GetTrack(iT);
+
+ if (!aodTrack){
+ continue;
+ }
+
+ Bool_t trkFlag = aodTrack->TestFilterBit(1);
+
+ if ((TMath::Abs(aodTrack->Eta()) > 0.8) || (aodTrack->Pt() < 0.2) || (aodTrack->GetTPCNcls() < fNcluster) || !trkFlag)
+ continue;
+
+ Double_t b[2] = {-99., -99.};
+ Double_t bCov[3] = {-99., -99., -99.};
+
+
+ AliAODTrack param(*aodTrack);
+ if (!param.PropagateToDCA(event->GetPrimaryVertex(), event->GetMagneticField(), 100., b, bCov)){
+ continue;
+ }
+
+ if ((TMath::Abs(b[0]) > 3.0) || (TMath::Abs(b[1]) > 2.4))
+ continue;
+
+ Qx2 += TMath::Cos(2*aodTrack->Phi());
+ Qy2 += TMath::Sin(2*aodTrack->Phi());
+ Qx3 += TMath::Cos(3*aodTrack->Phi());
+ Qy3 += TMath::Sin(3*aodTrack->Phi());
+
+ }
+
+ Float_t evPlAng2 = TMath::ATan2(Qy2, Qx2)/2.;
+ Float_t evPlAng3 = TMath::ATan2(Qy3, Qx3)/3.;
+
+ fgPsi2tpc = evPlAng2;
+ fgPsi3tpc = evPlAng3;
+
+ fPhiRPTPC->Fill(iC,evPlAng2);
+ fPhiRPTPCv3->Fill(iC,evPlAng3);
+
+
+
+//V0 info
+ AliAODVZERO* aodV0 = event->GetVZEROData();
+
+ for (Int_t iv0 = 0; iv0 < 64; iv0++) {
+ Double_t phiV0 = TMath::PiOver4()*(0.5 + iv0 % 8);
+ Float_t multv0 = aodV0->GetMultiplicity(iv0);
+
+ if(! fIsAfter2011){
+ if(fAnalysisType == "AOD"){//not for reco MC tracks, only for real data
+ if (iv0 < 32){ // V0C
+ Qxc2 += TMath::Cos(2*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1);
+ Qyc2 += TMath::Sin(2*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1);
+ Qxc3 += TMath::Cos(3*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1);
+ Qyc3 += TMath::Sin(3*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1);
+ } else { // V0A
+ Qxa2 += TMath::Cos(2*phiV0) * multv0*fV0Apol/fMultV0->GetBinContent(iv0+1);
+ Qya2 += TMath::Sin(2*phiV0) * multv0*fV0Apol/fMultV0->GetBinContent(iv0+1);
+ Qxa3 += TMath::Cos(3*phiV0) * multv0*fV0Apol/fMultV0->GetBinContent(iv0+1);
+ Qya3 += TMath::Sin(3*phiV0) * multv0*fV0Apol/fMultV0->GetBinContent(iv0+1);
+ }
+ }
+ else{
+ if (iv0 < 32){ // V0C
+ Qxc2 += TMath::Cos(2*phiV0) * multv0;//*fV0Cpol/fMultV0->GetBinContent(iv0+1);
+ Qyc2 += TMath::Sin(2*phiV0) * multv0;//*fV0Cpol/fMultV0->GetBinContent(iv0+1);
+ Qxc3 += TMath::Cos(3*phiV0) * multv0;//*fV0Cpol/fMultV0->GetBinContent(iv0+1);
+ Qyc3 += TMath::Sin(3*phiV0) * multv0;//*fV0Cpol/fMultV0->GetBinContent(iv0+1);
+ } else { // V0A
+ Qxa2 += TMath::Cos(2*phiV0) * multv0;//*fV0Apol/fMultV0->GetBinContent(iv0+1);
+ Qya2 += TMath::Sin(2*phiV0) * multv0;//*fV0Apol/fMultV0->GetBinContent(iv0+1);
+ Qxa3 += TMath::Cos(3*phiV0) * multv0;//*fV0Apol/fMultV0->GetBinContent(iv0+1);
+ Qya3 += TMath::Sin(3*phiV0) * multv0;//*fV0Apol/fMultV0->GetBinContent(iv0+1);
+ }
+ }
+ }
+ }
+ //grab for each centrality the proper histo with the Qx and Qy to do the recentering
+ Double_t Qxamean2 = fMeanQ[iCcal][1][0];
+ Double_t Qxarms2 = fWidthQ[iCcal][1][0];
+ Double_t Qyamean2 = fMeanQ[iCcal][1][1];
+ Double_t Qyarms2 = fWidthQ[iCcal][1][1];
+ Double_t Qxamean3 = fMeanQv3[iCcal][1][0];
+ Double_t Qxarms3 = fWidthQv3[iCcal][1][0];
+ Double_t Qyamean3 = fMeanQv3[iCcal][1][1];
+ Double_t Qyarms3 = fWidthQv3[iCcal][1][1];
+
+ Double_t Qxcmean2 = fMeanQ[iCcal][0][0];
+ Double_t Qxcrms2 = fWidthQ[iCcal][0][0];
+ Double_t Qycmean2 = fMeanQ[iCcal][0][1];
+ Double_t Qycrms2 = fWidthQ[iCcal][0][1];
+ Double_t Qxcmean3 = fMeanQv3[iCcal][0][0];
+ Double_t Qxcrms3 = fWidthQv3[iCcal][0][0];
+ Double_t Qycmean3 = fMeanQv3[iCcal][0][1];
+ Double_t Qycrms3 = fWidthQv3[iCcal][0][1];
+
+ Double_t QxaCor2 = (Qxa2 - Qxamean2)/Qxarms2;
+ Double_t QyaCor2 = (Qya2 - Qyamean2)/Qyarms2;
+ Double_t QxcCor2 = (Qxc2 - Qxcmean2)/Qxcrms2;
+ Double_t QycCor2 = (Qyc2 - Qycmean2)/Qycrms2;
+ Double_t QxaCor3 = (Qxa3 - Qxamean3)/Qxarms3;
+ Double_t QyaCor3 = (Qya3 - Qyamean3)/Qyarms3;
+ Double_t QxcCor3 = (Qxc3 - Qxcmean3)/Qxcrms3;
+ Double_t QycCor3 = (Qyc3 - Qycmean3)/Qycrms3;
+ /*
+ //to calculate 2nd order event plane with v0M
+ Double_t QxCor2 = (Qxa2 - Qxamean2 + Qxc2 - Qxcmean2)
+ /TMath::Sqrt(Qxarms2*Qxarms2 + Qxcrms2*Qxcrms2);
+ Double_t QyCor2 = (Qya2 - Qyamean2 + Qyc2 - Qycmean2)
+ /TMath::Sqrt(Qyarms2*Qyarms2 + Qycrms2*Qycrms2);
+
+ //here the calculated event plane is within -Pi to +Pi(delete it , no use here , only for definition)
+ Double_t psiV0A =(TMath::Pi() + TMath::ATan2(-QyaCor2, -QxaCor2))/2.;
+ Double_t psiV0C = (TMath::Pi() + TMath::ATan2(-QycCor2, -QxcCor2))/2.;
+ Double_t psiVZero = (TMath::Pi() + TMath::ATan2(-QyCor2, -QxCor2))/2.;
+
+ */
+
+ Float_t evPlAngV0ACor2=999.;
+ Float_t evPlAngV0CCor2=999.;
+ Float_t evPlAngV0ACor3=999.;
+ Float_t evPlAngV0CCor3=999.;
+
+ if(! fIsAfter2011){
+ if(fAnalysisType == "AOD"){
+ evPlAngV0ACor2 = TMath::ATan2(QyaCor2, QxaCor2)/2.;
+ evPlAngV0CCor2 = TMath::ATan2(QycCor2, QxcCor2)/2.;
+ evPlAngV0ACor3 = TMath::ATan2(QyaCor3, QxaCor3)/3.;
+ evPlAngV0CCor3 = TMath::ATan2(QycCor3, QxcCor3)/3.;
+ }
+ else{
+ evPlAngV0ACor2 = TMath::ATan2(Qya2, Qxa2)/2.;
+ evPlAngV0CCor2 = TMath::ATan2(Qyc2, Qxc2)/2.;
+ evPlAngV0ACor3 = TMath::ATan2(Qya3, Qxa3)/3.;
+ evPlAngV0CCor3 = TMath::ATan2(Qyc3, Qxc3)/3.;
+ }
+ }
+ else{
+ AliEventplane *ep = event->GetEventplane();
+ evPlAngV0ACor2 = ep->GetEventplane("V0A", event, 2);
+ evPlAngV0CCor2 = ep->GetEventplane("V0C", event, 2);
+ evPlAngV0ACor3 = ep->GetEventplane("V0A", event, 3);
+ evPlAngV0CCor3 = ep->GetEventplane("V0C", event, 3);
+ }
+
+ fgPsi2v0a = evPlAngV0ACor2;
+ fgPsi2v0c = evPlAngV0CCor2;
+ fgPsi3v0a = evPlAngV0ACor3;
+ fgPsi3v0c = evPlAngV0CCor3;
+
+ // Fill EP distribution histograms evPlAng
+
+ fPhiRPv0A->Fill(iC,evPlAngV0ACor2);
+ fPhiRPv0C->Fill(iC,evPlAngV0CCor2);
+
+ fPhiRPv0Av3->Fill(iC,evPlAngV0ACor3);
+ fPhiRPv0Cv3->Fill(iC,evPlAngV0CCor3);
+
+ // Fill histograms needed for resolution evaluation
+ fHResTPCv0A2->Fill(Double_t(iC), TMath::Cos(2*(evPlAng2 - evPlAngV0ACor2)));
+ fHResTPCv0C2->Fill(Double_t(iC), TMath::Cos(2*(evPlAng2 - evPlAngV0CCor2)));
+ fHResv0Cv0A2->Fill(Double_t(iC), TMath::Cos(2*(evPlAngV0ACor2 - evPlAngV0CCor2)));
+
+ fHResTPCv0A3->Fill(Double_t(iC), TMath::Cos(3*(evPlAng3 - evPlAngV0ACor3)));
+ fHResTPCv0C3->Fill(Double_t(iC), TMath::Cos(3*(evPlAng3 - evPlAngV0CCor3)));
+ fHResv0Cv0A3->Fill(Double_t(iC), TMath::Cos(3*(evPlAngV0ACor3 - evPlAngV0CCor3)));
+
+
+ /*
+ Float_t gVZEROEventPlane = -10.;
+ Float_t gReactionPlane = -10.;
+ Double_t qxTot = 0.0, qyTot = 0.0;
+
AliEventplane *ep = event->GetEventplane();
if(ep){
gVZEROEventPlane = ep->CalculateVZEROEventPlane(event,10,2,qxTot,qyTot);
//gReactionPlane = gVZEROEventPlane*TMath::RadToDeg();
gReactionPlane = gVZEROEventPlane;
}
+ */
}//AOD,ESD,ESDMC
- return gReactionPlane;
+ //return gReactionPlane;
+
+ //make the final 2nd order event plane within 0 to Pi
+ //using data and reco tracks only
+ if(fgPsi2v0a!=999. && fgPsi2v0a < 0.) fgPsi2v0a += TMath::Pi();
+ if(fgPsi2v0c!=999. && fgPsi2v0c < 0.) fgPsi2v0c += TMath::Pi();
+ if(fgPsi2tpc!=999. && fgPsi2tpc < 0.) fgPsi2tpc += TMath::Pi();
+ //using truth tracks only
+ if(evplaneMC!=999. && evplaneMC < 0.) evplaneMC += TMath::Pi();
+ if(fgPsi2v0aMC!=999. && fgPsi2v0aMC < 0.) fgPsi2v0aMC += TMath::Pi();
+ if(fgPsi2v0cMC!=999. && fgPsi2v0cMC < 0.) fgPsi2v0cMC += TMath::Pi();
+ if(fgPsi2tpcMC!=999. && fgPsi2tpcMC < 0.) fgPsi2tpcMC += TMath::Pi();
+ //for the time being leave the 3rd order event planes within -pi/3 t0 +pi/3
+
+ if(truth){//for truth MC
+ if(fV2 && fEPdet=="header")eventplane=evplaneMC;
+ if(fV2 && fEPdet=="V0A")eventplane=fgPsi2v0aMC;
+ if(fV2 && fEPdet=="V0C")eventplane=fgPsi2v0cMC;
+ if(fV2 && fEPdet=="TPC")eventplane=fgPsi2tpcMC;
+
+ if(fV3 && fEPdet=="V0A")eventplane=fgPsi3v0aMC;
+ if(fV3 && fEPdet=="V0C")eventplane=fgPsi3v0cMC;
+ if(fV3 && fEPdet=="TPC")eventplane=fgPsi3tpcMC;
}
+ else{//for data and recoMC
+ if(fV2 && fEPdet=="V0A")eventplane=fgPsi2v0a;
+ if(fV2 && fEPdet=="V0C")eventplane=fgPsi2v0c;
+ if(fV2 && fEPdet=="TPC")eventplane=fgPsi2tpc;
+ if(fV3 && fEPdet=="V0A")eventplane=fgPsi3v0a;
+ if(fV3 && fEPdet=="V0C")eventplane=fgPsi3v0c;
+ if(fV3 && fEPdet=="TPC")eventplane=fgPsi3tpc;
+ }
+
+
+ }//AOD/MCAOD condition
+
+ else eventplane=999.;
+
+ return eventplane;
+}
+//------------------------------------------------------------------------------------------------------------------
+void AliTwoParticlePIDCorr::OpenInfoCalbration(Int_t run){
+ TString oadbfilename = "$ALICE_ROOT/OADB/PWGCF/VZERO/VZEROcalibEP.root";
+ TFile *foadb = TFile::Open(oadbfilename.Data());
+
+ if(!foadb){
+ printf("OADB file %s cannot be opened\n",oadbfilename.Data());
+ return;
+ }
+ AliOADBContainer *cont = (AliOADBContainer*) foadb->Get("hMultV0BefCorr");
+ if(!cont){
+ printf("OADB object hMultV0BefCorr is not available in the file\n");
+ return;
+ }
+ if(!(cont->GetObject(run))){
+ printf("OADB object hMultV0BefCorr is not available for run %i (used run 137366)\n",run);
+ run = 137366;
+ }
+ fMultV0 = ((TH2F *) cont->GetObject(run))->ProfileX();
+
+ TF1 *fpol0 = new TF1("fpol0","pol0");
+ fMultV0->Fit(fpol0,"","",0,31);
+ fV0Cpol = fpol0->GetParameter(0);
+ fMultV0->Fit(fpol0,"","",32,64);
+ fV0Apol = fpol0->GetParameter(0);
+
+ for(Int_t iside=0;iside<2;iside++){
+ for(Int_t icoord=0;icoord<2;icoord++){
+ for(Int_t i=0;i < 9;i++){
+ char namecont[100];
+ if(iside==0 && icoord==0)
+ snprintf(namecont,100,"hQxc2_%i",i);
+ else if(iside==1 && icoord==0)
+ snprintf(namecont,100,"hQxa2_%i",i);
+ else if(iside==0 && icoord==1)
+ snprintf(namecont,100,"hQyc2_%i",i);
+ else if(iside==1 && icoord==1)
+ snprintf(namecont,100,"hQya2_%i",i);
+
+ cont = (AliOADBContainer*) foadb->Get(namecont);
+ if(!cont){
+ printf("OADB object %s is not available in the file\n",namecont);
+ return;
+ }
+
+ if(!(cont->GetObject(run))){
+ printf("OADB object %s is not available for run %i (used run 137366)\n",namecont,run);
+ run = 137366;
+ }
+ fMeanQ[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetMean();
+ fWidthQ[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetRMS();
+
+ //for v3
+ if(iside==0 && icoord==0)
+ snprintf(namecont,100,"hQxc3_%i",i);
+ else if(iside==1 && icoord==0)
+ snprintf(namecont,100,"hQxa3_%i",i);
+ else if(iside==0 && icoord==1)
+ snprintf(namecont,100,"hQyc3_%i",i);
+ else if(iside==1 && icoord==1)
+ snprintf(namecont,100,"hQya3_%i",i);
+
+ cont = (AliOADBContainer*) foadb->Get(namecont);
+ if(!cont){
+ printf("OADB object %s is not available in the file\n",namecont);
+ return;
+ }
+
+ if(!(cont->GetObject(run))){
+ printf("OADB object %s is not available for run %i (used run 137366)\n",namecont,run);
+ run = 137366;
+ }
+ fMeanQv3[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetMean();
+ fWidthQv3[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetRMS();
+ }
+ }
+ }
+}
//____________________________________________________________________
-void AliTwoParticlePIDCorr::Terminate(Option_t *)
+void AliTwoParticlePIDCorr::FillPIDEventPlane(Double_t centrality,Int_t par,Float_t trigphi,Float_t fReactionPlane)
{
- // 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; }
-
-
+
+ // Event plane (determine psi bin)
+ Double_t gPsiMinusPhi = 0.;
+ Double_t gPsiMinusPhiBin = -10.;
+if(fRequestEventPlane){
+ gPsiMinusPhi = TMath::Abs(trigphi - fReactionPlane);
+ //in-plane
+ if((gPsiMinusPhi <= 7.5*TMath::DegToRad())||
+ (gPsiMinusPhi >= 352.5*TMath::DegToRad())||
+ ((172.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 187.5*TMath::DegToRad())))
+ gPsiMinusPhiBin = 0.0;
+ //intermediate
+ else if(((37.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 52.5*TMath::DegToRad()))||
+ ((127.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 142.5*TMath::DegToRad()))||
+ ((217.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 232.5*TMath::DegToRad()))||
+ ((307.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 322.5*TMath::DegToRad())))
+ gPsiMinusPhiBin = 1.0;
+ //out of plane
+ else if(((82.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 97.5*TMath::DegToRad()))||
+ ((262.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 277.5*TMath::DegToRad())))
+ gPsiMinusPhiBin = 2.0;
+ //everything else
+ else
+ gPsiMinusPhiBin = 3.0;
+
+ fEventPlanePID->Fill(centrality,gPsiMinusPhiBin,(Float_t)par);
+ }
}
-//------------------------------------------------------------------
-/*
- // get centrality object and check quality
- Double_t cent_v0=0;
+//____________________________________________________________________________________________________
+
+TObjArray* AliTwoParticlePIDCorr::GetV0Particles(AliVEvent* event,Double_t Centrality)
+{
+ AliAODEvent* fAOD = dynamic_cast<AliAODEvent*>(event);
-if(fCentralityMethod=="V0M" || fCentralityMethod=="V0A" || fCentralityMethod=="V0C")//for PbPb, pPb, pp7TeV(still to be introduced)
- {
- AliCentrality *centralityObj=0;
- if(aod)
- AliAODHeader *header = (AliAODHeader*) aod->GetHeader();
- if(header){
- centralityObj = header->GetCentralityP();
- // if (centrality->GetQuality() != 0) return ;
+ //function to select v0's from AODs
+ trkVtx=fAOD->GetPrimaryVertex();
+ Float_t xv=trkVtx->GetX(), yv=trkVtx->GetY(), zv=trkVtx->GetZ();
+ Int_t nV0sTot = fAOD->GetNumberOfV0s();
- if(centralityObj)
- {
-if(fSampleType=="pPb" || fSampleType=="PbPb" || fSampleType=="pp") fHistCentStats->Fill(0.,centralityObj->GetCentralityPercentile("V0M"));//only available for LHC10d at present (Quantile info)
-if(fSampleType=="pPb" || fSampleType=="PbPb") fHistCentStats->Fill(2.,centralityObj->GetCentralityPercentile("V0A"));
-if(fSampleType=="pPb" || fSampleType=="PbPb") fHistCentStats->Fill(4.,centralityObj->GetCentralityPercentile("V0C"));
-if(fSampleType=="pPb" || fSampleType=="PbPb") fHistCentStats->Fill(6.,centralityObj->GetCentralityPercentile("CL1"));
-if(fSampleType=="pPb" || fSampleType=="PbPb") fHistCentStats->Fill(8.,centralityObj->GetCentralityPercentile("ZNA"));
+ TObjArray * selectedV0s = new TObjArray;
+ selectedV0s->SetOwner(kTRUE);
- cent_v0 = centralityObj->GetCentralityPercentile(fCentralityMethod);
- }
- else
+ for (Int_t iV0 = 0; iV0 < nV0sTot; iV0++)
{
- cent_v0= -1;
- }
- }//AOD header
- }//centralitymethod condition
+
+ AliAODv0 *v0=fAOD->GetV0(iV0);
+ if (!v0) continue;
+ if(!CheckStatusv0(v0)) continue;
-else if(fCentralityMethod=="V0M_MANUAL" || fCentralityMethod=="V0A_MANUAL" || fCentralityMethod=="V0C_MANUAL" || fCentralityMethod=="TRACKS_MANUAL")
- {
- cent_v0 = GetReferenceMultiplicityVZEROFromAOD(aod);
- fHistrefMultiplicity->Fill(cent_v0);
- }
- else cent_v0= -1;
+ AliAODTrack *ptrack=(AliAODTrack*)v0->GetDaughter(0);
+ AliAODTrack *ntrack=(AliAODTrack*)v0->GetDaughter(1);
+
+ Bool_t cutK0sPID=kFALSE;
+ Bool_t cutLambdaPID=kFALSE;
+ Bool_t cutAntiLambdaPID=kFALSE;
+ if(fUsev0DaughterPID)
+{
+ //use fHelperPID check PID of the daughter tracks
+ //v0 daughter PID may be helpful in distangling k0S and (Anti)Lamda
+
+ Int_t PIDptrack = GetParticle(ptrack,kFALSE);
+ Int_t PIDntrack = GetParticle(ntrack ,kFALSE);
-if(fSampleType=="pPb" || fSampleType=="PbPb") if (cent_v0 < 0) return;//for pp case it is the multiplicity
+ if(PIDptrack ==0 && PIDntrack == 0) cutK0sPID=kTRUE;
+ if(PIDptrack==2 && PIDntrack ==0) cutLambdaPID=kTRUE;
-
+ if (PIDptrack ==0 && PIDntrack == 2) cutAntiLambdaPID=kTRUE;
- 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?.
+ // effective mass calculations for each hypothesis(without daughter PID)
+ Double_t InvMassK0s = v0->MassK0Short();
+ Double_t InvMassAntiLambda = v0->MassAntiLambda();
+ Double_t InvMassLambda = v0->MassLambda();
+
+ Float_t v0Pt=TMath::Sqrt(v0->Pt2V0());
+ Float_t v0Eta=v0->Eta();
+ Float_t v0Phi=v0->Phi();
+
+ //This is simply raw v0 without any specialised cut
+ fHistRawPtCentInvK0s->Fill(InvMassK0s,v0Pt,Centrality);
+ fHistRawPtCentInvLambda->Fill(InvMassLambda,v0Pt,Centrality);
+ fHistRawPtCentInvAntiLambda->Fill(InvMassAntiLambda,v0Pt,Centrality);
+
+ // Decay vertex
+ Double_t xyz[3];
+ v0->GetSecondaryVtx(xyz);
+ Float_t dx,dy,dz;
+ dx=xyz[0]-xv, dy=xyz[1]-yv, dz=xyz[2]-zv;
+
+ // Float_t v0DecayRadius=TMath::Sqrt(dx*dx + dy*dy);
+ Float_t v0DecayLength=TMath::Sqrt(dx*dx + dy*dy + dz*dz);
+ // VO's main characteristics to check the reconstruction cuts
+ // Float_t DcaV0Daughters = v0->DcaV0Daughters();
+ Float_t V0cosPointAngle = v0->CosPointingAngle(trkVtx);
+ // Float_t DcaPosToPrimVertex = v0->DcaPosToPrimVertex();
+ //Float_t DcaNegToPrimVertex = v0->DcaNegToPrimVertex();
+ //Float_t Dcav0PVz = v0->DcaV0ToPrimVertex();
+ Float_t v0Pz=v0->Pz();
+ Float_t v0P= TMath::Sqrt(v0Pt*v0Pt + v0Pz*v0Pz);
+
+ Float_t ctauLambda =999.;
+ Float_t ctauAntiLambda = 999.;
+ Float_t ctauK0s = 999.;
+ if(v0P > 0.0)
{
- if (fAnalysisUtils && fAnalysisUtils->IsPileUpEvent(aod)) return;
-//count events after PileUP cut
- fEventCounter->Fill(3);
- }
+ ctauLambda = (v0DecayLength*InvMassLambda)/v0P;
+ ctauAntiLambda = (v0DecayLength*InvMassAntiLambda)/v0P;
+ ctauK0s = (v0DecayLength*InvMassK0s)/v0P;
+ }
+
+ Bool_t ctauCutK0s= ctauK0s < NCtau*fCutctauK0s ; //ctauK0s 2.68 cm, mean life time of K0s is 8.95 x10^(-11)
+ Bool_t ctauCutLambda = ctauLambda < NCtau*fCutctauLambda; //ctauLambda 7.8 cm ,mean life is 2.6 x10 ^(-10) ***** 3xctau is the accepted limit
+ Bool_t ctauCutAntiLambda= ctauAntiLambda < NCtau*fCutctauAntiLambda;
+ Bool_t RapCutK0s = v0->RapK0Short() < fRapCutK0s;
+ Bool_t RapCutLambda = v0->RapLambda() < fRapCutLambda;
+ Bool_t RapCutAntiLambda = v0->Y(-3122) < fRapCutLambda;
- //vertex selection(is it fine for PP?)
- if ( fVertextype==1){//for pPb basically if(!fAnalysisUtils->IsVertexSelected2013pA(aod)) return;
- 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(fVertextype==2) {//for pp and pb-pb case , taken from Jan's code
- 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;
+ Bool_t CPACut= V0cosPointAngle > fMinCPA; //cosine of pointing angle with v0 should be greater than 0.998
- // 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;
+ //Now we put a loose mass cut which will be tightened later
+ Bool_t MassCutLooseK0s=(TMath::Abs(InvMassK0s - 0.497614) < 0.1);
+ Bool_t MassCutLooseLambda=(TMath::Abs(InvMassLambda - 1.115683) < 0.1); // cut is same for Anti-Lambda
+ Bool_t MassCutLooseAntiLambda=(TMath::Abs(InvMassAntiLambda - 1.115683) < 0.1); // cut is same for Anti-Lambda
- }
- else if(fVertextype==0){//default case
- trkVtx = aod->GetPrimaryVertex();
- if (!trkVtx || trkVtx->GetNContributors()<=0) return;//proper number of contributors
- zvtx = trkVtx->GetZ();
- Double32_t fCov[6];
- trkVtx->GetCovarianceMatrix(fCov);
- if(fCov[5] == 0) return;//proper vertex resolution
- }
- else {
- AliInfo("Wrong Vertextype set for Primary-vertex Selection: event REJECTED ...");
- return;//as there is no proper sample type
- }
+ //Special Cut for Kshort arementeros podalanski plot
+ Bool_t ArmenterosCut =kFALSE;
+ if(ctauCutK0s && RapCutK0s && CPACut && MassCutLooseK0s)
+ {
- fHistQA[0]->Fill((trkVtx->GetX()));fHistQA[1]->Fill((trkVtx->GetY()));fHistQA[2]->Fill((trkVtx->GetZ())); //for trkVtx only before vertex cut |zvtx|<10 cm
+ Float_t lAlphaV0 = v0->AlphaV0();
+ Float_t lPtArmV0 = v0->PtArmV0();
-//count events having a proper vertex
- fEventCounter->Fill(5);
+ ArmenterosCut = lPtArmV0 > TMath::Abs(0.2*lAlphaV0);
+
+ }
- if (TMath::Abs(zvtx) > fzvtxcut) return;
+ Bool_t IskShortOk=(ctauCutK0s && RapCutK0s && CPACut && MassCutLooseK0s && ArmenterosCut);
-//count events after vertex cut
- fEventCounter->Fill(7);
+ Bool_t IsLambdaOk=(ctauCutLambda && RapCutLambda && CPACut && MassCutLooseLambda);
+ Bool_t IsAntiLambdaOk=(ctauCutAntiLambda && RapCutAntiLambda && CPACut && MassCutLooseAntiLambda);
- //if(!fAnalysisUtils->IsVertexSelected2013pA(aod)) return;
+//Difference on Lambda and Anti-Lambda can be made through daughter PID
+
+
+ Int_t particletype=999;
+
+ if( IskShortOk || cutK0sPID )
+ {
+ fHistFinalPtCentInvK0s->Fill(InvMassK0s,v0Pt,Centrality);
+ particletype=SpKs0;
+
+ Short_t chargeval=0;
+ Float_t effmatrix=1.0;
+ LRCParticlePID* copy1 = new LRCParticlePID(particletype,InvMassK0s,chargeval,v0Pt,v0Eta, v0Phi,effmatrix,ptrack->GetTPCSharedMapPtr(),ntrack->GetTPCSharedMapPtr());
+ copy1->SetUniqueID(eventno * 200000 + (Int_t)iV0);
+ selectedV0s->Add(copy1);
- fHistQA[3]->Fill((trkVtx->GetX()));fHistQA[4]->Fill((trkVtx->GetY()));fHistQA[5]->Fill((trkVtx->GetZ()));//after vertex cut,for trkVtx only
+ }
+
+
+ if(IsLambdaOk || cutLambdaPID)
+ {
+ fHistFinalPtCentInvLambda->Fill(InvMassLambda,v0Pt,Centrality);
+//Add in the LRCParticle and give Lambda a tag 5
+ particletype=SpLam;
+
+ Short_t chargeval=0;
+ Float_t effmatrix=1.0;
+ LRCParticlePID* copy1 = new LRCParticlePID(particletype,InvMassLambda,chargeval,v0Pt,v0Eta, v0Phi,effmatrix,ptrack->GetTPCSharedMapPtr(),ntrack->GetTPCSharedMapPtr());
+ copy1->SetUniqueID(eventno * 200000 + (Int_t)iV0);
+ selectedV0s->Add(copy1);
+ }
+
+ if(IsAntiLambdaOk || cutAntiLambdaPID)
+ {
+ fHistFinalPtCentInvLambda->Fill(InvMassAntiLambda,v0Pt,Centrality);
+//Add in the LRCParticle and give Lambda a tag 6
+ particletype=SpALam;
+ Short_t chargeval=0;
+ Float_t effmatrix=1.0;
+ LRCParticlePID* copy1 = new LRCParticlePID(particletype,InvMassAntiLambda,chargeval,v0Pt,v0Eta, v0Phi,effmatrix,ptrack->GetTPCSharedMapPtr(),ntrack->GetTPCSharedMapPtr());
+ copy1->SetUniqueID(eventno * 200000 + (Int_t)iV0);
+ selectedV0s->Add(copy1);
+ }
+
+
+ }//v0 loop
+ return selectedV0s;
+}
+
+//___________________________________________________________________
+ Bool_t AliTwoParticlePIDCorr :: CheckStatusv0Daughter(AliAODTrack *t1 ,AliAODTrack *t2)
+ {
+ if (!t1->IsOn(AliAODTrack::kTPCrefit) || !t2->IsOn(AliAODTrack::kTPCrefit)) return kFALSE;
+ // Float_t nCrossedRowsTPC = t->GetTPCClusterInfo(2,1);
+if(t1->GetTPCClusterInfo(2,1)<fDaugNClsTPC || t2->GetTPCClusterInfo(2,1)<fDaugNClsTPC) return kFALSE ;
- if(!aod) return; //for safety
+// ---------------- Fraction of TPC Shared Cluster
+ Float_t fracPosDaugTPCSharedMap = GetFractionTPCSharedCls(t1);
+ Float_t fracNegDaugTPCSharedMap = GetFractionTPCSharedCls(t2);
- Double_t frefMult=0;
+ if( (fracPosDaugTPCSharedMap > fFracTPCcls) || (fracNegDaugTPCSharedMap > fFracTPCcls) )
+ return kFALSE;
+
+ return kTRUE;
-//reference multiplicity for pp 7 TeV
- if ((fMultiplicityEstimator == "TRACKS_MANUAL") || (fMultiplicityEstimator == "V0M_MANUAL")|| (fMultiplicityEstimator == "V0A_MANUAL")||(fMultiplicityEstimator == "V0C_MANUAL")) {cent_v0=GetReferenceMultiplicityVZEROFromAOD(aod,bSign1);}
- else {frefMult=GetReferenceMultiplicityVZEROFromAOD(aod,bSign1);}
+ }
+//___________________________________________________________________________________________
+
+ Float_t AliTwoParticlePIDCorr :: GetFractionTPCSharedCls( AliAODTrack *track)
+{
+ // Rejects tracks with shared clusters after filling a control histogram
+ // This overload is used for primaries
+ // Get the shared maps
+ const TBits sharedMap = track->GetTPCSharedMap();
+ return 1.*sharedMap.CountBits()/track->GetTPCNclsF();
-// centrality weighting (optional for 2011 if central and semicentral triggers are used)
- if (fCentralityWeights && !AcceptEventCentralityWeight(cent_v0))
+}
+//______________________________________________________________________
+ Bool_t AliTwoParticlePIDCorr :: CheckStatusv0(AliAODv0 *v1)
{
- AliInfo(Form("Rejecting event because of centrality weighting: %f", cent_v0));
- return;
+
+ // Offline reconstructed V0 only
+ if (v1->GetOnFlyStatus()) return kFALSE;
+
+ AliAODTrack *ptrack=(AliAODTrack *)v1->GetDaughter(0);
+ AliAODTrack *ntrack=(AliAODTrack *)v1->GetDaughter(1);
+
+ if(!ptrack || !ntrack) return kFALSE;
+
+ if(ptrack->Charge()==-1 || ntrack->Charge()==1) return kFALSE; //remove wrongly identified charge pairs
+
+ if(ptrack->Charge()==0 || ntrack->Charge()==0) return kFALSE; //remove uncharged pairs
+
+ if(ptrack->Charge() == ntrack->Charge()) return kFALSE; //remove like sign pairs
+
+ if(!CheckStatusv0Daughter(ptrack,ntrack)) return kFALSE;//daughters need to pass some basic cuts
+
+ if(TMath::Abs(ptrack->Eta()) > fmaxeta || TMath::Abs(ntrack->Eta()) > fmaxeta) return kFALSE; // remove daughters beyond eta bound |0.8|
+
+ if(ptrack->Pt() < fMinPtDaughter || ntrack->Pt() < fMinPtDaughter) return kFALSE; // remove daughter tracks below minmum p |1.0 GeV/c|
+
+ if(ptrack->Pt() > fMaxPtDaughter || ntrack->Pt() > fMaxPtDaughter) return kFALSE; // remove daughter tracks above maximum p ** to make it compatiable with AliHelperPID**|4.0 GeV/C|
+
+ // Daughters: Impact parameter of daughter to prim vtx
+ Float_t xy = v1->DcaPosToPrimVertex();
+ if (TMath::Abs(xy)<fDCAToPrimVtx) return kFALSE; //0.1 cm
+ xy = v1->DcaNegToPrimVertex();
+ if (TMath::Abs(xy)<fDCAToPrimVtx) return kFALSE; //0.1 cm
+
+ // Daughters: DCA
+ Float_t dca = v1->DcaV0Daughters();
+ if (dca>fMaxDCADaughter) return kFALSE; //1.0 cm
+
+ // V0: Cosine of the pointing angle
+ Float_t cpa=v1->CosPointingAngle(trkVtx); //0.997
+ if (cpa<fMinCPA) return kFALSE;
+
+ // V0: Fiducial volume
+ Double_t xyz[3]; v1->GetSecondaryVtx(xyz);
+ Float_t r2=xyz[0]*xyz[0] + xyz[1]*xyz[1];
+ if (r2<5.*5.) return kFALSE;
+ if (r2>lMax*lMax) return kFALSE; //lmax=100 cm
+
+ return kTRUE;
+
+
}
+//__________________________________________________________________________________________
+Bool_t AliTwoParticlePIDCorr::IsTrackFromV0(AliAODEvent* fAOD,AliAODTrack* track)
+{
+//to check whether a daughter being taken as associated
+ Int_t assoID = track->GetID();
+
+ for(int i=0; i<fAOD->GetNumberOfV0s(); i++){ // loop over V0s
+ AliAODv0* aodV0 = fAOD->GetV0(i);
+
+ AliAODTrack *trackPos=(AliAODTrack *)(aodV0->GetDaughter(0));
+ AliAODTrack *trackNeg=(AliAODTrack *)(aodV0->GetDaughter(1));
+
+
+ Int_t negID = trackNeg->GetID();
+ Int_t posID = trackPos->GetID();
+
+ if ((TMath::Abs(negID)+1)==(TMath::Abs(assoID))){ return kTRUE;}
+ if ((TMath::Abs(posID)+1)==(TMath::Abs(assoID))){ return kTRUE;}
+ //----------------------------------
+ }
+ return kFALSE;
+}
-//count events after rejection due to centrality weighting
- fEventCounter->Fill(9);
-*/
+
+
+//----------------------------------------------------------
+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