-#include "TChain.h"\r
-#include "TList.h"\r
-#include "TCanvas.h"\r
-#include "TLorentzVector.h"\r
-#include "TGraphErrors.h"\r
-#include "TH1F.h"\r
-#include "TH2F.h"\r
-#include "TArrayF.h"\r
-#include "TF1.h"\r
-#include "TRandom.h"\r
-\r
-#include "AliAnalysisTaskSE.h"\r
-#include "AliAnalysisManager.h"\r
-\r
-#include "AliESDVertex.h"\r
-#include "AliESDEvent.h"\r
-#include "AliESDInputHandler.h"\r
-#include "AliAODEvent.h"\r
-#include "AliAODTrack.h"\r
-#include "AliAODInputHandler.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
-\r
-#include "TH2D.h" \r
-#include "AliPID.h" \r
-#include "AliPIDResponse.h" \r
-#include "AliPIDCombined.h" \r
-\r
-#include "AliAnalysisTaskBF.h"\r
-#include "AliBalance.h"\r
-\r
-\r
-// Analysis task for the BF code\r
-// Authors: Panos.Christakoglou@nikhef.nl\r
-\r
-ClassImp(AliAnalysisTaskBF)\r
-\r
-//________________________________________________________________________\r
-AliAnalysisTaskBF::AliAnalysisTaskBF(const char *name) \r
-: AliAnalysisTaskSE(name), \r
- fBalance(0),\r
- fRunShuffling(kFALSE),\r
- fShuffledBalance(0),\r
- fList(0),\r
- fListBF(0),\r
- fListBFS(0),\r
- fHistListPIDQA(0),\r
- fHistEventStats(0),\r
- fHistCentStats(0),\r
- fHistTriggerStats(0),\r
- fHistTrackStats(0),\r
- fHistVx(0),\r
- fHistVy(0),\r
- fHistVz(0),\r
- fHistClus(0),\r
- fHistDCA(0),\r
- fHistChi2(0),\r
- fHistPt(0),\r
- fHistEta(0),\r
- fHistRapidity(0),\r
- fHistPhi(0),\r
- fHistPhiBefore(0),\r
- fHistPhiAfter(0),\r
- fHistPhiPos(0),\r
- fHistPhiNeg(0),\r
- fHistV0M(0),\r
- fHistRefTracks(0),\r
- fHistdEdxVsPTPCbeforePID(NULL),\r
- fHistBetavsPTOFbeforePID(NULL), \r
- fHistProbTPCvsPtbeforePID(NULL), \r
- fHistProbTOFvsPtbeforePID(NULL), \r
- fHistProbTPCTOFvsPtbeforePID(NULL),\r
- fHistNSigmaTPCvsPtbeforePID(NULL), \r
- fHistNSigmaTOFvsPtbeforePID(NULL), \r
- fHistdEdxVsPTPCafterPID(NULL),\r
- fHistBetavsPTOFafterPID(NULL), \r
- fHistProbTPCvsPtafterPID(NULL), \r
- fHistProbTOFvsPtafterPID(NULL), \r
- fHistProbTPCTOFvsPtafterPID(NULL),\r
- fHistNSigmaTPCvsPtafterPID(NULL), \r
- fHistNSigmaTOFvsPtafterPID(NULL), \r
- fPIDResponse(0x0),\r
- fPIDCombined(0x0),\r
- fParticleOfInterest(kPion),\r
- fPidDetectorConfig(kTPCTOF),\r
- fUsePID(kFALSE),\r
- fUsePIDnSigma(kTRUE),\r
- fUsePIDPropabilities(kFALSE), \r
- fPIDNSigma(3.),\r
- fMinAcceptedPIDProbability(0.8),\r
- fESDtrackCuts(0),\r
- fCentralityEstimator("V0M"),\r
- fUseCentrality(kFALSE),\r
- fCentralityPercentileMin(0.), \r
- fCentralityPercentileMax(5.),\r
- fImpactParameterMin(0.),\r
- fImpactParameterMax(20.),\r
- fUseMultiplicity(kFALSE),\r
- fNumberOfAcceptedTracksMin(0),\r
- fNumberOfAcceptedTracksMax(10000),\r
- fHistNumberOfAcceptedTracks(0),\r
- fUseOfflineTrigger(kFALSE),\r
- fVxMax(0.3),\r
- fVyMax(0.3),\r
- fVzMax(10.),\r
- nAODtrackCutBit(128),\r
- fPtMin(0.3),\r
- fPtMax(1.5),\r
- fEtaMin(-0.8),\r
- fEtaMax(-0.8),\r
- fDCAxyCut(-1),\r
- fDCAzCut(-1),\r
- fTPCchi2Cut(-1),\r
- fNClustersTPCCut(-1),\r
- fAcceptanceParameterization(0),\r
- fDifferentialV2(0),\r
- fUseFlowAfterBurner(kFALSE),\r
- fExcludeResonancesInMC(kFALSE),\r
- fUseMCPdgCode(kFALSE),\r
- fPDGCodeToBeAnalyzed(-1) {\r
- // Constructor\r
- // Define input and output slots here\r
- // Input slot #0 works with a TChain\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
-}\r
-\r
-//________________________________________________________________________\r
-AliAnalysisTaskBF::~AliAnalysisTaskBF() {\r
-\r
- // delete fBalance; \r
- // delete fShuffledBalance; \r
- // delete fList;\r
- // delete fListBF; \r
- // delete fListBFS;\r
-\r
- // delete fHistEventStats; \r
- // delete fHistTrackStats; \r
- // delete fHistVx; \r
- // delete fHistVy; \r
- // delete fHistVz; \r
-\r
- // delete fHistClus;\r
- // delete fHistDCA;\r
- // delete fHistChi2;\r
- // delete fHistPt;\r
- // delete fHistEta;\r
- // delete fHistPhi;\r
- // delete fHistV0M;\r
-}\r
-\r
-//________________________________________________________________________\r
-void AliAnalysisTaskBF::UserCreateOutputObjects() {\r
- // Create histograms\r
- // Called once\r
- if(!fBalance) {\r
- fBalance = new AliBalance();\r
- fBalance->SetAnalysisLevel("ESD");\r
- //fBalance->SetNumberOfBins(-1,16);\r
- fBalance->SetInterval(-1,-0.8,0.8,16,0.,1.6);\r
- }\r
- if(fRunShuffling) {\r
- if(!fShuffledBalance) {\r
- fShuffledBalance = new AliBalance();\r
- fShuffledBalance->SetAnalysisLevel("ESD");\r
- //fShuffledBalance->SetNumberOfBins(-1,16);\r
- fShuffledBalance->SetInterval(-1,-0.8,0.8,16,0.,1.6);\r
- }\r
- }\r
-\r
- //QA list\r
- fList = new TList();\r
- fList->SetName("listQA");\r
- fList->SetOwner();\r
-\r
- //Balance Function list\r
- fListBF = new TList();\r
- fListBF->SetName("listBF");\r
- fListBF->SetOwner();\r
-\r
- if(fRunShuffling) {\r
- fListBFS = new TList();\r
- fListBFS->SetName("listBFShuffled");\r
- fListBFS->SetOwner();\r
- }\r
-\r
- //PID QA list\r
- if(fUsePID) {\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
- fHistEventStats = new TH1F("fHistEventStats",\r
- "Event statistics;;N_{events}",\r
- 4,0.5,4.5);\r
- for(Int_t i = 1; i <= 4; 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
- 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
- fHistCentStats->GetXaxis()->SetBinLabel(i,gCentName[i-1].Data());\r
- fList->Add(fHistCentStats);\r
-\r
- fHistTriggerStats = new TH1F("fHistTriggerStats","Trigger statistics;TriggerBit;N_{events}",130,0,130);\r
- fList->Add(fHistTriggerStats);\r
-\r
- fHistTrackStats = new TH1F("fHistTrackStats","Event statistics;TrackFilterBit;N_{events}",130,0,130);\r
- fList->Add(fHistTrackStats);\r
-\r
- fHistNumberOfAcceptedTracks = new TH1F("fHistNumberOfAcceptedTracks",";N_{acc.};Entries",4001,-0.5,4000.5);\r
- fList->Add(fHistNumberOfAcceptedTracks);\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
- fHistVy = new TH1F("fHistVy","Primary vertex distribution - y coordinate;V_{y} (cm);Entries",100,-0.5,0.5);\r
- fList->Add(fHistVy);\r
- fHistVz = new TH1F("fHistVz","Primary vertex distribution - z coordinate;V_{z} (cm);Entries",100,-20.,20.);\r
- fList->Add(fHistVz);\r
-\r
- // QA histograms\r
- fHistClus = new TH2F("fHistClus","# Cluster (TPC vs. ITS)",10,0,10,200,0,200);\r
- fList->Add(fHistClus);\r
- fHistChi2 = new TH1F("fHistChi2","Chi2/NDF distribution",200,0,10);\r
- fList->Add(fHistChi2);\r
- fHistDCA = new TH2F("fHistDCA","DCA (xy vs. z)",400,-5,5,400,-5,5); \r
- fList->Add(fHistDCA);\r
- fHistPt = new TH1F("fHistPt","p_{T} distribution",200,0,10);\r
- fList->Add(fHistPt);\r
- fHistEta = new TH1F("fHistEta","#eta distribution",200,-2,2);\r
- fList->Add(fHistEta);\r
- fHistRapidity = new TH1F("fHistRapidity","y distribution",200,-2,2);\r
- fList->Add(fHistRapidity);\r
- fHistPhi = new TH1F("fHistPhi","#phi distribution",200,-20,380);\r
- fList->Add(fHistPhi);\r
- fHistPhiBefore = new TH1F("fHistPhiBefore","#phi distribution",200,0.,2*TMath::Pi());\r
- fList->Add(fHistPhiBefore);\r
- fHistPhiAfter = new TH1F("fHistPhiAfter","#phi distribution",200,0.,2*TMath::Pi());\r
- fList->Add(fHistPhiAfter);\r
- fHistPhiPos = new TH1F("fHistPhiPos","#phi distribution for positive particles",200,-20,380);\r
- fList->Add(fHistPhiPos);\r
- fHistPhiNeg = new TH1F("fHistPhiNeg","#phi distribution for negative particles",200,-20,380);\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
- TString gRefTrackName[6] = {"tracks","tracksPos","tracksNeg","tracksTPConly","clusITS0","clusITS1"};\r
- fHistRefTracks = new TH2F("fHistRefTracks","Nr of Ref tracks/event vs. ref track estimator;;Nr of tracks",6, 0, 6, 400, 0, 20000);\r
- for(Int_t i = 1; i <= 6; i++)\r
- fHistRefTracks->GetXaxis()->SetBinLabel(i,gRefTrackName[i-1].Data());\r
- fList->Add(fHistRefTracks);\r
-\r
- // Balance function histograms\r
- // Initialize histograms if not done yet\r
- if(!fBalance->GetHistNp(0)){\r
- AliWarning("Histograms not yet initialized! --> Will be done now");\r
- AliWarning("--> Add 'gBalance->InitHistograms()' in your configBalanceFunction");\r
- fBalance->InitHistograms();\r
- }\r
-\r
- if(fRunShuffling) {\r
- if(!fShuffledBalance->GetHistNp(0)) {\r
- AliWarning("Histograms (shuffling) not yet initialized! --> Will be done now");\r
- AliWarning("--> Add 'gBalance->InitHistograms()' in your configBalanceFunction");\r
- fShuffledBalance->InitHistograms();\r
- }\r
- }\r
-\r
- for(Int_t a = 0; a < ANALYSIS_TYPES; a++){\r
- fListBF->Add(fBalance->GetHistNp(a));\r
- fListBF->Add(fBalance->GetHistNn(a));\r
- fListBF->Add(fBalance->GetHistNpn(a));\r
- fListBF->Add(fBalance->GetHistNnn(a));\r
- fListBF->Add(fBalance->GetHistNpp(a));\r
- fListBF->Add(fBalance->GetHistNnp(a));\r
-\r
- if(fRunShuffling) {\r
- fListBFS->Add(fShuffledBalance->GetHistNp(a));\r
- fListBFS->Add(fShuffledBalance->GetHistNn(a));\r
- fListBFS->Add(fShuffledBalance->GetHistNpn(a));\r
- fListBFS->Add(fShuffledBalance->GetHistNnn(a));\r
- fListBFS->Add(fShuffledBalance->GetHistNpp(a));\r
- fListBFS->Add(fShuffledBalance->GetHistNnp(a));\r
- } \r
- }\r
-\r
- if(fESDtrackCuts) fList->Add(fESDtrackCuts);\r
-\r
- //====================PID========================//\r
- if(fUsePID) {\r
- fPIDCombined = new AliPIDCombined();\r
- fPIDCombined->SetDefaultTPCPriors();\r
-\r
- fHistdEdxVsPTPCbeforePID = new TH2D ("dEdxVsPTPCbefore","dEdxVsPTPCbefore", 1000, -10.0, 10.0, 1000, 0, 1000); \r
- fHistListPIDQA->Add(fHistdEdxVsPTPCbeforePID); //addition \r
- \r
- fHistBetavsPTOFbeforePID = new TH2D ("BetavsPTOFbefore","BetavsPTOFbefore", 1000, -10.0, 10., 1000, 0, 1.2); \r
- fHistListPIDQA->Add(fHistBetavsPTOFbeforePID); //addition\r
- \r
- fHistProbTPCvsPtbeforePID = new TH2D ("ProbTPCvsPtbefore","ProbTPCvsPtbefore", 1000, -10.0,10.0, 1000, 0, 2.0); \r
- fHistListPIDQA->Add(fHistProbTPCvsPtbeforePID); //addition \r
- \r
- fHistProbTOFvsPtbeforePID = new TH2D ("ProbTOFvsPtbefore","ProbTOFvsPtbefore", 1000, -50, 50, 1000, 0, 2.0); \r
- fHistListPIDQA->Add(fHistProbTOFvsPtbeforePID); //addition \r
-\r
- fHistProbTPCTOFvsPtbeforePID =new TH2D ("ProbTPCTOFvsPtbefore","ProbTPCTOFvsPtbefore", 1000, -50, 50, 1000, 0, 2.0); \r
- fHistListPIDQA->Add(fHistProbTPCTOFvsPtbeforePID); //addition \r
- \r
- fHistNSigmaTPCvsPtbeforePID = new TH2D ("NSigmaTPCvsPtbefore","NSigmaTPCvsPtbefore", 1000, -10, 10, 1000, 0, 500); \r
- fHistListPIDQA->Add(fHistNSigmaTPCvsPtbeforePID); //addition \r
- \r
- fHistNSigmaTOFvsPtbeforePID = new TH2D ("NSigmaTOFvsPtbefore","NSigmaTOFvsPtbefore", 1000, -10, 10, 1000, 0, 500); \r
- fHistListPIDQA->Add(fHistNSigmaTOFvsPtbeforePID); //addition \r
- \r
- fHistdEdxVsPTPCafterPID = new TH2D ("dEdxVsPTPCafter","dEdxVsPTPCafter", 1000, -10, 10, 1000, 0, 1000); \r
- fHistListPIDQA->Add(fHistdEdxVsPTPCafterPID); //addition \r
- \r
- fHistBetavsPTOFafterPID = new TH2D ("BetavsPTOFafter","BetavsPTOFafter", 1000, -10, 10, 1000, 0, 1.2); \r
- fHistListPIDQA->Add(fHistBetavsPTOFafterPID); //addition \r
- \r
- fHistProbTPCvsPtafterPID = new TH2D ("ProbTPCvsPtafter","ProbTPCvsPtafter", 1000, -10, 10, 1000, 0, 2); \r
- fHistListPIDQA->Add(fHistProbTPCvsPtafterPID); //addition \r
- \r
- fHistProbTOFvsPtafterPID = new TH2D ("ProbTOFvsPtafter","ProbTOFvsPtafter", 1000, -10, 10, 1000, 0, 2); \r
- fHistListPIDQA->Add(fHistProbTOFvsPtafterPID); //addition \r
- \r
- fHistProbTPCTOFvsPtafterPID =new TH2D ("ProbTPCTOFvsPtafter","ProbTPCTOFvsPtafter", 1000, -50, 50, 1000, 0, 2.0); \r
- fHistListPIDQA->Add(fHistProbTPCTOFvsPtafterPID); //addition \r
-\r
- fHistNSigmaTPCvsPtafterPID = new TH2D ("NSigmaTPCvsPtafter","NSigmaTPCvsPtafter", 1000, -10, 10, 1000, 0, 500); \r
- fHistListPIDQA->Add(fHistNSigmaTPCvsPtafterPID); //addition \r
- \r
- fHistNSigmaTOFvsPtafterPID = new TH2D ("NSigmaTOFvsPtafter","NSigmaTOFvsPtafter", 1000, -10, 10, 1000, 0, 500); \r
- fHistListPIDQA->Add(fHistNSigmaTOFvsPtafterPID); //addition \r
- }\r
- //====================PID========================//\r
-\r
- // Post output data.\r
- PostData(1, fList);\r
- PostData(2, fListBF);\r
- if(fRunShuffling) PostData(3, fListBFS);\r
- if(fUsePID) PostData(4, fHistListPIDQA); //PID\r
-}\r
-\r
-//________________________________________________________________________\r
-void AliAnalysisTaskBF::UserExec(Option_t *) {\r
- // Main loop\r
- // Called for each event\r
- TString gAnalysisLevel = fBalance->GetAnalysisLevel();\r
-\r
- AliESDtrack *track_TPC = NULL;\r
-\r
- Int_t gNumberOfAcceptedTracks = 0;\r
- Float_t fCentrality = 0.;\r
-\r
- // for HBT like cuts need magnetic field sign\r
- Float_t bSign = 0; // only used in AOD so far\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
- }\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
- 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
- return;\r
- }\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); //all events\r
- Bool_t isSelected = kTRUE;\r
- if(fUseOfflineTrigger)\r
- isSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();\r
- if(isSelected) {\r
- fHistEventStats->Fill(2); //triggered events\r
-\r
- if(fUseCentrality) {\r
- //Centrality stuff\r
- AliCentrality *centrality = gESD->GetCentrality();\r
- \r
- fCentrality = centrality->GetCentralityPercentile(fCentralityEstimator.Data());\r
- \r
- // take only events inside centrality class\r
- if(!centrality->IsEventInCentralityClass(fCentralityPercentileMin,\r
- fCentralityPercentileMax,\r
- fCentralityEstimator.Data()))\r
- return;\r
- \r
- // centrality QA (V0M)\r
- fHistV0M->Fill(gESD->GetVZEROData()->GetMTotV0A(), gESD->GetVZEROData()->GetMTotV0C());\r
- }\r
- \r
- const AliESDVertex *vertex = gESD->GetPrimaryVertex();\r
- if(vertex) {\r
- if(vertex->GetNContributors() > 0) {\r
- if(vertex->GetZRes() != 0) {\r
- fHistEventStats->Fill(3); //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); //analayzed events\r
- fHistVx->Fill(vertex->GetXv());\r
- fHistVy->Fill(vertex->GetYv());\r
- fHistVz->Fill(vertex->GetZv());\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);\r
- fHistPt->Fill(v_pt);\r
- fHistEta->Fill(v_eta);\r
- fHistPhi->Fill(v_phi);\r
- fHistRapidity->Fill(v_y);\r
- if(v_charge > 0) fHistPhiPos->Fill(v_phi);\r
- else if(v_charge < 0) fHistPhiNeg->Fill(v_phi);\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
- \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
- }//ESD analysis\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
- }\r
-\r
- // for HBT like cuts need magnetic field sign\r
- bSign = (gAOD->GetMagneticField() > 0) ? 1 : -1;\r
-\r
- AliAODHeader *aodHeader = gAOD->GetHeader();\r
-\r
- // store offline trigger bits\r
- fHistTriggerStats->Fill(aodHeader->GetOfflineTrigger());\r
- \r
- // event selection done in AliAnalysisTaskSE::Exec() --> this is not used\r
- fHistEventStats->Fill(1); //all events\r
- Bool_t isSelected = kTRUE;\r
- if(fUseOfflineTrigger)\r
- isSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();\r
- if(isSelected) {\r
- fHistEventStats->Fill(2); //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
- \r
- // take only events inside centrality class\r
- if((fCentrality < fCentralityPercentileMin) || (fCentrality > fCentralityPercentileMax)) \r
- return;\r
- \r
- // centrality QA (V0M)\r
- fHistV0M->Fill(gAOD->GetVZEROData()->GetMTotV0A(), gAOD->GetVZEROData()->GetMTotV0C());\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
-\r
- const AliAODVertex *vertex = gAOD->GetPrimaryVertex();\r
- \r
- if(vertex) {\r
- Double32_t fCov[6];\r
- vertex->GetCovarianceMatrix(fCov);\r
- \r
- if(vertex->GetNContributors() > 0) {\r
- if(fCov[5] != 0) {\r
- fHistEventStats->Fill(3); //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); //analyzed events\r
- fHistVx->Fill(vertex->GetX());\r
- fHistVy->Fill(vertex->GetY());\r
- fHistVz->Fill(vertex->GetZ());\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
-\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());\r
- fHistPt->Fill(v_pt);\r
- fHistEta->Fill(v_eta);\r
- fHistPhi->Fill(v_phi);\r
- fHistRapidity->Fill(v_y);\r
- if(v_charge > 0) fHistPhiPos->Fill(v_phi);\r
- else if(v_charge < 0) fHistPhiNeg->Fill(v_phi);\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
- if (!mcEvent) {\r
- Printf("ERROR: mcEvent not available");\r
- return;\r
- }\r
-\r
- AliESDEvent* gESD = dynamic_cast<AliESDEvent*>(InputEvent()); // from TaskSE\r
- if (!gESD) {\r
- Printf("ERROR: gESD not available");\r
- return;\r
- }\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); //all events\r
- Bool_t isSelected = kTRUE;\r
- if(fUseOfflineTrigger)\r
- isSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();\r
- if(isSelected) {\r
- fHistEventStats->Fill(2); //triggered events\r
-\r
- if(fUseCentrality) {\r
- //Centrality stuff\r
- AliCentrality *centrality = gESD->GetCentrality();\r
-\r
- fCentrality = centrality->GetCentralityPercentile(fCentralityEstimator.Data());\r
- \r
- // take only events inside centrality class\r
- if(!centrality->IsEventInCentralityClass(fCentralityPercentileMin,\r
- fCentralityPercentileMax,\r
- fCentralityEstimator.Data()))\r
- return;\r
- \r
- // centrality QA (V0M)\r
- fHistV0M->Fill(gESD->GetVZEROData()->GetMTotV0A(), gESD->GetVZEROData()->GetMTotV0C());\r
- }\r
- \r
- const AliESDVertex *vertex = gESD->GetPrimaryVertex();\r
- if(vertex) {\r
- if(vertex->GetNContributors() > 0) {\r
- if(vertex->GetZRes() != 0) {\r
- fHistEventStats->Fill(3); //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); //analayzed events\r
- fHistVx->Fill(vertex->GetXv());\r
- fHistVy->Fill(vertex->GetYv());\r
- fHistVz->Fill(vertex->GetZv());\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);\r
- fHistPt->Fill(v_pt);\r
- fHistEta->Fill(v_eta);\r
- fHistPhi->Fill(v_phi);\r
- fHistRapidity->Fill(v_y);\r
- if(v_charge > 0) fHistPhiPos->Fill(v_phi);\r
- else if(v_charge < 0) fHistPhiNeg->Fill(v_phi);\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
-\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); //total events\r
- fHistEventStats->Fill(2); //offline trigger\r
-\r
- Double_t gReactionPlane = 0., 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
- }\r
- fCentrality = gImpactParameter;\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
- \r
- AliGenEventHeader *header = mcEvent->GenEventHeader();\r
- if(!header) return;\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); //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); //analayzed events\r
- fHistVx->Fill(gVertexArray.At(0));\r
- fHistVy->Fill(gVertexArray.At(1));\r
- fHistVz->Fill(gVertexArray.At(2));\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
-\r
- v_eta = track->Eta();\r
- v_pt = track->Pt();\r
- v_y = track->Y();\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
-\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
- }\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);\r
- fHistEta->Fill(v_eta);\r
- fHistPhi->Fill(v_phi);\r
- fHistRapidity->Fill(v_y);\r
- if(v_charge > 0) fHistPhiPos->Fill(v_phi);\r
- else if(v_charge < 0) fHistPhiNeg->Fill(v_phi);\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);\r
-\r
- Double_t v_DeltaphiAfter = v_phi - gReactionPlane;\r
- if(v_DeltaphiAfter < 0) v_DeltaphiAfter += 2*TMath::Pi();\r
- fHistPhiAfter->Fill(v_DeltaphiAfter);\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
- }//Vz cut\r
- }//Vy cut\r
- }//Vx cut\r
- }//MC analysis\r
- \r
- //multiplicity cut (used in pp)\r
- if(fUseMultiplicity) {\r
- if((gNumberOfAcceptedTracks < fNumberOfAcceptedTracksMin)||(gNumberOfAcceptedTracks > fNumberOfAcceptedTracksMax))\r
- return;\r
- }\r
- fHistNumberOfAcceptedTracks->Fill(gNumberOfAcceptedTracks);\r
-\r
- // calculate balance function\r
- if(fUseMultiplicity) \r
- fBalance->CalculateBalance(gNumberOfAcceptedTracks,chargeVector,bSign);\r
- else \r
- fBalance->CalculateBalance(fCentrality,chargeVector,bSign);\r
-\r
- if(fRunShuffling) {\r
- // shuffle charges\r
- random_shuffle( chargeVectorShuffle[0]->begin(), chargeVectorShuffle[0]->end() );\r
- if(fUseMultiplicity) \r
- fShuffledBalance->CalculateBalance(gNumberOfAcceptedTracks,chargeVectorShuffle,bSign);\r
- else \r
- fShuffledBalance->CalculateBalance(fCentrality,chargeVectorShuffle,bSign);\r
- }\r
-} \r
-\r
-//________________________________________________________________________\r
-void AliAnalysisTaskBF::FinishTaskOutput(){\r
- //Printf("END BF");\r
-\r
- if (!fBalance) {\r
- Printf("ERROR: fBalance not available");\r
- return;\r
- } \r
- if(fRunShuffling) {\r
- if (!fShuffledBalance) {\r
- Printf("ERROR: fShuffledBalance not available");\r
- return;\r
- }\r
- }\r
-\r
-}\r
-\r
-//________________________________________________________________________\r
-void AliAnalysisTaskBF::Terminate(Option_t *) {\r
- // Draw result to the screen\r
- // Called once at the end of the query\r
-\r
- // not implemented ...\r
-\r
-}\r
+#include "TChain.h"
+#include "TList.h"
+#include "TCanvas.h"
+#include "TLorentzVector.h"
+#include "TGraphErrors.h"
+#include "TH1F.h"
+#include "TH2F.h"
+#include "TArrayF.h"
+#include "TF1.h"
+#include "TRandom.h"
+
+#include "AliAnalysisTaskSE.h"
+#include "AliAnalysisManager.h"
+
+#include "AliESDVertex.h"
+#include "AliESDEvent.h"
+#include "AliESDInputHandler.h"
+#include "AliAODEvent.h"
+#include "AliAODTrack.h"
+#include "AliAODInputHandler.h"
+#include "AliGenEventHeader.h"
+#include "AliGenHijingEventHeader.h"
+#include "AliMCEventHandler.h"
+#include "AliMCEvent.h"
+#include "AliStack.h"
+#include "AliESDtrackCuts.h"
+
+#include "TH2D.h"
+#include "AliPID.h"
+#include "AliPIDResponse.h"
+#include "AliPIDCombined.h"
+
+#include "AliAnalysisTaskBF.h"
+#include "AliBalance.h"
+
+
+// Analysis task for the BF code
+// Authors: Panos.Christakoglou@nikhef.nl
+
+ClassImp(AliAnalysisTaskBF)
+
+//________________________________________________________________________
+AliAnalysisTaskBF::AliAnalysisTaskBF(const char *name)
+: AliAnalysisTaskSE(name),
+ fBalance(0),
+ fRunShuffling(kFALSE),
+ fShuffledBalance(0),
+ fList(0),
+ fListBF(0),
+ fListBFS(0),
+ fHistListPIDQA(0),
+ fHistEventStats(0),
+ fHistCentStats(0),
+ fHistTriggerStats(0),
+ fHistTrackStats(0),
+ fHistVx(0),
+ fHistVy(0),
+ fHistVz(0),
+ fHistClus(0),
+ fHistDCA(0),
+ fHistChi2(0),
+ fHistPt(0),
+ fHistEta(0),
+ fHistRapidity(0),
+ fHistPhi(0),
+ fHistPhiBefore(0),
+ fHistPhiAfter(0),
+ fHistPhiPos(0),
+ fHistPhiNeg(0),
+ fHistV0M(0),
+ fHistRefTracks(0),
+ fHistdEdxVsPTPCbeforePID(NULL),
+ fHistBetavsPTOFbeforePID(NULL),
+ fHistProbTPCvsPtbeforePID(NULL),
+ fHistProbTOFvsPtbeforePID(NULL),
+ fHistProbTPCTOFvsPtbeforePID(NULL),
+ fHistNSigmaTPCvsPtbeforePID(NULL),
+ fHistNSigmaTOFvsPtbeforePID(NULL),
+ fHistdEdxVsPTPCafterPID(NULL),
+ fHistBetavsPTOFafterPID(NULL),
+ fHistProbTPCvsPtafterPID(NULL),
+ fHistProbTOFvsPtafterPID(NULL),
+ fHistProbTPCTOFvsPtafterPID(NULL),
+ fHistNSigmaTPCvsPtafterPID(NULL),
+ fHistNSigmaTOFvsPtafterPID(NULL),
+ fPIDResponse(0x0),
+ fPIDCombined(0x0),
+ fParticleOfInterest(kPion),
+ fPidDetectorConfig(kTPCTOF),
+ fUsePID(kFALSE),
+ fUsePIDnSigma(kTRUE),
+ fUsePIDPropabilities(kFALSE),
+ fPIDNSigma(3.),
+ fMinAcceptedPIDProbability(0.8),
+ fESDtrackCuts(0),
+ fCentralityEstimator("V0M"),
+ fUseCentrality(kFALSE),
+ fCentralityPercentileMin(0.),
+ fCentralityPercentileMax(5.),
+ fImpactParameterMin(0.),
+ fImpactParameterMax(20.),
+ fUseMultiplicity(kFALSE),
+ fNumberOfAcceptedTracksMin(0),
+ fNumberOfAcceptedTracksMax(10000),
+ fHistNumberOfAcceptedTracks(0),
+ fUseOfflineTrigger(kFALSE),
+ fVxMax(0.3),
+ fVyMax(0.3),
+ fVzMax(10.),
+ nAODtrackCutBit(128),
+ fPtMin(0.3),
+ fPtMax(1.5),
+ fEtaMin(-0.8),
+ fEtaMax(-0.8),
+ fDCAxyCut(-1),
+ fDCAzCut(-1),
+ fTPCchi2Cut(-1),
+ fNClustersTPCCut(-1),
+ fAcceptanceParameterization(0),
+ fDifferentialV2(0),
+ fUseFlowAfterBurner(kFALSE),
+ fExcludeResonancesInMC(kFALSE),
+ fUseMCPdgCode(kFALSE),
+ fPDGCodeToBeAnalyzed(-1) {
+ // Constructor
+ // Define input and output slots here
+ // Input slot #0 works with a TChain
+ DefineInput(0, TChain::Class());
+ // Output slot #0 writes into a TH1 container
+ DefineOutput(1, TList::Class());
+ DefineOutput(2, TList::Class());
+ DefineOutput(3, TList::Class());
+ DefineOutput(4, TList::Class());
+}
+
+//________________________________________________________________________
+AliAnalysisTaskBF::~AliAnalysisTaskBF() {
+
+ // delete fBalance;
+ // delete fShuffledBalance;
+ // delete fList;
+ // delete fListBF;
+ // delete fListBFS;
+
+ // delete fHistEventStats;
+ // delete fHistTrackStats;
+ // delete fHistVx;
+ // delete fHistVy;
+ // delete fHistVz;
+
+ // delete fHistClus;
+ // delete fHistDCA;
+ // delete fHistChi2;
+ // delete fHistPt;
+ // delete fHistEta;
+ // delete fHistPhi;
+ // delete fHistV0M;
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskBF::UserCreateOutputObjects() {
+ // Create histograms
+ // Called once
+ if(!fBalance) {
+ fBalance = new AliBalance();
+ fBalance->SetAnalysisLevel("ESD");
+ //fBalance->SetNumberOfBins(-1,16);
+ fBalance->SetInterval(-1,-0.8,0.8,16,0.,1.6);
+ }
+ if(fRunShuffling) {
+ if(!fShuffledBalance) {
+ fShuffledBalance = new AliBalance();
+ fShuffledBalance->SetAnalysisLevel("ESD");
+ //fShuffledBalance->SetNumberOfBins(-1,16);
+ fShuffledBalance->SetInterval(-1,-0.8,0.8,16,0.,1.6);
+ }
+ }
+
+ //QA list
+ fList = new TList();
+ fList->SetName("listQA");
+ fList->SetOwner();
+
+ //Balance Function list
+ fListBF = new TList();
+ fListBF->SetName("listBF");
+ fListBF->SetOwner();
+
+ if(fRunShuffling) {
+ fListBFS = new TList();
+ fListBFS->SetName("listBFShuffled");
+ fListBFS->SetOwner();
+ }
+
+ //PID QA list
+ if(fUsePID) {
+ fHistListPIDQA = new TList();
+ fHistListPIDQA->SetName("listQAPID");
+ fHistListPIDQA->SetOwner();
+ }
+
+ //Event stats.
+ TString gCutName[4] = {"Total","Offline trigger",
+ "Vertex","Analyzed"};
+ fHistEventStats = new TH1F("fHistEventStats",
+ "Event statistics;;N_{events}",
+ 4,0.5,4.5);
+ for(Int_t i = 1; i <= 4; i++)
+ fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data());
+ fList->Add(fHistEventStats);
+
+ TString gCentName[9] = {"V0M","FMD","TRK","TKL","CL0","CL1","V0MvsFMD","TKLvsV0M","ZEMvsZDC"};
+ fHistCentStats = new TH2F("fHistCentStats",
+ "Centrality statistics;;Cent percentile",
+ 9,-0.5,8.5,220,-5,105);
+ for(Int_t i = 1; i <= 9; i++)
+ fHistCentStats->GetXaxis()->SetBinLabel(i,gCentName[i-1].Data());
+ fList->Add(fHistCentStats);
+
+ fHistTriggerStats = new TH1F("fHistTriggerStats","Trigger statistics;TriggerBit;N_{events}",130,0,130);
+ fList->Add(fHistTriggerStats);
+
+ fHistTrackStats = new TH1F("fHistTrackStats","Event statistics;TrackFilterBit;N_{events}",130,0,130);
+ fList->Add(fHistTrackStats);
+
+ fHistNumberOfAcceptedTracks = new TH1F("fHistNumberOfAcceptedTracks",";N_{acc.};Entries",4001,-0.5,4000.5);
+ fList->Add(fHistNumberOfAcceptedTracks);
+
+ // Vertex distributions
+ fHistVx = new TH1F("fHistVx","Primary vertex distribution - x coordinate;V_{x} (cm);Entries",100,-0.5,0.5);
+ fList->Add(fHistVx);
+ fHistVy = new TH1F("fHistVy","Primary vertex distribution - y coordinate;V_{y} (cm);Entries",100,-0.5,0.5);
+ fList->Add(fHistVy);
+ fHistVz = new TH1F("fHistVz","Primary vertex distribution - z coordinate;V_{z} (cm);Entries",100,-20.,20.);
+ fList->Add(fHistVz);
+
+ // QA histograms
+ fHistClus = new TH2F("fHistClus","# Cluster (TPC vs. ITS)",10,0,10,200,0,200);
+ fList->Add(fHistClus);
+ fHistChi2 = new TH1F("fHistChi2","Chi2/NDF distribution",200,0,10);
+ fList->Add(fHistChi2);
+ fHistDCA = new TH2F("fHistDCA","DCA (xy vs. z)",400,-5,5,400,-5,5);
+ fList->Add(fHistDCA);
+ fHistPt = new TH1F("fHistPt","p_{T} distribution",200,0,10);
+ fList->Add(fHistPt);
+ fHistEta = new TH1F("fHistEta","#eta distribution",200,-2,2);
+ fList->Add(fHistEta);
+ fHistRapidity = new TH1F("fHistRapidity","y distribution",200,-2,2);
+ fList->Add(fHistRapidity);
+ fHistPhi = new TH1F("fHistPhi","#phi distribution",200,-20,380);
+ fList->Add(fHistPhi);
+ fHistPhiBefore = new TH1F("fHistPhiBefore","#phi distribution",200,0.,2*TMath::Pi());
+ fList->Add(fHistPhiBefore);
+ fHistPhiAfter = new TH1F("fHistPhiAfter","#phi distribution",200,0.,2*TMath::Pi());
+ fList->Add(fHistPhiAfter);
+ fHistPhiPos = new TH1F("fHistPhiPos","#phi distribution for positive particles",200,-20,380);
+ fList->Add(fHistPhiPos);
+ fHistPhiNeg = new TH1F("fHistPhiNeg","#phi distribution for negative particles",200,-20,380);
+ fList->Add(fHistPhiNeg);
+ fHistV0M = new TH2F("fHistV0M","V0 Multiplicity C vs. A",500, 0, 20000, 500, 0, 20000);
+ fList->Add(fHistV0M);
+ TString gRefTrackName[6] = {"tracks","tracksPos","tracksNeg","tracksTPConly","clusITS0","clusITS1"};
+ fHistRefTracks = new TH2F("fHistRefTracks","Nr of Ref tracks/event vs. ref track estimator;;Nr of tracks",6, 0, 6, 400, 0, 20000);
+ for(Int_t i = 1; i <= 6; i++)
+ fHistRefTracks->GetXaxis()->SetBinLabel(i,gRefTrackName[i-1].Data());
+ fList->Add(fHistRefTracks);
+
+ // Balance function histograms
+ // Initialize histograms if not done yet
+ if(!fBalance->GetHistNp(0)){
+ AliWarning("Histograms not yet initialized! --> Will be done now");
+ AliWarning("--> Add 'gBalance->InitHistograms()' in your configBalanceFunction");
+ fBalance->InitHistograms();
+ }
+
+ if(fRunShuffling) {
+ if(!fShuffledBalance->GetHistNp(0)) {
+ AliWarning("Histograms (shuffling) not yet initialized! --> Will be done now");
+ AliWarning("--> Add 'gBalance->InitHistograms()' in your configBalanceFunction");
+ fShuffledBalance->InitHistograms();
+ }
+ }
+
+ for(Int_t a = 0; a < ANALYSIS_TYPES; a++){
+ fListBF->Add(fBalance->GetHistNp(a));
+ fListBF->Add(fBalance->GetHistNn(a));
+ fListBF->Add(fBalance->GetHistNpn(a));
+ fListBF->Add(fBalance->GetHistNnn(a));
+ fListBF->Add(fBalance->GetHistNpp(a));
+ fListBF->Add(fBalance->GetHistNnp(a));
+
+ if(fRunShuffling) {
+ fListBFS->Add(fShuffledBalance->GetHistNp(a));
+ fListBFS->Add(fShuffledBalance->GetHistNn(a));
+ fListBFS->Add(fShuffledBalance->GetHistNpn(a));
+ fListBFS->Add(fShuffledBalance->GetHistNnn(a));
+ fListBFS->Add(fShuffledBalance->GetHistNpp(a));
+ fListBFS->Add(fShuffledBalance->GetHistNnp(a));
+ }
+ }
+
+ if(fESDtrackCuts) fList->Add(fESDtrackCuts);
+
+ //====================PID========================//
+ if(fUsePID) {
+ fPIDCombined = new AliPIDCombined();
+ fPIDCombined->SetDefaultTPCPriors();
+
+ fHistdEdxVsPTPCbeforePID = new TH2D ("dEdxVsPTPCbefore","dEdxVsPTPCbefore", 1000, -10.0, 10.0, 1000, 0, 1000);
+ fHistListPIDQA->Add(fHistdEdxVsPTPCbeforePID); //addition
+
+ fHistBetavsPTOFbeforePID = new TH2D ("BetavsPTOFbefore","BetavsPTOFbefore", 1000, -10.0, 10., 1000, 0, 1.2);
+ fHistListPIDQA->Add(fHistBetavsPTOFbeforePID); //addition
+
+ fHistProbTPCvsPtbeforePID = new TH2D ("ProbTPCvsPtbefore","ProbTPCvsPtbefore", 1000, -10.0,10.0, 1000, 0, 2.0);
+ fHistListPIDQA->Add(fHistProbTPCvsPtbeforePID); //addition
+
+ fHistProbTOFvsPtbeforePID = new TH2D ("ProbTOFvsPtbefore","ProbTOFvsPtbefore", 1000, -50, 50, 1000, 0, 2.0);
+ fHistListPIDQA->Add(fHistProbTOFvsPtbeforePID); //addition
+
+ fHistProbTPCTOFvsPtbeforePID =new TH2D ("ProbTPCTOFvsPtbefore","ProbTPCTOFvsPtbefore", 1000, -50, 50, 1000, 0, 2.0);
+ fHistListPIDQA->Add(fHistProbTPCTOFvsPtbeforePID); //addition
+
+ fHistNSigmaTPCvsPtbeforePID = new TH2D ("NSigmaTPCvsPtbefore","NSigmaTPCvsPtbefore", 1000, -10, 10, 1000, 0, 500);
+ fHistListPIDQA->Add(fHistNSigmaTPCvsPtbeforePID); //addition
+
+ fHistNSigmaTOFvsPtbeforePID = new TH2D ("NSigmaTOFvsPtbefore","NSigmaTOFvsPtbefore", 1000, -10, 10, 1000, 0, 500);
+ fHistListPIDQA->Add(fHistNSigmaTOFvsPtbeforePID); //addition
+
+ fHistdEdxVsPTPCafterPID = new TH2D ("dEdxVsPTPCafter","dEdxVsPTPCafter", 1000, -10, 10, 1000, 0, 1000);
+ fHistListPIDQA->Add(fHistdEdxVsPTPCafterPID); //addition
+
+ fHistBetavsPTOFafterPID = new TH2D ("BetavsPTOFafter","BetavsPTOFafter", 1000, -10, 10, 1000, 0, 1.2);
+ fHistListPIDQA->Add(fHistBetavsPTOFafterPID); //addition
+
+ fHistProbTPCvsPtafterPID = new TH2D ("ProbTPCvsPtafter","ProbTPCvsPtafter", 1000, -10, 10, 1000, 0, 2);
+ fHistListPIDQA->Add(fHistProbTPCvsPtafterPID); //addition
+
+ fHistProbTOFvsPtafterPID = new TH2D ("ProbTOFvsPtafter","ProbTOFvsPtafter", 1000, -10, 10, 1000, 0, 2);
+ fHistListPIDQA->Add(fHistProbTOFvsPtafterPID); //addition
+
+ fHistProbTPCTOFvsPtafterPID =new TH2D ("ProbTPCTOFvsPtafter","ProbTPCTOFvsPtafter", 1000, -50, 50, 1000, 0, 2.0);
+ fHistListPIDQA->Add(fHistProbTPCTOFvsPtafterPID); //addition
+
+ fHistNSigmaTPCvsPtafterPID = new TH2D ("NSigmaTPCvsPtafter","NSigmaTPCvsPtafter", 1000, -10, 10, 1000, 0, 500);
+ fHistListPIDQA->Add(fHistNSigmaTPCvsPtafterPID); //addition
+
+ fHistNSigmaTOFvsPtafterPID = new TH2D ("NSigmaTOFvsPtafter","NSigmaTOFvsPtafter", 1000, -10, 10, 1000, 0, 500);
+ fHistListPIDQA->Add(fHistNSigmaTOFvsPtafterPID); //addition
+ }
+ //====================PID========================//
+
+ // Post output data.
+ PostData(1, fList);
+ PostData(2, fListBF);
+ if(fRunShuffling) PostData(3, fListBFS);
+ if(fUsePID) PostData(4, fHistListPIDQA); //PID
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskBF::UserExec(Option_t *) {
+ // Main loop
+ // Called for each event
+ TString gAnalysisLevel = fBalance->GetAnalysisLevel();
+
+ AliESDtrack *track_TPC = NULL;
+
+ Int_t gNumberOfAcceptedTracks = 0;
+ Float_t fCentrality = 0.;
+
+ // for HBT like cuts need magnetic field sign
+ Float_t bSign = 0; // only used in AOD so far
+
+ // vector holding the charges/kinematics of all tracks (charge,y,eta,phi,p0,p1,p2,pt,E)
+ vector<Double_t> *chargeVectorShuffle[9]; // this will be shuffled
+ vector<Double_t> *chargeVector[9]; // original charge
+ for(Int_t i = 0; i < 9; i++){
+ chargeVectorShuffle[i] = new vector<Double_t>;
+ chargeVector[i] = new vector<Double_t>;
+ }
+
+ Double_t v_charge;
+ Double_t v_y;
+ Double_t v_eta;
+ Double_t v_phi;
+ Double_t v_p[3];
+ Double_t v_pt;
+ Double_t v_E;
+
+ if(fUsePID) {
+ fPIDResponse = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetPIDResponse();
+ if (!fPIDResponse) AliFatal("This Task needs the PID response attached to the inputHandler");
+ }
+
+ //ESD analysis
+ if(gAnalysisLevel == "ESD") {
+ AliESDEvent* gESD = dynamic_cast<AliESDEvent*>(InputEvent()); // from TaskSE
+ if (!gESD) {
+ Printf("ERROR: gESD not available");
+ return;
+ }
+
+ // store offline trigger bits
+ fHistTriggerStats->Fill(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected());
+
+ // event selection done in AliAnalysisTaskSE::Exec() --> this is not used
+ fHistEventStats->Fill(1); //all events
+ Bool_t isSelected = kTRUE;
+ if(fUseOfflineTrigger)
+ isSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
+ if(isSelected) {
+ fHistEventStats->Fill(2); //triggered events
+
+ if(fUseCentrality) {
+ //Centrality stuff
+ AliCentrality *centrality = gESD->GetCentrality();
+
+ fCentrality = centrality->GetCentralityPercentile(fCentralityEstimator.Data());
+
+ // take only events inside centrality class
+ if(!centrality->IsEventInCentralityClass(fCentralityPercentileMin,
+ fCentralityPercentileMax,
+ fCentralityEstimator.Data()))
+ return;
+
+ // centrality QA (V0M)
+ fHistV0M->Fill(gESD->GetVZEROData()->GetMTotV0A(), gESD->GetVZEROData()->GetMTotV0C());
+ }
+
+ const AliESDVertex *vertex = gESD->GetPrimaryVertex();
+ if(vertex) {
+ if(vertex->GetNContributors() > 0) {
+ if(vertex->GetZRes() != 0) {
+ fHistEventStats->Fill(3); //events with a proper vertex
+ if(TMath::Abs(vertex->GetXv()) < fVxMax) {
+ if(TMath::Abs(vertex->GetYv()) < fVyMax) {
+ if(TMath::Abs(vertex->GetZv()) < fVzMax) {
+ fHistEventStats->Fill(4); //analayzed events
+ fHistVx->Fill(vertex->GetXv());
+ fHistVy->Fill(vertex->GetYv());
+ fHistVz->Fill(vertex->GetZv());
+
+ //Printf("There are %d tracks in this event", gESD->GetNumberOfTracks());
+ for (Int_t iTracks = 0; iTracks < gESD->GetNumberOfTracks(); iTracks++) {
+ AliESDtrack* track = dynamic_cast<AliESDtrack *>(gESD->GetTrack(iTracks));
+ if (!track) {
+ Printf("ERROR: Could not receive track %d", iTracks);
+ continue;
+ }
+
+ // take only TPC only tracks
+ track_TPC = new AliESDtrack();
+ if(!track->FillTPCOnlyTrack(*track_TPC)) continue;
+
+ //ESD track cuts
+ if(fESDtrackCuts)
+ if(!fESDtrackCuts->AcceptTrack(track_TPC)) continue;
+
+ // fill QA histograms
+ Float_t b[2];
+ Float_t bCov[3];
+ track_TPC->GetImpactParameters(b,bCov);
+ if (bCov[0]<=0 || bCov[2]<=0) {
+ AliDebug(1, "Estimated b resolution lower or equal zero!");
+ bCov[0]=0; bCov[2]=0;
+ }
+
+ Int_t nClustersTPC = -1;
+ nClustersTPC = track_TPC->GetTPCNclsIter1(); // TPC standalone
+ //nClustersTPC = track->GetTPCclusters(0); // global track
+ Float_t chi2PerClusterTPC = -1;
+ if (nClustersTPC!=0) {
+ chi2PerClusterTPC = track_TPC->GetTPCchi2Iter1()/Float_t(nClustersTPC); // TPC standalone
+ //chi2PerClusterTPC = track->GetTPCchi2()/Float_t(nClustersTPC); // global track
+ }
+
+ //===========================PID===============================//
+ if(fUsePID) {
+ Double_t prob[AliPID::kSPECIES]={0.};
+ Double_t probTPC[AliPID::kSPECIES]={0.};
+ Double_t probTOF[AliPID::kSPECIES]={0.};
+ Double_t probTPCTOF[AliPID::kSPECIES]={0.};
+
+ Double_t nSigma = 0.;
+ UInt_t detUsedTPC = 0;
+ UInt_t detUsedTOF = 0;
+ UInt_t detUsedTPCTOF = 0;
+
+ //Decide what detector configuration we want to use
+ switch(fPidDetectorConfig) {
+ case kTPCpid:
+ fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
+ nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)fParticleOfInterest));
+ detUsedTPC = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTPC);
+ for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
+ prob[iSpecies] = probTPC[iSpecies];
+ break;
+ case kTOFpid:
+ fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF);
+ nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)fParticleOfInterest));
+ detUsedTOF = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTOF);
+ for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
+ prob[iSpecies] = probTOF[iSpecies];
+ break;
+ case kTPCTOF:
+ fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC);
+ detUsedTPCTOF = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTPCTOF);
+ for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
+ prob[iSpecies] = probTPCTOF[iSpecies];
+ break;
+ default:
+ break;
+ }//end switch: define detector mask
+
+ //Filling the PID QA
+ Double_t tofTime = -999., length = 999., tof = -999.;
+ Double_t c = TMath::C()*1.E-9;// m/ns
+ Double_t beta = -999.;
+ Double_t nSigmaTOFForParticleOfInterest = -999.;
+ if ( (track->IsOn(AliESDtrack::kTOFin)) &&
+ (track->IsOn(AliESDtrack::kTIME)) ) {
+ tofTime = track->GetTOFsignal();//in ps
+ length = track->GetIntegratedLength();
+ tof = tofTime*1E-3; // ns
+
+ if (tof <= 0) {
+ //Printf("WARNING: track with negative TOF time found! Skipping this track for PID checks\n");
+ continue;
+ }
+ if (length <= 0){
+ //printf("WARNING: track with negative length found!Skipping this track for PID checks\n");
+ continue;
+ }
+
+ length = length*0.01; // in meters
+ tof = tof*c;
+ beta = length/tof;
+
+ nSigmaTOFForParticleOfInterest = fPIDResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)fParticleOfInterest);
+ fHistBetavsPTOFbeforePID ->Fill(track->P()*track->Charge(),beta);
+ fHistProbTOFvsPtbeforePID ->Fill(track->Pt(),probTOF[fParticleOfInterest]);
+ fHistNSigmaTOFvsPtbeforePID ->Fill(track->Pt(),nSigmaTOFForParticleOfInterest);
+ }//TOF signal
+
+
+ Double_t nSigmaTPCForParticleOfInterest = fPIDResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)fParticleOfInterest);
+ fHistdEdxVsPTPCbeforePID -> Fill(track->P()*track->Charge(),track->GetTPCsignal());
+ fHistProbTPCvsPtbeforePID -> Fill(track->Pt(),probTPC[fParticleOfInterest]);
+ fHistNSigmaTPCvsPtbeforePID -> Fill(track->Pt(),nSigmaTPCForParticleOfInterest);
+ fHistProbTPCTOFvsPtbeforePID -> Fill(track->Pt(),probTPCTOF[fParticleOfInterest]);
+ //end of QA-before pid
+
+ if ((detUsedTPC != 0)||(detUsedTOF != 0)||(detUsedTPCTOF != 0)) {
+ //Make the decision based on the n-sigma
+ if(fUsePIDnSigma) {
+ if(nSigma > fPIDNSigma) continue;}
+
+ //Make the decision based on the bayesian
+ else if(fUsePIDPropabilities) {
+ if(fParticleOfInterest != TMath::LocMax(AliPID::kSPECIES,prob)) continue;
+ if (prob[fParticleOfInterest] < fMinAcceptedPIDProbability) continue;
+ }
+
+ //Fill QA after the PID
+ fHistBetavsPTOFafterPID ->Fill(track->P()*track->Charge(),beta);
+ fHistProbTOFvsPtafterPID ->Fill(track->Pt(),probTOF[fParticleOfInterest]);
+ fHistNSigmaTOFvsPtafterPID ->Fill(track->Pt(),nSigmaTOFForParticleOfInterest);
+
+ fHistdEdxVsPTPCafterPID -> Fill(track->P()*track->Charge(),track->GetTPCsignal());
+ fHistProbTPCvsPtafterPID -> Fill(track->Pt(),probTPC[fParticleOfInterest]);
+ fHistProbTPCTOFvsPtafterPID -> Fill(track->Pt(),probTPCTOF[fParticleOfInterest]);
+ fHistNSigmaTPCvsPtafterPID -> Fill(track->Pt(),nSigmaTPCForParticleOfInterest);
+ }
+
+ PostData(4, fHistListPIDQA);
+ }
+ //===========================PID===============================//
+ v_charge = track_TPC->Charge();
+ v_y = track_TPC->Y();
+ v_eta = track_TPC->Eta();
+ v_phi = track_TPC->Phi() * TMath::RadToDeg();
+ v_E = track_TPC->E();
+ v_pt = track_TPC->Pt();
+ track_TPC->PxPyPz(v_p);
+ fHistClus->Fill(track_TPC->GetITSclusters(0),nClustersTPC);
+ fHistDCA->Fill(b[1],b[0]);
+ fHistChi2->Fill(chi2PerClusterTPC);
+ fHistPt->Fill(v_pt);
+ fHistEta->Fill(v_eta);
+ fHistPhi->Fill(v_phi);
+ fHistRapidity->Fill(v_y);
+ if(v_charge > 0) fHistPhiPos->Fill(v_phi);
+ else if(v_charge < 0) fHistPhiNeg->Fill(v_phi);
+
+ // fill charge vector
+ chargeVector[0]->push_back(v_charge);
+ chargeVector[1]->push_back(v_y);
+ chargeVector[2]->push_back(v_eta);
+ chargeVector[3]->push_back(v_phi);
+ chargeVector[4]->push_back(v_p[0]);
+ chargeVector[5]->push_back(v_p[1]);
+ chargeVector[6]->push_back(v_p[2]);
+ chargeVector[7]->push_back(v_pt);
+ chargeVector[8]->push_back(v_E);
+
+ if(fRunShuffling) {
+ chargeVectorShuffle[0]->push_back(v_charge);
+ chargeVectorShuffle[1]->push_back(v_y);
+ chargeVectorShuffle[2]->push_back(v_eta);
+ chargeVectorShuffle[3]->push_back(v_phi);
+ chargeVectorShuffle[4]->push_back(v_p[0]);
+ chargeVectorShuffle[5]->push_back(v_p[1]);
+ chargeVectorShuffle[6]->push_back(v_p[2]);
+ chargeVectorShuffle[7]->push_back(v_pt);
+ chargeVectorShuffle[8]->push_back(v_E);
+ }
+
+ delete track_TPC;
+
+ } //track loop
+ }//Vz cut
+ }//Vy cut
+ }//Vx cut
+ }//proper vertex resolution
+ }//proper number of contributors
+ }//vertex object valid
+ }//triggered event
+ }//ESD analysis
+
+ //AOD analysis (vertex and track cuts also here!!!!)
+ else if(gAnalysisLevel == "AOD") {
+ AliAODEvent* gAOD = dynamic_cast<AliAODEvent*>(InputEvent()); // from TaskSE
+ if(!gAOD) {
+ Printf("ERROR: gAOD not available");
+ return;
+ }
+
+ // for HBT like cuts need magnetic field sign
+ bSign = (gAOD->GetMagneticField() > 0) ? 1 : -1;
+
+ AliAODHeader *aodHeader = gAOD->GetHeader();
+
+ // store offline trigger bits
+ fHistTriggerStats->Fill(aodHeader->GetOfflineTrigger());
+
+ // event selection done in AliAnalysisTaskSE::Exec() --> this is not used
+ fHistEventStats->Fill(1); //all events
+ Bool_t isSelected = kTRUE;
+ if(fUseOfflineTrigger)
+ isSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
+ if(isSelected) {
+ fHistEventStats->Fill(2); //triggered events
+
+ //Centrality stuff (centrality in AOD header)
+ if(fUseCentrality) {
+ fCentrality = aodHeader->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());
+
+ // QA for centrality estimators
+ fHistCentStats->Fill(0.,aodHeader->GetCentralityP()->GetCentralityPercentile("V0M"));
+ fHistCentStats->Fill(1.,aodHeader->GetCentralityP()->GetCentralityPercentile("FMD"));
+ fHistCentStats->Fill(2.,aodHeader->GetCentralityP()->GetCentralityPercentile("TRK"));
+ fHistCentStats->Fill(3.,aodHeader->GetCentralityP()->GetCentralityPercentile("TKL"));
+ fHistCentStats->Fill(4.,aodHeader->GetCentralityP()->GetCentralityPercentile("CL0"));
+ fHistCentStats->Fill(5.,aodHeader->GetCentralityP()->GetCentralityPercentile("CL1"));
+ fHistCentStats->Fill(6.,aodHeader->GetCentralityP()->GetCentralityPercentile("V0MvsFMD"));
+ fHistCentStats->Fill(7.,aodHeader->GetCentralityP()->GetCentralityPercentile("TKLvsV0M"));
+ fHistCentStats->Fill(8.,aodHeader->GetCentralityP()->GetCentralityPercentile("ZEMvsZDC"));
+
+ // take only events inside centrality class
+ if((fCentrality < fCentralityPercentileMin) || (fCentrality > fCentralityPercentileMax))
+ return;
+
+ // centrality QA (V0M)
+ fHistV0M->Fill(gAOD->GetVZEROData()->GetMTotV0A(), gAOD->GetVZEROData()->GetMTotV0C());
+
+ // centrality QA (reference tracks)
+ fHistRefTracks->Fill(0.,aodHeader->GetRefMultiplicity());
+ fHistRefTracks->Fill(1.,aodHeader->GetRefMultiplicityPos());
+ fHistRefTracks->Fill(2.,aodHeader->GetRefMultiplicityNeg());
+ fHistRefTracks->Fill(3.,aodHeader->GetTPConlyRefMultiplicity());
+ fHistRefTracks->Fill(4.,aodHeader->GetNumberOfITSClusters(0));
+ fHistRefTracks->Fill(5.,aodHeader->GetNumberOfITSClusters(1));
+ fHistRefTracks->Fill(6.,aodHeader->GetNumberOfITSClusters(2));
+ fHistRefTracks->Fill(7.,aodHeader->GetNumberOfITSClusters(3));
+ fHistRefTracks->Fill(8.,aodHeader->GetNumberOfITSClusters(4));
+ }
+
+ const AliAODVertex *vertex = gAOD->GetPrimaryVertex();
+
+ if(vertex) {
+ Double32_t fCov[6];
+ vertex->GetCovarianceMatrix(fCov);
+
+ if(vertex->GetNContributors() > 0) {
+ if(fCov[5] != 0) {
+ fHistEventStats->Fill(3); //events with a proper vertex
+ if(TMath::Abs(vertex->GetX()) < fVxMax) {
+ if(TMath::Abs(vertex->GetY()) < fVyMax) {
+ if(TMath::Abs(vertex->GetZ()) < fVzMax) {
+ fHistEventStats->Fill(4); //analyzed events
+ fHistVx->Fill(vertex->GetX());
+ fHistVy->Fill(vertex->GetY());
+ fHistVz->Fill(vertex->GetZ());
+
+ //===========================================//
+ TExMap *trackMap = new TExMap();
+ for (Int_t iTracks = 0; iTracks < gAOD->GetNumberOfTracks(); iTracks++) {
+ AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(gAOD->GetTrack(iTracks));
+ if (!aodTrack) {
+ Printf("ERROR: Could not receive track %d", iTracks);
+ continue;
+ }
+ Int_t gID = aodTrack->GetID();
+ if (!aodTrack->TestFilterBit(nAODtrackCutBit)) trackMap->Add(gID, iTracks);
+ }
+ AliAODTrack* newAodTrack;
+ //===========================================//
+
+ //Printf("There are %d tracks in this event", gAOD->GetNumberOfTracks());
+ for (Int_t iTracks = 0; iTracks < gAOD->GetNumberOfTracks(); iTracks++) {
+ AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(gAOD->GetTrack(iTracks));
+ if (!aodTrack) {
+ Printf("ERROR: Could not receive track %d", iTracks);
+ continue;
+ }
+
+ // AOD track cuts
+
+ // For ESD Filter Information: ANALYSIS/macros/AddTaskESDfilter.C
+
+ //===========================================//
+ // take only TPC only tracks
+ fHistTrackStats->Fill(aodTrack->GetFilterMap());
+ if(!aodTrack->TestFilterBit(nAODtrackCutBit)) continue;
+ Int_t gID = aodTrack->GetID();
+ newAodTrack = gID >= 0 ? aodTrack : gAOD->GetTrack(trackMap->GetValue(-1-gID));
+ Printf("Label: %d - Pt: %lf (old) - %d - Pt: %lf(new)",gID,aodTrack->Pt(), newAodTrack->GetID(), newAodTrack->Pt());
+ //===========================================//
+
+ //fHistTrackStats->Fill(aodTrack->GetFilterMap());
+ //if(!aodTrack->TestFilterBit(nAODtrackCutBit)) continue;
+
+ v_charge = aodTrack->Charge();
+ v_y = aodTrack->Y();
+ v_eta = aodTrack->Eta();
+ v_phi = aodTrack->Phi() * TMath::RadToDeg();
+ v_E = aodTrack->E();
+ v_pt = aodTrack->Pt();
+ aodTrack->PxPyPz(v_p);
+
+ Float_t DCAxy = aodTrack->DCA(); // this is the DCA from global track (not exactly what is cut on)
+ Float_t DCAz = aodTrack->ZAtDCA(); // this is the DCA from global track (not exactly what is cut on)
+
+
+ // Kinematics cuts from ESD track cuts
+ if( v_pt < fPtMin || v_pt > fPtMax) continue;
+
+ if (!fUsePID) {
+ if( v_eta < fEtaMin || v_eta > fEtaMax) continue;
+ }
+
+ else if (fUsePID){
+ if( v_y < fEtaMin || v_y > fEtaMax) continue;
+ }
+
+ // Extra DCA cuts (for systematic studies [!= -1])
+ if( fDCAxyCut != -1 && fDCAzCut != -1){
+ if(TMath::Sqrt((DCAxy*DCAxy)/(fDCAxyCut*fDCAxyCut)+(DCAz*DCAz)/(fDCAzCut*fDCAzCut)) > 1 ){
+ continue; // 2D cut
+ }
+ }
+
+ // Extra TPC cuts (for systematic studies [!= -1])
+ if( fTPCchi2Cut != -1 && aodTrack->Chi2perNDF() > fTPCchi2Cut){
+ continue;
+ }
+ if( fNClustersTPCCut != -1 && aodTrack->GetTPCNcls() < fNClustersTPCCut){
+ continue;
+ }
+
+ //===============================================PID==================================//
+ if(fUsePID) {
+ Double_t prob[AliPID::kSPECIES]={0.};
+ Double_t probTPC[AliPID::kSPECIES]={0.};
+ Double_t probTOF[AliPID::kSPECIES]={0.};
+ Double_t probTPCTOF[AliPID::kSPECIES]={0.};
+
+ Double_t nSigma = 0.;
+ UInt_t detUsedTPC = 0;
+ UInt_t detUsedTOF = 0;
+ UInt_t detUsedTPCTOF = 0;
+
+ //Decide what detector configuration we want to use
+ switch(fPidDetectorConfig) {
+ case kTPCpid:
+ fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
+ nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(newAodTrack,(AliPID::EParticleType)fParticleOfInterest));
+ //detUsedTPC = fPIDCombined->ComputeProbabilities(aodTrack, fPIDResponse, probTPC);
+ detUsedTPC = (AliPIDResponse::kDetTPC);
+ for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
+ prob[iSpecies] = probTPC[iSpecies];
+ break;
+ case kTOFpid:
+ fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF);
+ nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(newAodTrack,(AliPID::EParticleType)fParticleOfInterest));
+ //detUsedTOF = fPIDCombined->ComputeProbabilities(aodTrack, fPIDResponse, probTOF);
+ detUsedTPC = (AliPIDResponse::kDetTPC);
+ for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
+ prob[iSpecies] = probTOF[iSpecies];
+ break;
+ case kTPCTOF:
+ fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC);
+ //detUsedTPCTOF = fPIDCombined->ComputeProbabilities(newAodTrack, fPIDResponse, probTPCTOF);
+ detUsedTPC = (AliPIDResponse::kDetTPC);
+ for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
+ prob[iSpecies] = probTPCTOF[iSpecies];
+ break;
+ default:
+ break;
+ }//end switch: define detector mask
+
+ //Filling the PID QA
+ Double_t tofTime = -999., tof = -999.; //length = 999., tof = -999.;
+ //Double_t c = TMath::C()*1.E-9;// m/ns
+ //Double_t beta = -999.;
+ Double_t nSigmaTOFForParticleOfInterest = -999.;
+ if ( (newAodTrack->IsOn(AliAODTrack::kTOFin)) &&
+ (newAodTrack->IsOn(AliAODTrack::kTIME)) ) {
+ tofTime = newAodTrack->GetTOFsignal();//in ps
+ //length = newAodTrack->GetIntegratedLength();
+ tof = tofTime*1E-3; // ns
+
+ if (tof <= 0) {
+ //Printf("WARNING: track with negative TOF time found! Skipping this track for PID checks\n");
+ continue;
+ }
+ //if (length <= 0){
+ //printf("WARNING: track with negative length found!Skipping this track for PID checks\n");
+ //continue;
+ //}
+
+ //length = length*0.01; // in meters
+ //tof = tof*c;
+ //beta = length/tof;
+
+ nSigmaTOFForParticleOfInterest = fPIDResponse->NumberOfSigmasTOF(newAodTrack,(AliPID::EParticleType)fParticleOfInterest);
+ //fHistBetavsPTOFbeforePID ->Fill(aodTrack->P()*aodTrack->Charge(),beta);
+ fHistProbTOFvsPtbeforePID ->Fill(newAodTrack->Pt(),probTOF[fParticleOfInterest]);
+ fHistNSigmaTOFvsPtbeforePID ->Fill(newAodTrack->Pt(),nSigmaTOFForParticleOfInterest);
+ }//TOF signal
+
+
+ Double_t nSigmaTPCForParticleOfInterest = fPIDResponse->NumberOfSigmasTPC(newAodTrack,(AliPID::EParticleType)fParticleOfInterest);
+ fHistdEdxVsPTPCbeforePID -> Fill(newAodTrack->P()*newAodTrack->Charge(),newAodTrack->GetTPCsignal());
+ fHistProbTPCvsPtbeforePID -> Fill(newAodTrack->Pt(),probTPC[fParticleOfInterest]);
+ fHistNSigmaTPCvsPtbeforePID -> Fill(newAodTrack->Pt(),nSigmaTPCForParticleOfInterest);
+ fHistProbTPCTOFvsPtbeforePID -> Fill(newAodTrack->Pt(),probTPCTOF[fParticleOfInterest]);
+ //end of QA-before pid
+
+ if ((detUsedTPC != 0)||(detUsedTOF != 0)||(detUsedTPCTOF != 0)) {
+ //Make the decision based on the n-sigma
+ if(fUsePIDnSigma) {
+ if(nSigma > fPIDNSigma) continue;}
+
+ //Make the decision based on the bayesian
+ else if(fUsePIDPropabilities) {
+ if(fParticleOfInterest != TMath::LocMax(AliPID::kSPECIES,prob)) continue;
+ if (prob[fParticleOfInterest] < fMinAcceptedPIDProbability) continue;
+ }
+
+ //Fill QA after the PID
+ //fHistBetavsPTOFafterPID ->Fill(newAodTrack->P()*newAodTrack->Charge(),beta);
+ fHistProbTOFvsPtafterPID ->Fill(newAodTrack->Pt(),probTOF[fParticleOfInterest]);
+ fHistNSigmaTOFvsPtafterPID ->Fill(newAodTrack->Pt(),nSigmaTOFForParticleOfInterest);
+
+ fHistdEdxVsPTPCafterPID -> Fill(newAodTrack->P()*newAodTrack->Charge(),newAodTrack->GetTPCsignal());
+ fHistProbTPCvsPtafterPID -> Fill(newAodTrack->Pt(),probTPC[fParticleOfInterest]);
+ fHistProbTPCTOFvsPtafterPID -> Fill(newAodTrack->Pt(),probTPCTOF[fParticleOfInterest]);
+ fHistNSigmaTPCvsPtafterPID -> Fill(newAodTrack->Pt(),nSigmaTPCForParticleOfInterest);
+ }
+
+ PostData(4, fHistListPIDQA);
+ }
+
+ //=========================================================PID=================================================================//
+
+ // fill QA histograms
+ fHistClus->Fill(aodTrack->GetITSNcls(),aodTrack->GetTPCNcls());
+ fHistDCA->Fill(DCAz,DCAxy);
+ fHistChi2->Fill(aodTrack->Chi2perNDF());
+ fHistPt->Fill(v_pt);
+ fHistEta->Fill(v_eta);
+ fHistPhi->Fill(v_phi);
+ fHistRapidity->Fill(v_y);
+ if(v_charge > 0) fHistPhiPos->Fill(v_phi);
+ else if(v_charge < 0) fHistPhiNeg->Fill(v_phi);
+
+ // fill charge vector
+ chargeVector[0]->push_back(v_charge);
+ chargeVector[1]->push_back(v_y);
+ chargeVector[2]->push_back(v_eta);
+ chargeVector[3]->push_back(v_phi);
+ chargeVector[4]->push_back(v_p[0]);
+ chargeVector[5]->push_back(v_p[1]);
+ chargeVector[6]->push_back(v_p[2]);
+ chargeVector[7]->push_back(v_pt);
+ chargeVector[8]->push_back(v_E);
+
+ if(fRunShuffling) {
+ chargeVectorShuffle[0]->push_back(v_charge);
+ chargeVectorShuffle[1]->push_back(v_y);
+ chargeVectorShuffle[2]->push_back(v_eta);
+ chargeVectorShuffle[3]->push_back(v_phi);
+ chargeVectorShuffle[4]->push_back(v_p[0]);
+ chargeVectorShuffle[5]->push_back(v_p[1]);
+ chargeVectorShuffle[6]->push_back(v_p[2]);
+ chargeVectorShuffle[7]->push_back(v_pt);
+ chargeVectorShuffle[8]->push_back(v_E);
+ }
+
+ gNumberOfAcceptedTracks += 1;
+
+ } //track loop
+ }//Vz cut
+ }//Vy cut
+ }//Vx cut
+ }//proper vertex resolution
+ }//proper number of contributors
+ }//vertex object valid
+ }//triggered event
+ }//AOD analysis
+
+ //MC-ESD analysis
+ if(gAnalysisLevel == "MCESD") {
+ AliMCEvent* mcEvent = MCEvent();
+ if (!mcEvent) {
+ Printf("ERROR: mcEvent not available");
+ return;
+ }
+
+ AliESDEvent* gESD = dynamic_cast<AliESDEvent*>(InputEvent()); // from TaskSE
+ if (!gESD) {
+ Printf("ERROR: gESD not available");
+ return;
+ }
+
+ // store offline trigger bits
+ fHistTriggerStats->Fill(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected());
+
+ // event selection done in AliAnalysisTaskSE::Exec() --> this is not used
+ fHistEventStats->Fill(1); //all events
+ Bool_t isSelected = kTRUE;
+ if(fUseOfflineTrigger)
+ isSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
+ if(isSelected) {
+ fHistEventStats->Fill(2); //triggered events
+
+ if(fUseCentrality) {
+ //Centrality stuff
+ AliCentrality *centrality = gESD->GetCentrality();
+
+ fCentrality = centrality->GetCentralityPercentile(fCentralityEstimator.Data());
+
+ // take only events inside centrality class
+ if(!centrality->IsEventInCentralityClass(fCentralityPercentileMin,
+ fCentralityPercentileMax,
+ fCentralityEstimator.Data()))
+ return;
+
+ // centrality QA (V0M)
+ fHistV0M->Fill(gESD->GetVZEROData()->GetMTotV0A(), gESD->GetVZEROData()->GetMTotV0C());
+ }
+
+ const AliESDVertex *vertex = gESD->GetPrimaryVertex();
+ if(vertex) {
+ if(vertex->GetNContributors() > 0) {
+ if(vertex->GetZRes() != 0) {
+ fHistEventStats->Fill(3); //events with a proper vertex
+ if(TMath::Abs(vertex->GetXv()) < fVxMax) {
+ if(TMath::Abs(vertex->GetYv()) < fVyMax) {
+ if(TMath::Abs(vertex->GetZv()) < fVzMax) {
+ fHistEventStats->Fill(4); //analayzed events
+ fHistVx->Fill(vertex->GetXv());
+ fHistVy->Fill(vertex->GetYv());
+ fHistVz->Fill(vertex->GetZv());
+
+ //Printf("There are %d tracks in this event", gESD->GetNumberOfTracks());
+ for (Int_t iTracks = 0; iTracks < gESD->GetNumberOfTracks(); iTracks++) {
+ AliESDtrack* track = dynamic_cast<AliESDtrack *>(gESD->GetTrack(iTracks));
+ if (!track) {
+ Printf("ERROR: Could not receive track %d", iTracks);
+ continue;
+ }
+
+ Int_t label = TMath::Abs(track->GetLabel());
+ if(label > mcEvent->GetNumberOfTracks()) continue;
+ if(label > mcEvent->GetNumberOfPrimaries()) continue;
+
+ AliMCParticle* mcTrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(label));
+ if(!mcTrack) continue;
+
+ // take only TPC only tracks
+ track_TPC = new AliESDtrack();
+ if(!track->FillTPCOnlyTrack(*track_TPC)) continue;
+
+ //ESD track cuts
+ if(fESDtrackCuts)
+ if(!fESDtrackCuts->AcceptTrack(track_TPC)) continue;
+
+ // fill QA histograms
+ Float_t b[2];
+ Float_t bCov[3];
+ track_TPC->GetImpactParameters(b,bCov);
+ if (bCov[0]<=0 || bCov[2]<=0) {
+ AliDebug(1, "Estimated b resolution lower or equal zero!");
+ bCov[0]=0; bCov[2]=0;
+ }
+
+ Int_t nClustersTPC = -1;
+ nClustersTPC = track_TPC->GetTPCNclsIter1(); // TPC standalone
+ //nClustersTPC = track->GetTPCclusters(0); // global track
+ Float_t chi2PerClusterTPC = -1;
+ if (nClustersTPC!=0) {
+ chi2PerClusterTPC = track_TPC->GetTPCchi2Iter1()/Float_t(nClustersTPC); // TPC standalone
+ //chi2PerClusterTPC = track->GetTPCchi2()/Float_t(nClustersTPC); // global track
+ }
+
+ v_charge = track_TPC->Charge();
+ v_y = track_TPC->Y();
+ v_eta = track_TPC->Eta();
+ v_phi = track_TPC->Phi() * TMath::RadToDeg();
+ v_E = track_TPC->E();
+ v_pt = track_TPC->Pt();
+ track_TPC->PxPyPz(v_p);
+
+ fHistClus->Fill(track_TPC->GetITSclusters(0),nClustersTPC);
+ fHistDCA->Fill(b[1],b[0]);
+ fHistChi2->Fill(chi2PerClusterTPC);
+ fHistPt->Fill(v_pt);
+ fHistEta->Fill(v_eta);
+ fHistPhi->Fill(v_phi);
+ fHistRapidity->Fill(v_y);
+ if(v_charge > 0) fHistPhiPos->Fill(v_phi);
+ else if(v_charge < 0) fHistPhiNeg->Fill(v_phi);
+
+ // fill charge vector
+ chargeVector[0]->push_back(v_charge);
+ chargeVector[1]->push_back(v_y);
+ chargeVector[2]->push_back(v_eta);
+ chargeVector[3]->push_back(v_phi);
+ chargeVector[4]->push_back(v_p[0]);
+ chargeVector[5]->push_back(v_p[1]);
+ chargeVector[6]->push_back(v_p[2]);
+ chargeVector[7]->push_back(v_pt);
+ chargeVector[8]->push_back(v_E);
+
+ if(fRunShuffling) {
+ chargeVectorShuffle[0]->push_back(v_charge);
+ chargeVectorShuffle[1]->push_back(v_y);
+ chargeVectorShuffle[2]->push_back(v_eta);
+ chargeVectorShuffle[3]->push_back(v_phi);
+ chargeVectorShuffle[4]->push_back(v_p[0]);
+ chargeVectorShuffle[5]->push_back(v_p[1]);
+ chargeVectorShuffle[6]->push_back(v_p[2]);
+ chargeVectorShuffle[7]->push_back(v_pt);
+ chargeVectorShuffle[8]->push_back(v_E);
+ }
+
+ delete track_TPC;
+ gNumberOfAcceptedTracks += 1;
+
+ } //track loop
+ }//Vz cut
+ }//Vy cut
+ }//Vx cut
+ }//proper vertex resolution
+ }//proper number of contributors
+ }//vertex object valid
+ }//triggered event
+ }//MC-ESD analysis
+
+ //MC analysis
+ else if(gAnalysisLevel == "MC") {
+ AliMCEvent* mcEvent = MCEvent();
+ if (!mcEvent) {
+ Printf("ERROR: mcEvent not available");
+ return;
+ }
+ fHistEventStats->Fill(1); //total events
+ fHistEventStats->Fill(2); //offline trigger
+
+ Double_t gReactionPlane = 0., gImpactParameter = 0.;
+ if(fUseCentrality) {
+ //Get the MC header
+ AliGenHijingEventHeader* headerH = dynamic_cast<AliGenHijingEventHeader*>(mcEvent->GenEventHeader());
+ if (headerH) {
+ //Printf("=====================================================");
+ //Printf("Reaction plane angle: %lf",headerH->ReactionPlaneAngle());
+ //Printf("Impact parameter: %lf",headerH->ImpactParameter());
+ //Printf("=====================================================");
+ gReactionPlane = headerH->ReactionPlaneAngle();
+ gImpactParameter = headerH->ImpactParameter();
+ fCentrality = gImpactParameter;
+ }
+ fCentrality = gImpactParameter;
+
+ // take only events inside centrality class (DIDN'T CHANGE THIS UP TO NOW)
+ if((fImpactParameterMin > gImpactParameter) || (fImpactParameterMax < gImpactParameter))
+ return;
+ }
+
+ AliGenEventHeader *header = mcEvent->GenEventHeader();
+ if(!header) return;
+
+ TArrayF gVertexArray;
+ header->PrimaryVertex(gVertexArray);
+ //Printf("Vertex: %lf (x) - %lf (y) - %lf (z)",
+ //gVertexArray.At(0),
+ //gVertexArray.At(1),
+ //gVertexArray.At(2));
+ fHistEventStats->Fill(3); //events with a proper vertex
+ if(TMath::Abs(gVertexArray.At(0)) < fVxMax) {
+ if(TMath::Abs(gVertexArray.At(1)) < fVyMax) {
+ if(TMath::Abs(gVertexArray.At(2)) < fVzMax) {
+ fHistEventStats->Fill(4); //analayzed events
+ fHistVx->Fill(gVertexArray.At(0));
+ fHistVy->Fill(gVertexArray.At(1));
+ fHistVz->Fill(gVertexArray.At(2));
+
+ Printf("There are %d tracks in this event", mcEvent->GetNumberOfPrimaries());
+ for (Int_t iTracks = 0; iTracks < mcEvent->GetNumberOfPrimaries(); iTracks++) {
+ AliMCParticle* track = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(iTracks));
+ if (!track) {
+ Printf("ERROR: Could not receive particle %d", iTracks);
+ continue;
+ }
+
+ //exclude non stable particles
+ if(!mcEvent->IsPhysicalPrimary(iTracks)) continue;
+
+ v_eta = track->Eta();
+ v_pt = track->Pt();
+ v_y = track->Y();
+
+ if( v_pt < fPtMin || v_pt > fPtMax)
+ continue;
+ if (!fUsePID) {
+ if( v_eta < fEtaMin || v_eta > fEtaMax) continue;
+ }
+ else if (fUsePID){
+ if( v_y < fEtaMin || v_y > fEtaMax) continue;
+ }
+
+ //analyze one set of particles
+ if(fUseMCPdgCode) {
+ TParticle *particle = track->Particle();
+ if(!particle) continue;
+
+ Int_t gPdgCode = particle->GetPdgCode();
+ if(TMath::Abs(fPDGCodeToBeAnalyzed) != TMath::Abs(gPdgCode))
+ continue;
+ }
+
+ //Use the acceptance parameterization
+ if(fAcceptanceParameterization) {
+ Double_t gRandomNumber = gRandom->Rndm();
+ if(gRandomNumber > fAcceptanceParameterization->Eval(track->Pt()))
+ continue;
+ }
+
+ //Exclude resonances
+ if(fExcludeResonancesInMC) {
+ TParticle *particle = track->Particle();
+ if(!particle) continue;
+
+ Bool_t kExcludeParticle = kFALSE;
+ Int_t gMotherIndex = particle->GetFirstMother();
+ if(gMotherIndex != -1) {
+ AliMCParticle* motherTrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(gMotherIndex));
+ if(motherTrack) {
+ TParticle *motherParticle = motherTrack->Particle();
+ if(motherParticle) {
+ Int_t pdgCodeOfMother = motherParticle->GetPdgCode();
+ //if((pdgCodeOfMother == 113)||(pdgCodeOfMother == 213)||(pdgCodeOfMother == 221)||(pdgCodeOfMother == 223)||(pdgCodeOfMother == 331)||(pdgCodeOfMother == 333)) {
+ if(pdgCodeOfMother == 113) {
+ kExcludeParticle = kTRUE;
+ }
+ }
+ }
+ }
+
+ //Exclude from the analysis decay products of rho0, rho+, eta, eta' and phi
+ if(kExcludeParticle) continue;
+ }
+
+ v_charge = track->Charge();
+ v_phi = track->Phi();
+ v_E = track->E();
+ track->PxPyPz(v_p);
+ //Printf("phi (before): %lf",v_phi);
+
+ fHistPt->Fill(v_pt);
+ fHistEta->Fill(v_eta);
+ fHistPhi->Fill(v_phi);
+ fHistRapidity->Fill(v_y);
+ if(v_charge > 0) fHistPhiPos->Fill(v_phi);
+ else if(v_charge < 0) fHistPhiNeg->Fill(v_phi);
+
+ //Flow after burner
+ if(fUseFlowAfterBurner) {
+ Double_t precisionPhi = 0.001;
+ Int_t maxNumberOfIterations = 100;
+
+ Double_t phi0 = v_phi;
+ Double_t gV2 = fDifferentialV2->Eval(v_pt);
+
+ for (Int_t j = 0; j < maxNumberOfIterations; j++) {
+ Double_t phiprev = v_phi;
+ Double_t fl = v_phi - phi0 + gV2*TMath::Sin(2.*(v_phi - gReactionPlane));
+ Double_t fp = 1.0 + 2.0*gV2*TMath::Cos(2.*(v_phi - gReactionPlane));
+ v_phi -= fl/fp;
+ if (TMath::AreEqualAbs(phiprev,v_phi,precisionPhi)) break;
+ }
+ //Printf("phi (after): %lf\n",v_phi);
+ Double_t v_DeltaphiBefore = phi0 - gReactionPlane;
+ if(v_DeltaphiBefore < 0) v_DeltaphiBefore += 2*TMath::Pi();
+ fHistPhiBefore->Fill(v_DeltaphiBefore);
+
+ Double_t v_DeltaphiAfter = v_phi - gReactionPlane;
+ if(v_DeltaphiAfter < 0) v_DeltaphiAfter += 2*TMath::Pi();
+ fHistPhiAfter->Fill(v_DeltaphiAfter);
+ }
+
+ v_phi *= TMath::RadToDeg();
+
+ // fill charge vector
+ chargeVector[0]->push_back(v_charge);
+ chargeVector[1]->push_back(v_y);
+ chargeVector[2]->push_back(v_eta);
+ chargeVector[3]->push_back(v_phi);
+ chargeVector[4]->push_back(v_p[0]);
+ chargeVector[5]->push_back(v_p[1]);
+ chargeVector[6]->push_back(v_p[2]);
+ chargeVector[7]->push_back(v_pt);
+ chargeVector[8]->push_back(v_E);
+
+ if(fRunShuffling) {
+ chargeVectorShuffle[0]->push_back(v_charge);
+ chargeVectorShuffle[1]->push_back(v_y);
+ chargeVectorShuffle[2]->push_back(v_eta);
+ chargeVectorShuffle[3]->push_back(v_phi);
+ chargeVectorShuffle[4]->push_back(v_p[0]);
+ chargeVectorShuffle[5]->push_back(v_p[1]);
+ chargeVectorShuffle[6]->push_back(v_p[2]);
+ chargeVectorShuffle[7]->push_back(v_pt);
+ chargeVectorShuffle[8]->push_back(v_E);
+ }
+ gNumberOfAcceptedTracks += 1;
+
+ } //track loop
+ }//Vz cut
+ }//Vy cut
+ }//Vx cut
+ }//MC analysis
+
+ //multiplicity cut (used in pp)
+ if(fUseMultiplicity) {
+ if((gNumberOfAcceptedTracks < fNumberOfAcceptedTracksMin)||(gNumberOfAcceptedTracks > fNumberOfAcceptedTracksMax))
+ return;
+ }
+ fHistNumberOfAcceptedTracks->Fill(gNumberOfAcceptedTracks);
+
+ // calculate balance function
+ if(fUseMultiplicity)
+ fBalance->CalculateBalance(gNumberOfAcceptedTracks,chargeVector,bSign);
+ else
+ fBalance->CalculateBalance(fCentrality,chargeVector,bSign);
+
+ if(fRunShuffling) {
+ // shuffle charges
+ random_shuffle( chargeVectorShuffle[0]->begin(), chargeVectorShuffle[0]->end() );
+ if(fUseMultiplicity)
+ fShuffledBalance->CalculateBalance(gNumberOfAcceptedTracks,chargeVectorShuffle,bSign);
+ else
+ fShuffledBalance->CalculateBalance(fCentrality,chargeVectorShuffle,bSign);
+ }
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskBF::FinishTaskOutput(){
+ //Printf("END BF");
+
+ if (!fBalance) {
+ Printf("ERROR: fBalance not available");
+ return;
+ }
+ if(fRunShuffling) {
+ if (!fShuffledBalance) {
+ Printf("ERROR: fShuffledBalance not available");
+ return;
+ }
+ }
+
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskBF::Terminate(Option_t *) {
+ // Draw result to the screen
+ // Called once at the end of the query
+
+ // not implemented ...
+
+}
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff