-#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 "AliMixInputEventHandler.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 "AliAnalysisTaskEventMixingBF.h"\r
-#include "AliBalanceEventMixing.h"\r
-\r
-\r
-// Analysis task for the EventMixingBF code\r
-// Authors: Panos.Christakoglou@nikhef.nl, m.weber@cern.ch\r
-\r
-ClassImp(AliAnalysisTaskEventMixingBF)\r
-\r
-//________________________________________________________________________\r
-AliAnalysisTaskEventMixingBF::AliAnalysisTaskEventMixingBF(const char *name) \r
-: AliAnalysisTaskSE(name), \r
- fBalance(0),\r
- fRunShuffling(kFALSE),\r
- fShuffledBalance(0),\r
- fList(0),\r
- fListEventMixingBF(0),\r
- fListEventMixingBFS(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
- fHistPhi(0),\r
- fHistPhiBefore(0),\r
- fHistPhiAfter(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
- fMainEvent(0x0),\r
- fMixEvent(0x0)\r
-{\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
-AliAnalysisTaskEventMixingBF::~AliAnalysisTaskEventMixingBF() {\r
-\r
- // delete fBalance; \r
- // delete fShuffledBalance; \r
- // delete fList;\r
- // delete fListEventMixingBF; \r
- // delete fListEventMixingBFS;\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 AliAnalysisTaskEventMixingBF::UserCreateOutputObjects() {\r
- // Create histograms\r
- // Called once\r
- if(!fBalance) {\r
- fBalance = new AliBalanceEventMixing();\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 AliBalanceEventMixing();\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
- fListEventMixingBF = new TList();\r
- fListEventMixingBF->SetName("listEventMixingBF");\r
- fListEventMixingBF->SetOwner();\r
-\r
- if(fRunShuffling) {\r
- fListEventMixingBFS = new TList();\r
- fListEventMixingBFS->SetName("listEventMixingBFShuffled");\r
- fListEventMixingBFS->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
- 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
- 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
- fListEventMixingBF->Add(fBalance->GetHistNp(a));\r
- fListEventMixingBF->Add(fBalance->GetHistNn(a));\r
- fListEventMixingBF->Add(fBalance->GetHistNpn(a));\r
- fListEventMixingBF->Add(fBalance->GetHistNnn(a));\r
- fListEventMixingBF->Add(fBalance->GetHistNpp(a));\r
- fListEventMixingBF->Add(fBalance->GetHistNnp(a));\r
-\r
- if(fRunShuffling) {\r
- fListEventMixingBFS->Add(fShuffledBalance->GetHistNp(a));\r
- fListEventMixingBFS->Add(fShuffledBalance->GetHistNn(a));\r
- fListEventMixingBFS->Add(fShuffledBalance->GetHistNpn(a));\r
- fListEventMixingBFS->Add(fShuffledBalance->GetHistNnn(a));\r
- fListEventMixingBFS->Add(fShuffledBalance->GetHistNpp(a));\r
- fListEventMixingBFS->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, fListEventMixingBF);\r
- if(fRunShuffling) PostData(3, fListEventMixingBFS);\r
- if(fUsePID) PostData(4, fHistListPIDQA); //PID\r
-}\r
-\r
-//________________________________________________________________________\r
-void AliAnalysisTaskEventMixingBF::UserExec(Option_t *) {\r
- // Main loop\r
- // Called for each event\r
- // NOTHING TO DO for event mixing!\r
-} \r
-\r
-//________________________________________________________________________\r
-void AliAnalysisTaskEventMixingBF::FinishTaskOutput(){\r
- //Printf("END EventMixingBF");\r
-\r
- if (!fBalance) {\r
- AliError("ERROR: fBalance not available");\r
- return;\r
- } \r
- if(fRunShuffling) {\r
- if (!fShuffledBalance) {\r
- AliError("ERROR: fShuffledBalance not available");\r
- return;\r
- }\r
- }\r
-\r
-}\r
-\r
-//________________________________________________________________________\r
-void AliAnalysisTaskEventMixingBF::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
-\r
-void AliAnalysisTaskEventMixingBF::UserExecMix(Option_t *)\r
-{\r
- // Main loop for event mixing\r
-\r
- TString gAnalysisLevel = fBalance->GetAnalysisLevel();\r
-\r
- AliMixInputEventHandler *mixIEH = SetupEventsForMixing();\r
-\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> *chargeVector[9]; // original charge\r
- for(Int_t i = 0; i < 9; i++){\r
- chargeVector[i] = new vector<Double_t>;\r
- }\r
- \r
- Double_t vCharge;\r
- Double_t vY;\r
- Double_t vEta;\r
- Double_t vPhi;\r
- Double_t vP[3];\r
- Double_t vPt;\r
- Double_t vE;\r
-\r
- Int_t iMainTrackUsed = -1;\r
-\r
- // ------------------------------------------------------------- \r
- // At the moment MIXING only for AODs\r
- if(mixIEH){\r
-\r
- //AOD analysis (vertex and track cuts also here!!!!)\r
- if(gAnalysisLevel == "AOD") {\r
- AliAODEvent* aodEventMain = dynamic_cast<AliAODEvent*>(fMainEvent); \r
- if(!aodEventMain) {\r
- AliError("ERROR: aodEventMain not available");\r
- return;\r
- }\r
- AliAODEvent *aodEventMix = dynamic_cast<AliAODEvent *>(fMixEvent); \r
- if(!aodEventMix) {\r
- AliError("ERROR: aodEventMix not available");\r
- return;\r
- }\r
-\r
- // for HBT like cuts need magnetic field sign\r
- bSign = (aodEventMain->GetMagneticField() > 0) ? 1 : -1;\r
- \r
- AliAODHeader *aodHeaderMain = aodEventMain->GetHeader(); \r
-\r
- // event selection done in AliAnalysisTaskSE::Exec() --> this is not used\r
- fHistEventStats->Fill(1); //all events\r
-\r
- // this is not needed (checked in mixing handler!)\r
- Bool_t isSelectedMain = kTRUE;\r
- Bool_t isSelectedMix = kTRUE;\r
- \r
- if(fUseOfflineTrigger){\r
- isSelectedMain = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();\r
- isSelectedMix = ((AliInputEventHandler*)((AliMultiInputEventHandler *)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetFirstMultiInputHandler())->IsEventSelected();\r
- }\r
- \r
- if(isSelectedMain && isSelectedMix) {\r
- fHistEventStats->Fill(2); //triggered events\r
- \r
- //Centrality stuff (centrality in AOD header)\r
- if(fUseCentrality) {\r
- fCentrality = aodHeaderMain->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());\r
- \r
- // QA for centrality estimators\r
- fHistCentStats->Fill(0.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("V0M"));\r
- fHistCentStats->Fill(1.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("FMD"));\r
- fHistCentStats->Fill(2.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("TRK"));\r
- fHistCentStats->Fill(3.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("TKL"));\r
- fHistCentStats->Fill(4.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("CL0"));\r
- fHistCentStats->Fill(5.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("CL1"));\r
- fHistCentStats->Fill(6.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("V0MvsFMD"));\r
- fHistCentStats->Fill(7.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("TKLvsV0M"));\r
- fHistCentStats->Fill(8.,aodHeaderMain->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(aodEventMain->GetVZEROData()->GetMTotV0A(), aodEventMain->GetVZEROData()->GetMTotV0C());\r
- \r
- // centrality QA (reference tracks)\r
- fHistRefTracks->Fill(0.,aodHeaderMain->GetRefMultiplicity());\r
- fHistRefTracks->Fill(1.,aodHeaderMain->GetRefMultiplicityPos());\r
- fHistRefTracks->Fill(2.,aodHeaderMain->GetRefMultiplicityNeg());\r
- fHistRefTracks->Fill(3.,aodHeaderMain->GetTPConlyRefMultiplicity());\r
- fHistRefTracks->Fill(4.,aodHeaderMain->GetNumberOfITSClusters(0));\r
- fHistRefTracks->Fill(5.,aodHeaderMain->GetNumberOfITSClusters(1));\r
- fHistRefTracks->Fill(6.,aodHeaderMain->GetNumberOfITSClusters(2));\r
- fHistRefTracks->Fill(7.,aodHeaderMain->GetNumberOfITSClusters(3));\r
- fHistRefTracks->Fill(8.,aodHeaderMain->GetNumberOfITSClusters(4));\r
- }\r
- \r
- // // this is crashing (Bug in ROOT (to be solved)) but not needed (checked in mixing handler)\r
- // const AliAODVertex *vertexMain = aodEventMain->GetPrimaryVertex();\r
- // const AliAODVertex *vertexMix = aodEventMix->GetPrimaryVertex();\r
- \r
- // if(vertexMain && vertexMix) {\r
- // Double32_t fCovMain[6];\r
- // Double32_t fCovMix[6];\r
- // vertexMain->GetCovarianceMatrix(fCovMain);\r
- // vertexMix->GetCovarianceMatrix(fCovMix);\r
- \r
- // if(vertexMain->GetNContributors() > 0 && vertexMix->GetNContributors() > 0) {\r
- // if(fCovMain[5] != 0 && fCovMix[5] != 0) {\r
- // fHistEventStats->Fill(3); //events with a proper vertex\r
- // if(TMath::Abs(vertexMain->GetX()) < fVxMax && TMath::Abs(vertexMix->GetX()) < fVxMax ) {\r
- // if(TMath::Abs(vertexMain->GetY()) < fVyMax && TMath::Abs(vertexMix->GetY()) < fVyMax) {\r
- // if(TMath::Abs(vertexMain->GetZ()) < fVzMax && TMath::Abs(vertexMix->GetZ()) < fVzMax) {\r
- // fHistEventStats->Fill(4); //analyzed events\r
- // fHistVx->Fill(vertexMain->GetX());\r
- // fHistVy->Fill(vertexMain->GetY());\r
- // fHistVz->Fill(vertexMain->GetZ());\r
-\r
- // Loop over tracks in main event\r
- for (Int_t iTracksMain = 0; iTracksMain < aodEventMain->GetNumberOfTracks(); iTracksMain++) {\r
- AliAODTrack* aodTrackMain = dynamic_cast<AliAODTrack *>(aodEventMain->GetTrack(iTracksMain));\r
- if (!aodTrackMain) {\r
- AliError(Form("ERROR: Could not receive track %d", iTracksMain));\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(aodTrackMain->GetFilterMap());\r
- if(!aodTrackMain->TestFilterBit(nAODtrackCutBit)) continue;\r
- \r
- vCharge = aodTrackMain->Charge();\r
- vY = aodTrackMain->Y();\r
- vEta = aodTrackMain->Eta();\r
- vPhi = aodTrackMain->Phi() * TMath::RadToDeg();\r
- vE = aodTrackMain->E();\r
- vPt = aodTrackMain->Pt();\r
- aodTrackMain->PxPyPz(vP);\r
- \r
- Float_t dcaXYMain = aodTrackMain->DCA(); // this is the DCA from global track (not exactly what is cut on)\r
- Float_t dcaZMain = aodTrackMain->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((dcaXYMain*dcaXYMain)/(fDCAxyCut*fDCAxyCut)+(dcaZMain*dcaZMain)/(fDCAzCut*fDCAzCut)) > 1 ){\r
- continue; // 2D cut\r
- }\r
- }\r
- \r
- // Extra TPC cuts (for systematic studies [!= -1])\r
- if( fTPCchi2Cut != -1 && aodTrackMain->Chi2perNDF() > fTPCchi2Cut){\r
- continue;\r
- }\r
- if( fNClustersTPCCut != -1 && aodTrackMain->GetTPCNcls() < fNClustersTPCCut){\r
- continue;\r
- }\r
- \r
- // fill QA histograms\r
- fHistClus->Fill(aodTrackMain->GetITSNcls(),aodTrackMain->GetTPCNcls());\r
- fHistDCA->Fill(dcaZMain,dcaXYMain);\r
- fHistChi2->Fill(aodTrackMain->Chi2perNDF());\r
- fHistPt->Fill(vPt);\r
- fHistEta->Fill(vEta);\r
- fHistPhi->Fill(vPhi);\r
- \r
- // fill charge vector\r
- chargeVector[0]->push_back(vCharge);\r
- chargeVector[1]->push_back(vY);\r
- chargeVector[2]->push_back(vEta);\r
- chargeVector[3]->push_back(vPhi);\r
- chargeVector[4]->push_back(vP[0]);\r
- chargeVector[5]->push_back(vP[1]);\r
- chargeVector[6]->push_back(vP[2]);\r
- chargeVector[7]->push_back(vPt);\r
- chargeVector[8]->push_back(vE);\r
-\r
- // ------------------------------------------------------------- \r
- // for each track in main event loop over all tracks in mix event\r
- for (Int_t iTracksMix = 0; iTracksMix < aodEventMix->GetNumberOfTracks(); iTracksMix++) {\r
-\r
- AliAODTrack* aodTrackMix = dynamic_cast<AliAODTrack *>(aodEventMix->GetTrack(iTracksMix));\r
- if (!aodTrackMix) {\r
- AliError(Form("ERROR: Could not receive track %d", iTracksMix));\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(aodTrackMix->GetFilterMap());\r
- if(!aodTrackMix->TestFilterBit(nAODtrackCutBit)) continue;\r
- \r
- vCharge = aodTrackMix->Charge();\r
- vY = aodTrackMix->Y();\r
- vEta = aodTrackMix->Eta();\r
- vPhi = aodTrackMix->Phi() * TMath::RadToDeg();\r
- vE = aodTrackMix->E();\r
- vPt = aodTrackMix->Pt();\r
- aodTrackMix->PxPyPz(vP);\r
- \r
- Float_t dcaXYMix = aodTrackMix->DCA(); // this is the DCA from global track (not exactly what is cut on)\r
- Float_t dcaZMix = aodTrackMix->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 && fDCAxyCut != -1){\r
- if(TMath::Sqrt((dcaXYMix*dcaXYMix)/(fDCAxyCut*fDCAxyCut)+(dcaZMix*dcaZMix)/(fDCAzCut*fDCAzCut)) > 1 ){\r
- continue; // 2D cut\r
- }\r
- }\r
- \r
- // Extra TPC cuts (for systematic studies [!= -1])\r
- if( fTPCchi2Cut != -1 && aodTrackMix->Chi2perNDF() > fTPCchi2Cut){\r
- continue;\r
- }\r
- if( fNClustersTPCCut != -1 && aodTrackMix->GetTPCNcls() < fNClustersTPCCut){\r
- continue;\r
- }\r
- \r
- // fill QA histograms\r
- fHistClus->Fill(aodTrackMix->GetITSNcls(),aodTrackMix->GetTPCNcls());\r
- fHistDCA->Fill(dcaZMix,dcaXYMix);\r
- fHistChi2->Fill(aodTrackMix->Chi2perNDF());\r
- fHistPt->Fill(vPt);\r
- fHistEta->Fill(vEta);\r
- fHistPhi->Fill(vPhi);\r
- \r
- // fill charge vector\r
- chargeVector[0]->push_back(vCharge);\r
- chargeVector[1]->push_back(vY);\r
- chargeVector[2]->push_back(vEta);\r
- chargeVector[3]->push_back(vPhi);\r
- chargeVector[4]->push_back(vP[0]);\r
- chargeVector[5]->push_back(vP[1]);\r
- chargeVector[6]->push_back(vP[2]);\r
- chargeVector[7]->push_back(vPt);\r
- chargeVector[8]->push_back(vE);\r
- \r
- \r
- \r
- } //mix track loop\r
-\r
- // calculate balance function for each track in main event\r
- iMainTrackUsed++; // is needed to do no double counting in Balance Function calculation \r
- if(iMainTrackUsed >= (Int_t)chargeVector[0]->size()) break; //do not allow more tracks than in mixed event!\r
- fBalance->CalculateBalance(fCentrality,chargeVector,iMainTrackUsed,bSign);\r
- // clean charge vector afterwards\r
- for(Int_t i = 0; i < 9; i++){ \r
- chargeVector[i]->clear();\r
- }\r
- \r
-\r
- } //main track loop\r
- // }//Vz cut\r
- // }//Vy cut\r
- // }//Vx cut\r
- // }//proper vertexresolution\r
- // }//proper number of contributors\r
- // }//vertex object valid\r
- }//triggered event \r
- }//AOD analysis\r
- }\r
-}\r
-\r
-AliMixInputEventHandler *AliAnalysisTaskEventMixingBF::SetupEventsForMixing() {\r
- //sets the input handlers for event mixing\r
-\r
- AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();\r
- AliMultiInputEventHandler *inEvHMain = dynamic_cast<AliMultiInputEventHandler *>(mgr->GetInputEventHandler());\r
- if (inEvHMain) {\r
-\r
- AliMixInputEventHandler *mixEH = dynamic_cast<AliMixInputEventHandler *>(inEvHMain->GetFirstMultiInputHandler());\r
- if (!mixEH) return nullptr;\r
- if (mixEH->CurrentBinIndex() < 0) {\r
- AliDebug(AliLog::kDebug + 1, "Current event mixEH->CurrentEntry() == -1");\r
- return nullptr;\r
- }\r
- AliDebug(AliLog::kDebug, Form("Mixing %lld %d [%lld,%lld] %d", mixEH->CurrentEntry(), mixEH->CurrentBinIndex(), mixEH->CurrentEntryMain(), mixEH->CurrentEntryMix(), mixEH->NumberMixed()));\r
-\r
- AliInputEventHandler *ihMainCurrent = inEvHMain->GetFirstInputEventHandler();\r
- fMainEvent = ihMainCurrent->GetEvent();\r
-\r
- AliMultiInputEventHandler *inEvHMixedCurrent = mixEH->GetFirstMultiInputHandler(); // for buffer = 1\r
- AliInputEventHandler *ihMixedCurrent = inEvHMixedCurrent->GetFirstInputEventHandler();\r
- fMixEvent = ihMixedCurrent->GetEvent();\r
- \r
- return mixEH;\r
- }\r
- return NULL;\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 "AliMixInputEventHandler.h"
+#include "AliStack.h"
+#include "AliESDtrackCuts.h"
+
+#include "TH2D.h"
+#include "AliPID.h"
+#include "AliPIDResponse.h"
+#include "AliPIDCombined.h"
+
+#include "AliAnalysisTaskEventMixingBF.h"
+#include "AliBalanceEventMixing.h"
+
+
+// Analysis task for the EventMixingBF code
+// Authors: Panos.Christakoglou@nikhef.nl, m.weber@cern.ch
+
+ClassImp(AliAnalysisTaskEventMixingBF)
+
+//________________________________________________________________________
+AliAnalysisTaskEventMixingBF::AliAnalysisTaskEventMixingBF(const char *name)
+: AliAnalysisTaskSE(name),
+ fBalance(0),
+ fRunShuffling(kFALSE),
+ fShuffledBalance(0),
+ fList(0),
+ fListEventMixingBF(0),
+ fListEventMixingBFS(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),
+ fHistPhi(0),
+ fHistPhiBefore(0),
+ fHistPhiAfter(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),
+ fMainEvent(0x0),
+ fMixEvent(0x0)
+{
+ // 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());
+}
+
+//________________________________________________________________________
+AliAnalysisTaskEventMixingBF::~AliAnalysisTaskEventMixingBF() {
+
+ // delete fBalance;
+ // delete fShuffledBalance;
+ // delete fList;
+ // delete fListEventMixingBF;
+ // delete fListEventMixingBFS;
+
+ // 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 AliAnalysisTaskEventMixingBF::UserCreateOutputObjects() {
+ // Create histograms
+ // Called once
+ if(!fBalance) {
+ fBalance = new AliBalanceEventMixing();
+ fBalance->SetAnalysisLevel("ESD");
+ //fBalance->SetNumberOfBins(-1,16);
+ fBalance->SetInterval(-1,-0.8,0.8,16,0.,1.6);
+ }
+ if(fRunShuffling) {
+ if(!fShuffledBalance) {
+ fShuffledBalance = new AliBalanceEventMixing();
+ 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
+ fListEventMixingBF = new TList();
+ fListEventMixingBF->SetName("listEventMixingBF");
+ fListEventMixingBF->SetOwner();
+
+ if(fRunShuffling) {
+ fListEventMixingBFS = new TList();
+ fListEventMixingBFS->SetName("listEventMixingBFShuffled");
+ fListEventMixingBFS->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);
+ 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);
+ 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++){
+ fListEventMixingBF->Add(fBalance->GetHistNp(a));
+ fListEventMixingBF->Add(fBalance->GetHistNn(a));
+ fListEventMixingBF->Add(fBalance->GetHistNpn(a));
+ fListEventMixingBF->Add(fBalance->GetHistNnn(a));
+ fListEventMixingBF->Add(fBalance->GetHistNpp(a));
+ fListEventMixingBF->Add(fBalance->GetHistNnp(a));
+
+ if(fRunShuffling) {
+ fListEventMixingBFS->Add(fShuffledBalance->GetHistNp(a));
+ fListEventMixingBFS->Add(fShuffledBalance->GetHistNn(a));
+ fListEventMixingBFS->Add(fShuffledBalance->GetHistNpn(a));
+ fListEventMixingBFS->Add(fShuffledBalance->GetHistNnn(a));
+ fListEventMixingBFS->Add(fShuffledBalance->GetHistNpp(a));
+ fListEventMixingBFS->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, fListEventMixingBF);
+ if(fRunShuffling) PostData(3, fListEventMixingBFS);
+ if(fUsePID) PostData(4, fHistListPIDQA); //PID
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskEventMixingBF::UserExec(Option_t *) {
+ // Main loop
+ // Called for each event
+ // NOTHING TO DO for event mixing!
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskEventMixingBF::FinishTaskOutput(){
+ //Printf("END EventMixingBF");
+
+ if (!fBalance) {
+ AliError("ERROR: fBalance not available");
+ return;
+ }
+ if(fRunShuffling) {
+ if (!fShuffledBalance) {
+ AliError("ERROR: fShuffledBalance not available");
+ return;
+ }
+ }
+
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskEventMixingBF::Terminate(Option_t *) {
+ // Draw result to the screen
+ // Called once at the end of the query
+
+ // not implemented ...
+
+}
+
+void AliAnalysisTaskEventMixingBF::UserExecMix(Option_t *)
+{
+ // Main loop for event mixing
+
+ TString gAnalysisLevel = fBalance->GetAnalysisLevel();
+
+ AliMixInputEventHandler *mixIEH = SetupEventsForMixing();
+
+ 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> *chargeVector[9]; // original charge
+ for(Int_t i = 0; i < 9; i++){
+ chargeVector[i] = new vector<Double_t>;
+ }
+
+ Double_t vCharge;
+ Double_t vY;
+ Double_t vEta;
+ Double_t vPhi;
+ Double_t vP[3];
+ Double_t vPt;
+ Double_t vE;
+
+ Int_t iMainTrackUsed = -1;
+
+ // -------------------------------------------------------------
+ // At the moment MIXING only for AODs
+ if(mixIEH){
+
+ //AOD analysis (vertex and track cuts also here!!!!)
+ if(gAnalysisLevel == "AOD") {
+ AliAODEvent* aodEventMain = dynamic_cast<AliAODEvent*>(fMainEvent);
+ if(!aodEventMain) {
+ AliError("ERROR: aodEventMain not available");
+ return;
+ }
+ AliAODEvent *aodEventMix = dynamic_cast<AliAODEvent *>(fMixEvent);
+ if(!aodEventMix) {
+ AliError("ERROR: aodEventMix not available");
+ return;
+ }
+
+ // for HBT like cuts need magnetic field sign
+ bSign = (aodEventMain->GetMagneticField() > 0) ? 1 : -1;
+
+ AliAODHeader *aodHeaderMain = dynamic_cast<AliAODHeader*>(aodEventMain->GetHeader());
+ if(!aodHeaderMain) AliFatal("Not a standard AOD");
+
+ // event selection done in AliAnalysisTaskSE::Exec() --> this is not used
+ fHistEventStats->Fill(1); //all events
+
+ // this is not needed (checked in mixing handler!)
+ Bool_t isSelectedMain = kTRUE;
+ Bool_t isSelectedMix = kTRUE;
+
+ if(fUseOfflineTrigger){
+ isSelectedMain = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
+ isSelectedMix = ((AliInputEventHandler*)((AliMultiInputEventHandler *)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetFirstMultiInputHandler())->IsEventSelected();
+ }
+
+ if(isSelectedMain && isSelectedMix) {
+ fHistEventStats->Fill(2); //triggered events
+
+ //Centrality stuff (centrality in AOD header)
+ if(fUseCentrality) {
+ fCentrality = aodHeaderMain->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());
+
+ // QA for centrality estimators
+ fHistCentStats->Fill(0.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("V0M"));
+ fHistCentStats->Fill(1.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("FMD"));
+ fHistCentStats->Fill(2.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("TRK"));
+ fHistCentStats->Fill(3.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("TKL"));
+ fHistCentStats->Fill(4.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("CL0"));
+ fHistCentStats->Fill(5.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("CL1"));
+ fHistCentStats->Fill(6.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("V0MvsFMD"));
+ fHistCentStats->Fill(7.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("TKLvsV0M"));
+ fHistCentStats->Fill(8.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("ZEMvsZDC"));
+
+ // take only events inside centrality class
+ if((fCentrality < fCentralityPercentileMin) || (fCentrality > fCentralityPercentileMax))
+ return;
+
+ // centrality QA (V0M)
+ fHistV0M->Fill(aodEventMain->GetVZEROData()->GetMTotV0A(), aodEventMain->GetVZEROData()->GetMTotV0C());
+
+ // centrality QA (reference tracks)
+ fHistRefTracks->Fill(0.,aodHeaderMain->GetRefMultiplicity());
+ fHistRefTracks->Fill(1.,aodHeaderMain->GetRefMultiplicityPos());
+ fHistRefTracks->Fill(2.,aodHeaderMain->GetRefMultiplicityNeg());
+ fHistRefTracks->Fill(3.,aodHeaderMain->GetTPConlyRefMultiplicity());
+ fHistRefTracks->Fill(4.,aodHeaderMain->GetNumberOfITSClusters(0));
+ fHistRefTracks->Fill(5.,aodHeaderMain->GetNumberOfITSClusters(1));
+ fHistRefTracks->Fill(6.,aodHeaderMain->GetNumberOfITSClusters(2));
+ fHistRefTracks->Fill(7.,aodHeaderMain->GetNumberOfITSClusters(3));
+ fHistRefTracks->Fill(8.,aodHeaderMain->GetNumberOfITSClusters(4));
+ }
+
+ // // this is crashing (Bug in ROOT (to be solved)) but not needed (checked in mixing handler)
+ // const AliAODVertex *vertexMain = aodEventMain->GetPrimaryVertex();
+ // const AliAODVertex *vertexMix = aodEventMix->GetPrimaryVertex();
+
+ // if(vertexMain && vertexMix) {
+ // Double32_t fCovMain[6];
+ // Double32_t fCovMix[6];
+ // vertexMain->GetCovarianceMatrix(fCovMain);
+ // vertexMix->GetCovarianceMatrix(fCovMix);
+
+ // if(vertexMain->GetNContributors() > 0 && vertexMix->GetNContributors() > 0) {
+ // if(fCovMain[5] != 0 && fCovMix[5] != 0) {
+ // fHistEventStats->Fill(3); //events with a proper vertex
+ // if(TMath::Abs(vertexMain->GetX()) < fVxMax && TMath::Abs(vertexMix->GetX()) < fVxMax ) {
+ // if(TMath::Abs(vertexMain->GetY()) < fVyMax && TMath::Abs(vertexMix->GetY()) < fVyMax) {
+ // if(TMath::Abs(vertexMain->GetZ()) < fVzMax && TMath::Abs(vertexMix->GetZ()) < fVzMax) {
+ // fHistEventStats->Fill(4); //analyzed events
+ // fHistVx->Fill(vertexMain->GetX());
+ // fHistVy->Fill(vertexMain->GetY());
+ // fHistVz->Fill(vertexMain->GetZ());
+
+ // Loop over tracks in main event
+ for (Int_t iTracksMain = 0; iTracksMain < aodEventMain->GetNumberOfTracks(); iTracksMain++) {
+ AliAODTrack* aodTrackMain = dynamic_cast<AliAODTrack *>(aodEventMain->GetTrack(iTracksMain));
+ if (!aodTrackMain) {
+ AliError(Form("ERROR: Could not receive track %d", iTracksMain));
+ continue;
+ }
+
+ // AOD track cuts
+
+ // For ESD Filter Information: ANALYSIS/macros/AddTaskESDfilter.C
+ // take only TPC only tracks
+ fHistTrackStats->Fill(aodTrackMain->GetFilterMap());
+ if(!aodTrackMain->TestFilterBit(nAODtrackCutBit)) continue;
+
+ vCharge = aodTrackMain->Charge();
+ vY = aodTrackMain->Y();
+ vEta = aodTrackMain->Eta();
+ vPhi = aodTrackMain->Phi() * TMath::RadToDeg();
+ vE = aodTrackMain->E();
+ vPt = aodTrackMain->Pt();
+ aodTrackMain->PxPyPz(vP);
+
+ Float_t dcaXYMain = aodTrackMain->DCA(); // this is the DCA from global track (not exactly what is cut on)
+ Float_t dcaZMain = aodTrackMain->ZAtDCA(); // this is the DCA from global track (not exactly what is cut on)
+
+
+ // Kinematics cuts from ESD track cuts
+ if( vPt < fPtMin || vPt > fPtMax) continue;
+ if( vEta < fEtaMin || vEta > fEtaMax) continue;
+
+ // Extra DCA cuts (for systematic studies [!= -1])
+ if( fDCAxyCut != -1 && fDCAzCut != -1){
+ if(TMath::Sqrt((dcaXYMain*dcaXYMain)/(fDCAxyCut*fDCAxyCut)+(dcaZMain*dcaZMain)/(fDCAzCut*fDCAzCut)) > 1 ){
+ continue; // 2D cut
+ }
+ }
+
+ // Extra TPC cuts (for systematic studies [!= -1])
+ if( fTPCchi2Cut != -1 && aodTrackMain->Chi2perNDF() > fTPCchi2Cut){
+ continue;
+ }
+ if( fNClustersTPCCut != -1 && aodTrackMain->GetTPCNcls() < fNClustersTPCCut){
+ continue;
+ }
+
+ // fill QA histograms
+ fHistClus->Fill(aodTrackMain->GetITSNcls(),aodTrackMain->GetTPCNcls());
+ fHistDCA->Fill(dcaZMain,dcaXYMain);
+ fHistChi2->Fill(aodTrackMain->Chi2perNDF());
+ fHistPt->Fill(vPt);
+ fHistEta->Fill(vEta);
+ fHistPhi->Fill(vPhi);
+
+ // fill charge vector
+ chargeVector[0]->push_back(vCharge);
+ chargeVector[1]->push_back(vY);
+ chargeVector[2]->push_back(vEta);
+ chargeVector[3]->push_back(vPhi);
+ chargeVector[4]->push_back(vP[0]);
+ chargeVector[5]->push_back(vP[1]);
+ chargeVector[6]->push_back(vP[2]);
+ chargeVector[7]->push_back(vPt);
+ chargeVector[8]->push_back(vE);
+
+ // -------------------------------------------------------------
+ // for each track in main event loop over all tracks in mix event
+ for (Int_t iTracksMix = 0; iTracksMix < aodEventMix->GetNumberOfTracks(); iTracksMix++) {
+
+ AliAODTrack* aodTrackMix = dynamic_cast<AliAODTrack *>(aodEventMix->GetTrack(iTracksMix));
+ if (!aodTrackMix) {
+ AliError(Form("ERROR: Could not receive track %d", iTracksMix));
+ continue;
+ }
+
+ // AOD track cuts
+
+ // For ESD Filter Information: ANALYSIS/macros/AddTaskESDfilter.C
+ // take only TPC only tracks
+ fHistTrackStats->Fill(aodTrackMix->GetFilterMap());
+ if(!aodTrackMix->TestFilterBit(nAODtrackCutBit)) continue;
+
+ vCharge = aodTrackMix->Charge();
+ vY = aodTrackMix->Y();
+ vEta = aodTrackMix->Eta();
+ vPhi = aodTrackMix->Phi() * TMath::RadToDeg();
+ vE = aodTrackMix->E();
+ vPt = aodTrackMix->Pt();
+ aodTrackMix->PxPyPz(vP);
+
+ Float_t dcaXYMix = aodTrackMix->DCA(); // this is the DCA from global track (not exactly what is cut on)
+ Float_t dcaZMix = aodTrackMix->ZAtDCA(); // this is the DCA from global track (not exactly what is cut on)
+
+
+ // Kinematics cuts from ESD track cuts
+ if( vPt < fPtMin || vPt > fPtMax) continue;
+ if( vEta < fEtaMin || vEta > fEtaMax) continue;
+
+ // Extra DCA cuts (for systematic studies [!= -1])
+ if( fDCAxyCut != -1 && fDCAxyCut != -1){
+ if(TMath::Sqrt((dcaXYMix*dcaXYMix)/(fDCAxyCut*fDCAxyCut)+(dcaZMix*dcaZMix)/(fDCAzCut*fDCAzCut)) > 1 ){
+ continue; // 2D cut
+ }
+ }
+
+ // Extra TPC cuts (for systematic studies [!= -1])
+ if( fTPCchi2Cut != -1 && aodTrackMix->Chi2perNDF() > fTPCchi2Cut){
+ continue;
+ }
+ if( fNClustersTPCCut != -1 && aodTrackMix->GetTPCNcls() < fNClustersTPCCut){
+ continue;
+ }
+
+ // fill QA histograms
+ fHistClus->Fill(aodTrackMix->GetITSNcls(),aodTrackMix->GetTPCNcls());
+ fHistDCA->Fill(dcaZMix,dcaXYMix);
+ fHistChi2->Fill(aodTrackMix->Chi2perNDF());
+ fHistPt->Fill(vPt);
+ fHistEta->Fill(vEta);
+ fHistPhi->Fill(vPhi);
+
+ // fill charge vector
+ chargeVector[0]->push_back(vCharge);
+ chargeVector[1]->push_back(vY);
+ chargeVector[2]->push_back(vEta);
+ chargeVector[3]->push_back(vPhi);
+ chargeVector[4]->push_back(vP[0]);
+ chargeVector[5]->push_back(vP[1]);
+ chargeVector[6]->push_back(vP[2]);
+ chargeVector[7]->push_back(vPt);
+ chargeVector[8]->push_back(vE);
+
+
+
+ } //mix track loop
+
+ // calculate balance function for each track in main event
+ iMainTrackUsed++; // is needed to do no double counting in Balance Function calculation
+ if(iMainTrackUsed >= (Int_t)chargeVector[0]->size()) break; //do not allow more tracks than in mixed event!
+ fBalance->CalculateBalance(fCentrality,chargeVector,iMainTrackUsed,bSign);
+ // clean charge vector afterwards
+ for(Int_t i = 0; i < 9; i++){
+ chargeVector[i]->clear();
+ }
+
+
+ } //main track loop
+ // }//Vz cut
+ // }//Vy cut
+ // }//Vx cut
+ // }//proper vertexresolution
+ // }//proper number of contributors
+ // }//vertex object valid
+ }//triggered event
+ }//AOD analysis
+ }
+}
+
+AliMixInputEventHandler *AliAnalysisTaskEventMixingBF::SetupEventsForMixing() {
+ //sets the input handlers for event mixing
+
+ AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+ AliMultiInputEventHandler *inEvHMain = dynamic_cast<AliMultiInputEventHandler *>(mgr->GetInputEventHandler());
+ if (inEvHMain) {
+
+ AliMixInputEventHandler *mixEH = dynamic_cast<AliMixInputEventHandler *>(inEvHMain->GetFirstMultiInputHandler());
+ if (!mixEH) return nullptr;
+ if (mixEH->CurrentBinIndex() < 0) {
+ AliDebug(AliLog::kDebug + 1, "Current event mixEH->CurrentEntry() == -1");
+ return nullptr;
+ }
+ AliDebug(AliLog::kDebug, Form("Mixing %lld %d [%lld,%lld] %d", mixEH->CurrentEntry(), mixEH->CurrentBinIndex(), mixEH->CurrentEntryMain(), mixEH->CurrentEntryMix(), mixEH->NumberMixed()));
+
+ AliInputEventHandler *ihMainCurrent = inEvHMain->GetFirstInputEventHandler();
+ fMainEvent = ihMainCurrent->GetEvent();
+
+ AliMultiInputEventHandler *inEvHMixedCurrent = mixEH->GetFirstMultiInputHandler(); // for buffer = 1
+ AliInputEventHandler *ihMixedCurrent = inEvHMixedCurrent->GetFirstInputEventHandler();
+ fMixEvent = ihMixedCurrent->GetEvent();
+
+ return mixEH;
+ }
+ return NULL;
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff