AliPhysics  b76e98e (b76e98e)
AliFMDMCTrackELoss.cxx
Go to the documentation of this file.
1 #include "AliFMDMCTrackELoss.h"
2 #include "AliESDFMD.h"
3 #include "AliTrackReference.h"
4 #include <TMath.h>
5 #include "AliFMDStripIndex.h"
6 #include "AliFMDEncodedEdx.h"
7 #include "AliLog.h"
8 #include "AliMCParticle.h"
9 #include <TH2D.h>
10 #include <TH1D.h>
11 #include <TList.h>
12 #include <TROOT.h>
13 #include <TTree.h>
14 #include <TLorentzVector.h>
15 #include <TParticle.h>
16 #include <TClonesArray.h>
17 #include <iostream>
18 
19 //====================================================================
20 void
22 {
23  angle = 0;
24  oldDetector = 0;
25  oldRing = '\0';
26  oldSector = 1024;
27  oldStrip = 1024;
28  startStrip = 1024;
29  nRefs = 0;
30  nStrips = 0;
31  longest = 0x0;
32  de = 0;
33  dx = 0;
34  if (alsoCount) count = 0;
35 }
36 
37 //____________________________________________________________________
40 {
41  if (&o == this) return *this;
42  angle = o.angle;
44  oldRing = o.oldRing;
45  oldSector = o.oldSector;
46  oldStrip = o.oldStrip;
48  nRefs = o.nRefs;
49  nStrips = o.nStrips;
50  count = o.count;
51  longest = o.longest;
52  de = o.de;
53  dx = o.dx;
54  return *this;
55 }
56 
57 //====================================================================
59  : fGamma(0), fBeta(0), fEta(0), fDe(0), fDx(0), fDetId(0), fPdg(0),
60  fPrimary(true), fDetector(0), fRing('\0'),
61  fSector(0xFFFF), fStrip(0xFFFF)
62 {}
63 //____________________________________________________________________
64 void
66 {
67  if (fDetector > 0) return;
69 }
70 //____________________________________________________________________
71 UInt_t
73 {
74  return TMath::Abs(fPdg);
75 }
76 
77 //====================================================================
80  fState(),
81  fEvent(),
83  fUseTree(false),
84  fHits(0),
85  fTree(0),
86  fNr(0),
87  fNt(0),
88  fNc(0),
89  fNcr(0),
90  fBetaGammadEdx(0),
91  fBetaGammaEta(0),
92  fDEdxEta(0),
93  fPrimaries(0),
94  fSecondaries(0),
95  fAll(0),
96  fEta(0)
97 {
98  // Default constructor
99 }
100 
101 //____________________________________________________________________
103  : AliBaseMCTrackDensity("fmdMCTrackELoss"),
104  fState(),
105  fEvent(),
107  fUseTree(false),
108  fHits(0),
109  fTree(0),
110  fNr(0),
111  fNt(0),
112  fNc(0),
113  fNcr(0),
114  fBetaGammadEdx(0),
115  fBetaGammaEta(0),
116  fDEdxEta(0),
117  fPrimaries(0),
118  fSecondaries(0),
119  fAll(0),
120  fEta(0)
121 {
122  // Normal constructor constructor
123  SetTitle("mcTrackELoss");
124  fHits = new TClonesArray("AliFMDMCTrackELoss::Hit");
125 }
126 
127 #if 0
128 //____________________________________________________________________
131  fState(o.fState),
133  fNr(o.fNr),
134  fNt(o.fNt),
135  fNc(o.fNc),
136  fNcr(o.fNcr),
137  fOutput(o.fOutput)
138 {
139  // Normal constructor constructor
140 }
141 
142 //____________________________________________________________________
145 {
146  // Assignment operator
147  if (&o == this) return *this;
150  fUseTree = o.fUseTree;
151  fNr = o.fNr;
152  fNt = o.fNt;
153  fNc = o.fNc;
154  fNcr = o.fNcr;
155  fState = o.fState;
156 
157  if (o.fTree)
158  fTree = static_cast<TTree*>(o.fTree->Clone());
159 
160  return *this;
161 }
162 #endif
163 
164 //____________________________________________________________________
165 void
167 {
169  TList* ll = static_cast<TList*>(l->FindObject(GetTitle()));
170  if (!ll) ll = l;
171 
172  fNr = new TH1D("clusterRefs", "# track references per cluster",
173  21, -.5, 20.5);
174  fNr->SetXTitle("# of track references/cluster");
175  fNr->SetFillColor(kRed+1);
176  fNr->SetFillStyle(3001);
177  fNr->SetDirectory(0);
178  ll->Add(fNr);
179 
180  fNt = new TH1D("clusterSize", "cluster length in strips", 21, -.5, 20.5);
181  fNt->SetXTitle("Cluster size [strips]");
182  fNt->SetFillColor(kBlue+1);
183  fNt->SetFillStyle(3001);
184  fNt->SetDirectory(0);
185  ll->Add(fNt);
186 
187  fNc = new TH1D("nClusters", "# clusters per track", 21, -.5, 20.5);
188  fNc->SetXTitle("# clusters");
189  fNc->SetFillColor(kGreen+1);
190  fNc->SetFillStyle(3001);
191  fNc->SetDirectory(0);
192  ll->Add(fNc);
193 
194  fNcr = new TH2D("clusterVsRefs", "# clusters vs. # refs",
195  21, -.5, 20.5, 21, -.5, 20.5);
196  fNcr->SetXTitle("# References");
197  fNcr->SetYTitle("# Clusters");
198  fNcr->SetOption("COLZ");
199  fNcr->SetDirectory(0);
200  ll->Add(fNcr);
201 
202 
203  TArrayD edepArray;
204  TArrayD bgArray(601);
206  AliForwardUtil::MakeLogScale(600, -2, 5, bgArray);
207  for (Int_t i = 0; i < edepArray.GetSize(); i++)
208  edepArray[i] *= 10;
209 
210  fBetaGammadEdx = new TH2D("betaGammadEdx", "Energy loss",
211  bgArray.GetSize()-1, bgArray.GetArray(),
212  edepArray.GetSize()-1, edepArray.GetArray());
213  fBetaGammadEdx->SetDirectory(0);
214  fBetaGammadEdx->SetXTitle("#it{#beta#gamma}");
215  fBetaGammadEdx->SetYTitle("d#it{#Delta}/d#it{x}");
216  fBetaGammadEdx->Sumw2();
217  ll->Add(fBetaGammadEdx);
218 
219  fBetaGammaEta = new TH2D("betaGammaEta", "#beta#gamma",
220  200, -4, 6, bgArray.GetSize()-1, bgArray.GetArray());
221  fBetaGammaEta->SetXTitle("#eta");
222  fBetaGammaEta->SetYTitle("#it{#beta#gamma}");
223  fBetaGammaEta->Sumw2();
224  ll->Add(fBetaGammaEta);
225 
226  fDEdxEta = new TH2D("dEdxEta", "d#it{#Delta}/d#it{x}",
227  200, -4, 6, edepArray.GetSize()-1, edepArray.GetArray());
228  fDEdxEta->SetXTitle("#eta");
229  fDEdxEta->SetYTitle("d#it{#Delta}/d#it{x}");
230  fDEdxEta->Sumw2();
231  ll->Add(fDEdxEta);
232 
233  if (!fUseTree) return;
234 
235  fTree = new TTree("tree", "Tree of hits");
236  fTree->Branch("event", &fEvent, "ipZ/D:cent");
237  fTree->Branch("hits", &fHits);
238 }
239 
240 //____________________________________________________________________
241 void
243 {
244  fPrimaries.Reset(0);
245  fSecondaries.Reset(0);
246  fAll.Reset(0);
247  fEta.Reset(AliESDFMD::kInvalidEta);
248  fHits->Clear();
249 }
250 
251 //____________________________________________________________________
252 Int_t
254 {
255  return AliTrackReference::kFMD;
256 }
257 
258 //____________________________________________________________________
259 void
261 {
262  fState.Clear(true);
263 }
264 
265 //____________________________________________________________________
266 void
268 {
269  fNc->Fill(fState.count);
270  fNcr->Fill(nRefs, fState.count);
271  fState.Clear(true);
272 }
273 
274 //____________________________________________________________________
275 AliTrackReference*
276 AliFMDMCTrackELoss::ProcessRef(AliMCParticle* particle,
277  const AliMCParticle* mother,
278  AliTrackReference* ref)
279 {
280  // Get the detector coordinates
281  UShort_t d, s, t;
282  Char_t r;
283  AliFMDStripIndex::Unpack(ref->UserId(), d, r, s, t);
284  Double_t edep, length, dEdep, dLength;
285  AliFMDEncodedEdx::Decode((ref->UserId() >> 19), edep, length, dEdep, dLength);
286 
287 
288 
289 
290  // Calculate distance of previous reference to base of cluster
291  UShort_t nT = TMath::Abs(t - fState.startStrip) + 1;
292 
293  // Now check if we should flush to output
294  Bool_t used = false;
295 
296  // If this is a new detector/ring, then reset the other one
297  // Check if we have a valid old detectorm ring, and sector
298  if (fState.oldDetector > 0 &&
299  fState.oldRing != '\0' &&
300  fState.oldSector != 1024) {
301  // New detector, new ring, or new sector
302  if (d != fState.oldDetector ||
303  r != fState.oldRing ||
304  s != fState.oldSector) {
305  if (fDebug) Info("Process", "New because new sector");
306  used = true;
307  }
308  else if (nT > fMaxConsequtiveStrips) {
309  if (fDebug) Info("Process", "New because too long: %d (%d,%d,%d)",
311  used = true;
312  }
313  }
314  if (used) {
315  if (fDebug)
316  Info("Process", "I=%p L=%p D=%d (was %d), R=%c (was %c), "
317  "S=%2d (was %2d) t=%3d (was %3d) nT=%3d/%4d",
318  ref, fState.longest,
319  d, fState.oldDetector,
320  r, fState.oldRing,
321  s, fState.oldSector,
322  t, fState.oldStrip,
324  // Int_t nnT = TMath::Abs(fState.oldStrip - fState.startStrip) + 1;
325  StoreParticle(particle, mother, fState.longest);
326  fState.Clear(false);
327  }
328 
329  // If base of cluster not set, set it here.
330  if (fState.startStrip == 1024) fState.startStrip = t;
331 
332  // Calculate distance of previous reference to base of cluster
333  fState.nStrips = TMath::Abs(t - fState.startStrip) + 1;
334 
335  // Count number of track refs in this sector
336  fState.nRefs++;
337 
338  fState.oldDetector = d;
339  fState.oldRing = r;
340  fState.oldSector = s;
341  fState.oldStrip = t;
342  fState.de += edep;
343  fState.dx += length;
344 
345  // Debug output
346  if (fDebug) {
347  if (t == fState.startStrip)
348  Info("Process", "New cluster starting at FMD%d%c[%2d,%3d]",
349  d, r, s, t);
350  else
351  Info("Process", "Adding to cluster starting at FMD%d%c[%2d,%3d], "
352  "length=%3d (now in %3d, previous %3d)",
354  }
355 
356  // The longest passage is determined through the angle
357  Double_t ang = GetTrackRefTheta(ref);
358  if (ang > fState.angle) {
359  fState.longest = ref;
360  fState.angle = ang;
361  }
362 
363  return fState.longest;
364 }
365 
366 //____________________________________________________________________
367 Double_t
368 AliFMDMCTrackELoss::StoreParticle(AliMCParticle* particle,
369  const AliMCParticle* mother,
370  AliTrackReference* ref) const
371 {
372  Double_t w =
373  AliBaseMCTrackDensity::StoreParticle(particle, mother, ref);
374  if (w <= 0) return w;
375 
376  // Get the detector coordinates
377  UShort_t d, s, t;
378  Char_t r;
379  AliFMDStripIndex::Unpack(ref->UserId(), d, r, s, t);
380 
381  Double_t de = fState.de;
382  Double_t dx = fState.dx;
383  Bool_t prim = particle == mother;
384  fAll(d,r,s,t) += de;
385  if (prim) fPrimaries(d,r,s,t) += de;
386  else fSecondaries(d,r,s,t) += de;
387 
388  // Fill histograms
389  fNr->Fill(fState.nRefs);
390  fNt->Fill(fState.nStrips);
391 
392  fState.count++;
393 
394  if (de <= 0 || dx <= 0) return w;
395 
396  TLorentzVector v;
397  particle->Particle()->Momentum(v);
398  if (v.E() <= 0) return w;
399  if (v.Beta() < 0 || v.Beta() > 1) return w;
400 
401  Double_t eta = fEta(d,r,s,t);
402  Double_t beta = v.Beta();
403  Double_t gamma = v.Gamma();
404  Double_t dEdx = de/dx;
405 
406 
407  fBetaGammadEdx->Fill(beta*gamma, dEdx);
408  fBetaGammaEta->Fill(eta, beta*gamma);
409  fDEdxEta->Fill(eta, dEdx);
410 
411  Int_t nHits = fHits->GetEntries();
412  Hit* hit = new((*fHits)[nHits]) Hit;
413  hit->fPrimary = prim;
414  hit->fEta = eta;
415  hit->fDe = de;
416  hit->fDx = dx;
417  hit->fDetId = ref->UserId() & AliFMDStripIndex::kIdMask;
418  hit->fBeta = beta;
419  hit->fGamma = gamma;
420  hit->fPdg = particle->PdgCode();
421 
422  return w;
423 }
424 
425 //____________________________________________________________________
426 Bool_t
428  const AliMCEvent& event,
429  const TVector3& ip,
430  Double_t cent)
431 {
432  //
433  // Filter the input kinematics and track references, using
434  // some of the ESD information
435  //
436  // Parameters:
437  // input Input ESD event
438  // event Input MC event
439  // vz Vertex position
440  // output Output ESD-like object
441  // primary Per-event histogram of primaries
442  //
443  // Return:
444  // True on succes, false otherwise
445  //
446  Clear();
447  fEvent.fCent = cent;
448  fEvent.fIpZ = ip.Z();
449 
450  // Copy eta values to output
451  for (UShort_t ed = 1; ed <= 3; ed++) {
452  UShort_t nq = (ed == 1 ? 1 : 2);
453  for (UShort_t eq = 0; eq < nq; eq++) {
454  Char_t er = (eq == 0 ? 'I' : 'O');
455  UShort_t ns = (eq == 0 ? 20 : 40);
456  UShort_t nt = (eq == 0 ? 512 : 256);
457  for (UShort_t es = 0; es < ns; es++)
458  for (UShort_t et = 0; et < nt; et++)
459  fEta(ed, er, es, et) = input.Eta(ed, er, es, et);
460  }
461  }
462 
463  Bool_t ret = ProcessTracks(event, ip, 0);
464 
465  if (fTree) fTree->Fill();
466 
467  return ret;
468 }
469 
470 #define PFV(N,VALUE) \
471  do { \
472  AliForwardUtil::PrintName(N); \
473  std::cout << (VALUE) << std::endl; } while(false)
474 //____________________________________________________________________
475 void
477 {
479  gROOT->IncreaseDirLevel();
480  PFV("Max cluster size", fMaxConsequtiveStrips);
481  gROOT->DecreaseDirLevel();
482 }
483 
484 //____________________________________________________________________
485 //
486 // EOF
487 //
double Double_t
Definition: External.C:58
AliFMDFloatMap fPrimaries
Double_t StoreParticle(AliMCParticle *particle, const AliMCParticle *mother, AliTrackReference *ref) const
void Clear(Bool_t alsoCount=false)
Longest track through.
virtual void CreateOutputObjects(TList *list)
AliTrackReference * longest
char Char_t
Definition: External.C:18
void Print(Option_t *option="") const
static const Spec & GetdEAxis()
struct AliFMDMCTrackELoss::State fState
State & operator=(const State &o)
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
unsigned int UInt_t
Definition: External.C:33
Definition: External.C:228
Definition: External.C:212
static void Decode(UInt_t bits, Double_t &edep, Double_t &length)
void FillBinArray(TArrayD &a) const
AliFMDMCTrackELoss & operator=(const AliFMDMCTrackELoss &o)
Class to encode the energy loss and path length of a particle into track reference bits...
AliBaseMCTrackDensity & operator=(const AliBaseMCTrackDensity &o)
void Clear(Option_t *opt="")
Bool_t Calculate(const AliESDFMD &esd, const AliMCEvent &event, const TVector3 &ip, Double_t cent)
void EndTrackRefs(Int_t nRefs)
static void MakeLogScale(Int_t nBins, Int_t minOrder, Int_t maxOrder, TArrayD &bins)
unsigned short UShort_t
Definition: External.C:28
const char Option_t
Definition: External.C:48
Int_t GetDetectorId() const
Bool_t ProcessTracks(const AliMCEvent &event, const TVector3 &ip, TH2D *primary)
bool Bool_t
Definition: External.C:53
AliFMDFloatMap fSecondaries
void CreateOutputObjects(TList *list)
#define PFV(N, VALUE)
AliTrackReference * ProcessRef(AliMCParticle *particle, const AliMCParticle *mother, AliTrackReference *ref)
static void Unpack(UInt_t id, UShort_t &det, Char_t &rng, UShort_t &sec, UShort_t &str)