bede67a042a829e65d77695cec71bebdb91f928d
[u/mrichter/AliRoot.git] / PWGCF / EBYE / macros / runBalanceFunctionOnHydro.C
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
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)
38 };
39 //========================================================//
40
41 //========================================================//
42 //Balance function analysis variables
43 Bool_t gRunShuffling=kFALSE;
44 Bool_t gRunMixing=kTRUE;
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
54 void runBalanceFunctionOnHydro(TString aEventDir = ".", 
55                                Int_t aEventFiles = 9) {
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();
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  
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
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;
166     
167     AliEventPoolManager *fPoolMgr = 0x0;
168     fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centbins, nVertexBins, vtxbins);
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
188   for(Int_t iFile = 1; iFile <= aEventFiles; iFile++) {
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   //========================================================//
199   //QA histograms
200   //Event stats.
201   TString gCutName[5] = {"Total","Offline trigger",
202                          "Vertex","Analyzed","sel. Centrality"};
203   TH1F *fHistEventStats = new TH1F("fHistEventStats",
204                                    "Event statistics;;N_{events}",
205                                    5,0.5,5.5);
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
211   TH1F *fHistNumberOfAcceptedTracks = new TH1F("fHistNumberOfAcceptedTracks",";N_{acc.};Entries",4001,-0.5,4000.5);
212   fList->Add(fHistNumberOfAcceptedTracks);
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);
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   
236   for(Int_t iEvent = 0; iEvent < nEvents; iEvent++) {
237     //for(Int_t iEvent = 0; iEvent < 1; iEvent++) {
238     eventChain->GetEntry(iEvent);
239
240     //========================================================//
241     //Create the TObjArray object to store the particles
242     TObjArray* tracksAccepted = new TObjArray;
243     tracksAccepted->SetOwner(kTRUE);
244     //========================================================//
245
246     Int_t nParticles = tStructEvents.entries;
247     if(((iEvent+1)%100)==0)
248       cout<<"Event: "<<iEvent+1<<" - ID: "<<tStructEvents.eventID<<" - Entries: "<<tStructEvents.entries<<endl;
249
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.;
260     Double_t gCharge = 0.;
261     //========================================================//
262     //loop over particles
263     for(Int_t iParticle = 0; iParticle < nParticles; iParticle++) {
264       particleChain->GetEntry(nTotalParticles+iParticle);
265       
266       //========================================================//
267       //consider only primordial particles
268       if(!IsPhysicalPrimary(tStructParticles.pid,
269                             tStructParticles.fathereid,
270                             gCharge)) continue;
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));
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
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       
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;
306     }//particle loop
307     
308     //========================================================//
309     // Event mixing (borrowed code from the task) 
310     if (gRunMixing) {
311       Int_t fMixingTracks = 50000;
312       AliEventPool* pool = fPoolMgr->GetEventPool(1.,0.,0.);
313       
314       if (!pool) {
315         AliFatal(Form("No pool found"));
316       }
317       else {
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         }
330       
331         // Update the Event pool
332         pool->UpdatePool(tracksAccepted);
333         //pool->PrintInfo();
334         
335       }//pool NULL check  
336     }//run mixing
337       //========================================================//
338
339     //========================================================//
340     // calculate balance function
341     bf->CalculateBalance(gReactionPlane,tracksAccepted,NULL,1.,0.);
342
343     fHistNumberOfAcceptedTracks->Fill(gNumberOfAcceptedParticles);
344     nTotalParticles += nParticles;
345   }
346
347   //========================================================//
348   //Output file
349   TFile *f = TFile::Open("AnalysisResults.root","recreate");
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);
367   f->Close();  
368   //========================================================//
369
370   // Print real and CPU time used for analysis:  
371   timer.Stop();
372   timer.Print();
373 }
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 }