-#include "TChain.h"\r
-#include "TList.h"\r
-#include "TCanvas.h"\r
-#include "TParticle.h"\r
-#include "TLorentzVector.h"\r
-#include "TGraphErrors.h"\r
-#include "TH1.h"\r
-#include "TH2.h"\r
-#include "TH3.h"\r
-#include "TArrayF.h"\r
-#include "TF1.h"\r
-#include "TRandom.h"\r
-#include "TFile.h"\r
-\r
-#include "AliAnalysisManager.h"\r
-#include "AliLog.h"\r
-\r
-#include "AliEventPoolManager.h" \r
-\r
-#include "AliAnalysisTaskToyModel.h"\r
-#include "AliBalance.h"\r
-#include "AliBalancePsi.h"\r
-#include "AliAnalysisTaskTriggeredBF.h"\r
-\r
-\r
-// Analysis task for the toy model analysis\r
-// Authors: Panos.Christakoglou@nikhef.nl\r
-\r
-using std::cout;\r
-using std::endl;\r
-\r
-ClassImp(AliAnalysisTaskToyModel)\r
-\r
-//________________________________________________________________________\r
-AliAnalysisTaskToyModel::AliAnalysisTaskToyModel() \r
-: TObject(),\r
- fUseDebug(kFALSE),\r
- fBalance(0),\r
- fRunShuffling(kFALSE), fShuffledBalance(0),\r
- fRunMixing(kFALSE), fMixedBalance(0), fPoolMgr(0),\r
- fList(0), fListBF(0), fListBFS(0), fListBFM(0),\r
- fHistEventStats(0),\r
- fHistNumberOfAcceptedParticles(0),\r
- fHistReactionPlane(0),\r
- fHistEtaTotal(0), fHistEta(0),\r
- fHistEtaPhiPos(0), fHistEtaPhiNeg(0),\r
- fHistRapidity(0),\r
- fHistRapidityPions(0), fHistRapidityKaons(0), fHistRapidityProtons(0),\r
- fHistPhi(0),\r
- fHistPhiPions(0), fHistPhiKaons(0), fHistPhiProtons(0),\r
- fHistPt(0),\r
- fHistPtPions(0), fHistPtKaons(0), fHistPtProtons(0),\r
- fTotalMultiplicityMean(100.), fTotalMultiplicitySigma(0.0),\r
- fNetChargeMean(0.0), fNetChargeSigma(0.0),\r
- fPtMin(0.0), fPtMax(100.0),\r
- fEtaMin(-1.0), fEtaMax(1.0),\r
- fUseAcceptanceParameterization(kFALSE), fAcceptanceParameterization(0),\r
- fSimulateDetectorEffects(kFALSE),\r
- fNumberOfInefficientSectors(0),\r
- fInefficiencyFactorInPhi(1.0),\r
- fNumberOfDeadSectors(0),\r
- fEfficiencyDropNearEtaEdges(kFALSE),\r
- fEfficiencyMatrix(0),\r
- fUseAllCharges(kFALSE), fParticleMass(0.0),\r
- fPtSpectraAllCharges(0), fTemperatureAllCharges(100.),\r
- fReactionPlane(0.0),\r
- fAzimuthalAngleAllCharges(0), fDirectedFlowAllCharges(0.0), \r
- fEllipticFlowAllCharges(0.0), fTriangularFlowAllCharges(0.0),\r
- fQuandrangularFlowAllCharges(0.0), fPentangularFlowAllCharges(0.0),\r
- fPionPercentage(0.8), fPionMass(0.0),\r
- fPtSpectraPions(0), fTemperaturePions(100.),\r
- fAzimuthalAnglePions(0), fDirectedFlowPions(0.0), \r
- fEllipticFlowPions(0.0), fTriangularFlowPions(0.0), \r
- fQuandrangularFlowPions(0.0), fPentangularFlowPions(0.0),\r
- fKaonPercentage(0.8), fKaonMass(0.0),\r
- fPtSpectraKaons(0), fTemperatureKaons(100.),\r
- fAzimuthalAngleKaons(0), fDirectedFlowKaons(0.0), \r
- fEllipticFlowKaons(0.0), fTriangularFlowKaons(0.0),\r
- fQuandrangularFlowKaons(0.0), fPentangularFlowKaons(0.0),\r
- fProtonPercentage(0.8), fProtonMass(0.0),\r
- fPtSpectraProtons(0), fTemperatureProtons(100.),\r
- fAzimuthalAngleProtons(0), fDirectedFlowProtons(0.0), \r
- fEllipticFlowProtons(0.0), fTriangularFlowProtons(0.0),\r
- fQuandrangularFlowProtons(0.0), fPentangularFlowProtons(0.0),\r
- fUseDynamicalCorrelations(kFALSE), fDynamicalCorrelationsPercentage(0.1),\r
- fUseJets(kFALSE), fPtAssoc(0) {\r
- // Constructor\r
-}\r
-\r
-//________________________________________________________________________\r
-AliAnalysisTaskToyModel::~AliAnalysisTaskToyModel() {\r
- //Destructor\r
- if(fUseAllCharges) {\r
- delete fPtSpectraAllCharges;\r
- delete fAzimuthalAngleAllCharges;\r
- }\r
- else {\r
- delete fPtSpectraPions;\r
- delete fAzimuthalAnglePions;\r
- delete fPtSpectraKaons;\r
- delete fAzimuthalAngleKaons;\r
- delete fPtSpectraProtons;\r
- delete fAzimuthalAngleProtons;\r
- }\r
- if(fUseJets) delete fPtAssoc;\r
-}\r
-\r
-//________________________________________________________________________\r
-void AliAnalysisTaskToyModel::Init() {\r
- //Initialize objects\r
- //==============gRandom Seed=======================//\r
- gRandom->SetSeed(0); //seed is set to a random value (depending on machine clock)\r
- //==============gRandom Seed=======================//\r
-\r
- //==============Particles and spectra==============//\r
- TParticle *pion = new TParticle();\r
- pion->SetPdgCode(211);\r
- fPionMass = pion->GetMass();\r
-\r
- TParticle *kaon = new TParticle();\r
- kaon->SetPdgCode(321);\r
- fKaonMass = kaon->GetMass();\r
- \r
- TParticle *proton = new TParticle();\r
- proton->SetPdgCode(2212);\r
- fProtonMass = proton->GetMass();\r
-\r
- if(fUseAllCharges) {\r
- fParticleMass = fPionMass;\r
- fPtSpectraAllCharges = new TF1("fPtSpectraAllCharges","x*TMath::Exp(-TMath::Power([0]*[0]+x*x,0.5)/[1])",0.,5.);\r
- fPtSpectraAllCharges->SetParName(0,"Mass");\r
- fPtSpectraAllCharges->SetParName(1,"Temperature");\r
- //fPtSpectraAllCharges = new TF1("fPtSpectraAllCharges","(x^2/TMath::Sqrt(TMath::Power(x,2) + TMath::Power(0.139,2)))*TMath::Power((1. + x/[0]),-[1])",0.,20.);\r
- //fPtSpectraAllCharges->SetParName(0,"pt0");\r
- //fPtSpectraAllCharges->SetParName(1,"b");\r
- }\r
- else {\r
- fPtSpectraPions = new TF1("fPtSpectraPions","x*TMath::Exp(-TMath::Power([0]*[0]+x*x,0.5)/[1])",0.,5.);\r
- fPtSpectraPions->SetParName(0,"Mass");\r
- fPtSpectraPions->SetParName(1,"Temperature");\r
- \r
- fPtSpectraKaons = new TF1("fPtSpectraKaons","x*TMath::Exp(-TMath::Power([0]*[0]+x*x,0.5)/[1])",0.,5.);\r
- fPtSpectraKaons->SetParName(0,"Mass");\r
- fPtSpectraKaons->SetParName(1,"Temperature");\r
- \r
- fPtSpectraProtons = new TF1("fPtSpectraProtons","x*TMath::Exp(-TMath::Power([0]*[0]+x*x,0.5)/[1])",0.,5.);\r
- fPtSpectraProtons->SetParName(0,"Mass");\r
- fPtSpectraProtons->SetParName(1,"Temperature");\r
- }\r
- //==============Particles and spectra==============//\r
-\r
- //==============Flow values==============//\r
- if(fUseAllCharges) {\r
- if(fUseDebug) {\r
- Printf("Directed flow: %lf",fDirectedFlowAllCharges);\r
- Printf("Elliptic flow: %lf",fEllipticFlowAllCharges);\r
- Printf("Triangular flow: %lf",fTriangularFlowAllCharges);\r
- Printf("Quandrangular flow: %lf",fQuandrangularFlowAllCharges);\r
- Printf("Pentangular flow: %lf",fPentangularFlowAllCharges);\r
- }\r
-\r
- fAzimuthalAngleAllCharges = new TF1("fAzimuthalAngleAllCharges","1+2.*[1]*TMath::Cos(x-[0])+2.*[2]*TMath::Cos(2*(x-[0]))+2.*[3]*TMath::Cos(3*(x-[0]))+2.*[4]*TMath::Cos(4*(x-[0]))+2.*[5]*TMath::Cos(5*(x-[0]))",0.,2.*TMath::Pi());\r
- fAzimuthalAngleAllCharges->SetParName(0,"Reaction Plane");\r
- fAzimuthalAngleAllCharges->SetParName(1,"Directed flow");\r
- fAzimuthalAngleAllCharges->SetParName(2,"Elliptic flow"); \r
- fAzimuthalAngleAllCharges->SetParName(3,"Triangular flow");\r
- fAzimuthalAngleAllCharges->SetParName(4,"Quandrangular flow");\r
- fAzimuthalAngleAllCharges->SetParName(5,"Pentangular flow");\r
- }\r
- else {\r
- fAzimuthalAnglePions = new TF1("fAzimuthalAnglePions","1+2.*[1]*TMath::Cos(x-[0])+2.*[2]*TMath::Cos(2*(x-[0]))+2.*[3]*TMath::Cos(3*(x-[0]))+2.*[4]*TMath::Cos(4*(x-[0]))+2.*[5]*TMath::Cos(5*(x-[0]))",0.,2.*TMath::Pi());\r
- fAzimuthalAnglePions->SetParName(0,"Reaction Plane");\r
- fAzimuthalAnglePions->SetParName(1,"Directed flow");\r
- fAzimuthalAnglePions->SetParName(2,"Elliptic flow"); \r
- fAzimuthalAnglePions->SetParName(3,"Triangular flow");\r
- fAzimuthalAnglePions->SetParName(4,"Quandrangular flow");\r
- fAzimuthalAnglePions->SetParName(5,"Pentangular flow");\r
- \r
- fAzimuthalAngleKaons = new TF1("fAzimuthalAngleKaons","1+2.*[1]*TMath::Cos(x-[0])+2.*[2]*TMath::Cos(2*(x-[0]))+2.*[3]*TMath::Cos(3*(x-[0]))+2.*[4]*TMath::Cos(4*(x-[0]))+2.*[5]*TMath::Cos(5*(x-[0]))",0.,2.*TMath::Pi());\r
- fAzimuthalAngleKaons->SetParName(0,"Reaction Plane");\r
- fAzimuthalAngleKaons->SetParName(1,"Directed flow");\r
- fAzimuthalAngleKaons->SetParName(2,"Elliptic flow"); \r
- fAzimuthalAngleKaons->SetParName(3,"Triangular flow");\r
- fAzimuthalAngleKaons->SetParName(4,"Quandrangular flow");\r
- fAzimuthalAngleKaons->SetParName(5,"Pentangular flow");\r
- \r
- fAzimuthalAngleProtons = new TF1("fAzimuthalAngleProtons","1+2.*[1]*TMath::Cos(x-[0])+2.*[2]*TMath::Cos(2*(x-[0]))+2.*[3]*TMath::Cos(3*(x-[0]))+2.*[4]*TMath::Cos(4*(x-[0]))+2.*[5]*TMath::Cos(5*(x-[0]))",0.,2.*TMath::Pi());\r
- fAzimuthalAngleProtons->SetParName(0,"Reaction Plane");\r
- fAzimuthalAngleProtons->SetParName(1,"Directed flow");\r
- fAzimuthalAngleProtons->SetParName(2,"Elliptic flow"); \r
- fAzimuthalAngleProtons->SetParName(3,"Triangular flow");\r
- fAzimuthalAngleProtons->SetParName(4,"Quandrangular flow");\r
- fAzimuthalAngleProtons->SetParName(5,"Pentangular flow");\r
- }\r
- //==============Flow values==============//\r
-\r
- //===================Jets===================//\r
- if(fUseJets) {\r
- fPtAssoc = new TF1("fPtAssoc","x*TMath::Exp(-TMath::Power([0]*[0]+x*x,0.5)/[1])",0.,20.);\r
- fPtAssoc->SetParName(0,"pt0");\r
- fPtAssoc->SetParName(1,"b");\r
- fPtAssoc->SetParameter(0,0.139);\r
- fPtAssoc->SetParameter(1,0.5);\r
- fPtAssoc->SetLineColor(1);\r
- }\r
-\r
- //===================Jets===================//\r
-\r
- //==============Efficiency matrix==============//\r
- if(fSimulateDetectorEffects) SetupEfficiencyMatrix();\r
- //==============Efficiency matrix==============//\r
-}\r
-\r
-//________________________________________________________________________\r
-void AliAnalysisTaskToyModel::SetupEfficiencyMatrix() {\r
- //Setup the efficiency matrix\r
- TH1F *hPt = new TH1F("hPt","",200,0.1,20.1);\r
- TH1F *hEta = new TH1F("hEta","",20,-0.95,0.95);\r
- TH1F *hPhi = new TH1F("hPhi","",72,0.,2.*TMath::Pi());\r
- fEfficiencyMatrix = new TH3F("fEfficiencyMatrix","",\r
- hEta->GetNbinsX(),\r
- hEta->GetXaxis()->GetXmin(),\r
- hEta->GetXaxis()->GetXmax(),\r
- hPt->GetNbinsX(),\r
- hPt->GetXaxis()->GetXmin(),\r
- hPt->GetXaxis()->GetXmax(),\r
- hPhi->GetNbinsX(),\r
- hPhi->GetXaxis()->GetXmin(),\r
- hPhi->GetXaxis()->GetXmax());\r
-\r
- //Efficiency in pt\r
- Double_t epsilon[20] = {0.3,0.6,0.77,0.79,0.80,0.80,0.80,0.80,0.80,0.80,0.80,0.80,0.80,0.80,0.80,0.80,0.80,0.80,0.80,0.80};\r
- for(Int_t i=1;i<=20;i++) {\r
- hPt->SetBinContent(i,epsilon[i-1]);\r
- hPt->SetBinError(i,0.01);\r
- }\r
- for(Int_t i=21;i<=200;i++) {\r
- hPt->SetBinContent(i,epsilon[19]);\r
- hPt->SetBinError(i,0.01);\r
- }\r
-\r
- //Efficiency in eta\r
- for(Int_t i=1;i<=hEta->GetNbinsX();i++) {\r
- hEta->SetBinContent(i,1.0);\r
- hEta->SetBinError(i,0.01);\r
- }\r
- if(fEfficiencyDropNearEtaEdges) {\r
- hEta->SetBinContent(1,0.7); hEta->SetBinContent(2,0.8);\r
- hEta->SetBinContent(3,0.9);\r
- hEta->SetBinContent(18,0.9); hEta->SetBinContent(19,0.8);\r
- hEta->SetBinContent(20,0.7);\r
- }\r
-\r
- //Efficiency in phi\r
- for(Int_t i=1;i<=hPhi->GetNbinsX();i++) {\r
- hPhi->SetBinContent(i,1.0);\r
- hPhi->SetBinError(i,0.01);\r
- }\r
- for(Int_t i=1;i<=fNumberOfInefficientSectors;i++)\r
- hPhi->SetBinContent(hPhi->FindBin(hPhi->GetRandom()),fInefficiencyFactorInPhi);\r
- for(Int_t i=1;i<=fNumberOfDeadSectors;i++)\r
- hPhi->SetBinContent(hPhi->FindBin(hPhi->GetRandom()),0.0);\r
- \r
- //Fill the 3D efficiency map\r
- for(Int_t iBinX = 1; iBinX<=fEfficiencyMatrix->GetXaxis()->GetNbins();iBinX++) {\r
- //cout<<"==================================="<<endl;\r
- for(Int_t iBinY = 1; iBinY<=fEfficiencyMatrix->GetYaxis()->GetNbins();iBinY++) {\r
- //cout<<"==================================="<<endl;\r
- for(Int_t iBinZ = 1; iBinZ<=fEfficiencyMatrix->GetZaxis()->GetNbins();iBinZ++) {\r
- fEfficiencyMatrix->SetBinContent(iBinX,iBinY,iBinZ,hEta->GetBinContent(iBinX)*hPt->GetBinContent(iBinY)*hPhi->GetBinContent(iBinZ));\r
- //cout<<"Eta: "<<hEta->GetBinCenter(iBinX)<<" - Pt: "<<hPt->GetBinCenter(iBinY)<<" - Phi: "<<hPhi->GetBinCenter(iBinZ)<<" - "<<hEta->GetBinContent(iBinX)<<" , "<<hPt->GetBinContent(iBinY)<<" , "<<hPhi->GetBinContent(iBinZ)<<" - Efficiency: "<<hEta->GetBinContent(iBinX)*hPt->GetBinContent(iBinY)*hPhi->GetBinContent(iBinZ)<<endl;\r
- }\r
- }\r
- }\r
-}\r
-\r
-//________________________________________________________________________\r
-void AliAnalysisTaskToyModel::CreateOutputObjects() {\r
- // Create histograms\r
- // Called once\r
-\r
- // global switch disabling the reference \r
- // (to avoid "Replacing existing TH1" if several wagons are created in train)\r
- Bool_t oldStatus = TH1::AddDirectoryStatus();\r
- TH1::AddDirectory(kFALSE);\r
-\r
- if(!fBalance) {\r
- fBalance = new AliBalancePsi();\r
- fBalance->SetDeltaEtaMax(2.0);\r
- }\r
- if(fRunShuffling) {\r
- if(!fShuffledBalance) {\r
- fShuffledBalance = new AliBalancePsi();\r
- fShuffledBalance->SetDeltaEtaMax(2.0);\r
- }\r
- }\r
- if(fRunMixing) {\r
- if(!fMixedBalance) {\r
- fMixedBalance = new AliBalancePsi();\r
- fMixedBalance->SetDeltaEtaMax(2.0);\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
- if(fRunMixing) {\r
- fListBFM = new TList();\r
- fListBFM->SetName("listBFMixed");\r
- fListBFM->SetOwner();\r
- }\r
-\r
- //==============QA================//\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
- fHistNumberOfAcceptedParticles = new TH1F("fHistNumberOfAcceptedParticles",";N_{acc.};Entries",10000,0,10000);\r
- fList->Add(fHistNumberOfAcceptedParticles);\r
- \r
- fHistReactionPlane = new TH1F("fHistReactionPlane","Reaction plane angle;#Psi [rad];Entries",1000,0.,2.*TMath::Pi());\r
- fList->Add(fHistReactionPlane);\r
-\r
- //Pseudo-rapidity\r
- fHistEtaTotal = new TH1F("fHistEtaTotal","Pseudo-rapidity (full phase space);#eta;Entries",1000,-15.,15.); \r
- fList->Add(fHistEtaTotal);\r
-\r
- fHistEta = new TH1F("fHistEta","Pseudo-rapidity (acceptance);#eta;Entries",1000,-1.5,1.5); \r
- fList->Add(fHistEta);\r
-\r
- fHistEtaPhiPos = new TH2F("fHistEtaPhiPos","#eta-#phi distribution (+);#eta;#varphi (rad)",80,-2.,2.,72,-TMath::Pi()/2.,3.*TMath::Pi()/2.);\r
- fList->Add(fHistEtaPhiPos);\r
- fHistEtaPhiNeg = new TH2F("fHistEtaPhiNeg","#eta-#phi distribution (-);#eta;#varphi (rad)",80,-2.,2.,72,-TMath::Pi()/2.,3.*TMath::Pi()/2.);\r
- fList->Add(fHistEtaPhiNeg);\r
-\r
- //Rapidity\r
- fHistRapidity = new TH1F("fHistRapidity","Rapidity (acceptance);y;Entries",1000,-1.5,1.5); \r
- fList->Add(fHistRapidity);\r
- fHistRapidityPions = new TH1F("fHistRapidityPions","Rapidity (acceptance - pions);y;Entries",1000,-1.5,1.5); \r
- fList->Add(fHistRapidityPions);\r
- fHistRapidityKaons = new TH1F("fHistRapidityKaons","Rapidity (acceptance - kaons);y;Entries",1000,-1.5,1.5); \r
- fList->Add(fHistRapidityKaons);\r
- fHistRapidityProtons = new TH1F("fHistRapidityProtons","Rapidity (acceptance - protons);y;Entries",1000,-1.5,1.5); \r
- fList->Add(fHistRapidityProtons);\r
-\r
- //Phi\r
- fHistPhi = new TH1F("fHistPhi","Phi (acceptance);#phi (rad);Entries",1000,0.,2*TMath::Pi());\r
- fList->Add(fHistPhi);\r
-\r
- fHistPhiPions = new TH1F("fHistPhiPions","Phi (acceptance - pions);#phi (rad);Entries",1000,0.,2*TMath::Pi());\r
- fList->Add(fHistPhiPions);\r
- fHistPhiKaons = new TH1F("fHistPhiKaons","Phi (acceptance - kaons);#phi (rad);Entries",1000,0.,2*TMath::Pi());\r
- fList->Add(fHistPhiKaons);\r
- fHistPhiProtons = new TH1F("fHistPhiProtons","Phi (acceptance - protons);#phi (rad);Entries",1000,0.,2*TMath::Pi());\r
- fList->Add(fHistPhiProtons);\r
-\r
- //Pt\r
- fHistPt = new TH1F("fHistPt","Pt (acceptance);p_{t} (GeV/c);Entries",1000,0.,10.);\r
- fList->Add(fHistPt);\r
-\r
- fHistPtPions = new TH1F("fHistPtPions","Pt (acceptance - pions);p_{t} (GeV/c);Entries",1000,0.,10.);\r
- fList->Add(fHistPtPions);\r
- fHistPtKaons = new TH1F("fHistPtKaons","Pt (acceptance - kaons);p_{t} (GeV/c);Entries",1000,0.,10.);\r
- fList->Add(fHistPtKaons);\r
- fHistPtProtons = new TH1F("fHistPtProtons","Pt (acceptance - protons);p_{t} (GeV/c);Entries",1000,0.,10.);\r
- fList->Add(fHistPtProtons);\r
-\r
- if(fEfficiencyMatrix) fList->Add(fEfficiencyMatrix);\r
-\r
- //==============Balance function histograms================//\r
- // Initialize histograms if not done yet\r
- if(!fBalance->GetHistNp()){\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()) {\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
- fListBF->Add(fBalance->GetHistNp());\r
- fListBF->Add(fBalance->GetHistNn());\r
- fListBF->Add(fBalance->GetHistNpn());\r
- fListBF->Add(fBalance->GetHistNnn());\r
- fListBF->Add(fBalance->GetHistNpp());\r
- fListBF->Add(fBalance->GetHistNnp());\r
-\r
- if(fRunShuffling) {\r
- fListBFS->Add(fShuffledBalance->GetHistNp());\r
- fListBFS->Add(fShuffledBalance->GetHistNn());\r
- fListBFS->Add(fShuffledBalance->GetHistNpn());\r
- fListBFS->Add(fShuffledBalance->GetHistNnn());\r
- fListBFS->Add(fShuffledBalance->GetHistNpp());\r
- fListBFS->Add(fShuffledBalance->GetHistNnp());\r
- }\r
-\r
- if(fRunMixing) {\r
- fListBFM->Add(fMixedBalance->GetHistNp());\r
- fListBFM->Add(fMixedBalance->GetHistNn());\r
- fListBFM->Add(fMixedBalance->GetHistNpn());\r
- fListBFM->Add(fMixedBalance->GetHistNnn());\r
- fListBFM->Add(fMixedBalance->GetHistNpp());\r
- fListBFM->Add(fMixedBalance->GetHistNnp());\r
- }\r
-\r
- // Event Mixing\r
- if(fRunMixing){\r
- Int_t fMixingTracks = 2000;\r
- Int_t trackDepth = fMixingTracks; \r
- Int_t poolsize = 1000; // Maximum number of events, ignored in the present implemented of AliEventPoolManager\r
- \r
- // centrality bins\r
- // Double_t centralityBins[] = {0.,1.,2.,3.,4.,5.,7.,10.,20.,30.,40.,50.,60.,70.,80.,100.}; // SHOULD BE DEDUCED FROM CREATED ALITHN!!!\r
- // Double_t* centbins = centralityBins;\r
- // Int_t nCentralityBins = sizeof(centralityBins) / sizeof(Double_t) - 1;\r
-\r
- // multiplicity bins\r
- Double_t multiplicityBins[] = {0,10,20,30,40,50,60,70,80,100,100000}; // SHOULD BE DEDUCED FROM CREATED ALITHN!!!\r
- Double_t* multbins = multiplicityBins;\r
- Int_t nMultiplicityBins = sizeof(multiplicityBins) / sizeof(Double_t) - 1;\r
- \r
- // Zvtx bins\r
- Double_t vertexBins[] = {-10., -7., -5., -3., -1., 1., 3., 5., 7., 10.}; // SHOULD BE DEDUCED FROM CREATED ALITHN!!!\r
- Double_t* vtxbins = vertexBins;\r
- Int_t nVertexBins = sizeof(vertexBins) / sizeof(Double_t) - 1;\r
- \r
- // Event plane angle (Psi) bins\r
- // Double_t psiBins[] = {0.,45.,135.,215.,305.,360.}; // SHOULD BE DEDUCED FROM CREATED ALITHN!!!\r
- // Double_t* psibins = psiBins;\r
- // Int_t nPsiBins = sizeof(psiBins) / sizeof(Double_t) - 1;\r
- \r
- // // run the event mixing also in bins of event plane (statistics!)\r
- // if(fRunMixingEventPlane){\r
- // if(fEventClass=="Multiplicity"){\r
- // fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nMultiplicityBins, multbins, nVertexBins, vtxbins, nPsiBins, psibins);\r
- // }\r
- // else{\r
- // fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centbins, nVertexBins, vtxbins, nPsiBins, psibins);\r
- // }\r
- // }\r
- // else{\r
- //if(fEventClass=="Multiplicity"){\r
- fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nMultiplicityBins, multbins, nVertexBins, vtxbins);\r
- //}\r
- //else{\r
- //fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centbins, nVertexBins, vtxbins);\r
- //}\r
- }\r
-\r
- TH1::AddDirectory(oldStatus);\r
-}\r
-\r
-//________________________________________________________________________\r
-void AliAnalysisTaskToyModel::Run(Int_t nEvents) {\r
- // Main loop\r
- // Called for each event\r
- Short_t vCharge = 0;\r
- Float_t vY = 0.0;\r
- Float_t vEta = 0.0;\r
- Float_t vPhi = 0.0;\r
- Float_t vP[3] = {0.,0.,0.};\r
- Float_t vPt = 0.0;\r
- Float_t vE = 0.0;\r
- Bool_t isPion = kFALSE, isKaon = kFALSE, isProton = kFALSE;\r
-\r
- Double_t gDecideCharge = 0.;\r
-\r
- if(fUseAllCharges) {\r
- //fPtSpectraAllCharges->SetParameter(0,fParticleMass);\r
- fPtSpectraAllCharges->SetParameter(0,1.05);\r
- fPtSpectraAllCharges->SetParameter(1,fTemperatureAllCharges);\r
-\r
- fAzimuthalAngleAllCharges->SetParameter(1,fDirectedFlowAllCharges);\r
- fAzimuthalAngleAllCharges->SetParameter(2,fEllipticFlowAllCharges);\r
- fAzimuthalAngleAllCharges->SetParameter(3,fTriangularFlowAllCharges);\r
- fAzimuthalAngleAllCharges->SetParameter(4,fQuandrangularFlowAllCharges);\r
- fAzimuthalAngleAllCharges->SetParameter(5,fPentangularFlowAllCharges);\r
- }\r
- else {\r
- fPtSpectraPions->SetParameter(0,fPionMass);\r
- fPtSpectraPions->SetParameter(1,fTemperaturePions);\r
- fPtSpectraKaons->SetParameter(0,fKaonMass);\r
- fPtSpectraKaons->SetParameter(1,fTemperatureKaons);\r
- fPtSpectraProtons->SetParameter(0,fProtonMass);\r
- fPtSpectraProtons->SetParameter(1,fTemperatureProtons);\r
-\r
- fAzimuthalAnglePions->SetParameter(1,fDirectedFlowPions);\r
- fAzimuthalAnglePions->SetParameter(2,fEllipticFlowPions);\r
- fAzimuthalAnglePions->SetParameter(3,fTriangularFlowPions);\r
- fAzimuthalAnglePions->SetParameter(4,fQuandrangularFlowPions);\r
- fAzimuthalAnglePions->SetParameter(5,fPentangularFlowPions);\r
-\r
- fAzimuthalAngleKaons->SetParameter(1,fDirectedFlowKaons);\r
- fAzimuthalAngleKaons->SetParameter(2,fEllipticFlowKaons);\r
- fAzimuthalAngleKaons->SetParameter(3,fTriangularFlowKaons);\r
- fAzimuthalAngleKaons->SetParameter(4,fQuandrangularFlowKaons);\r
- fAzimuthalAngleKaons->SetParameter(5,fPentangularFlowKaons);\r
-\r
- fAzimuthalAngleProtons->SetParameter(1,fDirectedFlowProtons);\r
- fAzimuthalAngleProtons->SetParameter(2,fEllipticFlowProtons);\r
- fAzimuthalAngleProtons->SetParameter(3,fTriangularFlowProtons);\r
- fAzimuthalAngleProtons->SetParameter(4,fQuandrangularFlowProtons);\r
- fAzimuthalAngleProtons->SetParameter(5,fPentangularFlowProtons);\r
- }\r
-\r
-\r
- for(Int_t iEvent = 0; iEvent < nEvents; iEvent++) {\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
- // TObjArray for the accepted particles \r
- // (has to be done here, otherwise mxing with event pool does not work, overwriting pointers!)\r
- TObjArray *tracksMain = new TObjArray();\r
- tracksMain->SetOwner(kTRUE);\r
- TObjArray *tracksMixing = 0x0;\r
- if(fRunMixing) {\r
- tracksMixing = new TObjArray();\r
- tracksMixing->SetOwner(kTRUE);\r
- }\r
-\r
- tracksMain->Clear();\r
- if(fRunMixing) tracksMixing->Clear();\r
- \r
- fHistEventStats->Fill(1);\r
- fHistEventStats->Fill(2);\r
- fHistEventStats->Fill(3);\r
-\r
- if((iEvent%1000) == 0) \r
- cout<<"Event: "<<iEvent<<"/"<<nEvents<<endl;\r
- \r
- //Multiplicities\r
- Int_t nMultiplicity = (Int_t)(gRandom->Gaus(fTotalMultiplicityMean,fTotalMultiplicitySigma));\r
- Int_t nNetCharge = (Int_t)(gRandom->Gaus(fNetChargeMean,fNetChargeSigma));\r
- \r
- Int_t nGeneratedPositive = (Int_t)((nMultiplicity/2) + nNetCharge);\r
- Int_t nGeneratedNegative = nMultiplicity - nGeneratedPositive;\r
- if(fUseDebug) \r
- Printf("Total multiplicity: %d - Generated positive: %d - Generated negative: %d",nMultiplicity,nGeneratedPositive,nGeneratedNegative);\r
- \r
- //Randomization of the reaction plane\r
- fReactionPlane = 2.0*TMath::Pi()*gRandom->Rndm();\r
- //fReactionPlane = 0.0;\r
- if(fUseAllCharges) \r
- fAzimuthalAngleAllCharges->SetParameter(0,fReactionPlane);\r
- else {\r
- fAzimuthalAnglePions->SetParameter(0,fReactionPlane);\r
- fAzimuthalAngleKaons->SetParameter(0,fReactionPlane);\r
- fAzimuthalAngleProtons->SetParameter(0,fReactionPlane);\r
- }\r
- \r
- Int_t gNumberOfAcceptedParticles = 0;\r
- Int_t gNumberOfAcceptedPositiveParticles = 0;\r
- Int_t gNumberOfAcceptedNegativeParticles = 0;\r
- Int_t nGeneratedPions = 0, nGeneratedKaons = 0, nGeneratedProtons = 0;\r
- \r
- //Generate particles\r
- for(Int_t iParticleCount = 0; iParticleCount < nMultiplicity; iParticleCount++) {\r
- isPion = kFALSE; isKaon = kFALSE; isProton = kFALSE;\r
- if(fUseDebug) \r
- Printf("Generating positive: %d(%d)",iParticleCount+1,nGeneratedPositive);\r
-\r
- //Pseudo-rapidity sampled from a Gaussian centered @ 0\r
- vEta = gRandom->Gaus(0.0,4.0);\r
-\r
- //Fill QA histograms (full phase space)\r
- fHistEtaTotal->Fill(vEta);\r
-\r
- //Decide the charge\r
- gDecideCharge = gRandom->Rndm();\r
- if(gDecideCharge <= 1.*nGeneratedPositive/nMultiplicity) \r
- vCharge = 1;\r
- else \r
- vCharge = -1;\r
- \r
- //Acceptance\r
- if((vEta < fEtaMin) || (vEta > fEtaMax)) continue;\r
-\r
- if(!fUseAllCharges) {\r
- //Decide the specie\r
- Double_t randomNumberSpecies = gRandom->Rndm();\r
- if((randomNumberSpecies >= 0.0)&&(randomNumberSpecies < fPionPercentage)) {\r
- nGeneratedPions += 1;\r
- vPt = fPtSpectraPions->GetRandom();\r
- vPhi = fAzimuthalAnglePions->GetRandom();\r
- fParticleMass = fPionMass;\r
- isPion = kTRUE;\r
- }\r
- else if((randomNumberSpecies >= fPionPercentage)&&(randomNumberSpecies < fPionPercentage + fKaonPercentage)) {\r
- nGeneratedKaons += 1;\r
- vPt = fPtSpectraKaons->GetRandom();\r
- vPhi = fAzimuthalAngleKaons->GetRandom();\r
- fParticleMass = fKaonMass;\r
- isKaon = kTRUE;\r
- }\r
- else if((randomNumberSpecies >= fPionPercentage + fKaonPercentage)&&(randomNumberSpecies < fPionPercentage + fKaonPercentage + fProtonPercentage)) {\r
- nGeneratedProtons += 1;\r
- vPt = fPtSpectraProtons->GetRandom();\r
- vPhi = fAzimuthalAngleProtons->GetRandom();\r
- fParticleMass = fProtonMass;\r
- isProton = kTRUE;\r
- }\r
- }\r
- else {\r
- vPt = fPtSpectraAllCharges->GetRandom();\r
- vPhi = fAzimuthalAngleAllCharges->GetRandom();\r
- }\r
- \r
- vP[0] = vPt*TMath::Cos(vPhi);\r
- vP[1] = vPt*TMath::Sin(vPhi);\r
- vP[2] = vPt*TMath::SinH(vEta);\r
- vE = TMath::Sqrt(TMath::Power(fParticleMass,2) +\r
- TMath::Power(vP[0],2) +\r
- TMath::Power(vP[1],2) +\r
- TMath::Power(vP[2],2));\r
- \r
- vY = 0.5*TMath::Log((vE + vP[2])/(vE - vP[2]));\r
- \r
- //pt coverage\r
- if((vPt < fPtMin) || (vPt > fPtMax)) continue;\r
- //Printf("pt: %lf - mins: %lf - max: %lf",vPt,fPtMin,fPtMax);\r
-\r
- //acceptance filter\r
- if(fUseAcceptanceParameterization) {\r
- Double_t gRandomNumberForAcceptance = gRandom->Rndm();\r
- if(gRandomNumberForAcceptance > fAcceptanceParameterization->Eval(vPt)) \r
- continue;\r
- }\r
-\r
- //Detector effects\r
- if(fSimulateDetectorEffects) {\r
- Double_t randomNumber = gRandom->Rndm();\r
- if(randomNumber > fEfficiencyMatrix->GetBinContent(fEfficiencyMatrix->FindBin(vEta,vPt,vPhi)))\r
- continue;\r
- }\r
-\r
- //Fill QA histograms (acceptance);\r
- if(vCharge > 0) {\r
- gNumberOfAcceptedPositiveParticles += 1;\r
- if(vPhi > 3.*TMath::Pi()/2.)\r
- fHistEtaPhiPos->Fill(vEta,vPhi-2.*TMath::Pi());\r
- else\r
- fHistEtaPhiPos->Fill(vEta,vPhi);\r
- }\r
- else {\r
- gNumberOfAcceptedNegativeParticles += 1;\r
- if(vPhi > 3.*TMath::Pi()/2.)\r
- fHistEtaPhiNeg->Fill(vEta,vPhi-2.*TMath::Pi());\r
- else\r
- fHistEtaPhiNeg->Fill(vEta,vPhi);\r
- }\r
-\r
- fHistEta->Fill(vEta);\r
- fHistRapidity->Fill(vY);\r
- fHistPhi->Fill(vPhi);\r
- fHistPt->Fill(vPt);\r
- if(isPion) {\r
- fHistRapidityPions->Fill(vY);\r
- fHistPhiPions->Fill(vPhi);\r
- fHistPtPions->Fill(vPt);\r
- }\r
- else if(isKaon) {\r
- fHistRapidityKaons->Fill(vY);\r
- fHistPhiKaons->Fill(vPhi);\r
- fHistPtKaons->Fill(vPt);\r
- }\r
- else if(isProton) {\r
- fHistRapidityProtons->Fill(vY);\r
- fHistPhiProtons->Fill(vPhi);\r
- fHistPtProtons->Fill(vPt);\r
- }\r
- gNumberOfAcceptedParticles += 1;\r
-\r
- // add the track to the TObjArray\r
- tracksMain->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, 1.0)); \r
- if(fRunMixing) \r
- tracksMixing->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, 1.0)); \r
- }//generated positive particle loop\r
-\r
- //Jets\r
- if(fUseJets) {\r
- const Int_t nAssociated = 3;\r
-\r
- Double_t gPtTrig1 = 0., gPtTrig2 = 0., gPtAssoc = 0.;\r
- Double_t gPhiTrig1 = 0., gPhiTrig2 = 0., gPhiAssoc = 0.;\r
- Double_t gEtaTrig1 = 0., gEtaTrig2 = 0., gEtaAssoc = 0.;\r
- Short_t gChargeTrig1 = 0, gChargeTrig2 = 0, gChargeAssoc = 0;\r
- \r
- Double_t gJetCone = 0.2;\r
-\r
- //First leading particle\r
- gPtTrig1 = gRandom->Uniform(3.,5.);\r
- gEtaTrig1 = gRandom->Uniform(-0.8,0.8);\r
- gPhiTrig1 = gRandom->Uniform(0.,TMath::TwoPi());\r
-\r
- //Decide the charge\r
- gDecideCharge = gRandom->Rndm();\r
- if(gDecideCharge <= 0.5)\r
- gChargeTrig1 = 1;\r
- else \r
- gChargeTrig1 = -1;\r
- \r
- //Acceptance\r
- if((gEtaTrig1 < fEtaMin) || (gEtaTrig1 > fEtaMax)) continue;\r
- //pt coverage\r
- if((gPtTrig1 < fPtMin) || (gPtTrig1 > fPtMax)) continue;\r
- //Printf("pt: %lf - mins: %lf - max: %lf",vPt,fPtMin,fPtMax);\r
-\r
- //acceptance filter\r
- if(fUseAcceptanceParameterization) {\r
- Double_t gRandomNumberForAcceptance = gRandom->Rndm();\r
- if(gRandomNumberForAcceptance > fAcceptanceParameterization->Eval(gPtTrig1)) \r
- continue;\r
- }\r
-\r
- //Detector effects\r
- if(fSimulateDetectorEffects) {\r
- Double_t randomNumber = gRandom->Rndm();\r
- if(randomNumber > fEfficiencyMatrix->GetBinContent(fEfficiencyMatrix->FindBin(gEtaTrig1,gPtTrig1,gPhiTrig1)))\r
- continue;\r
- }\r
-\r
- gNumberOfAcceptedParticles += 1;\r
-\r
- // add the track to the TObjArray\r
- tracksMain->Add(new AliBFBasicParticle(gEtaTrig1, gPhiTrig1, gPtTrig1, gChargeTrig1, 1.0)); \r
- if(fRunMixing) \r
- tracksMixing->Add(new AliBFBasicParticle(gEtaTrig1, gPhiTrig1, gPtTrig1, gChargeTrig1, 1.0)); \r
-\r
- Int_t iAssociated = 0; \r
- while(iAssociated < nAssociated) {\r
- gPtAssoc = fPtAssoc->GetRandom();\r
- if(gPtAssoc < gPtTrig1) {\r
- gEtaAssoc = gRandom->Uniform(gEtaTrig1 - gJetCone/2.,gEtaTrig1 + gJetCone/2.);\r
- gPhiAssoc = gRandom->Uniform(gPhiTrig1 - gJetCone/2.,gPhiTrig1 + gJetCone/2.);\r
- if(gPhiAssoc < 0.) gPhiAssoc += TMath::TwoPi();\r
- else if(gPhiAssoc > TMath::TwoPi()) gPhiAssoc -= TMath::TwoPi();\r
- \r
- iAssociated += 1;\r
-\r
- gDecideCharge = gRandom->Rndm();\r
- if(gDecideCharge <= 0.5)\r
- gChargeAssoc = 1;\r
- else \r
- gChargeAssoc = -1;\r
- \r
- //Acceptance\r
- if((gEtaAssoc < fEtaMin) || (gEtaAssoc > fEtaMax)) continue;\r
- //pt coverage\r
- if((gPtAssoc < fPtMin) || (gPtAssoc > fPtMax)) continue;\r
- //Printf("pt: %lf - mins: %lf - max: %lf",vPt,fPtMin,fPtMax);\r
-\r
- //acceptance filter\r
- if(fUseAcceptanceParameterization) {\r
- Double_t gRandomNumberForAcceptance = gRandom->Rndm();\r
- if(gRandomNumberForAcceptance > fAcceptanceParameterization->Eval(gPtAssoc)) \r
- continue;\r
- }\r
-\r
- //Detector effects\r
- if(fSimulateDetectorEffects) {\r
- Double_t randomNumber = gRandom->Rndm();\r
- if(randomNumber > fEfficiencyMatrix->GetBinContent(fEfficiencyMatrix->FindBin(gEtaAssoc,gPtAssoc,gPhiAssoc)))\r
- continue;\r
- }\r
-\r
- gNumberOfAcceptedParticles += 1;\r
-\r
- // add the track to the TObjArray\r
- tracksMain->Add(new AliBFBasicParticle(gEtaAssoc, gPhiAssoc, gPtAssoc, gChargeAssoc, 1.0)); \r
- if(fRunMixing) \r
- tracksMixing->Add(new AliBFBasicParticle(gEtaAssoc, gPhiAssoc, gPtAssoc, gChargeAssoc, 1.0)); \r
- }//pt,assoc < pt,trig\r
- }//associated\r
- \r
- //back2back\r
- gPtTrig2 = gPtTrig1;\r
- gEtaTrig2 = -gEtaTrig1;\r
- gPhiTrig2 = TMath::Pi() + gPhiTrig1;\r
- if(gPhiTrig2 < 0.) gPhiTrig2 += TMath::TwoPi();\r
- else if(gPhiTrig2 > TMath::TwoPi()) gPhiTrig2 -= TMath::TwoPi();\r
-\r
- //Decide the charge\r
- gDecideCharge = gRandom->Rndm();\r
- if(gDecideCharge <= 0.5)\r
- gChargeTrig2 = 1;\r
- else \r
- gChargeTrig2 = -1;\r
- \r
- //Acceptance\r
- if((gEtaTrig2 < fEtaMin) || (gEtaTrig2 > fEtaMax)) continue;\r
- //pt coverage\r
- if((gPtTrig2 < fPtMin) || (gPtTrig2 > fPtMax)) continue;\r
- //Printf("pt: %lf - mins: %lf - max: %lf",vPt,fPtMin,fPtMax);\r
-\r
- //acceptance filter\r
- if(fUseAcceptanceParameterization) {\r
- Double_t gRandomNumberForAcceptance = gRandom->Rndm();\r
- if(gRandomNumberForAcceptance > fAcceptanceParameterization->Eval(gPtTrig2)) \r
- continue;\r
- }\r
-\r
- //Detector effects\r
- if(fSimulateDetectorEffects) {\r
- Double_t randomNumber = gRandom->Rndm();\r
- if(randomNumber > fEfficiencyMatrix->GetBinContent(fEfficiencyMatrix->FindBin(gEtaTrig2,gPtTrig2,gPhiTrig2)))\r
- continue;\r
- }\r
-\r
- gNumberOfAcceptedParticles += 1;\r
-\r
- // add the track to the TObjArray\r
- tracksMain->Add(new AliBFBasicParticle(gEtaTrig2, gPhiTrig2, gPtTrig2, gChargeTrig2, 1.0)); \r
- if(fRunMixing) \r
- tracksMixing->Add(new AliBFBasicParticle(gEtaTrig2, gPhiTrig2, gPtTrig2, gChargeTrig2, 1.0)); \r
-\r
- iAssociated = 0; \r
- while(iAssociated < nAssociated) {\r
- gPtAssoc = fPtAssoc->GetRandom();\r
- if(gPtAssoc < gPtTrig2) {\r
- gEtaAssoc = gRandom->Uniform(gEtaTrig2 - gJetCone/2.,gEtaTrig2 + gJetCone/2.);\r
- gPhiAssoc = gRandom->Uniform(gPhiTrig2 - gJetCone/2.,gPhiTrig2 + gJetCone/2.);\r
- if(gPhiAssoc < 0.) gPhiAssoc += TMath::TwoPi();\r
- else if(gPhiAssoc > TMath::TwoPi()) gPhiAssoc -= TMath::TwoPi();\r
- \r
- iAssociated += 1;\r
-\r
- gDecideCharge = gRandom->Rndm();\r
- if(gDecideCharge <= 0.5)\r
- gChargeAssoc = 1;\r
- else \r
- gChargeAssoc = -1;\r
- \r
- //Acceptance\r
- if((gEtaAssoc < fEtaMin) || (gEtaAssoc > fEtaMax)) continue;\r
- //pt coverage\r
- if((gPtAssoc < fPtMin) || (gPtAssoc > fPtMax)) continue;\r
- //Printf("pt: %lf - mins: %lf - max: %lf",vPt,fPtMin,fPtMax);\r
-\r
- //acceptance filter\r
- if(fUseAcceptanceParameterization) {\r
- Double_t gRandomNumberForAcceptance = gRandom->Rndm();\r
- if(gRandomNumberForAcceptance > fAcceptanceParameterization->Eval(gPtAssoc)) \r
- continue;\r
- }\r
-\r
- //Detector effects\r
- if(fSimulateDetectorEffects) {\r
- Double_t randomNumber = gRandom->Rndm();\r
- if(randomNumber > fEfficiencyMatrix->GetBinContent(fEfficiencyMatrix->FindBin(gEtaAssoc,gPtAssoc,gPhiAssoc)))\r
- continue;\r
- }\r
-\r
- gNumberOfAcceptedParticles += 1;\r
-\r
- // add the track to the TObjArray\r
- tracksMain->Add(new AliBFBasicParticle(gEtaAssoc, gPhiAssoc, gPtAssoc, gChargeAssoc, 1.0)); \r
- if(fRunMixing) \r
- tracksMixing->Add(new AliBFBasicParticle(gEtaAssoc, gPhiAssoc, gPtAssoc, gChargeAssoc, 1.0)); \r
- }//pt,assoc < pt,trig\r
- }//associated\r
- }//Jet usage\r
-\r
-\r
-\r
-\r
-\r
-\r
-\r
-\r
- \r
- //Dynamical correlations\r
- Int_t nGeneratedPositiveDynamicalCorrelations = 0;\r
- Int_t nGeneratedNegativeDynamicalCorrelations = 0;\r
- /*Double_t vChargePrime = 0;\r
- Double_t vYPrime = 0.0;\r
- Double_t vEtaPrime = 0.0;\r
- Double_t vPhiPrime = 0.0;\r
- Double_t vPPrime[3] = {0.,0.,0.};\r
- Double_t vPtPrime = 0.0;\r
- Double_t vEPrime = 0.0;\r
- //Generate "correlated" particles \r
- if(fUseDynamicalCorrelations) {\r
- Int_t gNumberOfDynamicalCorrelations = (Int_t)(0.5*gNumberOfAcceptedParticles*fDynamicalCorrelationsPercentage);\r
- for(Int_t iDynamicalCorrelations = 0; iDynamicalCorrelations < gNumberOfDynamicalCorrelations; iDynamicalCorrelations++) {\r
- isPion = kFALSE; isKaon = kFALSE; isProton = kFALSE;\r
- \r
- //Pseudo-rapidity sampled from a Gaussian centered @ 0\r
- vEta = gRandom->Gaus(0.0,0.1);\r
- vCharge = 1.0;\r
- nGeneratedPositiveDynamicalCorrelations += 1;\r
- \r
- vEtaPrime = -vEta;\r
- vChargePrime = -1.0;\r
- nGeneratedNegativeDynamicalCorrelations += 1;\r
- \r
- //Acceptance\r
- if((vEta < fEtaMin) || (vEta > fEtaMax)) continue;\r
- if((vEtaPrime < fEtaMin) || (vEtaPrime > fEtaMax)) continue;\r
- \r
- if(!fUseAllCharges) {\r
- //Decide the specie\r
- Double_t randomNumberSpecies = gRandom->Rndm();\r
- if((randomNumberSpecies >= 0.0)&&(randomNumberSpecies < fPionPercentage)) {\r
- nGeneratedPions += 1;\r
- vPt = fPtSpectraPions->GetRandom();\r
- vPhi = fAzimuthalAnglePions->GetRandom();\r
- fParticleMass = fPionMass;\r
- isPion = kTRUE;\r
- }\r
- else if((randomNumberSpecies >= fPionPercentage)&&(randomNumberSpecies < fPionPercentage + fKaonPercentage)) {\r
- nGeneratedKaons += 1;\r
- vPt = fPtSpectraKaons->GetRandom();\r
- vPhi = fAzimuthalAngleKaons->GetRandom();\r
- fParticleMass = fKaonMass;\r
- isKaon = kTRUE;\r
- }\r
- else if((randomNumberSpecies >= fPionPercentage + fKaonPercentage)&&(randomNumberSpecies < fPionPercentage + fKaonPercentage + fProtonPercentage)) {\r
- nGeneratedProtons += 1;\r
- vPt = fPtSpectraProtons->GetRandom();\r
- vPtPrime = vPt;\r
- vPhi = fAzimuthalAngleProtons->GetRandom();\r
- fParticleMass = fProtonMass;\r
- isProton = kTRUE;\r
- }\r
- }\r
- else {\r
- vPt = fPtSpectraAllCharges->GetRandom();\r
- vPhi = fAzimuthalAngleAllCharges->GetRandom();\r
- }\r
- vPtPrime = vPt;\r
- vPhiPrime = vPhi;\r
-\r
- vP[0] = vPt*TMath::Cos(vPhi);\r
- vP[1] = vPt*TMath::Sin(vPhi);\r
- vP[2] = vPt*TMath::SinH(vEta);\r
- vE = TMath::Sqrt(TMath::Power(fParticleMass,2) +\r
- TMath::Power(vP[0],2) +\r
- TMath::Power(vP[1],2) +\r
- TMath::Power(vP[2],2));\r
- \r
- vY = 0.5*TMath::Log((vE + vP[2])/(vE - vP[2]));\r
-\r
- vPPrime[0] = vPtPrime*TMath::Cos(vPhiPrime);\r
- vPPrime[1] = vPtPrime*TMath::Sin(vPhiPrime);\r
- vPPrime[2] = vPtPrime*TMath::SinH(vEtaPrime);\r
- vEPrime = TMath::Sqrt(TMath::Power(fParticleMass,2) +\r
- TMath::Power(vPPrime[0],2) +\r
- TMath::Power(vPPrime[1],2) +\r
- TMath::Power(vPPrime[2],2));\r
- \r
- vYPrime = 0.5*TMath::Log((vEPrime + vPPrime[2])/(vEPrime - vPPrime[2]));\r
- \r
- //pt coverage\r
- if((vPt < fPtMin) || (vPt > fPtMax)) continue;\r
- if((vPtPrime < fPtMin) || (vPtPrime > fPtMax)) continue;\r
-\r
- //acceptance filter\r
- if(fUseAcceptanceParameterization) {\r
- Double_t gRandomNumberForAcceptance = gRandom->Rndm();\r
- if(gRandomNumberForAcceptance > fAcceptanceParameterization->Eval(vPt)) \r
- continue;\r
- \r
- Double_t gRandomNumberForAcceptancePrime = gRandom->Rndm();\r
- if(gRandomNumberForAcceptancePrime > fAcceptanceParameterization->Eval(vPtPrime)) \r
- continue;\r
- }\r
- \r
- // fill charge vector (positive)\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
- if(fRunShuffling) {\r
- chargeVectorShuffle[0]->push_back(vCharge);\r
- chargeVectorShuffle[1]->push_back(vY);\r
- chargeVectorShuffle[2]->push_back(vEta);\r
- chargeVectorShuffle[3]->push_back(vPhi);\r
- chargeVectorShuffle[4]->push_back(vP[0]);\r
- chargeVectorShuffle[5]->push_back(vP[1]);\r
- chargeVectorShuffle[6]->push_back(vP[2]);\r
- chargeVectorShuffle[7]->push_back(vPt);\r
- chargeVectorShuffle[8]->push_back(vE);\r
- }\r
-\r
- // fill charge vector (negative)\r
- chargeVector[0]->push_back(vChargePrime);\r
- chargeVector[1]->push_back(vYPrime);\r
- chargeVector[2]->push_back(vEtaPrime);\r
- chargeVector[3]->push_back(vPhiPrime);\r
- chargeVector[4]->push_back(vPPrime[0]);\r
- chargeVector[5]->push_back(vPPrime[1]);\r
- chargeVector[6]->push_back(vPPrime[2]);\r
- chargeVector[7]->push_back(vPtPrime);\r
- chargeVector[8]->push_back(vEPrime);\r
- \r
- if(fRunShuffling) {\r
- chargeVectorShuffle[0]->push_back(vChargePrime);\r
- chargeVectorShuffle[1]->push_back(vYPrime);\r
- chargeVectorShuffle[2]->push_back(vEtaPrime);\r
- chargeVectorShuffle[3]->push_back(vPhiPrime);\r
- chargeVectorShuffle[4]->push_back(vPPrime[0]);\r
- chargeVectorShuffle[5]->push_back(vPPrime[1]);\r
- chargeVectorShuffle[6]->push_back(vPPrime[2]);\r
- chargeVectorShuffle[7]->push_back(vPtPrime);\r
- chargeVectorShuffle[8]->push_back(vEPrime);\r
- }\r
-\r
- gNumberOfAcceptedParticles += 2;\r
- }//loop over the dynamical correlations\r
- }*/// usage of dynamical correlations\r
-\r
- if(fUseDebug) {\r
- Printf("=======================================================");\r
- Printf("Total: %d - Total positive: %d - Total negative: %d",nMultiplicity,nGeneratedPositive,nGeneratedNegative);\r
- Printf("Accepted positive: %d - Accepted negative: %d",gNumberOfAcceptedPositiveParticles,gNumberOfAcceptedNegativeParticles);\r
- Printf("Correlations: %d - Correlations positive: %d - Correlations negative: %d",nGeneratedPositiveDynamicalCorrelations+nGeneratedNegativeDynamicalCorrelations,nGeneratedPositiveDynamicalCorrelations,nGeneratedNegativeDynamicalCorrelations);\r
- Printf("Number of accepted particles: %d",gNumberOfAcceptedParticles);\r
- if(!fUseAllCharges)\r
- Printf("Pions: %lf - Kaons: %lf - Protons: %lf",1.*nGeneratedPions/nMultiplicity,1.*nGeneratedKaons/nMultiplicity,1.*nGeneratedProtons/nMultiplicity);\r
- //Printf("Calculating the balance function for %d particles",chargeVector[0]->size());\r
- }\r
-\r
- fHistEventStats->Fill(4);\r
- fHistNumberOfAcceptedParticles->Fill(gNumberOfAcceptedParticles);\r
- fHistReactionPlane->Fill(fReactionPlane);\r
-\r
- // Event mixing \r
- if (fRunMixing)\r
- {\r
- // 1. First get an event pool corresponding in mult (cent) and\r
- // zvertex to the current event. Once initialized, the pool\r
- // should contain nMix (reduced) events. This routine does not\r
- // pre-scan the chain. The first several events of every chain\r
- // will be skipped until the needed pools are filled to the\r
- // specified depth. If the pool categories are not too rare, this\r
- // should not be a problem. If they are rare, you could lose`\r
- // statistics.\r
- \r
- // 2. Collect the whole pool's content of tracks into one TObjArray\r
- // (bgTracks), which is effectively a single background super-event.\r
- \r
- // 3. The reduced and bgTracks arrays must both be passed into\r
- // FillCorrelations(). Also nMix should be passed in, so a weight\r
- // of 1./nMix can be applied.\r
- \r
- AliEventPool* pool = fPoolMgr->GetEventPool(1., 0.,fReactionPlane);\r
- \r
- if (!pool){\r
- AliFatal(Form("No pool found for centrality = %f, zVtx = %f, psi = %f", 1., 0.,fReactionPlane));\r
- }\r
- else{\r
- \r
- //pool->SetDebug(1);\r
- \r
- if (pool->IsReady() || pool->NTracksInPool() > 50000 / 10 || pool->GetCurrentNEvents() >= 5){ \r
- \r
- \r
- Int_t nMix = pool->GetCurrentNEvents();\r
- //cout << "nMix = " << nMix << " tracks in pool = " << pool->NTracksInPool() << ", tracks in this event = "<<gNumberOfAcceptedParticles <<", event Plane = "<<fReactionPlane<<endl;\r
- \r
- //((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(2);\r
- //((TH2F*) fListOfHistos->FindObject("mixedDist"))->Fill(centrality, pool->NTracksInPool());\r
- //if (pool->IsReady())\r
- //((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(3);\r
- \r
- // Fill mixed-event histos here \r
- for (Int_t jMix=0; jMix<nMix; jMix++) \r
- {\r
- TObjArray* tracksMixed = pool->GetEvent(jMix);\r
- fMixedBalance->CalculateBalance(fReactionPlane,tracksMain,tracksMixed,1,1.,0.);\r
- }\r
- }\r
- \r
- // Update the Event pool\r
- pool->UpdatePool(tracksMain);\r
- //pool->PrintInfo();\r
- \r
- }//pool NULL check \r
- }//run mixing\r
- \r
- //Calculate the balance function\r
- fBalance->CalculateBalance(fReactionPlane,tracksMain,NULL,1,1.,0.);\r
- //if(fRunMixing)\r
- // fMixedBalance->CalculateBalance(fReactionPlane,tracksMixing,NULL,1,1.,0.); \r
-\r
- }//event loop\r
-} \r
-\r
-//________________________________________________________________________\r
-void AliAnalysisTaskToyModel::FinishOutput() {\r
- //Printf("END BF");\r
-\r
- TFile *gOutput = TFile::Open("outputToyModel.root","recreate");\r
- TDirectoryFile *dir = new TDirectoryFile("PWGCFEbyE.outputBalanceFunctionPsiAnalysis","PWGCFEbyE.outputBalanceFunctionPsiAnalysis");\r
- \r
- fList->SetName("listQA");\r
- fList->SetOwner(kTRUE);\r
- dir->Add(fList);\r
- //fList->Write("listQA",TObject::kSingleKey);\r
-\r
- if (!fBalance) {\r
- AliError("ERROR: fBalance not available");\r
- return;\r
- } \r
- fListBF->SetName("listBF");\r
- fListBF->SetOwner(kTRUE);\r
- dir->Add(fListBF);\r
- //fListBF->Write("listBF",TObject::kSingleKey);\r
-\r
- if(fRunShuffling) {\r
- if (!fShuffledBalance) {\r
- AliError("ERROR: fShuffledBalance not available");\r
- return;\r
- }\r
- fListBFS->SetName("listBFShuffled");\r
- fListBFS->SetOwner(kTRUE);\r
- dir->Add(fListBFS);\r
- //fListBFS->Write("listBFShuffled",TObject::kSingleKey);\r
- }\r
-\r
- if(fRunMixing) {\r
- if (!fMixedBalance) {\r
- AliError("ERROR: fMixedBalance not available");\r
- return;\r
- }\r
- fListBFM->SetName("listBFMixed");\r
- fListBFM->SetOwner(kTRUE);\r
- dir->Add(fListBFM);\r
- //fListBFM->Write("listBFMixed",TObject::kSingleKey);\r
- }\r
-\r
- dir->Write(dir->GetName(),TObject::kSingleKey);\r
- gOutput->Close();\r
-}\r
-\r
+#include "TChain.h"
+#include "TList.h"
+#include "TCanvas.h"
+#include "TParticle.h"
+#include "TLorentzVector.h"
+#include "TGraphErrors.h"
+#include "TH1.h"
+#include "TH2.h"
+#include "TH3.h"
+#include "TArrayF.h"
+#include "TF1.h"
+#include "TRandom.h"
+#include "TFile.h"
+
+#include "AliAnalysisManager.h"
+#include "AliLog.h"
+
+#include "AliEventPoolManager.h"
+
+#include "AliAnalysisTaskToyModel.h"
+#include "AliBalance.h"
+#include "AliBalancePsi.h"
+#include "AliAnalysisTaskTriggeredBF.h"
+
+
+// Analysis task for the toy model analysis
+// Authors: Panos.Christakoglou@nikhef.nl
+
+using std::cout;
+using std::endl;
+
+ClassImp(AliAnalysisTaskToyModel)
+
+//________________________________________________________________________
+AliAnalysisTaskToyModel::AliAnalysisTaskToyModel()
+: TObject(),
+ fUseDebug(kFALSE),
+ fBalance(0),
+ fRunShuffling(kFALSE), fShuffledBalance(0),
+ fRunMixing(kFALSE), fMixedBalance(0), fPoolMgr(0),
+ fList(0), fListBF(0), fListBFS(0), fListBFM(0),
+ fHistEventStats(0),
+ fHistNumberOfAcceptedParticles(0),
+ fHistReactionPlane(0),
+ fHistEtaTotal(0), fHistEta(0),
+ fHistEtaPhiPos(0), fHistEtaPhiNeg(0),
+ fHistRapidity(0),
+ fHistRapidityPions(0), fHistRapidityKaons(0), fHistRapidityProtons(0),
+ fHistPhi(0),
+ fHistPhiPions(0), fHistPhiKaons(0), fHistPhiProtons(0),
+ fHistPt(0),
+ fHistPtPions(0), fHistPtKaons(0), fHistPtProtons(0),
+ fTotalMultiplicityMean(100.), fTotalMultiplicitySigma(0.0),
+ fNetChargeMean(0.0), fNetChargeSigma(0.0),
+ fPtMin(0.0), fPtMax(100.0),
+ fEtaMin(-1.0), fEtaMax(1.0),
+ fUseAcceptanceParameterization(kFALSE), fAcceptanceParameterization(0),
+ fSimulateDetectorEffects(kFALSE),
+ fNumberOfInefficientSectors(0),
+ fInefficiencyFactorInPhi(1.0),
+ fNumberOfDeadSectors(0),
+ fEfficiencyDropNearEtaEdges(kFALSE),
+ fEfficiencyMatrix(0),
+ fUseAllCharges(kFALSE), fParticleMass(0.0),
+ fPtSpectraAllCharges(0), fTemperatureAllCharges(100.),
+ fReactionPlane(0.0),
+ fAzimuthalAngleAllCharges(0), fDirectedFlowAllCharges(0.0),
+ fEllipticFlowAllCharges(0.0), fTriangularFlowAllCharges(0.0),
+ fQuandrangularFlowAllCharges(0.0), fPentangularFlowAllCharges(0.0),
+ fPionPercentage(0.8), fPionMass(0.0),
+ fPtSpectraPions(0), fTemperaturePions(100.),
+ fAzimuthalAnglePions(0), fDirectedFlowPions(0.0),
+ fEllipticFlowPions(0.0), fTriangularFlowPions(0.0),
+ fQuandrangularFlowPions(0.0), fPentangularFlowPions(0.0),
+ fKaonPercentage(0.8), fKaonMass(0.0),
+ fPtSpectraKaons(0), fTemperatureKaons(100.),
+ fAzimuthalAngleKaons(0), fDirectedFlowKaons(0.0),
+ fEllipticFlowKaons(0.0), fTriangularFlowKaons(0.0),
+ fQuandrangularFlowKaons(0.0), fPentangularFlowKaons(0.0),
+ fProtonPercentage(0.8), fProtonMass(0.0),
+ fPtSpectraProtons(0), fTemperatureProtons(100.),
+ fAzimuthalAngleProtons(0), fDirectedFlowProtons(0.0),
+ fEllipticFlowProtons(0.0), fTriangularFlowProtons(0.0),
+ fQuandrangularFlowProtons(0.0), fPentangularFlowProtons(0.0),
+ fUseDynamicalCorrelations(kFALSE), fDynamicalCorrelationsPercentage(0.1),
+ fUseJets(kFALSE), fPtAssoc(0) {
+ // Constructor
+}
+
+//________________________________________________________________________
+AliAnalysisTaskToyModel::~AliAnalysisTaskToyModel() {
+ //Destructor
+ if(fUseAllCharges) {
+ delete fPtSpectraAllCharges;
+ delete fAzimuthalAngleAllCharges;
+ }
+ else {
+ delete fPtSpectraPions;
+ delete fAzimuthalAnglePions;
+ delete fPtSpectraKaons;
+ delete fAzimuthalAngleKaons;
+ delete fPtSpectraProtons;
+ delete fAzimuthalAngleProtons;
+ }
+ if(fUseJets) delete fPtAssoc;
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskToyModel::Init() {
+ //Initialize objects
+ //==============gRandom Seed=======================//
+ gRandom->SetSeed(0); //seed is set to a random value (depending on machine clock)
+ //==============gRandom Seed=======================//
+
+ //==============Particles and spectra==============//
+ TParticle *pion = new TParticle();
+ pion->SetPdgCode(211);
+ fPionMass = pion->GetMass();
+
+ TParticle *kaon = new TParticle();
+ kaon->SetPdgCode(321);
+ fKaonMass = kaon->GetMass();
+
+ TParticle *proton = new TParticle();
+ proton->SetPdgCode(2212);
+ fProtonMass = proton->GetMass();
+
+ if(fUseAllCharges) {
+ fParticleMass = fPionMass;
+ fPtSpectraAllCharges = new TF1("fPtSpectraAllCharges","x*TMath::Exp(-TMath::Power([0]*[0]+x*x,0.5)/[1])",0.,5.);
+ fPtSpectraAllCharges->SetParName(0,"Mass");
+ fPtSpectraAllCharges->SetParName(1,"Temperature");
+ //fPtSpectraAllCharges = new TF1("fPtSpectraAllCharges","(x^2/TMath::Sqrt(TMath::Power(x,2) + TMath::Power(0.139,2)))*TMath::Power((1. + x/[0]),-[1])",0.,20.);
+ //fPtSpectraAllCharges->SetParName(0,"pt0");
+ //fPtSpectraAllCharges->SetParName(1,"b");
+ }
+ else {
+ fPtSpectraPions = new TF1("fPtSpectraPions","x*TMath::Exp(-TMath::Power([0]*[0]+x*x,0.5)/[1])",0.,5.);
+ fPtSpectraPions->SetParName(0,"Mass");
+ fPtSpectraPions->SetParName(1,"Temperature");
+
+ fPtSpectraKaons = new TF1("fPtSpectraKaons","x*TMath::Exp(-TMath::Power([0]*[0]+x*x,0.5)/[1])",0.,5.);
+ fPtSpectraKaons->SetParName(0,"Mass");
+ fPtSpectraKaons->SetParName(1,"Temperature");
+
+ fPtSpectraProtons = new TF1("fPtSpectraProtons","x*TMath::Exp(-TMath::Power([0]*[0]+x*x,0.5)/[1])",0.,5.);
+ fPtSpectraProtons->SetParName(0,"Mass");
+ fPtSpectraProtons->SetParName(1,"Temperature");
+ }
+ //==============Particles and spectra==============//
+
+ //==============Flow values==============//
+ if(fUseAllCharges) {
+ if(fUseDebug) {
+ Printf("Directed flow: %lf",fDirectedFlowAllCharges);
+ Printf("Elliptic flow: %lf",fEllipticFlowAllCharges);
+ Printf("Triangular flow: %lf",fTriangularFlowAllCharges);
+ Printf("Quandrangular flow: %lf",fQuandrangularFlowAllCharges);
+ Printf("Pentangular flow: %lf",fPentangularFlowAllCharges);
+ }
+
+ fAzimuthalAngleAllCharges = new TF1("fAzimuthalAngleAllCharges","1+2.*[1]*TMath::Cos(x-[0])+2.*[2]*TMath::Cos(2*(x-[0]))+2.*[3]*TMath::Cos(3*(x-[0]))+2.*[4]*TMath::Cos(4*(x-[0]))+2.*[5]*TMath::Cos(5*(x-[0]))",0.,2.*TMath::Pi());
+ fAzimuthalAngleAllCharges->SetParName(0,"Reaction Plane");
+ fAzimuthalAngleAllCharges->SetParName(1,"Directed flow");
+ fAzimuthalAngleAllCharges->SetParName(2,"Elliptic flow");
+ fAzimuthalAngleAllCharges->SetParName(3,"Triangular flow");
+ fAzimuthalAngleAllCharges->SetParName(4,"Quandrangular flow");
+ fAzimuthalAngleAllCharges->SetParName(5,"Pentangular flow");
+ }
+ else {
+ fAzimuthalAnglePions = new TF1("fAzimuthalAnglePions","1+2.*[1]*TMath::Cos(x-[0])+2.*[2]*TMath::Cos(2*(x-[0]))+2.*[3]*TMath::Cos(3*(x-[0]))+2.*[4]*TMath::Cos(4*(x-[0]))+2.*[5]*TMath::Cos(5*(x-[0]))",0.,2.*TMath::Pi());
+ fAzimuthalAnglePions->SetParName(0,"Reaction Plane");
+ fAzimuthalAnglePions->SetParName(1,"Directed flow");
+ fAzimuthalAnglePions->SetParName(2,"Elliptic flow");
+ fAzimuthalAnglePions->SetParName(3,"Triangular flow");
+ fAzimuthalAnglePions->SetParName(4,"Quandrangular flow");
+ fAzimuthalAnglePions->SetParName(5,"Pentangular flow");
+
+ fAzimuthalAngleKaons = new TF1("fAzimuthalAngleKaons","1+2.*[1]*TMath::Cos(x-[0])+2.*[2]*TMath::Cos(2*(x-[0]))+2.*[3]*TMath::Cos(3*(x-[0]))+2.*[4]*TMath::Cos(4*(x-[0]))+2.*[5]*TMath::Cos(5*(x-[0]))",0.,2.*TMath::Pi());
+ fAzimuthalAngleKaons->SetParName(0,"Reaction Plane");
+ fAzimuthalAngleKaons->SetParName(1,"Directed flow");
+ fAzimuthalAngleKaons->SetParName(2,"Elliptic flow");
+ fAzimuthalAngleKaons->SetParName(3,"Triangular flow");
+ fAzimuthalAngleKaons->SetParName(4,"Quandrangular flow");
+ fAzimuthalAngleKaons->SetParName(5,"Pentangular flow");
+
+ fAzimuthalAngleProtons = new TF1("fAzimuthalAngleProtons","1+2.*[1]*TMath::Cos(x-[0])+2.*[2]*TMath::Cos(2*(x-[0]))+2.*[3]*TMath::Cos(3*(x-[0]))+2.*[4]*TMath::Cos(4*(x-[0]))+2.*[5]*TMath::Cos(5*(x-[0]))",0.,2.*TMath::Pi());
+ fAzimuthalAngleProtons->SetParName(0,"Reaction Plane");
+ fAzimuthalAngleProtons->SetParName(1,"Directed flow");
+ fAzimuthalAngleProtons->SetParName(2,"Elliptic flow");
+ fAzimuthalAngleProtons->SetParName(3,"Triangular flow");
+ fAzimuthalAngleProtons->SetParName(4,"Quandrangular flow");
+ fAzimuthalAngleProtons->SetParName(5,"Pentangular flow");
+ }
+ //==============Flow values==============//
+
+ //===================Jets===================//
+ if(fUseJets) {
+ fPtAssoc = new TF1("fPtAssoc","x*TMath::Exp(-TMath::Power([0]*[0]+x*x,0.5)/[1])",0.,20.);
+ fPtAssoc->SetParName(0,"pt0");
+ fPtAssoc->SetParName(1,"b");
+ fPtAssoc->SetParameter(0,0.139);
+ fPtAssoc->SetParameter(1,0.5);
+ fPtAssoc->SetLineColor(1);
+ }
+
+ //===================Jets===================//
+
+ //==============Efficiency matrix==============//
+ if(fSimulateDetectorEffects) SetupEfficiencyMatrix();
+ //==============Efficiency matrix==============//
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskToyModel::SetupEfficiencyMatrix() {
+ //Setup the efficiency matrix
+ TH1F *hPt = new TH1F("hPt","",200,0.1,20.1);
+ TH1F *hEta = new TH1F("hEta","",20,-0.95,0.95);
+ TH1F *hPhi = new TH1F("hPhi","",72,0.,2.*TMath::Pi());
+ fEfficiencyMatrix = new TH3F("fEfficiencyMatrix","",
+ hEta->GetNbinsX(),
+ hEta->GetXaxis()->GetXmin(),
+ hEta->GetXaxis()->GetXmax(),
+ hPt->GetNbinsX(),
+ hPt->GetXaxis()->GetXmin(),
+ hPt->GetXaxis()->GetXmax(),
+ hPhi->GetNbinsX(),
+ hPhi->GetXaxis()->GetXmin(),
+ hPhi->GetXaxis()->GetXmax());
+
+ //Efficiency in pt
+ Double_t epsilon[20] = {0.3,0.6,0.77,0.79,0.80,0.80,0.80,0.80,0.80,0.80,0.80,0.80,0.80,0.80,0.80,0.80,0.80,0.80,0.80,0.80};
+ for(Int_t i=1;i<=20;i++) {
+ hPt->SetBinContent(i,epsilon[i-1]);
+ hPt->SetBinError(i,0.01);
+ }
+ for(Int_t i=21;i<=200;i++) {
+ hPt->SetBinContent(i,epsilon[19]);
+ hPt->SetBinError(i,0.01);
+ }
+
+ //Efficiency in eta
+ for(Int_t i=1;i<=hEta->GetNbinsX();i++) {
+ hEta->SetBinContent(i,1.0);
+ hEta->SetBinError(i,0.01);
+ }
+ if(fEfficiencyDropNearEtaEdges) {
+ hEta->SetBinContent(1,0.7); hEta->SetBinContent(2,0.8);
+ hEta->SetBinContent(3,0.9);
+ hEta->SetBinContent(18,0.9); hEta->SetBinContent(19,0.8);
+ hEta->SetBinContent(20,0.7);
+ }
+
+ //Efficiency in phi
+ for(Int_t i=1;i<=hPhi->GetNbinsX();i++) {
+ hPhi->SetBinContent(i,1.0);
+ hPhi->SetBinError(i,0.01);
+ }
+ for(Int_t i=1;i<=fNumberOfInefficientSectors;i++)
+ hPhi->SetBinContent(hPhi->FindBin(hPhi->GetRandom()),fInefficiencyFactorInPhi);
+ for(Int_t i=1;i<=fNumberOfDeadSectors;i++)
+ hPhi->SetBinContent(hPhi->FindBin(hPhi->GetRandom()),0.0);
+
+ //Fill the 3D efficiency map
+ for(Int_t iBinX = 1; iBinX<=fEfficiencyMatrix->GetXaxis()->GetNbins();iBinX++) {
+ //cout<<"==================================="<<endl;
+ for(Int_t iBinY = 1; iBinY<=fEfficiencyMatrix->GetYaxis()->GetNbins();iBinY++) {
+ //cout<<"==================================="<<endl;
+ for(Int_t iBinZ = 1; iBinZ<=fEfficiencyMatrix->GetZaxis()->GetNbins();iBinZ++) {
+ fEfficiencyMatrix->SetBinContent(iBinX,iBinY,iBinZ,hEta->GetBinContent(iBinX)*hPt->GetBinContent(iBinY)*hPhi->GetBinContent(iBinZ));
+ //cout<<"Eta: "<<hEta->GetBinCenter(iBinX)<<" - Pt: "<<hPt->GetBinCenter(iBinY)<<" - Phi: "<<hPhi->GetBinCenter(iBinZ)<<" - "<<hEta->GetBinContent(iBinX)<<" , "<<hPt->GetBinContent(iBinY)<<" , "<<hPhi->GetBinContent(iBinZ)<<" - Efficiency: "<<hEta->GetBinContent(iBinX)*hPt->GetBinContent(iBinY)*hPhi->GetBinContent(iBinZ)<<endl;
+ }
+ }
+ }
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskToyModel::CreateOutputObjects() {
+ // Create histograms
+ // Called once
+
+ // global switch disabling the reference
+ // (to avoid "Replacing existing TH1" if several wagons are created in train)
+ Bool_t oldStatus = TH1::AddDirectoryStatus();
+ TH1::AddDirectory(kFALSE);
+
+ if(!fBalance) {
+ fBalance = new AliBalancePsi();
+ fBalance->SetDeltaEtaMax(2.0);
+ }
+ if(fRunShuffling) {
+ if(!fShuffledBalance) {
+ fShuffledBalance = new AliBalancePsi();
+ fShuffledBalance->SetDeltaEtaMax(2.0);
+ }
+ }
+ if(fRunMixing) {
+ if(!fMixedBalance) {
+ fMixedBalance = new AliBalancePsi();
+ fMixedBalance->SetDeltaEtaMax(2.0);
+ }
+ }
+
+ //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();
+ }
+
+ if(fRunMixing) {
+ fListBFM = new TList();
+ fListBFM->SetName("listBFMixed");
+ fListBFM->SetOwner();
+ }
+
+ //==============QA================//
+ //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);
+
+ fHistNumberOfAcceptedParticles = new TH1F("fHistNumberOfAcceptedParticles",";N_{acc.};Entries",10000,0,10000);
+ fList->Add(fHistNumberOfAcceptedParticles);
+
+ fHistReactionPlane = new TH1F("fHistReactionPlane","Reaction plane angle;#Psi [rad];Entries",1000,0.,2.*TMath::Pi());
+ fList->Add(fHistReactionPlane);
+
+ //Pseudo-rapidity
+ fHistEtaTotal = new TH1F("fHistEtaTotal","Pseudo-rapidity (full phase space);#eta;Entries",1000,-15.,15.);
+ fList->Add(fHistEtaTotal);
+
+ fHistEta = new TH1F("fHistEta","Pseudo-rapidity (acceptance);#eta;Entries",1000,-1.5,1.5);
+ fList->Add(fHistEta);
+
+ fHistEtaPhiPos = new TH2F("fHistEtaPhiPos","#eta-#phi distribution (+);#eta;#varphi (rad)",80,-2.,2.,72,-TMath::Pi()/2.,3.*TMath::Pi()/2.);
+ fList->Add(fHistEtaPhiPos);
+ fHistEtaPhiNeg = new TH2F("fHistEtaPhiNeg","#eta-#phi distribution (-);#eta;#varphi (rad)",80,-2.,2.,72,-TMath::Pi()/2.,3.*TMath::Pi()/2.);
+ fList->Add(fHistEtaPhiNeg);
+
+ //Rapidity
+ fHistRapidity = new TH1F("fHistRapidity","Rapidity (acceptance);y;Entries",1000,-1.5,1.5);
+ fList->Add(fHistRapidity);
+ fHistRapidityPions = new TH1F("fHistRapidityPions","Rapidity (acceptance - pions);y;Entries",1000,-1.5,1.5);
+ fList->Add(fHistRapidityPions);
+ fHistRapidityKaons = new TH1F("fHistRapidityKaons","Rapidity (acceptance - kaons);y;Entries",1000,-1.5,1.5);
+ fList->Add(fHistRapidityKaons);
+ fHistRapidityProtons = new TH1F("fHistRapidityProtons","Rapidity (acceptance - protons);y;Entries",1000,-1.5,1.5);
+ fList->Add(fHistRapidityProtons);
+
+ //Phi
+ fHistPhi = new TH1F("fHistPhi","Phi (acceptance);#phi (rad);Entries",1000,0.,2*TMath::Pi());
+ fList->Add(fHistPhi);
+
+ fHistPhiPions = new TH1F("fHistPhiPions","Phi (acceptance - pions);#phi (rad);Entries",1000,0.,2*TMath::Pi());
+ fList->Add(fHistPhiPions);
+ fHistPhiKaons = new TH1F("fHistPhiKaons","Phi (acceptance - kaons);#phi (rad);Entries",1000,0.,2*TMath::Pi());
+ fList->Add(fHistPhiKaons);
+ fHistPhiProtons = new TH1F("fHistPhiProtons","Phi (acceptance - protons);#phi (rad);Entries",1000,0.,2*TMath::Pi());
+ fList->Add(fHistPhiProtons);
+
+ //Pt
+ fHistPt = new TH1F("fHistPt","Pt (acceptance);p_{t} (GeV/c);Entries",1000,0.,10.);
+ fList->Add(fHistPt);
+
+ fHistPtPions = new TH1F("fHistPtPions","Pt (acceptance - pions);p_{t} (GeV/c);Entries",1000,0.,10.);
+ fList->Add(fHistPtPions);
+ fHistPtKaons = new TH1F("fHistPtKaons","Pt (acceptance - kaons);p_{t} (GeV/c);Entries",1000,0.,10.);
+ fList->Add(fHistPtKaons);
+ fHistPtProtons = new TH1F("fHistPtProtons","Pt (acceptance - protons);p_{t} (GeV/c);Entries",1000,0.,10.);
+ fList->Add(fHistPtProtons);
+
+ if(fEfficiencyMatrix) fList->Add(fEfficiencyMatrix);
+
+ //==============Balance function histograms================//
+ // Initialize histograms if not done yet
+ if(!fBalance->GetHistNp()){
+ AliWarning("Histograms not yet initialized! --> Will be done now");
+ AliWarning("--> Add 'gBalance->InitHistograms()' in your configBalanceFunction");
+ fBalance->InitHistograms();
+ }
+
+ if(fRunShuffling) {
+ if(!fShuffledBalance->GetHistNp()) {
+ AliWarning("Histograms (shuffling) not yet initialized! --> Will be done now");
+ AliWarning("--> Add 'gBalance->InitHistograms()' in your configBalanceFunction");
+ fShuffledBalance->InitHistograms();
+ }
+ }
+
+ fListBF->Add(fBalance->GetHistNp());
+ fListBF->Add(fBalance->GetHistNn());
+ fListBF->Add(fBalance->GetHistNpn());
+ fListBF->Add(fBalance->GetHistNnn());
+ fListBF->Add(fBalance->GetHistNpp());
+ fListBF->Add(fBalance->GetHistNnp());
+
+ if(fRunShuffling) {
+ fListBFS->Add(fShuffledBalance->GetHistNp());
+ fListBFS->Add(fShuffledBalance->GetHistNn());
+ fListBFS->Add(fShuffledBalance->GetHistNpn());
+ fListBFS->Add(fShuffledBalance->GetHistNnn());
+ fListBFS->Add(fShuffledBalance->GetHistNpp());
+ fListBFS->Add(fShuffledBalance->GetHistNnp());
+ }
+
+ if(fRunMixing) {
+ fListBFM->Add(fMixedBalance->GetHistNp());
+ fListBFM->Add(fMixedBalance->GetHistNn());
+ fListBFM->Add(fMixedBalance->GetHistNpn());
+ fListBFM->Add(fMixedBalance->GetHistNnn());
+ fListBFM->Add(fMixedBalance->GetHistNpp());
+ fListBFM->Add(fMixedBalance->GetHistNnp());
+ }
+
+ // Event Mixing
+ if(fRunMixing){
+ Int_t fMixingTracks = 2000;
+ Int_t trackDepth = fMixingTracks;
+ Int_t poolsize = 1000; // Maximum number of events, ignored in the present implemented of AliEventPoolManager
+
+ // centrality bins
+ // Double_t centralityBins[] = {0.,1.,2.,3.,4.,5.,7.,10.,20.,30.,40.,50.,60.,70.,80.,100.}; // SHOULD BE DEDUCED FROM CREATED ALITHN!!!
+ // Double_t* centbins = centralityBins;
+ // Int_t nCentralityBins = sizeof(centralityBins) / sizeof(Double_t) - 1;
+
+ // multiplicity bins
+ Double_t multiplicityBins[] = {0,10,20,30,40,50,60,70,80,100,100000}; // SHOULD BE DEDUCED FROM CREATED ALITHN!!!
+ Double_t* multbins = multiplicityBins;
+ Int_t nMultiplicityBins = sizeof(multiplicityBins) / sizeof(Double_t) - 1;
+
+ // Zvtx bins
+ Double_t vertexBins[] = {-10., -7., -5., -3., -1., 1., 3., 5., 7., 10.}; // SHOULD BE DEDUCED FROM CREATED ALITHN!!!
+ Double_t* vtxbins = vertexBins;
+ Int_t nVertexBins = sizeof(vertexBins) / sizeof(Double_t) - 1;
+
+ // Event plane angle (Psi) bins
+ // Double_t psiBins[] = {0.,45.,135.,215.,305.,360.}; // SHOULD BE DEDUCED FROM CREATED ALITHN!!!
+ // Double_t* psibins = psiBins;
+ // Int_t nPsiBins = sizeof(psiBins) / sizeof(Double_t) - 1;
+
+ // // run the event mixing also in bins of event plane (statistics!)
+ // if(fRunMixingEventPlane){
+ // if(fEventClass=="Multiplicity"){
+ // fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nMultiplicityBins, multbins, nVertexBins, vtxbins, nPsiBins, psibins);
+ // }
+ // else{
+ // fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centbins, nVertexBins, vtxbins, nPsiBins, psibins);
+ // }
+ // }
+ // else{
+ //if(fEventClass=="Multiplicity"){
+ fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nMultiplicityBins, multbins, nVertexBins, vtxbins);
+ //}
+ //else{
+ //fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centbins, nVertexBins, vtxbins);
+ //}
+ }
+
+ TH1::AddDirectory(oldStatus);
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskToyModel::Run(Int_t nEvents) {
+ // Main loop
+ // Called for each event
+ Short_t vCharge = 0;
+ Float_t vY = 0.0;
+ Float_t vEta = 0.0;
+ Float_t vPhi = 0.0;
+ Float_t vP[3] = {0.,0.,0.};
+ Float_t vPt = 0.0;
+ Float_t vE = 0.0;
+ Bool_t isPion = kFALSE, isKaon = kFALSE, isProton = kFALSE;
+
+ Double_t gDecideCharge = 0.;
+
+ if(fUseAllCharges) {
+ //fPtSpectraAllCharges->SetParameter(0,fParticleMass);
+ fPtSpectraAllCharges->SetParameter(0,1.05);
+ fPtSpectraAllCharges->SetParameter(1,fTemperatureAllCharges);
+
+ fAzimuthalAngleAllCharges->SetParameter(1,fDirectedFlowAllCharges);
+ fAzimuthalAngleAllCharges->SetParameter(2,fEllipticFlowAllCharges);
+ fAzimuthalAngleAllCharges->SetParameter(3,fTriangularFlowAllCharges);
+ fAzimuthalAngleAllCharges->SetParameter(4,fQuandrangularFlowAllCharges);
+ fAzimuthalAngleAllCharges->SetParameter(5,fPentangularFlowAllCharges);
+ }
+ else {
+ fPtSpectraPions->SetParameter(0,fPionMass);
+ fPtSpectraPions->SetParameter(1,fTemperaturePions);
+ fPtSpectraKaons->SetParameter(0,fKaonMass);
+ fPtSpectraKaons->SetParameter(1,fTemperatureKaons);
+ fPtSpectraProtons->SetParameter(0,fProtonMass);
+ fPtSpectraProtons->SetParameter(1,fTemperatureProtons);
+
+ fAzimuthalAnglePions->SetParameter(1,fDirectedFlowPions);
+ fAzimuthalAnglePions->SetParameter(2,fEllipticFlowPions);
+ fAzimuthalAnglePions->SetParameter(3,fTriangularFlowPions);
+ fAzimuthalAnglePions->SetParameter(4,fQuandrangularFlowPions);
+ fAzimuthalAnglePions->SetParameter(5,fPentangularFlowPions);
+
+ fAzimuthalAngleKaons->SetParameter(1,fDirectedFlowKaons);
+ fAzimuthalAngleKaons->SetParameter(2,fEllipticFlowKaons);
+ fAzimuthalAngleKaons->SetParameter(3,fTriangularFlowKaons);
+ fAzimuthalAngleKaons->SetParameter(4,fQuandrangularFlowKaons);
+ fAzimuthalAngleKaons->SetParameter(5,fPentangularFlowKaons);
+
+ fAzimuthalAngleProtons->SetParameter(1,fDirectedFlowProtons);
+ fAzimuthalAngleProtons->SetParameter(2,fEllipticFlowProtons);
+ fAzimuthalAngleProtons->SetParameter(3,fTriangularFlowProtons);
+ fAzimuthalAngleProtons->SetParameter(4,fQuandrangularFlowProtons);
+ fAzimuthalAngleProtons->SetParameter(5,fPentangularFlowProtons);
+ }
+
+
+ for(Int_t iEvent = 0; iEvent < nEvents; iEvent++) {
+ // 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>;
+ //}
+
+ // TObjArray for the accepted particles
+ // (has to be done here, otherwise mxing with event pool does not work, overwriting pointers!)
+ TObjArray *tracksMain = new TObjArray();
+ tracksMain->SetOwner(kTRUE);
+ TObjArray *tracksMixing = 0x0;
+ if(fRunMixing) {
+ tracksMixing = new TObjArray();
+ tracksMixing->SetOwner(kTRUE);
+ }
+
+ tracksMain->Clear();
+ if(fRunMixing) tracksMixing->Clear();
+
+ fHistEventStats->Fill(1);
+ fHistEventStats->Fill(2);
+ fHistEventStats->Fill(3);
+
+ if((iEvent%1000) == 0)
+ cout<<"Event: "<<iEvent<<"/"<<nEvents<<endl;
+
+ //Multiplicities
+ Int_t nMultiplicity = (Int_t)(gRandom->Gaus(fTotalMultiplicityMean,fTotalMultiplicitySigma));
+ Int_t nNetCharge = (Int_t)(gRandom->Gaus(fNetChargeMean,fNetChargeSigma));
+
+ Int_t nGeneratedPositive = (Int_t)((nMultiplicity/2) + nNetCharge);
+ Int_t nGeneratedNegative = nMultiplicity - nGeneratedPositive;
+ if(fUseDebug)
+ Printf("Total multiplicity: %d - Generated positive: %d - Generated negative: %d",nMultiplicity,nGeneratedPositive,nGeneratedNegative);
+
+ //Randomization of the reaction plane
+ fReactionPlane = 2.0*TMath::Pi()*gRandom->Rndm();
+ //fReactionPlane = 0.0;
+ if(fUseAllCharges)
+ fAzimuthalAngleAllCharges->SetParameter(0,fReactionPlane);
+ else {
+ fAzimuthalAnglePions->SetParameter(0,fReactionPlane);
+ fAzimuthalAngleKaons->SetParameter(0,fReactionPlane);
+ fAzimuthalAngleProtons->SetParameter(0,fReactionPlane);
+ }
+
+ Int_t gNumberOfAcceptedParticles = 0;
+ Int_t gNumberOfAcceptedPositiveParticles = 0;
+ Int_t gNumberOfAcceptedNegativeParticles = 0;
+ Int_t nGeneratedPions = 0, nGeneratedKaons = 0, nGeneratedProtons = 0;
+
+ //Generate particles
+ for(Int_t iParticleCount = 0; iParticleCount < nMultiplicity; iParticleCount++) {
+ isPion = kFALSE; isKaon = kFALSE; isProton = kFALSE;
+ if(fUseDebug)
+ Printf("Generating positive: %d(%d)",iParticleCount+1,nGeneratedPositive);
+
+ //Pseudo-rapidity sampled from a Gaussian centered @ 0
+ vEta = gRandom->Gaus(0.0,4.0);
+
+ //Fill QA histograms (full phase space)
+ fHistEtaTotal->Fill(vEta);
+
+ //Decide the charge
+ gDecideCharge = gRandom->Rndm();
+ if(gDecideCharge <= 1.*nGeneratedPositive/nMultiplicity)
+ vCharge = 1;
+ else
+ vCharge = -1;
+
+ //Acceptance
+ if((vEta < fEtaMin) || (vEta > fEtaMax)) continue;
+
+ if(!fUseAllCharges) {
+ //Decide the specie
+ Double_t randomNumberSpecies = gRandom->Rndm();
+ if((randomNumberSpecies >= 0.0)&&(randomNumberSpecies < fPionPercentage)) {
+ nGeneratedPions += 1;
+ vPt = fPtSpectraPions->GetRandom();
+ vPhi = fAzimuthalAnglePions->GetRandom();
+ fParticleMass = fPionMass;
+ isPion = kTRUE;
+ }
+ else if((randomNumberSpecies >= fPionPercentage)&&(randomNumberSpecies < fPionPercentage + fKaonPercentage)) {
+ nGeneratedKaons += 1;
+ vPt = fPtSpectraKaons->GetRandom();
+ vPhi = fAzimuthalAngleKaons->GetRandom();
+ fParticleMass = fKaonMass;
+ isKaon = kTRUE;
+ }
+ else if((randomNumberSpecies >= fPionPercentage + fKaonPercentage)&&(randomNumberSpecies < fPionPercentage + fKaonPercentage + fProtonPercentage)) {
+ nGeneratedProtons += 1;
+ vPt = fPtSpectraProtons->GetRandom();
+ vPhi = fAzimuthalAngleProtons->GetRandom();
+ fParticleMass = fProtonMass;
+ isProton = kTRUE;
+ }
+ }
+ else {
+ vPt = fPtSpectraAllCharges->GetRandom();
+ vPhi = fAzimuthalAngleAllCharges->GetRandom();
+ }
+
+ vP[0] = vPt*TMath::Cos(vPhi);
+ vP[1] = vPt*TMath::Sin(vPhi);
+ vP[2] = vPt*TMath::SinH(vEta);
+ vE = TMath::Sqrt(TMath::Power(fParticleMass,2) +
+ TMath::Power(vP[0],2) +
+ TMath::Power(vP[1],2) +
+ TMath::Power(vP[2],2));
+
+ vY = 0.5*TMath::Log((vE + vP[2])/(vE - vP[2]));
+
+ //pt coverage
+ if((vPt < fPtMin) || (vPt > fPtMax)) continue;
+ //Printf("pt: %lf - mins: %lf - max: %lf",vPt,fPtMin,fPtMax);
+
+ //acceptance filter
+ if(fUseAcceptanceParameterization) {
+ Double_t gRandomNumberForAcceptance = gRandom->Rndm();
+ if(gRandomNumberForAcceptance > fAcceptanceParameterization->Eval(vPt))
+ continue;
+ }
+
+ //Detector effects
+ if(fSimulateDetectorEffects) {
+ Double_t randomNumber = gRandom->Rndm();
+ if(randomNumber > fEfficiencyMatrix->GetBinContent(fEfficiencyMatrix->FindBin(vEta,vPt,vPhi)))
+ continue;
+ }
+
+ //Fill QA histograms (acceptance);
+ if(vCharge > 0) {
+ gNumberOfAcceptedPositiveParticles += 1;
+ if(vPhi > 3.*TMath::Pi()/2.)
+ fHistEtaPhiPos->Fill(vEta,vPhi-2.*TMath::Pi());
+ else
+ fHistEtaPhiPos->Fill(vEta,vPhi);
+ }
+ else {
+ gNumberOfAcceptedNegativeParticles += 1;
+ if(vPhi > 3.*TMath::Pi()/2.)
+ fHistEtaPhiNeg->Fill(vEta,vPhi-2.*TMath::Pi());
+ else
+ fHistEtaPhiNeg->Fill(vEta,vPhi);
+ }
+
+ fHistEta->Fill(vEta);
+ fHistRapidity->Fill(vY);
+ fHistPhi->Fill(vPhi);
+ fHistPt->Fill(vPt);
+ if(isPion) {
+ fHistRapidityPions->Fill(vY);
+ fHistPhiPions->Fill(vPhi);
+ fHistPtPions->Fill(vPt);
+ }
+ else if(isKaon) {
+ fHistRapidityKaons->Fill(vY);
+ fHistPhiKaons->Fill(vPhi);
+ fHistPtKaons->Fill(vPt);
+ }
+ else if(isProton) {
+ fHistRapidityProtons->Fill(vY);
+ fHistPhiProtons->Fill(vPhi);
+ fHistPtProtons->Fill(vPt);
+ }
+ gNumberOfAcceptedParticles += 1;
+
+ // add the track to the TObjArray
+ tracksMain->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, 1.0));
+ if(fRunMixing)
+ tracksMixing->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, 1.0));
+ }//generated positive particle loop
+
+ //Jets
+ if(fUseJets) {
+ const Int_t nAssociated = 3;
+
+ Double_t gPtTrig1 = 0., gPtTrig2 = 0., gPtAssoc = 0.;
+ Double_t gPhiTrig1 = 0., gPhiTrig2 = 0., gPhiAssoc = 0.;
+ Double_t gEtaTrig1 = 0., gEtaTrig2 = 0., gEtaAssoc = 0.;
+ Short_t gChargeTrig1 = 0, gChargeTrig2 = 0, gChargeAssoc = 0;
+
+ Double_t gJetCone = 0.2;
+
+ //First leading particle
+ gPtTrig1 = gRandom->Uniform(3.,5.);
+ gEtaTrig1 = gRandom->Uniform(-0.8,0.8);
+ gPhiTrig1 = gRandom->Uniform(0.,TMath::TwoPi());
+
+ //Decide the charge
+ gDecideCharge = gRandom->Rndm();
+ if(gDecideCharge <= 0.5)
+ gChargeTrig1 = 1;
+ else
+ gChargeTrig1 = -1;
+
+ //Acceptance
+ if((gEtaTrig1 < fEtaMin) || (gEtaTrig1 > fEtaMax)) continue;
+ //pt coverage
+ if((gPtTrig1 < fPtMin) || (gPtTrig1 > fPtMax)) continue;
+ //Printf("pt: %lf - mins: %lf - max: %lf",vPt,fPtMin,fPtMax);
+
+ //acceptance filter
+ if(fUseAcceptanceParameterization) {
+ Double_t gRandomNumberForAcceptance = gRandom->Rndm();
+ if(gRandomNumberForAcceptance > fAcceptanceParameterization->Eval(gPtTrig1))
+ continue;
+ }
+
+ //Detector effects
+ if(fSimulateDetectorEffects) {
+ Double_t randomNumber = gRandom->Rndm();
+ if(randomNumber > fEfficiencyMatrix->GetBinContent(fEfficiencyMatrix->FindBin(gEtaTrig1,gPtTrig1,gPhiTrig1)))
+ continue;
+ }
+
+ gNumberOfAcceptedParticles += 1;
+
+ // add the track to the TObjArray
+ tracksMain->Add(new AliBFBasicParticle(gEtaTrig1, gPhiTrig1, gPtTrig1, gChargeTrig1, 1.0));
+ if(fRunMixing)
+ tracksMixing->Add(new AliBFBasicParticle(gEtaTrig1, gPhiTrig1, gPtTrig1, gChargeTrig1, 1.0));
+
+ Int_t iAssociated = 0;
+ while(iAssociated < nAssociated) {
+ gPtAssoc = fPtAssoc->GetRandom();
+ if(gPtAssoc < gPtTrig1) {
+ gEtaAssoc = gRandom->Uniform(gEtaTrig1 - gJetCone/2.,gEtaTrig1 + gJetCone/2.);
+ gPhiAssoc = gRandom->Uniform(gPhiTrig1 - gJetCone/2.,gPhiTrig1 + gJetCone/2.);
+ if(gPhiAssoc < 0.) gPhiAssoc += TMath::TwoPi();
+ else if(gPhiAssoc > TMath::TwoPi()) gPhiAssoc -= TMath::TwoPi();
+
+ iAssociated += 1;
+
+ gDecideCharge = gRandom->Rndm();
+ if(gDecideCharge <= 0.5)
+ gChargeAssoc = 1;
+ else
+ gChargeAssoc = -1;
+
+ //Acceptance
+ if((gEtaAssoc < fEtaMin) || (gEtaAssoc > fEtaMax)) continue;
+ //pt coverage
+ if((gPtAssoc < fPtMin) || (gPtAssoc > fPtMax)) continue;
+ //Printf("pt: %lf - mins: %lf - max: %lf",vPt,fPtMin,fPtMax);
+
+ //acceptance filter
+ if(fUseAcceptanceParameterization) {
+ Double_t gRandomNumberForAcceptance = gRandom->Rndm();
+ if(gRandomNumberForAcceptance > fAcceptanceParameterization->Eval(gPtAssoc))
+ continue;
+ }
+
+ //Detector effects
+ if(fSimulateDetectorEffects) {
+ Double_t randomNumber = gRandom->Rndm();
+ if(randomNumber > fEfficiencyMatrix->GetBinContent(fEfficiencyMatrix->FindBin(gEtaAssoc,gPtAssoc,gPhiAssoc)))
+ continue;
+ }
+
+ gNumberOfAcceptedParticles += 1;
+
+ // add the track to the TObjArray
+ tracksMain->Add(new AliBFBasicParticle(gEtaAssoc, gPhiAssoc, gPtAssoc, gChargeAssoc, 1.0));
+ if(fRunMixing)
+ tracksMixing->Add(new AliBFBasicParticle(gEtaAssoc, gPhiAssoc, gPtAssoc, gChargeAssoc, 1.0));
+ }//pt,assoc < pt,trig
+ }//associated
+
+ //back2back
+ gPtTrig2 = gPtTrig1;
+ gEtaTrig2 = -gEtaTrig1;
+ gPhiTrig2 = TMath::Pi() + gPhiTrig1;
+ if(gPhiTrig2 < 0.) gPhiTrig2 += TMath::TwoPi();
+ else if(gPhiTrig2 > TMath::TwoPi()) gPhiTrig2 -= TMath::TwoPi();
+
+ //Decide the charge
+ gDecideCharge = gRandom->Rndm();
+ if(gDecideCharge <= 0.5)
+ gChargeTrig2 = 1;
+ else
+ gChargeTrig2 = -1;
+
+ //Acceptance
+ if((gEtaTrig2 < fEtaMin) || (gEtaTrig2 > fEtaMax)) continue;
+ //pt coverage
+ if((gPtTrig2 < fPtMin) || (gPtTrig2 > fPtMax)) continue;
+ //Printf("pt: %lf - mins: %lf - max: %lf",vPt,fPtMin,fPtMax);
+
+ //acceptance filter
+ if(fUseAcceptanceParameterization) {
+ Double_t gRandomNumberForAcceptance = gRandom->Rndm();
+ if(gRandomNumberForAcceptance > fAcceptanceParameterization->Eval(gPtTrig2))
+ continue;
+ }
+
+ //Detector effects
+ if(fSimulateDetectorEffects) {
+ Double_t randomNumber = gRandom->Rndm();
+ if(randomNumber > fEfficiencyMatrix->GetBinContent(fEfficiencyMatrix->FindBin(gEtaTrig2,gPtTrig2,gPhiTrig2)))
+ continue;
+ }
+
+ gNumberOfAcceptedParticles += 1;
+
+ // add the track to the TObjArray
+ tracksMain->Add(new AliBFBasicParticle(gEtaTrig2, gPhiTrig2, gPtTrig2, gChargeTrig2, 1.0));
+ if(fRunMixing)
+ tracksMixing->Add(new AliBFBasicParticle(gEtaTrig2, gPhiTrig2, gPtTrig2, gChargeTrig2, 1.0));
+
+ iAssociated = 0;
+ while(iAssociated < nAssociated) {
+ gPtAssoc = fPtAssoc->GetRandom();
+ if(gPtAssoc < gPtTrig2) {
+ gEtaAssoc = gRandom->Uniform(gEtaTrig2 - gJetCone/2.,gEtaTrig2 + gJetCone/2.);
+ gPhiAssoc = gRandom->Uniform(gPhiTrig2 - gJetCone/2.,gPhiTrig2 + gJetCone/2.);
+ if(gPhiAssoc < 0.) gPhiAssoc += TMath::TwoPi();
+ else if(gPhiAssoc > TMath::TwoPi()) gPhiAssoc -= TMath::TwoPi();
+
+ iAssociated += 1;
+
+ gDecideCharge = gRandom->Rndm();
+ if(gDecideCharge <= 0.5)
+ gChargeAssoc = 1;
+ else
+ gChargeAssoc = -1;
+
+ //Acceptance
+ if((gEtaAssoc < fEtaMin) || (gEtaAssoc > fEtaMax)) continue;
+ //pt coverage
+ if((gPtAssoc < fPtMin) || (gPtAssoc > fPtMax)) continue;
+ //Printf("pt: %lf - mins: %lf - max: %lf",vPt,fPtMin,fPtMax);
+
+ //acceptance filter
+ if(fUseAcceptanceParameterization) {
+ Double_t gRandomNumberForAcceptance = gRandom->Rndm();
+ if(gRandomNumberForAcceptance > fAcceptanceParameterization->Eval(gPtAssoc))
+ continue;
+ }
+
+ //Detector effects
+ if(fSimulateDetectorEffects) {
+ Double_t randomNumber = gRandom->Rndm();
+ if(randomNumber > fEfficiencyMatrix->GetBinContent(fEfficiencyMatrix->FindBin(gEtaAssoc,gPtAssoc,gPhiAssoc)))
+ continue;
+ }
+
+ gNumberOfAcceptedParticles += 1;
+
+ // add the track to the TObjArray
+ tracksMain->Add(new AliBFBasicParticle(gEtaAssoc, gPhiAssoc, gPtAssoc, gChargeAssoc, 1.0));
+ if(fRunMixing)
+ tracksMixing->Add(new AliBFBasicParticle(gEtaAssoc, gPhiAssoc, gPtAssoc, gChargeAssoc, 1.0));
+ }//pt,assoc < pt,trig
+ }//associated
+ }//Jet usage
+
+
+
+
+
+
+
+
+
+ //Dynamical correlations
+ Int_t nGeneratedPositiveDynamicalCorrelations = 0;
+ Int_t nGeneratedNegativeDynamicalCorrelations = 0;
+ /*Double_t vChargePrime = 0;
+ Double_t vYPrime = 0.0;
+ Double_t vEtaPrime = 0.0;
+ Double_t vPhiPrime = 0.0;
+ Double_t vPPrime[3] = {0.,0.,0.};
+ Double_t vPtPrime = 0.0;
+ Double_t vEPrime = 0.0;
+ //Generate "correlated" particles
+ if(fUseDynamicalCorrelations) {
+ Int_t gNumberOfDynamicalCorrelations = (Int_t)(0.5*gNumberOfAcceptedParticles*fDynamicalCorrelationsPercentage);
+ for(Int_t iDynamicalCorrelations = 0; iDynamicalCorrelations < gNumberOfDynamicalCorrelations; iDynamicalCorrelations++) {
+ isPion = kFALSE; isKaon = kFALSE; isProton = kFALSE;
+
+ //Pseudo-rapidity sampled from a Gaussian centered @ 0
+ vEta = gRandom->Gaus(0.0,0.1);
+ vCharge = 1.0;
+ nGeneratedPositiveDynamicalCorrelations += 1;
+
+ vEtaPrime = -vEta;
+ vChargePrime = -1.0;
+ nGeneratedNegativeDynamicalCorrelations += 1;
+
+ //Acceptance
+ if((vEta < fEtaMin) || (vEta > fEtaMax)) continue;
+ if((vEtaPrime < fEtaMin) || (vEtaPrime > fEtaMax)) continue;
+
+ if(!fUseAllCharges) {
+ //Decide the specie
+ Double_t randomNumberSpecies = gRandom->Rndm();
+ if((randomNumberSpecies >= 0.0)&&(randomNumberSpecies < fPionPercentage)) {
+ nGeneratedPions += 1;
+ vPt = fPtSpectraPions->GetRandom();
+ vPhi = fAzimuthalAnglePions->GetRandom();
+ fParticleMass = fPionMass;
+ isPion = kTRUE;
+ }
+ else if((randomNumberSpecies >= fPionPercentage)&&(randomNumberSpecies < fPionPercentage + fKaonPercentage)) {
+ nGeneratedKaons += 1;
+ vPt = fPtSpectraKaons->GetRandom();
+ vPhi = fAzimuthalAngleKaons->GetRandom();
+ fParticleMass = fKaonMass;
+ isKaon = kTRUE;
+ }
+ else if((randomNumberSpecies >= fPionPercentage + fKaonPercentage)&&(randomNumberSpecies < fPionPercentage + fKaonPercentage + fProtonPercentage)) {
+ nGeneratedProtons += 1;
+ vPt = fPtSpectraProtons->GetRandom();
+ vPtPrime = vPt;
+ vPhi = fAzimuthalAngleProtons->GetRandom();
+ fParticleMass = fProtonMass;
+ isProton = kTRUE;
+ }
+ }
+ else {
+ vPt = fPtSpectraAllCharges->GetRandom();
+ vPhi = fAzimuthalAngleAllCharges->GetRandom();
+ }
+ vPtPrime = vPt;
+ vPhiPrime = vPhi;
+
+ vP[0] = vPt*TMath::Cos(vPhi);
+ vP[1] = vPt*TMath::Sin(vPhi);
+ vP[2] = vPt*TMath::SinH(vEta);
+ vE = TMath::Sqrt(TMath::Power(fParticleMass,2) +
+ TMath::Power(vP[0],2) +
+ TMath::Power(vP[1],2) +
+ TMath::Power(vP[2],2));
+
+ vY = 0.5*TMath::Log((vE + vP[2])/(vE - vP[2]));
+
+ vPPrime[0] = vPtPrime*TMath::Cos(vPhiPrime);
+ vPPrime[1] = vPtPrime*TMath::Sin(vPhiPrime);
+ vPPrime[2] = vPtPrime*TMath::SinH(vEtaPrime);
+ vEPrime = TMath::Sqrt(TMath::Power(fParticleMass,2) +
+ TMath::Power(vPPrime[0],2) +
+ TMath::Power(vPPrime[1],2) +
+ TMath::Power(vPPrime[2],2));
+
+ vYPrime = 0.5*TMath::Log((vEPrime + vPPrime[2])/(vEPrime - vPPrime[2]));
+
+ //pt coverage
+ if((vPt < fPtMin) || (vPt > fPtMax)) continue;
+ if((vPtPrime < fPtMin) || (vPtPrime > fPtMax)) continue;
+
+ //acceptance filter
+ if(fUseAcceptanceParameterization) {
+ Double_t gRandomNumberForAcceptance = gRandom->Rndm();
+ if(gRandomNumberForAcceptance > fAcceptanceParameterization->Eval(vPt))
+ continue;
+
+ Double_t gRandomNumberForAcceptancePrime = gRandom->Rndm();
+ if(gRandomNumberForAcceptancePrime > fAcceptanceParameterization->Eval(vPtPrime))
+ continue;
+ }
+
+ // fill charge vector (positive)
+ 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);
+
+ if(fRunShuffling) {
+ chargeVectorShuffle[0]->push_back(vCharge);
+ chargeVectorShuffle[1]->push_back(vY);
+ chargeVectorShuffle[2]->push_back(vEta);
+ chargeVectorShuffle[3]->push_back(vPhi);
+ chargeVectorShuffle[4]->push_back(vP[0]);
+ chargeVectorShuffle[5]->push_back(vP[1]);
+ chargeVectorShuffle[6]->push_back(vP[2]);
+ chargeVectorShuffle[7]->push_back(vPt);
+ chargeVectorShuffle[8]->push_back(vE);
+ }
+
+ // fill charge vector (negative)
+ chargeVector[0]->push_back(vChargePrime);
+ chargeVector[1]->push_back(vYPrime);
+ chargeVector[2]->push_back(vEtaPrime);
+ chargeVector[3]->push_back(vPhiPrime);
+ chargeVector[4]->push_back(vPPrime[0]);
+ chargeVector[5]->push_back(vPPrime[1]);
+ chargeVector[6]->push_back(vPPrime[2]);
+ chargeVector[7]->push_back(vPtPrime);
+ chargeVector[8]->push_back(vEPrime);
+
+ if(fRunShuffling) {
+ chargeVectorShuffle[0]->push_back(vChargePrime);
+ chargeVectorShuffle[1]->push_back(vYPrime);
+ chargeVectorShuffle[2]->push_back(vEtaPrime);
+ chargeVectorShuffle[3]->push_back(vPhiPrime);
+ chargeVectorShuffle[4]->push_back(vPPrime[0]);
+ chargeVectorShuffle[5]->push_back(vPPrime[1]);
+ chargeVectorShuffle[6]->push_back(vPPrime[2]);
+ chargeVectorShuffle[7]->push_back(vPtPrime);
+ chargeVectorShuffle[8]->push_back(vEPrime);
+ }
+
+ gNumberOfAcceptedParticles += 2;
+ }//loop over the dynamical correlations
+ }*/// usage of dynamical correlations
+
+ if(fUseDebug) {
+ Printf("=======================================================");
+ Printf("Total: %d - Total positive: %d - Total negative: %d",nMultiplicity,nGeneratedPositive,nGeneratedNegative);
+ Printf("Accepted positive: %d - Accepted negative: %d",gNumberOfAcceptedPositiveParticles,gNumberOfAcceptedNegativeParticles);
+ Printf("Correlations: %d - Correlations positive: %d - Correlations negative: %d",nGeneratedPositiveDynamicalCorrelations+nGeneratedNegativeDynamicalCorrelations,nGeneratedPositiveDynamicalCorrelations,nGeneratedNegativeDynamicalCorrelations);
+ Printf("Number of accepted particles: %d",gNumberOfAcceptedParticles);
+ if(!fUseAllCharges)
+ Printf("Pions: %lf - Kaons: %lf - Protons: %lf",1.*nGeneratedPions/nMultiplicity,1.*nGeneratedKaons/nMultiplicity,1.*nGeneratedProtons/nMultiplicity);
+ //Printf("Calculating the balance function for %d particles",chargeVector[0]->size());
+ }
+
+ fHistEventStats->Fill(4);
+ fHistNumberOfAcceptedParticles->Fill(gNumberOfAcceptedParticles);
+ fHistReactionPlane->Fill(fReactionPlane);
+
+ // Event mixing
+ if (fRunMixing)
+ {
+ // 1. First get an event pool corresponding in mult (cent) and
+ // zvertex to the current event. Once initialized, the pool
+ // should contain nMix (reduced) events. This routine does not
+ // pre-scan the chain. The first several events of every chain
+ // will be skipped until the needed pools are filled to the
+ // specified depth. If the pool categories are not too rare, this
+ // should not be a problem. If they are rare, you could lose`
+ // statistics.
+
+ // 2. Collect the whole pool's content of tracks into one TObjArray
+ // (bgTracks), which is effectively a single background super-event.
+
+ // 3. The reduced and bgTracks arrays must both be passed into
+ // FillCorrelations(). Also nMix should be passed in, so a weight
+ // of 1./nMix can be applied.
+
+ AliEventPool* pool = fPoolMgr->GetEventPool(1., 0.,fReactionPlane);
+
+ if (!pool){
+ AliFatal(Form("No pool found for centrality = %f, zVtx = %f, psi = %f", 1., 0.,fReactionPlane));
+ }
+ else{
+
+ //pool->SetDebug(1);
+
+ if (pool->IsReady() || pool->NTracksInPool() > 50000 / 10 || pool->GetCurrentNEvents() >= 5){
+
+
+ Int_t nMix = pool->GetCurrentNEvents();
+ //cout << "nMix = " << nMix << " tracks in pool = " << pool->NTracksInPool() << ", tracks in this event = "<<gNumberOfAcceptedParticles <<", event Plane = "<<fReactionPlane<<endl;
+
+ //((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(2);
+ //((TH2F*) fListOfHistos->FindObject("mixedDist"))->Fill(centrality, pool->NTracksInPool());
+ //if (pool->IsReady())
+ //((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(3);
+
+ // Fill mixed-event histos here
+ for (Int_t jMix=0; jMix<nMix; jMix++)
+ {
+ TObjArray* tracksMixed = pool->GetEvent(jMix);
+ fMixedBalance->CalculateBalance(fReactionPlane,tracksMain,tracksMixed,1,1.,0.);
+ }
+ }
+
+ // Update the Event pool
+ pool->UpdatePool(tracksMain);
+ //pool->PrintInfo();
+
+ }//pool NULL check
+ }//run mixing
+
+ //Calculate the balance function
+ fBalance->CalculateBalance(fReactionPlane,tracksMain,NULL,1,1.,0.);
+ //if(fRunMixing)
+ // fMixedBalance->CalculateBalance(fReactionPlane,tracksMixing,NULL,1,1.,0.);
+
+ }//event loop
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskToyModel::FinishOutput() {
+ //Printf("END BF");
+
+ TFile *gOutput = TFile::Open("outputToyModel.root","recreate");
+ TDirectoryFile *dir = new TDirectoryFile("PWGCFEbyE.outputBalanceFunctionPsiAnalysis","PWGCFEbyE.outputBalanceFunctionPsiAnalysis");
+
+ fList->SetName("listQA");
+ fList->SetOwner(kTRUE);
+ dir->Add(fList);
+ //fList->Write("listQA",TObject::kSingleKey);
+
+ if (!fBalance) {
+ AliError("ERROR: fBalance not available");
+ return;
+ }
+ fListBF->SetName("listBF");
+ fListBF->SetOwner(kTRUE);
+ dir->Add(fListBF);
+ //fListBF->Write("listBF",TObject::kSingleKey);
+
+ if(fRunShuffling) {
+ if (!fShuffledBalance) {
+ AliError("ERROR: fShuffledBalance not available");
+ return;
+ }
+ fListBFS->SetName("listBFShuffled");
+ fListBFS->SetOwner(kTRUE);
+ dir->Add(fListBFS);
+ //fListBFS->Write("listBFShuffled",TObject::kSingleKey);
+ }
+
+ if(fRunMixing) {
+ if (!fMixedBalance) {
+ AliError("ERROR: fMixedBalance not available");
+ return;
+ }
+ fListBFM->SetName("listBFMixed");
+ fListBFM->SetOwner(kTRUE);
+ dir->Add(fListBFM);
+ //fListBFM->Write("listBFMixed",TObject::kSingleKey);
+ }
+
+ dir->Write(dir->GetName(),TObject::kSingleKey);
+ gOutput->Close();
+}
+
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff