42cf1adcbc49672e92400c173f8aedfad9b73e8b
[u/mrichter/AliRoot.git] / PWGLF / FORWARD / analysis2 / AliBaseMCTrackDensity.cxx
1 #include "AliBaseMCTrackDensity.h"
2 #include <AliMCEvent.h>
3 #include <AliTrackReference.h>
4 #include <AliStack.h>
5 #include <TMath.h>
6 #include <AliLog.h>
7 #include <TH2D.h>
8 #include <TH1D.h>
9 #include <TList.h>
10 #include <TROOT.h>
11 #include <iostream>
12 #include "AliCollisionGeometry.h"
13 #include "AliGenEventHeader.h"
14 #include "AliForwardUtil.h"
15 #include <TF1.h>
16 #include <TGraph.h>
17
18 //____________________________________________________________________
19 AliBaseMCTrackDensity::AliBaseMCTrackDensity()
20   : TNamed(), 
21     fUseOnlyPrimary(false), 
22     fUseFlowWeights(false),
23     fBinFlow(0), 
24     fEtaBinFlow(0),
25     fPhiBinFlow(0),
26     fNRefs(0),
27     fV2Eta(0), 
28     fV22Pt(0), 
29     fV24Pt(0), 
30     fV2B(0),
31     fVz(0), 
32     fB(0),
33     fPhiR(0),
34     fDebug(false)
35 {
36   // Default constructor 
37   DGUARD(0,0,"MC track density default construction");
38 }
39
40 //____________________________________________________________________
41 AliBaseMCTrackDensity::AliBaseMCTrackDensity(const char* name)
42   : TNamed(name,"mcTrackDensity"), 
43     fUseOnlyPrimary(false), 
44     fUseFlowWeights(false),
45     fBinFlow(0), 
46     fEtaBinFlow(0),
47     fPhiBinFlow(0),
48     fNRefs(0),
49     fV2Eta(0), 
50     fV22Pt(0), 
51     fV24Pt(0), 
52     fV2B(0), 
53     fVz(0), 
54     fB(0),
55     fPhiR(0),
56     fDebug(false)
57 {
58   // Normal constructor constructor 
59   DGUARD(0,0,"MC track density named construction: %s", name);
60 }
61
62 //____________________________________________________________________
63 AliBaseMCTrackDensity::AliBaseMCTrackDensity(const AliBaseMCTrackDensity& o)
64   : TNamed(o),
65     fUseOnlyPrimary(o.fUseOnlyPrimary), 
66     fUseFlowWeights(o.fUseFlowWeights),
67     fBinFlow(o.fBinFlow), 
68     fEtaBinFlow(o.fEtaBinFlow),
69     fPhiBinFlow(o.fPhiBinFlow),
70     fNRefs(o.fNRefs),
71     fV2Eta(o.fV2Eta), 
72     fV22Pt(o.fV22Pt), 
73     fV24Pt(o.fV24Pt), 
74     fV2B(o.fV2B), 
75     fVz(o.fVz), 
76     fB(o.fB),
77     fPhiR(o.fPhiR),
78     fDebug(o.fDebug)
79 {
80   // Normal constructor constructor 
81   DGUARD(0,0,"MC track density copy construction");
82 }
83
84 //____________________________________________________________________
85 AliBaseMCTrackDensity&
86 AliBaseMCTrackDensity::operator=(const AliBaseMCTrackDensity& o)
87 {
88   // Assignment operator 
89   DGUARD(fDebug,3,"MC track density assignmetn");
90   if (&o == this) return *this; 
91   TNamed::operator=(o);
92   fUseOnlyPrimary       = o.fUseOnlyPrimary;
93   fBinFlow              = o.fBinFlow;
94   fEtaBinFlow           = o.fEtaBinFlow;
95   fPhiBinFlow           = o.fPhiBinFlow;
96   fNRefs                = o.fNRefs;
97   fDebug                = o.fDebug;
98   fUseFlowWeights       = o.fUseFlowWeights;
99   fV2Eta                = o.fV2Eta;
100   fV22Pt                = o.fV22Pt;
101   fV24Pt                = o.fV24Pt;
102   fV2B                  = o.fV2B;
103   fVz                   = o.fVz;
104   fB                    = o.fB;
105   fPhiR                 = o.fPhiR;
106   return *this;
107 }
108
109 //____________________________________________________________________
110 void
111 AliBaseMCTrackDensity::DefineOutput(TList* l)
112 {
113   DGUARD(fDebug,1,"MC track defines output");
114   TList* ll = new TList;
115   ll->SetName(GetTitle());
116   ll->SetOwner();
117   l->Add(ll);
118   
119   fBinFlow = new TH2D("binFlow", "#eta and #varphi bin flow", 
120                       200, -5, 5, 40, -180, 180);
121   fBinFlow->SetXTitle("#Delta#eta");
122   fBinFlow->SetYTitle("#Delta#varphi");
123   fBinFlow->SetOption("colz");
124   fBinFlow->SetDirectory(0);
125   ll->Add(fBinFlow);
126
127   fEtaBinFlow = new TH2D("binFlowEta", "#eta bin flow vs #eta", 
128                          200, -4, 6, 200, -5, 5);
129   fEtaBinFlow->SetXTitle("#eta");
130   fEtaBinFlow->SetYTitle("#Delta#eta");
131   fEtaBinFlow->SetOption("colz");
132   fEtaBinFlow->SetDirectory(0);
133   ll->Add(fEtaBinFlow);
134
135   fPhiBinFlow = new TH2D("binFlowPhi", "#phi bin flow vs #phi", 
136                          40, 0, 360, 40, -180, 180);
137   fPhiBinFlow->SetXTitle("#varphi");
138   fPhiBinFlow->SetYTitle("#Delta#varphi");
139   fPhiBinFlow->SetOption("colz");
140   fPhiBinFlow->SetDirectory(0);
141   ll->Add(fPhiBinFlow);
142
143   fNRefs = new TH1D("nRefs", "# references per track", 21, -.5, 20.5);
144   fNRefs->SetXTitle("# references");
145   fNRefs->SetFillColor(kMagenta+1);
146   fNRefs->SetFillStyle(3001);
147   fNRefs->SetDirectory(0);
148   ll->Add(fNRefs);
149
150   SetupWeights(ll);
151 }
152
153 //____________________________________________________________________
154 void 
155 AliBaseMCTrackDensity::SetupWeights(TList* l)
156 {
157   DGUARD(fDebug,2,"MC track setup phi weights");
158   fV2Eta = new TF1("v2eta", "gaus", -6, 6);
159   fV2Eta->SetParameters(20 * 0.1, 0, 9);
160   l->Add(fV2Eta);
161
162   Int_t          ptN     = 19;
163   const Double_t ptX[]   = {0.00,     0.25,     0.350,    0.45, 
164                             0.55,     0.650,    0.75,     0.85, 
165                             0.950,    1.10,     1.30,     1.500,
166                             1.70,     1.90,     2.250,    2.75, 
167                             3.25,     3.750,    4.50};
168   { 
169     // v2{2} dependence on pt
170     const Double_t y[] = {0.00000,  0.043400, 0.059911, 0.073516,
171                           0.089756, 0.105486, 0.117391, 0.128199,
172                           0.138013, 0.158271, 0.177726, 0.196383,
173                           0.208277, 0.216648, 0.242954, 0.249961,
174                           0.240131, 0.269006, 0.207796};
175     
176     fV22Pt = new TGraph(ptN, ptX, y);
177     fV22Pt->SetName("v2_2_vs_pt");
178     fV22Pt->SetMarkerStyle(20);
179     fV22Pt->SetMarkerColor(kRed+1);
180     l->Add(fV22Pt);
181   }
182
183   {
184     const Double_t y[] = {0.000000, 0.038646, 0.049824, 0.066662,
185                           0.075856, 0.081583, 0.099778, 0.104674,
186                           0.118545, 0.131874, 0.152959, 0.155348,
187                           0.169751, 0.179052, 0.178532, 0.198851,
188                           0.185737, 0.239901, 0.186098};
189
190     // v2{4} dependence on pt 
191     fV24Pt = new TGraph(ptN, ptX, y);
192     fV24Pt->SetName("v2_4_vs_pt");
193     fV24Pt->SetMarkerStyle(20);
194     fV24Pt->SetMarkerColor(kBlue+1);
195     l->Add(fV24Pt);
196   }
197   {
198     // V2 dependence on centrality
199     Int_t n            = 8;
200     const Double_t x[] = {1.75,     4.225,    5.965,    7.765,
201                           9.215,    10.46,    11.565,   12.575};
202     const Double_t y[] = {0.017855, 0.032440, 0.055818, 0.073137,
203                           0.083898, 0.086690, 0.082040, 0.077777};
204     fV2B = new TGraph(n, x, y);
205     fV2B->SetName("v2_vs_b");
206     fV2B->SetMarkerStyle(20);
207     fV2B->SetMarkerColor(kGreen+1);
208     l->Add(fV2B);
209   }
210 }
211
212
213 //____________________________________________________________________
214 Double_t
215 AliBaseMCTrackDensity::StoreParticle(AliMCParticle*       particle, 
216                                      const AliMCParticle* mother, 
217                                      AliTrackReference*   ref) const
218 {
219   DGUARD(fDebug,3,"MC track density store particle");
220   // Store a particle. 
221   if (!ref) return 0;
222
223   Double_t weight = 1;
224   if (fUseFlowWeights) {
225     Double_t phi = (mother ? mother->Phi() : particle->Phi());
226     Double_t eta = (mother ? mother->Eta() : particle->Eta());
227     Double_t pt  = (mother ? mother->Pt() : particle->Pt());
228     Int_t    id  = (mother ? mother->PdgCode() : 2212);
229     weight       = CalculateWeight(eta, pt, phi, id);
230   }
231
232   // Get track-reference stuff 
233   Double_t x      = ref->X();
234   Double_t y      = ref->Y();
235   Double_t z      = ref->Z()-fVz;
236   Double_t rr     = TMath::Sqrt(x*x+y*y);
237   Double_t thetaR = TMath::ATan2(rr,z);
238   Double_t phiR   = TMath::ATan2(y,x);
239   Double_t etaR   = -TMath::Log(TMath::Tan(thetaR/2));
240   
241   // Correct angle and convert to degrees 
242   if (thetaR < 0) thetaR += 2*TMath::Pi();
243   thetaR *= 180. / TMath::Pi();
244   if (phiR < 0) phiR += 2*TMath::Pi();
245   phiR *= 180. / TMath::Pi();
246
247   const AliMCParticle* mp = (mother ? mother : particle);
248   Double_t dEta = mp->Eta() - etaR;
249   Double_t dPhi = mp->Phi() * 180 / TMath::Pi() - phiR;
250   if (dPhi >  180) dPhi -= 360;
251   if (dPhi < -180) dPhi += 360;
252   fBinFlow->Fill(dEta, dPhi);
253   fEtaBinFlow->Fill(etaR, dEta);
254   fPhiBinFlow->Fill(phiR, dPhi);
255
256   return weight;
257 }
258
259 //____________________________________________________________________
260 Double_t
261 AliBaseMCTrackDensity::GetTrackRefTheta(const AliTrackReference* ref) const
262 {
263   // Get the incidient angle of the track reference. 
264   Double_t x    = ref->X();
265   Double_t y    = ref->Y();
266   Double_t z    = ref->Z()-fVz;
267   Double_t rr   = TMath::Sqrt(x*x+y*y);
268   Double_t theta= TMath::ATan2(rr,z);
269   Double_t ang  = TMath::Abs(TMath::Pi()-theta);
270   return ang;
271 }
272                                     
273 //____________________________________________________________________
274 const AliMCParticle*
275 AliBaseMCTrackDensity::GetMother(Int_t     iTr,
276                                 const AliMCEvent& event) const
277 {
278   // 
279   // Track down primary mother 
280   // 
281   Int_t i  = iTr;
282   do { 
283     const AliMCParticle* p = static_cast<AliMCParticle*>(event.GetTrack(i));
284     if (const_cast<AliMCEvent&>(event).Stack()->IsPhysicalPrimary(i)) return p;
285     
286     i = p->GetMother();
287   } while (i > 0);
288
289   return 0;
290 }  
291
292 //____________________________________________________________________
293 Bool_t
294 AliBaseMCTrackDensity::GetCollisionParameters(const AliMCEvent& event)
295
296   DGUARD(fDebug,3,"MC track density get collision parameters");
297   AliCollisionGeometry* hd = 
298     dynamic_cast<AliCollisionGeometry*>(event.GenEventHeader());
299   fPhiR = (hd ? hd->ReactionPlaneAngle() : 0.);
300   fB    = (hd ? hd->ImpactParameter() : -1 );
301   return hd != 0;
302 }
303
304 //____________________________________________________________________
305 Bool_t
306 AliBaseMCTrackDensity::ProcessTrack(AliMCParticle* particle, 
307                                     const AliMCParticle* mother)
308 {
309   // Check the returned particle 
310   DGUARD(fDebug,3,"MC track density Process a track");
311   if (!particle) return false;
312     
313   Int_t              nTrRef = particle->GetNumberOfTrackReferences();
314   AliTrackReference* store  = 0;
315
316   BeginTrackRefs();
317
318   // Double_t oTheta= 0;
319   Int_t nRefs = 0;
320   for (Int_t iTrRef = 0; iTrRef < nTrRef; iTrRef++) { 
321     AliTrackReference* ref = particle->GetTrackReference(iTrRef);
322       
323     // Check existence 
324     if (!ref) continue;
325
326     // Check that we hit an Base element 
327     if (ref->DetectorId() != GetDetectorId()) continue;
328     if (!CheckTrackRef(ref)) continue;
329
330     nRefs++;
331
332     AliTrackReference* test = ProcessRef(particle, mother, ref);
333     if (test) store = test;
334
335   } // Loop over track references
336   if (!store) return true; // Nothing found
337   
338   StoreParticle(particle, mother, store);
339
340   fNRefs->Fill(nRefs);
341
342   EndTrackRefs(nRefs);
343
344   return true;
345 }
346
347
348 //____________________________________________________________________
349 Bool_t
350 AliBaseMCTrackDensity::ProcessTracks(const AliMCEvent& event,
351                                      Double_t          vz,
352                                      TH2D*             primary)
353 {
354   // 
355   // Filter the input kinematics and track references, using 
356   // some of the ESD information
357   // 
358   // Parameters:
359   //    input   Input ESD event
360   //    event   Input MC event
361   //    vz      Vertex position 
362   //    output  Output ESD-like object
363   //    primary Per-event histogram of primaries 
364   //
365   // Return:
366   //    True on succes, false otherwise 
367   //
368   DGUARD(fDebug,3,"MC track density Process a tracks");
369   fVz = vz;
370   GetCollisionParameters(event);
371   
372   AliStack* stack = const_cast<AliMCEvent&>(event).Stack();
373   Int_t nTracks   = stack->GetNtrack();//event.GetNumberOfTracks();
374   Int_t nPrim     = stack->GetNprimary();//event.GetNumberOfPrimary();
375   for (Int_t iTr = 0; iTr < nTracks; iTr++) { 
376     AliMCParticle* particle = 
377       static_cast<AliMCParticle*>(event.GetTrack(iTr));
378     
379     // Check if this charged and a primary 
380     if (!particle->Charge() != 0) continue;
381     
382     Bool_t isPrimary = stack->IsPhysicalPrimary(iTr) && iTr < nPrim;
383     
384     // Fill 'dn/deta' histogram 
385     if (isPrimary && primary) 
386       primary->Fill(particle->Eta(), particle->Phi());
387
388     // Bail out if we're only processing primaries - perhaps we should
389     // track back to the original primary?
390     if (fUseOnlyPrimary && !isPrimary) continue;
391
392     const AliMCParticle* mother = GetMother(iTr, event);
393     ProcessTrack(particle, mother);
394
395   } // Loop over tracks
396   return kTRUE;
397 }
398
399   
400 //____________________________________________________________________
401 Double_t
402 AliBaseMCTrackDensity::CalculateWeight(Double_t eta, Double_t pt, 
403                                        Double_t phi, Int_t id) const
404 {
405   Double_t w = (fV2Eta->Eval(eta) * 0.5 * (fV22Pt->Eval(pt) + fV24Pt->Eval(pt)) 
406                 * fV2B->Eval(fB) / fV2B->Eval(10.46) 
407                 * 2 * TMath::Cos(2* (phi - fPhiR)));
408   UInt_t aid = TMath::Abs(id);
409   switch (aid) { 
410   case 211:  w *= 1.3; break; // pions 
411   case 2212: w *= 1.0; break; // protons 
412   default:   w *= 0.7; break;
413   }
414   return w;
415 }
416 //____________________________________________________________________
417 void
418 AliBaseMCTrackDensity::Print(Option_t* /*option*/) const 
419 {
420   char ind[gROOT->GetDirLevel()+1];
421   for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
422   ind[gROOT->GetDirLevel()] = '\0';
423   std::cout << ind << ClassName() << ": " << GetName() << '\n'
424             << std::boolalpha 
425             << ind << " Only primary tracks:    " << fUseOnlyPrimary << '\n'
426             << ind << " Use flow after burner:  " << fUseFlowWeights 
427             << std::noboolalpha << std::endl;
428   
429 }
430
431 //____________________________________________________________________
432 //
433 // EOF
434 //