Moved calculations of flow weights (after burner) into a
[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     fWeights(),
28     fVz(0), 
29     fB(0),
30     fPhiR(0),
31     fDebug(false)
32 {
33   // Default constructor 
34   DGUARD(0,0,"MC track density default construction");
35 }
36
37 //____________________________________________________________________
38 AliBaseMCTrackDensity::AliBaseMCTrackDensity(const char* name)
39   : TNamed(name,"mcTrackDensity"), 
40     fUseOnlyPrimary(false), 
41     fUseFlowWeights(false),
42     fBinFlow(0), 
43     fEtaBinFlow(0),
44     fPhiBinFlow(0),
45     fNRefs(0),
46     fWeights(),
47     fVz(0), 
48     fB(0),
49     fPhiR(0),
50     fDebug(false)
51 {
52   // Normal constructor constructor 
53   DGUARD(0,0,"MC track density named construction: %s", name);
54 }
55
56 //____________________________________________________________________
57 AliBaseMCTrackDensity::AliBaseMCTrackDensity(const AliBaseMCTrackDensity& o)
58   : TNamed(o),
59     fUseOnlyPrimary(o.fUseOnlyPrimary), 
60     fUseFlowWeights(o.fUseFlowWeights),
61     fBinFlow(o.fBinFlow), 
62     fEtaBinFlow(o.fEtaBinFlow),
63     fPhiBinFlow(o.fPhiBinFlow),
64     fNRefs(o.fNRefs),
65     fWeights(o.fWeights),
66     fVz(o.fVz), 
67     fB(o.fB),
68     fPhiR(o.fPhiR),
69     fDebug(o.fDebug)
70 {
71   // Normal constructor constructor 
72   DGUARD(0,0,"MC track density copy construction");
73 }
74
75 //____________________________________________________________________
76 AliBaseMCTrackDensity&
77 AliBaseMCTrackDensity::operator=(const AliBaseMCTrackDensity& o)
78 {
79   // Assignment operator 
80   DGUARD(fDebug,3,"MC track density assignmetn");
81   if (&o == this) return *this; 
82   TNamed::operator=(o);
83   fUseOnlyPrimary       = o.fUseOnlyPrimary;
84   fBinFlow              = o.fBinFlow;
85   fEtaBinFlow           = o.fEtaBinFlow;
86   fPhiBinFlow           = o.fPhiBinFlow;
87   fNRefs                = o.fNRefs;
88   fDebug                = o.fDebug;
89   fUseFlowWeights       = o.fUseFlowWeights;
90   fWeights              = o.fWeights;
91   fVz                   = o.fVz;
92   fB                    = o.fB;
93   fPhiR                 = o.fPhiR;
94   return *this;
95 }
96
97 //____________________________________________________________________
98 void
99 AliBaseMCTrackDensity::DefineOutput(TList* l)
100 {
101   DGUARD(fDebug,1,"MC track defines output");
102   TList* ll = new TList;
103   ll->SetName(GetTitle());
104   ll->SetOwner();
105   l->Add(ll);
106   
107   fBinFlow = new TH2D("binFlow", "#eta and #varphi bin flow", 
108                       200, -5, 5, 40, -180, 180);
109   fBinFlow->SetXTitle("#Delta#eta");
110   fBinFlow->SetYTitle("#Delta#varphi");
111   fBinFlow->SetOption("colz");
112   fBinFlow->SetDirectory(0);
113   ll->Add(fBinFlow);
114
115   fEtaBinFlow = new TH2D("binFlowEta", "#eta bin flow vs #eta", 
116                          200, -4, 6, 200, -5, 5);
117   fEtaBinFlow->SetXTitle("#eta");
118   fEtaBinFlow->SetYTitle("#Delta#eta");
119   fEtaBinFlow->SetOption("colz");
120   fEtaBinFlow->SetDirectory(0);
121   ll->Add(fEtaBinFlow);
122
123   fPhiBinFlow = new TH2D("binFlowPhi", "#phi bin flow vs #phi", 
124                          40, 0, 360, 40, -180, 180);
125   fPhiBinFlow->SetXTitle("#varphi");
126   fPhiBinFlow->SetYTitle("#Delta#varphi");
127   fPhiBinFlow->SetOption("colz");
128   fPhiBinFlow->SetDirectory(0);
129   ll->Add(fPhiBinFlow);
130
131   fNRefs = new TH1D("nRefs", "# references per track", 21, -.5, 20.5);
132   fNRefs->SetXTitle("# references");
133   fNRefs->SetFillColor(kMagenta+1);
134   fNRefs->SetFillStyle(3001);
135   fNRefs->SetDirectory(0);
136   ll->Add(fNRefs);
137   
138   fWeights.Init(ll);
139 }
140
141
142 //____________________________________________________________________
143 Double_t
144 AliBaseMCTrackDensity::StoreParticle(AliMCParticle*       particle, 
145                                      const AliMCParticle* mother, 
146                                      AliTrackReference*   ref) const
147 {
148   DGUARD(fDebug,3,"MC track density store particle");
149   // Store a particle. 
150   if (!ref) return 0;
151
152   Double_t weight = 1;
153   if (fUseFlowWeights) {
154     Double_t phi = (mother ? mother->Phi() : particle->Phi());
155     Double_t eta = (mother ? mother->Eta() : particle->Eta());
156     Double_t pt  = (mother ? mother->Pt() : particle->Pt());
157     Int_t    id  = (mother ? mother->PdgCode() : 2212);
158     weight       = CalculateWeight(eta, pt, phi, id);
159   }
160
161   // Get track-reference stuff 
162   Double_t x      = ref->X();
163   Double_t y      = ref->Y();
164   Double_t z      = ref->Z()-fVz;
165   Double_t rr     = TMath::Sqrt(x*x+y*y);
166   Double_t thetaR = TMath::ATan2(rr,z);
167   Double_t phiR   = TMath::ATan2(y,x);
168   Double_t etaR   = -TMath::Log(TMath::Tan(thetaR/2));
169   
170   // Correct angle and convert to degrees 
171   if (thetaR < 0) thetaR += 2*TMath::Pi();
172   thetaR *= 180. / TMath::Pi();
173   if (phiR < 0) phiR += 2*TMath::Pi();
174   phiR *= 180. / TMath::Pi();
175
176   const AliMCParticle* mp = (mother ? mother : particle);
177   Double_t dEta = mp->Eta() - etaR;
178   Double_t dPhi = mp->Phi() * 180 / TMath::Pi() - phiR;
179   if (dPhi >  180) dPhi -= 360;
180   if (dPhi < -180) dPhi += 360;
181   fBinFlow->Fill(dEta, dPhi);
182   fEtaBinFlow->Fill(etaR, dEta);
183   fPhiBinFlow->Fill(phiR, dPhi);
184
185   return weight;
186 }
187
188 //____________________________________________________________________
189 Double_t
190 AliBaseMCTrackDensity::GetTrackRefTheta(const AliTrackReference* ref) const
191 {
192   // Get the incidient angle of the track reference. 
193   Double_t x    = ref->X();
194   Double_t y    = ref->Y();
195   Double_t z    = ref->Z()-fVz;
196   Double_t rr   = TMath::Sqrt(x*x+y*y);
197   Double_t theta= TMath::ATan2(rr,z);
198   Double_t ang  = TMath::Abs(TMath::Pi()-theta);
199   return ang;
200 }
201                                     
202 //____________________________________________________________________
203 const AliMCParticle*
204 AliBaseMCTrackDensity::GetMother(Int_t     iTr,
205                                 const AliMCEvent& event) const
206 {
207   // 
208   // Track down primary mother 
209   // 
210   Int_t i  = iTr;
211   do { 
212     const AliMCParticle* p = static_cast<AliMCParticle*>(event.GetTrack(i));
213     if (const_cast<AliMCEvent&>(event).Stack()->IsPhysicalPrimary(i)) return p;
214     
215     i = p->GetMother();
216   } while (i > 0);
217
218   return 0;
219 }  
220
221 //____________________________________________________________________
222 Bool_t
223 AliBaseMCTrackDensity::GetCollisionParameters(const AliMCEvent& event)
224
225   DGUARD(fDebug,3,"MC track density get collision parameters");
226   AliCollisionGeometry* hd = 
227     dynamic_cast<AliCollisionGeometry*>(event.GenEventHeader());
228   fPhiR = (hd ? hd->ReactionPlaneAngle() : 0.);
229   fB    = (hd ? hd->ImpactParameter() : -1 );
230   return hd != 0;
231 }
232
233 //____________________________________________________________________
234 Bool_t
235 AliBaseMCTrackDensity::ProcessTrack(AliMCParticle* particle, 
236                                     const AliMCParticle* mother)
237 {
238   // Check the returned particle 
239   DGUARD(fDebug,3,"MC track density Process a track");
240   if (!particle) return false;
241     
242   Int_t              nTrRef = particle->GetNumberOfTrackReferences();
243   AliTrackReference* store  = 0;
244
245   BeginTrackRefs();
246
247   // Double_t oTheta= 0;
248   Int_t nRefs = 0;
249   for (Int_t iTrRef = 0; iTrRef < nTrRef; iTrRef++) { 
250     AliTrackReference* ref = particle->GetTrackReference(iTrRef);
251       
252     // Check existence 
253     if (!ref) continue;
254
255     // Check that we hit an Base element 
256     if (ref->DetectorId() != GetDetectorId()) continue;
257     if (!CheckTrackRef(ref)) continue;
258
259     nRefs++;
260
261     AliTrackReference* test = ProcessRef(particle, mother, ref);
262     if (test) store = test;
263
264   } // Loop over track references
265   if (!store) return true; // Nothing found
266   
267   StoreParticle(particle, mother, store);
268
269   fNRefs->Fill(nRefs);
270
271   EndTrackRefs(nRefs);
272
273   return true;
274 }
275
276
277 //____________________________________________________________________
278 Bool_t
279 AliBaseMCTrackDensity::ProcessTracks(const AliMCEvent& event,
280                                      Double_t          vz,
281                                      TH2D*             primary)
282 {
283   // 
284   // Filter the input kinematics and track references, using 
285   // some of the ESD information
286   // 
287   // Parameters:
288   //    input   Input ESD event
289   //    event   Input MC event
290   //    vz      Vertex position 
291   //    output  Output ESD-like object
292   //    primary Per-event histogram of primaries 
293   //
294   // Return:
295   //    True on succes, false otherwise 
296   //
297   DGUARD(fDebug,3,"MC track density Process a tracks");
298   fVz = vz;
299   GetCollisionParameters(event);
300   
301   AliStack* stack = const_cast<AliMCEvent&>(event).Stack();
302   Int_t nTracks   = stack->GetNtrack();//event.GetNumberOfTracks();
303   Int_t nPrim     = stack->GetNprimary();//event.GetNumberOfPrimary();
304   for (Int_t iTr = 0; iTr < nTracks; iTr++) { 
305     AliMCParticle* particle = 
306       static_cast<AliMCParticle*>(event.GetTrack(iTr));
307     
308     // Check if this charged and a primary 
309     if (!particle->Charge() != 0) continue;
310     
311     Bool_t isPrimary = stack->IsPhysicalPrimary(iTr) && iTr < nPrim;
312     
313     // Fill 'dn/deta' histogram 
314     if (isPrimary && primary) 
315       primary->Fill(particle->Eta(), particle->Phi());
316
317     // Bail out if we're only processing primaries - perhaps we should
318     // track back to the original primary?
319     if (fUseOnlyPrimary && !isPrimary) continue;
320
321     const AliMCParticle* mother = GetMother(iTr, event);
322     ProcessTrack(particle, mother);
323
324   } // Loop over tracks
325   return kTRUE;
326 }
327
328   
329 //____________________________________________________________________
330 Double_t
331 AliBaseMCTrackDensity::CalculateWeight(Double_t eta, Double_t pt, 
332                                        Double_t phi, Int_t id) const
333 {
334   return fWeights.CalcWeight(eta, pt, phi, id, fPhiR, fB);
335 }
336 //____________________________________________________________________
337 void
338 AliBaseMCTrackDensity::Print(Option_t* /*option*/) const 
339 {
340   char ind[gROOT->GetDirLevel()+1];
341   for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
342   ind[gROOT->GetDirLevel()] = '\0';
343   std::cout << ind << ClassName() << ": " << GetName() << '\n'
344             << std::boolalpha 
345             << ind << " Only primary tracks:    " << fUseOnlyPrimary << '\n'
346             << ind << " Use flow after burner:  " << fUseFlowWeights 
347             << std::noboolalpha << std::endl;
348   
349 }
350
351 //____________________________________________________________________
352 //
353 // EOF
354 //