AliPhysics  58f3d52 (58f3d52)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliBaseMCTrackDensity.cxx
Go to the documentation of this file.
3 #include <AliMCEvent.h>
4 #include <AliTrackReference.h>
5 //#include <AliStack.h>
6 #include <TMath.h>
7 #include <AliLog.h>
8 #include <TH2D.h>
9 #include <TH1D.h>
10 #include <TList.h>
11 #include <TROOT.h>
12 #include <TVector3.h>
13 #include <iostream>
14 #include "AliCollisionGeometry.h"
15 #include "AliGenEventHeader.h"
16 #include "AliForwardUtil.h"
17 #include <TF1.h>
18 #include <TGraph.h>
19 
20 //____________________________________________________________________
22  : TNamed(),
23  fUseOnlyPrimary(false),
24  fBinFlow(0),
25  fEtaBinFlow(0),
26  fPhiBinFlow(0),
27  fNRefs(0),
28  fWeights(0),
29  fTruthWeights(0),
30  fIP(0,0,0),
31  fB(0),
32  fPhiR(0),
33  fDebug(false),
34  fTrackGammaToPi0(false)
35 {
36  // Default constructor
37  DGUARD(fDebug, 3,"Default CTOR of AliBasMCTrackDensity");
38 }
39 
40 //____________________________________________________________________
42  : TNamed(name,"mcTrackDensity"),
43  fUseOnlyPrimary(false),
44  fBinFlow(0),
45  fEtaBinFlow(0),
46  fPhiBinFlow(0),
47  fNRefs(0),
48  fWeights(0),
49  fTruthWeights(0),
50  fIP(0,0,0),
51  fB(0),
52  fPhiR(0),
53  fDebug(false),
54  fTrackGammaToPi0(false)
55 {
56  // Normal constructor constructor
57  DGUARD(fDebug, 3,"Named CTOR of AliBasMCTrackDensity: %s", name);
58 }
59 
60 //____________________________________________________________________
62  : TNamed(o),
63  fUseOnlyPrimary(o.fUseOnlyPrimary),
64  fBinFlow(o.fBinFlow),
65  fEtaBinFlow(o.fEtaBinFlow),
66  fPhiBinFlow(o.fPhiBinFlow),
67  fNRefs(o.fNRefs),
68  fWeights(o.fWeights),
69  fTruthWeights(o.fTruthWeights),
70  fIP(o.fIP),
71  fB(o.fB),
72  fPhiR(o.fPhiR),
73  fDebug(o.fDebug),
74  fTrackGammaToPi0(o.fTrackGammaToPi0)
75 {
76  // Normal constructor constructor
77  DGUARD(fDebug, 3,"Copy CTOR of AliBasMCTrackDensity");
78 }
79 
80 //____________________________________________________________________
83 {
84  // Assignment operator
85  DGUARD(fDebug,3,"MC track density assignmetn");
86  if (&o == this) return *this;
87  TNamed::operator=(o);
89  fBinFlow = o.fBinFlow;
92  fNRefs = o.fNRefs;
93  fDebug = o.fDebug;
94  fWeights = o.fWeights;
96  fIP.SetXYZ(o.fIP.X(), o.fIP.Y(), o.fIP.Z());
97  fB = o.fB;
98  fPhiR = o.fPhiR;
100  return *this;
101 }
102 
103 //____________________________________________________________________
104 void
106 {
107  if (fWeights) {
108  delete fWeights;
109  fWeights = 0;
110  }
111  fWeights = w;
112 }
113 //____________________________________________________________________
114 void
116 {
117  if (fTruthWeights) {
118  delete fTruthWeights;
119  fTruthWeights = 0;
120  }
121  fTruthWeights = w;
122 }
123 //____________________________________________________________________
124 void
126 {
127  if (!use) return;
128 
130 }
131 //____________________________________________________________________
132 void
134 {
135  DGUARD(fDebug,1,"MC track defines output");
136  TList* ll = new TList;
137  ll->SetName(GetTitle());
138  ll->SetOwner();
139  l->Add(ll);
140 
141  fBinFlow = new TH2D("binFlow", "#eta and #varphi bin flow",
142  200, -5, 5, 40, -180, 180);
143  fBinFlow->SetXTitle("#Delta#eta");
144  fBinFlow->SetYTitle("#Delta#varphi");
145  fBinFlow->SetOption("colz");
146  fBinFlow->SetDirectory(0);
147  ll->Add(fBinFlow);
148 
149  fEtaBinFlow = new TH2D("binFlowEta", "#eta bin flow vs #eta",
150  200, -4, 6, 200, -5, 5);
151  fEtaBinFlow->SetXTitle("#eta");
152  fEtaBinFlow->SetYTitle("#Delta#eta");
153  fEtaBinFlow->SetOption("colz");
154  fEtaBinFlow->SetDirectory(0);
155  ll->Add(fEtaBinFlow);
156 
157  fPhiBinFlow = new TH2D("binFlowPhi", "#phi bin flow vs #phi",
158  40, 0, 360, 40, -180, 180);
159  fPhiBinFlow->SetXTitle("#varphi");
160  fPhiBinFlow->SetYTitle("#Delta#varphi");
161  fPhiBinFlow->SetOption("colz");
162  fPhiBinFlow->SetDirectory(0);
163  ll->Add(fPhiBinFlow);
164 
165  fNRefs = new TH1D("nRefs", "# references per track", 21, -.5, 20.5);
166  fNRefs->SetXTitle("# references");
167  fNRefs->SetFillColor(kMagenta+1);
168  fNRefs->SetFillStyle(3001);
169  fNRefs->SetDirectory(0);
170  ll->Add(fNRefs);
171 
172  if(fWeights) fWeights->Init(ll);
173 }
174 
175 
176 //____________________________________________________________________
177 Double_t
178 AliBaseMCTrackDensity::StoreParticle(AliMCParticle* particle,
179  const AliMCParticle* mother,
180  AliTrackReference* ref) const
181 {
182  // Store this particle
183  //
184  // Note: If particle refers to a primary, then particle and mother
185  // refers to the same particle (the address are the same)
186  //
187  DGUARD(fDebug,3,"MC track density store particle");
188  // Store a particle.
189  if (!ref) return 0;
190 
191  Double_t weight = 1;
192  if (fWeights)
193  weight = CalculateWeight(mother, (mother == particle));
194 
195 
196  // Get track-reference stuff
197  Double_t x = ref->X()-fIP.X();
198  Double_t y = ref->Y()-fIP.Y();
199  Double_t z = ref->Z()-fIP.Z();
200  Double_t rr = TMath::Sqrt(x*x+y*y);
201  Double_t thetaR = TMath::ATan2(rr,z);
202  Double_t phiR = TMath::ATan2(y,x);
203  Double_t etaR = -TMath::Log(TMath::Tan(thetaR/2));
204 
205  // Correct angle and convert to degrees
206  if (thetaR < 0) thetaR += 2*TMath::Pi();
207  thetaR *= 180. / TMath::Pi();
208  if (phiR < 0) phiR += 2*TMath::Pi();
209  phiR *= 180. / TMath::Pi();
210 
211  const AliMCParticle* mp = (mother ? mother : particle);
212  Double_t dEta = mp->Eta() - etaR;
213  Double_t dPhi = mp->Phi() * 180 / TMath::Pi() - phiR;
214  if (dPhi > 180) dPhi -= 360;
215  if (dPhi < -180) dPhi += 360;
216  fBinFlow->Fill(dEta, dPhi);
217  fEtaBinFlow->Fill(etaR, dEta);
218  fPhiBinFlow->Fill(phiR, dPhi);
219  return weight;
220 }
221 
222 //____________________________________________________________________
223 Double_t
224 AliBaseMCTrackDensity::GetTrackRefTheta(const AliTrackReference* ref) const
225 {
226  // Get the incidient angle of the track reference.
227  Double_t x = ref->X()-fIP.X();
228  Double_t y = ref->Y()-fIP.Y();
229  Double_t z = ref->Z()-fIP.Z();
230  Double_t rr = TMath::Sqrt(x*x+y*y);
231  Double_t theta= TMath::ATan2(rr,z);
232  Double_t ang = TMath::Abs(TMath::Pi()-theta);
233  return ang;
234 }
235 
236 //____________________________________________________________________
237 const AliMCParticle*
238 AliBaseMCTrackDensity::GetMother(Int_t iTr, const AliMCEvent& event) const
239 {
240  //
241  // Track down primary mother
242  //
243  Int_t i = iTr;
244  Bool_t gammaSeen = false;
245  const AliMCParticle* candidate = 0;
246  do {
247  const AliMCParticle* p = static_cast<AliMCParticle*>(event.GetTrack(i));
248  if (!p) break;
249  if (gammaSeen && TMath::Abs(p->PdgCode()) == 111)
250  // If we're looking for a mother pi0 of gamma, and we find it
251  // here, we return it - irrespective of whether it's flagged as
252  // a primary or not.
253  return p;
254 
255  if (const_cast<AliMCEvent&>(event).IsPhysicalPrimary(i)) {
256  candidate = p;
257  if (fTrackGammaToPi0 && TMath::Abs(p->PdgCode()) == 22)
258  // If we want to track gammas back to a possible pi0, we flag
259  // the gamma seen, and store it as a candidate in case we do
260  // not find a pi0 in the stack
261  gammaSeen = true;
262  else
263  break;
264  }
265 
266  // We get here if the current track isn't a primary, or it was a
267  // primary gamma and we want to track back to a pi0.
268  i = p->GetMother();
269  } while (i > 0);
270 
271  // Return our candidate (gamma) if we find no mother pi0. Note, we
272  // should never get here with a null pointer, so we issue a warning
273  // in that case.
274  if (!candidate)
275  AliWarningF("No mother found for track # %d", iTr);
276  return candidate;
277 }
278 
279 //____________________________________________________________________
280 Bool_t
282 {
283  DGUARD(fDebug,3,"MC track density get collision parameters");
284  AliCollisionGeometry* hd =
285  dynamic_cast<AliCollisionGeometry*>(event.GenEventHeader());
286  fPhiR = (hd ? hd->ReactionPlaneAngle() : 0.);
287  fB = (hd ? hd->ImpactParameter() : -1 );
288  return hd != 0;
289 }
290 
291 //____________________________________________________________________
292 Bool_t
293 AliBaseMCTrackDensity::ProcessTrack(AliMCParticle* particle,
294  const AliMCParticle* mother)
295 {
296  // Check the returned particle
297  //
298  // Note: If particle refers to a primary, then particle and mother
299  // refers to the same particle (the address are the same)
300  //
301  DGUARD(fDebug,3,"MC track density Process a track");
302  if (!particle) return false;
303 
304  Int_t nTrRef = particle->GetNumberOfTrackReferences();
305  AliTrackReference* store = 0;
306 
307  BeginTrackRefs();
308 
309  // Double_t oTheta= 0;
310  Int_t nRefs = 0;
311  for (Int_t iTrRef = 0; iTrRef < nTrRef; iTrRef++) {
312  AliTrackReference* ref = particle->GetTrackReference(iTrRef);
313 
314  // Check existence
315  if (!ref) continue;
316 
317  // Check that we hit an Base element
318  if (ref->DetectorId() != GetDetectorId()) continue;
319  if (!CheckTrackRef(ref)) continue;
320 
321  nRefs++;
322 
323  AliTrackReference* test = ProcessRef(particle, mother, ref);
324  if (test) store = test;
325 
326  } // Loop over track references
327  if (!store) return true; // Nothing found
328 
329  StoreParticle(particle, mother, store);
330 
331  fNRefs->Fill(nRefs);
332 
333  EndTrackRefs(nRefs);
334 
335  return true;
336 }
337 
338 
339 //____________________________________________________________________
340 Bool_t
341 AliBaseMCTrackDensity::ProcessTracks(const AliMCEvent& event,
342  const TVector3& ip,
343  TH2D* primary)
344 {
345  //
346  // Filter the input kinematics and track references, using
347  // some of the ESD information
348  //
349  // Parameters:
350  // input Input ESD event
351  // event Input MC event
352  // vz Vertex position
353  // output Output ESD-like object
354  // primary Per-event histogram of primaries
355  //
356  // Return:
357  // True on succes, false otherwise
358  //
359  DGUARD(fDebug,3,"MC track density Process a tracks");
360  fIP.SetXYZ(ip.X(), ip.Y(), ip.Z());
361  GetCollisionParameters(event);
362 
363  // AliStack* stack = const_cast<AliMCEvent&>(event).Stack();
364  // if (!stack) return kFALSE;
365 
366 
367  Int_t nTracks = /*stack->GetNtrack();*/event.GetNumberOfTracks();
368  // temporary remove constness, since event.GetNumberOfPrimaries() was not const
369  Int_t nPrim = /*stack->GetNprimary();*/((AliMCEvent&)event).GetNumberOfPrimaries();
370  for (Int_t iTr = 0; iTr < nTracks; iTr++) {
371  AliMCParticle* particle =
372  static_cast<AliMCParticle*>(event.GetTrack(iTr));
373 
374  // Check if this charged and a primary
375  if (particle->Charge() == 0) continue;
376 
377  Bool_t isPrimary = event.IsPhysicalPrimary(iTr) && iTr < nPrim;
378 
379  // Fill 'dn/deta' histogram
380  if (isPrimary && primary) {
381  Double_t w = 1;
382  if (fTruthWeights) w = CalculateTruthWeight(particle);
383  primary->Fill(particle->Eta(), particle->Phi(), w);
384  }
385 
386  // Bail out if we're only processing primaries - perhaps we should
387  // track back to the original primary?
388  if (fUseOnlyPrimary && !isPrimary) continue;
389 
390  const AliMCParticle* mother = isPrimary ? particle : GetMother(iTr, event);
391  // IF the track corresponds to a primary, pass that as both
392  // arguments.
393  ProcessTrack(particle, mother);
394 
395  } // Loop over tracks
396  return kTRUE;
397 }
398 
399 
400 //____________________________________________________________________
401 Double_t
403  Bool_t isPrimary) const
404 {
405  // Note, it is always the ultimate mother that is passed here.
406  if (!p) return 0;
407  return fWeights->CalcWeight(p, isPrimary, fPhiR, fB);
408 }
409 //____________________________________________________________________
410 Double_t
411 AliBaseMCTrackDensity::CalculateTruthWeight(const AliMCParticle* p) const
412 {
413  // Note, it is always the ultimate mother that is passed here.
414  return fTruthWeights->CalcWeight(p, true, fPhiR, fB);
415 }
416 
417 #define PF(N,V,...) \
418  AliForwardUtil::PrintField(N,V, ## __VA_ARGS__)
419 #define PFB(N,FLAG) \
420  do { \
421  AliForwardUtil::PrintName(N); \
422  std::cout << std::boolalpha << (FLAG) << std::noboolalpha << std::endl; \
423  } while(false)
424 #define PFV(N,VALUE) \
425  do { \
426  AliForwardUtil::PrintName(N); \
427  std::cout << (VALUE) << std::endl; } while(false)
428 
429 //____________________________________________________________________
430 void
432 {
434  gROOT->IncreaseDirLevel();
435  PFB("Only primary tracks", fUseOnlyPrimary);
436  PFB("Use weights", fWeights);
437  if (fWeights) fWeights->Print(option);
438  gROOT->DecreaseDirLevel();
439 }
440 
441 //____________________________________________________________________
442 //
443 // EOF
444 //
void SetUseFlowWeights(Bool_t use)
double Double_t
Definition: External.C:58
void SetTruthWeights(AliBaseMCWeights *weights)
virtual void Print(Option_t *option="") const
const AliMCParticle * GetMother(Int_t iTr, const AliMCEvent &event) const
virtual void CreateOutputObjects(TList *list)
void SetWeights(AliBaseMCWeights *weights)
virtual Bool_t CheckTrackRef(AliTrackReference *) const
virtual Double_t StoreParticle(AliMCParticle *particle, const AliMCParticle *mother, AliTrackReference *ref) const
virtual void Print(Option_t *option="") const
Double_t GetTrackRefTheta(const AliTrackReference *ref) const
int Int_t
Definition: External.C:63
Various utilities used in PWGLF/FORWARD.
virtual AliTrackReference * ProcessRef(AliMCParticle *particle, const AliMCParticle *mother, AliTrackReference *ref)=0
Definition: External.C:228
Definition: External.C:212
Double_t CalculateWeight(const AliMCParticle *p, Bool_t isPrimary) const
AliBaseMCWeights * fWeights
AliBaseMCTrackDensity & operator=(const AliBaseMCTrackDensity &o)
virtual void Init(TList *l)
virtual Int_t GetDetectorId() const =0
#define DGUARD(L, N, F,...)
virtual void EndTrackRefs(Int_t)
static void PrintTask(const TObject &o)
#define PFB(N, FLAG)
Double_t CalculateTruthWeight(const AliMCParticle *p) const
const char Option_t
Definition: External.C:48
void test(int runnumber=195345)
Bool_t ProcessTracks(const AliMCEvent &event, const TVector3 &ip, TH2D *primary)
bool Bool_t
Definition: External.C:53
AliBaseMCWeights * fTruthWeights
virtual Double_t CalcWeight(const AliMCParticle *p, Bool_t isPrimary, Double_t phiR, Double_t b) const
Bool_t GetCollisionParameters(const AliMCEvent &event)
Bool_t ProcessTrack(AliMCParticle *particle, const AliMCParticle *mother)