end-of-line normalization
[u/mrichter/AliRoot.git] / PWGCF / EBYE / BalanceFunctions / AliAnalysisTaskToyModel.cxx
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
29 using std::cout;
30 using std::endl;
31
32 ClassImp(AliAnalysisTaskToyModel)
33
34 //________________________________________________________________________
35 AliAnalysisTaskToyModel::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 //________________________________________________________________________
91 AliAnalysisTaskToyModel::~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 //________________________________________________________________________
109 void 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 //________________________________________________________________________
215 void 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 //________________________________________________________________________
278 void 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 //________________________________________________________________________
478 void 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 //________________________________________________________________________
1125 void  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