#include "TGraphErrors.h"\r
#include "TH1F.h"\r
#include "TH2F.h"\r
+#include "TH3F.h" \r
#include "TH2D.h" \r
#include "TH3D.h"\r
#include "TArrayF.h"\r
#include "AliAODEvent.h"\r
#include "AliAODTrack.h"\r
#include "AliAODInputHandler.h"\r
+#include "AliAODMCParticle.h" \r
+#include "AliCollisionGeometry.h"\r
#include "AliGenEventHeader.h"\r
-#include "AliGenHijingEventHeader.h"\r
#include "AliMCEventHandler.h"\r
#include "AliMCEvent.h"\r
#include "AliStack.h"\r
#include "AliESDtrackCuts.h"\r
#include "AliEventplane.h"\r
-#include "AliTHn.h" \r
+#include "AliTHn.h" \r
+#include "AliLog.h"\r
+#include "AliAnalysisUtils.h"\r
+\r
+#include "AliEventPoolManager.h" \r
\r
#include "AliPID.h" \r
#include "AliPIDResponse.h" \r
\r
#include "AliAnalysisTaskBFPsi.h"\r
#include "AliBalancePsi.h"\r
+#include "AliAnalysisTaskTriggeredBF.h"\r
\r
\r
// Analysis task for the BF vs Psi code\r
\r
//________________________________________________________________________\r
AliAnalysisTaskBFPsi::AliAnalysisTaskBFPsi(const char *name) \r
-: AliAnalysisTaskSE(name), \r
+: AliAnalysisTaskSE(name),\r
+ fDebugLevel(kFALSE),\r
+ fArrayMC(0), //+++++++++++++\r
fBalance(0),\r
fRunShuffling(kFALSE),\r
fShuffledBalance(0),\r
+ fRunMixing(kFALSE),\r
+ fRunMixingEventPlane(kFALSE),\r
+ fMixingTracks(50000),\r
+ fMixedBalance(0),\r
+ fPoolMgr(0),\r
fList(0),\r
fListBF(0),\r
fListBFS(0),\r
+ fListBFM(0),\r
fHistListPIDQA(0),\r
fHistEventStats(0),\r
fHistCentStats(0),\r
+ fHistCentStatsUsed(0),\r
fHistTriggerStats(0),\r
fHistTrackStats(0),\r
fHistVx(0),\r
fHistVy(0),\r
fHistVz(0),\r
+ fHistTPCvsVZEROMultiplicity(0),\r
+ fHistVZEROSignal(0),\r
fHistEventPlane(0),\r
fHistClus(0),\r
fHistDCA(0),\r
fHistEta(0),\r
fHistRapidity(0),\r
fHistPhi(0),\r
+ fHistEtaPhiPos(0), \r
+ fHistEtaPhiNeg(0), \r
fHistPhiBefore(0),\r
fHistPhiAfter(0),\r
fHistPhiPos(0),\r
fHistProbTPCTOFvsPtbeforePID(NULL),\r
fHistNSigmaTPCvsPtbeforePID(NULL), \r
fHistNSigmaTOFvsPtbeforePID(NULL), \r
+ fHistBetaVsdEdXbeforePID(NULL), //+++++++ \r
+ fHistNSigmaTPCTOFvsPtbeforePID(NULL), //++++++\r
fHistdEdxVsPTPCafterPID(NULL),\r
fHistBetavsPTOFafterPID(NULL), \r
fHistProbTPCvsPtafterPID(NULL), \r
fHistProbTPCTOFvsPtafterPID(NULL),\r
fHistNSigmaTPCvsPtafterPID(NULL), \r
fHistNSigmaTOFvsPtafterPID(NULL), \r
+ fHistBetaVsdEdXafterPID(NULL), //+++++++ \r
+ fHistNSigmaTPCTOFvsPtafterPID(NULL), //+++++++\r
+ fHistdEdxVsPTPCbeforePIDelectron(NULL), //+++++++\r
+ fHistNSigmaTPCvsPtbeforePIDelectron(NULL), //+++++++\r
+ fHistdEdxVsPTPCafterPIDelectron(NULL), //+++++++\r
+ fHistNSigmaTPCvsPtafterPIDelectron(NULL), //+++++++\r
+ fCentralityArrayBinsForCorrections(kCENTRALITY),\r
fPIDResponse(0x0),\r
fPIDCombined(0x0),\r
fParticleOfInterest(kPion),\r
fUsePIDPropabilities(kFALSE), \r
fPIDNSigma(3.),\r
fMinAcceptedPIDProbability(0.8),\r
+ fElectronRejection(kFALSE),\r
+ fElectronOnlyRejection(kFALSE),\r
+ fElectronRejectionNSigma(-1.),\r
+ fElectronRejectionMinPt(0.),\r
+ fElectronRejectionMaxPt(1000.),\r
fESDtrackCuts(0),\r
fCentralityEstimator("V0M"),\r
fUseCentrality(kFALSE),\r
fCentralityPercentileMax(5.),\r
fImpactParameterMin(0.),\r
fImpactParameterMax(20.),\r
+ fMultiplicityEstimator("V0A"),\r
fUseMultiplicity(kFALSE),\r
fNumberOfAcceptedTracksMin(0),\r
fNumberOfAcceptedTracksMax(10000),\r
fHistNumberOfAcceptedTracks(0),\r
+ fHistMultiplicity(0),\r
fUseOfflineTrigger(kFALSE),\r
+ fCheckFirstEventInChunk(kFALSE),\r
+ fCheckPileUp(kFALSE),\r
+ fUseMCforKinematics(kFALSE),\r
fVxMax(0.3),\r
fVyMax(0.3),\r
fVzMax(10.),\r
- nAODtrackCutBit(128),\r
+ fnAODtrackCutBit(128),\r
fPtMin(0.3),\r
fPtMax(1.5),\r
+ fPtMinForCorrections(0.3),//=================================correction\r
+ fPtMaxForCorrections(1.5),//=================================correction\r
+ fPtBinForCorrections(36), //=================================correction\r
fEtaMin(-0.8),\r
fEtaMax(-0.8),\r
+ fEtaMinForCorrections(-0.8),//=================================correction\r
+ fEtaMaxForCorrections(0.8),//=================================correction\r
+ fEtaBinForCorrections(16), //=================================correction\r
+ fPhiMin(0.),\r
+ fPhiMax(360.),\r
+ fPhiMinForCorrections(0.),//=================================correction\r
+ fPhiMaxForCorrections(360.),//=================================correction\r
+ fPhiBinForCorrections(100), //=================================correction\r
fDCAxyCut(-1),\r
fDCAzCut(-1),\r
fTPCchi2Cut(-1),\r
fDifferentialV2(0),\r
fUseFlowAfterBurner(kFALSE),\r
fExcludeResonancesInMC(kFALSE),\r
+ fExcludeElectronsInMC(kFALSE),\r
fUseMCPdgCode(kFALSE),\r
- fPDGCodeToBeAnalyzed(-1) {\r
+ fPDGCodeToBeAnalyzed(-1),\r
+ fEventClass("EventPlane"), \r
+ fCustomBinning(""),\r
+ fHistVZEROAGainEqualizationMap(0),\r
+ fHistVZEROCGainEqualizationMap(0),\r
+ fHistVZEROChannelGainEqualizationMap(0) {\r
// Constructor\r
// Define input and output slots here\r
// Input slot #0 works with a TChain\r
+\r
+ //======================================================correction\r
+ for (Int_t i=0; i<kCENTRALITY; i++){\r
+ fHistCorrectionPlus[i] = NULL; \r
+ fHistCorrectionMinus[i] = NULL; \r
+ fCentralityArrayForCorrections[i] = -1.;\r
+ }\r
+ //=====================================================correction\r
+\r
DefineInput(0, TChain::Class());\r
// Output slot #0 writes into a TH1 container\r
DefineOutput(1, TList::Class());\r
DefineOutput(2, TList::Class());\r
DefineOutput(3, TList::Class());\r
DefineOutput(4, TList::Class());\r
+ DefineOutput(5, TList::Class());\r
}\r
\r
//________________________________________________________________________\r
// delete fHistPt;\r
// delete fHistEta;\r
// delete fHistPhi;\r
+ // delete fHistEtaPhiPos; \r
+ // delete fHistEtaPhiNeg;\r
// delete fHistV0M;\r
}\r
\r
void AliAnalysisTaskBFPsi::UserCreateOutputObjects() {\r
// Create histograms\r
// Called once\r
+\r
+ // global switch disabling the reference \r
+ // (to avoid "Replacing existing TH1" if several wagons are created in train)\r
+ Bool_t oldStatus = TH1::AddDirectoryStatus();\r
+ TH1::AddDirectory(kFALSE);\r
+\r
if(!fBalance) {\r
fBalance = new AliBalancePsi();\r
fBalance->SetAnalysisLevel("ESD");\r
+ fBalance->SetEventClass(fEventClass);\r
//fBalance->SetNumberOfBins(-1,16);\r
//fBalance->SetInterval(-1,-0.8,0.8,16,0.,1.6,15.);\r
}\r
if(!fShuffledBalance) {\r
fShuffledBalance = new AliBalancePsi();\r
fShuffledBalance->SetAnalysisLevel("ESD");\r
+ fShuffledBalance->SetEventClass(fEventClass);\r
//fShuffledBalance->SetNumberOfBins(-1,16);\r
//fShuffledBalance->SetInterval(-1,-0.8,0.8,16,0.,1.6,15.);\r
}\r
}\r
+ if(fRunMixing) {\r
+ if(!fMixedBalance) {\r
+ fMixedBalance = new AliBalancePsi();\r
+ fMixedBalance->SetAnalysisLevel("ESD");\r
+ fMixedBalance->SetEventClass(fEventClass);\r
+ }\r
+ }\r
\r
//QA list\r
fList = new TList();\r
fListBFS->SetOwner();\r
}\r
\r
+ if(fRunMixing) {\r
+ fListBFM = new TList();\r
+ fListBFM->SetName("listTriggeredBFMixed");\r
+ fListBFM->SetOwner();\r
+ }\r
+\r
//PID QA list\r
- if(fUsePID) {\r
+ if(fUsePID || fElectronRejection) {\r
fHistListPIDQA = new TList();\r
fHistListPIDQA->SetName("listQAPID");\r
fHistListPIDQA->SetOwner();\r
}\r
\r
//Event stats.\r
- TString gCutName[4] = {"Total","Offline trigger",\r
- "Vertex","Analyzed"};\r
+ TString gCutName[7] = {"Total","Offline trigger",\r
+ "Vertex","Analyzed","sel. Centrality","Not1stEvInChunk","No Pile-Up"};\r
fHistEventStats = new TH2F("fHistEventStats",\r
"Event statistics;;Centrality percentile;N_{events}",\r
- 4,0.5,4.5,220,-5,105);\r
- for(Int_t i = 1; i <= 4; i++)\r
+ 7,0.5,7.5,220,-5,105);\r
+ for(Int_t i = 1; i <= 7; i++)\r
fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data());\r
fList->Add(fHistEventStats);\r
\r
- TString gCentName[9] = {"V0M","FMD","TRK","TKL","CL0","CL1","V0MvsFMD","TKLvsV0M","ZEMvsZDC"};\r
+ TString gCentName[13] = {"V0M","V0A","V0C","FMD","TRK","TKL","CL0","CL1","ZNA","ZPA","V0MvsFMD","TKLvsV0M","ZEMvsZDC"};\r
fHistCentStats = new TH2F("fHistCentStats",\r
"Centrality statistics;;Cent percentile",\r
- 9,-0.5,8.5,220,-5,105);\r
- for(Int_t i = 1; i <= 9; i++)\r
+ 13,-0.5,12.5,220,-5,105);\r
+ for(Int_t i = 1; i <= 13; i++){\r
fHistCentStats->GetXaxis()->SetBinLabel(i,gCentName[i-1].Data());\r
+ //fHistCentStatsUsed->GetXaxis()->SetBinLabel(i,gCentName[i-1].Data());\r
+ }\r
fList->Add(fHistCentStats);\r
\r
- fHistTriggerStats = new TH1F("fHistTriggerStats","Trigger statistics;TriggerBit;N_{events}",130,0,130);\r
+ fHistCentStatsUsed = new TH2F("fHistCentStatsUsed","Centrality statistics;;Cent percentile", 1,-0.5,0.5,220,-5,105);\r
+ fHistCentStatsUsed->GetXaxis()->SetBinLabel(1,fCentralityEstimator.Data());\r
+ fList->Add(fHistCentStatsUsed);\r
+\r
+ fHistTriggerStats = new TH1F("fHistTriggerStats","Trigger statistics;TriggerBit;N_{events}",1025,0,1025);\r
fList->Add(fHistTriggerStats);\r
\r
- fHistTrackStats = new TH1F("fHistTrackStats","Event statistics;TrackFilterBit;N_{events}",130,0,130);\r
+ fHistTrackStats = new TH1F("fHistTrackStats","Event statistics;TrackFilterBit;N_{events}",16,0,16);\r
fList->Add(fHistTrackStats);\r
\r
fHistNumberOfAcceptedTracks = new TH2F("fHistNumberOfAcceptedTracks",";N_{acc.};Centrality percentile;Entries",4001,-0.5,4000.5,220,-5,105);\r
fList->Add(fHistNumberOfAcceptedTracks);\r
\r
+ fHistMultiplicity = new TH1F("fHistMultiplicity",";N_{ch.};Entries",30001,-0.5,30000.5);\r
+ fList->Add(fHistMultiplicity);\r
+\r
// Vertex distributions\r
fHistVx = new TH1F("fHistVx","Primary vertex distribution - x coordinate;V_{x} (cm);Entries",100,-0.5,0.5);\r
fList->Add(fHistVx);\r
fHistVz = new TH2F("fHistVz","Primary vertex distribution - z coordinate;V_{z} (cm);Centrality percentile;Entries",100,-20.,20.,220,-5,105);\r
fList->Add(fHistVz);\r
\r
+ //TPC vs VZERO multiplicity\r
+ fHistTPCvsVZEROMultiplicity = new TH2F("fHistTPCvsVZEROMultiplicity","VZERO vs TPC multiplicity",3001,-0.5,30000.5,4001,-0.5,4000.5);\r
+ if(fMultiplicityEstimator == "V0A") \r
+ fHistTPCvsVZEROMultiplicity->GetXaxis()->SetTitle("VZERO-A multiplicity (a.u.)");\r
+ else if(fMultiplicityEstimator == "V0C") \r
+ fHistTPCvsVZEROMultiplicity->GetXaxis()->SetTitle("VZERO-C multiplicity (a.u.)");\r
+ else \r
+ fHistTPCvsVZEROMultiplicity->GetXaxis()->SetTitle("VZERO multiplicity (a.u.)");\r
+ fList->Add(fHistTPCvsVZEROMultiplicity);\r
+\r
+ fHistVZEROSignal = new TH2F("fHistVZEROSignal","VZERO signal vs VZERO channel;VZERO channel; Signal (a.u.)",64,0.5,64.5,3001,-0.5,30000.5);\r
+ fList->Add(fHistVZEROSignal);\r
+\r
//Event plane\r
fHistEventPlane = new TH2F("fHistEventPlane",";#Psi_{2} [deg.];Centrality percentile;Counts",100,0,360.,220,-5,105);\r
fList->Add(fHistEventPlane);\r
fList->Add(fHistEta);\r
fHistRapidity = new TH2F("fHistRapidity","y distribution;y;Centrality percentile",200,-2,2,220,-5,105);\r
fList->Add(fHistRapidity);\r
- fHistPhi = new TH2F("fHistPhi","#phi distribution;#phi;Centrality percentile",200,-20,380,220,-5,105);\r
+ fHistPhi = new TH2F("fHistPhi","#phi distribution;#phi (rad);Centrality percentile",200,0.0,2.*TMath::Pi(),220,-5,105);\r
fList->Add(fHistPhi);\r
+ fHistEtaPhiPos = new TH3F("fHistEtaPhiPos","#eta-#phi distribution (+);#eta;#phi (rad);Centrality percentile",80,-2,2,72,0.0,2.*TMath::Pi(),220,-5,105); \r
+ fList->Add(fHistEtaPhiPos); \r
+ fHistEtaPhiNeg = new TH3F("fHistEtaPhiNeg","#eta-#phi distribution (-);#eta;#phi (rad);Centrality percentile",80,-2,2,72,0.0,2.*TMath::Pi(),220,-5,105); \r
+ fList->Add(fHistEtaPhiNeg);\r
fHistPhiBefore = new TH2F("fHistPhiBefore","#phi distribution;#phi;Centrality percentile",200,0.,2*TMath::Pi(),220,-5,105);\r
fList->Add(fHistPhiBefore);\r
fHistPhiAfter = new TH2F("fHistPhiAfter","#phi distribution;#phi;Centrality percentile",200,0.,2*TMath::Pi(),220,-5,105);\r
fList->Add(fHistPhiAfter);\r
- fHistPhiPos = new TH2F("fHistPhiPos","#phi distribution for positive particles;#phi;Centrality percentile",200,-20,380,220,-5,105);\r
+ fHistPhiPos = new TH2F("fHistPhiPos","#phi distribution for positive particles;#phi;Centrality percentile",200,0.,2*TMath::Pi(),220,-5,105);\r
fList->Add(fHistPhiPos);\r
- fHistPhiNeg = new TH2F("fHistPhiNeg","#phi distribution for negative particles;#phi;Centrality percentile",200,-20,380,220,-5,105);\r
+ fHistPhiNeg = new TH2F("fHistPhiNeg","#phi distribution for negative particles;#phi;Centrality percentile",200,0.,2.*TMath::Pi(),220,-5,105);\r
fList->Add(fHistPhiNeg);\r
fHistV0M = new TH2F("fHistV0M","V0 Multiplicity C vs. A",500, 0, 20000, 500, 0, 20000);\r
fList->Add(fHistV0M);\r
fList->Add(fHistRefTracks);\r
\r
// Balance function histograms\r
- // Initialize histograms if not done yet\r
+ // Initialize histograms if not done yet (including the custom binning)\r
if(!fBalance->GetHistNp()){\r
- AliWarning("Histograms not yet initialized! --> Will be done now");\r
- AliWarning("--> Add 'gBalance->InitHistograms()' in your configBalanceFunction");\r
+ AliInfo("Histograms not yet initialized! --> Will be done now");\r
+ fBalance->SetCustomBinning(fCustomBinning);\r
fBalance->InitHistograms();\r
}\r
\r
if(fRunShuffling) {\r
if(!fShuffledBalance->GetHistNp()) {\r
- AliWarning("Histograms (shuffling) not yet initialized! --> Will be done now");\r
- AliWarning("--> Add 'gBalance->InitHistograms()' in your configBalanceFunction");\r
+ AliInfo("Histograms (shuffling) not yet initialized! --> Will be done now");\r
+ fShuffledBalance->SetCustomBinning(fCustomBinning);\r
fShuffledBalance->InitHistograms();\r
}\r
}\r
\r
+ if(fRunMixing) {\r
+ if(!fMixedBalance->GetHistNp()) {\r
+ AliInfo("Histograms (mixing) not yet initialized! --> Will be done now");\r
+ fMixedBalance->SetCustomBinning(fCustomBinning);\r
+ fMixedBalance->InitHistograms();\r
+ }\r
+ }\r
+\r
+ // QA histograms for different cuts\r
+ fList->Add(fBalance->GetQAHistHBTbefore());\r
+ fList->Add(fBalance->GetQAHistHBTafter());\r
+ fList->Add(fBalance->GetQAHistConversionbefore());\r
+ fList->Add(fBalance->GetQAHistConversionafter());\r
+ fList->Add(fBalance->GetQAHistPsiMinusPhi());\r
+ fList->Add(fBalance->GetQAHistResonancesBefore());\r
+ fList->Add(fBalance->GetQAHistResonancesRho());\r
+ fList->Add(fBalance->GetQAHistResonancesK0());\r
+ fList->Add(fBalance->GetQAHistResonancesLambda());\r
+ fList->Add(fBalance->GetQAHistQbefore());\r
+ fList->Add(fBalance->GetQAHistQafter());\r
+\r
//for(Int_t a = 0; a < ANALYSIS_TYPES; a++){\r
fListBF->Add(fBalance->GetHistNp());\r
fListBF->Add(fBalance->GetHistNn());\r
fListBFS->Add(fShuffledBalance->GetHistNpp());\r
fListBFS->Add(fShuffledBalance->GetHistNnp());\r
} \r
+\r
+ if(fRunMixing) {\r
+ fListBFM->Add(fMixedBalance->GetHistNp());\r
+ fListBFM->Add(fMixedBalance->GetHistNn());\r
+ fListBFM->Add(fMixedBalance->GetHistNpn());\r
+ fListBFM->Add(fMixedBalance->GetHistNnn());\r
+ fListBFM->Add(fMixedBalance->GetHistNpp());\r
+ fListBFM->Add(fMixedBalance->GetHistNnp());\r
+ }\r
//}\r
\r
+\r
+ // Event Mixing\r
+ if(fRunMixing){\r
+ Int_t trackDepth = fMixingTracks; \r
+ Int_t poolsize = 1000; // Maximum number of events, ignored in the present implemented of AliEventPoolManager\r
+ \r
+ // centrality bins\r
+ Double_t* centbins;\r
+ Int_t nCentralityBins;\r
+ if(fBalance->IsUseVertexBinning()){\r
+ centbins = fBalance->GetBinning(fBalance->GetBinningString(), "centralityVertex", nCentralityBins);\r
+ }\r
+ else{\r
+ centbins = fBalance->GetBinning(fBalance->GetBinningString(), "centrality", nCentralityBins);\r
+ }\r
+ \r
+ // multiplicity bins\r
+ Double_t* multbins;\r
+ Int_t nMultiplicityBins;\r
+ multbins = fBalance->GetBinning(fBalance->GetBinningString(), "multiplicity", nMultiplicityBins);\r
+ \r
+ // Zvtx bins\r
+ Double_t* vtxbins; \r
+ Int_t nVertexBins;\r
+ if(fBalance->IsUseVertexBinning()){\r
+ vtxbins = fBalance->GetBinning(fBalance->GetBinningString(), "vertexVertex", nVertexBins);\r
+ }\r
+ else{\r
+ vtxbins = fBalance->GetBinning(fBalance->GetBinningString(), "vertex", nVertexBins);\r
+ }\r
+\r
+ // Event plane angle (Psi) bins\r
+ Double_t* psibins;\r
+ Int_t nPsiBins; \r
+ psibins = fBalance->GetBinning(fBalance->GetBinningString(), "eventPlane", nPsiBins);\r
+\r
+ \r
+ // run the event mixing also in bins of event plane (statistics!)\r
+ if(fRunMixingEventPlane){\r
+ if(fEventClass=="Multiplicity"){\r
+ fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nMultiplicityBins, multbins, nVertexBins, vtxbins, nPsiBins, psibins);\r
+ }\r
+ else{\r
+ fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centbins, nVertexBins, vtxbins, nPsiBins, psibins);\r
+ }\r
+ }\r
+ else{\r
+ if(fEventClass=="Multiplicity"){\r
+ fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nMultiplicityBins, multbins, nVertexBins, vtxbins);\r
+ }\r
+ else{\r
+ fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centbins, nVertexBins, vtxbins);\r
+ }\r
+ }\r
+ }\r
+\r
if(fESDtrackCuts) fList->Add(fESDtrackCuts);\r
\r
//====================PID========================//\r
fPIDCombined->SetDefaultTPCPriors();\r
\r
fHistdEdxVsPTPCbeforePID = new TH2D ("dEdxVsPTPCbefore","dEdxVsPTPCbefore", 1000, -10.0, 10.0, 1000, 0, 1000); \r
- fHistListPIDQA->Add(fHistdEdxVsPTPCbeforePID); //addition \r
+ fHistListPIDQA->Add(fHistdEdxVsPTPCbeforePID);\r
\r
fHistBetavsPTOFbeforePID = new TH2D ("BetavsPTOFbefore","BetavsPTOFbefore", 1000, -10.0, 10., 1000, 0, 1.2); \r
- fHistListPIDQA->Add(fHistBetavsPTOFbeforePID); //addition\r
+ fHistListPIDQA->Add(fHistBetavsPTOFbeforePID); \r
\r
fHistProbTPCvsPtbeforePID = new TH2D ("ProbTPCvsPtbefore","ProbTPCvsPtbefore", 1000, -10.0,10.0, 1000, 0, 2.0); \r
- fHistListPIDQA->Add(fHistProbTPCvsPtbeforePID); //addition \r
+ fHistListPIDQA->Add(fHistProbTPCvsPtbeforePID); \r
\r
fHistProbTOFvsPtbeforePID = new TH2D ("ProbTOFvsPtbefore","ProbTOFvsPtbefore", 1000, -50, 50, 1000, 0, 2.0); \r
- fHistListPIDQA->Add(fHistProbTOFvsPtbeforePID); //addition \r
+ fHistListPIDQA->Add(fHistProbTOFvsPtbeforePID);\r
\r
fHistProbTPCTOFvsPtbeforePID =new TH2D ("ProbTPCTOFvsPtbefore","ProbTPCTOFvsPtbefore", 1000, -50, 50, 1000, 0, 2.0); \r
- fHistListPIDQA->Add(fHistProbTPCTOFvsPtbeforePID); //addition \r
+ fHistListPIDQA->Add(fHistProbTPCTOFvsPtbeforePID);\r
\r
fHistNSigmaTPCvsPtbeforePID = new TH2D ("NSigmaTPCvsPtbefore","NSigmaTPCvsPtbefore", 1000, -10, 10, 1000, 0, 500); \r
- fHistListPIDQA->Add(fHistNSigmaTPCvsPtbeforePID); //addition \r
+ fHistListPIDQA->Add(fHistNSigmaTPCvsPtbeforePID);\r
\r
fHistNSigmaTOFvsPtbeforePID = new TH2D ("NSigmaTOFvsPtbefore","NSigmaTOFvsPtbefore", 1000, -10, 10, 1000, 0, 500); \r
- fHistListPIDQA->Add(fHistNSigmaTOFvsPtbeforePID); //addition \r
+ fHistListPIDQA->Add(fHistNSigmaTOFvsPtbeforePID); \r
+\r
+ fHistBetaVsdEdXbeforePID = new TH2D ("BetaVsdEdXbefore","BetaVsdEdXbefore", 1000, 0., 1000, 1000, 0, 1.2); \r
+ fHistListPIDQA->Add(fHistBetaVsdEdXbeforePID); //++++++++\r
+ \r
+ fHistNSigmaTPCTOFvsPtbeforePID = new TH2D ("NSigmaTPCTOFvsPtbefore","NSigmaTPCTOFvsPtbefore", 1000, -10., 10., 1000, 0, 500); \r
+ fHistListPIDQA->Add(fHistNSigmaTPCTOFvsPtbeforePID); //++++++++\r
\r
fHistdEdxVsPTPCafterPID = new TH2D ("dEdxVsPTPCafter","dEdxVsPTPCafter", 1000, -10, 10, 1000, 0, 1000); \r
- fHistListPIDQA->Add(fHistdEdxVsPTPCafterPID); //addition \r
+ fHistListPIDQA->Add(fHistdEdxVsPTPCafterPID);\r
\r
fHistBetavsPTOFafterPID = new TH2D ("BetavsPTOFafter","BetavsPTOFafter", 1000, -10, 10, 1000, 0, 1.2); \r
- fHistListPIDQA->Add(fHistBetavsPTOFafterPID); //addition \r
+ fHistListPIDQA->Add(fHistBetavsPTOFafterPID); \r
\r
fHistProbTPCvsPtafterPID = new TH2D ("ProbTPCvsPtafter","ProbTPCvsPtafter", 1000, -10, 10, 1000, 0, 2); \r
- fHistListPIDQA->Add(fHistProbTPCvsPtafterPID); //addition \r
+ fHistListPIDQA->Add(fHistProbTPCvsPtafterPID);\r
\r
fHistProbTOFvsPtafterPID = new TH2D ("ProbTOFvsPtafter","ProbTOFvsPtafter", 1000, -10, 10, 1000, 0, 2); \r
- fHistListPIDQA->Add(fHistProbTOFvsPtafterPID); //addition \r
+ fHistListPIDQA->Add(fHistProbTOFvsPtafterPID); \r
\r
fHistProbTPCTOFvsPtafterPID =new TH2D ("ProbTPCTOFvsPtafter","ProbTPCTOFvsPtafter", 1000, -50, 50, 1000, 0, 2.0); \r
- fHistListPIDQA->Add(fHistProbTPCTOFvsPtafterPID); //addition \r
+ fHistListPIDQA->Add(fHistProbTPCTOFvsPtafterPID);\r
\r
fHistNSigmaTPCvsPtafterPID = new TH2D ("NSigmaTPCvsPtafter","NSigmaTPCvsPtafter", 1000, -10, 10, 1000, 0, 500); \r
- fHistListPIDQA->Add(fHistNSigmaTPCvsPtafterPID); //addition \r
+ fHistListPIDQA->Add(fHistNSigmaTPCvsPtafterPID);\r
\r
fHistNSigmaTOFvsPtafterPID = new TH2D ("NSigmaTOFvsPtafter","NSigmaTOFvsPtafter", 1000, -10, 10, 1000, 0, 500); \r
- fHistListPIDQA->Add(fHistNSigmaTOFvsPtafterPID); //addition \r
+ fHistListPIDQA->Add(fHistNSigmaTOFvsPtafterPID);\r
+\r
+ fHistBetaVsdEdXafterPID = new TH2D ("BetaVsdEdXafter","BetaVsdEdXafter", 1000, 0., 1000, 1000, 0, 1.2); \r
+ fHistListPIDQA->Add(fHistBetaVsdEdXafterPID); //++++++++\r
+\r
+ fHistNSigmaTPCTOFvsPtafterPID = new TH2D ("NSigmaTPCTOFvsPtafter","NSigmaTPCTOFvsPtafter", 1000, -10., 10., 1000, 0, 500); \r
+ fHistListPIDQA->Add(fHistNSigmaTPCTOFvsPtafterPID); //++++++++\r
+ }\r
+\r
+ // for electron rejection only TPC nsigma histograms\r
+ if(fElectronRejection) {\r
+ \r
+ fHistdEdxVsPTPCbeforePIDelectron = new TH2D ("dEdxVsPTPCbeforeelectron","dEdxVsPTPCbeforeelectron", 1000, -10.0, 10.0, 1000, 0, 1000); \r
+ fHistListPIDQA->Add(fHistdEdxVsPTPCbeforePIDelectron);\r
+ \r
+ fHistNSigmaTPCvsPtbeforePIDelectron = new TH2D ("NSigmaTPCvsPtbeforeelectron","NSigmaTPCvsPtbeforeelectron", 1000, -10, 10, 1000, 0, 500); \r
+ fHistListPIDQA->Add(fHistNSigmaTPCvsPtbeforePIDelectron);\r
+ \r
+ fHistdEdxVsPTPCafterPIDelectron = new TH2D ("dEdxVsPTPCafterelectron","dEdxVsPTPCafterelectron", 1000, -10, 10, 1000, 0, 1000); \r
+ fHistListPIDQA->Add(fHistdEdxVsPTPCafterPIDelectron);\r
+\r
+ fHistNSigmaTPCvsPtafterPIDelectron = new TH2D ("NSigmaTPCvsPtafterelectron","NSigmaTPCvsPtafterelectron", 1000, -10, 10, 1000, 0, 500); \r
+ fHistListPIDQA->Add(fHistNSigmaTPCvsPtafterPIDelectron); \r
}\r
//====================PID========================//\r
\r
PostData(1, fList);\r
PostData(2, fListBF);\r
if(fRunShuffling) PostData(3, fListBFS);\r
- if(fUsePID) PostData(4, fHistListPIDQA); //PID\r
+ if(fRunMixing) PostData(4, fListBFM);\r
+ if(fUsePID || fElectronRejection) PostData(5, fHistListPIDQA); //PID\r
+\r
+ AliInfo("Finished setting up the Output");\r
+\r
+ TH1::AddDirectory(oldStatus);\r
+}\r
+\r
+\r
+//________________________________________________________________________\r
+void AliAnalysisTaskBFPsi::SetInputCorrection(TString filename, \r
+ Int_t nCentralityBins, \r
+ Double_t *centralityArrayForCorrections) {\r
+ //Open files that will be used for correction\r
+ fCentralityArrayBinsForCorrections = nCentralityBins;\r
+ for (Int_t i=0; i<nCentralityBins; i++)\r
+ fCentralityArrayForCorrections[i] = centralityArrayForCorrections[i];\r
+\r
+ // No file specified -> run without corrections\r
+ if(!filename.Contains(".root")) {\r
+ AliInfo(Form("No correction file specified (= %s) --> run without corrections",filename.Data()));\r
+ return;\r
+ }\r
+\r
+ //Open the input file\r
+ TFile *f = TFile::Open(filename);\r
+ if(!f->IsOpen()) {\r
+ AliInfo(Form("File %s not found --> run without corrections",filename.Data()));\r
+ return;\r
+ }\r
+ \r
+ //TString listEffName = "";\r
+ for (Int_t iCent = 0; iCent < fCentralityArrayBinsForCorrections-1; iCent++) { \r
+ //Printf("iCent %d:",iCent); \r
+ TString histoName = "fHistCorrectionPlus";\r
+ histoName += Form("%d-%d",(Int_t)(fCentralityArrayForCorrections[iCent]),(Int_t)(fCentralityArrayForCorrections[iCent+1]));\r
+ fHistCorrectionPlus[iCent]= dynamic_cast<TH3F *>(f->Get(histoName.Data()));\r
+ if(!fHistCorrectionPlus[iCent]) {\r
+ AliError(Form("fHist %s not found!!!",histoName.Data()));\r
+ return;\r
+ }\r
+ \r
+ histoName = "fHistCorrectionMinus";\r
+ histoName += Form("%d-%d",(Int_t)(fCentralityArrayForCorrections[iCent]),(Int_t)(fCentralityArrayForCorrections[iCent+1]));\r
+ fHistCorrectionMinus[iCent] = dynamic_cast<TH3F *>(f->Get(histoName.Data())); \r
+ if(!fHistCorrectionMinus[iCent]) {\r
+ AliError(Form("fHist %s not found!!!",histoName.Data()));\r
+ return; \r
+ }\r
+ }//loop over centralities: ONLY the PbPb case is covered\r
+\r
+ if(fHistCorrectionPlus[0]){\r
+ fEtaMinForCorrections = fHistCorrectionPlus[0]->GetXaxis()->GetXmin();\r
+ fEtaMaxForCorrections = fHistCorrectionPlus[0]->GetXaxis()->GetXmax();\r
+ fEtaBinForCorrections = fHistCorrectionPlus[0]->GetNbinsX();\r
+ \r
+ fPtMinForCorrections = fHistCorrectionPlus[0]->GetYaxis()->GetXmin();\r
+ fPtMaxForCorrections = fHistCorrectionPlus[0]->GetYaxis()->GetXmax();\r
+ fPtBinForCorrections = fHistCorrectionPlus[0]->GetNbinsY();\r
+ \r
+ fPhiMinForCorrections = fHistCorrectionPlus[0]->GetZaxis()->GetXmin();\r
+ fPhiMaxForCorrections = fHistCorrectionPlus[0]->GetZaxis()->GetXmax();\r
+ fPhiBinForCorrections = fHistCorrectionPlus[0]->GetNbinsZ();\r
+ }\r
}\r
\r
//________________________________________________________________________\r
void AliAnalysisTaskBFPsi::UserExec(Option_t *) {\r
// Main loop\r
// Called for each event\r
- TString gAnalysisLevel = fBalance->GetAnalysisLevel();\r
-\r
- AliESDtrack *track_TPC = NULL;\r
\r
+ TString gAnalysisLevel = fBalance->GetAnalysisLevel();\r
Int_t gNumberOfAcceptedTracks = 0;\r
- Float_t fCentrality = 0.;\r
- Double_t gReactionPlane = 0.;\r
-\r
- // vector holding the charges/kinematics of all tracks (charge,y,eta,phi,p0,p1,p2,pt,E)\r
- vector<Double_t> *chargeVectorShuffle[9]; // this will be shuffled\r
- vector<Double_t> *chargeVector[9]; // original charge\r
- for(Int_t i = 0; i < 9; i++){\r
- chargeVectorShuffle[i] = new vector<Double_t>;\r
- chargeVector[i] = new vector<Double_t>;\r
+ Double_t lMultiplicityVar = -999.; //-1\r
+ Double_t gReactionPlane = -1.; \r
+ Float_t bSign = 0.;\r
+ \r
+ // get the event (for generator level: MCEvent())\r
+ AliVEvent* eventMain = NULL;\r
+ if(gAnalysisLevel == "MC") {\r
+ eventMain = dynamic_cast<AliVEvent*>(MCEvent()); \r
}\r
-\r
- Double_t v_charge;\r
- Double_t v_y;\r
- Double_t v_eta;\r
- Double_t v_phi;\r
- Double_t v_p[3];\r
- Double_t v_pt;\r
- Double_t v_E;\r
-\r
- if(fUsePID) {\r
+ else{\r
+ eventMain = dynamic_cast<AliVEvent*>(InputEvent()); \r
+ // for HBT like cuts need magnetic field sign\r
+ bSign = (eventMain->GetMagneticField() > 0) ? 1 : -1;\r
+ }\r
+ if(!eventMain) {\r
+ AliError("eventMain not available");\r
+ return;\r
+ }\r
+ \r
+ // PID Response task active?\r
+ if(fUsePID || fElectronRejection) {\r
fPIDResponse = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetPIDResponse();\r
if (!fPIDResponse) AliFatal("This Task needs the PID response attached to the inputHandler");\r
}\r
\r
- //ESD analysis\r
- if(gAnalysisLevel == "ESD") {\r
- AliESDEvent* gESD = dynamic_cast<AliESDEvent*>(InputEvent()); // from TaskSE\r
- if (!gESD) {\r
- Printf("ERROR: gESD not available");\r
+ // check event cuts and fill event histograms\r
+ if((lMultiplicityVar = IsEventAccepted(eventMain)) < 0){ \r
+ return;\r
+ }\r
+ // get the reaction plane\r
+ if(fEventClass != "Multiplicity") {\r
+ gReactionPlane = GetEventPlane(eventMain);\r
+ fHistEventPlane->Fill(gReactionPlane,lMultiplicityVar);\r
+ if(gReactionPlane < 0){\r
return;\r
}\r
+ }\r
+ \r
+ // get the accepted tracks in main event\r
+ TObjArray *tracksMain = GetAcceptedTracks(eventMain,lMultiplicityVar,gReactionPlane);\r
+ gNumberOfAcceptedTracks = tracksMain->GetEntriesFast();\r
\r
- // store offline trigger bits\r
- fHistTriggerStats->Fill(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected());\r
-\r
- // event selection done in AliAnalysisTaskSE::Exec() --> this is not used\r
- fHistEventStats->Fill(1,fCentrality); //all events\r
- Bool_t isSelected = kTRUE;\r
- if(fUseOfflineTrigger)\r
- isSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();\r
- if(isSelected) {\r
- fHistEventStats->Fill(2,fCentrality); //triggered events\r
+ //multiplicity cut (used in pp)\r
+ fHistNumberOfAcceptedTracks->Fill(gNumberOfAcceptedTracks,lMultiplicityVar);\r
\r
- if(fUseCentrality) {\r
- //Centrality stuff\r
- AliCentrality *centrality = gESD->GetCentrality();\r
+ // store charges of all accepted tracks,shuffle and reassign(two extra loops)\r
+ TObjArray* tracksShuffled = NULL;\r
+ if(fRunShuffling){\r
+ tracksShuffled = GetShuffledTracks(tracksMain,lMultiplicityVar);\r
+ }\r
+ \r
+ // Event mixing \r
+ if (fRunMixing)\r
+ {\r
+ // 1. First get an event pool corresponding in mult (cent) and\r
+ // zvertex to the current event. Once initialized, the pool\r
+ // should contain nMix (reduced) events. This routine does not\r
+ // pre-scan the chain. The first several events of every chain\r
+ // will be skipped until the needed pools are filled to the\r
+ // specified depth. If the pool categories are not too rare, this\r
+ // should not be a problem. If they are rare, you could lose`\r
+ // statistics.\r
+ \r
+ // 2. Collect the whole pool's content of tracks into one TObjArray\r
+ // (bgTracks), which is effectively a single background super-event.\r
+ \r
+ // 3. The reduced and bgTracks arrays must both be passed into\r
+ // FillCorrelations(). Also nMix should be passed in, so a weight\r
+ // of 1./nMix can be applied.\r
+ \r
+ AliEventPool* pool = fPoolMgr->GetEventPool(lMultiplicityVar, eventMain->GetPrimaryVertex()->GetZ(),gReactionPlane);\r
+ \r
+ if (!pool){\r
+ AliFatal(Form("No pool found for centrality = %f, zVtx = %f, psi = %f", lMultiplicityVar, eventMain->GetPrimaryVertex()->GetZ(),gReactionPlane));\r
+ }\r
+ else{\r
\r
- fCentrality = centrality->GetCentralityPercentile(fCentralityEstimator.Data());\r
+ //pool->SetDebug(1);\r
+\r
+ if (pool->IsReady() || pool->NTracksInPool() > fMixingTracks / 10 || pool->GetCurrentNEvents() >= 5){ \r
+ \r
+ \r
+ Int_t nMix = pool->GetCurrentNEvents();\r
+ //cout << "nMix = " << nMix << " tracks in pool = " << pool->NTracksInPool() << endl;\r
+ \r
+ //((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(2);\r
+ //((TH2F*) fListOfHistos->FindObject("mixedDist"))->Fill(centrality, pool->NTracksInPool());\r
+ //if (pool->IsReady())\r
+ //((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(3);\r
+ \r
+ // Fill mixed-event histos here \r
+ for (Int_t jMix=0; jMix<nMix; jMix++) \r
+ {\r
+ TObjArray* tracksMixed = pool->GetEvent(jMix);\r
+ fMixedBalance->CalculateBalance(gReactionPlane,tracksMain,tracksMixed,bSign,lMultiplicityVar,eventMain->GetPrimaryVertex()->GetZ());\r
+ }\r
+ }\r
\r
- // take only events inside centrality class\r
- if(!centrality->IsEventInCentralityClass(fCentralityPercentileMin,\r
- fCentralityPercentileMax,\r
- fCentralityEstimator.Data()))\r
- return;\r
+ // Update the Event pool\r
+ pool->UpdatePool(tracksMain);\r
+ //pool->PrintInfo();\r
\r
- // centrality QA (V0M)\r
- fHistV0M->Fill(gESD->GetVZEROData()->GetMTotV0A(), gESD->GetVZEROData()->GetMTotV0C());\r
+ }//pool NULL check \r
+ }//run mixing\r
+ \r
+ // calculate balance function\r
+ fBalance->CalculateBalance(gReactionPlane,tracksMain,NULL,bSign,lMultiplicityVar,eventMain->GetPrimaryVertex()->GetZ());\r
+ \r
+ // calculate shuffled balance function\r
+ if(fRunShuffling && tracksShuffled != NULL) {\r
+ fShuffledBalance->CalculateBalance(gReactionPlane,tracksShuffled,NULL,bSign,lMultiplicityVar,eventMain->GetPrimaryVertex()->GetZ());\r
+ }\r
+} \r
+\r
+//________________________________________________________________________\r
+Double_t AliAnalysisTaskBFPsi::IsEventAccepted(AliVEvent *event){\r
+ // Checks the Event cuts\r
+ // Fills Event statistics histograms\r
+ \r
+ Bool_t isSelectedMain = kTRUE;\r
+ Float_t gRefMultiplicity = -1.;\r
+ TString gAnalysisLevel = fBalance->GetAnalysisLevel();\r
+\r
+ AliMCEvent *mcevent = dynamic_cast<AliMCEvent*>(event);\r
+\r
+ fHistEventStats->Fill(1,gRefMultiplicity); //all events\r
+\r
+ // check first event in chunk (is not needed for new reconstructions)\r
+ if(fCheckFirstEventInChunk){\r
+ AliAnalysisUtils ut;\r
+ if(ut.IsFirstEventInChunk(event)) \r
+ return -1.;\r
+ fHistEventStats->Fill(6,gRefMultiplicity); \r
+ }\r
+ // check for pile-up event\r
+ if(fCheckPileUp){\r
+ AliAnalysisUtils ut;\r
+ ut.SetUseMVPlpSelection(kTRUE);\r
+ ut.SetUseOutOfBunchPileUp(kTRUE);\r
+ if(ut.IsPileUpEvent(event))\r
+ return -1.;\r
+ fHistEventStats->Fill(7,gRefMultiplicity); \r
+ }\r
+\r
+ // Event trigger bits\r
+ fHistTriggerStats->Fill(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected());\r
+ if(fUseOfflineTrigger)\r
+ isSelectedMain = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();\r
+ \r
+ if(isSelectedMain) {\r
+ fHistEventStats->Fill(2,gRefMultiplicity); //triggered events\r
+ \r
+ // Event Vertex MC\r
+ if(gAnalysisLevel == "MC") {\r
+ if(!event) {\r
+ AliError("mcEvent not available");\r
+ return 0x0;\r
}\r
- \r
- const AliESDVertex *vertex = gESD->GetPrimaryVertex();\r
+ \r
+ if(mcevent){\r
+ AliGenEventHeader *header = dynamic_cast<AliGenEventHeader*>(mcevent->GenEventHeader());\r
+ if(header){ \r
+ TArrayF gVertexArray;\r
+ header->PrimaryVertex(gVertexArray);\r
+ //Printf("Vertex: %lf (x) - %lf (y) - %lf (z)",\r
+ //gVertexArray.At(0),\r
+ //gVertexArray.At(1),\r
+ //gVertexArray.At(2));\r
+ fHistEventStats->Fill(3,gRefMultiplicity); //events with a proper vertex\r
+ if(TMath::Abs(gVertexArray.At(0)) < fVxMax) {\r
+ if(TMath::Abs(gVertexArray.At(1)) < fVyMax) {\r
+ if(TMath::Abs(gVertexArray.At(2)) < fVzMax) {\r
+ fHistEventStats->Fill(4,gRefMultiplicity);//analyzed events\r
+\r
+ // get the reference multiplicty or centrality\r
+ gRefMultiplicity = GetRefMultiOrCentrality(event);\r
+\r
+ fHistVx->Fill(gVertexArray.At(0));\r
+ fHistVy->Fill(gVertexArray.At(1));\r
+ fHistVz->Fill(gVertexArray.At(2),gRefMultiplicity);\r
+ \r
+ // take only events inside centrality class\r
+ if(fUseCentrality) {\r
+ if((fImpactParameterMin < gRefMultiplicity) && (fImpactParameterMax > gRefMultiplicity)){\r
+ fHistEventStats->Fill(5,gRefMultiplicity); //events with correct centrality\r
+ return gRefMultiplicity; \r
+ }//centrality class\r
+ }\r
+ // take events only within the same multiplicity class\r
+ else if(fUseMultiplicity) {\r
+ if((gRefMultiplicity > fNumberOfAcceptedTracksMin) && (gRefMultiplicity < fNumberOfAcceptedTracksMax)) {\r
+ fHistEventStats->Fill(5,gRefMultiplicity); //events with correct multiplicity\r
+ return gRefMultiplicity;\r
+ }\r
+ }//multiplicity range\r
+ }//Vz cut\r
+ }//Vy cut\r
+ }//Vx cut\r
+ }//header \r
+ }//MC event object\r
+ }//MC\r
+ \r
+ // Event Vertex AOD, ESD, ESDMC\r
+ else{\r
+ const AliVVertex *vertex = event->GetPrimaryVertex();\r
+ \r
if(vertex) {\r
+ Double32_t fCov[6];\r
+ vertex->GetCovarianceMatrix(fCov);\r
if(vertex->GetNContributors() > 0) {\r
- if(vertex->GetZRes() != 0) {\r
- fHistEventStats->Fill(3,fCentrality); //events with a proper vertex\r
- if(TMath::Abs(vertex->GetXv()) < fVxMax) {\r
- if(TMath::Abs(vertex->GetYv()) < fVyMax) {\r
- if(TMath::Abs(vertex->GetZv()) < fVzMax) {\r
- fHistEventStats->Fill(4,fCentrality); //analayzed events\r
- fHistVx->Fill(vertex->GetXv());\r
- fHistVy->Fill(vertex->GetYv());\r
- fHistVz->Fill(vertex->GetZv(),fCentrality);\r
+ if(fCov[5] != 0) {\r
+ fHistEventStats->Fill(3,gRefMultiplicity); //proper vertex\r
+ if(TMath::Abs(vertex->GetX()) < fVxMax) {\r
+ if(TMath::Abs(vertex->GetY()) < fVyMax) {\r
+ if(TMath::Abs(vertex->GetZ()) < fVzMax) {\r
+ fHistEventStats->Fill(4,gRefMultiplicity);//analyzed events\r
+\r
+ // get the reference multiplicty or centrality\r
+ gRefMultiplicity = GetRefMultiOrCentrality(event);\r
\r
- //========Get the VZERO event plane========//\r
- Double_t gVZEROEventPlane = -10.0;\r
- Double_t qxTot = 0.0, qyTot = 0.0;\r
- AliEventplane *ep = gESD->GetEventplane();\r
- if(ep) \r
- gVZEROEventPlane = ep->CalculateVZEROEventPlane(gESD,10,2,qxTot,qyTot);\r
- if(gVZEROEventPlane < 0.) gVZEROEventPlane += TMath::Pi();\r
- gReactionPlane = gVZEROEventPlane*TMath::RadToDeg();\r
- fHistEventPlane->Fill(gReactionPlane,fCentrality);\r
- //========Get the VZERO event plane========//\r
-\r
- //Printf("There are %d tracks in this event", gESD->GetNumberOfTracks());\r
- for (Int_t iTracks = 0; iTracks < gESD->GetNumberOfTracks(); iTracks++) {\r
- AliESDtrack* track = dynamic_cast<AliESDtrack *>(gESD->GetTrack(iTracks));\r
- if (!track) {\r
- Printf("ERROR: Could not receive track %d", iTracks);\r
- continue;\r
- } \r
- \r
- // take only TPC only tracks\r
- track_TPC = new AliESDtrack();\r
- if(!track->FillTPCOnlyTrack(*track_TPC)) continue;\r
- \r
- //ESD track cuts\r
- if(fESDtrackCuts) \r
- if(!fESDtrackCuts->AcceptTrack(track_TPC)) continue;\r
- \r
- // fill QA histograms\r
- Float_t b[2];\r
- Float_t bCov[3];\r
- track_TPC->GetImpactParameters(b,bCov);\r
- if (bCov[0]<=0 || bCov[2]<=0) {\r
- AliDebug(1, "Estimated b resolution lower or equal zero!");\r
- bCov[0]=0; bCov[2]=0;\r
- }\r
- \r
- Int_t nClustersTPC = -1;\r
- nClustersTPC = track_TPC->GetTPCNclsIter1(); // TPC standalone\r
- //nClustersTPC = track->GetTPCclusters(0); // global track\r
- Float_t chi2PerClusterTPC = -1;\r
- if (nClustersTPC!=0) {\r
- chi2PerClusterTPC = track_TPC->GetTPCchi2Iter1()/Float_t(nClustersTPC); // TPC standalone\r
- //chi2PerClusterTPC = track->GetTPCchi2()/Float_t(nClustersTPC); // global track\r
- }\r
-\r
- //===========================PID===============================// \r
- if(fUsePID) {\r
- Double_t prob[AliPID::kSPECIES]={0.};\r
- Double_t probTPC[AliPID::kSPECIES]={0.};\r
- Double_t probTOF[AliPID::kSPECIES]={0.};\r
- Double_t probTPCTOF[AliPID::kSPECIES]={0.};\r
-\r
- Double_t nSigma = 0.;\r
- UInt_t detUsedTPC = 0;\r
- UInt_t detUsedTOF = 0;\r
- UInt_t detUsedTPCTOF = 0;\r
-\r
- //Decide what detector configuration we want to use\r
- switch(fPidDetectorConfig) {\r
- case kTPCpid:\r
- fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);\r
- nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)fParticleOfInterest));\r
- detUsedTPC = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTPC);\r
- for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)\r
- prob[iSpecies] = probTPC[iSpecies];\r
- break;\r
- case kTOFpid:\r
- fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF);\r
- nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)fParticleOfInterest));\r
- detUsedTOF = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTOF);\r
- for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)\r
- prob[iSpecies] = probTOF[iSpecies];\r
- break;\r
- case kTPCTOF:\r
- fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC);\r
- detUsedTPCTOF = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTPCTOF);\r
- for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)\r
- prob[iSpecies] = probTPCTOF[iSpecies];\r
- break;\r
- default:\r
- break;\r
- }//end switch: define detector mask\r
- \r
- //Filling the PID QA\r
- Double_t tofTime = -999., length = 999., tof = -999.;\r
- Double_t c = TMath::C()*1.E-9;// m/ns\r
- Double_t beta = -999.;\r
- Double_t nSigmaTOFForParticleOfInterest = -999.;\r
- if ( (track->IsOn(AliESDtrack::kTOFin)) &&\r
- (track->IsOn(AliESDtrack::kTIME)) ) { \r
- tofTime = track->GetTOFsignal();//in ps\r
- length = track->GetIntegratedLength();\r
- tof = tofTime*1E-3; // ns \r
- \r
- if (tof <= 0) {\r
- //Printf("WARNING: track with negative TOF time found! Skipping this track for PID checks\n");\r
- continue;\r
- }\r
- if (length <= 0){\r
- //printf("WARNING: track with negative length found!Skipping this track for PID checks\n");\r
- continue;\r
- }\r
- \r
- length = length*0.01; // in meters\r
- tof = tof*c;\r
- beta = length/tof;\r
- \r
- nSigmaTOFForParticleOfInterest = fPIDResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)fParticleOfInterest);\r
- fHistBetavsPTOFbeforePID ->Fill(track->P()*track->Charge(),beta);\r
- fHistProbTOFvsPtbeforePID ->Fill(track->Pt(),probTOF[fParticleOfInterest]);\r
- fHistNSigmaTOFvsPtbeforePID ->Fill(track->Pt(),nSigmaTOFForParticleOfInterest);\r
- }//TOF signal \r
- \r
- \r
- Double_t nSigmaTPCForParticleOfInterest = fPIDResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)fParticleOfInterest);\r
- fHistdEdxVsPTPCbeforePID -> Fill(track->P()*track->Charge(),track->GetTPCsignal());\r
- fHistProbTPCvsPtbeforePID -> Fill(track->Pt(),probTPC[fParticleOfInterest]); \r
- fHistNSigmaTPCvsPtbeforePID -> Fill(track->Pt(),nSigmaTPCForParticleOfInterest); \r
- fHistProbTPCTOFvsPtbeforePID -> Fill(track->Pt(),probTPCTOF[fParticleOfInterest]);\r
- //end of QA-before pid\r
- \r
- if ((detUsedTPC != 0)||(detUsedTOF != 0)||(detUsedTPCTOF != 0)) {\r
- //Make the decision based on the n-sigma\r
- if(fUsePIDnSigma) {\r
- if(nSigma > fPIDNSigma) continue;}\r
- \r
- //Make the decision based on the bayesian\r
- else if(fUsePIDPropabilities) {\r
- if(fParticleOfInterest != TMath::LocMax(AliPID::kSPECIES,prob)) continue;\r
- if (prob[fParticleOfInterest] < fMinAcceptedPIDProbability) continue; \r
- }\r
- \r
- //Fill QA after the PID\r
- fHistBetavsPTOFafterPID ->Fill(track->P()*track->Charge(),beta);\r
- fHistProbTOFvsPtafterPID ->Fill(track->Pt(),probTOF[fParticleOfInterest]);\r
- fHistNSigmaTOFvsPtafterPID ->Fill(track->Pt(),nSigmaTOFForParticleOfInterest);\r
- \r
- fHistdEdxVsPTPCafterPID -> Fill(track->P()*track->Charge(),track->GetTPCsignal());\r
- fHistProbTPCvsPtafterPID -> Fill(track->Pt(),probTPC[fParticleOfInterest]); \r
- fHistProbTPCTOFvsPtafterPID -> Fill(track->Pt(),probTPCTOF[fParticleOfInterest]);\r
- fHistNSigmaTPCvsPtafterPID -> Fill(track->Pt(),nSigmaTPCForParticleOfInterest); \r
- }\r
- \r
- PostData(4, fHistListPIDQA);\r
- }\r
- //===========================PID===============================//\r
- v_charge = track_TPC->Charge();\r
- v_y = track_TPC->Y();\r
- v_eta = track_TPC->Eta();\r
- v_phi = track_TPC->Phi() * TMath::RadToDeg();\r
- v_E = track_TPC->E();\r
- v_pt = track_TPC->Pt();\r
- track_TPC->PxPyPz(v_p);\r
- fHistClus->Fill(track_TPC->GetITSclusters(0),nClustersTPC);\r
- fHistDCA->Fill(b[1],b[0]);\r
- fHistChi2->Fill(chi2PerClusterTPC,fCentrality);\r
- fHistPt->Fill(v_pt,fCentrality);\r
- fHistEta->Fill(v_eta,fCentrality);\r
- fHistPhi->Fill(v_phi,fCentrality);\r
- fHistRapidity->Fill(v_y,fCentrality);\r
- if(v_charge > 0) fHistPhiPos->Fill(v_phi,fCentrality);\r
- else if(v_charge < 0) fHistPhiNeg->Fill(v_phi,fCentrality);\r
-\r
- // fill charge vector\r
- chargeVector[0]->push_back(v_charge);\r
- chargeVector[1]->push_back(v_y);\r
- chargeVector[2]->push_back(v_eta);\r
- chargeVector[3]->push_back(v_phi);\r
- chargeVector[4]->push_back(v_p[0]);\r
- chargeVector[5]->push_back(v_p[1]);\r
- chargeVector[6]->push_back(v_p[2]);\r
- chargeVector[7]->push_back(v_pt);\r
- chargeVector[8]->push_back(v_E);\r
-\r
- if(fRunShuffling) {\r
- chargeVectorShuffle[0]->push_back(v_charge);\r
- chargeVectorShuffle[1]->push_back(v_y);\r
- chargeVectorShuffle[2]->push_back(v_eta);\r
- chargeVectorShuffle[3]->push_back(v_phi);\r
- chargeVectorShuffle[4]->push_back(v_p[0]);\r
- chargeVectorShuffle[5]->push_back(v_p[1]);\r
- chargeVectorShuffle[6]->push_back(v_p[2]);\r
- chargeVectorShuffle[7]->push_back(v_pt);\r
- chargeVectorShuffle[8]->push_back(v_E);\r
+ fHistVx->Fill(vertex->GetX());\r
+ fHistVy->Fill(vertex->GetY());\r
+ fHistVz->Fill(vertex->GetZ(),gRefMultiplicity);\r
+ \r
+ // take only events inside centrality class\r
+ if(fUseCentrality) {\r
+ if((gRefMultiplicity > fCentralityPercentileMin) && (gRefMultiplicity < fCentralityPercentileMax)){\r
+ fHistEventStats->Fill(5,gRefMultiplicity); //events with correct centrality\r
+ return gRefMultiplicity; \r
+ }//centrality class\r
+ }\r
+ // take events only within the same multiplicity class\r
+ else if(fUseMultiplicity) {\r
+ //if(fDebugLevel) \r
+ //Printf("N(min): %.0f, N(max): %.0f - N(ref): %.0f",fNumberOfAcceptedTracksMin,\r
+ //fNumberOfAcceptedTracksMax,gRefMultiplicity);\r
+\r
+ if((gRefMultiplicity > fNumberOfAcceptedTracksMin) && (gRefMultiplicity < fNumberOfAcceptedTracksMax)) {\r
+ fHistEventStats->Fill(5,gRefMultiplicity); //events with correct multiplicity\r
+ return gRefMultiplicity;\r
}\r
- \r
- delete track_TPC;\r
- \r
- } //track loop\r
+ }//multiplicity range\r
}//Vz cut\r
}//Vy cut\r
}//Vx cut\r
}//proper number of contributors\r
}//vertex object valid\r
}//triggered event \r
- }//ESD analysis\r
+ }//AOD,ESD,ESDMC\r
\r
- //AOD analysis (vertex and track cuts also here!!!!)\r
- else if(gAnalysisLevel == "AOD") {\r
- AliAODEvent* gAOD = dynamic_cast<AliAODEvent*>(InputEvent()); // from TaskSE\r
- if(!gAOD) {\r
- Printf("ERROR: gAOD not available");\r
- return;\r
+ // in all other cases return -1 (event not accepted)\r
+ return -1;\r
+}\r
+\r
+\r
+//________________________________________________________________________\r
+Double_t AliAnalysisTaskBFPsi::GetRefMultiOrCentrality(AliVEvent *event){\r
+ // Checks the Event cuts\r
+ // Fills Event statistics histograms\r
+ \r
+ Float_t gCentrality = -1.;\r
+ Double_t gMultiplicity = -1.;\r
+ TString gAnalysisLevel = fBalance->GetAnalysisLevel();\r
+\r
+\r
+ // calculate centrality always (not only in centrality mode)\r
+ if(gAnalysisLevel == "AOD"|| gAnalysisLevel == "MCAOD" || gAnalysisLevel == "MCAODrec" ) { //centrality in AOD header //++++++++++++++\r
+ AliAODHeader *header = (AliAODHeader*) event->GetHeader();\r
+ if(header){\r
+ gCentrality = header->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());\r
+\r
+ // QA for centrality estimators\r
+ fHistCentStats->Fill(0.,header->GetCentralityP()->GetCentralityPercentile("V0M"));\r
+ fHistCentStats->Fill(1.,header->GetCentralityP()->GetCentralityPercentile("V0A"));\r
+ fHistCentStats->Fill(2.,header->GetCentralityP()->GetCentralityPercentile("V0C"));\r
+ fHistCentStats->Fill(3.,header->GetCentralityP()->GetCentralityPercentile("FMD"));\r
+ fHistCentStats->Fill(4.,header->GetCentralityP()->GetCentralityPercentile("TRK"));\r
+ fHistCentStats->Fill(5.,header->GetCentralityP()->GetCentralityPercentile("TKL")); \r
+ fHistCentStats->Fill(6.,header->GetCentralityP()->GetCentralityPercentile("CL0"));\r
+ fHistCentStats->Fill(7.,header->GetCentralityP()->GetCentralityPercentile("CL1"));\r
+ fHistCentStats->Fill(8.,header->GetCentralityP()->GetCentralityPercentile("ZNA"));\r
+ fHistCentStats->Fill(9.,header->GetCentralityP()->GetCentralityPercentile("ZPA"));\r
+ fHistCentStats->Fill(10.,header->GetCentralityP()->GetCentralityPercentile("V0MvsFMD"));\r
+ fHistCentStats->Fill(11.,header->GetCentralityP()->GetCentralityPercentile("TKLvsV0M"));\r
+ fHistCentStats->Fill(12.,header->GetCentralityP()->GetCentralityPercentile("ZEMvsZDC"));\r
+ \r
+ // Centrality estimator USED ++++++++++++++++++++++++++++++\r
+ fHistCentStatsUsed->Fill(0.,header->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data()));\r
+ \r
+ // centrality QA (V0M)\r
+ fHistV0M->Fill(event->GetVZEROData()->GetMTotV0A(), event->GetVZEROData()->GetMTotV0C());\r
+ \r
+ // centrality QA (reference tracks)\r
+ fHistRefTracks->Fill(0.,header->GetRefMultiplicity());\r
+ fHistRefTracks->Fill(1.,header->GetRefMultiplicityPos());\r
+ fHistRefTracks->Fill(2.,header->GetRefMultiplicityNeg());\r
+ fHistRefTracks->Fill(3.,header->GetTPConlyRefMultiplicity());\r
+ fHistRefTracks->Fill(4.,header->GetNumberOfITSClusters(0));\r
+ fHistRefTracks->Fill(5.,header->GetNumberOfITSClusters(1));\r
+ fHistRefTracks->Fill(6.,header->GetNumberOfITSClusters(2));\r
+ fHistRefTracks->Fill(7.,header->GetNumberOfITSClusters(3));\r
+ fHistRefTracks->Fill(8.,header->GetNumberOfITSClusters(4));\r
+\r
+ }//AOD header\r
+ }//AOD\r
+ \r
+ else if(gAnalysisLevel == "ESD" || gAnalysisLevel == "MCESD"){ // centrality class for ESDs or MC-ESDs\r
+ AliCentrality *centrality = event->GetCentrality();\r
+ gCentrality = centrality->GetCentralityPercentile(fCentralityEstimator.Data());\r
+\r
+ // QA for centrality estimators\r
+ fHistCentStats->Fill(0.,centrality->GetCentralityPercentile("V0M"));\r
+ fHistCentStats->Fill(1.,centrality->GetCentralityPercentile("V0A"));\r
+ fHistCentStats->Fill(2.,centrality->GetCentralityPercentile("V0C"));\r
+ fHistCentStats->Fill(3.,centrality->GetCentralityPercentile("FMD"));\r
+ fHistCentStats->Fill(4.,centrality->GetCentralityPercentile("TRK"));\r
+ fHistCentStats->Fill(5.,centrality->GetCentralityPercentile("TKL"));\r
+ fHistCentStats->Fill(6.,centrality->GetCentralityPercentile("CL0"));\r
+ fHistCentStats->Fill(7.,centrality->GetCentralityPercentile("CL1"));\r
+ fHistCentStats->Fill(8.,centrality->GetCentralityPercentile("ZNA"));\r
+ fHistCentStats->Fill(9.,centrality->GetCentralityPercentile("ZPA"));\r
+ fHistCentStats->Fill(10.,centrality->GetCentralityPercentile("V0MvsFMD"));\r
+ fHistCentStats->Fill(11.,centrality->GetCentralityPercentile("TKLvsV0M"));\r
+ fHistCentStats->Fill(12.,centrality->GetCentralityPercentile("ZEMvsZDC"));\r
+ \r
+ // Centrality estimator USED ++++++++++++++++++++++++++++++\r
+ fHistCentStatsUsed->Fill(0.,centrality->GetCentralityPercentile(fCentralityEstimator.Data()));\r
+ \r
+ // centrality QA (V0M)\r
+ fHistV0M->Fill(event->GetVZEROData()->GetMTotV0A(), event->GetVZEROData()->GetMTotV0C());\r
+ }//ESD\r
+\r
+ else if(gAnalysisLevel == "MC"){\r
+ Double_t gImpactParameter = 0.;\r
+ AliMCEvent *gMCEvent = dynamic_cast<AliMCEvent*>(event);\r
+ if(gMCEvent){\r
+ AliCollisionGeometry* headerH = dynamic_cast<AliCollisionGeometry*>(gMCEvent->GenEventHeader()); \r
+ if(headerH){\r
+ gImpactParameter = headerH->ImpactParameter();\r
+ gCentrality = gImpactParameter;\r
+ }//MC header\r
+ }//MC event cast\r
+ }//MC\r
+\r
+ else{\r
+ gCentrality = -1.;\r
+ }\r
+ \r
+ // calculate reference multiplicity always (not only in multiplicity mode)\r
+ if(gAnalysisLevel == "ESD" || gAnalysisLevel == "MCESD"){\r
+ AliESDEvent* gESDEvent = dynamic_cast<AliESDEvent*>(event);\r
+ if(gESDEvent){\r
+ gMultiplicity = fESDtrackCuts->GetReferenceMultiplicity(gESDEvent, AliESDtrackCuts::kTrackletsITSTPC,0.5);\r
+ fHistMultiplicity->Fill(gMultiplicity);\r
+ }//AliESDevent cast\r
+ }//ESD mode\r
+\r
+ else if(gAnalysisLevel == "AOD"|| gAnalysisLevel == "MCAOD" || gAnalysisLevel == "MCAODrec" ){\r
+ AliAODHeader *header = (AliAODHeader*) event->GetHeader();\r
+ if ((fMultiplicityEstimator == "V0M")||\r
+ (fMultiplicityEstimator == "V0A")||\r
+ (fMultiplicityEstimator == "V0C") ||\r
+ (fMultiplicityEstimator == "TPC")) {\r
+ gMultiplicity = GetReferenceMultiplicityFromAOD(event);\r
+ if(fDebugLevel) Printf("Reference multiplicity (calculated): %.0f",gMultiplicity);\r
+ }\r
+ else {\r
+ if(header)\r
+ gMultiplicity = header->GetRefMultiplicity();\r
+ if(fDebugLevel) Printf("Reference multiplicity (AOD header): %.0f",gMultiplicity);\r
}\r
+ fHistMultiplicity->Fill(gMultiplicity);\r
+ }//AOD mode\r
+ else if(gAnalysisLevel == "MC") {\r
+ AliMCEvent* gMCEvent = dynamic_cast<AliMCEvent*>(event);\r
+ //Calculating the multiplicity as the number of charged primaries\r
+ //within \pm 0.8 in eta and pT > 0.1 GeV/c\r
+ for(Int_t iParticle = 0; iParticle < gMCEvent->GetNumberOfPrimaries(); iParticle++) {\r
+ AliMCParticle* track = dynamic_cast<AliMCParticle *>(gMCEvent->GetTrack(iParticle));\r
+ if (!track) {\r
+ AliError(Form("Could not receive particle %d", iParticle));\r
+ continue;\r
+ }\r
+ \r
+ //exclude non stable particles\r
+ if(!(gMCEvent->IsPhysicalPrimary(iParticle))) continue;\r
+ \r
+ //++++++++++++++++\r
+ if (fMultiplicityEstimator == "V0M") {\r
+ if((track->Eta() > 5.1 || track->Eta() < 2.8)&&(track->Eta() < -3.7 || track->Eta() > -1.7)) \r
+ continue;}\r
+ else if (fMultiplicityEstimator == "V0A") {\r
+ if(track->Eta() > 5.1 || track->Eta() < 2.8) continue;}\r
+ else if (fMultiplicityEstimator == "V0C") {\r
+ if(track->Eta() > -1.7 || track->Eta() < -3.7) continue;}\r
+ else if (fMultiplicityEstimator == "TPC") {\r
+ if(track->Eta() < fEtaMin || track->Eta() > fEtaMax) continue;\r
+ if(track->Pt() < fPtMin || track->Pt() > fPtMax) continue;\r
+ }\r
+ else{\r
+ if(track->Pt() < fPtMin || track->Pt() > fPtMax) continue;\r
+ if(track->Eta() < fEtaMin || track->Eta() > fEtaMax) continue;\r
+ }\r
+ //++++++++++++++++\r
+ \r
+ if(track->Charge() == 0) continue;\r
+ \r
+ gMultiplicity += 1;\r
+ }//loop over primaries\r
+ fHistMultiplicity->Fill(gMultiplicity);\r
+ }//MC mode\r
+ else{\r
+ gMultiplicity = -1;\r
+ }\r
+ \r
\r
- AliAODHeader *aodHeader = gAOD->GetHeader();\r
+ // decide what should be returned only here\r
+ Double_t lReturnVal = -100;\r
+ if(fEventClass=="Multiplicity"){\r
+ lReturnVal = gMultiplicity;\r
+ }else if(fEventClass=="Centrality"){\r
+ lReturnVal = gCentrality;\r
+ }\r
+ return lReturnVal;\r
+}\r
\r
- // store offline trigger bits\r
- fHistTriggerStats->Fill(aodHeader->GetOfflineTrigger());\r
+//________________________________________________________________________\r
+Double_t AliAnalysisTaskBFPsi::GetReferenceMultiplicityFromAOD(AliVEvent *event){\r
+ //Function that returns the reference multiplicity from AODs (data or reco MC)\r
+ //Different ref. mult. implemented: V0M, V0A, V0C, TPC\r
+ Double_t gRefMultiplicity = 0., gRefMultiplicityTPC = 0.;\r
+ Double_t gRefMultiplicityVZERO = 0., gRefMultiplicityVZEROA = 0., gRefMultiplicityVZEROC = 0.;\r
+\r
+ AliAODHeader *header = dynamic_cast<AliAODHeader *>(event->GetHeader());\r
+ if(!header) {\r
+ Printf("ERROR: AOD header not available");\r
+ return -999;\r
+ }\r
+ Int_t gRunNumber = header->GetRunNumber();\r
+\r
+ // Loop over tracks in event\r
+ for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {\r
+ AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(event->GetTrack(iTracks));\r
+ if (!aodTrack) {\r
+ AliError(Form("Could not receive track %d", iTracks));\r
+ continue;\r
+ }\r
\r
- // event selection done in AliAnalysisTaskSE::Exec() --> this is not used\r
- fHistEventStats->Fill(1,fCentrality); //all events\r
- Bool_t isSelected = kTRUE;\r
- if(fUseOfflineTrigger)\r
- isSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();\r
- if(isSelected) {\r
- fHistEventStats->Fill(2,fCentrality); //triggered events\r
- \r
- //Centrality stuff (centrality in AOD header)\r
- if(fUseCentrality) {\r
- fCentrality = aodHeader->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());\r
- \r
- // QA for centrality estimators\r
- fHistCentStats->Fill(0.,aodHeader->GetCentralityP()->GetCentralityPercentile("V0M"));\r
- fHistCentStats->Fill(1.,aodHeader->GetCentralityP()->GetCentralityPercentile("FMD"));\r
- fHistCentStats->Fill(2.,aodHeader->GetCentralityP()->GetCentralityPercentile("TRK"));\r
- fHistCentStats->Fill(3.,aodHeader->GetCentralityP()->GetCentralityPercentile("TKL"));\r
- fHistCentStats->Fill(4.,aodHeader->GetCentralityP()->GetCentralityPercentile("CL0"));\r
- fHistCentStats->Fill(5.,aodHeader->GetCentralityP()->GetCentralityPercentile("CL1"));\r
- fHistCentStats->Fill(6.,aodHeader->GetCentralityP()->GetCentralityPercentile("V0MvsFMD"));\r
- fHistCentStats->Fill(7.,aodHeader->GetCentralityP()->GetCentralityPercentile("TKLvsV0M"));\r
- fHistCentStats->Fill(8.,aodHeader->GetCentralityP()->GetCentralityPercentile("ZEMvsZDC"));\r
+ // AOD track cuts\r
+ if(!aodTrack->TestFilterBit(fnAODtrackCutBit)) continue;\r
+ \r
+ if(aodTrack->Charge() == 0) continue;\r
+ // Kinematics cuts from ESD track cuts\r
+ if( aodTrack->Pt() < fPtMin || aodTrack->Pt() > fPtMax) continue;\r
+ if( aodTrack->Eta() < fEtaMin || aodTrack->Eta() > fEtaMax) continue;\r
+ \r
+ //=================PID (so far only for electron rejection)==========================//\r
+ if(fElectronRejection) {\r
+ // get the electron nsigma\r
+ Double_t nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kElectron));\r
+ \r
+ // check only for given momentum range\r
+ if( aodTrack->Pt() > fElectronRejectionMinPt && aodTrack->Pt() < fElectronRejectionMaxPt ){\r
+ //look only at electron nsigma\r
+ if(!fElectronOnlyRejection) {\r
+ //Make the decision based on the n-sigma of electrons\r
+ if(nSigma < fElectronRejectionNSigma) continue;\r
+ }\r
+ else {\r
+ Double_t nSigmaPions = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kPion));\r
+ Double_t nSigmaKaons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kKaon));\r
+ Double_t nSigmaProtons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kProton));\r
+ \r
+ //Make the decision based on the n-sigma of electrons exclusively ( = track not in nsigma region for other species)\r
+ if(nSigma < fElectronRejectionNSigma\r
+ && nSigmaPions > fElectronRejectionNSigma\r
+ && nSigmaKaons > fElectronRejectionNSigma\r
+ && nSigmaProtons > fElectronRejectionNSigma ) continue;\r
+ }\r
+ }\r
+ }//electron rejection\r
+ \r
+ gRefMultiplicityTPC += 1.0;\r
+ }// track loop\r
+ \r
+ //VZERO segmentation in two detectors (0-31: VZERO-C, 32-63: VZERO-A)\r
+ for(Int_t iChannel = 0; iChannel < 64; iChannel++) {\r
+ fHistVZEROSignal->Fill(iChannel,event->GetVZEROEqMultiplicity(iChannel));\r
+ \r
+ if(iChannel < 32) \r
+ gRefMultiplicityVZEROC += event->GetVZEROEqMultiplicity(iChannel);\r
+ else if(iChannel >= 32) \r
+ gRefMultiplicityVZEROA += event->GetVZEROEqMultiplicity(iChannel);\r
+ }//loop over PMTs\r
+ \r
+ //Equalization of gain\r
+ Double_t gFactorA = GetEqualizationFactor(gRunNumber,"A");\r
+ if(gFactorA != 0)\r
+ gRefMultiplicityVZEROA /= gFactorA;\r
+ Double_t gFactorC = GetEqualizationFactor(gRunNumber,"C");\r
+ if(gFactorC != 0)\r
+ gRefMultiplicityVZEROC /= gFactorC;\r
+ if((gFactorA != 0)&&(gFactorC != 0)) \r
+ gRefMultiplicityVZERO = (gRefMultiplicityVZEROA/gFactorA)+(gRefMultiplicityVZEROC/gFactorC);\r
+ \r
+ if(fDebugLevel) \r
+ Printf("VZERO multiplicity: %.0f - TPC multiplicity: %.0f",gRefMultiplicityVZERO,gRefMultiplicityTPC);\r
+\r
+ fHistTPCvsVZEROMultiplicity->Fill(gRefMultiplicityVZERO,gRefMultiplicityTPC);\r
+\r
+ if(fMultiplicityEstimator == "TPC") \r
+ gRefMultiplicity = gRefMultiplicityTPC;\r
+ else if(fMultiplicityEstimator == "V0M")\r
+ gRefMultiplicity = gRefMultiplicityVZERO;\r
+ else if(fMultiplicityEstimator == "V0A")\r
+ gRefMultiplicity = gRefMultiplicityVZEROA;\r
+ else if(fMultiplicityEstimator == "V0C")\r
+ gRefMultiplicity = gRefMultiplicityVZEROC;\r
+ \r
+ return gRefMultiplicity;\r
+}\r
+\r
+//________________________________________________________________________\r
+Double_t AliAnalysisTaskBFPsi::GetEventPlane(AliVEvent *event){\r
+ // Get the event plane\r
+\r
+ TString gAnalysisLevel = fBalance->GetAnalysisLevel();\r
+\r
+ Float_t gVZEROEventPlane = -10.;\r
+ Float_t gReactionPlane = -10.;\r
+ Double_t qxTot = 0.0, qyTot = 0.0;\r
+\r
+ //MC: from reaction plane\r
+ if(gAnalysisLevel == "MC"){\r
+ if(!event) {\r
+ AliError("mcEvent not available");\r
+ return 0x0;\r
+ }\r
+\r
+ AliMCEvent *gMCEvent = dynamic_cast<AliMCEvent*>(event);\r
+ if(gMCEvent){\r
+ AliCollisionGeometry* headerH = dynamic_cast<AliCollisionGeometry*>(gMCEvent->GenEventHeader()); \r
+ if (headerH) {\r
+ gReactionPlane = headerH->ReactionPlaneAngle();\r
+ //gReactionPlane *= TMath::RadToDeg();\r
+ }//MC header\r
+ }//MC event cast\r
+ }//MC\r
+ \r
+ // AOD,ESD,ESDMC: from VZERO Event Plane\r
+ else{\r
+ \r
+ AliEventplane *ep = event->GetEventplane();\r
+ if(ep){ \r
+ gVZEROEventPlane = ep->CalculateVZEROEventPlane(event,10,2,qxTot,qyTot);\r
+ if(gVZEROEventPlane < 0.) gVZEROEventPlane += TMath::Pi();\r
+ //gReactionPlane = gVZEROEventPlane*TMath::RadToDeg();\r
+ gReactionPlane = gVZEROEventPlane;\r
+ }\r
+ }//AOD,ESD,ESDMC\r
+\r
+ return gReactionPlane;\r
+}\r
+\r
+//________________________________________________________________________\r
+Double_t AliAnalysisTaskBFPsi::GetTrackbyTrackCorrectionMatrix( Double_t vEta, \r
+ Double_t vPhi, \r
+ Double_t vPt, \r
+ Short_t vCharge, \r
+ Double_t gCentrality) {\r
+ // -- Get efficiency correction of particle dependent on (eta, phi, pt, charge, centrality) \r
+\r
+ Double_t correction = 1.;\r
+ Int_t binEta = 0, binPt = 0, binPhi = 0;\r
+\r
+ //Printf("EtaMAx: %lf - EtaMin: %lf - EtaBin: %lf", fEtaMaxForCorrections,fEtaMinForCorrections,fEtaBinForCorrections);\r
+ if(fEtaBinForCorrections != 0) {\r
+ Double_t widthEta = (fEtaMaxForCorrections - fEtaMinForCorrections)/fEtaBinForCorrections;\r
+ if(fEtaMaxForCorrections != fEtaMinForCorrections) \r
+ binEta = (Int_t)((vEta-fEtaMinForCorrections)/widthEta)+1;\r
+ }\r
+\r
+ if(fPtBinForCorrections != 0) {\r
+ Double_t widthPt = (fPtMaxForCorrections - fPtMinForCorrections)/fPtBinForCorrections;\r
+ if(fPtMaxForCorrections != fPtMinForCorrections) \r
+ binPt = (Int_t)((vPt-fPtMinForCorrections)/widthPt) + 1;\r
+ }\r
+ \r
+ if(fPhiBinForCorrections != 0) {\r
+ Double_t widthPhi = (fPhiMaxForCorrections - fPhiMinForCorrections)/fPhiBinForCorrections;\r
+ if(fPhiMaxForCorrections != fPhiMinForCorrections) \r
+ binPhi = (Int_t)((vPhi-fPhiMinForCorrections)/widthPhi)+ 1;\r
+ }\r
+\r
+ Int_t gCentralityInt = -1;\r
+ for (Int_t i=0; i<fCentralityArrayBinsForCorrections-1; i++){\r
+ if((fCentralityArrayForCorrections[i] <= gCentrality)&&(gCentrality <= fCentralityArrayForCorrections[i+1])){\r
+ gCentralityInt = i;\r
+ break;\r
+ }\r
+ } \r
+\r
+ // centrality not in array --> no correction\r
+ if(gCentralityInt < 0){\r
+ correction = 1.;\r
+ }\r
+ else{\r
+ \r
+ //Printf("//=============CENTRALITY=============// %d:",gCentralityInt);\r
+ \r
+ if(fHistCorrectionPlus[gCentralityInt]){\r
+ if (vCharge > 0) {\r
+ correction = fHistCorrectionPlus[gCentralityInt]->GetBinContent(fHistCorrectionPlus[gCentralityInt]->GetBin(binEta, binPt, binPhi));\r
+ //Printf("CORRECTIONplus: %.2f | Centrality %d",correction,gCentralityInt); \r
+ }\r
+ if (vCharge < 0) {\r
+ correction = fHistCorrectionMinus[gCentralityInt]->GetBinContent(fHistCorrectionMinus[gCentralityInt]->GetBin(binEta, binPt, binPhi));\r
+ //Printf("CORRECTIONminus: %.2f | Centrality %d",correction,gCentralityInt);\r
+ }\r
+ }\r
+ else {\r
+ correction = 1.;\r
+ }\r
+ }//centrality in array\r
+ \r
+ if (correction == 0.) { \r
+ AliError(Form("Should not happen : bin content = 0. >> eta: %.2f | phi : %.2f | pt : %.2f | cent %d",vEta, vPhi, vPt, gCentralityInt)); \r
+ correction = 1.; \r
+ } \r
+ \r
+ return correction;\r
+}\r
+\r
+//________________________________________________________________________\r
+TObjArray* AliAnalysisTaskBFPsi::GetAcceptedTracks(AliVEvent *event, Double_t gCentrality, Double_t gReactionPlane){\r
+ // Returns TObjArray with tracks after all track cuts (only for AOD!)\r
+ // Fills QA histograms\r
+\r
+ TString gAnalysisLevel = fBalance->GetAnalysisLevel();\r
+\r
+ //output TObjArray holding all good tracks\r
+ TObjArray* tracksAccepted = new TObjArray;\r
+ tracksAccepted->SetOwner(kTRUE);\r
+\r
+ Short_t vCharge;\r
+ Double_t vEta;\r
+ Double_t vY;\r
+ Double_t vPhi;\r
+ Double_t vPt;\r
+\r
+\r
+ if(gAnalysisLevel == "AOD") { // handling of TPC only tracks different in AOD and ESD\r
+ // Loop over tracks in event\r
+ \r
+ for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {\r
+ AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(event->GetTrack(iTracks));\r
+ if (!aodTrack) {\r
+ AliError(Form("Could not receive track %d", iTracks));\r
+ continue;\r
+ }\r
+ \r
+ // AOD track cuts\r
+ \r
+ // For ESD Filter Information: ANALYSIS/macros/AddTaskESDfilter.C\r
+ // take only TPC only tracks \r
+ //fHistTrackStats->Fill(aodTrack->GetFilterMap());\r
+ for(Int_t iTrackBit = 0; iTrackBit < 16; iTrackBit++){\r
+ fHistTrackStats->Fill(iTrackBit,aodTrack->TestFilterBit(1<<iTrackBit));\r
+ }\r
+\r
+ if(!aodTrack->TestFilterBit(fnAODtrackCutBit)) continue;\r
+ \r
+ vCharge = aodTrack->Charge();\r
+ vEta = aodTrack->Eta();\r
+ vY = aodTrack->Y();\r
+ vPhi = aodTrack->Phi();// * TMath::RadToDeg();\r
+ vPt = aodTrack->Pt();\r
+ \r
+ //===========================PID (so far only for electron rejection)===============================// \r
+ if(fElectronRejection) {\r
+\r
+ // get the electron nsigma\r
+ Double_t nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kElectron));\r
\r
- // take only events inside centrality class\r
- if((fCentrality < fCentralityPercentileMin) || (fCentrality > fCentralityPercentileMax)) \r
- return;\r
+ //Fill QA before the PID\r
+ fHistdEdxVsPTPCbeforePIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());\r
+ fHistNSigmaTPCvsPtbeforePIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),nSigma); \r
+ //end of QA-before pid\r
\r
- // centrality QA (V0M)\r
- fHistV0M->Fill(gAOD->GetVZEROData()->GetMTotV0A(), gAOD->GetVZEROData()->GetMTotV0C());\r
+ // check only for given momentum range\r
+ if( vPt > fElectronRejectionMinPt && vPt < fElectronRejectionMaxPt ){\r
+ \r
+ //look only at electron nsigma\r
+ if(!fElectronOnlyRejection){\r
+ \r
+ //Make the decision based on the n-sigma of electrons\r
+ if(nSigma < fElectronRejectionNSigma) continue;\r
+ }\r
+ else{\r
+ \r
+ Double_t nSigmaPions = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kPion));\r
+ Double_t nSigmaKaons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kKaon));\r
+ Double_t nSigmaProtons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kProton));\r
+ \r
+ //Make the decision based on the n-sigma of electrons exclusively ( = track not in nsigma region for other species)\r
+ if(nSigma < fElectronRejectionNSigma\r
+ && nSigmaPions > fElectronRejectionNSigma\r
+ && nSigmaKaons > fElectronRejectionNSigma\r
+ && nSigmaProtons > fElectronRejectionNSigma ) continue;\r
+ }\r
+ }\r
+ \r
+ //Fill QA after the PID\r
+ fHistdEdxVsPTPCafterPIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());\r
+ fHistNSigmaTPCvsPtafterPIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),nSigma); \r
\r
- // centrality QA (reference tracks)\r
- fHistRefTracks->Fill(0.,aodHeader->GetRefMultiplicity());\r
- fHistRefTracks->Fill(1.,aodHeader->GetRefMultiplicityPos());\r
- fHistRefTracks->Fill(2.,aodHeader->GetRefMultiplicityNeg());\r
- fHistRefTracks->Fill(3.,aodHeader->GetTPConlyRefMultiplicity());\r
- fHistRefTracks->Fill(4.,aodHeader->GetNumberOfITSClusters(0));\r
- fHistRefTracks->Fill(5.,aodHeader->GetNumberOfITSClusters(1));\r
- fHistRefTracks->Fill(6.,aodHeader->GetNumberOfITSClusters(2));\r
- fHistRefTracks->Fill(7.,aodHeader->GetNumberOfITSClusters(3));\r
- fHistRefTracks->Fill(8.,aodHeader->GetNumberOfITSClusters(4));\r
}\r
+ //===========================end of PID (so far only for electron rejection)===============================//\r
+\r
+ //+++++++++++++++++++++++++++++//\r
+ //===========================PID===============================// \r
+ if(fUsePID) {\r
+ Double_t prob[AliPID::kSPECIES]={0.};\r
+ Double_t probTPC[AliPID::kSPECIES]={0.};\r
+ Double_t probTOF[AliPID::kSPECIES]={0.};\r
+ Double_t probTPCTOF[AliPID::kSPECIES]={0.};\r
+ \r
+ Double_t nSigma = 0.;\r
+ Double_t nSigmaTPC = 0.; //++++\r
+ Double_t nSigmaTOF = 0.; //++++\r
+ Double_t nSigmaTPCTOF = 0.; //++++\r
+ UInt_t detUsedTPC = 0;\r
+ UInt_t detUsedTOF = 0;\r
+ UInt_t detUsedTPCTOF = 0;\r
+ \r
+ nSigmaTPC = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)fParticleOfInterest));\r
+ nSigmaTOF = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(aodTrack,(AliPID::EParticleType)fParticleOfInterest));\r
+ nSigmaTPCTOF = TMath::Sqrt(nSigmaTPC*nSigmaTPC + nSigmaTOF*nSigmaTOF);\r
+\r
+ //Decide what detector configuration we want to use\r
+ switch(fPidDetectorConfig) {\r
+ case kTPCpid:\r
+ fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);\r
+ nSigma = nSigmaTPC;\r
+ detUsedTPC = fPIDCombined->ComputeProbabilities(aodTrack, fPIDResponse, probTPC);\r
+ for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)\r
+ prob[iSpecies] = probTPC[iSpecies];\r
+ break;\r
+ case kTOFpid:\r
+ fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF);\r
+ nSigma = nSigmaTOF;\r
+ detUsedTOF = fPIDCombined->ComputeProbabilities(aodTrack, fPIDResponse, probTOF);\r
+ for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)\r
+ prob[iSpecies] = probTOF[iSpecies];\r
+ break;\r
+ case kTPCTOF:\r
+ fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC);\r
+ nSigma = nSigmaTPCTOF;\r
+ detUsedTPCTOF = fPIDCombined->ComputeProbabilities(aodTrack, fPIDResponse, probTPCTOF);\r
+ for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)\r
+ prob[iSpecies] = probTPCTOF[iSpecies];\r
+ break;\r
+ default:\r
+ break;\r
+ }//end switch: define detector mask\r
+ \r
+ //Filling the PID QA\r
+ Double_t tofTime = -999., length = 999., tof = -999.;\r
+ Double_t c = TMath::C()*1.E-9;// m/ns\r
+ Double_t beta = -999.;\r
+ if ( (aodTrack->IsOn(AliAODTrack::kTOFin)) &&\r
+ (aodTrack->IsOn(AliAODTrack::kTIME)) ) { \r
+ tofTime = aodTrack->GetTOFsignal();//in ps\r
+ length = aodTrack->GetIntegratedLength();\r
+ tof = tofTime*1E-3; // ns \r
+ \r
+ if (tof <= 0) {\r
+ //Printf("WARNING: track with negative TOF time found! Skipping this track for PID checks\n");\r
+ continue;\r
+ }\r
+ if (length <= 0){\r
+ //printf("WARNING: track with negative length found!Skipping this track for PID checks\n");\r
+ continue;\r
+ }\r
+ \r
+ length = length*0.01; // in meters\r
+ tof = tof*c;\r
+ beta = length/tof;\r
+ \r
+ fHistBetavsPTOFbeforePID ->Fill(aodTrack->P()*aodTrack->Charge(),beta);\r
+ fHistProbTOFvsPtbeforePID ->Fill(aodTrack->Pt(),probTOF[fParticleOfInterest]);\r
+ fHistNSigmaTOFvsPtbeforePID ->Fill(aodTrack->Pt(),nSigmaTOF);\r
+ }//TOF signal \r
+ \r
+ fHistdEdxVsPTPCbeforePID -> Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());\r
+ fHistProbTPCvsPtbeforePID -> Fill(aodTrack->Pt(),probTPC[fParticleOfInterest]); \r
+ fHistNSigmaTPCvsPtbeforePID -> Fill(aodTrack->Pt(),nSigmaTPC);\r
+ \r
+ fHistProbTPCTOFvsPtbeforePID -> Fill(aodTrack->Pt(),probTPCTOF[fParticleOfInterest]);\r
\r
- const AliAODVertex *vertex = gAOD->GetPrimaryVertex();\r
- \r
- if(vertex) {\r
- Double32_t fCov[6];\r
- vertex->GetCovarianceMatrix(fCov);\r
+ //combined TPC&TOF\r
+ fHistBetaVsdEdXbeforePID->Fill(aodTrack->GetTPCsignal(),beta); //+++++++++ \r
+ fHistNSigmaTPCTOFvsPtbeforePID -> Fill(aodTrack->Pt(),nSigmaTPCTOF);\r
+ Printf("NSIGMA: %lf - nSigmaTOF: %lf - nSigmaTPC: %lf - nSigmaTPCTOF: %lf",nSigma,nSigmaTOF,nSigmaTPC,nSigmaTPCTOF);\r
\r
- if(vertex->GetNContributors() > 0) {\r
- if(fCov[5] != 0) {\r
- fHistEventStats->Fill(3,fCentrality); //events with a proper vertex\r
- if(TMath::Abs(vertex->GetX()) < fVxMax) {\r
- if(TMath::Abs(vertex->GetY()) < fVyMax) {\r
- if(TMath::Abs(vertex->GetZ()) < fVzMax) {\r
- fHistEventStats->Fill(4,fCentrality); //analyzed events\r
- fHistVx->Fill(vertex->GetX());\r
- fHistVy->Fill(vertex->GetY());\r
- fHistVz->Fill(vertex->GetZ(),fCentrality);\r
+ //end of QA-before pid\r
\r
- //========Get the VZERO event plane========//\r
- Double_t gVZEROEventPlane = -10.0;\r
- Double_t qxTot = 0.0, qyTot = 0.0;\r
- AliEventplane *ep = gAOD->GetEventplane();\r
- if(ep) \r
- gVZEROEventPlane = ep->CalculateVZEROEventPlane(gAOD,10,2,qxTot,qyTot);\r
- if(gVZEROEventPlane < 0.) gVZEROEventPlane += TMath::Pi();\r
- gReactionPlane = gVZEROEventPlane*TMath::RadToDeg();\r
- fHistEventPlane->Fill(gReactionPlane,fCentrality);\r
- //========Get the VZERO event plane========//\r
-\r
- //Printf("There are %d tracks in this event", gAOD->GetNumberOfTracks());\r
- for (Int_t iTracks = 0; iTracks < gAOD->GetNumberOfTracks(); iTracks++) {\r
- AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(gAOD->GetTrack(iTracks));\r
- if (!aodTrack) {\r
- Printf("ERROR: Could not receive track %d", iTracks);\r
- continue;\r
- }\r
- \r
- // AOD track cuts\r
- \r
- // For ESD Filter Information: ANALYSIS/macros/AddTaskESDfilter.C\r
- // take only TPC only tracks \r
- fHistTrackStats->Fill(aodTrack->GetFilterMap());\r
- if(!aodTrack->TestFilterBit(nAODtrackCutBit)) continue;\r
- \r
- v_charge = aodTrack->Charge();\r
- v_y = aodTrack->Y();\r
- v_eta = aodTrack->Eta();\r
- v_phi = aodTrack->Phi() * TMath::RadToDeg();\r
- v_E = aodTrack->E();\r
- v_pt = aodTrack->Pt();\r
- aodTrack->PxPyPz(v_p);\r
- \r
- Float_t DCAxy = aodTrack->DCA(); // this is the DCA from global track (not exactly what is cut on)\r
- Float_t DCAz = aodTrack->ZAtDCA(); // this is the DCA from global track (not exactly what is cut on)\r
- \r
- \r
- // Kinematics cuts from ESD track cuts\r
- if( v_pt < fPtMin || v_pt > fPtMax) continue;\r
- if (!fUsePID) {\r
- if( v_eta < fEtaMin || v_eta > fEtaMax) continue;\r
- }\r
- else if (fUsePID){\r
- if( v_y < fEtaMin || v_y > fEtaMax) continue;\r
- }\r
+ if ((detUsedTPC != 0)||(detUsedTOF != 0)||(detUsedTPCTOF != 0)) {\r
+ //Make the decision based on the n-sigma\r
+ if(fUsePIDnSigma) {\r
+ if(nSigma > fPIDNSigma) continue; \r
+ \r
+ fHistNSigmaTOFvsPtafterPID ->Fill(aodTrack->Pt(),nSigmaTOF);\r
+ fHistNSigmaTPCvsPtafterPID ->Fill(aodTrack->Pt(),nSigmaTPC);\r
+ fHistNSigmaTPCTOFvsPtafterPID ->Fill(aodTrack->Pt(),nSigmaTPCTOF);\r
+ }\r
+ //Make the decision based on the bayesian\r
+ else if(fUsePIDPropabilities) {\r
+ if(fParticleOfInterest != TMath::LocMax(AliPID::kSPECIES,prob)) continue;\r
+ if (prob[fParticleOfInterest] < fMinAcceptedPIDProbability) continue; \r
+ \r
+ fHistProbTOFvsPtafterPID ->Fill(aodTrack->Pt(),probTOF[fParticleOfInterest]);\r
+ fHistProbTPCvsPtafterPID ->Fill(aodTrack->Pt(),probTPC[fParticleOfInterest]); \r
+ fHistProbTPCTOFvsPtafterPID ->Fill(aodTrack->Pt(),probTPCTOF[fParticleOfInterest]);\r
+ \r
+ }\r
+ \r
+ //Fill QA after the PID\r
+ fHistBetavsPTOFafterPID ->Fill(aodTrack->P()*aodTrack->Charge(),beta);\r
+ fHistdEdxVsPTPCafterPID ->Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());\r
+ fHistBetaVsdEdXafterPID->Fill(aodTrack->GetTPCsignal(),beta); //+++++++++ \r
+ }\r
+ }\r
+ //===========================PID===============================//\r
+ //+++++++++++++++++++++++++++++//\r
\r
- // Extra DCA cuts (for systematic studies [!= -1])\r
- if( fDCAxyCut != -1 && fDCAzCut != -1){\r
- if(TMath::Sqrt((DCAxy*DCAxy)/(fDCAxyCut*fDCAxyCut)+(DCAz*DCAz)/(fDCAzCut*fDCAzCut)) > 1 ){\r
- continue; // 2D cut\r
- }\r
- }\r
- \r
- // Extra TPC cuts (for systematic studies [!= -1])\r
- if( fTPCchi2Cut != -1 && aodTrack->Chi2perNDF() > fTPCchi2Cut){\r
- continue;\r
- }\r
- if( fNClustersTPCCut != -1 && aodTrack->GetTPCNcls() < fNClustersTPCCut){\r
- continue;\r
- }\r
- \r
- // fill QA histograms\r
- fHistClus->Fill(aodTrack->GetITSNcls(),aodTrack->GetTPCNcls());\r
- fHistDCA->Fill(DCAz,DCAxy);\r
- fHistChi2->Fill(aodTrack->Chi2perNDF(),fCentrality);\r
- fHistPt->Fill(v_pt,fCentrality);\r
- fHistEta->Fill(v_eta,fCentrality);\r
- fHistPhi->Fill(v_phi,fCentrality);\r
- fHistRapidity->Fill(v_y,fCentrality);\r
- if(v_charge > 0) fHistPhiPos->Fill(v_phi,fCentrality);\r
- else if(v_charge < 0) fHistPhiNeg->Fill(v_phi,fCentrality);\r
-\r
- // fill charge vector\r
- chargeVector[0]->push_back(v_charge);\r
- chargeVector[1]->push_back(v_y);\r
- chargeVector[2]->push_back(v_eta);\r
- chargeVector[3]->push_back(v_phi);\r
- chargeVector[4]->push_back(v_p[0]);\r
- chargeVector[5]->push_back(v_p[1]);\r
- chargeVector[6]->push_back(v_p[2]);\r
- chargeVector[7]->push_back(v_pt);\r
- chargeVector[8]->push_back(v_E);\r
-\r
- if(fRunShuffling) {\r
- chargeVectorShuffle[0]->push_back(v_charge);\r
- chargeVectorShuffle[1]->push_back(v_y);\r
- chargeVectorShuffle[2]->push_back(v_eta);\r
- chargeVectorShuffle[3]->push_back(v_phi);\r
- chargeVectorShuffle[4]->push_back(v_p[0]);\r
- chargeVectorShuffle[5]->push_back(v_p[1]);\r
- chargeVectorShuffle[6]->push_back(v_p[2]);\r
- chargeVectorShuffle[7]->push_back(v_pt);\r
- chargeVectorShuffle[8]->push_back(v_E);\r
- }\r
- \r
- gNumberOfAcceptedTracks += 1;\r
- \r
- } //track loop\r
- }//Vz cut\r
- }//Vy cut\r
- }//Vx cut\r
- }//proper vertex resolution\r
- }//proper number of contributors\r
- }//vertex object valid\r
- }//triggered event \r
- }//AOD analysis\r
\r
- //MC-ESD analysis\r
- if(gAnalysisLevel == "MCESD") {\r
- AliMCEvent* mcEvent = MCEvent(); \r
+ Float_t dcaXY = aodTrack->DCA(); // this is the DCA from global track (not exactly what is cut on)\r
+ Float_t dcaZ = aodTrack->ZAtDCA(); // this is the DCA from global track (not exactly what is cut on)\r
+ \r
+ \r
+ // Kinematics cuts from ESD track cuts\r
+ if( vPt < fPtMin || vPt > fPtMax) continue;\r
+ if( vEta < fEtaMin || vEta > fEtaMax) continue;\r
+ \r
+ // Extra DCA cuts (for systematic studies [!= -1])\r
+ if( fDCAxyCut != -1 && fDCAzCut != -1){\r
+ if(TMath::Sqrt((dcaXY*dcaXY)/(fDCAxyCut*fDCAxyCut)+(dcaZ*dcaZ)/(fDCAzCut*fDCAzCut)) > 1 ){\r
+ continue; // 2D cut\r
+ }\r
+ }\r
+ \r
+ // Extra TPC cuts (for systematic studies [!= -1])\r
+ if( fTPCchi2Cut != -1 && aodTrack->Chi2perNDF() > fTPCchi2Cut){\r
+ continue;\r
+ }\r
+ if( fNClustersTPCCut != -1 && aodTrack->GetTPCNcls() < fNClustersTPCCut){\r
+ continue;\r
+ }\r
+ \r
+ // fill QA histograms\r
+ fHistClus->Fill(aodTrack->GetITSNcls(),aodTrack->GetTPCNcls());\r
+ fHistDCA->Fill(dcaZ,dcaXY);\r
+ fHistChi2->Fill(aodTrack->Chi2perNDF(),gCentrality);\r
+ fHistPt->Fill(vPt,gCentrality);\r
+ fHistEta->Fill(vEta,gCentrality);\r
+ fHistRapidity->Fill(vY,gCentrality);\r
+ if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);\r
+ else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);\r
+ fHistPhi->Fill(vPhi,gCentrality);\r
+ if(vCharge > 0) fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality); \r
+ else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);\r
+ \r
+ //=======================================correction\r
+ Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality); \r
+ //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);\r
+ \r
+ // add the track to the TObjArray\r
+ tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction)); \r
+ }//track loop\r
+ }// AOD analysis\r
+\r
+ //==============================================================================================================\r
+ else if(gAnalysisLevel == "MCAOD") {\r
+ \r
+ AliMCEvent* mcEvent = MCEvent();\r
if (!mcEvent) {\r
- Printf("ERROR: mcEvent not available");\r
- return;\r
+ AliError("ERROR: Could not retrieve MC event");\r
}\r
+ else{\r
+ \r
+ for (Int_t iTracks = 0; iTracks < mcEvent->GetNumberOfTracks(); iTracks++) {\r
+ AliAODMCParticle *aodTrack = (AliAODMCParticle*) mcEvent->GetTrack(iTracks); \r
+ if (!aodTrack) {\r
+ AliError(Form("ERROR: Could not receive track %d (mc loop)", iTracks));\r
+ continue;\r
+ }\r
+ \r
+ if(!aodTrack->IsPhysicalPrimary()) continue; \r
+ \r
+ vCharge = aodTrack->Charge();\r
+ vEta = aodTrack->Eta();\r
+ vY = aodTrack->Y();\r
+ vPhi = aodTrack->Phi();// * TMath::RadToDeg();\r
+ vPt = aodTrack->Pt();\r
+ \r
+ // Kinematics cuts from ESD track cuts\r
+ if( vPt < fPtMin || vPt > fPtMax) continue;\r
+ if( vEta < fEtaMin || vEta > fEtaMax) continue;\r
+ \r
+ // Remove neutral tracks\r
+ if( vCharge == 0 ) continue;\r
+ \r
+ //Exclude resonances\r
+ if(fExcludeResonancesInMC) {\r
+ \r
+ Bool_t kExcludeParticle = kFALSE;\r
+ Int_t gMotherIndex = aodTrack->GetMother();\r
+ if(gMotherIndex != -1) {\r
+ AliAODMCParticle* motherTrack = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(gMotherIndex));\r
+ if(motherTrack) {\r
+ Int_t pdgCodeOfMother = motherTrack->GetPdgCode();\r
+ //if((pdgCodeOfMother == 113)||(pdgCodeOfMother == 213)||(pdgCodeOfMother == 221)||(pdgCodeOfMother == 223)||(pdgCodeOfMother == 331)||(pdgCodeOfMother == 333)) {\r
+ //if(pdgCodeOfMother == 113) {\r
+ if(pdgCodeOfMother == 113 // rho0\r
+ || pdgCodeOfMother == 213 || pdgCodeOfMother == -213 // rho+\r
+ // || pdgCodeOfMother == 221 // eta\r
+ // || pdgCodeOfMother == 331 // eta'\r
+ // || pdgCodeOfMother == 223 // omega\r
+ // || pdgCodeOfMother == 333 // phi\r
+ || pdgCodeOfMother == 311 || pdgCodeOfMother == -311 // K0\r
+ // || pdgCodeOfMother == 313 || pdgCodeOfMother == -313 // K0*\r
+ // || pdgCodeOfMother == 323 || pdgCodeOfMother == -323 // K+*\r
+ || pdgCodeOfMother == 3122 || pdgCodeOfMother == -3122 // Lambda\r
+ || pdgCodeOfMother == 111 // pi0 Dalitz\r
+ || pdgCodeOfMother == 22 // photon\r
+ ) {\r
+ kExcludeParticle = kTRUE;\r
+ }\r
+ }\r
+ }\r
+ \r
+ //Exclude from the analysis decay products of rho0, rho+, eta, eta' and phi\r
+ if(kExcludeParticle) continue;\r
+ }\r
\r
- AliESDEvent* gESD = dynamic_cast<AliESDEvent*>(InputEvent()); // from TaskSE\r
- if (!gESD) {\r
- Printf("ERROR: gESD not available");\r
+ //Exclude electrons with PDG\r
+ if(fExcludeElectronsInMC) {\r
+ \r
+ if(TMath::Abs(aodTrack->GetPdgCode()) == 11) continue;\r
+ \r
+ }\r
+ \r
+ // fill QA histograms\r
+ fHistPt->Fill(vPt,gCentrality);\r
+ fHistEta->Fill(vEta,gCentrality);\r
+ fHistRapidity->Fill(vY,gCentrality);\r
+ if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);\r
+ else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);\r
+ fHistPhi->Fill(vPhi,gCentrality);\r
+ if(vCharge > 0) fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality); \r
+ else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);\r
+ \r
+ //=======================================correction\r
+ Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality); \r
+ //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality); \r
+ \r
+ // add the track to the TObjArray\r
+ tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction)); \r
+ }//aodTracks\r
+ }//MC event\r
+ }//MCAOD\r
+ //==============================================================================================================\r
+\r
+ //==============================================================================================================\r
+ else if(gAnalysisLevel == "MCAODrec") {\r
+ \r
+ /* fAOD = dynamic_cast<AliAODEvent*>(InputEvent());\r
+ if (!fAOD) {\r
+ printf("ERROR: fAOD not available\n");\r
return;\r
+ }*/\r
+ \r
+ fArrayMC = dynamic_cast<TClonesArray*>(event->FindListObject(AliAODMCParticle::StdBranchName())); \r
+ if (!fArrayMC) { \r
+ AliError("No array of MC particles found !!!");\r
}\r
\r
- // store offline trigger bits\r
- fHistTriggerStats->Fill(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected());\r
+ AliMCEvent* mcEvent = MCEvent();\r
+ if (!mcEvent) {\r
+ AliError("ERROR: Could not retrieve MC event");\r
+ return tracksAccepted;\r
+ }\r
+ \r
+ for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {\r
+ AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(event->GetTrack(iTracks));\r
+ if (!aodTrack) {\r
+ AliError(Form("Could not receive track %d", iTracks));\r
+ continue;\r
+ }\r
+ \r
+ for(Int_t iTrackBit = 0; iTrackBit < 16; iTrackBit++){\r
+ fHistTrackStats->Fill(iTrackBit,aodTrack->TestFilterBit(1<<iTrackBit));\r
+ }\r
+ if(!aodTrack->TestFilterBit(fnAODtrackCutBit)) continue;\r
+ \r
+ vCharge = aodTrack->Charge();\r
+ vEta = aodTrack->Eta();\r
+ vY = aodTrack->Y();\r
+ vPhi = aodTrack->Phi();// * TMath::RadToDeg();\r
+ vPt = aodTrack->Pt();\r
+ \r
+ //===========================use MC information for Kinematics===============================// \r
+ if(fUseMCforKinematics){\r
+\r
+ Int_t label = TMath::Abs(aodTrack->GetLabel());\r
+ AliAODMCParticle *AODmcTrack = (AliAODMCParticle*) fArrayMC->At(label);\r
+\r
+ if(AODmcTrack){\r
+ vCharge = AODmcTrack->Charge();\r
+ vEta = AODmcTrack->Eta();\r
+ vY = AODmcTrack->Y();\r
+ vPhi = AODmcTrack->Phi();// * TMath::RadToDeg();\r
+ vPt = AODmcTrack->Pt();\r
+ }\r
+ else{\r
+ AliDebug(1, "no MC particle for this track"); \r
+ continue;\r
+ }\r
+ }\r
+ //===========================end of use MC information for Kinematics========================// \r
+\r
\r
- // event selection done in AliAnalysisTaskSE::Exec() --> this is not used\r
- fHistEventStats->Fill(1,fCentrality); //all events\r
- Bool_t isSelected = kTRUE;\r
- if(fUseOfflineTrigger)\r
- isSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();\r
- if(isSelected) {\r
- fHistEventStats->Fill(2,fCentrality); //triggered events\r
+ //===========================PID (so far only for electron rejection)===============================// \r
+ if(fElectronRejection) {\r
\r
- if(fUseCentrality) {\r
- //Centrality stuff\r
- AliCentrality *centrality = gESD->GetCentrality();\r
+ // get the electron nsigma\r
+ Double_t nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kElectron));\r
\r
- fCentrality = centrality->GetCentralityPercentile(fCentralityEstimator.Data());\r
+ //Fill QA before the PID\r
+ fHistdEdxVsPTPCbeforePIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());\r
+ fHistNSigmaTPCvsPtbeforePIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),nSigma); \r
+ //end of QA-before pid\r
\r
- // take only events inside centrality class\r
- if(!centrality->IsEventInCentralityClass(fCentralityPercentileMin,\r
- fCentralityPercentileMax,\r
- fCentralityEstimator.Data()))\r
- return;\r
+ // check only for given momentum range\r
+ if( vPt > fElectronRejectionMinPt && vPt < fElectronRejectionMaxPt ){\r
+ \r
+ //look only at electron nsigma\r
+ if(!fElectronOnlyRejection){\r
+ \r
+ //Make the decision based on the n-sigma of electrons\r
+ if(nSigma < fElectronRejectionNSigma) continue;\r
+ }\r
+ else{\r
+ \r
+ Double_t nSigmaPions = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kPion));\r
+ Double_t nSigmaKaons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kKaon));\r
+ Double_t nSigmaProtons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kProton));\r
+ \r
+ //Make the decision based on the n-sigma of electrons exclusively ( = track not in nsigma region for other species)\r
+ if(nSigma < fElectronRejectionNSigma\r
+ && nSigmaPions > fElectronRejectionNSigma\r
+ && nSigmaKaons > fElectronRejectionNSigma\r
+ && nSigmaProtons > fElectronRejectionNSigma ) continue;\r
+ }\r
+ }\r
+ \r
+ //Fill QA after the PID\r
+ fHistdEdxVsPTPCafterPIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());\r
+ fHistNSigmaTPCvsPtafterPIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),nSigma); \r
\r
- // centrality QA (V0M)\r
- fHistV0M->Fill(gESD->GetVZEROData()->GetMTotV0A(), gESD->GetVZEROData()->GetMTotV0C());\r
}\r
+ //===========================end of PID (so far only for electron rejection)===============================//\r
+ \r
+ Float_t dcaXY = aodTrack->DCA(); // this is the DCA from global track (not exactly what is cut on)\r
+ Float_t dcaZ = aodTrack->ZAtDCA(); // this is the DCA from global track (not exactly what is cut on)\r
+ \r
+ // Kinematics cuts from ESD track cuts\r
+ if( vPt < fPtMin || vPt > fPtMax) continue;\r
+ if( vEta < fEtaMin || vEta > fEtaMax) continue;\r
+ \r
+ // Extra DCA cuts (for systematic studies [!= -1])\r
+ if( fDCAxyCut != -1 && fDCAzCut != -1){\r
+ if(TMath::Sqrt((dcaXY*dcaXY)/(fDCAxyCut*fDCAxyCut)+(dcaZ*dcaZ)/(fDCAzCut*fDCAzCut)) > 1 ){\r
+ continue; // 2D cut\r
+ }\r
+ }\r
+ \r
+ // Extra TPC cuts (for systematic studies [!= -1])\r
+ if( fTPCchi2Cut != -1 && aodTrack->Chi2perNDF() > fTPCchi2Cut){\r
+ continue;\r
+ }\r
+ if( fNClustersTPCCut != -1 && aodTrack->GetTPCNcls() < fNClustersTPCCut){\r
+ continue;\r
+ }\r
+\r
+ //Exclude resonances\r
+ if(fExcludeResonancesInMC) {\r
\r
- const AliESDVertex *vertex = gESD->GetPrimaryVertex();\r
- if(vertex) {\r
- if(vertex->GetNContributors() > 0) {\r
- if(vertex->GetZRes() != 0) {\r
- fHistEventStats->Fill(3,fCentrality); //events with a proper vertex\r
- if(TMath::Abs(vertex->GetXv()) < fVxMax) {\r
- if(TMath::Abs(vertex->GetYv()) < fVyMax) {\r
- if(TMath::Abs(vertex->GetZv()) < fVzMax) {\r
- fHistEventStats->Fill(4,fCentrality); //analayzed events\r
- fHistVx->Fill(vertex->GetXv());\r
- fHistVy->Fill(vertex->GetYv());\r
- fHistVz->Fill(vertex->GetZv(),fCentrality);\r
- \r
- //========Get the VZERO event plane========//\r
- Double_t gVZEROEventPlane = -10.0;\r
- Double_t qxTot = 0.0, qyTot = 0.0;\r
- AliEventplane *ep = gESD->GetEventplane();\r
- if(ep) \r
- gVZEROEventPlane = ep->CalculateVZEROEventPlane(gESD,10,2,qxTot,qyTot);\r
- if(gVZEROEventPlane < 0.) gVZEROEventPlane += TMath::Pi();\r
- gReactionPlane = gVZEROEventPlane*TMath::RadToDeg();\r
- fHistEventPlane->Fill(gReactionPlane,fCentrality);\r
- //========Get the VZERO event plane========//\r
-\r
- //Printf("There are %d tracks in this event", gESD->GetNumberOfTracks());\r
- for (Int_t iTracks = 0; iTracks < gESD->GetNumberOfTracks(); iTracks++) {\r
- AliESDtrack* track = dynamic_cast<AliESDtrack *>(gESD->GetTrack(iTracks));\r
- if (!track) {\r
- Printf("ERROR: Could not receive track %d", iTracks);\r
- continue;\r
- } \r
- \r
- Int_t label = TMath::Abs(track->GetLabel());\r
- if(label > mcEvent->GetNumberOfTracks()) continue;\r
- if(label > mcEvent->GetNumberOfPrimaries()) continue;\r
- \r
- AliMCParticle* mcTrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(label));\r
- if(!mcTrack) continue;\r
-\r
- // take only TPC only tracks\r
- track_TPC = new AliESDtrack();\r
- if(!track->FillTPCOnlyTrack(*track_TPC)) continue;\r
- \r
- //ESD track cuts\r
- if(fESDtrackCuts) \r
- if(!fESDtrackCuts->AcceptTrack(track_TPC)) continue;\r
- \r
- // fill QA histograms\r
- Float_t b[2];\r
- Float_t bCov[3];\r
- track_TPC->GetImpactParameters(b,bCov);\r
- if (bCov[0]<=0 || bCov[2]<=0) {\r
- AliDebug(1, "Estimated b resolution lower or equal zero!");\r
- bCov[0]=0; bCov[2]=0;\r
- }\r
- \r
- Int_t nClustersTPC = -1;\r
- nClustersTPC = track_TPC->GetTPCNclsIter1(); // TPC standalone\r
- //nClustersTPC = track->GetTPCclusters(0); // global track\r
- Float_t chi2PerClusterTPC = -1;\r
- if (nClustersTPC!=0) {\r
- chi2PerClusterTPC = track_TPC->GetTPCchi2Iter1()/Float_t(nClustersTPC); // TPC standalone\r
- //chi2PerClusterTPC = track->GetTPCchi2()/Float_t(nClustersTPC); // global track\r
- }\r
- \r
- v_charge = track_TPC->Charge();\r
- v_y = track_TPC->Y();\r
- v_eta = track_TPC->Eta();\r
- v_phi = track_TPC->Phi() * TMath::RadToDeg();\r
- v_E = track_TPC->E();\r
- v_pt = track_TPC->Pt();\r
- track_TPC->PxPyPz(v_p);\r
-\r
- fHistClus->Fill(track_TPC->GetITSclusters(0),nClustersTPC);\r
- fHistDCA->Fill(b[1],b[0]);\r
- fHistChi2->Fill(chi2PerClusterTPC,fCentrality);\r
- fHistPt->Fill(v_pt,fCentrality);\r
- fHistEta->Fill(v_eta,fCentrality);\r
- fHistPhi->Fill(v_phi,fCentrality);\r
- fHistRapidity->Fill(v_y,fCentrality);\r
- if(v_charge > 0) fHistPhiPos->Fill(v_phi,fCentrality);\r
- else if(v_charge < 0) fHistPhiNeg->Fill(v_phi,fCentrality);\r
-\r
- // fill charge vector\r
- chargeVector[0]->push_back(v_charge);\r
- chargeVector[1]->push_back(v_y);\r
- chargeVector[2]->push_back(v_eta);\r
- chargeVector[3]->push_back(v_phi);\r
- chargeVector[4]->push_back(v_p[0]);\r
- chargeVector[5]->push_back(v_p[1]);\r
- chargeVector[6]->push_back(v_p[2]);\r
- chargeVector[7]->push_back(v_pt);\r
- chargeVector[8]->push_back(v_E);\r
-\r
- if(fRunShuffling) {\r
- chargeVectorShuffle[0]->push_back(v_charge);\r
- chargeVectorShuffle[1]->push_back(v_y);\r
- chargeVectorShuffle[2]->push_back(v_eta);\r
- chargeVectorShuffle[3]->push_back(v_phi);\r
- chargeVectorShuffle[4]->push_back(v_p[0]);\r
- chargeVectorShuffle[5]->push_back(v_p[1]);\r
- chargeVectorShuffle[6]->push_back(v_p[2]);\r
- chargeVectorShuffle[7]->push_back(v_pt);\r
- chargeVectorShuffle[8]->push_back(v_E);\r
- }\r
- \r
- delete track_TPC;\r
- gNumberOfAcceptedTracks += 1;\r
- \r
- } //track loop\r
- }//Vz cut\r
- }//Vy cut\r
- }//Vx cut\r
- }//proper vertex resolution\r
- }//proper number of contributors\r
- }//vertex object valid\r
- }//triggered event \r
- }//MC-ESD analysis\r
+ Bool_t kExcludeParticle = kFALSE;\r
\r
- //MC analysis\r
- else if(gAnalysisLevel == "MC") {\r
- AliMCEvent* mcEvent = MCEvent(); \r
- if (!mcEvent) {\r
- Printf("ERROR: mcEvent not available");\r
- return;\r
- }\r
- fHistEventStats->Fill(1,fCentrality); //total events\r
- fHistEventStats->Fill(2,fCentrality); //offline trigger\r
+ Int_t label = TMath::Abs(aodTrack->GetLabel());\r
+ AliAODMCParticle *AODmcTrack = (AliAODMCParticle*) fArrayMC->At(label);\r
+ \r
+ if (AODmcTrack){ \r
+ //if (AODmcTrack->IsPhysicalPrimary()){\r
+ \r
+ Int_t gMotherIndex = AODmcTrack->GetMother();\r
+ if(gMotherIndex != -1) {\r
+ AliAODMCParticle* motherTrack = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(gMotherIndex));\r
+ if(motherTrack) {\r
+ Int_t pdgCodeOfMother = motherTrack->GetPdgCode();\r
+ if(pdgCodeOfMother == 113 // rho0\r
+ || pdgCodeOfMother == 213 || pdgCodeOfMother == -213 // rho+\r
+ // || pdgCodeOfMother == 221 // eta\r
+ // || pdgCodeOfMother == 331 // eta'\r
+ // || pdgCodeOfMother == 223 // omega\r
+ // || pdgCodeOfMother == 333 // phi\r
+ || pdgCodeOfMother == 311 || pdgCodeOfMother == -311 // K0\r
+ // || pdgCodeOfMother == 313 || pdgCodeOfMother == -313 // K0*\r
+ // || pdgCodeOfMother == 323 || pdgCodeOfMother == -323 // K+*\r
+ || pdgCodeOfMother == 3122 || pdgCodeOfMother == -3122 // Lambda\r
+ || pdgCodeOfMother == 111 // pi0 Dalitz\r
+ || pdgCodeOfMother == 22 // photon\r
+ ) {\r
+ kExcludeParticle = kTRUE;\r
+ }\r
+ }\r
+ }\r
+ } \r
+ //Exclude from the analysis decay products of rho0, rho+, eta, eta' and phi\r
+ if(kExcludeParticle) continue;\r
+ }\r
\r
- Double_t gImpactParameter = 0.;\r
- if(fUseCentrality) {\r
- //Get the MC header\r
- AliGenHijingEventHeader* headerH = dynamic_cast<AliGenHijingEventHeader*>(mcEvent->GenEventHeader());\r
- if (headerH) {\r
- //Printf("=====================================================");\r
- //Printf("Reaction plane angle: %lf",headerH->ReactionPlaneAngle());\r
- //Printf("Impact parameter: %lf",headerH->ImpactParameter());\r
- //Printf("=====================================================");\r
- gReactionPlane = headerH->ReactionPlaneAngle();\r
- gImpactParameter = headerH->ImpactParameter();\r
- fCentrality = gImpactParameter;\r
+ //Exclude electrons with PDG\r
+ if(fExcludeElectronsInMC) {\r
+ \r
+ Int_t label = TMath::Abs(aodTrack->GetLabel());\r
+ AliAODMCParticle *AODmcTrack = (AliAODMCParticle*) fArrayMC->At(label);\r
+ \r
+ if (AODmcTrack){ \r
+ if(TMath::Abs(AODmcTrack->GetPdgCode()) == 11) continue;\r
+ }\r
}\r
- fCentrality = gImpactParameter;\r
- fHistEventPlane->Fill(gReactionPlane*TMath::RadToDeg(),fCentrality);\r
+ \r
+ // fill QA histograms\r
+ fHistClus->Fill(aodTrack->GetITSNcls(),aodTrack->GetTPCNcls());\r
+ fHistDCA->Fill(dcaZ,dcaXY);\r
+ fHistChi2->Fill(aodTrack->Chi2perNDF(),gCentrality);\r
+ fHistPt->Fill(vPt,gCentrality);\r
+ fHistEta->Fill(vEta,gCentrality);\r
+ fHistRapidity->Fill(vY,gCentrality);\r
+ if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);\r
+ else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);\r
+ fHistPhi->Fill(vPhi,gCentrality);\r
+ if(vCharge > 0) fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality); \r
+ else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);\r
+ \r
+ //=======================================correction\r
+ Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality); \r
+ //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);\r
+ \r
+ // add the track to the TObjArray\r
+ tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction)); \r
+ }//track loop\r
+ }//MCAODrec\r
+ //==============================================================================================================\r
\r
- // take only events inside centrality class (DIDN'T CHANGE THIS UP TO NOW)\r
- if((fImpactParameterMin > gImpactParameter) || (fImpactParameterMax < gImpactParameter))\r
- return;\r
- }\r
+ else if(gAnalysisLevel == "ESD" || gAnalysisLevel == "MCESD") { // handling of TPC only tracks different in AOD and ESD\r
+\r
+ AliESDtrack *trackTPC = NULL;\r
+ AliMCParticle *mcTrack = NULL;\r
+\r
+ AliMCEvent* mcEvent = NULL;\r
\r
- AliGenEventHeader *header = mcEvent->GenEventHeader();\r
- if(!header) return;\r
+ // for MC ESDs use also MC information\r
+ if(gAnalysisLevel == "MCESD"){\r
+ mcEvent = MCEvent(); \r
+ if (!mcEvent) {\r
+ AliError("mcEvent not available");\r
+ return tracksAccepted;\r
+ }\r
+ }\r
\r
- TArrayF gVertexArray;\r
- header->PrimaryVertex(gVertexArray);\r
- //Printf("Vertex: %lf (x) - %lf (y) - %lf (z)",\r
- //gVertexArray.At(0),\r
- //gVertexArray.At(1),\r
- //gVertexArray.At(2));\r
- fHistEventStats->Fill(3,fCentrality); //events with a proper vertex\r
- if(TMath::Abs(gVertexArray.At(0)) < fVxMax) {\r
- if(TMath::Abs(gVertexArray.At(1)) < fVyMax) {\r
- if(TMath::Abs(gVertexArray.At(2)) < fVzMax) {\r
- fHistEventStats->Fill(4,fCentrality); //analayzed events\r
- fHistVx->Fill(gVertexArray.At(0));\r
- fHistVy->Fill(gVertexArray.At(1));\r
- fHistVz->Fill(gVertexArray.At(2),fCentrality);\r
+ // Loop over tracks in event\r
+ for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {\r
+ AliESDtrack* track = dynamic_cast<AliESDtrack *>(event->GetTrack(iTracks));\r
+ if (!track) {\r
+ AliError(Form("Could not receive track %d", iTracks));\r
+ continue;\r
+ } \r
+\r
+ // for MC ESDs use also MC information --> MC track not used further???\r
+ if(gAnalysisLevel == "MCESD"){\r
+ Int_t label = TMath::Abs(track->GetLabel());\r
+ if(label > mcEvent->GetNumberOfTracks()) continue;\r
+ if(label > mcEvent->GetNumberOfPrimaries()) continue;\r
+ \r
+ mcTrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(label));\r
+ if(!mcTrack) continue;\r
+ }\r
+\r
+ // take only TPC only tracks\r
+ trackTPC = new AliESDtrack();\r
+ if(!track->FillTPCOnlyTrack(*trackTPC)) continue;\r
+ \r
+ //ESD track cuts\r
+ if(fESDtrackCuts) \r
+ if(!fESDtrackCuts->AcceptTrack(trackTPC)) continue;\r
+ \r
+ // fill QA histograms\r
+ Float_t b[2];\r
+ Float_t bCov[3];\r
+ trackTPC->GetImpactParameters(b,bCov);\r
+ if (bCov[0]<=0 || bCov[2]<=0) {\r
+ AliDebug(1, "Estimated b resolution lower or equal zero!");\r
+ bCov[0]=0; bCov[2]=0;\r
+ }\r
+ \r
+ Int_t nClustersTPC = -1;\r
+ nClustersTPC = trackTPC->GetTPCNclsIter1(); // TPC standalone\r
+ //nClustersTPC = track->GetTPCclusters(0); // global track\r
+ Float_t chi2PerClusterTPC = -1;\r
+ if (nClustersTPC!=0) {\r
+ chi2PerClusterTPC = trackTPC->GetTPCchi2Iter1()/Float_t(nClustersTPC); // TPC standalone\r
+ //chi2PerClusterTPC = track->GetTPCchi2()/Float_t(nClustersTPC); // global track\r
+ }\r
+ \r
+ //===========================PID===============================// \r
+ if(fUsePID) {\r
+ Double_t prob[AliPID::kSPECIES]={0.};\r
+ Double_t probTPC[AliPID::kSPECIES]={0.};\r
+ Double_t probTOF[AliPID::kSPECIES]={0.};\r
+ Double_t probTPCTOF[AliPID::kSPECIES]={0.};\r
+ \r
+ Double_t nSigma = 0.;\r
+ UInt_t detUsedTPC = 0;\r
+ UInt_t detUsedTOF = 0;\r
+ UInt_t detUsedTPCTOF = 0;\r
+ \r
+ //Decide what detector configuration we want to use\r
+ switch(fPidDetectorConfig) {\r
+ case kTPCpid:\r
+ fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);\r
+ nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)fParticleOfInterest));\r
+ detUsedTPC = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTPC);\r
+ for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)\r
+ prob[iSpecies] = probTPC[iSpecies];\r
+ break;\r
+ case kTOFpid:\r
+ fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF);\r
+ nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)fParticleOfInterest));\r
+ detUsedTOF = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTOF);\r
+ for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)\r
+ prob[iSpecies] = probTOF[iSpecies];\r
+ break;\r
+ case kTPCTOF:\r
+ fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC);\r
+ detUsedTPCTOF = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTPCTOF);\r
+ for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)\r
+ prob[iSpecies] = probTPCTOF[iSpecies];\r
+ break;\r
+ default:\r
+ break;\r
+ }//end switch: define detector mask\r
+ \r
+ //Filling the PID QA\r
+ Double_t tofTime = -999., length = 999., tof = -999.;\r
+ Double_t c = TMath::C()*1.E-9;// m/ns\r
+ Double_t beta = -999.;\r
+ Double_t nSigmaTOFForParticleOfInterest = -999.;\r
+ if ( (track->IsOn(AliESDtrack::kTOFin)) &&\r
+ (track->IsOn(AliESDtrack::kTIME)) ) { \r
+ tofTime = track->GetTOFsignal();//in ps\r
+ length = track->GetIntegratedLength();\r
+ tof = tofTime*1E-3; // ns \r
\r
- Printf("There are %d tracks in this event", mcEvent->GetNumberOfPrimaries());\r
- for (Int_t iTracks = 0; iTracks < mcEvent->GetNumberOfPrimaries(); iTracks++) {\r
- AliMCParticle* track = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(iTracks));\r
- if (!track) {\r
- Printf("ERROR: Could not receive particle %d", iTracks);\r
- continue;\r
- }\r
- \r
- //exclude non stable particles\r
- if(!mcEvent->IsPhysicalPrimary(iTracks)) continue;\r
+ if (tof <= 0) {\r
+ //Printf("WARNING: track with negative TOF time found! Skipping this track for PID checks\n");\r
+ continue;\r
+ }\r
+ if (length <= 0){\r
+ //printf("WARNING: track with negative length found!Skipping this track for PID checks\n");\r
+ continue;\r
+ }\r
+ \r
+ length = length*0.01; // in meters\r
+ tof = tof*c;\r
+ beta = length/tof;\r
+ \r
+ nSigmaTOFForParticleOfInterest = fPIDResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)fParticleOfInterest);\r
+ fHistBetavsPTOFbeforePID ->Fill(track->P()*track->Charge(),beta);\r
+ fHistProbTOFvsPtbeforePID ->Fill(track->Pt(),probTOF[fParticleOfInterest]);\r
+ fHistNSigmaTOFvsPtbeforePID ->Fill(track->Pt(),nSigmaTOFForParticleOfInterest);\r
+ }//TOF signal \r
+ \r
+ \r
+ Double_t nSigmaTPCForParticleOfInterest = fPIDResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)fParticleOfInterest);\r
+ fHistdEdxVsPTPCbeforePID -> Fill(track->P()*track->Charge(),track->GetTPCsignal());\r
+ fHistProbTPCvsPtbeforePID -> Fill(track->Pt(),probTPC[fParticleOfInterest]); \r
+ fHistNSigmaTPCvsPtbeforePID -> Fill(track->Pt(),nSigmaTPCForParticleOfInterest); \r
+ fHistProbTPCTOFvsPtbeforePID -> Fill(track->Pt(),probTPCTOF[fParticleOfInterest]);\r
+ //end of QA-before pid\r
+ \r
+ if ((detUsedTPC != 0)||(detUsedTOF != 0)||(detUsedTPCTOF != 0)) {\r
+ //Make the decision based on the n-sigma\r
+ if(fUsePIDnSigma) {\r
+ if(nSigma > fPIDNSigma) continue;}\r
+ \r
+ //Make the decision based on the bayesian\r
+ else if(fUsePIDPropabilities) {\r
+ if(fParticleOfInterest != TMath::LocMax(AliPID::kSPECIES,prob)) continue;\r
+ if (prob[fParticleOfInterest] < fMinAcceptedPIDProbability) continue; \r
+ }\r
+ \r
+ //Fill QA after the PID\r
+ fHistBetavsPTOFafterPID ->Fill(track->P()*track->Charge(),beta);\r
+ fHistProbTOFvsPtafterPID ->Fill(track->Pt(),probTOF[fParticleOfInterest]);\r
+ fHistNSigmaTOFvsPtafterPID ->Fill(track->Pt(),nSigmaTOFForParticleOfInterest);\r
+ \r
+ fHistdEdxVsPTPCafterPID -> Fill(track->P()*track->Charge(),track->GetTPCsignal());\r
+ fHistProbTPCvsPtafterPID -> Fill(track->Pt(),probTPC[fParticleOfInterest]); \r
+ fHistProbTPCTOFvsPtafterPID -> Fill(track->Pt(),probTPCTOF[fParticleOfInterest]);\r
+ fHistNSigmaTPCvsPtafterPID -> Fill(track->Pt(),nSigmaTPCForParticleOfInterest); \r
+ }\r
+ }\r
+ //===========================PID===============================//\r
+ vCharge = trackTPC->Charge();\r
+ vY = trackTPC->Y();\r
+ vEta = trackTPC->Eta();\r
+ vPhi = trackTPC->Phi();// * TMath::RadToDeg();\r
+ vPt = trackTPC->Pt();\r
+ fHistClus->Fill(trackTPC->GetITSclusters(0),nClustersTPC);\r
+ fHistDCA->Fill(b[1],b[0]);\r
+ fHistChi2->Fill(chi2PerClusterTPC,gCentrality);\r
+ fHistPt->Fill(vPt,gCentrality);\r
+ fHistEta->Fill(vEta,gCentrality);\r
+ fHistPhi->Fill(vPhi,gCentrality);\r
+ if(vCharge > 0) fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);\r
+ else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);\r
+ fHistRapidity->Fill(vY,gCentrality);\r
+ if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);\r
+ else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);\r
+ \r
+ //=======================================correction\r
+ Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality); \r
+ //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);\r
\r
- v_eta = track->Eta();\r
- v_pt = track->Pt();\r
- v_y = track->Y();\r
+ // add the track to the TObjArray\r
+ tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction)); \r
\r
- if( v_pt < fPtMin || v_pt > fPtMax) \r
- continue;\r
- if (!fUsePID) {\r
- if( v_eta < fEtaMin || v_eta > fEtaMax) continue;\r
- }\r
- else if (fUsePID){\r
- if( v_y < fEtaMin || v_y > fEtaMax) continue;\r
- }\r
+ delete trackTPC;\r
+ }//track loop\r
+ }// ESD analysis\r
\r
- //analyze one set of particles\r
- if(fUseMCPdgCode) {\r
- TParticle *particle = track->Particle();\r
- if(!particle) continue;\r
- \r
- Int_t gPdgCode = particle->GetPdgCode();\r
- if(TMath::Abs(fPDGCodeToBeAnalyzed) != TMath::Abs(gPdgCode)) \r
- continue;\r
- }\r
- \r
- //Use the acceptance parameterization\r
- if(fAcceptanceParameterization) {\r
- Double_t gRandomNumber = gRandom->Rndm();\r
- if(gRandomNumber > fAcceptanceParameterization->Eval(track->Pt())) \r
- continue;\r
- }\r
- \r
- //Exclude resonances\r
- if(fExcludeResonancesInMC) {\r
- TParticle *particle = track->Particle();\r
- if(!particle) continue;\r
- \r
- Bool_t kExcludeParticle = kFALSE;\r
- Int_t gMotherIndex = particle->GetFirstMother();\r
- if(gMotherIndex != -1) {\r
- AliMCParticle* motherTrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(gMotherIndex));\r
- if(motherTrack) {\r
- TParticle *motherParticle = motherTrack->Particle();\r
- if(motherParticle) {\r
- Int_t pdgCodeOfMother = motherParticle->GetPdgCode();\r
- //if((pdgCodeOfMother == 113)||(pdgCodeOfMother == 213)||(pdgCodeOfMother == 221)||(pdgCodeOfMother == 223)||(pdgCodeOfMother == 331)||(pdgCodeOfMother == 333)) {\r
- if(pdgCodeOfMother == 113) {\r
- kExcludeParticle = kTRUE;\r
- }\r
- }\r
+ else if(gAnalysisLevel == "MC"){\r
+ if(!event) {\r
+ AliError("mcEvent not available");\r
+ return 0x0;\r
+ }\r
+\r
+ AliMCEvent *gMCEvent = dynamic_cast<AliMCEvent*>(event);\r
+ if(gMCEvent) {\r
+ // Loop over tracks in event\r
+ for (Int_t iTracks = 0; iTracks < gMCEvent->GetNumberOfPrimaries(); iTracks++) {\r
+ AliMCParticle* track = dynamic_cast<AliMCParticle *>(gMCEvent->GetTrack(iTracks));\r
+ if (!track) {\r
+ AliError(Form("Could not receive particle %d", iTracks));\r
+ continue;\r
+ }\r
+ \r
+ //exclude non stable particles\r
+ if(!(gMCEvent->IsPhysicalPrimary(iTracks))) continue;\r
+ \r
+ vCharge = track->Charge();\r
+ vEta = track->Eta();\r
+ vPt = track->Pt();\r
+ vY = track->Y();\r
+ \r
+ if( vPt < fPtMin || vPt > fPtMax) \r
+ continue;\r
+ if (!fUsePID) {\r
+ if( vEta < fEtaMin || vEta > fEtaMax) continue;\r
+ }\r
+ else if (fUsePID){\r
+ if( vY < fEtaMin || vY > fEtaMax) continue;\r
+ }\r
+\r
+ // Remove neutral tracks\r
+ if( vCharge == 0 ) continue;\r
+ \r
+ //analyze one set of particles\r
+ if(fUseMCPdgCode) {\r
+ TParticle *particle = track->Particle();\r
+ if(!particle) continue;\r
+ \r
+ Int_t gPdgCode = particle->GetPdgCode();\r
+ if(TMath::Abs(fPDGCodeToBeAnalyzed) != TMath::Abs(gPdgCode)) \r
+ continue;\r
+ }\r
+ \r
+ //Use the acceptance parameterization\r
+ if(fAcceptanceParameterization) {\r
+ Double_t gRandomNumber = gRandom->Rndm();\r
+ if(gRandomNumber > fAcceptanceParameterization->Eval(track->Pt())) \r
+ continue;\r
+ }\r
+ \r
+ //Exclude resonances\r
+ if(fExcludeResonancesInMC) {\r
+ TParticle *particle = track->Particle();\r
+ if(!particle) continue;\r
+ \r
+ Bool_t kExcludeParticle = kFALSE;\r
+ Int_t gMotherIndex = particle->GetFirstMother();\r
+ if(gMotherIndex != -1) {\r
+ AliMCParticle* motherTrack = dynamic_cast<AliMCParticle *>(event->GetTrack(gMotherIndex));\r
+ if(motherTrack) {\r
+ TParticle *motherParticle = motherTrack->Particle();\r
+ if(motherParticle) {\r
+ Int_t pdgCodeOfMother = motherParticle->GetPdgCode();\r
+ //if((pdgCodeOfMother == 113)||(pdgCodeOfMother == 213)||(pdgCodeOfMother == 221)||(pdgCodeOfMother == 223)||(pdgCodeOfMother == 331)||(pdgCodeOfMother == 333)) {\r
+ if(pdgCodeOfMother == 113 // rho0\r
+ || pdgCodeOfMother == 213 || pdgCodeOfMother == -213 // rho+\r
+ // || pdgCodeOfMother == 221 // eta\r
+ // || pdgCodeOfMother == 331 // eta'\r
+ // || pdgCodeOfMother == 223 // omega\r
+ // || pdgCodeOfMother == 333 // phi\r
+ || pdgCodeOfMother == 311 || pdgCodeOfMother == -311 // K0\r
+ // || pdgCodeOfMother == 313 || pdgCodeOfMother == -313 // K0*\r
+ // || pdgCodeOfMother == 323 || pdgCodeOfMother == -323 // K+*\r
+ || pdgCodeOfMother == 3122 || pdgCodeOfMother == -3122 // Lambda\r
+ || pdgCodeOfMother == 111 // pi0 Dalitz\r
+ ) {\r
+ kExcludeParticle = kTRUE;\r
}\r
}\r
- \r
- //Exclude from the analysis decay products of rho0, rho+, eta, eta' and phi\r
- if(kExcludeParticle) continue;\r
}\r
+ }\r
+ \r
+ //Exclude from the analysis decay products of rho0, rho+, eta, eta' and phi\r
+ if(kExcludeParticle) continue;\r
+ }\r
\r
- v_charge = track->Charge();\r
- v_phi = track->Phi();\r
- v_E = track->E();\r
- track->PxPyPz(v_p);\r
- //Printf("phi (before): %lf",v_phi);\r
-\r
- fHistPt->Fill(v_pt,fCentrality);\r
- fHistEta->Fill(v_eta,fCentrality);\r
- fHistPhi->Fill(v_phi*TMath::RadToDeg(),fCentrality);\r
- fHistRapidity->Fill(v_y,fCentrality);\r
- if(v_charge > 0) fHistPhiPos->Fill(v_phi*TMath::RadToDeg(),fCentrality);\r
- else if(v_charge < 0) fHistPhiNeg->Fill(v_phi*TMath::RadToDeg(),fCentrality);\r
-\r
- //Flow after burner\r
- if(fUseFlowAfterBurner) {\r
- Double_t precisionPhi = 0.001;\r
- Int_t maxNumberOfIterations = 100;\r
-\r
- Double_t phi0 = v_phi;\r
- Double_t gV2 = fDifferentialV2->Eval(v_pt);\r
-\r
- for (Int_t j = 0; j < maxNumberOfIterations; j++) {\r
- Double_t phiprev = v_phi;\r
- Double_t fl = v_phi - phi0 + gV2*TMath::Sin(2.*(v_phi - gReactionPlane));\r
- Double_t fp = 1.0 + 2.0*gV2*TMath::Cos(2.*(v_phi - gReactionPlane)); \r
- v_phi -= fl/fp;\r
- if (TMath::AreEqualAbs(phiprev,v_phi,precisionPhi)) break;\r
- }\r
- //Printf("phi (after): %lf\n",v_phi);\r
- Double_t v_DeltaphiBefore = phi0 - gReactionPlane;\r
- if(v_DeltaphiBefore < 0) v_DeltaphiBefore += 2*TMath::Pi();\r
- fHistPhiBefore->Fill(v_DeltaphiBefore,fCentrality);\r
-\r
- Double_t v_DeltaphiAfter = v_phi - gReactionPlane;\r
- if(v_DeltaphiAfter < 0) v_DeltaphiAfter += 2*TMath::Pi();\r
- fHistPhiAfter->Fill(v_DeltaphiAfter,fCentrality);\r
- }\r
- \r
- v_phi *= TMath::RadToDeg();\r
-\r
- // fill charge vector\r
- chargeVector[0]->push_back(v_charge);\r
- chargeVector[1]->push_back(v_y);\r
- chargeVector[2]->push_back(v_eta);\r
- chargeVector[3]->push_back(v_phi);\r
- chargeVector[4]->push_back(v_p[0]);\r
- chargeVector[5]->push_back(v_p[1]);\r
- chargeVector[6]->push_back(v_p[2]);\r
- chargeVector[7]->push_back(v_pt);\r
- chargeVector[8]->push_back(v_E);\r
- \r
- if(fRunShuffling) {\r
- chargeVectorShuffle[0]->push_back(v_charge);\r
- chargeVectorShuffle[1]->push_back(v_y);\r
- chargeVectorShuffle[2]->push_back(v_eta);\r
- chargeVectorShuffle[3]->push_back(v_phi);\r
- chargeVectorShuffle[4]->push_back(v_p[0]);\r
- chargeVectorShuffle[5]->push_back(v_p[1]);\r
- chargeVectorShuffle[6]->push_back(v_p[2]);\r
- chargeVectorShuffle[7]->push_back(v_pt);\r
- chargeVectorShuffle[8]->push_back(v_E);\r
- }\r
- gNumberOfAcceptedTracks += 1;\r
- \r
- } //track loop\r
- gReactionPlane *= TMath::RadToDeg();\r
- }//Vz cut\r
- }//Vy cut\r
- }//Vx cut\r
- }//MC analysis\r
+ //Exclude electrons with PDG\r
+ if(fExcludeElectronsInMC) {\r
+ \r
+ TParticle *particle = track->Particle();\r
+ \r
+ if (particle){ \r
+ if(TMath::Abs(particle->GetPdgCode()) == 11) continue;\r
+ }\r
+ }\r
+ \r
+ vPhi = track->Phi();\r
+ //Printf("phi (before): %lf",vPhi);\r
+ \r
+ fHistPt->Fill(vPt,gCentrality);\r
+ fHistEta->Fill(vEta,gCentrality);\r
+ fHistPhi->Fill(vPhi,gCentrality);\r
+ if(vCharge > 0) fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);\r
+ else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);\r
+ //fHistPhi->Fill(vPhi*TMath::RadToDeg(),gCentrality);\r
+ fHistRapidity->Fill(vY,gCentrality);\r
+ //if(vCharge > 0) fHistPhiPos->Fill(vPhi*TMath::RadToDeg(),gCentrality);\r
+ //else if(vCharge < 0) fHistPhiNeg->Fill(vPhi*TMath::RadToDeg(),gCentrality);\r
+ if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);\r
+ else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);\r
+ \r
+ //Flow after burner\r
+ if(fUseFlowAfterBurner) {\r
+ Double_t precisionPhi = 0.001;\r
+ Int_t maxNumberOfIterations = 100;\r
+ \r
+ Double_t phi0 = vPhi;\r
+ Double_t gV2 = fDifferentialV2->Eval(vPt);\r
+ \r
+ for (Int_t j = 0; j < maxNumberOfIterations; j++) {\r
+ Double_t phiprev = vPhi;\r
+ Double_t fl = vPhi - phi0 + gV2*TMath::Sin(2.*(vPhi - gReactionPlane*TMath::DegToRad()));\r
+ Double_t fp = 1.0 + 2.0*gV2*TMath::Cos(2.*(vPhi - gReactionPlane*TMath::DegToRad())); \r
+ vPhi -= fl/fp;\r
+ if (TMath::AreEqualAbs(phiprev,vPhi,precisionPhi)) break;\r
+ }\r
+ //Printf("phi (after): %lf\n",vPhi);\r
+ Double_t vDeltaphiBefore = phi0 - gReactionPlane*TMath::DegToRad();\r
+ if(vDeltaphiBefore < 0) vDeltaphiBefore += 2*TMath::Pi();\r
+ fHistPhiBefore->Fill(vDeltaphiBefore,gCentrality);\r
+ \r
+ Double_t vDeltaphiAfter = vPhi - gReactionPlane*TMath::DegToRad();\r
+ if(vDeltaphiAfter < 0) vDeltaphiAfter += 2*TMath::Pi();\r
+ fHistPhiAfter->Fill(vDeltaphiAfter,gCentrality);\r
+ \r
+ }\r
+ \r
+ //vPhi *= TMath::RadToDeg();\r
+ \r
+ //=======================================correction\r
+ Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality); \r
+ //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);\r
+\r
+ tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction)); \r
+ } //track loop\r
+ }//MC event object\r
+ }//MC\r
\r
- //multiplicity cut (used in pp)\r
- if(fUseMultiplicity) {\r
- if((gNumberOfAcceptedTracks < fNumberOfAcceptedTracksMin)||(gNumberOfAcceptedTracks > fNumberOfAcceptedTracksMax))\r
- return;\r
- }\r
- fHistNumberOfAcceptedTracks->Fill(gNumberOfAcceptedTracks,fCentrality);\r
+ return tracksAccepted; \r
+}\r
+\r
+//________________________________________________________________________\r
+TObjArray* AliAnalysisTaskBFPsi::GetShuffledTracks(TObjArray *tracks, Double_t gCentrality){\r
+ // Clones TObjArray and returns it with tracks after shuffling the charges\r
+\r
+ TObjArray* tracksShuffled = new TObjArray;\r
+ tracksShuffled->SetOwner(kTRUE);\r
+\r
+ vector<Short_t> *chargeVector = new vector<Short_t>; //original charge of accepted tracks \r
+\r
+ for (Int_t i=0; i<tracks->GetEntriesFast(); i++)\r
+ {\r
+ AliVParticle* track = (AliVParticle*) tracks->At(i);\r
+ chargeVector->push_back(track->Charge());\r
+ } \r
+ \r
+ random_shuffle(chargeVector->begin(), chargeVector->end());\r
\r
- // calculate balance function\r
- if(fUseMultiplicity) \r
- //fBalance->CalculateBalance(gNumberOfAcceptedTracks,gReactionPlane,chargeVector);\r
- fBalance->CalculateBalance(gReactionPlane,chargeVector);\r
- else \r
- fBalance->CalculateBalance(gReactionPlane,chargeVector);\r
+ for(Int_t i = 0; i < tracks->GetEntriesFast(); i++){\r
+ AliVParticle* track = (AliVParticle*) tracks->At(i);\r
+ //==============================correction\r
+ Double_t correction = GetTrackbyTrackCorrectionMatrix(track->Eta(), track->Phi(),track->Pt(), chargeVector->at(i), gCentrality);\r
+ //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);\r
+ tracksShuffled->Add(new AliBFBasicParticle(track->Eta(), track->Phi(), track->Pt(),chargeVector->at(i), correction));\r
+ }\r
\r
- if(fRunShuffling) {\r
- // shuffle charges\r
- random_shuffle( chargeVectorShuffle[0]->begin(), chargeVectorShuffle[0]->end() );\r
- if(fUseMultiplicity) \r
- fShuffledBalance->CalculateBalance(gReactionPlane,chargeVectorShuffle);\r
- else \r
- fShuffledBalance->CalculateBalance(gReactionPlane,chargeVectorShuffle);\r
+ delete chargeVector;\r
+ \r
+ return tracksShuffled;\r
+}\r
+\r
+//________________________________________________________________________\r
+void AliAnalysisTaskBFPsi::SetVZEROCalibrationFile(const char* filename,\r
+ const char* lhcPeriod) {\r
+ //Function to setup the VZERO gain equalization\r
+ //============Get the equilization map============//\r
+ TFile *calibrationFile = TFile::Open(filename);\r
+ if((!calibrationFile)||(!calibrationFile->IsOpen())) {\r
+ Printf("No calibration file found!!!");\r
+ return;\r
+ }\r
+\r
+ TList *list = dynamic_cast<TList *>(calibrationFile->Get(lhcPeriod));\r
+ if(!list) {\r
+ Printf("Calibration TList not found!!!");\r
+ return;\r
+ }\r
+\r
+ fHistVZEROAGainEqualizationMap = dynamic_cast<TH1F *>(list->FindObject("gHistVZEROAGainEqualizationMap"));\r
+ if(!fHistVZEROAGainEqualizationMap) {\r
+ Printf("VZERO-A calibration object not found!!!");\r
+ return;\r
+ }\r
+ fHistVZEROCGainEqualizationMap = dynamic_cast<TH1F *>(list->FindObject("gHistVZEROCGainEqualizationMap"));\r
+ if(!fHistVZEROCGainEqualizationMap) {\r
+ Printf("VZERO-C calibration object not found!!!");\r
+ return;\r
}\r
-} \r
+\r
+ fHistVZEROChannelGainEqualizationMap = dynamic_cast<TH2F *>(list->FindObject("gHistVZEROChannelGainEqualizationMap"));\r
+ if(!fHistVZEROChannelGainEqualizationMap) {\r
+ Printf("VZERO channel calibration object not found!!!");\r
+ return;\r
+ }\r
+}\r
+\r
+//________________________________________________________________________\r
+Double_t AliAnalysisTaskBFPsi::GetChannelEqualizationFactor(Int_t run, \r
+ Int_t channel) {\r
+ //\r
+ if(!fHistVZEROAGainEqualizationMap) return 1.0;\r
+\r
+ for(Int_t iBinX = 1; iBinX <= fHistVZEROChannelGainEqualizationMap->GetNbinsX(); iBinX++) {\r
+ Int_t gRunNumber = atoi(fHistVZEROChannelGainEqualizationMap->GetXaxis()->GetBinLabel(iBinX));\r
+ if(gRunNumber == run)\r
+ return fHistVZEROChannelGainEqualizationMap->GetBinContent(iBinX,channel+1);\r
+ }\r
+\r
+ return 1.0;\r
+}\r
+\r
+//________________________________________________________________________\r
+Double_t AliAnalysisTaskBFPsi::GetEqualizationFactor(Int_t run, \r
+ const char* side) {\r
+ //\r
+ if(!fHistVZEROAGainEqualizationMap) return 1.0;\r
+\r
+ TString gVZEROSide = side;\r
+ for(Int_t iBinX = 1; iBinX < fHistVZEROAGainEqualizationMap->GetNbinsX(); iBinX++) {\r
+ Int_t gRunNumber = atoi(fHistVZEROAGainEqualizationMap->GetXaxis()->GetBinLabel(iBinX));\r
+ //cout<<"Looking for run "<<run<<" - current run: "<<gRunNumber<<endl;\r
+ if(gRunNumber == run) {\r
+ if(gVZEROSide == "A") \r
+ return fHistVZEROAGainEqualizationMap->GetBinContent(iBinX);\r
+ else if(gVZEROSide == "C") \r
+ return fHistVZEROCGainEqualizationMap->GetBinContent(iBinX);\r
+ }\r
+ }\r
+\r
+ return 1.0;\r
+}\r
\r
//________________________________________________________________________\r
void AliAnalysisTaskBFPsi::FinishTaskOutput(){\r
//Printf("END BF");\r
\r
if (!fBalance) {\r
- Printf("ERROR: fBalance not available");\r
+ AliError("fBalance not available");\r
return;\r
} \r
if(fRunShuffling) {\r
if (!fShuffledBalance) {\r
- Printf("ERROR: fShuffledBalance not available");\r
+ AliError("fShuffledBalance not available");\r
return;\r
}\r
}\r
// not implemented ...\r
\r
}\r
+\r
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff