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