4 #include "TParticle.h"
\r
5 #include "TLorentzVector.h"
\r
6 #include "TGraphErrors.h"
\r
11 #include "TRandom.h"
\r
13 #include "AliAnalysisTaskSE.h"
\r
14 #include "AliAnalysisManager.h"
\r
17 #include "AliAnalysisTaskToyModel.h"
\r
18 #include "AliBalance.h"
\r
21 // Analysis task for the toy model analysis
\r
22 // Authors: Panos.Christakoglou@nikhef.nl
\r
24 ClassImp(AliAnalysisTaskToyModel)
\r
26 //________________________________________________________________________
\r
27 AliAnalysisTaskToyModel::AliAnalysisTaskToyModel(const char *name)
\r
28 : AliAnalysisTaskSE(name),
\r
30 fRunShuffling(kFALSE), fShuffledBalance(0),
\r
31 fList(0), fListBF(0), fListBFS(0),
\r
33 fTotalMultiplicityMean(100.), fTotalMultiplicitySigma(0.0),
\r
34 fNetChargeMean(0.0), fNetChargeSigma(0.0),
\r
35 fPtMin(0.0), fPtMax(100.0),
\r
36 fEtaMin(-1.0), fEtaMax(1.0),
\r
37 fUseAcceptanceParameterization(kFALSE), fAcceptanceParameterization(0),
\r
38 fPtSpectraAllCharges(0), fTemperatureAllCharges(100.),
\r
39 fReactionPlane(0.0),
\r
40 fAzimuthalAngleAllCharges(0), fDirectedFlowAllCharges(0.0),
\r
41 fEllipticFlowAllCharges(0.0), fTriangularFlowAllCharges(0.0),
\r
42 fQuandrangularFlowAllCharges(0.0), fPentangularFlowAllCharges(0.0),
\r
43 fPionPercentage(0.8),
\r
44 fPtSpectraPions(0), fTemperaturePions(100.),
\r
45 fAzimuthalAnglePions(0), fDirectedFlowPions(0.0),
\r
46 fEllipticFlowPions(0.0), fTriangularFlowPions(0.0),
\r
47 fQuandrangularFlowPions(0.0), fPentangularFlowPions(0.0),
\r
48 fKaonPercentage(0.8),
\r
49 fPtSpectraKaons(0), fTemperatureKaons(100.),
\r
50 fAzimuthalAngleKaons(0), fDirectedFlowKaons(0.0),
\r
51 fEllipticFlowKaons(0.0), fTriangularFlowKaons(0.0),
\r
52 fQuandrangularFlowKaons(0.0), fPentangularFlowKaons(0.0),
\r
53 fProtonPercentage(0.8),
\r
54 fPtSpectraProtons(0), fTemperatureProtons(100.),
\r
55 fAzimuthalAngleProtons(0), fDirectedFlowProtons(0.0),
\r
56 fEllipticFlowProtons(0.0), fTriangularFlowProtons(0.0),
\r
57 fQuandrangularFlowProtons(0.0), fPentangularFlowProtons(0.0) {
\r
60 // Define input and output slots here
\r
61 // Input slot #0 works with a TChain
\r
62 DefineInput(0, TChain::Class());
\r
63 // Output slot #0 writes into a TH1 container
\r
64 DefineOutput(1, TList::Class());
\r
65 DefineOutput(2, TList::Class());
\r
66 DefineOutput(3, TList::Class());
\r
69 //________________________________________________________________________
\r
70 AliAnalysisTaskToyModel::~AliAnalysisTaskToyModel() {
\r
72 delete fPtSpectraAllCharges;
\r
73 delete fAzimuthalAngleAllCharges;
\r
74 delete fPtSpectraPions;
\r
75 delete fAzimuthalAnglePions;
\r
76 delete fPtSpectraKaons;
\r
77 delete fAzimuthalAngleKaons;
\r
78 delete fPtSpectraProtons;
\r
79 delete fAzimuthalAngleProtons;
\r
82 //________________________________________________________________________
\r
83 void AliAnalysisTaskToyModel::Init() {
\r
84 //Initialize objects
\r
85 //==============Particles and spectra==============//
\r
86 TParticle *pion = new TParticle();
\r
87 pion->SetPdgCode(211);
\r
88 Double_t gPionMass = pion->GetMass();
\r
90 fPtSpectraAllCharges = new TF1("fPtSpectraAllCharges","x*TMath::Exp(-TMath::Power([0]*[0]+x*x,0.5)/[1])",0.,100.);
\r
91 fPtSpectraAllCharges->SetParName(0,"Mass");
\r
92 fPtSpectraAllCharges->SetParameter(0,gPionMass);
\r
93 fPtSpectraAllCharges->SetParName(1,"Temperature");
\r
94 fPtSpectraAllCharges->SetParameter(1,fTemperatureAllCharges);
\r
96 fPtSpectraPions = new TF1("fPtSpectraPions","x*TMath::Exp(-TMath::Power([0]*[0]+x*x,0.5)/[1])",0.,100.);
\r
97 fPtSpectraPions->SetParName(0,"Mass");
\r
98 fPtSpectraPions->SetParameter(0,gPionMass);
\r
99 fPtSpectraPions->SetParName(1,"Temperature");
\r
100 fPtSpectraPions->SetParameter(1,fTemperaturePions);
\r
102 TParticle *kaon = new TParticle();
\r
103 kaon->SetPdgCode(321);
\r
104 Double_t gKaonMass = kaon->GetMass();
\r
105 fPtSpectraKaons = new TF1("fPtSpectraKaons","x*TMath::Exp(-TMath::Power([0]*[0]+x*x,0.5)/[1])",0.,100.);
\r
106 fPtSpectraKaons->SetParName(0,"Mass");
\r
107 fPtSpectraKaons->SetParameter(0,gKaonMass);
\r
108 fPtSpectraKaons->SetParName(1,"Temperature");
\r
109 fPtSpectraKaons->SetParameter(1,fTemperatureKaons);
\r
111 TParticle *proton = new TParticle();
\r
112 proton->SetPdgCode(2212);
\r
113 Double_t gProtonMass = proton->GetMass();
\r
114 fPtSpectraProtons = new TF1("fPtSpectraProtons","x*TMath::Exp(-TMath::Power([0]*[0]+x*x,0.5)/[1])",0.,5.);
\r
115 fPtSpectraProtons->SetParName(0,"Mass");
\r
116 fPtSpectraProtons->SetParameter(0,gProtonMass);
\r
117 fPtSpectraProtons->SetParName(1,"Temperature");
\r
118 fPtSpectraProtons->SetParameter(1,fTemperatureProtons);
\r
119 //==============Particles and spectra==============//
\r
121 //==============Flow values==============//
\r
122 fAzimuthalAngleAllCharges = new TF1("fAzimuthalAngleAllCharges","1+2.*[1]*TMath::Cos(x-[0])+2.*[2]*TMath::Cos(2*(x-[0]))+2.*[3]*TMath::Cos(3*(x-[0]))+2.*[4]*TMath::Cos(4*(x-[0]))+2.*[5]*TMath::Cos(5*(x-[0]))",0.,2.*TMath::Pi());
\r
123 fAzimuthalAngleAllCharges->SetParName(0,"Reaction Plane");
\r
124 fAzimuthalAngleAllCharges->SetParameter(0,fReactionPlane);
\r
125 fAzimuthalAngleAllCharges->SetParName(1,"Directed flow");
\r
126 fAzimuthalAngleAllCharges->SetParameter(1,fDirectedFlowAllCharges);
\r
127 fAzimuthalAngleAllCharges->SetParName(2,"Elliptic flow");
\r
128 fAzimuthalAngleAllCharges->SetParameter(2,fEllipticFlowAllCharges);
\r
129 fAzimuthalAngleAllCharges->SetParName(3,"Triangular flow");
\r
130 fAzimuthalAngleAllCharges->SetParameter(3,fTriangularFlowAllCharges);
\r
131 fAzimuthalAngleAllCharges->SetParName(4,"Quandrangular flow");
\r
132 fAzimuthalAngleAllCharges->SetParameter(4,fQuandrangularFlowAllCharges);
\r
133 fAzimuthalAngleAllCharges->SetParName(5,"Pentangular flow");
\r
134 fAzimuthalAngleAllCharges->SetParameter(5,fPentangularFlowAllCharges);
\r
136 fAzimuthalAnglePions = new TF1("fAzimuthalAnglePions","1+2.*[1]*TMath::Cos(x-[0])+2.*[2]*TMath::Cos(2*(x-[0]))+2.*[3]*TMath::Cos(3*(x-[0]))+2.*[4]*TMath::Cos(4*(x-[0]))+2.*[5]*TMath::Cos(5*(x-[0]))",0.,2.*TMath::Pi());
\r
137 fAzimuthalAnglePions->SetParName(0,"Reaction Plane");
\r
138 fAzimuthalAnglePions->SetParameter(0,fReactionPlane);
\r
139 fAzimuthalAnglePions->SetParName(1,"Directed flow");
\r
140 fAzimuthalAnglePions->SetParameter(1,fDirectedFlowPions);
\r
141 fAzimuthalAnglePions->SetParName(2,"Elliptic flow");
\r
142 fAzimuthalAnglePions->SetParameter(2,fEllipticFlowPions);
\r
143 fAzimuthalAnglePions->SetParName(3,"Triangular flow");
\r
144 fAzimuthalAnglePions->SetParameter(3,fTriangularFlowPions);
\r
145 fAzimuthalAnglePions->SetParName(4,"Quandrangular flow");
\r
146 fAzimuthalAnglePions->SetParameter(4,fQuandrangularFlowPions);
\r
147 fAzimuthalAnglePions->SetParName(5,"Pentangular flow");
\r
148 fAzimuthalAnglePions->SetParameter(5,fPentangularFlowPions);
\r
150 fAzimuthalAngleKaons = new TF1("fAzimuthalAngleKaons","1+2.*[1]*TMath::Cos(x-[0])+2.*[2]*TMath::Cos(2*(x-[0]))+2.*[3]*TMath::Cos(3*(x-[0]))+2.*[4]*TMath::Cos(4*(x-[0]))+2.*[5]*TMath::Cos(5*(x-[0]))",0.,2.*TMath::Pi());
\r
151 fAzimuthalAngleKaons->SetParName(0,"Reaction Plane");
\r
152 fAzimuthalAngleKaons->SetParameter(0,fReactionPlane);
\r
153 fAzimuthalAngleKaons->SetParName(1,"Directed flow");
\r
154 fAzimuthalAngleKaons->SetParameter(1,fDirectedFlowKaons);
\r
155 fAzimuthalAngleKaons->SetParName(2,"Elliptic flow");
\r
156 fAzimuthalAngleKaons->SetParameter(2,fEllipticFlowKaons);
\r
157 fAzimuthalAngleKaons->SetParName(3,"Triangular flow");
\r
158 fAzimuthalAngleKaons->SetParameter(3,fTriangularFlowKaons);
\r
159 fAzimuthalAngleKaons->SetParName(4,"Quandrangular flow");
\r
160 fAzimuthalAngleKaons->SetParameter(4,fQuandrangularFlowKaons);
\r
161 fAzimuthalAngleKaons->SetParName(5,"Pentangular flow");
\r
162 fAzimuthalAngleKaons->SetParameter(5,fPentangularFlowKaons);
\r
164 fAzimuthalAngleProtons = new TF1("fAzimuthalAngleProtons","1+2.*[1]*TMath::Cos(x-[0])+2.*[2]*TMath::Cos(2*(x-[0]))+2.*[3]*TMath::Cos(3*(x-[0]))+2.*[4]*TMath::Cos(4*(x-[0]))+2.*[5]*TMath::Cos(5*(x-[0]))",0.,2.*TMath::Pi());
\r
165 fAzimuthalAngleProtons->SetParName(0,"Reaction Plane");
\r
166 fAzimuthalAngleProtons->SetParameter(0,fReactionPlane);
\r
167 fAzimuthalAngleProtons->SetParName(1,"Directed flow");
\r
168 fAzimuthalAngleProtons->SetParameter(1,fDirectedFlowProtons);
\r
169 fAzimuthalAngleProtons->SetParName(2,"Elliptic flow");
\r
170 fAzimuthalAngleProtons->SetParameter(2,fEllipticFlowProtons);
\r
171 fAzimuthalAngleProtons->SetParName(3,"Triangular flow");
\r
172 fAzimuthalAngleProtons->SetParameter(3,fTriangularFlowProtons);
\r
173 fAzimuthalAngleProtons->SetParName(4,"Quandrangular flow");
\r
174 fAzimuthalAngleProtons->SetParameter(4,fQuandrangularFlowProtons);
\r
175 fAzimuthalAngleProtons->SetParName(5,"Pentangular flow");
\r
176 fAzimuthalAngleProtons->SetParameter(5,fPentangularFlowProtons);
\r
177 //==============Flow values==============//
\r
180 //________________________________________________________________________
\r
181 void AliAnalysisTaskToyModel::UserCreateOutputObjects() {
\r
182 // Create histograms
\r
185 fBalance = new AliBalance();
\r
186 fBalance->SetInterval(-1,-0.8,0.8,16,0.,1.6);
\r
188 if(fRunShuffling) {
\r
189 if(!fShuffledBalance) {
\r
190 fShuffledBalance = new AliBalance();
\r
191 fShuffledBalance->SetInterval(-1,-0.8,0.8,16,0.,1.6);
\r
196 fList = new TList();
\r
197 fList->SetName("listQA");
\r
200 //Balance Function list
\r
201 fListBF = new TList();
\r
202 fListBF->SetName("listBF");
\r
203 fListBF->SetOwner();
\r
205 if(fRunShuffling) {
\r
206 fListBFS = new TList();
\r
207 fListBFS->SetName("listBFShuffled");
\r
208 fListBFS->SetOwner();
\r
212 TString gCutName[4] = {"Total","Offline trigger",
\r
213 "Vertex","Analyzed"};
\r
214 fHistEventStats = new TH1F("fHistEventStats",
\r
215 "Event statistics;;N_{events}",
\r
217 for(Int_t i = 1; i <= 4; i++)
\r
218 fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data());
\r
219 fList->Add(fHistEventStats);
\r
221 // Balance function histograms
\r
222 // Initialize histograms if not done yet
\r
223 if(!fBalance->GetHistNp(0)){
\r
224 AliWarning("Histograms not yet initialized! --> Will be done now");
\r
225 AliWarning("--> Add 'gBalance->InitHistograms()' in your configBalanceFunction");
\r
226 fBalance->InitHistograms();
\r
229 if(fRunShuffling) {
\r
230 if(!fShuffledBalance->GetHistNp(0)) {
\r
231 AliWarning("Histograms (shuffling) not yet initialized! --> Will be done now");
\r
232 AliWarning("--> Add 'gBalance->InitHistograms()' in your configBalanceFunction");
\r
233 fShuffledBalance->InitHistograms();
\r
237 for(Int_t a = 0; a < ANALYSIS_TYPES; a++){
\r
238 fListBF->Add(fBalance->GetHistNp(a));
\r
239 fListBF->Add(fBalance->GetHistNn(a));
\r
240 fListBF->Add(fBalance->GetHistNpn(a));
\r
241 fListBF->Add(fBalance->GetHistNnn(a));
\r
242 fListBF->Add(fBalance->GetHistNpp(a));
\r
243 fListBF->Add(fBalance->GetHistNnp(a));
\r
245 if(fRunShuffling) {
\r
246 fListBFS->Add(fShuffledBalance->GetHistNp(a));
\r
247 fListBFS->Add(fShuffledBalance->GetHistNn(a));
\r
248 fListBFS->Add(fShuffledBalance->GetHistNpn(a));
\r
249 fListBFS->Add(fShuffledBalance->GetHistNnn(a));
\r
250 fListBFS->Add(fShuffledBalance->GetHistNpp(a));
\r
251 fListBFS->Add(fShuffledBalance->GetHistNnp(a));
\r
255 // Post output data.
\r
256 PostData(1, fList);
\r
257 PostData(2, fListBF);
\r
258 if(fRunShuffling) PostData(3, fListBFS);
\r
261 //________________________________________________________________________
\r
262 void AliAnalysisTaskToyModel::UserExec(Option_t *) {
\r
264 // Called for each event
\r
266 // vector holding the charges/kinematics of all tracks (charge,y,eta,phi,p0,p1,p2,pt,E)
\r
267 vector<Double_t> *chargeVectorShuffle[9]; // this will be shuffled
\r
268 vector<Double_t> *chargeVector[9]; // original charge
\r
269 for(Int_t i = 0; i < 9; i++){
\r
270 chargeVectorShuffle[i] = new vector<Double_t>;
\r
271 chargeVector[i] = new vector<Double_t>;
\r
283 Int_t nMultiplicity = (Int_t)(gRandom->Gaus(fTotalMultiplicityMean,fTotalMultiplicitySigma));
\r
284 Int_t nNetCharge = (Int_t)(gRandom->Gaus(fNetChargeMean,fNetChargeSigma));
\r
285 Int_t nGeneratedPositive = 0.5*(nMultiplicity + nNetCharge);
\r
286 Int_t nGeneratedNegative = nMultiplicity - nGeneratedPositive;
\r
287 Int_t nGeneratedPositivePions = (Int_t)(fPionPercentage*nGeneratedPositive);
\r
288 Int_t nGeneratedNegativePions = (Int_t)(fPionPercentage*nGeneratedNegative);
\r
289 Int_t nGeneratedPositiveKaons = (Int_t)(fKaonPercentage*nGeneratedPositive);
\r
290 Int_t nGeneratedNegativeKaons = (Int_t)(fKaonPercentage*nGeneratedNegative);
\r
291 Int_t nGeneratedPositiveProtons = (Int_t)(fProtonPercentage*nGeneratedPositive);
\r
292 Int_t nGeneratedNegativeProtons = (Int_t)(fProtonPercentage*nGeneratedNegative);
\r
294 Printf("Total positive: %d - Total negative: %d",nGeneratedPositive,nGeneratedNegative);
\r
295 Printf("Positive pions: %d - Negative pions: %d",nGeneratedPositivePions,nGeneratedNegativePions);
\r
296 Printf("Positive kaons: %d - Negative kaons: %d",nGeneratedPositiveKaons,nGeneratedNegativeKaons);
\r
297 Printf("Positive protons: %d - Negative protons: %d",nGeneratedPositiveProtons,nGeneratedNegativeProtons);
\r
299 Int_t gNumberOfAcceptedParticles = 0;
\r
300 //positive particles
\r
301 for(Int_t iPosCount = 0; iPosCount < nGeneratedPositive; iPosCount++) {
\r
304 v_pt = fPtSpectraAllCharges->GetRandom();
\r
305 v_phi = fAzimuthalAngleAllCharges->GetRandom();
\r
306 v_eta = gRandom->Gaus(0.0,4.0);
\r
308 v_p[0] = v_pt*TMath::Cos(v_phi);
\r
309 v_p[1] = v_pt*TMath::Sin(v_phi);
\r
310 v_p[2] = v_pt*TMath::SinH(v_eta);
\r
311 v_E = TMath::Sqrt(TMath::Power(0.139,2) +
\r
312 TMath::Power(v_p[0],2) +
\r
313 TMath::Power(v_p[1],2) +
\r
314 TMath::Power(v_p[2],2));
\r
316 v_y = 0.5*TMath::Log((v_E + v_p[2])/(v_E - v_p[2]));
\r
319 if((v_eta < fEtaMin) || (v_eta > fEtaMax)) continue;
\r
320 if((v_pt < fPtMin) || (v_pt > fPtMax)) continue;
\r
322 // fill charge vector
\r
323 chargeVector[0]->push_back(v_charge);
\r
324 chargeVector[1]->push_back(v_y);
\r
325 chargeVector[2]->push_back(v_eta);
\r
326 chargeVector[3]->push_back(v_phi);
\r
327 chargeVector[4]->push_back(v_p[0]);
\r
328 chargeVector[5]->push_back(v_p[1]);
\r
329 chargeVector[6]->push_back(v_p[2]);
\r
330 chargeVector[7]->push_back(v_pt);
\r
331 chargeVector[8]->push_back(v_E);
\r
333 if(fRunShuffling) {
\r
334 chargeVectorShuffle[0]->push_back(v_charge);
\r
335 chargeVectorShuffle[1]->push_back(v_y);
\r
336 chargeVectorShuffle[2]->push_back(v_eta);
\r
337 chargeVectorShuffle[3]->push_back(v_phi);
\r
338 chargeVectorShuffle[4]->push_back(v_p[0]);
\r
339 chargeVectorShuffle[5]->push_back(v_p[1]);
\r
340 chargeVectorShuffle[6]->push_back(v_p[2]);
\r
341 chargeVectorShuffle[7]->push_back(v_pt);
\r
342 chargeVectorShuffle[8]->push_back(v_E);
\r
344 gNumberOfAcceptedParticles += 1;
\r
346 }//positive particle loop
\r
348 //negative particles
\r
349 for(Int_t iNegCount = 0; iNegCount < nGeneratedNegative; iNegCount++) {
\r
352 v_pt = fPtSpectraAllCharges->GetRandom();
\r
353 v_phi = fAzimuthalAngleAllCharges->GetRandom();
\r
354 v_eta = gRandom->Gaus(0.0,4.0);
\r
356 v_p[0] = v_pt*TMath::Cos(v_phi);
\r
357 v_p[1] = v_pt*TMath::Sin(v_phi);
\r
358 v_p[2] = v_pt*TMath::SinH(v_eta);
\r
359 v_E = TMath::Sqrt(TMath::Power(0.139,2) +
\r
360 TMath::Power(v_p[0],2) +
\r
361 TMath::Power(v_p[1],2) +
\r
362 TMath::Power(v_p[2],2));
\r
364 v_y = 0.5*TMath::Log((v_E + v_p[2])/(v_E - v_p[2]));
\r
367 if((v_eta < fEtaMin) || (v_eta > fEtaMax)) continue;
\r
368 if((v_pt < fPtMin) || (v_pt > fPtMax)) continue;
\r
370 // fill charge vector
\r
371 chargeVector[0]->push_back(v_charge);
\r
372 chargeVector[1]->push_back(v_y);
\r
373 chargeVector[2]->push_back(v_eta);
\r
374 chargeVector[3]->push_back(v_phi);
\r
375 chargeVector[4]->push_back(v_p[0]);
\r
376 chargeVector[5]->push_back(v_p[1]);
\r
377 chargeVector[6]->push_back(v_p[2]);
\r
378 chargeVector[7]->push_back(v_pt);
\r
379 chargeVector[8]->push_back(v_E);
\r
381 if(fRunShuffling) {
\r
382 chargeVectorShuffle[0]->push_back(v_charge);
\r
383 chargeVectorShuffle[1]->push_back(v_y);
\r
384 chargeVectorShuffle[2]->push_back(v_eta);
\r
385 chargeVectorShuffle[3]->push_back(v_phi);
\r
386 chargeVectorShuffle[4]->push_back(v_p[0]);
\r
387 chargeVectorShuffle[5]->push_back(v_p[1]);
\r
388 chargeVectorShuffle[6]->push_back(v_p[2]);
\r
389 chargeVectorShuffle[7]->push_back(v_pt);
\r
390 chargeVectorShuffle[8]->push_back(v_E);
\r
392 gNumberOfAcceptedParticles += 1;
\r
394 }//negative particle loop
\r
396 fBalance->CalculateBalance(gNumberOfAcceptedParticles,chargeVector);
\r
398 if(fRunShuffling) {
\r
400 random_shuffle( chargeVectorShuffle[0]->begin(), chargeVectorShuffle[0]->end() );
\r
401 fShuffledBalance->CalculateBalance(gNumberOfAcceptedParticles,chargeVectorShuffle);
\r
405 //________________________________________________________________________
\r
406 void AliAnalysisTaskToyModel::FinishTaskOutput() {
\r
407 //Printf("END BF");
\r
410 Printf("ERROR: fBalance not available");
\r
413 if(fRunShuffling) {
\r
414 if (!fShuffledBalance) {
\r
415 Printf("ERROR: fShuffledBalance not available");
\r
421 //________________________________________________________________________
\r
422 void AliAnalysisTaskToyModel::Terminate(Option_t *) {
\r
423 // Draw result to the screen
\r
424 // Called once at the end of the query
\r
426 // not implemented ...
\r