]> 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
Updates in the balance function code
[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
211b716d 180\r
181 // global switch disabling the reference \r
182 // (to avoid "Replacing existing TH1" if several wagons are created in train)\r
183 Bool_t oldStatus = TH1::AddDirectoryStatus();\r
184 TH1::AddDirectory(kFALSE);\r
185\r
63d15f4b 186 if(!fBalance) {\r
187 fBalance = new AliBalance();\r
188 fBalance->SetInterval(-1,-0.8,0.8,16,0.,1.6);\r
189 }\r
190 if(fRunShuffling) {\r
191 if(!fShuffledBalance) {\r
192 fShuffledBalance = new AliBalance();\r
193 fShuffledBalance->SetInterval(-1,-0.8,0.8,16,0.,1.6);\r
194 }\r
195 }\r
196\r
197 //QA list\r
198 fList = new TList();\r
199 fList->SetName("listQA");\r
200 fList->SetOwner();\r
201\r
202 //Balance Function list\r
203 fListBF = new TList();\r
204 fListBF->SetName("listBF");\r
205 fListBF->SetOwner();\r
206\r
207 if(fRunShuffling) {\r
208 fListBFS = new TList();\r
209 fListBFS->SetName("listBFShuffled");\r
210 fListBFS->SetOwner();\r
211 }\r
212\r
0ed078e1 213 //==============QA================//\r
63d15f4b 214 //Event stats.\r
215 TString gCutName[4] = {"Total","Offline trigger",\r
216 "Vertex","Analyzed"};\r
217 fHistEventStats = new TH1F("fHistEventStats",\r
218 "Event statistics;;N_{events}",\r
219 4,0.5,4.5);\r
220 for(Int_t i = 1; i <= 4; i++)\r
221 fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data());\r
222 fList->Add(fHistEventStats);\r
223\r
0ed078e1 224 fHistNumberOfAcceptedParticles = new TH1F("fHistNumberOfAcceptedParticles",";N_{acc.};Entries",10000,0,10000);\r
225 fList->Add(fHistNumberOfAcceptedParticles);\r
226 \r
227 fHistReactionPlane = new TH1F("fHistReactionPlane","Reaction plane angle;#Psi [rad];Entries",1000,0.,2.*TMath::Pi());\r
228 fList->Add(fHistReactionPlane);\r
229\r
230 //Pseudo-rapidity\r
231 fHistEtaTotal = new TH1F("fHistEtaTotal","Pseudo-rapidity (full phase space);#eta;Entries",1000,-15.,15.); \r
232 fList->Add(fHistEtaTotal);\r
233\r
234 fHistEta = new TH1F("fHistEta","Pseudo-rapidity (acceptance);#eta;Entries",1000,-1.5,1.5); \r
235 fList->Add(fHistEta);\r
236\r
237 //Rapidity\r
0ed078e1 238 fHistRapidity = new TH1F("fHistRapidity","Rapidity (acceptance);y;Entries",1000,-1.5,1.5); \r
239 fList->Add(fHistRapidity);\r
240 fHistRapidityPions = new TH1F("fHistRapidityPions","Rapidity (acceptance - pions);y;Entries",1000,-1.5,1.5); \r
241 fList->Add(fHistRapidityPions);\r
242 fHistRapidityKaons = new TH1F("fHistRapidityKaons","Rapidity (acceptance - kaons);y;Entries",1000,-1.5,1.5); \r
243 fList->Add(fHistRapidityKaons);\r
244 fHistRapidityProtons = new TH1F("fHistRapidityProtons","Rapidity (acceptance - protons);y;Entries",1000,-1.5,1.5); \r
245 fList->Add(fHistRapidityProtons);\r
246\r
247 //Phi\r
0ed078e1 248 fHistPhi = new TH1F("fHistPhi","Phi (acceptance);#phi (rad);Entries",1000,0.,2*TMath::Pi());\r
249 fList->Add(fHistPhi);\r
250\r
251 fHistPhiPions = new TH1F("fHistPhiPions","Phi (acceptance - pions);#phi (rad);Entries",1000,0.,2*TMath::Pi());\r
252 fList->Add(fHistPhiPions);\r
253 fHistPhiKaons = new TH1F("fHistPhiKaons","Phi (acceptance - kaons);#phi (rad);Entries",1000,0.,2*TMath::Pi());\r
254 fList->Add(fHistPhiKaons);\r
255 fHistPhiProtons = new TH1F("fHistPhiProtons","Phi (acceptance - protons);#phi (rad);Entries",1000,0.,2*TMath::Pi());\r
256 fList->Add(fHistPhiProtons);\r
257\r
258\r
259 //Pt\r
0ed078e1 260 fHistPt = new TH1F("fHistPt","Pt (acceptance);p_{t} (GeV/c);Entries",1000,0.,10.);\r
261 fList->Add(fHistPt);\r
262\r
263 fHistPtPions = new TH1F("fHistPtPions","Pt (acceptance - pions);p_{t} (GeV/c);Entries",1000,0.,10.);\r
264 fList->Add(fHistPtPions);\r
265 fHistPtKaons = new TH1F("fHistPtKaons","Pt (acceptance - kaons);p_{t} (GeV/c);Entries",1000,0.,10.);\r
266 fList->Add(fHistPtKaons);\r
267 fHistPtProtons = new TH1F("fHistPtProtons","Pt (acceptance - protons);p_{t} (GeV/c);Entries",1000,0.,10.);\r
268 fList->Add(fHistPtProtons);\r
269\r
270 //==============Balance function histograms================//\r
63d15f4b 271 // Initialize histograms if not done yet\r
272 if(!fBalance->GetHistNp(0)){\r
273 AliWarning("Histograms not yet initialized! --> Will be done now");\r
274 AliWarning("--> Add 'gBalance->InitHistograms()' in your configBalanceFunction");\r
275 fBalance->InitHistograms();\r
276 }\r
277\r
278 if(fRunShuffling) {\r
279 if(!fShuffledBalance->GetHistNp(0)) {\r
280 AliWarning("Histograms (shuffling) not yet initialized! --> Will be done now");\r
281 AliWarning("--> Add 'gBalance->InitHistograms()' in your configBalanceFunction");\r
282 fShuffledBalance->InitHistograms();\r
283 }\r
284 }\r
285\r
286 for(Int_t a = 0; a < ANALYSIS_TYPES; a++){\r
287 fListBF->Add(fBalance->GetHistNp(a));\r
288 fListBF->Add(fBalance->GetHistNn(a));\r
289 fListBF->Add(fBalance->GetHistNpn(a));\r
290 fListBF->Add(fBalance->GetHistNnn(a));\r
291 fListBF->Add(fBalance->GetHistNpp(a));\r
292 fListBF->Add(fBalance->GetHistNnp(a));\r
293\r
294 if(fRunShuffling) {\r
295 fListBFS->Add(fShuffledBalance->GetHistNp(a));\r
296 fListBFS->Add(fShuffledBalance->GetHistNn(a));\r
297 fListBFS->Add(fShuffledBalance->GetHistNpn(a));\r
298 fListBFS->Add(fShuffledBalance->GetHistNnn(a));\r
299 fListBFS->Add(fShuffledBalance->GetHistNpp(a));\r
300 fListBFS->Add(fShuffledBalance->GetHistNnp(a));\r
301 } \r
302 }\r
303\r
304 // Post output data.\r
0ed078e1 305 //PostData(1, fList);\r
306 //PostData(2, fListBF);\r
307 //if(fRunShuffling) PostData(3, fListBFS);\r
211b716d 308\r
309 TH1::AddDirectory(oldStatus);\r
310\r
63d15f4b 311}\r
312\r
313//________________________________________________________________________\r
0ed078e1 314void AliAnalysisTaskToyModel::Run(Int_t nEvents) {\r
63d15f4b 315 // Main loop\r
316 // Called for each event\r
9fd4b54e 317 Double_t vCharge = 0;\r
318 Double_t vY = 0.0;\r
319 Double_t vEta = 0.0;\r
320 Double_t vPhi = 0.0;\r
321 Double_t vP[3] = {0.,0.,0.};\r
322 Double_t vPt = 0.0;\r
323 Double_t vE = 0.0;\r
0ed078e1 324 Bool_t isPion = kFALSE, isKaon = kFALSE, isProton = kFALSE;\r
63d15f4b 325\r
0ed078e1 326 if(fUseAllCharges) {\r
327 fPtSpectraAllCharges->SetParameter(0,fParticleMass);\r
328 fPtSpectraAllCharges->SetParameter(1,fTemperatureAllCharges);\r
b4665e0f 329\r
0ed078e1 330 fAzimuthalAngleAllCharges->SetParameter(1,fDirectedFlowAllCharges);\r
331 fAzimuthalAngleAllCharges->SetParameter(2,fEllipticFlowAllCharges);\r
332 fAzimuthalAngleAllCharges->SetParameter(3,fTriangularFlowAllCharges);\r
333 fAzimuthalAngleAllCharges->SetParameter(4,fQuandrangularFlowAllCharges);\r
334 fAzimuthalAngleAllCharges->SetParameter(5,fPentangularFlowAllCharges);\r
335 }\r
b4665e0f 336 else {\r
0ed078e1 337 fPtSpectraPions->SetParameter(0,fPionMass);\r
338 fPtSpectraPions->SetParameter(1,fTemperaturePions);\r
339 fPtSpectraKaons->SetParameter(0,fKaonMass);\r
340 fPtSpectraKaons->SetParameter(1,fTemperatureKaons);\r
341 fPtSpectraProtons->SetParameter(0,fProtonMass);\r
342 fPtSpectraProtons->SetParameter(1,fTemperatureProtons);\r
343\r
344 fAzimuthalAnglePions->SetParameter(1,fDirectedFlowPions);\r
345 fAzimuthalAnglePions->SetParameter(2,fEllipticFlowPions);\r
346 fAzimuthalAnglePions->SetParameter(3,fTriangularFlowPions);\r
347 fAzimuthalAnglePions->SetParameter(4,fQuandrangularFlowPions);\r
348 fAzimuthalAnglePions->SetParameter(5,fPentangularFlowPions);\r
349\r
350 fAzimuthalAngleKaons->SetParameter(1,fDirectedFlowKaons);\r
351 fAzimuthalAngleKaons->SetParameter(2,fEllipticFlowKaons);\r
352 fAzimuthalAngleKaons->SetParameter(3,fTriangularFlowKaons);\r
353 fAzimuthalAngleKaons->SetParameter(4,fQuandrangularFlowKaons);\r
354 fAzimuthalAngleKaons->SetParameter(5,fPentangularFlowKaons);\r
355\r
356 fAzimuthalAngleProtons->SetParameter(1,fDirectedFlowProtons);\r
357 fAzimuthalAngleProtons->SetParameter(2,fEllipticFlowProtons);\r
358 fAzimuthalAngleProtons->SetParameter(3,fTriangularFlowProtons);\r
359 fAzimuthalAngleProtons->SetParameter(4,fQuandrangularFlowProtons);\r
360 fAzimuthalAngleProtons->SetParameter(5,fPentangularFlowProtons);\r
b4665e0f 361 }\r
63d15f4b 362\r
0ed078e1 363 for(Int_t iEvent = 0; iEvent < nEvents; iEvent++) {\r
364 // vector holding the charges/kinematics of all tracks (charge,y,eta,phi,p0,p1,p2,pt,E)\r
365 vector<Double_t> *chargeVectorShuffle[9]; // this will be shuffled\r
366 vector<Double_t> *chargeVector[9]; // original charge\r
367 for(Int_t i = 0; i < 9; i++){\r
368 chargeVectorShuffle[i] = new vector<Double_t>;\r
369 chargeVector[i] = new vector<Double_t>;\r
b4665e0f 370 }\r
0ed078e1 371\r
372 fHistEventStats->Fill(1);\r
373 fHistEventStats->Fill(2);\r
374 fHistEventStats->Fill(3);\r
375\r
376 if((iEvent%1000) == 0) \r
377 cout<<"Event: "<<iEvent<<"/"<<nEvents<<endl;\r
378 \r
379 //Multiplicities\r
380 Int_t nMultiplicity = (Int_t)(gRandom->Gaus(fTotalMultiplicityMean,fTotalMultiplicitySigma));\r
381 Int_t nNetCharge = (Int_t)(gRandom->Gaus(fNetChargeMean,fNetChargeSigma));\r
6b3bc65c 382 \r
383 Int_t nGeneratedPositive = (Int_t)((nMultiplicity/2) + nNetCharge);\r
384 Int_t nGeneratedNegative = nMultiplicity - nGeneratedPositive;\r
385 if(fUseDebug) \r
386 Printf("Total multiplicity: %d - Generated positive: %d - Generated negative: %d",nMultiplicity,nGeneratedPositive,nGeneratedNegative);\r
387\r
388 //Int_t nGeneratedPositive = 0, nGeneratedNegative = 0;\r
0ed078e1 389 Int_t nGeneratedPions = 0, nGeneratedKaons = 0, nGeneratedProtons = 0;\r
390 \r
391 //Randomization of the reaction plane\r
392 //fReactionPlane = 2.0*TMath::Pi()*gRandom->Rndm();\r
393 fReactionPlane = 0.0;\r
394 if(fUseAllCharges) \r
395 fAzimuthalAngleAllCharges->SetParameter(0,fReactionPlane);\r
b4665e0f 396 else {\r
0ed078e1 397 fAzimuthalAnglePions->SetParameter(0,fReactionPlane);\r
398 fAzimuthalAngleKaons->SetParameter(0,fReactionPlane);\r
399 fAzimuthalAngleProtons->SetParameter(0,fReactionPlane);\r
63d15f4b 400 }\r
0ed078e1 401 \r
402 Int_t gNumberOfAcceptedParticles = 0;\r
6b3bc65c 403 Int_t gNumberOfAcceptedPositiveParticles = 0;\r
404 Int_t gNumberOfAcceptedNegativeParticles = 0;\r
405 \r
406 //Generate positive particles\r
407 for(Int_t iParticleCount = 0; iParticleCount < nGeneratedPositive; iParticleCount++) {\r
0ed078e1 408 isPion = kFALSE; isKaon = kFALSE; isProton = kFALSE;\r
6b3bc65c 409 if(fUseDebug) \r
410 Printf("Generating positive: %d(%d)",iParticleCount+1,nGeneratedPositive);\r
0ed078e1 411\r
6b3bc65c 412 //Pseudo-rapidity sampled from a Gaussian centered @ 0\r
9fd4b54e 413 vEta = gRandom->Gaus(0.0,4.0);\r
6b3bc65c 414\r
415 //Fill QA histograms (full phase space)\r
9fd4b54e 416 fHistEtaTotal->Fill(vEta);\r
6b3bc65c 417\r
9fd4b54e 418 vCharge = 1.0;\r
6b3bc65c 419 //nGeneratedPositive += 1;\r
0ed078e1 420 \r
6b3bc65c 421 //Acceptance\r
9fd4b54e 422 if((vEta < fEtaMin) || (vEta > fEtaMax)) continue;\r
6b3bc65c 423\r
0ed078e1 424 if(!fUseAllCharges) {\r
425 //Decide the specie\r
426 Double_t randomNumberSpecies = gRandom->Rndm();\r
427 if((randomNumberSpecies >= 0.0)&&(randomNumberSpecies < fPionPercentage)) {\r
428 nGeneratedPions += 1;\r
9fd4b54e 429 vPt = fPtSpectraPions->GetRandom();\r
430 vPhi = fAzimuthalAnglePions->GetRandom();\r
0ed078e1 431 fParticleMass = fPionMass;\r
432 isPion = kTRUE;\r
433 }\r
434 else if((randomNumberSpecies >= fPionPercentage)&&(randomNumberSpecies < fPionPercentage + fKaonPercentage)) {\r
435 nGeneratedKaons += 1;\r
9fd4b54e 436 vPt = fPtSpectraKaons->GetRandom();\r
437 vPhi = fAzimuthalAngleKaons->GetRandom();\r
0ed078e1 438 fParticleMass = fKaonMass;\r
439 isKaon = kTRUE;\r
440 }\r
441 else if((randomNumberSpecies >= fPionPercentage + fKaonPercentage)&&(randomNumberSpecies < fPionPercentage + fKaonPercentage + fProtonPercentage)) {\r
442 nGeneratedProtons += 1;\r
9fd4b54e 443 vPt = fPtSpectraProtons->GetRandom();\r
444 vPhi = fAzimuthalAngleProtons->GetRandom();\r
0ed078e1 445 fParticleMass = fProtonMass;\r
446 isProton = kTRUE;\r
447 }\r
b4665e0f 448 }\r
0ed078e1 449 else {\r
9fd4b54e 450 vPt = fPtSpectraAllCharges->GetRandom();\r
451 vPhi = fAzimuthalAngleAllCharges->GetRandom();\r
0ed078e1 452 }\r
0ed078e1 453 \r
9fd4b54e 454 vP[0] = vPt*TMath::Cos(vPhi);\r
455 vP[1] = vPt*TMath::Sin(vPhi);\r
456 vP[2] = vPt*TMath::SinH(vEta);\r
457 vE = TMath::Sqrt(TMath::Power(fParticleMass,2) +\r
458 TMath::Power(vP[0],2) +\r
459 TMath::Power(vP[1],2) +\r
460 TMath::Power(vP[2],2));\r
0ed078e1 461 \r
9fd4b54e 462 vY = 0.5*TMath::Log((vE + vP[2])/(vE - vP[2]));\r
0ed078e1 463 \r
6b3bc65c 464 //pt coverage\r
9fd4b54e 465 if((vPt < fPtMin) || (vPt > fPtMax)) continue;\r
466 //Printf("pt: %lf - mins: %lf - max: %lf",vPt,fPtMin,fPtMax);\r
6b3bc65c 467\r
468 //acceptance filter\r
469 if(fUseAcceptanceParameterization) {\r
470 Double_t gRandomNumberForAcceptance = gRandom->Rndm();\r
9fd4b54e 471 if(gRandomNumberForAcceptance > fAcceptanceParameterization->Eval(vPt)) \r
6b3bc65c 472 continue;\r
473 }\r
474\r
475 gNumberOfAcceptedPositiveParticles += 1;\r
476\r
477 //Fill QA histograms (acceptance);\r
9fd4b54e 478 fHistEta->Fill(vEta);\r
479 fHistRapidity->Fill(vY);\r
480 fHistPhi->Fill(vPhi);\r
481 fHistPt->Fill(vPt);\r
6b3bc65c 482 if(isPion) {\r
9fd4b54e 483 fHistRapidityPions->Fill(vY);\r
484 fHistPhiPions->Fill(vPhi);\r
485 fHistPtPions->Fill(vPt);\r
6b3bc65c 486 }\r
487 else if(isKaon) {\r
9fd4b54e 488 fHistRapidityKaons->Fill(vY);\r
489 fHistPhiKaons->Fill(vPhi);\r
490 fHistPtKaons->Fill(vPt);\r
6b3bc65c 491 }\r
492 else if(isProton) {\r
9fd4b54e 493 fHistRapidityProtons->Fill(vY);\r
494 fHistPhiProtons->Fill(vPhi);\r
495 fHistPtProtons->Fill(vPt);\r
6b3bc65c 496 }\r
497\r
498 // fill charge vector\r
9fd4b54e 499 chargeVector[0]->push_back(vCharge);\r
500 chargeVector[1]->push_back(vY);\r
501 chargeVector[2]->push_back(vEta);\r
502 chargeVector[3]->push_back(TMath::RadToDeg()*vPhi);\r
503 chargeVector[4]->push_back(vP[0]);\r
504 chargeVector[5]->push_back(vP[1]);\r
505 chargeVector[6]->push_back(vP[2]);\r
506 chargeVector[7]->push_back(vPt);\r
507 chargeVector[8]->push_back(vE);\r
6b3bc65c 508 \r
509 if(fRunShuffling) {\r
9fd4b54e 510 chargeVectorShuffle[0]->push_back(vCharge);\r
511 chargeVectorShuffle[1]->push_back(vY);\r
512 chargeVectorShuffle[2]->push_back(vEta);\r
513 chargeVectorShuffle[3]->push_back(TMath::RadToDeg()*vPhi);\r
514 chargeVectorShuffle[4]->push_back(vP[0]);\r
515 chargeVectorShuffle[5]->push_back(vP[1]);\r
516 chargeVectorShuffle[6]->push_back(vP[2]);\r
517 chargeVectorShuffle[7]->push_back(vPt);\r
518 chargeVectorShuffle[8]->push_back(vE);\r
6b3bc65c 519 }\r
520 gNumberOfAcceptedParticles += 1;\r
521 }//generated positive particle loop\r
522 \r
523 //Generate negative particles\r
524 for(Int_t iParticleCount = 0; iParticleCount < nGeneratedNegative; iParticleCount++) {\r
525 isPion = kFALSE; isKaon = kFALSE; isProton = kFALSE;\r
526 if(fUseDebug) \r
527 Printf("Generating negative: %d(%d)",iParticleCount+1,nGeneratedNegative);\r
528\r
529 //Pseudo-rapidity sampled from a Gaussian centered @ 0\r
9fd4b54e 530 vEta = gRandom->Gaus(0.0,4.0);\r
6b3bc65c 531\r
0ed078e1 532 //Fill QA histograms (full phase space)\r
9fd4b54e 533 fHistEtaTotal->Fill(vEta);\r
0ed078e1 534\r
9fd4b54e 535 vCharge = -1.0;\r
6b3bc65c 536 //nGeneratedNegative += 1;\r
537 \r
0ed078e1 538 //Acceptance\r
9fd4b54e 539 if((vEta < fEtaMin) || (vEta > fEtaMax)) continue;\r
6b3bc65c 540\r
541 if(!fUseAllCharges) {\r
542 //Decide the specie\r
543 Double_t randomNumberSpecies = gRandom->Rndm();\r
544 if((randomNumberSpecies >= 0.0)&&(randomNumberSpecies < fPionPercentage)) {\r
545 nGeneratedPions += 1;\r
9fd4b54e 546 vPt = fPtSpectraPions->GetRandom();\r
547 vPhi = fAzimuthalAnglePions->GetRandom();\r
6b3bc65c 548 fParticleMass = fPionMass;\r
549 isPion = kTRUE;\r
550 }\r
551 else if((randomNumberSpecies >= fPionPercentage)&&(randomNumberSpecies < fPionPercentage + fKaonPercentage)) {\r
552 nGeneratedKaons += 1;\r
9fd4b54e 553 vPt = fPtSpectraKaons->GetRandom();\r
554 vPhi = fAzimuthalAngleKaons->GetRandom();\r
6b3bc65c 555 fParticleMass = fKaonMass;\r
556 isKaon = kTRUE;\r
557 }\r
558 else if((randomNumberSpecies >= fPionPercentage + fKaonPercentage)&&(randomNumberSpecies < fPionPercentage + fKaonPercentage + fProtonPercentage)) {\r
559 nGeneratedProtons += 1;\r
9fd4b54e 560 vPt = fPtSpectraProtons->GetRandom();\r
561 vPhi = fAzimuthalAngleProtons->GetRandom();\r
6b3bc65c 562 fParticleMass = fProtonMass;\r
563 isProton = kTRUE;\r
564 }\r
565 }\r
566 else {\r
9fd4b54e 567 vPt = fPtSpectraAllCharges->GetRandom();\r
568 vPhi = fAzimuthalAngleAllCharges->GetRandom();\r
6b3bc65c 569 }\r
570 \r
9fd4b54e 571 vP[0] = vPt*TMath::Cos(vPhi);\r
572 vP[1] = vPt*TMath::Sin(vPhi);\r
573 vP[2] = vPt*TMath::SinH(vEta);\r
574 vE = TMath::Sqrt(TMath::Power(fParticleMass,2) +\r
575 TMath::Power(vP[0],2) +\r
576 TMath::Power(vP[1],2) +\r
577 TMath::Power(vP[2],2));\r
6b3bc65c 578 \r
9fd4b54e 579 vY = 0.5*TMath::Log((vE + vP[2])/(vE - vP[2]));\r
6b3bc65c 580 \r
581 //pt coverage\r
9fd4b54e 582 if((vPt < fPtMin) || (vPt > fPtMax)) continue;\r
583 //Printf("pt: %lf - mins: %lf - max: %lf",vPt,fPtMin,fPtMax);\r
0ed078e1 584\r
6b3bc65c 585 //acceptance filter\r
586 if(fUseAcceptanceParameterization) {\r
587 Double_t gRandomNumberForAcceptance = gRandom->Rndm();\r
9fd4b54e 588 if(gRandomNumberForAcceptance > fAcceptanceParameterization->Eval(vPt)) \r
6b3bc65c 589 continue;\r
590 }\r
591\r
592 gNumberOfAcceptedNegativeParticles += 1;\r
593\r
0ed078e1 594 //Fill QA histograms (acceptance);\r
9fd4b54e 595 fHistEta->Fill(vEta);\r
596 fHistRapidity->Fill(vY);\r
597 fHistPhi->Fill(vPhi);\r
598 fHistPt->Fill(vPt);\r
0ed078e1 599 if(isPion) {\r
9fd4b54e 600 fHistRapidityPions->Fill(vY);\r
601 fHistPhiPions->Fill(vPhi);\r
602 fHistPtPions->Fill(vPt);\r
0ed078e1 603 }\r
604 else if(isKaon) {\r
9fd4b54e 605 fHistRapidityKaons->Fill(vY);\r
606 fHistPhiKaons->Fill(vPhi);\r
607 fHistPtKaons->Fill(vPt);\r
0ed078e1 608 }\r
609 else if(isProton) {\r
9fd4b54e 610 fHistRapidityProtons->Fill(vY);\r
611 fHistPhiProtons->Fill(vPhi);\r
612 fHistPtProtons->Fill(vPt);\r
0ed078e1 613 }\r
614\r
615 // fill charge vector\r
9fd4b54e 616 chargeVector[0]->push_back(vCharge);\r
617 chargeVector[1]->push_back(vY);\r
618 chargeVector[2]->push_back(vEta);\r
619 chargeVector[3]->push_back(TMath::RadToDeg()*vPhi);\r
620 chargeVector[4]->push_back(vP[0]);\r
621 chargeVector[5]->push_back(vP[1]);\r
622 chargeVector[6]->push_back(vP[2]);\r
623 chargeVector[7]->push_back(vPt);\r
624 chargeVector[8]->push_back(vE);\r
0ed078e1 625 \r
626 if(fRunShuffling) {\r
9fd4b54e 627 chargeVectorShuffle[0]->push_back(vCharge);\r
628 chargeVectorShuffle[1]->push_back(vY);\r
629 chargeVectorShuffle[2]->push_back(vEta);\r
630 chargeVectorShuffle[3]->push_back(TMath::RadToDeg()*vPhi);\r
631 chargeVectorShuffle[4]->push_back(vP[0]);\r
632 chargeVectorShuffle[5]->push_back(vP[1]);\r
633 chargeVectorShuffle[6]->push_back(vP[2]);\r
634 chargeVectorShuffle[7]->push_back(vPt);\r
635 chargeVectorShuffle[8]->push_back(vE);\r
0ed078e1 636 }\r
637 gNumberOfAcceptedParticles += 1;\r
6b3bc65c 638 }//generated negative particle loop\r
78a6252b 639 \r
6b3bc65c 640 //Dynamical correlations\r
9fd4b54e 641 Double_t vChargePrime = 0;\r
642 Double_t vYPrime = 0.0;\r
643 Double_t vEtaPrime = 0.0;\r
644 Double_t vPhiPrime = 0.0;\r
645 Double_t vPPrime[3] = {0.,0.,0.};\r
646 Double_t vPtPrime = 0.0;\r
647 Double_t vEPrime = 0.0;\r
6b3bc65c 648 Int_t nGeneratedPositiveDynamicalCorrelations = 0;\r
649 Int_t nGeneratedNegativeDynamicalCorrelations = 0;\r
650 //Generate "correlated" particles \r
651 if(fUseDynamicalCorrelations) {\r
652 Int_t gNumberOfDynamicalCorrelations = (Int_t)(0.5*gNumberOfAcceptedParticles*fDynamicalCorrelationsPercentage);\r
653 for(Int_t iDynamicalCorrelations = 0; iDynamicalCorrelations < gNumberOfDynamicalCorrelations; iDynamicalCorrelations++) {\r
654 isPion = kFALSE; isKaon = kFALSE; isProton = kFALSE;\r
655 \r
656 //Pseudo-rapidity sampled from a Gaussian centered @ 0\r
9fd4b54e 657 vEta = gRandom->Gaus(0.0,0.1);\r
658 vCharge = 1.0;\r
6b3bc65c 659 nGeneratedPositiveDynamicalCorrelations += 1;\r
660 \r
9fd4b54e 661 vEtaPrime = -vEta;\r
662 vChargePrime = -1.0;\r
6b3bc65c 663 nGeneratedNegativeDynamicalCorrelations += 1;\r
664 \r
665 //Acceptance\r
9fd4b54e 666 if((vEta < fEtaMin) || (vEta > fEtaMax)) continue;\r
667 if((vEtaPrime < fEtaMin) || (vEtaPrime > fEtaMax)) continue;\r
0ed078e1 668 \r
6b3bc65c 669 if(!fUseAllCharges) {\r
670 //Decide the specie\r
671 Double_t randomNumberSpecies = gRandom->Rndm();\r
672 if((randomNumberSpecies >= 0.0)&&(randomNumberSpecies < fPionPercentage)) {\r
673 nGeneratedPions += 1;\r
9fd4b54e 674 vPt = fPtSpectraPions->GetRandom();\r
675 vPhi = fAzimuthalAnglePions->GetRandom();\r
6b3bc65c 676 fParticleMass = fPionMass;\r
677 isPion = kTRUE;\r
678 }\r
679 else if((randomNumberSpecies >= fPionPercentage)&&(randomNumberSpecies < fPionPercentage + fKaonPercentage)) {\r
680 nGeneratedKaons += 1;\r
9fd4b54e 681 vPt = fPtSpectraKaons->GetRandom();\r
682 vPhi = fAzimuthalAngleKaons->GetRandom();\r
6b3bc65c 683 fParticleMass = fKaonMass;\r
684 isKaon = kTRUE;\r
685 }\r
686 else if((randomNumberSpecies >= fPionPercentage + fKaonPercentage)&&(randomNumberSpecies < fPionPercentage + fKaonPercentage + fProtonPercentage)) {\r
687 nGeneratedProtons += 1;\r
9fd4b54e 688 vPt = fPtSpectraProtons->GetRandom();\r
689 vPtPrime = vPt;\r
690 vPhi = fAzimuthalAngleProtons->GetRandom();\r
6b3bc65c 691 fParticleMass = fProtonMass;\r
692 isProton = kTRUE;\r
693 }\r
694 }\r
695 else {\r
9fd4b54e 696 vPt = fPtSpectraAllCharges->GetRandom();\r
697 vPhi = fAzimuthalAngleAllCharges->GetRandom();\r
6b3bc65c 698 }\r
9fd4b54e 699 vPtPrime = vPt;\r
700 vPhiPrime = vPhi;\r
701\r
702 vP[0] = vPt*TMath::Cos(vPhi);\r
703 vP[1] = vPt*TMath::Sin(vPhi);\r
704 vP[2] = vPt*TMath::SinH(vEta);\r
705 vE = TMath::Sqrt(TMath::Power(fParticleMass,2) +\r
706 TMath::Power(vP[0],2) +\r
707 TMath::Power(vP[1],2) +\r
708 TMath::Power(vP[2],2));\r
6b3bc65c 709 \r
9fd4b54e 710 vY = 0.5*TMath::Log((vE + vP[2])/(vE - vP[2]));\r
711\r
712 vPPrime[0] = vPtPrime*TMath::Cos(vPhiPrime);\r
713 vPPrime[1] = vPtPrime*TMath::Sin(vPhiPrime);\r
714 vPPrime[2] = vPtPrime*TMath::SinH(vEtaPrime);\r
715 vEPrime = TMath::Sqrt(TMath::Power(fParticleMass,2) +\r
716 TMath::Power(vPPrime[0],2) +\r
717 TMath::Power(vPPrime[1],2) +\r
718 TMath::Power(vPPrime[2],2));\r
6b3bc65c 719 \r
9fd4b54e 720 vYPrime = 0.5*TMath::Log((vEPrime + vPPrime[2])/(vEPrime - vPPrime[2]));\r
6b3bc65c 721 \r
722 //pt coverage\r
9fd4b54e 723 if((vPt < fPtMin) || (vPt > fPtMax)) continue;\r
724 if((vPtPrime < fPtMin) || (vPtPrime > fPtMax)) continue;\r
6b3bc65c 725\r
726 //acceptance filter\r
727 if(fUseAcceptanceParameterization) {\r
728 Double_t gRandomNumberForAcceptance = gRandom->Rndm();\r
9fd4b54e 729 if(gRandomNumberForAcceptance > fAcceptanceParameterization->Eval(vPt)) \r
6b3bc65c 730 continue;\r
731 \r
732 Double_t gRandomNumberForAcceptancePrime = gRandom->Rndm();\r
9fd4b54e 733 if(gRandomNumberForAcceptancePrime > fAcceptanceParameterization->Eval(vPtPrime)) \r
6b3bc65c 734 continue;\r
735 }\r
736 \r
737 // fill charge vector (positive)\r
9fd4b54e 738 chargeVector[0]->push_back(vCharge);\r
739 chargeVector[1]->push_back(vY);\r
740 chargeVector[2]->push_back(vEta);\r
741 chargeVector[3]->push_back(TMath::RadToDeg()*vPhi);\r
742 chargeVector[4]->push_back(vP[0]);\r
743 chargeVector[5]->push_back(vP[1]);\r
744 chargeVector[6]->push_back(vP[2]);\r
745 chargeVector[7]->push_back(vPt);\r
746 chargeVector[8]->push_back(vE);\r
6b3bc65c 747 \r
748 if(fRunShuffling) {\r
9fd4b54e 749 chargeVectorShuffle[0]->push_back(vCharge);\r
750 chargeVectorShuffle[1]->push_back(vY);\r
751 chargeVectorShuffle[2]->push_back(vEta);\r
752 chargeVectorShuffle[3]->push_back(TMath::RadToDeg()*vPhi);\r
753 chargeVectorShuffle[4]->push_back(vP[0]);\r
754 chargeVectorShuffle[5]->push_back(vP[1]);\r
755 chargeVectorShuffle[6]->push_back(vP[2]);\r
756 chargeVectorShuffle[7]->push_back(vPt);\r
757 chargeVectorShuffle[8]->push_back(vE);\r
6b3bc65c 758 }\r
759\r
760 // fill charge vector (negative)\r
9fd4b54e 761 chargeVector[0]->push_back(vChargePrime);\r
762 chargeVector[1]->push_back(vYPrime);\r
763 chargeVector[2]->push_back(vEtaPrime);\r
764 chargeVector[3]->push_back(TMath::RadToDeg()*vPhiPrime);\r
765 chargeVector[4]->push_back(vPPrime[0]);\r
766 chargeVector[5]->push_back(vPPrime[1]);\r
767 chargeVector[6]->push_back(vPPrime[2]);\r
768 chargeVector[7]->push_back(vPtPrime);\r
769 chargeVector[8]->push_back(vEPrime);\r
6b3bc65c 770 \r
771 if(fRunShuffling) {\r
9fd4b54e 772 chargeVectorShuffle[0]->push_back(vChargePrime);\r
773 chargeVectorShuffle[1]->push_back(vYPrime);\r
774 chargeVectorShuffle[2]->push_back(vEtaPrime);\r
775 chargeVectorShuffle[3]->push_back(TMath::RadToDeg()*vPhiPrime);\r
776 chargeVectorShuffle[4]->push_back(vPPrime[0]);\r
777 chargeVectorShuffle[5]->push_back(vPPrime[1]);\r
778 chargeVectorShuffle[6]->push_back(vPPrime[2]);\r
779 chargeVectorShuffle[7]->push_back(vPtPrime);\r
780 chargeVectorShuffle[8]->push_back(vEPrime);\r
6b3bc65c 781 }\r
782\r
783 gNumberOfAcceptedParticles += 2;\r
784 }//loop over the dynamical correlations\r
785 }// usage of dynamical correlations\r
786\r
0ed078e1 787 if(fUseDebug) {\r
788 Printf("=======================================================");\r
789 Printf("Total: %d - Total positive: %d - Total negative: %d",nMultiplicity,nGeneratedPositive,nGeneratedNegative);\r
6b3bc65c 790 Printf("Accepted positive: %d - Accepted negative: %d",gNumberOfAcceptedPositiveParticles,gNumberOfAcceptedNegativeParticles);\r
791 Printf("Correlations: %d - Correlations positive: %d - Correlations negative: %d",nGeneratedPositiveDynamicalCorrelations+nGeneratedNegativeDynamicalCorrelations,nGeneratedPositiveDynamicalCorrelations,nGeneratedNegativeDynamicalCorrelations);\r
792 Printf("Number of accepted particles: %d",gNumberOfAcceptedParticles);\r
0ed078e1 793 if(!fUseAllCharges)\r
794 Printf("Pions: %lf - Kaons: %lf - Protons: %lf",1.*nGeneratedPions/nMultiplicity,1.*nGeneratedKaons/nMultiplicity,1.*nGeneratedProtons/nMultiplicity);\r
6b3bc65c 795 //Printf("Calculating the balance function for %d particles",chargeVector[0]->size());\r
63d15f4b 796 }\r
63d15f4b 797\r
0ed078e1 798 fHistEventStats->Fill(4);\r
799 fHistNumberOfAcceptedParticles->Fill(gNumberOfAcceptedParticles);\r
800 fHistReactionPlane->Fill(fReactionPlane);\r
63d15f4b 801\r
0ed078e1 802 //Calculate the balance function\r
803 fBalance->CalculateBalance(gNumberOfAcceptedParticles,chargeVector);\r
804 \r
805 if(fRunShuffling) {\r
806 // shuffle charges\r
807 random_shuffle( chargeVectorShuffle[0]->begin(), chargeVectorShuffle[0]->end() );\r
808 fShuffledBalance->CalculateBalance(gNumberOfAcceptedParticles,chargeVectorShuffle);\r
809 }\r
63d15f4b 810 }\r
811} \r
812\r
813//________________________________________________________________________\r
0ed078e1 814void AliAnalysisTaskToyModel::FinishOutput() {\r
63d15f4b 815 //Printf("END BF");\r
816\r
0ed078e1 817 TFile *gOutput = TFile::Open("outputToyModel.root","recreate");\r
818 fList->Write();\r
819\r
63d15f4b 820 if (!fBalance) {\r
821 Printf("ERROR: fBalance not available");\r
822 return;\r
823 } \r
0ed078e1 824 fListBF->Write();\r
825\r
63d15f4b 826 if(fRunShuffling) {\r
827 if (!fShuffledBalance) {\r
828 Printf("ERROR: fShuffledBalance not available");\r
829 return;\r
830 }\r
0ed078e1 831 fListBFS->Write();\r
63d15f4b 832 }\r
0ed078e1 833 gOutput->Close();\r
63d15f4b 834}\r
835\r