1 //========================================================//
2 //Macro that reads the output of the hydro calculations
3 //(Therminator) using the format that can be found at
5 //Author: Panos.Christakoglou@nikhef.nl
6 //========================================================//
8 //========================================================//
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
15 //========================================================//
17 //========================================================//
19 //Primordial particles: fathereid==-1 ==>
20 //pid==fatherpid==rootpid
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)
39 //========================================================//
41 //========================================================//
42 //Balance function analysis variables
43 Bool_t gRunShuffling=kFALSE;
44 Bool_t gRunMixing=kTRUE;
45 Bool_t gRunMixingWithEventPlane=kFALSE;
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 //========================================================//
54 void runBalanceFunctionOnHydro(TString aEventDir = ".",
55 Int_t aEventFiles = 9) {
56 //Macro that reads the themrinator events
57 //Author: Panos.Christakoglou@nikhef.nl
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 //========================================================//
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));
84 AliBalancePsi *bfs = 0x0;
86 bfs = new AliBalancePsi();
87 bfs->SetAnalysisLevel("MC");
88 bfs->SetEventClass("EventPlane");
89 bfs->SetDeltaEtaMax(TMath::Abs(gEtaMax-gEtaMin));
90 bfs->InitHistograms();
93 AliBalancePsi *bfm = 0x0;
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();
102 //========================================================//
104 //========================================================//
108 fList->SetName("listQA");
111 //Balance Function list
112 TList *fListBF = new TList();
113 fListBF->SetName("listBF");
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());
122 //Balance function list: shuffling
123 TList *fListBFS = 0x0;
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());
136 //Balance function list: event mixing
137 TList *fListBFM = 0x0;
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());
149 //========================================================//
151 //========================================================//
154 Int_t trackDepth = 50000;
155 Int_t poolsize = 1000; // Maximum number of events, ignored in the present implemented of AliEventPoolManager
158 Double_t centralityBins[] = {0.,100.};
159 Double_t* centbins = centralityBins;
160 Int_t nCentralityBins = sizeof(centralityBins) / sizeof(Double_t) - 1;
163 Double_t vertexBins[] = {-10., 10.};
164 Double_t* vtxbins = vertexBins;
165 Int_t nVertexBins = sizeof(vertexBins) / sizeof(Double_t) - 1;
167 AliEventPoolManager *fPoolMgr = 0x0;
168 fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centbins, nVertexBins, vtxbins);
170 //========================================================//
172 //========================================================//
173 //Create the TChain object: events
174 TChain *eventChain = new TChain("events");
175 StructEvent tStructEvents;
176 eventChain->SetBranchAddress("event",&tStructEvents);
177 //========================================================//
179 //========================================================//
180 //Create the TChain object: particles
181 TChain *particleChain = new TChain("particles");
182 StructParticle tStructParticles;
183 particleChain->SetBranchAddress("particle",&tStructParticles);
184 //========================================================//
186 //========================================================//
187 //Fill the TChain with the files in the directory
188 for(Int_t iFile = 1; iFile <= aEventFiles; iFile++) {
189 TString filename = aEventDir.Data();
190 filename += "/event00"; filename += iFile;
192 cout<<"Adding file "<<filename.Data()<<" to the chain..."<<endl;
193 eventChain->Add(filename.Data());
194 particleChain->Add(filename.Data());
196 //========================================================//
198 //========================================================//
201 TString gCutName[5] = {"Total","Offline trigger",
202 "Vertex","Analyzed","sel. Centrality"};
203 TH1F *fHistEventStats = new TH1F("fHistEventStats",
204 "Event statistics;;N_{events}",
206 for(Int_t i = 1; i <= 5; i++)
207 fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data());
208 fList->Add(fHistEventStats);
210 //Number of accepted particles
211 TH1F *fHistNumberOfAcceptedTracks = new TH1F("fHistNumberOfAcceptedTracks",";N_{acc.};Entries",4001,-0.5,4000.5);
212 fList->Add(fHistNumberOfAcceptedTracks);
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);
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);
225 //========================================================//
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;
236 for(Int_t iEvent = 0; iEvent < nEvents; iEvent++) {
237 //for(Int_t iEvent = 0; iEvent < 1; iEvent++) {
238 eventChain->GetEntry(iEvent);
240 //========================================================//
241 //Create the TObjArray object to store the particles
242 TObjArray* tracksAccepted = new TObjArray;
243 tracksAccepted->SetOwner(kTRUE);
244 //========================================================//
246 Int_t nParticles = tStructEvents.entries;
247 if(((iEvent+1)%100)==0)
248 cout<<"Event: "<<iEvent+1<<" - ID: "<<tStructEvents.eventID<<" - Entries: "<<tStructEvents.entries<<endl;
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);
258 Int_t gNumberOfAcceptedParticles = 0;
259 Double_t gReactionPlane = 0.;
260 Double_t gCharge = 0.;
261 //========================================================//
262 //loop over particles
263 for(Int_t iParticle = 0; iParticle < nParticles; iParticle++) {
264 particleChain->GetEntry(nTotalParticles+iParticle);
266 //========================================================//
267 //consider only primordial particles
268 if(!IsPhysicalPrimary(tStructParticles.pid,
269 tStructParticles.fathereid,
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));
282 Double_t gPhi = TMath::Pi()+TMath::ATan2(-tStructParticles.py,-tStructParticles.px);
284 //========================================================//
287 fHistEtaPhiPositive->Fill(gEta,gPhi);
288 fHistPtPositive->Fill(gPt);
290 else if(gCharge < 0) {
291 fHistEtaPhiNegative->Fill(gEta,gPhi);
292 fHistPtNegative->Fill(gPt);
294 //========================================================//
296 //========================================================//
298 if((gEta > gEtaMax)||(gEta < gEtaMin)) continue;
299 if((gPt > gPtMax)||(gPt < gPtMin)) continue;
301 tracksAccepted->Add(new AliBFBasicParticle(gEta,gPhi,gPt,gCharge, 1.));
302 gNumberOfAcceptedParticles += 1;
304 iParticleCounter += 1;
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;
308 //========================================================//
309 // Event mixing (borrowed code from the task)
311 Int_t fMixingTracks = 50000;
312 AliEventPool* pool = fPoolMgr->GetEventPool(1.,0.,0.);
315 AliFatal(Form("No pool found"));
319 if (pool->IsReady() || pool->NTracksInPool() > fMixingTracks / 10 || pool->GetCurrentNEvents() >= 5){
321 Int_t nMix = pool->GetCurrentNEvents();
322 //cout << "nMix = " << nMix << " tracks in pool = " << pool->NTracksInPool() << endl;
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.);
331 // Update the Event pool
332 pool->UpdatePool(tracksAccepted);
337 //========================================================//
339 //========================================================//
340 // calculate balance function
341 bf->CalculateBalance(gReactionPlane,tracksAccepted,NULL,1.,0.);
343 fHistNumberOfAcceptedTracks->Fill(gNumberOfAcceptedParticles);
344 nTotalParticles += nParticles;
347 //========================================================//
349 TFile *f = TFile::Open("AnalysisResults.root","recreate");
350 TDirectoryFile *dir = new TDirectoryFile("PWGCFEbyE.outputBalanceFunctionPsiAnalysis","PWGCFEbyE.outputBalanceFunctionPsiAnalysis");
352 fList->SetName("listQA");
353 fList->SetOwner(kTRUE);
356 fListBF->SetName("listBF");
357 fListBF->SetOwner(kTRUE);
361 fListBFM->SetName("listBFMixed");
362 fListBFM->SetOwner(kTRUE);
366 dir->Write(dir->GetName(),TObject::kSingleKey);
368 //========================================================//
370 // Print real and CPU time used for analysis:
375 //______________________________________________________________//
376 Bool_t IsPhysicalPrimary(Int_t gPDGCode,
379 //Check whether the primordial particle belongs to the list
381 Bool_t kStatus = kFALSE;
385 //List of stable particles
386 const Int_t kNstable = 3;
387 Int_t pdgStable[kNstable] = {
393 if(gPDGCode < 0) gCharge = -1;
394 else if(gPDGCode > 0) gCharge = 1;
396 for(Int_t iParticle = 0; iParticle < kNstable; iParticle++) {
397 if((TMath::Abs(gPDGCode) == pdgStable[iParticle])) {