]>
Commit | Line | Data |
---|---|---|
ee97f33d | 1 | //========================================================// |
2 | //Macro that reads the output of the hydro calculations | |
3 | //(Therminator) using the format that can be found at | |
4 | //arXiv:1102.0273 | |
5 | //Author: Panos.Christakoglou@nikhef.nl | |
6 | //========================================================// | |
7 | ||
8 | //========================================================// | |
9 | //Event structure | |
10 | struct StructEvent { | |
11 | UInt_t eventID;//unique identifier of the event (random number) | |
12 | UInt_t entries;//number of particles of each event | |
13 | UInt_t entriesprev;//set to 0 by default | |
14 | }; | |
15 | //========================================================// | |
16 | ||
17 | //========================================================// | |
18 | //Particle structure | |
a7f9d873 | 19 | //Primordial particles: fathereid==-1 ==> |
20 | //pid==fatherpid==rootpid | |
ee97f33d | 21 | struct StructParticle { |
22 | Float_t mass;//the mass of the particle (in GeV) | |
23 | Float_t t;//the time coordinate of the particle in fm/c | |
24 | Float_t x;//the spacial coordinate x in fm/c | |
25 | Float_t y;//the spacial coordinate y in fm/c | |
26 | Float_t z;//the spacial coordinate z in fm/c | |
27 | Float_t e;//the energy in GeV | |
28 | Float_t px;//the x-coordinate of the particle's momentum in GeV | |
29 | Float_t py;//the y-coordinate of the particle's momentum in GeV | |
30 | Float_t pz;//the z-coordinate of the particle's momentum in GeV | |
31 | Int_t decayed;//decay flag (1: decayed, 0: no decay) | |
32 | Int_t pid;//PDG of particle | |
33 | Int_t fatherpid;//PDG of parent | |
34 | Int_t rootpid;//root (primordial) particle PDG number | |
35 | Int_t eid;//sequence number in the event | |
36 | Int_t fathereid;//parent sequence number in the event | |
37 | UInt_t eventid;//unique identifier of the event (random number) | |
38 | }; | |
39 | //========================================================// | |
40 | ||
41 | //========================================================// | |
42 | //Balance function analysis variables | |
43 | Bool_t gRunShuffling=kFALSE; | |
eaae9cd6 | 44 | Bool_t gRunMixing=kTRUE; |
ee97f33d | 45 | Bool_t gRunMixingWithEventPlane=kFALSE; |
46 | ||
47 | Double_t gEtaMin = -0.8; | |
48 | Double_t gEtaMax = 0.8; | |
49 | Double_t gPtMin = 0.2; | |
50 | Double_t gPtMax = 20.0; | |
51 | //========================================================// | |
52 | ||
53 | ||
8e81f4e0 | 54 | void runBalanceFunctionOnHydro(TString aEventDir = ".", |
55 | Int_t aEventFiles = 9) { | |
ee97f33d | 56 | //Macro that reads the themrinator events |
57 | //Author: Panos.Christakoglou@nikhef.nl | |
58 | // Time: | |
59 | TStopwatch timer; | |
60 | timer.Start(); | |
61 | ||
62 | //========================================================// | |
63 | //Load the aliroot libraries | |
64 | gSystem->Load("libSTEERBase.so"); | |
65 | gSystem->Load("libESD.so"); | |
66 | gSystem->Load("libAOD.so"); | |
67 | gSystem->Load("libANALYSIS.so"); | |
68 | gSystem->Load("libANALYSISalice.so"); | |
69 | gSystem->Load("libEventMixing.so"); | |
70 | gSystem->Load("libCORRFW.so"); | |
71 | gSystem->Load("libPWGTools.so"); | |
72 | gSystem->Load("libPWGCFebye.so"); | |
73 | //========================================================// | |
74 | ||
75 | //========================================================// | |
76 | //Configure the bf objects | |
77 | AliBalancePsi *bf = new AliBalancePsi(); | |
78 | bf->SetAnalysisLevel("MC"); | |
79 | bf->SetShuffle(gRunShuffling); | |
80 | bf->SetEventClass("EventPlane"); | |
81 | bf->SetDeltaEtaMax(TMath::Abs(gEtaMax-gEtaMin)); | |
82 | bf->InitHistograms(); | |
83 | ||
84 | AliBalancePsi *bfs = 0x0; | |
85 | if(gRunShuffling) { | |
86 | bfs = new AliBalancePsi(); | |
87 | bfs->SetAnalysisLevel("MC"); | |
88 | bfs->SetEventClass("EventPlane"); | |
89 | bfs->SetDeltaEtaMax(TMath::Abs(gEtaMax-gEtaMin)); | |
90 | bfs->InitHistograms(); | |
91 | } | |
92 | ||
93 | AliBalancePsi *bfm = 0x0; | |
94 | if(gRunMixing) { | |
95 | bfm = new AliBalancePsi(); | |
96 | bfm->SetAnalysisLevel("MC"); | |
97 | bfm->SetShuffle(gRunShuffling); | |
98 | bfm->SetEventClass("EventPlane"); | |
99 | bfm->SetDeltaEtaMax(TMath::Abs(gEtaMax-gEtaMin)); | |
100 | bfm->InitHistograms(); | |
101 | } | |
102 | //========================================================// | |
103 | ||
104 | //========================================================// | |
105 | //Output objects | |
106 | //QA list | |
107 | fList = new TList(); | |
108 | fList->SetName("listQA"); | |
109 | fList->SetOwner(); | |
110 | ||
111 | //Balance Function list | |
112 | TList *fListBF = new TList(); | |
113 | fListBF->SetName("listBF"); | |
114 | fListBF->SetOwner(); | |
eaae9cd6 | 115 | fListBF->Add(bf->GetHistNp()); |
116 | fListBF->Add(bf->GetHistNn()); | |
117 | fListBF->Add(bf->GetHistNpn()); | |
118 | fListBF->Add(bf->GetHistNnn()); | |
119 | fListBF->Add(bf->GetHistNpp()); | |
120 | fListBF->Add(bf->GetHistNnp()); | |
121 | ||
ee97f33d | 122 | //Balance function list: shuffling |
123 | TList *fListBFS = 0x0; | |
124 | if(gRunShuffling) { | |
125 | fListBFS = new TList(); | |
126 | fListBFS->SetName("listBFShuffled"); | |
127 | fListBFS->SetOwner(); | |
128 | fListBFS->Add(bfs->GetHistNp()); | |
129 | fListBFS->Add(bfs->GetHistNn()); | |
130 | fListBFS->Add(bfs->GetHistNpn()); | |
131 | fListBFS->Add(bfs->GetHistNnn()); | |
132 | fListBFS->Add(bfs->GetHistNpp()); | |
133 | fListBFS->Add(bfs->GetHistNnp()); | |
134 | } | |
135 | ||
136 | //Balance function list: event mixing | |
137 | TList *fListBFM = 0x0; | |
138 | if(gRunMixing) { | |
139 | fListBFM = new TList(); | |
140 | fListBFM->SetName("listBFMixed"); | |
141 | fListBFM->SetOwner(); | |
142 | fListBFM->Add(bfm->GetHistNp()); | |
143 | fListBFM->Add(bfm->GetHistNn()); | |
144 | fListBFM->Add(bfm->GetHistNpn()); | |
145 | fListBFM->Add(bfm->GetHistNnn()); | |
146 | fListBFM->Add(bfm->GetHistNpp()); | |
147 | fListBFM->Add(bfm->GetHistNnp()); | |
148 | } | |
149 | //========================================================// | |
150 | ||
151 | //========================================================// | |
152 | //Event Mixing | |
153 | if(gRunMixing){ | |
154 | Int_t trackDepth = 50000; | |
155 | Int_t poolsize = 1000; // Maximum number of events, ignored in the present implemented of AliEventPoolManager | |
eaae9cd6 | 156 | |
157 | // centrality bins | |
158 | Double_t centralityBins[] = {0.,100.}; | |
159 | Double_t* centbins = centralityBins; | |
160 | Int_t nCentralityBins = sizeof(centralityBins) / sizeof(Double_t) - 1; | |
161 | ||
162 | // Zvtx bins | |
163 | Double_t vertexBins[] = {-10., 10.}; | |
164 | Double_t* vtxbins = vertexBins; | |
165 | Int_t nVertexBins = sizeof(vertexBins) / sizeof(Double_t) - 1; | |
ee97f33d | 166 | |
ee97f33d | 167 | AliEventPoolManager *fPoolMgr = 0x0; |
eaae9cd6 | 168 | fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centbins, nVertexBins, vtxbins); |
ee97f33d | 169 | } |
170 | //========================================================// | |
171 | ||
172 | //========================================================// | |
173 | //Create the TChain object: events | |
174 | TChain *eventChain = new TChain("events"); | |
175 | StructEvent tStructEvents; | |
176 | eventChain->SetBranchAddress("event",&tStructEvents); | |
177 | //========================================================// | |
178 | ||
179 | //========================================================// | |
180 | //Create the TChain object: particles | |
181 | TChain *particleChain = new TChain("particles"); | |
182 | StructParticle tStructParticles; | |
183 | particleChain->SetBranchAddress("particle",&tStructParticles); | |
184 | //========================================================// | |
185 | ||
186 | //========================================================// | |
187 | //Fill the TChain with the files in the directory | |
8e81f4e0 | 188 | for(Int_t iFile = 1; iFile <= aEventFiles; iFile++) { |
ee97f33d | 189 | TString filename = aEventDir.Data(); |
190 | filename += "/event00"; filename += iFile; | |
191 | filename += ".root"; | |
192 | cout<<"Adding file "<<filename.Data()<<" to the chain..."<<endl; | |
193 | eventChain->Add(filename.Data()); | |
194 | particleChain->Add(filename.Data()); | |
195 | } | |
196 | //========================================================// | |
197 | ||
198 | //========================================================// | |
eaae9cd6 | 199 | //QA histograms |
ee97f33d | 200 | //Event stats. |
201 | TString gCutName[5] = {"Total","Offline trigger", | |
202 | "Vertex","Analyzed","sel. Centrality"}; | |
a7f9d873 | 203 | TH1F *fHistEventStats = new TH1F("fHistEventStats", |
204 | "Event statistics;;N_{events}", | |
205 | 5,0.5,5.5); | |
ee97f33d | 206 | for(Int_t i = 1; i <= 5; i++) |
207 | fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data()); | |
208 | fList->Add(fHistEventStats); | |
209 | ||
210 | //Number of accepted particles | |
a7f9d873 | 211 | TH1F *fHistNumberOfAcceptedTracks = new TH1F("fHistNumberOfAcceptedTracks",";N_{acc.};Entries",4001,-0.5,4000.5); |
ee97f33d | 212 | fList->Add(fHistNumberOfAcceptedTracks); |
eaae9cd6 | 213 | |
214 | //Particle level: eta-phi | |
215 | TH2F *fHistEtaPhiPositive = new TH2F("fHistEtaPhiPositive",";#eta;#varphi (rad)",80,gEtaMin,gEtaMax,72,0.0,2.*TMath::Pi()); | |
216 | fList->Add(fHistEtaPhiPositive); | |
217 | TH2F *fHistEtaPhiNegative = new TH2F("fHistEtaPhiNegative",";#eta;#varphi (rad)",80,gEtaMin,gEtaMax,72,0.0,2.*TMath::Pi()); | |
218 | fList->Add(fHistEtaPhiNegative); | |
219 | ||
220 | //Particle level: pt | |
221 | TH1F *fHistPtPositive = new TH1F("fHistPtPositive",";p_{T} (GeV/c)",100,0.01,30.01); | |
222 | fList->Add(fHistPtPositive); | |
223 | TH1F *fHistPtNegative = new TH1F("fHistPtNegative",";p_{T} (GeV/c)",100,0.01,30.01); | |
224 | fList->Add(fHistPtNegative); | |
ee97f33d | 225 | //========================================================// |
226 | ||
227 | //========================================================// | |
228 | //loop over the events | |
229 | Int_t nEvents = eventChain->GetEntries(); | |
230 | cout<<"========================================="<<endl; | |
231 | cout<<"Number of events in the chain: "<<nEvents<<endl; | |
232 | cout<<"========================================="<<endl; | |
233 | Int_t iParticleCounter = 0; | |
234 | Int_t nTotalParticles = 0; | |
235 | ||
eaae9cd6 | 236 | for(Int_t iEvent = 0; iEvent < nEvents; iEvent++) { |
237 | //for(Int_t iEvent = 0; iEvent < 1; iEvent++) { | |
ee97f33d | 238 | eventChain->GetEntry(iEvent); |
239 | ||
a7f9d873 | 240 | //========================================================// |
241 | //Create the TObjArray object to store the particles | |
242 | TObjArray* tracksAccepted = new TObjArray; | |
243 | tracksAccepted->SetOwner(kTRUE); | |
244 | //========================================================// | |
245 | ||
ee97f33d | 246 | Int_t nParticles = tStructEvents.entries; |
a7f9d873 | 247 | if(((iEvent+1)%100)==0) |
248 | cout<<"Event: "<<iEvent+1<<" - ID: "<<tStructEvents.eventID<<" - Entries: "<<tStructEvents.entries<<endl; | |
ee97f33d | 249 | |
a7f9d873 | 250 | //========================================================// |
251 | //Filling event level histos | |
252 | fHistEventStats->Fill(1); | |
253 | fHistEventStats->Fill(2); | |
254 | fHistEventStats->Fill(3); | |
255 | fHistEventStats->Fill(4); | |
256 | fHistEventStats->Fill(5); | |
257 | ||
258 | Int_t gNumberOfAcceptedParticles = 0; | |
259 | Double_t gReactionPlane = 0.; | |
eaae9cd6 | 260 | Double_t gCharge = 0.; |
ee97f33d | 261 | //========================================================// |
262 | //loop over particles | |
263 | for(Int_t iParticle = 0; iParticle < nParticles; iParticle++) { | |
264 | particleChain->GetEntry(nTotalParticles+iParticle); | |
265 | ||
a7f9d873 | 266 | //========================================================// |
267 | //consider only primordial particles | |
eaae9cd6 | 268 | if(!IsPhysicalPrimary(tStructParticles.pid, |
269 | tStructParticles.fathereid, | |
270 | gCharge)) continue; | |
a7f9d873 | 271 | |
272 | //========================================================// | |
273 | //Calculate kinematic variables | |
274 | Double_t gPt = TMath::Sqrt(TMath::Power(tStructParticles.px,2) + | |
275 | TMath::Power(tStructParticles.py,2)); | |
276 | Double_t gP = TMath::Sqrt(TMath::Power(tStructParticles.px,2) + | |
277 | TMath::Power(tStructParticles.py,2) + | |
278 | TMath::Power(tStructParticles.pz,2) ); | |
279 | Double_t gEta = -100.; | |
280 | if(gP != tStructParticles.pz) | |
281 | gEta = 0.5*TMath::Log((gP + tStructParticles.pz)/(gP - tStructParticles.pz)); | |
eaae9cd6 | 282 | Double_t gPhi = TMath::Pi()+TMath::ATan2(-tStructParticles.py,-tStructParticles.px); |
283 | ||
284 | //========================================================// | |
285 | //Fill QA | |
286 | if(gCharge > 0) { | |
287 | fHistEtaPhiPositive->Fill(gEta,gPhi); | |
288 | fHistPtPositive->Fill(gPt); | |
289 | } | |
290 | else if(gCharge < 0) { | |
291 | fHistEtaPhiNegative->Fill(gEta,gPhi); | |
292 | fHistPtNegative->Fill(gPt); | |
293 | } | |
294 | //========================================================// | |
295 | ||
a7f9d873 | 296 | //========================================================// |
297 | //Apply cuts | |
298 | if((gEta > gEtaMax)||(gEta < gEtaMin)) continue; | |
299 | if((gPt > gPtMax)||(gPt < gPtMin)) continue; | |
300 | ||
301 | tracksAccepted->Add(new AliBFBasicParticle(gEta,gPhi,gPt,gCharge, 1.)); | |
302 | gNumberOfAcceptedParticles += 1; | |
303 | ||
ee97f33d | 304 | iParticleCounter += 1; |
a7f9d873 | 305 | //cout<<"\t Particle counter: "<<iParticleCounter<<" - Particle in the event: "<<iParticle+1<<" - eventID: "<<tStructParticles.eventid<<" - pid: "<<tStructParticles.pid<<" - fatherpid: "<<tStructParticles.fatherpid<<" - rootpid: "<<tStructParticles.rootpid<<" - fathereid: "<<tStructParticles.fathereid<<" - eid: "<<tStructParticles.eid<<endl; |
ee97f33d | 306 | }//particle loop |
a7f9d873 | 307 | |
308 | //========================================================// | |
309 | // Event mixing (borrowed code from the task) | |
eaae9cd6 | 310 | if (gRunMixing) { |
311 | Int_t fMixingTracks = 50000; | |
312 | AliEventPool* pool = fPoolMgr->GetEventPool(1.,0.,0.); | |
a7f9d873 | 313 | |
314 | if (!pool) { | |
eaae9cd6 | 315 | AliFatal(Form("No pool found")); |
a7f9d873 | 316 | } |
317 | else { | |
eaae9cd6 | 318 | //pool->SetDebug(1); |
319 | if (pool->IsReady() || pool->NTracksInPool() > fMixingTracks / 10 || pool->GetCurrentNEvents() >= 5){ | |
320 | ||
321 | Int_t nMix = pool->GetCurrentNEvents(); | |
322 | //cout << "nMix = " << nMix << " tracks in pool = " << pool->NTracksInPool() << endl; | |
323 | ||
324 | // Fill mixed-event histos here | |
325 | for (Int_t jMix=0; jMix<nMix; jMix++) { | |
326 | TObjArray* tracksMixed = pool->GetEvent(jMix); | |
327 | bfm->CalculateBalance(gReactionPlane,tracksAccepted,tracksMixed,1.,0.); | |
328 | } | |
329 | } | |
a7f9d873 | 330 | |
eaae9cd6 | 331 | // Update the Event pool |
332 | pool->UpdatePool(tracksAccepted); | |
333 | //pool->PrintInfo(); | |
334 | ||
a7f9d873 | 335 | }//pool NULL check |
eaae9cd6 | 336 | }//run mixing |
a7f9d873 | 337 | //========================================================// |
338 | ||
339 | //========================================================// | |
340 | // calculate balance function | |
eaae9cd6 | 341 | bf->CalculateBalance(gReactionPlane,tracksAccepted,NULL,1.,0.); |
a7f9d873 | 342 | |
343 | fHistNumberOfAcceptedTracks->Fill(gNumberOfAcceptedParticles); | |
ee97f33d | 344 | nTotalParticles += nParticles; |
345 | } | |
346 | ||
347 | //========================================================// | |
348 | //Output file | |
a7f9d873 | 349 | TFile *f = TFile::Open("AnalysisResults.root","recreate"); |
8e81f4e0 | 350 | TDirectoryFile *dir = new TDirectoryFile("PWGCFEbyE.outputBalanceFunctionPsiAnalysis","PWGCFEbyE.outputBalanceFunctionPsiAnalysis"); |
351 | ||
352 | fList->SetName("listQA"); | |
353 | fList->SetOwner(kTRUE); | |
354 | dir->Add(fList); | |
355 | ||
356 | fListBF->SetName("listBF"); | |
357 | fListBF->SetOwner(kTRUE); | |
358 | dir->Add(fListBF); | |
359 | ||
360 | if(gRunMixing) { | |
361 | fListBFM->SetName("listBFMixed"); | |
362 | fListBFM->SetOwner(kTRUE); | |
363 | dir->Add(fListBFM); | |
364 | } | |
365 | ||
366 | dir->Write(dir->GetName(),TObject::kSingleKey); | |
ee97f33d | 367 | f->Close(); |
368 | //========================================================// | |
369 | ||
370 | // Print real and CPU time used for analysis: | |
371 | timer.Stop(); | |
372 | timer.Print(); | |
373 | } | |
eaae9cd6 | 374 | |
375 | //______________________________________________________________// | |
376 | Bool_t IsPhysicalPrimary(Int_t gPDGCode, | |
377 | Int_t gFathereid, | |
378 | Double_t &gCharge) { | |
379 | //Check whether the primordial particle belongs to the list | |
380 | //of known particles | |
381 | Bool_t kStatus = kFALSE; | |
382 | if(gFathereid != -1) | |
383 | return kStatus; | |
384 | ||
385 | //List of stable particles | |
386 | const Int_t kNstable = 3; | |
387 | Int_t pdgStable[kNstable] = { | |
388 | 211, // Pion | |
389 | 321, // Kaon | |
390 | 2212, // Proton | |
391 | }; | |
392 | ||
393 | if(gPDGCode < 0) gCharge = -1; | |
394 | else if(gPDGCode > 0) gCharge = 1; | |
395 | ||
396 | for(Int_t iParticle = 0; iParticle < kNstable; iParticle++) { | |
397 | if((TMath::Abs(gPDGCode) == pdgStable[iParticle])) { | |
398 | kStatus = kTRUE; | |
399 | return kStatus; | |
400 | } | |
401 | } | |
402 | ||
403 | return kFALSE; | |
404 | } |