]> 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
end-of-line normalization
[u/mrichter/AliRoot.git] / PWGCF / EBYE / BalanceFunctions / AliAnalysisTaskToyModel.cxx
CommitLineData
7fac8669 1#include "TChain.h"
2#include "TList.h"
3#include "TCanvas.h"
4#include "TParticle.h"
5#include "TLorentzVector.h"
6#include "TGraphErrors.h"
7#include "TH1.h"
8#include "TH2.h"
9#include "TH3.h"
10#include "TArrayF.h"
11#include "TF1.h"
12#include "TRandom.h"
13#include "TFile.h"
14
15#include "AliAnalysisManager.h"
16#include "AliLog.h"
17
18#include "AliEventPoolManager.h"
19
20#include "AliAnalysisTaskToyModel.h"
21#include "AliBalance.h"
22#include "AliBalancePsi.h"
23#include "AliAnalysisTaskTriggeredBF.h"
24
25
26// Analysis task for the toy model analysis
27// Authors: Panos.Christakoglou@nikhef.nl
28
29using std::cout;
30using std::endl;
31
32ClassImp(AliAnalysisTaskToyModel)
33
34//________________________________________________________________________
35AliAnalysisTaskToyModel::AliAnalysisTaskToyModel()
36: TObject(),
37 fUseDebug(kFALSE),
38 fBalance(0),
39 fRunShuffling(kFALSE), fShuffledBalance(0),
40 fRunMixing(kFALSE), fMixedBalance(0), fPoolMgr(0),
41 fList(0), fListBF(0), fListBFS(0), fListBFM(0),
42 fHistEventStats(0),
43 fHistNumberOfAcceptedParticles(0),
44 fHistReactionPlane(0),
45 fHistEtaTotal(0), fHistEta(0),
46 fHistEtaPhiPos(0), fHistEtaPhiNeg(0),
47 fHistRapidity(0),
48 fHistRapidityPions(0), fHistRapidityKaons(0), fHistRapidityProtons(0),
49 fHistPhi(0),
50 fHistPhiPions(0), fHistPhiKaons(0), fHistPhiProtons(0),
51 fHistPt(0),
52 fHistPtPions(0), fHistPtKaons(0), fHistPtProtons(0),
53 fTotalMultiplicityMean(100.), fTotalMultiplicitySigma(0.0),
54 fNetChargeMean(0.0), fNetChargeSigma(0.0),
55 fPtMin(0.0), fPtMax(100.0),
56 fEtaMin(-1.0), fEtaMax(1.0),
57 fUseAcceptanceParameterization(kFALSE), fAcceptanceParameterization(0),
58 fSimulateDetectorEffects(kFALSE),
59 fNumberOfInefficientSectors(0),
60 fInefficiencyFactorInPhi(1.0),
61 fNumberOfDeadSectors(0),
62 fEfficiencyDropNearEtaEdges(kFALSE),
63 fEfficiencyMatrix(0),
64 fUseAllCharges(kFALSE), fParticleMass(0.0),
65 fPtSpectraAllCharges(0), fTemperatureAllCharges(100.),
66 fReactionPlane(0.0),
67 fAzimuthalAngleAllCharges(0), fDirectedFlowAllCharges(0.0),
68 fEllipticFlowAllCharges(0.0), fTriangularFlowAllCharges(0.0),
69 fQuandrangularFlowAllCharges(0.0), fPentangularFlowAllCharges(0.0),
70 fPionPercentage(0.8), fPionMass(0.0),
71 fPtSpectraPions(0), fTemperaturePions(100.),
72 fAzimuthalAnglePions(0), fDirectedFlowPions(0.0),
73 fEllipticFlowPions(0.0), fTriangularFlowPions(0.0),
74 fQuandrangularFlowPions(0.0), fPentangularFlowPions(0.0),
75 fKaonPercentage(0.8), fKaonMass(0.0),
76 fPtSpectraKaons(0), fTemperatureKaons(100.),
77 fAzimuthalAngleKaons(0), fDirectedFlowKaons(0.0),
78 fEllipticFlowKaons(0.0), fTriangularFlowKaons(0.0),
79 fQuandrangularFlowKaons(0.0), fPentangularFlowKaons(0.0),
80 fProtonPercentage(0.8), fProtonMass(0.0),
81 fPtSpectraProtons(0), fTemperatureProtons(100.),
82 fAzimuthalAngleProtons(0), fDirectedFlowProtons(0.0),
83 fEllipticFlowProtons(0.0), fTriangularFlowProtons(0.0),
84 fQuandrangularFlowProtons(0.0), fPentangularFlowProtons(0.0),
85 fUseDynamicalCorrelations(kFALSE), fDynamicalCorrelationsPercentage(0.1),
86 fUseJets(kFALSE), fPtAssoc(0) {
87 // Constructor
88}
89
90//________________________________________________________________________
91AliAnalysisTaskToyModel::~AliAnalysisTaskToyModel() {
92 //Destructor
93 if(fUseAllCharges) {
94 delete fPtSpectraAllCharges;
95 delete fAzimuthalAngleAllCharges;
96 }
97 else {
98 delete fPtSpectraPions;
99 delete fAzimuthalAnglePions;
100 delete fPtSpectraKaons;
101 delete fAzimuthalAngleKaons;
102 delete fPtSpectraProtons;
103 delete fAzimuthalAngleProtons;
104 }
105 if(fUseJets) delete fPtAssoc;
106}
107
108//________________________________________________________________________
109void AliAnalysisTaskToyModel::Init() {
110 //Initialize objects
111 //==============gRandom Seed=======================//
112 gRandom->SetSeed(0); //seed is set to a random value (depending on machine clock)
113 //==============gRandom Seed=======================//
114
115 //==============Particles and spectra==============//
116 TParticle *pion = new TParticle();
117 pion->SetPdgCode(211);
118 fPionMass = pion->GetMass();
119
120 TParticle *kaon = new TParticle();
121 kaon->SetPdgCode(321);
122 fKaonMass = kaon->GetMass();
123
124 TParticle *proton = new TParticle();
125 proton->SetPdgCode(2212);
126 fProtonMass = proton->GetMass();
127
128 if(fUseAllCharges) {
129 fParticleMass = fPionMass;
130 fPtSpectraAllCharges = new TF1("fPtSpectraAllCharges","x*TMath::Exp(-TMath::Power([0]*[0]+x*x,0.5)/[1])",0.,5.);
131 fPtSpectraAllCharges->SetParName(0,"Mass");
132 fPtSpectraAllCharges->SetParName(1,"Temperature");
133 //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.);
134 //fPtSpectraAllCharges->SetParName(0,"pt0");
135 //fPtSpectraAllCharges->SetParName(1,"b");
136 }
137 else {
138 fPtSpectraPions = new TF1("fPtSpectraPions","x*TMath::Exp(-TMath::Power([0]*[0]+x*x,0.5)/[1])",0.,5.);
139 fPtSpectraPions->SetParName(0,"Mass");
140 fPtSpectraPions->SetParName(1,"Temperature");
141
142 fPtSpectraKaons = new TF1("fPtSpectraKaons","x*TMath::Exp(-TMath::Power([0]*[0]+x*x,0.5)/[1])",0.,5.);
143 fPtSpectraKaons->SetParName(0,"Mass");
144 fPtSpectraKaons->SetParName(1,"Temperature");
145
146 fPtSpectraProtons = new TF1("fPtSpectraProtons","x*TMath::Exp(-TMath::Power([0]*[0]+x*x,0.5)/[1])",0.,5.);
147 fPtSpectraProtons->SetParName(0,"Mass");
148 fPtSpectraProtons->SetParName(1,"Temperature");
149 }
150 //==============Particles and spectra==============//
151
152 //==============Flow values==============//
153 if(fUseAllCharges) {
154 if(fUseDebug) {
155 Printf("Directed flow: %lf",fDirectedFlowAllCharges);
156 Printf("Elliptic flow: %lf",fEllipticFlowAllCharges);
157 Printf("Triangular flow: %lf",fTriangularFlowAllCharges);
158 Printf("Quandrangular flow: %lf",fQuandrangularFlowAllCharges);
159 Printf("Pentangular flow: %lf",fPentangularFlowAllCharges);
160 }
161
162 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());
163 fAzimuthalAngleAllCharges->SetParName(0,"Reaction Plane");
164 fAzimuthalAngleAllCharges->SetParName(1,"Directed flow");
165 fAzimuthalAngleAllCharges->SetParName(2,"Elliptic flow");
166 fAzimuthalAngleAllCharges->SetParName(3,"Triangular flow");
167 fAzimuthalAngleAllCharges->SetParName(4,"Quandrangular flow");
168 fAzimuthalAngleAllCharges->SetParName(5,"Pentangular flow");
169 }
170 else {
171 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());
172 fAzimuthalAnglePions->SetParName(0,"Reaction Plane");
173 fAzimuthalAnglePions->SetParName(1,"Directed flow");
174 fAzimuthalAnglePions->SetParName(2,"Elliptic flow");
175 fAzimuthalAnglePions->SetParName(3,"Triangular flow");
176 fAzimuthalAnglePions->SetParName(4,"Quandrangular flow");
177 fAzimuthalAnglePions->SetParName(5,"Pentangular flow");
178
179 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());
180 fAzimuthalAngleKaons->SetParName(0,"Reaction Plane");
181 fAzimuthalAngleKaons->SetParName(1,"Directed flow");
182 fAzimuthalAngleKaons->SetParName(2,"Elliptic flow");
183 fAzimuthalAngleKaons->SetParName(3,"Triangular flow");
184 fAzimuthalAngleKaons->SetParName(4,"Quandrangular flow");
185 fAzimuthalAngleKaons->SetParName(5,"Pentangular flow");
186
187 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());
188 fAzimuthalAngleProtons->SetParName(0,"Reaction Plane");
189 fAzimuthalAngleProtons->SetParName(1,"Directed flow");
190 fAzimuthalAngleProtons->SetParName(2,"Elliptic flow");
191 fAzimuthalAngleProtons->SetParName(3,"Triangular flow");
192 fAzimuthalAngleProtons->SetParName(4,"Quandrangular flow");
193 fAzimuthalAngleProtons->SetParName(5,"Pentangular flow");
194 }
195 //==============Flow values==============//
196
197 //===================Jets===================//
198 if(fUseJets) {
199 fPtAssoc = new TF1("fPtAssoc","x*TMath::Exp(-TMath::Power([0]*[0]+x*x,0.5)/[1])",0.,20.);
200 fPtAssoc->SetParName(0,"pt0");
201 fPtAssoc->SetParName(1,"b");
202 fPtAssoc->SetParameter(0,0.139);
203 fPtAssoc->SetParameter(1,0.5);
204 fPtAssoc->SetLineColor(1);
205 }
206
207 //===================Jets===================//
208
209 //==============Efficiency matrix==============//
210 if(fSimulateDetectorEffects) SetupEfficiencyMatrix();
211 //==============Efficiency matrix==============//
212}
213
214//________________________________________________________________________
215void AliAnalysisTaskToyModel::SetupEfficiencyMatrix() {
216 //Setup the efficiency matrix
217 TH1F *hPt = new TH1F("hPt","",200,0.1,20.1);
218 TH1F *hEta = new TH1F("hEta","",20,-0.95,0.95);
219 TH1F *hPhi = new TH1F("hPhi","",72,0.,2.*TMath::Pi());
220 fEfficiencyMatrix = new TH3F("fEfficiencyMatrix","",
221 hEta->GetNbinsX(),
222 hEta->GetXaxis()->GetXmin(),
223 hEta->GetXaxis()->GetXmax(),
224 hPt->GetNbinsX(),
225 hPt->GetXaxis()->GetXmin(),
226 hPt->GetXaxis()->GetXmax(),
227 hPhi->GetNbinsX(),
228 hPhi->GetXaxis()->GetXmin(),
229 hPhi->GetXaxis()->GetXmax());
230
231 //Efficiency in pt
232 Double_t epsilon[20] = {0.3,0.6,0.77,0.79,0.80,0.80,0.80,0.80,0.80,0.80,0.80,0.80,0.80,0.80,0.80,0.80,0.80,0.80,0.80,0.80};
233 for(Int_t i=1;i<=20;i++) {
234 hPt->SetBinContent(i,epsilon[i-1]);
235 hPt->SetBinError(i,0.01);
236 }
237 for(Int_t i=21;i<=200;i++) {
238 hPt->SetBinContent(i,epsilon[19]);
239 hPt->SetBinError(i,0.01);
240 }
241
242 //Efficiency in eta
243 for(Int_t i=1;i<=hEta->GetNbinsX();i++) {
244 hEta->SetBinContent(i,1.0);
245 hEta->SetBinError(i,0.01);
246 }
247 if(fEfficiencyDropNearEtaEdges) {
248 hEta->SetBinContent(1,0.7); hEta->SetBinContent(2,0.8);
249 hEta->SetBinContent(3,0.9);
250 hEta->SetBinContent(18,0.9); hEta->SetBinContent(19,0.8);
251 hEta->SetBinContent(20,0.7);
252 }
253
254 //Efficiency in phi
255 for(Int_t i=1;i<=hPhi->GetNbinsX();i++) {
256 hPhi->SetBinContent(i,1.0);
257 hPhi->SetBinError(i,0.01);
258 }
259 for(Int_t i=1;i<=fNumberOfInefficientSectors;i++)
260 hPhi->SetBinContent(hPhi->FindBin(hPhi->GetRandom()),fInefficiencyFactorInPhi);
261 for(Int_t i=1;i<=fNumberOfDeadSectors;i++)
262 hPhi->SetBinContent(hPhi->FindBin(hPhi->GetRandom()),0.0);
263
264 //Fill the 3D efficiency map
265 for(Int_t iBinX = 1; iBinX<=fEfficiencyMatrix->GetXaxis()->GetNbins();iBinX++) {
266 //cout<<"==================================="<<endl;
267 for(Int_t iBinY = 1; iBinY<=fEfficiencyMatrix->GetYaxis()->GetNbins();iBinY++) {
268 //cout<<"==================================="<<endl;
269 for(Int_t iBinZ = 1; iBinZ<=fEfficiencyMatrix->GetZaxis()->GetNbins();iBinZ++) {
270 fEfficiencyMatrix->SetBinContent(iBinX,iBinY,iBinZ,hEta->GetBinContent(iBinX)*hPt->GetBinContent(iBinY)*hPhi->GetBinContent(iBinZ));
271 //cout<<"Eta: "<<hEta->GetBinCenter(iBinX)<<" - Pt: "<<hPt->GetBinCenter(iBinY)<<" - Phi: "<<hPhi->GetBinCenter(iBinZ)<<" - "<<hEta->GetBinContent(iBinX)<<" , "<<hPt->GetBinContent(iBinY)<<" , "<<hPhi->GetBinContent(iBinZ)<<" - Efficiency: "<<hEta->GetBinContent(iBinX)*hPt->GetBinContent(iBinY)*hPhi->GetBinContent(iBinZ)<<endl;
272 }
273 }
274 }
275}
276
277//________________________________________________________________________
278void AliAnalysisTaskToyModel::CreateOutputObjects() {
279 // Create histograms
280 // Called once
281
282 // global switch disabling the reference
283 // (to avoid "Replacing existing TH1" if several wagons are created in train)
284 Bool_t oldStatus = TH1::AddDirectoryStatus();
285 TH1::AddDirectory(kFALSE);
286
287 if(!fBalance) {
288 fBalance = new AliBalancePsi();
289 fBalance->SetDeltaEtaMax(2.0);
290 }
291 if(fRunShuffling) {
292 if(!fShuffledBalance) {
293 fShuffledBalance = new AliBalancePsi();
294 fShuffledBalance->SetDeltaEtaMax(2.0);
295 }
296 }
297 if(fRunMixing) {
298 if(!fMixedBalance) {
299 fMixedBalance = new AliBalancePsi();
300 fMixedBalance->SetDeltaEtaMax(2.0);
301 }
302 }
303
304 //QA list
305 fList = new TList();
306 fList->SetName("listQA");
307 fList->SetOwner();
308
309 //Balance Function list
310 fListBF = new TList();
311 fListBF->SetName("listBF");
312 fListBF->SetOwner();
313
314 if(fRunShuffling) {
315 fListBFS = new TList();
316 fListBFS->SetName("listBFShuffled");
317 fListBFS->SetOwner();
318 }
319
320 if(fRunMixing) {
321 fListBFM = new TList();
322 fListBFM->SetName("listBFMixed");
323 fListBFM->SetOwner();
324 }
325
326 //==============QA================//
327 //Event stats.
328 TString gCutName[4] = {"Total","Offline trigger",
329 "Vertex","Analyzed"};
330 fHistEventStats = new TH1F("fHistEventStats",
331 "Event statistics;;N_{events}",
332 4,0.5,4.5);
333 for(Int_t i = 1; i <= 4; i++)
334 fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data());
335 fList->Add(fHistEventStats);
336
337 fHistNumberOfAcceptedParticles = new TH1F("fHistNumberOfAcceptedParticles",";N_{acc.};Entries",10000,0,10000);
338 fList->Add(fHistNumberOfAcceptedParticles);
339
340 fHistReactionPlane = new TH1F("fHistReactionPlane","Reaction plane angle;#Psi [rad];Entries",1000,0.,2.*TMath::Pi());
341 fList->Add(fHistReactionPlane);
342
343 //Pseudo-rapidity
344 fHistEtaTotal = new TH1F("fHistEtaTotal","Pseudo-rapidity (full phase space);#eta;Entries",1000,-15.,15.);
345 fList->Add(fHistEtaTotal);
346
347 fHistEta = new TH1F("fHistEta","Pseudo-rapidity (acceptance);#eta;Entries",1000,-1.5,1.5);
348 fList->Add(fHistEta);
349
350 fHistEtaPhiPos = new TH2F("fHistEtaPhiPos","#eta-#phi distribution (+);#eta;#varphi (rad)",80,-2.,2.,72,-TMath::Pi()/2.,3.*TMath::Pi()/2.);
351 fList->Add(fHistEtaPhiPos);
352 fHistEtaPhiNeg = new TH2F("fHistEtaPhiNeg","#eta-#phi distribution (-);#eta;#varphi (rad)",80,-2.,2.,72,-TMath::Pi()/2.,3.*TMath::Pi()/2.);
353 fList->Add(fHistEtaPhiNeg);
354
355 //Rapidity
356 fHistRapidity = new TH1F("fHistRapidity","Rapidity (acceptance);y;Entries",1000,-1.5,1.5);
357 fList->Add(fHistRapidity);
358 fHistRapidityPions = new TH1F("fHistRapidityPions","Rapidity (acceptance - pions);y;Entries",1000,-1.5,1.5);
359 fList->Add(fHistRapidityPions);
360 fHistRapidityKaons = new TH1F("fHistRapidityKaons","Rapidity (acceptance - kaons);y;Entries",1000,-1.5,1.5);
361 fList->Add(fHistRapidityKaons);
362 fHistRapidityProtons = new TH1F("fHistRapidityProtons","Rapidity (acceptance - protons);y;Entries",1000,-1.5,1.5);
363 fList->Add(fHistRapidityProtons);
364
365 //Phi
366 fHistPhi = new TH1F("fHistPhi","Phi (acceptance);#phi (rad);Entries",1000,0.,2*TMath::Pi());
367 fList->Add(fHistPhi);
368
369 fHistPhiPions = new TH1F("fHistPhiPions","Phi (acceptance - pions);#phi (rad);Entries",1000,0.,2*TMath::Pi());
370 fList->Add(fHistPhiPions);
371 fHistPhiKaons = new TH1F("fHistPhiKaons","Phi (acceptance - kaons);#phi (rad);Entries",1000,0.,2*TMath::Pi());
372 fList->Add(fHistPhiKaons);
373 fHistPhiProtons = new TH1F("fHistPhiProtons","Phi (acceptance - protons);#phi (rad);Entries",1000,0.,2*TMath::Pi());
374 fList->Add(fHistPhiProtons);
375
376 //Pt
377 fHistPt = new TH1F("fHistPt","Pt (acceptance);p_{t} (GeV/c);Entries",1000,0.,10.);
378 fList->Add(fHistPt);
379
380 fHistPtPions = new TH1F("fHistPtPions","Pt (acceptance - pions);p_{t} (GeV/c);Entries",1000,0.,10.);
381 fList->Add(fHistPtPions);
382 fHistPtKaons = new TH1F("fHistPtKaons","Pt (acceptance - kaons);p_{t} (GeV/c);Entries",1000,0.,10.);
383 fList->Add(fHistPtKaons);
384 fHistPtProtons = new TH1F("fHistPtProtons","Pt (acceptance - protons);p_{t} (GeV/c);Entries",1000,0.,10.);
385 fList->Add(fHistPtProtons);
386
387 if(fEfficiencyMatrix) fList->Add(fEfficiencyMatrix);
388
389 //==============Balance function histograms================//
390 // Initialize histograms if not done yet
391 if(!fBalance->GetHistNp()){
392 AliWarning("Histograms not yet initialized! --> Will be done now");
393 AliWarning("--> Add 'gBalance->InitHistograms()' in your configBalanceFunction");
394 fBalance->InitHistograms();
395 }
396
397 if(fRunShuffling) {
398 if(!fShuffledBalance->GetHistNp()) {
399 AliWarning("Histograms (shuffling) not yet initialized! --> Will be done now");
400 AliWarning("--> Add 'gBalance->InitHistograms()' in your configBalanceFunction");
401 fShuffledBalance->InitHistograms();
402 }
403 }
404
405 fListBF->Add(fBalance->GetHistNp());
406 fListBF->Add(fBalance->GetHistNn());
407 fListBF->Add(fBalance->GetHistNpn());
408 fListBF->Add(fBalance->GetHistNnn());
409 fListBF->Add(fBalance->GetHistNpp());
410 fListBF->Add(fBalance->GetHistNnp());
411
412 if(fRunShuffling) {
413 fListBFS->Add(fShuffledBalance->GetHistNp());
414 fListBFS->Add(fShuffledBalance->GetHistNn());
415 fListBFS->Add(fShuffledBalance->GetHistNpn());
416 fListBFS->Add(fShuffledBalance->GetHistNnn());
417 fListBFS->Add(fShuffledBalance->GetHistNpp());
418 fListBFS->Add(fShuffledBalance->GetHistNnp());
419 }
420
421 if(fRunMixing) {
422 fListBFM->Add(fMixedBalance->GetHistNp());
423 fListBFM->Add(fMixedBalance->GetHistNn());
424 fListBFM->Add(fMixedBalance->GetHistNpn());
425 fListBFM->Add(fMixedBalance->GetHistNnn());
426 fListBFM->Add(fMixedBalance->GetHistNpp());
427 fListBFM->Add(fMixedBalance->GetHistNnp());
428 }
429
430 // Event Mixing
431 if(fRunMixing){
432 Int_t fMixingTracks = 2000;
433 Int_t trackDepth = fMixingTracks;
434 Int_t poolsize = 1000; // Maximum number of events, ignored in the present implemented of AliEventPoolManager
435
436 // centrality bins
437 // Double_t centralityBins[] = {0.,1.,2.,3.,4.,5.,7.,10.,20.,30.,40.,50.,60.,70.,80.,100.}; // SHOULD BE DEDUCED FROM CREATED ALITHN!!!
438 // Double_t* centbins = centralityBins;
439 // Int_t nCentralityBins = sizeof(centralityBins) / sizeof(Double_t) - 1;
440
441 // multiplicity bins
442 Double_t multiplicityBins[] = {0,10,20,30,40,50,60,70,80,100,100000}; // SHOULD BE DEDUCED FROM CREATED ALITHN!!!
443 Double_t* multbins = multiplicityBins;
444 Int_t nMultiplicityBins = sizeof(multiplicityBins) / sizeof(Double_t) - 1;
445
446 // Zvtx bins
447 Double_t vertexBins[] = {-10., -7., -5., -3., -1., 1., 3., 5., 7., 10.}; // SHOULD BE DEDUCED FROM CREATED ALITHN!!!
448 Double_t* vtxbins = vertexBins;
449 Int_t nVertexBins = sizeof(vertexBins) / sizeof(Double_t) - 1;
450
451 // Event plane angle (Psi) bins
452 // Double_t psiBins[] = {0.,45.,135.,215.,305.,360.}; // SHOULD BE DEDUCED FROM CREATED ALITHN!!!
453 // Double_t* psibins = psiBins;
454 // Int_t nPsiBins = sizeof(psiBins) / sizeof(Double_t) - 1;
455
456 // // run the event mixing also in bins of event plane (statistics!)
457 // if(fRunMixingEventPlane){
458 // if(fEventClass=="Multiplicity"){
459 // fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nMultiplicityBins, multbins, nVertexBins, vtxbins, nPsiBins, psibins);
460 // }
461 // else{
462 // fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centbins, nVertexBins, vtxbins, nPsiBins, psibins);
463 // }
464 // }
465 // else{
466 //if(fEventClass=="Multiplicity"){
467 fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nMultiplicityBins, multbins, nVertexBins, vtxbins);
468 //}
469 //else{
470 //fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centbins, nVertexBins, vtxbins);
471 //}
472 }
473
474 TH1::AddDirectory(oldStatus);
475}
476
477//________________________________________________________________________
478void AliAnalysisTaskToyModel::Run(Int_t nEvents) {
479 // Main loop
480 // Called for each event
481 Short_t vCharge = 0;
482 Float_t vY = 0.0;
483 Float_t vEta = 0.0;
484 Float_t vPhi = 0.0;
485 Float_t vP[3] = {0.,0.,0.};
486 Float_t vPt = 0.0;
487 Float_t vE = 0.0;
488 Bool_t isPion = kFALSE, isKaon = kFALSE, isProton = kFALSE;
489
490 Double_t gDecideCharge = 0.;
491
492 if(fUseAllCharges) {
493 //fPtSpectraAllCharges->SetParameter(0,fParticleMass);
494 fPtSpectraAllCharges->SetParameter(0,1.05);
495 fPtSpectraAllCharges->SetParameter(1,fTemperatureAllCharges);
496
497 fAzimuthalAngleAllCharges->SetParameter(1,fDirectedFlowAllCharges);
498 fAzimuthalAngleAllCharges->SetParameter(2,fEllipticFlowAllCharges);
499 fAzimuthalAngleAllCharges->SetParameter(3,fTriangularFlowAllCharges);
500 fAzimuthalAngleAllCharges->SetParameter(4,fQuandrangularFlowAllCharges);
501 fAzimuthalAngleAllCharges->SetParameter(5,fPentangularFlowAllCharges);
502 }
503 else {
504 fPtSpectraPions->SetParameter(0,fPionMass);
505 fPtSpectraPions->SetParameter(1,fTemperaturePions);
506 fPtSpectraKaons->SetParameter(0,fKaonMass);
507 fPtSpectraKaons->SetParameter(1,fTemperatureKaons);
508 fPtSpectraProtons->SetParameter(0,fProtonMass);
509 fPtSpectraProtons->SetParameter(1,fTemperatureProtons);
510
511 fAzimuthalAnglePions->SetParameter(1,fDirectedFlowPions);
512 fAzimuthalAnglePions->SetParameter(2,fEllipticFlowPions);
513 fAzimuthalAnglePions->SetParameter(3,fTriangularFlowPions);
514 fAzimuthalAnglePions->SetParameter(4,fQuandrangularFlowPions);
515 fAzimuthalAnglePions->SetParameter(5,fPentangularFlowPions);
516
517 fAzimuthalAngleKaons->SetParameter(1,fDirectedFlowKaons);
518 fAzimuthalAngleKaons->SetParameter(2,fEllipticFlowKaons);
519 fAzimuthalAngleKaons->SetParameter(3,fTriangularFlowKaons);
520 fAzimuthalAngleKaons->SetParameter(4,fQuandrangularFlowKaons);
521 fAzimuthalAngleKaons->SetParameter(5,fPentangularFlowKaons);
522
523 fAzimuthalAngleProtons->SetParameter(1,fDirectedFlowProtons);
524 fAzimuthalAngleProtons->SetParameter(2,fEllipticFlowProtons);
525 fAzimuthalAngleProtons->SetParameter(3,fTriangularFlowProtons);
526 fAzimuthalAngleProtons->SetParameter(4,fQuandrangularFlowProtons);
527 fAzimuthalAngleProtons->SetParameter(5,fPentangularFlowProtons);
528 }
529
530
531 for(Int_t iEvent = 0; iEvent < nEvents; iEvent++) {
532 // vector holding the charges/kinematics of all tracks (charge,y,eta,phi,p0,p1,p2,pt,E)
533 //vector<Double_t> *chargeVectorShuffle[9]; // this will be shuffled
534 //vector<Double_t> *chargeVector[9]; // original charge
535 //for(Int_t i = 0; i < 9; i++){
536 //chargeVectorShuffle[i] = new vector<Double_t>;
537 //chargeVector[i] = new vector<Double_t>;
538 //}
539
540 // TObjArray for the accepted particles
541 // (has to be done here, otherwise mxing with event pool does not work, overwriting pointers!)
542 TObjArray *tracksMain = new TObjArray();
543 tracksMain->SetOwner(kTRUE);
544 TObjArray *tracksMixing = 0x0;
545 if(fRunMixing) {
546 tracksMixing = new TObjArray();
547 tracksMixing->SetOwner(kTRUE);
548 }
549
550 tracksMain->Clear();
551 if(fRunMixing) tracksMixing->Clear();
552
553 fHistEventStats->Fill(1);
554 fHistEventStats->Fill(2);
555 fHistEventStats->Fill(3);
556
557 if((iEvent%1000) == 0)
558 cout<<"Event: "<<iEvent<<"/"<<nEvents<<endl;
559
560 //Multiplicities
561 Int_t nMultiplicity = (Int_t)(gRandom->Gaus(fTotalMultiplicityMean,fTotalMultiplicitySigma));
562 Int_t nNetCharge = (Int_t)(gRandom->Gaus(fNetChargeMean,fNetChargeSigma));
563
564 Int_t nGeneratedPositive = (Int_t)((nMultiplicity/2) + nNetCharge);
565 Int_t nGeneratedNegative = nMultiplicity - nGeneratedPositive;
566 if(fUseDebug)
567 Printf("Total multiplicity: %d - Generated positive: %d - Generated negative: %d",nMultiplicity,nGeneratedPositive,nGeneratedNegative);
568
569 //Randomization of the reaction plane
570 fReactionPlane = 2.0*TMath::Pi()*gRandom->Rndm();
571 //fReactionPlane = 0.0;
572 if(fUseAllCharges)
573 fAzimuthalAngleAllCharges->SetParameter(0,fReactionPlane);
574 else {
575 fAzimuthalAnglePions->SetParameter(0,fReactionPlane);
576 fAzimuthalAngleKaons->SetParameter(0,fReactionPlane);
577 fAzimuthalAngleProtons->SetParameter(0,fReactionPlane);
578 }
579
580 Int_t gNumberOfAcceptedParticles = 0;
581 Int_t gNumberOfAcceptedPositiveParticles = 0;
582 Int_t gNumberOfAcceptedNegativeParticles = 0;
583 Int_t nGeneratedPions = 0, nGeneratedKaons = 0, nGeneratedProtons = 0;
584
585 //Generate particles
586 for(Int_t iParticleCount = 0; iParticleCount < nMultiplicity; iParticleCount++) {
587 isPion = kFALSE; isKaon = kFALSE; isProton = kFALSE;
588 if(fUseDebug)
589 Printf("Generating positive: %d(%d)",iParticleCount+1,nGeneratedPositive);
590
591 //Pseudo-rapidity sampled from a Gaussian centered @ 0
592 vEta = gRandom->Gaus(0.0,4.0);
593
594 //Fill QA histograms (full phase space)
595 fHistEtaTotal->Fill(vEta);
596
597 //Decide the charge
598 gDecideCharge = gRandom->Rndm();
599 if(gDecideCharge <= 1.*nGeneratedPositive/nMultiplicity)
600 vCharge = 1;
601 else
602 vCharge = -1;
603
604 //Acceptance
605 if((vEta < fEtaMin) || (vEta > fEtaMax)) continue;
606
607 if(!fUseAllCharges) {
608 //Decide the specie
609 Double_t randomNumberSpecies = gRandom->Rndm();
610 if((randomNumberSpecies >= 0.0)&&(randomNumberSpecies < fPionPercentage)) {
611 nGeneratedPions += 1;
612 vPt = fPtSpectraPions->GetRandom();
613 vPhi = fAzimuthalAnglePions->GetRandom();
614 fParticleMass = fPionMass;
615 isPion = kTRUE;
616 }
617 else if((randomNumberSpecies >= fPionPercentage)&&(randomNumberSpecies < fPionPercentage + fKaonPercentage)) {
618 nGeneratedKaons += 1;
619 vPt = fPtSpectraKaons->GetRandom();
620 vPhi = fAzimuthalAngleKaons->GetRandom();
621 fParticleMass = fKaonMass;
622 isKaon = kTRUE;
623 }
624 else if((randomNumberSpecies >= fPionPercentage + fKaonPercentage)&&(randomNumberSpecies < fPionPercentage + fKaonPercentage + fProtonPercentage)) {
625 nGeneratedProtons += 1;
626 vPt = fPtSpectraProtons->GetRandom();
627 vPhi = fAzimuthalAngleProtons->GetRandom();
628 fParticleMass = fProtonMass;
629 isProton = kTRUE;
630 }
631 }
632 else {
633 vPt = fPtSpectraAllCharges->GetRandom();
634 vPhi = fAzimuthalAngleAllCharges->GetRandom();
635 }
636
637 vP[0] = vPt*TMath::Cos(vPhi);
638 vP[1] = vPt*TMath::Sin(vPhi);
639 vP[2] = vPt*TMath::SinH(vEta);
640 vE = TMath::Sqrt(TMath::Power(fParticleMass,2) +
641 TMath::Power(vP[0],2) +
642 TMath::Power(vP[1],2) +
643 TMath::Power(vP[2],2));
644
645 vY = 0.5*TMath::Log((vE + vP[2])/(vE - vP[2]));
646
647 //pt coverage
648 if((vPt < fPtMin) || (vPt > fPtMax)) continue;
649 //Printf("pt: %lf - mins: %lf - max: %lf",vPt,fPtMin,fPtMax);
650
651 //acceptance filter
652 if(fUseAcceptanceParameterization) {
653 Double_t gRandomNumberForAcceptance = gRandom->Rndm();
654 if(gRandomNumberForAcceptance > fAcceptanceParameterization->Eval(vPt))
655 continue;
656 }
657
658 //Detector effects
659 if(fSimulateDetectorEffects) {
660 Double_t randomNumber = gRandom->Rndm();
661 if(randomNumber > fEfficiencyMatrix->GetBinContent(fEfficiencyMatrix->FindBin(vEta,vPt,vPhi)))
662 continue;
663 }
664
665 //Fill QA histograms (acceptance);
666 if(vCharge > 0) {
667 gNumberOfAcceptedPositiveParticles += 1;
668 if(vPhi > 3.*TMath::Pi()/2.)
669 fHistEtaPhiPos->Fill(vEta,vPhi-2.*TMath::Pi());
670 else
671 fHistEtaPhiPos->Fill(vEta,vPhi);
672 }
673 else {
674 gNumberOfAcceptedNegativeParticles += 1;
675 if(vPhi > 3.*TMath::Pi()/2.)
676 fHistEtaPhiNeg->Fill(vEta,vPhi-2.*TMath::Pi());
677 else
678 fHistEtaPhiNeg->Fill(vEta,vPhi);
679 }
680
681 fHistEta->Fill(vEta);
682 fHistRapidity->Fill(vY);
683 fHistPhi->Fill(vPhi);
684 fHistPt->Fill(vPt);
685 if(isPion) {
686 fHistRapidityPions->Fill(vY);
687 fHistPhiPions->Fill(vPhi);
688 fHistPtPions->Fill(vPt);
689 }
690 else if(isKaon) {
691 fHistRapidityKaons->Fill(vY);
692 fHistPhiKaons->Fill(vPhi);
693 fHistPtKaons->Fill(vPt);
694 }
695 else if(isProton) {
696 fHistRapidityProtons->Fill(vY);
697 fHistPhiProtons->Fill(vPhi);
698 fHistPtProtons->Fill(vPt);
699 }
700 gNumberOfAcceptedParticles += 1;
701
702 // add the track to the TObjArray
703 tracksMain->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, 1.0));
704 if(fRunMixing)
705 tracksMixing->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, 1.0));
706 }//generated positive particle loop
707
708 //Jets
709 if(fUseJets) {
710 const Int_t nAssociated = 3;
711
712 Double_t gPtTrig1 = 0., gPtTrig2 = 0., gPtAssoc = 0.;
713 Double_t gPhiTrig1 = 0., gPhiTrig2 = 0., gPhiAssoc = 0.;
714 Double_t gEtaTrig1 = 0., gEtaTrig2 = 0., gEtaAssoc = 0.;
715 Short_t gChargeTrig1 = 0, gChargeTrig2 = 0, gChargeAssoc = 0;
716
717 Double_t gJetCone = 0.2;
718
719 //First leading particle
720 gPtTrig1 = gRandom->Uniform(3.,5.);
721 gEtaTrig1 = gRandom->Uniform(-0.8,0.8);
722 gPhiTrig1 = gRandom->Uniform(0.,TMath::TwoPi());
723
724 //Decide the charge
725 gDecideCharge = gRandom->Rndm();
726 if(gDecideCharge <= 0.5)
727 gChargeTrig1 = 1;
728 else
729 gChargeTrig1 = -1;
730
731 //Acceptance
732 if((gEtaTrig1 < fEtaMin) || (gEtaTrig1 > fEtaMax)) continue;
733 //pt coverage
734 if((gPtTrig1 < fPtMin) || (gPtTrig1 > fPtMax)) continue;
735 //Printf("pt: %lf - mins: %lf - max: %lf",vPt,fPtMin,fPtMax);
736
737 //acceptance filter
738 if(fUseAcceptanceParameterization) {
739 Double_t gRandomNumberForAcceptance = gRandom->Rndm();
740 if(gRandomNumberForAcceptance > fAcceptanceParameterization->Eval(gPtTrig1))
741 continue;
742 }
743
744 //Detector effects
745 if(fSimulateDetectorEffects) {
746 Double_t randomNumber = gRandom->Rndm();
747 if(randomNumber > fEfficiencyMatrix->GetBinContent(fEfficiencyMatrix->FindBin(gEtaTrig1,gPtTrig1,gPhiTrig1)))
748 continue;
749 }
750
751 gNumberOfAcceptedParticles += 1;
752
753 // add the track to the TObjArray
754 tracksMain->Add(new AliBFBasicParticle(gEtaTrig1, gPhiTrig1, gPtTrig1, gChargeTrig1, 1.0));
755 if(fRunMixing)
756 tracksMixing->Add(new AliBFBasicParticle(gEtaTrig1, gPhiTrig1, gPtTrig1, gChargeTrig1, 1.0));
757
758 Int_t iAssociated = 0;
759 while(iAssociated < nAssociated) {
760 gPtAssoc = fPtAssoc->GetRandom();
761 if(gPtAssoc < gPtTrig1) {
762 gEtaAssoc = gRandom->Uniform(gEtaTrig1 - gJetCone/2.,gEtaTrig1 + gJetCone/2.);
763 gPhiAssoc = gRandom->Uniform(gPhiTrig1 - gJetCone/2.,gPhiTrig1 + gJetCone/2.);
764 if(gPhiAssoc < 0.) gPhiAssoc += TMath::TwoPi();
765 else if(gPhiAssoc > TMath::TwoPi()) gPhiAssoc -= TMath::TwoPi();
766
767 iAssociated += 1;
768
769 gDecideCharge = gRandom->Rndm();
770 if(gDecideCharge <= 0.5)
771 gChargeAssoc = 1;
772 else
773 gChargeAssoc = -1;
774
775 //Acceptance
776 if((gEtaAssoc < fEtaMin) || (gEtaAssoc > fEtaMax)) continue;
777 //pt coverage
778 if((gPtAssoc < fPtMin) || (gPtAssoc > fPtMax)) continue;
779 //Printf("pt: %lf - mins: %lf - max: %lf",vPt,fPtMin,fPtMax);
780
781 //acceptance filter
782 if(fUseAcceptanceParameterization) {
783 Double_t gRandomNumberForAcceptance = gRandom->Rndm();
784 if(gRandomNumberForAcceptance > fAcceptanceParameterization->Eval(gPtAssoc))
785 continue;
786 }
787
788 //Detector effects
789 if(fSimulateDetectorEffects) {
790 Double_t randomNumber = gRandom->Rndm();
791 if(randomNumber > fEfficiencyMatrix->GetBinContent(fEfficiencyMatrix->FindBin(gEtaAssoc,gPtAssoc,gPhiAssoc)))
792 continue;
793 }
794
795 gNumberOfAcceptedParticles += 1;
796
797 // add the track to the TObjArray
798 tracksMain->Add(new AliBFBasicParticle(gEtaAssoc, gPhiAssoc, gPtAssoc, gChargeAssoc, 1.0));
799 if(fRunMixing)
800 tracksMixing->Add(new AliBFBasicParticle(gEtaAssoc, gPhiAssoc, gPtAssoc, gChargeAssoc, 1.0));
801 }//pt,assoc < pt,trig
802 }//associated
803
804 //back2back
805 gPtTrig2 = gPtTrig1;
806 gEtaTrig2 = -gEtaTrig1;
807 gPhiTrig2 = TMath::Pi() + gPhiTrig1;
808 if(gPhiTrig2 < 0.) gPhiTrig2 += TMath::TwoPi();
809 else if(gPhiTrig2 > TMath::TwoPi()) gPhiTrig2 -= TMath::TwoPi();
810
811 //Decide the charge
812 gDecideCharge = gRandom->Rndm();
813 if(gDecideCharge <= 0.5)
814 gChargeTrig2 = 1;
815 else
816 gChargeTrig2 = -1;
817
818 //Acceptance
819 if((gEtaTrig2 < fEtaMin) || (gEtaTrig2 > fEtaMax)) continue;
820 //pt coverage
821 if((gPtTrig2 < fPtMin) || (gPtTrig2 > fPtMax)) continue;
822 //Printf("pt: %lf - mins: %lf - max: %lf",vPt,fPtMin,fPtMax);
823
824 //acceptance filter
825 if(fUseAcceptanceParameterization) {
826 Double_t gRandomNumberForAcceptance = gRandom->Rndm();
827 if(gRandomNumberForAcceptance > fAcceptanceParameterization->Eval(gPtTrig2))
828 continue;
829 }
830
831 //Detector effects
832 if(fSimulateDetectorEffects) {
833 Double_t randomNumber = gRandom->Rndm();
834 if(randomNumber > fEfficiencyMatrix->GetBinContent(fEfficiencyMatrix->FindBin(gEtaTrig2,gPtTrig2,gPhiTrig2)))
835 continue;
836 }
837
838 gNumberOfAcceptedParticles += 1;
839
840 // add the track to the TObjArray
841 tracksMain->Add(new AliBFBasicParticle(gEtaTrig2, gPhiTrig2, gPtTrig2, gChargeTrig2, 1.0));
842 if(fRunMixing)
843 tracksMixing->Add(new AliBFBasicParticle(gEtaTrig2, gPhiTrig2, gPtTrig2, gChargeTrig2, 1.0));
844
845 iAssociated = 0;
846 while(iAssociated < nAssociated) {
847 gPtAssoc = fPtAssoc->GetRandom();
848 if(gPtAssoc < gPtTrig2) {
849 gEtaAssoc = gRandom->Uniform(gEtaTrig2 - gJetCone/2.,gEtaTrig2 + gJetCone/2.);
850 gPhiAssoc = gRandom->Uniform(gPhiTrig2 - gJetCone/2.,gPhiTrig2 + gJetCone/2.);
851 if(gPhiAssoc < 0.) gPhiAssoc += TMath::TwoPi();
852 else if(gPhiAssoc > TMath::TwoPi()) gPhiAssoc -= TMath::TwoPi();
853
854 iAssociated += 1;
855
856 gDecideCharge = gRandom->Rndm();
857 if(gDecideCharge <= 0.5)
858 gChargeAssoc = 1;
859 else
860 gChargeAssoc = -1;
861
862 //Acceptance
863 if((gEtaAssoc < fEtaMin) || (gEtaAssoc > fEtaMax)) continue;
864 //pt coverage
865 if((gPtAssoc < fPtMin) || (gPtAssoc > fPtMax)) continue;
866 //Printf("pt: %lf - mins: %lf - max: %lf",vPt,fPtMin,fPtMax);
867
868 //acceptance filter
869 if(fUseAcceptanceParameterization) {
870 Double_t gRandomNumberForAcceptance = gRandom->Rndm();
871 if(gRandomNumberForAcceptance > fAcceptanceParameterization->Eval(gPtAssoc))
872 continue;
873 }
874
875 //Detector effects
876 if(fSimulateDetectorEffects) {
877 Double_t randomNumber = gRandom->Rndm();
878 if(randomNumber > fEfficiencyMatrix->GetBinContent(fEfficiencyMatrix->FindBin(gEtaAssoc,gPtAssoc,gPhiAssoc)))
879 continue;
880 }
881
882 gNumberOfAcceptedParticles += 1;
883
884 // add the track to the TObjArray
885 tracksMain->Add(new AliBFBasicParticle(gEtaAssoc, gPhiAssoc, gPtAssoc, gChargeAssoc, 1.0));
886 if(fRunMixing)
887 tracksMixing->Add(new AliBFBasicParticle(gEtaAssoc, gPhiAssoc, gPtAssoc, gChargeAssoc, 1.0));
888 }//pt,assoc < pt,trig
889 }//associated
890 }//Jet usage
891
892
893
894
895
896
897
898
899
900 //Dynamical correlations
901 Int_t nGeneratedPositiveDynamicalCorrelations = 0;
902 Int_t nGeneratedNegativeDynamicalCorrelations = 0;
903 /*Double_t vChargePrime = 0;
904 Double_t vYPrime = 0.0;
905 Double_t vEtaPrime = 0.0;
906 Double_t vPhiPrime = 0.0;
907 Double_t vPPrime[3] = {0.,0.,0.};
908 Double_t vPtPrime = 0.0;
909 Double_t vEPrime = 0.0;
910 //Generate "correlated" particles
911 if(fUseDynamicalCorrelations) {
912 Int_t gNumberOfDynamicalCorrelations = (Int_t)(0.5*gNumberOfAcceptedParticles*fDynamicalCorrelationsPercentage);
913 for(Int_t iDynamicalCorrelations = 0; iDynamicalCorrelations < gNumberOfDynamicalCorrelations; iDynamicalCorrelations++) {
914 isPion = kFALSE; isKaon = kFALSE; isProton = kFALSE;
915
916 //Pseudo-rapidity sampled from a Gaussian centered @ 0
917 vEta = gRandom->Gaus(0.0,0.1);
918 vCharge = 1.0;
919 nGeneratedPositiveDynamicalCorrelations += 1;
920
921 vEtaPrime = -vEta;
922 vChargePrime = -1.0;
923 nGeneratedNegativeDynamicalCorrelations += 1;
924
925 //Acceptance
926 if((vEta < fEtaMin) || (vEta > fEtaMax)) continue;
927 if((vEtaPrime < fEtaMin) || (vEtaPrime > fEtaMax)) continue;
928
929 if(!fUseAllCharges) {
930 //Decide the specie
931 Double_t randomNumberSpecies = gRandom->Rndm();
932 if((randomNumberSpecies >= 0.0)&&(randomNumberSpecies < fPionPercentage)) {
933 nGeneratedPions += 1;
934 vPt = fPtSpectraPions->GetRandom();
935 vPhi = fAzimuthalAnglePions->GetRandom();
936 fParticleMass = fPionMass;
937 isPion = kTRUE;
938 }
939 else if((randomNumberSpecies >= fPionPercentage)&&(randomNumberSpecies < fPionPercentage + fKaonPercentage)) {
940 nGeneratedKaons += 1;
941 vPt = fPtSpectraKaons->GetRandom();
942 vPhi = fAzimuthalAngleKaons->GetRandom();
943 fParticleMass = fKaonMass;
944 isKaon = kTRUE;
945 }
946 else if((randomNumberSpecies >= fPionPercentage + fKaonPercentage)&&(randomNumberSpecies < fPionPercentage + fKaonPercentage + fProtonPercentage)) {
947 nGeneratedProtons += 1;
948 vPt = fPtSpectraProtons->GetRandom();
949 vPtPrime = vPt;
950 vPhi = fAzimuthalAngleProtons->GetRandom();
951 fParticleMass = fProtonMass;
952 isProton = kTRUE;
953 }
954 }
955 else {
956 vPt = fPtSpectraAllCharges->GetRandom();
957 vPhi = fAzimuthalAngleAllCharges->GetRandom();
958 }
959 vPtPrime = vPt;
960 vPhiPrime = vPhi;
961
962 vP[0] = vPt*TMath::Cos(vPhi);
963 vP[1] = vPt*TMath::Sin(vPhi);
964 vP[2] = vPt*TMath::SinH(vEta);
965 vE = TMath::Sqrt(TMath::Power(fParticleMass,2) +
966 TMath::Power(vP[0],2) +
967 TMath::Power(vP[1],2) +
968 TMath::Power(vP[2],2));
969
970 vY = 0.5*TMath::Log((vE + vP[2])/(vE - vP[2]));
971
972 vPPrime[0] = vPtPrime*TMath::Cos(vPhiPrime);
973 vPPrime[1] = vPtPrime*TMath::Sin(vPhiPrime);
974 vPPrime[2] = vPtPrime*TMath::SinH(vEtaPrime);
975 vEPrime = TMath::Sqrt(TMath::Power(fParticleMass,2) +
976 TMath::Power(vPPrime[0],2) +
977 TMath::Power(vPPrime[1],2) +
978 TMath::Power(vPPrime[2],2));
979
980 vYPrime = 0.5*TMath::Log((vEPrime + vPPrime[2])/(vEPrime - vPPrime[2]));
981
982 //pt coverage
983 if((vPt < fPtMin) || (vPt > fPtMax)) continue;
984 if((vPtPrime < fPtMin) || (vPtPrime > fPtMax)) continue;
985
986 //acceptance filter
987 if(fUseAcceptanceParameterization) {
988 Double_t gRandomNumberForAcceptance = gRandom->Rndm();
989 if(gRandomNumberForAcceptance > fAcceptanceParameterization->Eval(vPt))
990 continue;
991
992 Double_t gRandomNumberForAcceptancePrime = gRandom->Rndm();
993 if(gRandomNumberForAcceptancePrime > fAcceptanceParameterization->Eval(vPtPrime))
994 continue;
995 }
996
997 // fill charge vector (positive)
998 chargeVector[0]->push_back(vCharge);
999 chargeVector[1]->push_back(vY);
1000 chargeVector[2]->push_back(vEta);
1001 chargeVector[3]->push_back(vPhi);
1002 chargeVector[4]->push_back(vP[0]);
1003 chargeVector[5]->push_back(vP[1]);
1004 chargeVector[6]->push_back(vP[2]);
1005 chargeVector[7]->push_back(vPt);
1006 chargeVector[8]->push_back(vE);
1007
1008 if(fRunShuffling) {
1009 chargeVectorShuffle[0]->push_back(vCharge);
1010 chargeVectorShuffle[1]->push_back(vY);
1011 chargeVectorShuffle[2]->push_back(vEta);
1012 chargeVectorShuffle[3]->push_back(vPhi);
1013 chargeVectorShuffle[4]->push_back(vP[0]);
1014 chargeVectorShuffle[5]->push_back(vP[1]);
1015 chargeVectorShuffle[6]->push_back(vP[2]);
1016 chargeVectorShuffle[7]->push_back(vPt);
1017 chargeVectorShuffle[8]->push_back(vE);
1018 }
1019
1020 // fill charge vector (negative)
1021 chargeVector[0]->push_back(vChargePrime);
1022 chargeVector[1]->push_back(vYPrime);
1023 chargeVector[2]->push_back(vEtaPrime);
1024 chargeVector[3]->push_back(vPhiPrime);
1025 chargeVector[4]->push_back(vPPrime[0]);
1026 chargeVector[5]->push_back(vPPrime[1]);
1027 chargeVector[6]->push_back(vPPrime[2]);
1028 chargeVector[7]->push_back(vPtPrime);
1029 chargeVector[8]->push_back(vEPrime);
1030
1031 if(fRunShuffling) {
1032 chargeVectorShuffle[0]->push_back(vChargePrime);
1033 chargeVectorShuffle[1]->push_back(vYPrime);
1034 chargeVectorShuffle[2]->push_back(vEtaPrime);
1035 chargeVectorShuffle[3]->push_back(vPhiPrime);
1036 chargeVectorShuffle[4]->push_back(vPPrime[0]);
1037 chargeVectorShuffle[5]->push_back(vPPrime[1]);
1038 chargeVectorShuffle[6]->push_back(vPPrime[2]);
1039 chargeVectorShuffle[7]->push_back(vPtPrime);
1040 chargeVectorShuffle[8]->push_back(vEPrime);
1041 }
1042
1043 gNumberOfAcceptedParticles += 2;
1044 }//loop over the dynamical correlations
1045 }*/// usage of dynamical correlations
1046
1047 if(fUseDebug) {
1048 Printf("=======================================================");
1049 Printf("Total: %d - Total positive: %d - Total negative: %d",nMultiplicity,nGeneratedPositive,nGeneratedNegative);
1050 Printf("Accepted positive: %d - Accepted negative: %d",gNumberOfAcceptedPositiveParticles,gNumberOfAcceptedNegativeParticles);
1051 Printf("Correlations: %d - Correlations positive: %d - Correlations negative: %d",nGeneratedPositiveDynamicalCorrelations+nGeneratedNegativeDynamicalCorrelations,nGeneratedPositiveDynamicalCorrelations,nGeneratedNegativeDynamicalCorrelations);
1052 Printf("Number of accepted particles: %d",gNumberOfAcceptedParticles);
1053 if(!fUseAllCharges)
1054 Printf("Pions: %lf - Kaons: %lf - Protons: %lf",1.*nGeneratedPions/nMultiplicity,1.*nGeneratedKaons/nMultiplicity,1.*nGeneratedProtons/nMultiplicity);
1055 //Printf("Calculating the balance function for %d particles",chargeVector[0]->size());
1056 }
1057
1058 fHistEventStats->Fill(4);
1059 fHistNumberOfAcceptedParticles->Fill(gNumberOfAcceptedParticles);
1060 fHistReactionPlane->Fill(fReactionPlane);
1061
1062 // Event mixing
1063 if (fRunMixing)
1064 {
1065 // 1. First get an event pool corresponding in mult (cent) and
1066 // zvertex to the current event. Once initialized, the pool
1067 // should contain nMix (reduced) events. This routine does not
1068 // pre-scan the chain. The first several events of every chain
1069 // will be skipped until the needed pools are filled to the
1070 // specified depth. If the pool categories are not too rare, this
1071 // should not be a problem. If they are rare, you could lose`
1072 // statistics.
1073
1074 // 2. Collect the whole pool's content of tracks into one TObjArray
1075 // (bgTracks), which is effectively a single background super-event.
1076
1077 // 3. The reduced and bgTracks arrays must both be passed into
1078 // FillCorrelations(). Also nMix should be passed in, so a weight
1079 // of 1./nMix can be applied.
1080
1081 AliEventPool* pool = fPoolMgr->GetEventPool(1., 0.,fReactionPlane);
1082
1083 if (!pool){
1084 AliFatal(Form("No pool found for centrality = %f, zVtx = %f, psi = %f", 1., 0.,fReactionPlane));
1085 }
1086 else{
1087
1088 //pool->SetDebug(1);
1089
1090 if (pool->IsReady() || pool->NTracksInPool() > 50000 / 10 || pool->GetCurrentNEvents() >= 5){
1091
1092
1093 Int_t nMix = pool->GetCurrentNEvents();
1094 //cout << "nMix = " << nMix << " tracks in pool = " << pool->NTracksInPool() << ", tracks in this event = "<<gNumberOfAcceptedParticles <<", event Plane = "<<fReactionPlane<<endl;
1095
1096 //((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(2);
1097 //((TH2F*) fListOfHistos->FindObject("mixedDist"))->Fill(centrality, pool->NTracksInPool());
1098 //if (pool->IsReady())
1099 //((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(3);
1100
1101 // Fill mixed-event histos here
1102 for (Int_t jMix=0; jMix<nMix; jMix++)
1103 {
1104 TObjArray* tracksMixed = pool->GetEvent(jMix);
1105 fMixedBalance->CalculateBalance(fReactionPlane,tracksMain,tracksMixed,1,1.,0.);
1106 }
1107 }
1108
1109 // Update the Event pool
1110 pool->UpdatePool(tracksMain);
1111 //pool->PrintInfo();
1112
1113 }//pool NULL check
1114 }//run mixing
1115
1116 //Calculate the balance function
1117 fBalance->CalculateBalance(fReactionPlane,tracksMain,NULL,1,1.,0.);
1118 //if(fRunMixing)
1119 // fMixedBalance->CalculateBalance(fReactionPlane,tracksMixing,NULL,1,1.,0.);
1120
1121 }//event loop
1122}
1123
1124//________________________________________________________________________
1125void AliAnalysisTaskToyModel::FinishOutput() {
1126 //Printf("END BF");
1127
1128 TFile *gOutput = TFile::Open("outputToyModel.root","recreate");
1129 TDirectoryFile *dir = new TDirectoryFile("PWGCFEbyE.outputBalanceFunctionPsiAnalysis","PWGCFEbyE.outputBalanceFunctionPsiAnalysis");
1130
1131 fList->SetName("listQA");
1132 fList->SetOwner(kTRUE);
1133 dir->Add(fList);
1134 //fList->Write("listQA",TObject::kSingleKey);
1135
1136 if (!fBalance) {
1137 AliError("ERROR: fBalance not available");
1138 return;
1139 }
1140 fListBF->SetName("listBF");
1141 fListBF->SetOwner(kTRUE);
1142 dir->Add(fListBF);
1143 //fListBF->Write("listBF",TObject::kSingleKey);
1144
1145 if(fRunShuffling) {
1146 if (!fShuffledBalance) {
1147 AliError("ERROR: fShuffledBalance not available");
1148 return;
1149 }
1150 fListBFS->SetName("listBFShuffled");
1151 fListBFS->SetOwner(kTRUE);
1152 dir->Add(fListBFS);
1153 //fListBFS->Write("listBFShuffled",TObject::kSingleKey);
1154 }
1155
1156 if(fRunMixing) {
1157 if (!fMixedBalance) {
1158 AliError("ERROR: fMixedBalance not available");
1159 return;
1160 }
1161 fListBFM->SetName("listBFMixed");
1162 fListBFM->SetOwner(kTRUE);
1163 dir->Add(fListBFM);
1164 //fListBFM->Write("listBFMixed",TObject::kSingleKey);
1165 }
1166
1167 dir->Write(dir->GetName(),TObject::kSingleKey);
1168 gOutput->Close();
1169}
1170