Compatibility with the Root trunk
[u/mrichter/AliRoot.git] / PWGCF / EBYE / BalanceFunctions / AliAnalysisTaskToyModel.cxx
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
12 #include "TFile.h"\r
13 \r
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
24 using std::cout;\r
25 using std::endl;\r
26 \r
27 ClassImp(AliAnalysisTaskToyModel)\r
28 \r
29 //________________________________________________________________________\r
30 AliAnalysisTaskToyModel::AliAnalysisTaskToyModel() \r
31 : TObject(),\r
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
73   // Constructor\r
74 }\r
75 \r
76 //________________________________________________________________________\r
77 AliAnalysisTaskToyModel::~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
90 void AliAnalysisTaskToyModel::Init() {\r
91   //Initialize objects\r
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
96   //==============Particles and spectra==============//\r
97   TParticle *pion = new TParticle();\r
98   pion->SetPdgCode(211);\r
99   fPionMass = pion->GetMass();\r
100 \r
101   TParticle *kaon = new TParticle();\r
102   kaon->SetPdgCode(321);\r
103   fKaonMass = kaon->GetMass();\r
104   \r
105   TParticle *proton = new TParticle();\r
106   proton->SetPdgCode(2212);\r
107   fProtonMass = proton->GetMass();\r
108 \r
109   if(fUseAllCharges) {\r
110     fParticleMass = fPionMass;\r
111     fPtSpectraAllCharges = new TF1("fPtSpectraAllCharges","x*TMath::Exp(-TMath::Power([0]*[0]+x*x,0.5)/[1])",0.,5.);\r
112     fPtSpectraAllCharges->SetParName(0,"Mass");\r
113     fPtSpectraAllCharges->SetParName(1,"Temperature");\r
114   }\r
115   else {\r
116     fPtSpectraPions = new TF1("fPtSpectraPions","x*TMath::Exp(-TMath::Power([0]*[0]+x*x,0.5)/[1])",0.,5.);\r
117     fPtSpectraPions->SetParName(0,"Mass");\r
118     fPtSpectraPions->SetParName(1,"Temperature");\r
119     \r
120     fPtSpectraKaons = new TF1("fPtSpectraKaons","x*TMath::Exp(-TMath::Power([0]*[0]+x*x,0.5)/[1])",0.,5.);\r
121     fPtSpectraKaons->SetParName(0,"Mass");\r
122     fPtSpectraKaons->SetParName(1,"Temperature");\r
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
126     fPtSpectraProtons->SetParName(1,"Temperature");\r
127   }\r
128   //==============Particles and spectra==============//\r
129 \r
130   //==============Flow values==============//\r
131   if(fUseAllCharges) {\r
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
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
143     fAzimuthalAngleAllCharges->SetParName(2,"Elliptic flow"); \r
144     fAzimuthalAngleAllCharges->SetParName(3,"Triangular flow");\r
145     fAzimuthalAngleAllCharges->SetParName(4,"Quandrangular flow");\r
146     fAzimuthalAngleAllCharges->SetParName(5,"Pentangular flow");\r
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
152     fAzimuthalAnglePions->SetParName(2,"Elliptic flow"); \r
153     fAzimuthalAnglePions->SetParName(3,"Triangular flow");\r
154     fAzimuthalAnglePions->SetParName(4,"Quandrangular flow");\r
155     fAzimuthalAnglePions->SetParName(5,"Pentangular flow");\r
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
160     fAzimuthalAngleKaons->SetParName(2,"Elliptic flow"); \r
161     fAzimuthalAngleKaons->SetParName(3,"Triangular flow");\r
162     fAzimuthalAngleKaons->SetParName(4,"Quandrangular flow");\r
163     fAzimuthalAngleKaons->SetParName(5,"Pentangular flow");\r
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
168     fAzimuthalAngleProtons->SetParName(2,"Elliptic flow"); \r
169     fAzimuthalAngleProtons->SetParName(3,"Triangular flow");\r
170     fAzimuthalAngleProtons->SetParName(4,"Quandrangular flow");\r
171     fAzimuthalAngleProtons->SetParName(5,"Pentangular flow");\r
172   }\r
173   //==============Flow values==============//\r
174 }\r
175 \r
176 //________________________________________________________________________\r
177 void AliAnalysisTaskToyModel::CreateOutputObjects() {\r
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
207   //==============QA================//\r
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
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
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
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
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
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
299   //PostData(1, fList);\r
300   //PostData(2, fListBF);\r
301   //if(fRunShuffling) PostData(3, fListBFS);\r
302 }\r
303 \r
304 //________________________________________________________________________\r
305 void AliAnalysisTaskToyModel::Run(Int_t nEvents) {\r
306   // Main loop\r
307   // Called for each event\r
308   Double_t v_charge = 0;\r
309   Double_t v_y = 0.0;\r
310   Double_t v_eta = 0.0;\r
311   Double_t v_phi = 0.0;\r
312   Double_t v_p[3] = {0.,0.,0.};\r
313   Double_t v_pt = 0.0;\r
314   Double_t v_E = 0.0;\r
315   Bool_t isPion = kFALSE, isKaon = kFALSE, isProton = kFALSE;\r
316 \r
317   if(fUseAllCharges) {\r
318     fPtSpectraAllCharges->SetParameter(0,fParticleMass);\r
319     fPtSpectraAllCharges->SetParameter(1,fTemperatureAllCharges);\r
320 \r
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
327   else {\r
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
352   }\r
353 \r
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
361     }\r
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
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
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
387     else {\r
388       fAzimuthalAnglePions->SetParameter(0,fReactionPlane);\r
389       fAzimuthalAngleKaons->SetParameter(0,fReactionPlane);\r
390       fAzimuthalAngleProtons->SetParameter(0,fReactionPlane);\r
391     }\r
392     \r
393     Int_t gNumberOfAcceptedParticles = 0;\r
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
399       isPion = kFALSE; isKaon = kFALSE; isProton = kFALSE;\r
400       if(fUseDebug) \r
401         Printf("Generating positive: %d(%d)",iParticleCount+1,nGeneratedPositive);\r
402 \r
403       //Pseudo-rapidity sampled from a Gaussian centered @ 0\r
404       v_eta = gRandom->Gaus(0.0,4.0);\r
405 \r
406       //Fill QA histograms (full phase space)\r
407       fHistEtaTotal->Fill(v_eta);\r
408 \r
409       v_charge = 1.0;\r
410       //nGeneratedPositive += 1;\r
411       \r
412       //Acceptance\r
413       if((v_eta < fEtaMin) || (v_eta > fEtaMax)) continue;\r
414 \r
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
420           v_pt = fPtSpectraPions->GetRandom();\r
421           v_phi = fAzimuthalAnglePions->GetRandom();\r
422           fParticleMass = fPionMass;\r
423           isPion = kTRUE;\r
424         }\r
425         else if((randomNumberSpecies >= fPionPercentage)&&(randomNumberSpecies < fPionPercentage + fKaonPercentage)) {\r
426           nGeneratedKaons += 1;\r
427           v_pt = fPtSpectraKaons->GetRandom();\r
428           v_phi = fAzimuthalAngleKaons->GetRandom();\r
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
434           v_pt = fPtSpectraProtons->GetRandom();\r
435           v_phi = fAzimuthalAngleProtons->GetRandom();\r
436           fParticleMass = fProtonMass;\r
437           isProton = kTRUE;\r
438         }\r
439       }\r
440       else {\r
441         v_pt = fPtSpectraAllCharges->GetRandom();\r
442         v_phi = fAzimuthalAngleAllCharges->GetRandom();\r
443       }\r
444       \r
445       v_p[0] = v_pt*TMath::Cos(v_phi);\r
446       v_p[1] = v_pt*TMath::Sin(v_phi);\r
447       v_p[2] = v_pt*TMath::SinH(v_eta);\r
448       v_E = TMath::Sqrt(TMath::Power(fParticleMass,2) +\r
449                         TMath::Power(v_p[0],2) +\r
450                         TMath::Power(v_p[1],2) +\r
451                         TMath::Power(v_p[2],2));\r
452       \r
453       v_y = 0.5*TMath::Log((v_E + v_p[2])/(v_E - v_p[2]));\r
454       \r
455       //pt coverage\r
456       if((v_pt < fPtMin) || (v_pt > fPtMax)) continue;\r
457       //Printf("pt: %lf - mins: %lf - max: %lf",v_pt,fPtMin,fPtMax);\r
458 \r
459       //acceptance filter\r
460       if(fUseAcceptanceParameterization) {\r
461         Double_t gRandomNumberForAcceptance = gRandom->Rndm();\r
462         if(gRandomNumberForAcceptance > fAcceptanceParameterization->Eval(v_pt)) \r
463           continue;\r
464       }\r
465 \r
466       gNumberOfAcceptedPositiveParticles += 1;\r
467 \r
468       //Fill QA histograms (acceptance);\r
469       fHistEta->Fill(v_eta);\r
470       fHistRapidity->Fill(v_y);\r
471       fHistPhi->Fill(v_phi);\r
472       fHistPt->Fill(v_pt);\r
473       if(isPion) {\r
474         fHistRapidityPions->Fill(v_y);\r
475         fHistPhiPions->Fill(v_phi);\r
476         fHistPtPions->Fill(v_pt);\r
477       }\r
478       else if(isKaon) {\r
479         fHistRapidityKaons->Fill(v_y);\r
480         fHistPhiKaons->Fill(v_phi);\r
481         fHistPtKaons->Fill(v_pt);\r
482       }\r
483       else if(isProton) {\r
484         fHistRapidityProtons->Fill(v_y);\r
485         fHistPhiProtons->Fill(v_phi);\r
486         fHistPtProtons->Fill(v_pt);\r
487       }\r
488 \r
489       // fill charge vector\r
490       chargeVector[0]->push_back(v_charge);\r
491       chargeVector[1]->push_back(v_y);\r
492       chargeVector[2]->push_back(v_eta);\r
493       chargeVector[3]->push_back(TMath::RadToDeg()*v_phi);\r
494       chargeVector[4]->push_back(v_p[0]);\r
495       chargeVector[5]->push_back(v_p[1]);\r
496       chargeVector[6]->push_back(v_p[2]);\r
497       chargeVector[7]->push_back(v_pt);\r
498       chargeVector[8]->push_back(v_E);\r
499       \r
500       if(fRunShuffling) {\r
501         chargeVectorShuffle[0]->push_back(v_charge);\r
502         chargeVectorShuffle[1]->push_back(v_y);\r
503         chargeVectorShuffle[2]->push_back(v_eta);\r
504         chargeVectorShuffle[3]->push_back(TMath::RadToDeg()*v_phi);\r
505         chargeVectorShuffle[4]->push_back(v_p[0]);\r
506         chargeVectorShuffle[5]->push_back(v_p[1]);\r
507         chargeVectorShuffle[6]->push_back(v_p[2]);\r
508         chargeVectorShuffle[7]->push_back(v_pt);\r
509         chargeVectorShuffle[8]->push_back(v_E);\r
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
521       v_eta = gRandom->Gaus(0.0,4.0);\r
522 \r
523       //Fill QA histograms (full phase space)\r
524       fHistEtaTotal->Fill(v_eta);\r
525 \r
526       v_charge = -1.0;\r
527       //nGeneratedNegative += 1;\r
528       \r
529       //Acceptance\r
530       if((v_eta < fEtaMin) || (v_eta > fEtaMax)) continue;\r
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
537           v_pt = fPtSpectraPions->GetRandom();\r
538           v_phi = fAzimuthalAnglePions->GetRandom();\r
539           fParticleMass = fPionMass;\r
540           isPion = kTRUE;\r
541         }\r
542         else if((randomNumberSpecies >= fPionPercentage)&&(randomNumberSpecies < fPionPercentage + fKaonPercentage)) {\r
543           nGeneratedKaons += 1;\r
544           v_pt = fPtSpectraKaons->GetRandom();\r
545           v_phi = fAzimuthalAngleKaons->GetRandom();\r
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
551           v_pt = fPtSpectraProtons->GetRandom();\r
552           v_phi = fAzimuthalAngleProtons->GetRandom();\r
553           fParticleMass = fProtonMass;\r
554           isProton = kTRUE;\r
555         }\r
556       }\r
557       else {\r
558         v_pt = fPtSpectraAllCharges->GetRandom();\r
559         v_phi = fAzimuthalAngleAllCharges->GetRandom();\r
560       }\r
561       \r
562       v_p[0] = v_pt*TMath::Cos(v_phi);\r
563       v_p[1] = v_pt*TMath::Sin(v_phi);\r
564       v_p[2] = v_pt*TMath::SinH(v_eta);\r
565       v_E = TMath::Sqrt(TMath::Power(fParticleMass,2) +\r
566                         TMath::Power(v_p[0],2) +\r
567                         TMath::Power(v_p[1],2) +\r
568                         TMath::Power(v_p[2],2));\r
569       \r
570       v_y = 0.5*TMath::Log((v_E + v_p[2])/(v_E - v_p[2]));\r
571       \r
572       //pt coverage\r
573       if((v_pt < fPtMin) || (v_pt > fPtMax)) continue;\r
574       //Printf("pt: %lf - mins: %lf - max: %lf",v_pt,fPtMin,fPtMax);\r
575 \r
576      //acceptance filter\r
577       if(fUseAcceptanceParameterization) {\r
578         Double_t gRandomNumberForAcceptance = gRandom->Rndm();\r
579         if(gRandomNumberForAcceptance > fAcceptanceParameterization->Eval(v_pt)) \r
580           continue;\r
581       }\r
582 \r
583       gNumberOfAcceptedNegativeParticles += 1;\r
584 \r
585       //Fill QA histograms (acceptance);\r
586       fHistEta->Fill(v_eta);\r
587       fHistRapidity->Fill(v_y);\r
588       fHistPhi->Fill(v_phi);\r
589       fHistPt->Fill(v_pt);\r
590       if(isPion) {\r
591         fHistRapidityPions->Fill(v_y);\r
592         fHistPhiPions->Fill(v_phi);\r
593         fHistPtPions->Fill(v_pt);\r
594       }\r
595       else if(isKaon) {\r
596         fHistRapidityKaons->Fill(v_y);\r
597         fHistPhiKaons->Fill(v_phi);\r
598         fHistPtKaons->Fill(v_pt);\r
599       }\r
600       else if(isProton) {\r
601         fHistRapidityProtons->Fill(v_y);\r
602         fHistPhiProtons->Fill(v_phi);\r
603         fHistPtProtons->Fill(v_pt);\r
604       }\r
605 \r
606       // fill charge vector\r
607       chargeVector[0]->push_back(v_charge);\r
608       chargeVector[1]->push_back(v_y);\r
609       chargeVector[2]->push_back(v_eta);\r
610       chargeVector[3]->push_back(TMath::RadToDeg()*v_phi);\r
611       chargeVector[4]->push_back(v_p[0]);\r
612       chargeVector[5]->push_back(v_p[1]);\r
613       chargeVector[6]->push_back(v_p[2]);\r
614       chargeVector[7]->push_back(v_pt);\r
615       chargeVector[8]->push_back(v_E);\r
616       \r
617       if(fRunShuffling) {\r
618         chargeVectorShuffle[0]->push_back(v_charge);\r
619         chargeVectorShuffle[1]->push_back(v_y);\r
620         chargeVectorShuffle[2]->push_back(v_eta);\r
621         chargeVectorShuffle[3]->push_back(TMath::RadToDeg()*v_phi);\r
622         chargeVectorShuffle[4]->push_back(v_p[0]);\r
623         chargeVectorShuffle[5]->push_back(v_p[1]);\r
624         chargeVectorShuffle[6]->push_back(v_p[2]);\r
625         chargeVectorShuffle[7]->push_back(v_pt);\r
626         chargeVectorShuffle[8]->push_back(v_E);\r
627       }\r
628       gNumberOfAcceptedParticles += 1;\r
629     }//generated negative particle loop\r
630     \r
631     //Dynamical correlations\r
632     Double_t v_chargePrime = 0;\r
633     Double_t v_yPrime = 0.0;\r
634     Double_t v_etaPrime = 0.0;\r
635     Double_t v_phiPrime = 0.0;\r
636     Double_t v_pPrime[3] = {0.,0.,0.};\r
637     Double_t v_ptPrime = 0.0;\r
638     Double_t v_EPrime = 0.0;\r
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
648         v_eta = gRandom->Gaus(0.0,0.1);\r
649         v_charge = 1.0;\r
650         nGeneratedPositiveDynamicalCorrelations += 1;\r
651         \r
652         v_etaPrime = -v_eta;\r
653         v_chargePrime = -1.0;\r
654         nGeneratedNegativeDynamicalCorrelations += 1;\r
655           \r
656         //Acceptance\r
657         if((v_eta < fEtaMin) || (v_eta > fEtaMax)) continue;\r
658         if((v_etaPrime < fEtaMin) || (v_etaPrime > fEtaMax)) continue;\r
659       \r
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
665             v_pt = fPtSpectraPions->GetRandom();\r
666             v_phi = fAzimuthalAnglePions->GetRandom();\r
667             fParticleMass = fPionMass;\r
668             isPion = kTRUE;\r
669           }\r
670           else if((randomNumberSpecies >= fPionPercentage)&&(randomNumberSpecies < fPionPercentage + fKaonPercentage)) {\r
671             nGeneratedKaons += 1;\r
672             v_pt = fPtSpectraKaons->GetRandom();\r
673             v_phi = fAzimuthalAngleKaons->GetRandom();\r
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
679             v_pt = fPtSpectraProtons->GetRandom();\r
680             v_ptPrime = v_pt;\r
681             v_phi = fAzimuthalAngleProtons->GetRandom();\r
682             fParticleMass = fProtonMass;\r
683             isProton = kTRUE;\r
684           }\r
685         }\r
686         else {\r
687           v_pt = fPtSpectraAllCharges->GetRandom();\r
688           v_phi = fAzimuthalAngleAllCharges->GetRandom();\r
689         }\r
690         v_ptPrime = v_pt;\r
691         v_phiPrime = v_phi;\r
692 \r
693         v_p[0] = v_pt*TMath::Cos(v_phi);\r
694         v_p[1] = v_pt*TMath::Sin(v_phi);\r
695         v_p[2] = v_pt*TMath::SinH(v_eta);\r
696         v_E = TMath::Sqrt(TMath::Power(fParticleMass,2) +\r
697                           TMath::Power(v_p[0],2) +\r
698                           TMath::Power(v_p[1],2) +\r
699                           TMath::Power(v_p[2],2));\r
700         \r
701         v_y = 0.5*TMath::Log((v_E + v_p[2])/(v_E - v_p[2]));\r
702 \r
703         v_pPrime[0] = v_ptPrime*TMath::Cos(v_phiPrime);\r
704         v_pPrime[1] = v_ptPrime*TMath::Sin(v_phiPrime);\r
705         v_pPrime[2] = v_ptPrime*TMath::SinH(v_etaPrime);\r
706         v_EPrime = TMath::Sqrt(TMath::Power(fParticleMass,2) +\r
707                           TMath::Power(v_pPrime[0],2) +\r
708                           TMath::Power(v_pPrime[1],2) +\r
709                           TMath::Power(v_pPrime[2],2));\r
710         \r
711         v_yPrime = 0.5*TMath::Log((v_EPrime + v_pPrime[2])/(v_EPrime - v_pPrime[2]));\r
712       \r
713         //pt coverage\r
714         if((v_pt < fPtMin) || (v_pt > fPtMax)) continue;\r
715         if((v_ptPrime < fPtMin) || (v_ptPrime > fPtMax)) continue;\r
716 \r
717         //acceptance filter\r
718         if(fUseAcceptanceParameterization) {\r
719           Double_t gRandomNumberForAcceptance = gRandom->Rndm();\r
720           if(gRandomNumberForAcceptance > fAcceptanceParameterization->Eval(v_pt)) \r
721             continue;\r
722           \r
723           Double_t gRandomNumberForAcceptancePrime = gRandom->Rndm();\r
724           if(gRandomNumberForAcceptancePrime > fAcceptanceParameterization->Eval(v_ptPrime)) \r
725             continue;\r
726         }\r
727       \r
728         // fill charge vector (positive)\r
729         chargeVector[0]->push_back(v_charge);\r
730         chargeVector[1]->push_back(v_y);\r
731         chargeVector[2]->push_back(v_eta);\r
732         chargeVector[3]->push_back(TMath::RadToDeg()*v_phi);\r
733         chargeVector[4]->push_back(v_p[0]);\r
734         chargeVector[5]->push_back(v_p[1]);\r
735         chargeVector[6]->push_back(v_p[2]);\r
736         chargeVector[7]->push_back(v_pt);\r
737         chargeVector[8]->push_back(v_E);\r
738         \r
739         if(fRunShuffling) {\r
740           chargeVectorShuffle[0]->push_back(v_charge);\r
741           chargeVectorShuffle[1]->push_back(v_y);\r
742           chargeVectorShuffle[2]->push_back(v_eta);\r
743           chargeVectorShuffle[3]->push_back(TMath::RadToDeg()*v_phi);\r
744           chargeVectorShuffle[4]->push_back(v_p[0]);\r
745           chargeVectorShuffle[5]->push_back(v_p[1]);\r
746           chargeVectorShuffle[6]->push_back(v_p[2]);\r
747           chargeVectorShuffle[7]->push_back(v_pt);\r
748           chargeVectorShuffle[8]->push_back(v_E);\r
749         }\r
750 \r
751         // fill charge vector (negative)\r
752         chargeVector[0]->push_back(v_chargePrime);\r
753         chargeVector[1]->push_back(v_yPrime);\r
754         chargeVector[2]->push_back(v_etaPrime);\r
755         chargeVector[3]->push_back(TMath::RadToDeg()*v_phiPrime);\r
756         chargeVector[4]->push_back(v_pPrime[0]);\r
757         chargeVector[5]->push_back(v_pPrime[1]);\r
758         chargeVector[6]->push_back(v_pPrime[2]);\r
759         chargeVector[7]->push_back(v_ptPrime);\r
760         chargeVector[8]->push_back(v_EPrime);\r
761         \r
762         if(fRunShuffling) {\r
763           chargeVectorShuffle[0]->push_back(v_chargePrime);\r
764           chargeVectorShuffle[1]->push_back(v_yPrime);\r
765           chargeVectorShuffle[2]->push_back(v_etaPrime);\r
766           chargeVectorShuffle[3]->push_back(TMath::RadToDeg()*v_phiPrime);\r
767           chargeVectorShuffle[4]->push_back(v_pPrime[0]);\r
768           chargeVectorShuffle[5]->push_back(v_pPrime[1]);\r
769           chargeVectorShuffle[6]->push_back(v_pPrime[2]);\r
770           chargeVectorShuffle[7]->push_back(v_ptPrime);\r
771           chargeVectorShuffle[8]->push_back(v_EPrime);\r
772         }\r
773 \r
774         gNumberOfAcceptedParticles += 2;\r
775       }//loop over the dynamical correlations\r
776     }// usage of dynamical correlations\r
777 \r
778     if(fUseDebug) {\r
779       Printf("=======================================================");\r
780       Printf("Total: %d - Total positive: %d - Total negative: %d",nMultiplicity,nGeneratedPositive,nGeneratedNegative);\r
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
784       if(!fUseAllCharges)\r
785         Printf("Pions: %lf - Kaons: %lf - Protons: %lf",1.*nGeneratedPions/nMultiplicity,1.*nGeneratedKaons/nMultiplicity,1.*nGeneratedProtons/nMultiplicity);\r
786       //Printf("Calculating the balance function for %d particles",chargeVector[0]->size());\r
787     }\r
788 \r
789     fHistEventStats->Fill(4);\r
790     fHistNumberOfAcceptedParticles->Fill(gNumberOfAcceptedParticles);\r
791     fHistReactionPlane->Fill(fReactionPlane);\r
792 \r
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
801   }\r
802 }      \r
803 \r
804 //________________________________________________________________________\r
805 void  AliAnalysisTaskToyModel::FinishOutput() {\r
806   //Printf("END BF");\r
807 \r
808   TFile *gOutput = TFile::Open("outputToyModel.root","recreate");\r
809   fList->Write();\r
810 \r
811   if (!fBalance) {\r
812     Printf("ERROR: fBalance not available");\r
813     return;\r
814   }  \r
815   fListBF->Write();\r
816 \r
817   if(fRunShuffling) {\r
818     if (!fShuffledBalance) {\r
819       Printf("ERROR: fShuffledBalance not available");\r
820       return;\r
821     }\r
822     fListBFS->Write();\r
823   }\r
824   gOutput->Close();\r
825 }\r
826 \r