]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGCF/EBYE/BalanceFunctions/AliAnalysisTaskToyModel.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Simulating detector effects (more later)
[u/mrichter/AliRoot.git] / PWGCF / EBYE / BalanceFunctions / AliAnalysisTaskToyModel.cxx
CommitLineData
63d15f4b 1#include "TChain.h"\r
2#include "TList.h"\r
3#include "TCanvas.h"\r
4#include "TParticle.h"\r
5#include "TLorentzVector.h"\r
6#include "TGraphErrors.h"\r
59d226d3 7#include "TH1.h"\r
8#include "TH2.h"\r
9#include "TH3.h"\r
63d15f4b 10#include "TArrayF.h"\r
11#include "TF1.h"\r
12#include "TRandom.h"\r
0ed078e1 13#include "TFile.h"\r
63d15f4b 14\r
63d15f4b 15#include "AliAnalysisManager.h"\r
16#include "AliLog.h"\r
17\r
18#include "AliAnalysisTaskToyModel.h"\r
19#include "AliBalance.h"\r
c5914570 20#include "AliBalancePsi.h"\r
21#include "AliAnalysisTaskTriggeredBF.h"\r
63d15f4b 22\r
23\r
24// Analysis task for the toy model analysis\r
25// Authors: Panos.Christakoglou@nikhef.nl\r
26\r
c64cb1f6 27using std::cout;\r
28using std::endl;\r
29\r
63d15f4b 30ClassImp(AliAnalysisTaskToyModel)\r
31\r
32//________________________________________________________________________\r
0ed078e1 33AliAnalysisTaskToyModel::AliAnalysisTaskToyModel() \r
34: TObject(),\r
6b3bc65c 35 fUseDebug(kFALSE),\r
36 fBalance(0),\r
37 fRunShuffling(kFALSE), fShuffledBalance(0),\r
38 fList(0), fListBF(0), fListBFS(0),\r
39 fHistEventStats(0),\r
40 fHistNumberOfAcceptedParticles(0),\r
41 fHistReactionPlane(0),\r
42 fHistEtaTotal(0), fHistEta(0),\r
c5914570 43 fHistEtaPhiPos(0), fHistEtaPhiNeg(0),\r
6b3bc65c 44 fHistRapidity(0),\r
45 fHistRapidityPions(0), fHistRapidityKaons(0), fHistRapidityProtons(0),\r
46 fHistPhi(0),\r
47 fHistPhiPions(0), fHistPhiKaons(0), fHistPhiProtons(0),\r
48 fHistPt(0),\r
49 fHistPtPions(0), fHistPtKaons(0), fHistPtProtons(0),\r
50 fTotalMultiplicityMean(100.), fTotalMultiplicitySigma(0.0),\r
51 fNetChargeMean(0.0), fNetChargeSigma(0.0),\r
52 fPtMin(0.0), fPtMax(100.0),\r
53 fEtaMin(-1.0), fEtaMax(1.0),\r
54 fUseAcceptanceParameterization(kFALSE), fAcceptanceParameterization(0),\r
59d226d3 55 fSimulateDetectorEffects(kFALSE),\r
56 fNumberOfInefficientSectors(0),\r
57 fInefficiencyFactorInPhi(1.0),\r
58 fNumberOfDeadSectors(0),\r
59 fEfficiencyDropNearEtaEdges(kFALSE),\r
60 fEfficiencyMatrix(0),\r
6b3bc65c 61 fUseAllCharges(kFALSE), fParticleMass(0.0),\r
62 fPtSpectraAllCharges(0), fTemperatureAllCharges(100.),\r
63 fReactionPlane(0.0),\r
64 fAzimuthalAngleAllCharges(0), fDirectedFlowAllCharges(0.0), \r
65 fEllipticFlowAllCharges(0.0), fTriangularFlowAllCharges(0.0),\r
66 fQuandrangularFlowAllCharges(0.0), fPentangularFlowAllCharges(0.0),\r
67 fPionPercentage(0.8), fPionMass(0.0),\r
68 fPtSpectraPions(0), fTemperaturePions(100.),\r
69 fAzimuthalAnglePions(0), fDirectedFlowPions(0.0), \r
70 fEllipticFlowPions(0.0), fTriangularFlowPions(0.0), \r
71 fQuandrangularFlowPions(0.0), fPentangularFlowPions(0.0),\r
72 fKaonPercentage(0.8), fKaonMass(0.0),\r
73 fPtSpectraKaons(0), fTemperatureKaons(100.),\r
74 fAzimuthalAngleKaons(0), fDirectedFlowKaons(0.0), \r
75 fEllipticFlowKaons(0.0), fTriangularFlowKaons(0.0),\r
76 fQuandrangularFlowKaons(0.0), fPentangularFlowKaons(0.0),\r
77 fProtonPercentage(0.8), fProtonMass(0.0),\r
78 fPtSpectraProtons(0), fTemperatureProtons(100.),\r
79 fAzimuthalAngleProtons(0), fDirectedFlowProtons(0.0), \r
80 fEllipticFlowProtons(0.0), fTriangularFlowProtons(0.0),\r
81 fQuandrangularFlowProtons(0.0), fPentangularFlowProtons(0.0),\r
82 fUseDynamicalCorrelations(kFALSE), fDynamicalCorrelationsPercentage(0.1) {\r
63d15f4b 83 // Constructor\r
63d15f4b 84}\r
85\r
86//________________________________________________________________________\r
87AliAnalysisTaskToyModel::~AliAnalysisTaskToyModel() {\r
88 //Destructor\r
89 delete fPtSpectraAllCharges;\r
90 delete fAzimuthalAngleAllCharges;\r
91 delete fPtSpectraPions;\r
92 delete fAzimuthalAnglePions;\r
93 delete fPtSpectraKaons;\r
94 delete fAzimuthalAngleKaons;\r
95 delete fPtSpectraProtons;\r
96 delete fAzimuthalAngleProtons;\r
97}\r
98\r
99//________________________________________________________________________\r
100void AliAnalysisTaskToyModel::Init() {\r
101 //Initialize objects\r
78a6252b 102 //==============gRandom Seed=======================//\r
103 gRandom->SetSeed(0); //seed is set to a random value (depending on machine clock)\r
104 //==============gRandom Seed=======================//\r
105\r
63d15f4b 106 //==============Particles and spectra==============//\r
107 TParticle *pion = new TParticle();\r
108 pion->SetPdgCode(211);\r
b4665e0f 109 fPionMass = pion->GetMass();\r
63d15f4b 110\r
111 TParticle *kaon = new TParticle();\r
112 kaon->SetPdgCode(321);\r
b4665e0f 113 fKaonMass = kaon->GetMass();\r
114 \r
63d15f4b 115 TParticle *proton = new TParticle();\r
116 proton->SetPdgCode(2212);\r
b4665e0f 117 fProtonMass = proton->GetMass();\r
118\r
119 if(fUseAllCharges) {\r
120 fParticleMass = fPionMass;\r
c5914570 121 //fPtSpectraAllCharges = new TF1("fPtSpectraAllCharges","x*TMath::Exp(-TMath::Power([0]*[0]+x*x,0.5)/[1])",0.,5.);\r
122 //fPtSpectraAllCharges->SetParName(0,"Mass");\r
123 //fPtSpectraAllCharges->SetParName(1,"Temperature");\r
124 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
125 fPtSpectraAllCharges->SetParName(0,"pt0");\r
126 fPtSpectraAllCharges->SetParName(1,"b");\r
b4665e0f 127 }\r
128 else {\r
0ed078e1 129 fPtSpectraPions = new TF1("fPtSpectraPions","x*TMath::Exp(-TMath::Power([0]*[0]+x*x,0.5)/[1])",0.,5.);\r
b4665e0f 130 fPtSpectraPions->SetParName(0,"Mass");\r
b4665e0f 131 fPtSpectraPions->SetParName(1,"Temperature");\r
b4665e0f 132 \r
0ed078e1 133 fPtSpectraKaons = new TF1("fPtSpectraKaons","x*TMath::Exp(-TMath::Power([0]*[0]+x*x,0.5)/[1])",0.,5.);\r
b4665e0f 134 fPtSpectraKaons->SetParName(0,"Mass");\r
b4665e0f 135 fPtSpectraKaons->SetParName(1,"Temperature");\r
b4665e0f 136 \r
137 fPtSpectraProtons = new TF1("fPtSpectraProtons","x*TMath::Exp(-TMath::Power([0]*[0]+x*x,0.5)/[1])",0.,5.);\r
138 fPtSpectraProtons->SetParName(0,"Mass");\r
b4665e0f 139 fPtSpectraProtons->SetParName(1,"Temperature");\r
b4665e0f 140 }\r
63d15f4b 141 //==============Particles and spectra==============//\r
142\r
143 //==============Flow values==============//\r
b4665e0f 144 if(fUseAllCharges) {\r
0ed078e1 145 if(fUseDebug) {\r
146 Printf("Directed flow: %lf",fDirectedFlowAllCharges);\r
147 Printf("Elliptic flow: %lf",fEllipticFlowAllCharges);\r
148 Printf("Triangular flow: %lf",fTriangularFlowAllCharges);\r
149 Printf("Quandrangular flow: %lf",fQuandrangularFlowAllCharges);\r
150 Printf("Pentangular flow: %lf",fPentangularFlowAllCharges);\r
151 }\r
152\r
b4665e0f 153 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
154 fAzimuthalAngleAllCharges->SetParName(0,"Reaction Plane");\r
155 fAzimuthalAngleAllCharges->SetParName(1,"Directed flow");\r
b4665e0f 156 fAzimuthalAngleAllCharges->SetParName(2,"Elliptic flow"); \r
b4665e0f 157 fAzimuthalAngleAllCharges->SetParName(3,"Triangular flow");\r
b4665e0f 158 fAzimuthalAngleAllCharges->SetParName(4,"Quandrangular flow");\r
b4665e0f 159 fAzimuthalAngleAllCharges->SetParName(5,"Pentangular flow");\r
b4665e0f 160 }\r
161 else {\r
162 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
163 fAzimuthalAnglePions->SetParName(0,"Reaction Plane");\r
164 fAzimuthalAnglePions->SetParName(1,"Directed flow");\r
b4665e0f 165 fAzimuthalAnglePions->SetParName(2,"Elliptic flow"); \r
b4665e0f 166 fAzimuthalAnglePions->SetParName(3,"Triangular flow");\r
b4665e0f 167 fAzimuthalAnglePions->SetParName(4,"Quandrangular flow");\r
b4665e0f 168 fAzimuthalAnglePions->SetParName(5,"Pentangular flow");\r
b4665e0f 169 \r
170 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
171 fAzimuthalAngleKaons->SetParName(0,"Reaction Plane");\r
172 fAzimuthalAngleKaons->SetParName(1,"Directed flow");\r
b4665e0f 173 fAzimuthalAngleKaons->SetParName(2,"Elliptic flow"); \r
b4665e0f 174 fAzimuthalAngleKaons->SetParName(3,"Triangular flow");\r
b4665e0f 175 fAzimuthalAngleKaons->SetParName(4,"Quandrangular flow");\r
b4665e0f 176 fAzimuthalAngleKaons->SetParName(5,"Pentangular flow");\r
b4665e0f 177 \r
178 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
179 fAzimuthalAngleProtons->SetParName(0,"Reaction Plane");\r
180 fAzimuthalAngleProtons->SetParName(1,"Directed flow");\r
b4665e0f 181 fAzimuthalAngleProtons->SetParName(2,"Elliptic flow"); \r
b4665e0f 182 fAzimuthalAngleProtons->SetParName(3,"Triangular flow");\r
b4665e0f 183 fAzimuthalAngleProtons->SetParName(4,"Quandrangular flow");\r
b4665e0f 184 fAzimuthalAngleProtons->SetParName(5,"Pentangular flow");\r
b4665e0f 185 }\r
63d15f4b 186 //==============Flow values==============//\r
59d226d3 187\r
188 //==============Efficiency matrix==============//\r
189 if(fSimulateDetectorEffects) SetupEfficiencyMatrix();\r
190 //==============Efficiency matrix==============//\r
191}\r
192\r
193//________________________________________________________________________\r
194void AliAnalysisTaskToyModel::SetupEfficiencyMatrix() {\r
195 //Setup the efficiency matrix\r
196 TH1F *hPt = new TH1F("hPt","",200,0.1,20.1);\r
197 TH1F *hEta = new TH1F("hEta","",20,-0.95,0.95);\r
198 TH1F *hPhi = new TH1F("hPhi","",72,0.,2.*TMath::Pi());\r
199 fEfficiencyMatrix = new TH3F("fEfficiencyMatrix","",\r
200 hEta->GetNbinsX(),\r
201 hEta->GetXaxis()->GetXmin(),\r
202 hEta->GetXaxis()->GetXmax(),\r
203 hPt->GetNbinsX(),\r
204 hPt->GetXaxis()->GetXmin(),\r
205 hPt->GetXaxis()->GetXmax(),\r
206 hPhi->GetNbinsX(),\r
207 hPhi->GetXaxis()->GetXmin(),\r
208 hPhi->GetXaxis()->GetXmax());\r
209\r
210 //Efficiency in pt\r
211 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
212 for(Int_t i=1;i<=20;i++) {\r
213 hPt->SetBinContent(i,epsilon[i-1]);\r
214 hPt->SetBinError(i,0.01);\r
215 }\r
216 for(Int_t i=21;i<=200;i++) {\r
217 hPt->SetBinContent(i,epsilon[19]);\r
218 hPt->SetBinError(i,0.01);\r
219 }\r
220\r
221 //Efficiency in eta\r
222 for(Int_t i=1;i<=hEta->GetNbinsX();i++) {\r
223 hEta->SetBinContent(i,1.0);\r
224 hEta->SetBinError(i,0.01);\r
225 }\r
226 if(fEfficiencyDropNearEtaEdges) {\r
227 hEta->SetBinContent(1,0.7); hEta->SetBinContent(2,0.8);\r
228 hEta->SetBinContent(3,0.9);\r
229 hEta->SetBinContent(18,0.9); hEta->SetBinContent(19,0.8);\r
230 hEta->SetBinContent(20,0.7);\r
231 }\r
232\r
233 //Efficiency in phi\r
234 for(Int_t i=1;i<=hPhi->GetNbinsX();i++) {\r
235 hPhi->SetBinContent(i,1.0);\r
236 hPhi->SetBinError(i,0.01);\r
237 }\r
238 for(Int_t i=1;i<=fNumberOfInefficientSectors;i++)\r
239 hPhi->SetBinContent(hPhi->FindBin(hPhi->GetRandom()),fInefficiencyFactorInPhi);\r
240 for(Int_t i=1;i<=fNumberOfDeadSectors;i++)\r
241 hPhi->SetBinContent(hPhi->FindBin(hPhi->GetRandom()),0.0);\r
242 \r
243 //Fill the 3D efficiency map\r
244 for(Int_t iBinX = 1; iBinX<=fEfficiencyMatrix->GetXaxis()->GetNbins();iBinX++) {\r
245 //cout<<"==================================="<<endl;\r
246 for(Int_t iBinY = 1; iBinY<=fEfficiencyMatrix->GetYaxis()->GetNbins();iBinY++) {\r
247 //cout<<"==================================="<<endl;\r
248 for(Int_t iBinZ = 1; iBinZ<=fEfficiencyMatrix->GetZaxis()->GetNbins();iBinZ++) {\r
249 fEfficiencyMatrix->SetBinContent(iBinX,iBinY,iBinZ,hEta->GetBinContent(iBinX)*hPt->GetBinContent(iBinY)*hPhi->GetBinContent(iBinZ));\r
250 //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
251 }\r
252 }\r
253 }\r
63d15f4b 254}\r
255\r
256//________________________________________________________________________\r
0ed078e1 257void AliAnalysisTaskToyModel::CreateOutputObjects() {\r
63d15f4b 258 // Create histograms\r
259 // Called once\r
211b716d 260\r
261 // global switch disabling the reference \r
262 // (to avoid "Replacing existing TH1" if several wagons are created in train)\r
263 Bool_t oldStatus = TH1::AddDirectoryStatus();\r
264 TH1::AddDirectory(kFALSE);\r
265\r
63d15f4b 266 if(!fBalance) {\r
c5914570 267 fBalance = new AliBalancePsi();\r
268 fBalance->SetDeltaEtaMax(2.0);\r
269 //fBalance->SetInterval(-1,-0.8,0.8,16,0.,1.6);\r
63d15f4b 270 }\r
271 if(fRunShuffling) {\r
272 if(!fShuffledBalance) {\r
c5914570 273 fShuffledBalance = new AliBalancePsi();\r
274 fShuffledBalance->SetDeltaEtaMax(2.0);\r
275 //fShuffledBalance->SetInterval(-1,-0.8,0.8,16,0.,1.6);\r
63d15f4b 276 }\r
277 }\r
278\r
279 //QA list\r
280 fList = new TList();\r
281 fList->SetName("listQA");\r
282 fList->SetOwner();\r
283\r
284 //Balance Function list\r
285 fListBF = new TList();\r
286 fListBF->SetName("listBF");\r
287 fListBF->SetOwner();\r
288\r
289 if(fRunShuffling) {\r
290 fListBFS = new TList();\r
291 fListBFS->SetName("listBFShuffled");\r
292 fListBFS->SetOwner();\r
293 }\r
294\r
0ed078e1 295 //==============QA================//\r
63d15f4b 296 //Event stats.\r
297 TString gCutName[4] = {"Total","Offline trigger",\r
298 "Vertex","Analyzed"};\r
299 fHistEventStats = new TH1F("fHistEventStats",\r
300 "Event statistics;;N_{events}",\r
301 4,0.5,4.5);\r
302 for(Int_t i = 1; i <= 4; i++)\r
303 fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data());\r
304 fList->Add(fHistEventStats);\r
305\r
0ed078e1 306 fHistNumberOfAcceptedParticles = new TH1F("fHistNumberOfAcceptedParticles",";N_{acc.};Entries",10000,0,10000);\r
307 fList->Add(fHistNumberOfAcceptedParticles);\r
308 \r
309 fHistReactionPlane = new TH1F("fHistReactionPlane","Reaction plane angle;#Psi [rad];Entries",1000,0.,2.*TMath::Pi());\r
310 fList->Add(fHistReactionPlane);\r
311\r
312 //Pseudo-rapidity\r
313 fHistEtaTotal = new TH1F("fHistEtaTotal","Pseudo-rapidity (full phase space);#eta;Entries",1000,-15.,15.); \r
314 fList->Add(fHistEtaTotal);\r
315\r
316 fHistEta = new TH1F("fHistEta","Pseudo-rapidity (acceptance);#eta;Entries",1000,-1.5,1.5); \r
317 fList->Add(fHistEta);\r
318\r
c5914570 319 fHistEtaPhiPos = new TH2F("fHistEtaPhiPos","#eta-#phi distribution (+);#eta;#varphi (rad)",80,-2.,2.,72,0.,2.*TMath::Pi());\r
320 fList->Add(fHistEtaPhiPos);\r
321 fHistEtaPhiNeg = new TH2F("fHistEtaPhiNeg","#eta-#phi distribution (-);#eta;#varphi (rad)",80,-2.,2.,72,0.,2.*TMath::Pi());\r
322 fList->Add(fHistEtaPhiNeg);\r
323\r
0ed078e1 324 //Rapidity\r
0ed078e1 325 fHistRapidity = new TH1F("fHistRapidity","Rapidity (acceptance);y;Entries",1000,-1.5,1.5); \r
326 fList->Add(fHistRapidity);\r
327 fHistRapidityPions = new TH1F("fHistRapidityPions","Rapidity (acceptance - pions);y;Entries",1000,-1.5,1.5); \r
328 fList->Add(fHistRapidityPions);\r
329 fHistRapidityKaons = new TH1F("fHistRapidityKaons","Rapidity (acceptance - kaons);y;Entries",1000,-1.5,1.5); \r
330 fList->Add(fHistRapidityKaons);\r
331 fHistRapidityProtons = new TH1F("fHistRapidityProtons","Rapidity (acceptance - protons);y;Entries",1000,-1.5,1.5); \r
332 fList->Add(fHistRapidityProtons);\r
333\r
334 //Phi\r
0ed078e1 335 fHistPhi = new TH1F("fHistPhi","Phi (acceptance);#phi (rad);Entries",1000,0.,2*TMath::Pi());\r
336 fList->Add(fHistPhi);\r
337\r
338 fHistPhiPions = new TH1F("fHistPhiPions","Phi (acceptance - pions);#phi (rad);Entries",1000,0.,2*TMath::Pi());\r
339 fList->Add(fHistPhiPions);\r
340 fHistPhiKaons = new TH1F("fHistPhiKaons","Phi (acceptance - kaons);#phi (rad);Entries",1000,0.,2*TMath::Pi());\r
341 fList->Add(fHistPhiKaons);\r
342 fHistPhiProtons = new TH1F("fHistPhiProtons","Phi (acceptance - protons);#phi (rad);Entries",1000,0.,2*TMath::Pi());\r
343 fList->Add(fHistPhiProtons);\r
344\r
345\r
346 //Pt\r
0ed078e1 347 fHistPt = new TH1F("fHistPt","Pt (acceptance);p_{t} (GeV/c);Entries",1000,0.,10.);\r
348 fList->Add(fHistPt);\r
349\r
350 fHistPtPions = new TH1F("fHistPtPions","Pt (acceptance - pions);p_{t} (GeV/c);Entries",1000,0.,10.);\r
351 fList->Add(fHistPtPions);\r
352 fHistPtKaons = new TH1F("fHistPtKaons","Pt (acceptance - kaons);p_{t} (GeV/c);Entries",1000,0.,10.);\r
353 fList->Add(fHistPtKaons);\r
354 fHistPtProtons = new TH1F("fHistPtProtons","Pt (acceptance - protons);p_{t} (GeV/c);Entries",1000,0.,10.);\r
355 fList->Add(fHistPtProtons);\r
356\r
357 //==============Balance function histograms================//\r
63d15f4b 358 // Initialize histograms if not done yet\r
c5914570 359 if(!fBalance->GetHistNp()){\r
63d15f4b 360 AliWarning("Histograms not yet initialized! --> Will be done now");\r
361 AliWarning("--> Add 'gBalance->InitHistograms()' in your configBalanceFunction");\r
362 fBalance->InitHistograms();\r
363 }\r
364\r
365 if(fRunShuffling) {\r
c5914570 366 if(!fShuffledBalance->GetHistNp()) {\r
63d15f4b 367 AliWarning("Histograms (shuffling) not yet initialized! --> Will be done now");\r
368 AliWarning("--> Add 'gBalance->InitHistograms()' in your configBalanceFunction");\r
369 fShuffledBalance->InitHistograms();\r
370 }\r
371 }\r
372\r
c5914570 373 fListBF->Add(fBalance->GetHistNp());\r
374 fListBF->Add(fBalance->GetHistNn());\r
375 fListBF->Add(fBalance->GetHistNpn());\r
376 fListBF->Add(fBalance->GetHistNnn());\r
377 fListBF->Add(fBalance->GetHistNpp());\r
378 fListBF->Add(fBalance->GetHistNnp());\r
379\r
380 if(fRunShuffling) {\r
381 fListBFS->Add(fShuffledBalance->GetHistNp());\r
382 fListBFS->Add(fShuffledBalance->GetHistNn());\r
383 fListBFS->Add(fShuffledBalance->GetHistNpn());\r
384 fListBFS->Add(fShuffledBalance->GetHistNnn());\r
385 fListBFS->Add(fShuffledBalance->GetHistNpp());\r
386 fListBFS->Add(fShuffledBalance->GetHistNnp());\r
63d15f4b 387 }\r
388\r
389 // Post output data.\r
0ed078e1 390 //PostData(1, fList);\r
391 //PostData(2, fListBF);\r
392 //if(fRunShuffling) PostData(3, fListBFS);\r
211b716d 393\r
394 TH1::AddDirectory(oldStatus);\r
395\r
63d15f4b 396}\r
397\r
398//________________________________________________________________________\r
0ed078e1 399void AliAnalysisTaskToyModel::Run(Int_t nEvents) {\r
63d15f4b 400 // Main loop\r
401 // Called for each event\r
c5914570 402 Short_t vCharge = 0;\r
403 Float_t vY = 0.0;\r
404 Float_t vEta = 0.0;\r
405 Float_t vPhi = 0.0;\r
406 Float_t vP[3] = {0.,0.,0.};\r
407 Float_t vPt = 0.0;\r
408 Float_t vE = 0.0;\r
0ed078e1 409 Bool_t isPion = kFALSE, isKaon = kFALSE, isProton = kFALSE;\r
63d15f4b 410\r
0ed078e1 411 if(fUseAllCharges) {\r
c5914570 412 //fPtSpectraAllCharges->SetParameter(0,fParticleMass);\r
413 fPtSpectraAllCharges->SetParameter(0,1.05);\r
0ed078e1 414 fPtSpectraAllCharges->SetParameter(1,fTemperatureAllCharges);\r
b4665e0f 415\r
0ed078e1 416 fAzimuthalAngleAllCharges->SetParameter(1,fDirectedFlowAllCharges);\r
417 fAzimuthalAngleAllCharges->SetParameter(2,fEllipticFlowAllCharges);\r
418 fAzimuthalAngleAllCharges->SetParameter(3,fTriangularFlowAllCharges);\r
419 fAzimuthalAngleAllCharges->SetParameter(4,fQuandrangularFlowAllCharges);\r
420 fAzimuthalAngleAllCharges->SetParameter(5,fPentangularFlowAllCharges);\r
421 }\r
b4665e0f 422 else {\r
0ed078e1 423 fPtSpectraPions->SetParameter(0,fPionMass);\r
424 fPtSpectraPions->SetParameter(1,fTemperaturePions);\r
425 fPtSpectraKaons->SetParameter(0,fKaonMass);\r
426 fPtSpectraKaons->SetParameter(1,fTemperatureKaons);\r
427 fPtSpectraProtons->SetParameter(0,fProtonMass);\r
428 fPtSpectraProtons->SetParameter(1,fTemperatureProtons);\r
429\r
430 fAzimuthalAnglePions->SetParameter(1,fDirectedFlowPions);\r
431 fAzimuthalAnglePions->SetParameter(2,fEllipticFlowPions);\r
432 fAzimuthalAnglePions->SetParameter(3,fTriangularFlowPions);\r
433 fAzimuthalAnglePions->SetParameter(4,fQuandrangularFlowPions);\r
434 fAzimuthalAnglePions->SetParameter(5,fPentangularFlowPions);\r
435\r
436 fAzimuthalAngleKaons->SetParameter(1,fDirectedFlowKaons);\r
437 fAzimuthalAngleKaons->SetParameter(2,fEllipticFlowKaons);\r
438 fAzimuthalAngleKaons->SetParameter(3,fTriangularFlowKaons);\r
439 fAzimuthalAngleKaons->SetParameter(4,fQuandrangularFlowKaons);\r
440 fAzimuthalAngleKaons->SetParameter(5,fPentangularFlowKaons);\r
441\r
442 fAzimuthalAngleProtons->SetParameter(1,fDirectedFlowProtons);\r
443 fAzimuthalAngleProtons->SetParameter(2,fEllipticFlowProtons);\r
444 fAzimuthalAngleProtons->SetParameter(3,fTriangularFlowProtons);\r
445 fAzimuthalAngleProtons->SetParameter(4,fQuandrangularFlowProtons);\r
446 fAzimuthalAngleProtons->SetParameter(5,fPentangularFlowProtons);\r
b4665e0f 447 }\r
63d15f4b 448\r
c5914570 449 //TObjArray for the accepted particles\r
450 TObjArray *tracksMain = new TObjArray();\r
451 tracksMain->SetOwner(kTRUE);\r
452\r
0ed078e1 453 for(Int_t iEvent = 0; iEvent < nEvents; iEvent++) {\r
454 // vector holding the charges/kinematics of all tracks (charge,y,eta,phi,p0,p1,p2,pt,E)\r
c5914570 455 //vector<Double_t> *chargeVectorShuffle[9]; // this will be shuffled\r
456 //vector<Double_t> *chargeVector[9]; // original charge\r
457 //for(Int_t i = 0; i < 9; i++){\r
458 //chargeVectorShuffle[i] = new vector<Double_t>;\r
459 //chargeVector[i] = new vector<Double_t>;\r
460 //}\r
461 tracksMain->Clear();\r
462 \r
0ed078e1 463 fHistEventStats->Fill(1);\r
464 fHistEventStats->Fill(2);\r
465 fHistEventStats->Fill(3);\r
466\r
467 if((iEvent%1000) == 0) \r
468 cout<<"Event: "<<iEvent<<"/"<<nEvents<<endl;\r
469 \r
470 //Multiplicities\r
471 Int_t nMultiplicity = (Int_t)(gRandom->Gaus(fTotalMultiplicityMean,fTotalMultiplicitySigma));\r
472 Int_t nNetCharge = (Int_t)(gRandom->Gaus(fNetChargeMean,fNetChargeSigma));\r
6b3bc65c 473 \r
474 Int_t nGeneratedPositive = (Int_t)((nMultiplicity/2) + nNetCharge);\r
475 Int_t nGeneratedNegative = nMultiplicity - nGeneratedPositive;\r
476 if(fUseDebug) \r
477 Printf("Total multiplicity: %d - Generated positive: %d - Generated negative: %d",nMultiplicity,nGeneratedPositive,nGeneratedNegative);\r
478\r
479 //Int_t nGeneratedPositive = 0, nGeneratedNegative = 0;\r
0ed078e1 480 Int_t nGeneratedPions = 0, nGeneratedKaons = 0, nGeneratedProtons = 0;\r
481 \r
482 //Randomization of the reaction plane\r
c5914570 483 fReactionPlane = 2.0*TMath::Pi()*gRandom->Rndm();\r
484 //fReactionPlane = 0.0;\r
0ed078e1 485 if(fUseAllCharges) \r
486 fAzimuthalAngleAllCharges->SetParameter(0,fReactionPlane);\r
b4665e0f 487 else {\r
0ed078e1 488 fAzimuthalAnglePions->SetParameter(0,fReactionPlane);\r
489 fAzimuthalAngleKaons->SetParameter(0,fReactionPlane);\r
490 fAzimuthalAngleProtons->SetParameter(0,fReactionPlane);\r
63d15f4b 491 }\r
0ed078e1 492 \r
493 Int_t gNumberOfAcceptedParticles = 0;\r
6b3bc65c 494 Int_t gNumberOfAcceptedPositiveParticles = 0;\r
495 Int_t gNumberOfAcceptedNegativeParticles = 0;\r
496 \r
497 //Generate positive particles\r
498 for(Int_t iParticleCount = 0; iParticleCount < nGeneratedPositive; iParticleCount++) {\r
0ed078e1 499 isPion = kFALSE; isKaon = kFALSE; isProton = kFALSE;\r
6b3bc65c 500 if(fUseDebug) \r
501 Printf("Generating positive: %d(%d)",iParticleCount+1,nGeneratedPositive);\r
0ed078e1 502\r
6b3bc65c 503 //Pseudo-rapidity sampled from a Gaussian centered @ 0\r
9fd4b54e 504 vEta = gRandom->Gaus(0.0,4.0);\r
6b3bc65c 505\r
506 //Fill QA histograms (full phase space)\r
9fd4b54e 507 fHistEtaTotal->Fill(vEta);\r
6b3bc65c 508\r
c5914570 509 vCharge = 1;\r
6b3bc65c 510 //nGeneratedPositive += 1;\r
0ed078e1 511 \r
6b3bc65c 512 //Acceptance\r
9fd4b54e 513 if((vEta < fEtaMin) || (vEta > fEtaMax)) continue;\r
6b3bc65c 514\r
0ed078e1 515 if(!fUseAllCharges) {\r
516 //Decide the specie\r
517 Double_t randomNumberSpecies = gRandom->Rndm();\r
518 if((randomNumberSpecies >= 0.0)&&(randomNumberSpecies < fPionPercentage)) {\r
519 nGeneratedPions += 1;\r
9fd4b54e 520 vPt = fPtSpectraPions->GetRandom();\r
521 vPhi = fAzimuthalAnglePions->GetRandom();\r
0ed078e1 522 fParticleMass = fPionMass;\r
523 isPion = kTRUE;\r
524 }\r
525 else if((randomNumberSpecies >= fPionPercentage)&&(randomNumberSpecies < fPionPercentage + fKaonPercentage)) {\r
526 nGeneratedKaons += 1;\r
9fd4b54e 527 vPt = fPtSpectraKaons->GetRandom();\r
528 vPhi = fAzimuthalAngleKaons->GetRandom();\r
0ed078e1 529 fParticleMass = fKaonMass;\r
530 isKaon = kTRUE;\r
531 }\r
532 else if((randomNumberSpecies >= fPionPercentage + fKaonPercentage)&&(randomNumberSpecies < fPionPercentage + fKaonPercentage + fProtonPercentage)) {\r
533 nGeneratedProtons += 1;\r
9fd4b54e 534 vPt = fPtSpectraProtons->GetRandom();\r
535 vPhi = fAzimuthalAngleProtons->GetRandom();\r
0ed078e1 536 fParticleMass = fProtonMass;\r
537 isProton = kTRUE;\r
538 }\r
b4665e0f 539 }\r
0ed078e1 540 else {\r
9fd4b54e 541 vPt = fPtSpectraAllCharges->GetRandom();\r
542 vPhi = fAzimuthalAngleAllCharges->GetRandom();\r
0ed078e1 543 }\r
0ed078e1 544 \r
9fd4b54e 545 vP[0] = vPt*TMath::Cos(vPhi);\r
546 vP[1] = vPt*TMath::Sin(vPhi);\r
547 vP[2] = vPt*TMath::SinH(vEta);\r
548 vE = TMath::Sqrt(TMath::Power(fParticleMass,2) +\r
549 TMath::Power(vP[0],2) +\r
550 TMath::Power(vP[1],2) +\r
551 TMath::Power(vP[2],2));\r
0ed078e1 552 \r
9fd4b54e 553 vY = 0.5*TMath::Log((vE + vP[2])/(vE - vP[2]));\r
0ed078e1 554 \r
6b3bc65c 555 //pt coverage\r
9fd4b54e 556 if((vPt < fPtMin) || (vPt > fPtMax)) continue;\r
557 //Printf("pt: %lf - mins: %lf - max: %lf",vPt,fPtMin,fPtMax);\r
6b3bc65c 558\r
559 //acceptance filter\r
560 if(fUseAcceptanceParameterization) {\r
561 Double_t gRandomNumberForAcceptance = gRandom->Rndm();\r
9fd4b54e 562 if(gRandomNumberForAcceptance > fAcceptanceParameterization->Eval(vPt)) \r
6b3bc65c 563 continue;\r
564 }\r
565\r
566 gNumberOfAcceptedPositiveParticles += 1;\r
567\r
568 //Fill QA histograms (acceptance);\r
c5914570 569 fHistEtaPhiPos->Fill(vEta,vPhi);\r
9fd4b54e 570 fHistEta->Fill(vEta);\r
571 fHistRapidity->Fill(vY);\r
572 fHistPhi->Fill(vPhi);\r
573 fHistPt->Fill(vPt);\r
6b3bc65c 574 if(isPion) {\r
9fd4b54e 575 fHistRapidityPions->Fill(vY);\r
576 fHistPhiPions->Fill(vPhi);\r
577 fHistPtPions->Fill(vPt);\r
6b3bc65c 578 }\r
579 else if(isKaon) {\r
9fd4b54e 580 fHistRapidityKaons->Fill(vY);\r
581 fHistPhiKaons->Fill(vPhi);\r
582 fHistPtKaons->Fill(vPt);\r
6b3bc65c 583 }\r
584 else if(isProton) {\r
9fd4b54e 585 fHistRapidityProtons->Fill(vY);\r
586 fHistPhiProtons->Fill(vPhi);\r
587 fHistPtProtons->Fill(vPt);\r
6b3bc65c 588 }\r
589\r
590 // fill charge vector\r
c5914570 591 /*chargeVector[0]->push_back(vCharge);\r
9fd4b54e 592 chargeVector[1]->push_back(vY);\r
593 chargeVector[2]->push_back(vEta);\r
c5914570 594 chargeVector[3]->push_back(vPhi);\r
9fd4b54e 595 chargeVector[4]->push_back(vP[0]);\r
596 chargeVector[5]->push_back(vP[1]);\r
597 chargeVector[6]->push_back(vP[2]);\r
598 chargeVector[7]->push_back(vPt);\r
599 chargeVector[8]->push_back(vE);\r
6b3bc65c 600 \r
601 if(fRunShuffling) {\r
9fd4b54e 602 chargeVectorShuffle[0]->push_back(vCharge);\r
603 chargeVectorShuffle[1]->push_back(vY);\r
604 chargeVectorShuffle[2]->push_back(vEta);\r
c5914570 605 chargeVectorShuffle[3]->push_back(vPhi);\r
9fd4b54e 606 chargeVectorShuffle[4]->push_back(vP[0]);\r
607 chargeVectorShuffle[5]->push_back(vP[1]);\r
608 chargeVectorShuffle[6]->push_back(vP[2]);\r
609 chargeVectorShuffle[7]->push_back(vPt);\r
610 chargeVectorShuffle[8]->push_back(vE);\r
c5914570 611 }*/\r
6b3bc65c 612 gNumberOfAcceptedParticles += 1;\r
c5914570 613\r
614 // add the track to the TObjArray\r
615 tracksMain->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, 1.0)); \r
6b3bc65c 616 }//generated positive particle loop\r
617 \r
618 //Generate negative particles\r
619 for(Int_t iParticleCount = 0; iParticleCount < nGeneratedNegative; iParticleCount++) {\r
620 isPion = kFALSE; isKaon = kFALSE; isProton = kFALSE;\r
621 if(fUseDebug) \r
622 Printf("Generating negative: %d(%d)",iParticleCount+1,nGeneratedNegative);\r
623\r
624 //Pseudo-rapidity sampled from a Gaussian centered @ 0\r
9fd4b54e 625 vEta = gRandom->Gaus(0.0,4.0);\r
6b3bc65c 626\r
0ed078e1 627 //Fill QA histograms (full phase space)\r
9fd4b54e 628 fHistEtaTotal->Fill(vEta);\r
0ed078e1 629\r
c5914570 630 vCharge = -1;\r
6b3bc65c 631 //nGeneratedNegative += 1;\r
632 \r
0ed078e1 633 //Acceptance\r
9fd4b54e 634 if((vEta < fEtaMin) || (vEta > fEtaMax)) continue;\r
6b3bc65c 635\r
636 if(!fUseAllCharges) {\r
637 //Decide the specie\r
638 Double_t randomNumberSpecies = gRandom->Rndm();\r
639 if((randomNumberSpecies >= 0.0)&&(randomNumberSpecies < fPionPercentage)) {\r
640 nGeneratedPions += 1;\r
9fd4b54e 641 vPt = fPtSpectraPions->GetRandom();\r
642 vPhi = fAzimuthalAnglePions->GetRandom();\r
6b3bc65c 643 fParticleMass = fPionMass;\r
644 isPion = kTRUE;\r
645 }\r
646 else if((randomNumberSpecies >= fPionPercentage)&&(randomNumberSpecies < fPionPercentage + fKaonPercentage)) {\r
647 nGeneratedKaons += 1;\r
9fd4b54e 648 vPt = fPtSpectraKaons->GetRandom();\r
649 vPhi = fAzimuthalAngleKaons->GetRandom();\r
6b3bc65c 650 fParticleMass = fKaonMass;\r
651 isKaon = kTRUE;\r
652 }\r
653 else if((randomNumberSpecies >= fPionPercentage + fKaonPercentage)&&(randomNumberSpecies < fPionPercentage + fKaonPercentage + fProtonPercentage)) {\r
654 nGeneratedProtons += 1;\r
9fd4b54e 655 vPt = fPtSpectraProtons->GetRandom();\r
656 vPhi = fAzimuthalAngleProtons->GetRandom();\r
6b3bc65c 657 fParticleMass = fProtonMass;\r
658 isProton = kTRUE;\r
659 }\r
660 }\r
661 else {\r
9fd4b54e 662 vPt = fPtSpectraAllCharges->GetRandom();\r
663 vPhi = fAzimuthalAngleAllCharges->GetRandom();\r
6b3bc65c 664 }\r
665 \r
9fd4b54e 666 vP[0] = vPt*TMath::Cos(vPhi);\r
667 vP[1] = vPt*TMath::Sin(vPhi);\r
668 vP[2] = vPt*TMath::SinH(vEta);\r
669 vE = TMath::Sqrt(TMath::Power(fParticleMass,2) +\r
670 TMath::Power(vP[0],2) +\r
671 TMath::Power(vP[1],2) +\r
672 TMath::Power(vP[2],2));\r
6b3bc65c 673 \r
9fd4b54e 674 vY = 0.5*TMath::Log((vE + vP[2])/(vE - vP[2]));\r
6b3bc65c 675 \r
676 //pt coverage\r
9fd4b54e 677 if((vPt < fPtMin) || (vPt > fPtMax)) continue;\r
678 //Printf("pt: %lf - mins: %lf - max: %lf",vPt,fPtMin,fPtMax);\r
0ed078e1 679\r
6b3bc65c 680 //acceptance filter\r
681 if(fUseAcceptanceParameterization) {\r
682 Double_t gRandomNumberForAcceptance = gRandom->Rndm();\r
9fd4b54e 683 if(gRandomNumberForAcceptance > fAcceptanceParameterization->Eval(vPt)) \r
6b3bc65c 684 continue;\r
685 }\r
686\r
687 gNumberOfAcceptedNegativeParticles += 1;\r
c5914570 688 \r
0ed078e1 689 //Fill QA histograms (acceptance);\r
c5914570 690 fHistEtaPhiNeg->Fill(vEta,vPhi);\r
9fd4b54e 691 fHistEta->Fill(vEta);\r
692 fHistRapidity->Fill(vY);\r
693 fHistPhi->Fill(vPhi);\r
694 fHistPt->Fill(vPt);\r
0ed078e1 695 if(isPion) {\r
9fd4b54e 696 fHistRapidityPions->Fill(vY);\r
697 fHistPhiPions->Fill(vPhi);\r
698 fHistPtPions->Fill(vPt);\r
0ed078e1 699 }\r
700 else if(isKaon) {\r
9fd4b54e 701 fHistRapidityKaons->Fill(vY);\r
702 fHistPhiKaons->Fill(vPhi);\r
703 fHistPtKaons->Fill(vPt);\r
0ed078e1 704 }\r
705 else if(isProton) {\r
9fd4b54e 706 fHistRapidityProtons->Fill(vY);\r
707 fHistPhiProtons->Fill(vPhi);\r
708 fHistPtProtons->Fill(vPt);\r
0ed078e1 709 }\r
710\r
711 // fill charge vector\r
c5914570 712 /*chargeVector[0]->push_back(vCharge);\r
9fd4b54e 713 chargeVector[1]->push_back(vY);\r
714 chargeVector[2]->push_back(vEta);\r
c5914570 715 chargeVector[3]->push_back(vPhi);\r
9fd4b54e 716 chargeVector[4]->push_back(vP[0]);\r
717 chargeVector[5]->push_back(vP[1]);\r
718 chargeVector[6]->push_back(vP[2]);\r
719 chargeVector[7]->push_back(vPt);\r
720 chargeVector[8]->push_back(vE);\r
0ed078e1 721 \r
722 if(fRunShuffling) {\r
9fd4b54e 723 chargeVectorShuffle[0]->push_back(vCharge);\r
724 chargeVectorShuffle[1]->push_back(vY);\r
725 chargeVectorShuffle[2]->push_back(vEta);\r
c5914570 726 chargeVectorShuffle[3]->push_back(vPhi);\r
9fd4b54e 727 chargeVectorShuffle[4]->push_back(vP[0]);\r
728 chargeVectorShuffle[5]->push_back(vP[1]);\r
729 chargeVectorShuffle[6]->push_back(vP[2]);\r
730 chargeVectorShuffle[7]->push_back(vPt);\r
731 chargeVectorShuffle[8]->push_back(vE);\r
c5914570 732 }*/\r
0ed078e1 733 gNumberOfAcceptedParticles += 1;\r
c5914570 734\r
735 // add the track to the TObjArray\r
736 tracksMain->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, 1.0)); \r
6b3bc65c 737 }//generated negative particle loop\r
78a6252b 738 \r
6b3bc65c 739 //Dynamical correlations\r
c5914570 740 Int_t nGeneratedPositiveDynamicalCorrelations = 0;\r
741 Int_t nGeneratedNegativeDynamicalCorrelations = 0;\r
742 /*Double_t vChargePrime = 0;\r
9fd4b54e 743 Double_t vYPrime = 0.0;\r
744 Double_t vEtaPrime = 0.0;\r
745 Double_t vPhiPrime = 0.0;\r
746 Double_t vPPrime[3] = {0.,0.,0.};\r
747 Double_t vPtPrime = 0.0;\r
748 Double_t vEPrime = 0.0;\r
6b3bc65c 749 //Generate "correlated" particles \r
750 if(fUseDynamicalCorrelations) {\r
751 Int_t gNumberOfDynamicalCorrelations = (Int_t)(0.5*gNumberOfAcceptedParticles*fDynamicalCorrelationsPercentage);\r
752 for(Int_t iDynamicalCorrelations = 0; iDynamicalCorrelations < gNumberOfDynamicalCorrelations; iDynamicalCorrelations++) {\r
753 isPion = kFALSE; isKaon = kFALSE; isProton = kFALSE;\r
754 \r
755 //Pseudo-rapidity sampled from a Gaussian centered @ 0\r
9fd4b54e 756 vEta = gRandom->Gaus(0.0,0.1);\r
757 vCharge = 1.0;\r
6b3bc65c 758 nGeneratedPositiveDynamicalCorrelations += 1;\r
759 \r
9fd4b54e 760 vEtaPrime = -vEta;\r
761 vChargePrime = -1.0;\r
6b3bc65c 762 nGeneratedNegativeDynamicalCorrelations += 1;\r
763 \r
764 //Acceptance\r
9fd4b54e 765 if((vEta < fEtaMin) || (vEta > fEtaMax)) continue;\r
766 if((vEtaPrime < fEtaMin) || (vEtaPrime > fEtaMax)) continue;\r
0ed078e1 767 \r
6b3bc65c 768 if(!fUseAllCharges) {\r
769 //Decide the specie\r
770 Double_t randomNumberSpecies = gRandom->Rndm();\r
771 if((randomNumberSpecies >= 0.0)&&(randomNumberSpecies < fPionPercentage)) {\r
772 nGeneratedPions += 1;\r
9fd4b54e 773 vPt = fPtSpectraPions->GetRandom();\r
774 vPhi = fAzimuthalAnglePions->GetRandom();\r
6b3bc65c 775 fParticleMass = fPionMass;\r
776 isPion = kTRUE;\r
777 }\r
778 else if((randomNumberSpecies >= fPionPercentage)&&(randomNumberSpecies < fPionPercentage + fKaonPercentage)) {\r
779 nGeneratedKaons += 1;\r
9fd4b54e 780 vPt = fPtSpectraKaons->GetRandom();\r
781 vPhi = fAzimuthalAngleKaons->GetRandom();\r
6b3bc65c 782 fParticleMass = fKaonMass;\r
783 isKaon = kTRUE;\r
784 }\r
785 else if((randomNumberSpecies >= fPionPercentage + fKaonPercentage)&&(randomNumberSpecies < fPionPercentage + fKaonPercentage + fProtonPercentage)) {\r
786 nGeneratedProtons += 1;\r
9fd4b54e 787 vPt = fPtSpectraProtons->GetRandom();\r
788 vPtPrime = vPt;\r
789 vPhi = fAzimuthalAngleProtons->GetRandom();\r
6b3bc65c 790 fParticleMass = fProtonMass;\r
791 isProton = kTRUE;\r
792 }\r
793 }\r
794 else {\r
9fd4b54e 795 vPt = fPtSpectraAllCharges->GetRandom();\r
796 vPhi = fAzimuthalAngleAllCharges->GetRandom();\r
6b3bc65c 797 }\r
9fd4b54e 798 vPtPrime = vPt;\r
799 vPhiPrime = vPhi;\r
800\r
801 vP[0] = vPt*TMath::Cos(vPhi);\r
802 vP[1] = vPt*TMath::Sin(vPhi);\r
803 vP[2] = vPt*TMath::SinH(vEta);\r
804 vE = TMath::Sqrt(TMath::Power(fParticleMass,2) +\r
805 TMath::Power(vP[0],2) +\r
806 TMath::Power(vP[1],2) +\r
807 TMath::Power(vP[2],2));\r
6b3bc65c 808 \r
9fd4b54e 809 vY = 0.5*TMath::Log((vE + vP[2])/(vE - vP[2]));\r
810\r
811 vPPrime[0] = vPtPrime*TMath::Cos(vPhiPrime);\r
812 vPPrime[1] = vPtPrime*TMath::Sin(vPhiPrime);\r
813 vPPrime[2] = vPtPrime*TMath::SinH(vEtaPrime);\r
814 vEPrime = TMath::Sqrt(TMath::Power(fParticleMass,2) +\r
815 TMath::Power(vPPrime[0],2) +\r
816 TMath::Power(vPPrime[1],2) +\r
817 TMath::Power(vPPrime[2],2));\r
6b3bc65c 818 \r
9fd4b54e 819 vYPrime = 0.5*TMath::Log((vEPrime + vPPrime[2])/(vEPrime - vPPrime[2]));\r
6b3bc65c 820 \r
821 //pt coverage\r
9fd4b54e 822 if((vPt < fPtMin) || (vPt > fPtMax)) continue;\r
823 if((vPtPrime < fPtMin) || (vPtPrime > fPtMax)) continue;\r
6b3bc65c 824\r
825 //acceptance filter\r
826 if(fUseAcceptanceParameterization) {\r
827 Double_t gRandomNumberForAcceptance = gRandom->Rndm();\r
9fd4b54e 828 if(gRandomNumberForAcceptance > fAcceptanceParameterization->Eval(vPt)) \r
6b3bc65c 829 continue;\r
830 \r
831 Double_t gRandomNumberForAcceptancePrime = gRandom->Rndm();\r
9fd4b54e 832 if(gRandomNumberForAcceptancePrime > fAcceptanceParameterization->Eval(vPtPrime)) \r
6b3bc65c 833 continue;\r
834 }\r
835 \r
836 // fill charge vector (positive)\r
9fd4b54e 837 chargeVector[0]->push_back(vCharge);\r
838 chargeVector[1]->push_back(vY);\r
839 chargeVector[2]->push_back(vEta);\r
c5914570 840 chargeVector[3]->push_back(vPhi);\r
9fd4b54e 841 chargeVector[4]->push_back(vP[0]);\r
842 chargeVector[5]->push_back(vP[1]);\r
843 chargeVector[6]->push_back(vP[2]);\r
844 chargeVector[7]->push_back(vPt);\r
845 chargeVector[8]->push_back(vE);\r
6b3bc65c 846 \r
847 if(fRunShuffling) {\r
9fd4b54e 848 chargeVectorShuffle[0]->push_back(vCharge);\r
849 chargeVectorShuffle[1]->push_back(vY);\r
850 chargeVectorShuffle[2]->push_back(vEta);\r
c5914570 851 chargeVectorShuffle[3]->push_back(vPhi);\r
9fd4b54e 852 chargeVectorShuffle[4]->push_back(vP[0]);\r
853 chargeVectorShuffle[5]->push_back(vP[1]);\r
854 chargeVectorShuffle[6]->push_back(vP[2]);\r
855 chargeVectorShuffle[7]->push_back(vPt);\r
856 chargeVectorShuffle[8]->push_back(vE);\r
6b3bc65c 857 }\r
858\r
859 // fill charge vector (negative)\r
9fd4b54e 860 chargeVector[0]->push_back(vChargePrime);\r
861 chargeVector[1]->push_back(vYPrime);\r
862 chargeVector[2]->push_back(vEtaPrime);\r
c5914570 863 chargeVector[3]->push_back(vPhiPrime);\r
9fd4b54e 864 chargeVector[4]->push_back(vPPrime[0]);\r
865 chargeVector[5]->push_back(vPPrime[1]);\r
866 chargeVector[6]->push_back(vPPrime[2]);\r
867 chargeVector[7]->push_back(vPtPrime);\r
868 chargeVector[8]->push_back(vEPrime);\r
6b3bc65c 869 \r
870 if(fRunShuffling) {\r
9fd4b54e 871 chargeVectorShuffle[0]->push_back(vChargePrime);\r
872 chargeVectorShuffle[1]->push_back(vYPrime);\r
873 chargeVectorShuffle[2]->push_back(vEtaPrime);\r
c5914570 874 chargeVectorShuffle[3]->push_back(vPhiPrime);\r
9fd4b54e 875 chargeVectorShuffle[4]->push_back(vPPrime[0]);\r
876 chargeVectorShuffle[5]->push_back(vPPrime[1]);\r
877 chargeVectorShuffle[6]->push_back(vPPrime[2]);\r
878 chargeVectorShuffle[7]->push_back(vPtPrime);\r
879 chargeVectorShuffle[8]->push_back(vEPrime);\r
6b3bc65c 880 }\r
881\r
882 gNumberOfAcceptedParticles += 2;\r
c5914570 883 }//loop over the dynamical correlations\r
884 }*/// usage of dynamical correlations\r
6b3bc65c 885\r
0ed078e1 886 if(fUseDebug) {\r
887 Printf("=======================================================");\r
888 Printf("Total: %d - Total positive: %d - Total negative: %d",nMultiplicity,nGeneratedPositive,nGeneratedNegative);\r
6b3bc65c 889 Printf("Accepted positive: %d - Accepted negative: %d",gNumberOfAcceptedPositiveParticles,gNumberOfAcceptedNegativeParticles);\r
890 Printf("Correlations: %d - Correlations positive: %d - Correlations negative: %d",nGeneratedPositiveDynamicalCorrelations+nGeneratedNegativeDynamicalCorrelations,nGeneratedPositiveDynamicalCorrelations,nGeneratedNegativeDynamicalCorrelations);\r
891 Printf("Number of accepted particles: %d",gNumberOfAcceptedParticles);\r
0ed078e1 892 if(!fUseAllCharges)\r
893 Printf("Pions: %lf - Kaons: %lf - Protons: %lf",1.*nGeneratedPions/nMultiplicity,1.*nGeneratedKaons/nMultiplicity,1.*nGeneratedProtons/nMultiplicity);\r
6b3bc65c 894 //Printf("Calculating the balance function for %d particles",chargeVector[0]->size());\r
63d15f4b 895 }\r
63d15f4b 896\r
0ed078e1 897 fHistEventStats->Fill(4);\r
898 fHistNumberOfAcceptedParticles->Fill(gNumberOfAcceptedParticles);\r
899 fHistReactionPlane->Fill(fReactionPlane);\r
63d15f4b 900\r
0ed078e1 901 //Calculate the balance function\r
c5914570 902 //fBalance->CalculateBalance(gNumberOfAcceptedParticles,chargeVector);\r
903 fBalance->CalculateBalance(fReactionPlane,tracksMain,NULL,1,1.,0.);\r
904 /*if(fRunShuffling) {\r
0ed078e1 905 // shuffle charges\r
906 random_shuffle( chargeVectorShuffle[0]->begin(), chargeVectorShuffle[0]->end() );\r
907 fShuffledBalance->CalculateBalance(gNumberOfAcceptedParticles,chargeVectorShuffle);\r
c5914570 908 }*/\r
909 }//event loop\r
63d15f4b 910} \r
911\r
912//________________________________________________________________________\r
0ed078e1 913void AliAnalysisTaskToyModel::FinishOutput() {\r
63d15f4b 914 //Printf("END BF");\r
915\r
0ed078e1 916 TFile *gOutput = TFile::Open("outputToyModel.root","recreate");\r
917 fList->Write();\r
918\r
63d15f4b 919 if (!fBalance) {\r
920 Printf("ERROR: fBalance not available");\r
921 return;\r
922 } \r
0ed078e1 923 fListBF->Write();\r
924\r
63d15f4b 925 if(fRunShuffling) {\r
926 if (!fShuffledBalance) {\r
927 Printf("ERROR: fShuffledBalance not available");\r
928 return;\r
929 }\r
0ed078e1 930 fListBFS->Write();\r
63d15f4b 931 }\r
0ed078e1 932 gOutput->Close();\r
63d15f4b 933}\r
934\r