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