AliPhysics  d497547 (d497547)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliFMDMCHitEnergyFitter.cxx
Go to the documentation of this file.
2 #include "AliMCEvent.h"
3 #include "AliESDEvent.h"
4 #include "AliStack.h"
5 #include "AliFMDHit.h"
6 #include "AliESDFMD.h"
7 #include "AliMCAuxHandler.h"
8 #include "AliGenEventHeader.h"
9 #include "AliLog.h"
10 #include "AliFMDEncodedEdx.h"
11 #include <TH1.h>
12 #include <TH2.h>
13 #include <TNtuple.h>
14 #include <TMath.h>
15 
16 
17 //____________________________________________________________________
20  fSumPrimary(0),
21  fSumSecondary(0),
22  fIp(),
23  fNTrack(0),
24  fNPrimary(0),
25  fTuple(0)
26 {
27 }
28 
29 //____________________________________________________________________
31  Bool_t useTuple)
32  : AliFMDEnergyFitter(title),
33  fSumPrimary(0),
34  fSumSecondary(0),
35  fIp(),
36  fNTrack(0),
37  fNPrimary(0),
38  fTuple(0)
39 {
40  // Some defaults
41  fUseIncreasingBins = true;
42  fDoFits = true;
43  fDoMakeObject = false;
45 
46  if (!useTuple) return;
47 
48  fTuple = new TNtuple("hits", "Hits",
49  "primary:eta:phi:edep:length:"
50  "detector:ring:sector:strip:ipz");
51 
52 }
53 //____________________________________________________________________
55 {}
56 
57 //____________________________________________________________________
58 void
60 {
61  DGUARD(fDebug, 2, "Create output objects w/MC hits (%p)", dir);
63  // TList* d = static_cast<TList*>(l->FindObject(GetName()));
64 }
65 //____________________________________________________________________
66 Bool_t
67 AliFMDMCHitEnergyFitter::PreEvent(const AliMCEvent& mcInput)
68 {
69  DGUARD(fDebug, 5, "Reset for event w/MC hits (%p)", &mcInput);
70  fSumPrimary.Reset(0);
71  fSumSecondary.Reset(0);
72 
73  AliMCEvent& mc = const_cast<AliMCEvent&>(mcInput);
74  AliHeader* header = mc.Header();
75  AliStack* stack = mc.Stack();
76  AliGenEventHeader* genHeader = header->GenEventHeader();
77 
78  genHeader->PrimaryVertex(fIp);
79  fNTrack = stack->GetNtrack();
80  fNPrimary = stack->GetNprimary();
81 
82  return true;
83 }
84 
85 //____________________________________________________________________
86 Bool_t
88  const AliMCEvent& mcInput,
89  AliMCAuxHandler& handler)
90 {
91  DGUARD(fDebug, 5, "Process an event for MC hits (%p,%p,%p)",
92  &esdInput, &mcInput, &handler);
93  PreEvent(mcInput);
94 
95  Int_t nEntries = handler.GetNEntry();
96  DMSG(fDebug,5, "We got a total of %d particles", nEntries);
97 
98  for (Int_t i = 0; i < nEntries; i++) {
99  TClonesArray* hits = handler.GetEntryArray(i);
100  if (!hits) continue;
101 
102  AccumulateHits(mcInput, *hits);
103  }
104  return PostEvent(esdInput);
105 }
106 
107 //____________________________________________________________________
108 Bool_t
109 AliFMDMCHitEnergyFitter::AccumulateHits(const AliMCEvent& mcInput,
110  const TClonesArray& hits)
111 {
112  DGUARD(fDebug, 5, "Accumulate MC hit energy loss (%p,%p,%d)",
113  &mcInput, &hits, hits.GetEntries());
114 
115  AliStack* stack = const_cast<AliMCEvent&>(mcInput).Stack();
116  Float_t cache[10];
117  Int_t nHit = hits.GetEntries();
118  Int_t oldTr = -1;
119  UShort_t oldD = 0;
120  Char_t oldR = '\0';
121  for (Int_t j = 0; j < nHit; j++) {
122 
123  // Check the hit
124  AliFMDHit* hit = static_cast<AliFMDHit*>(hits.At(j));
125  if (!hit) continue;
126 
127  // Zero fill array
128  if (fTuple) for (Int_t k = 0; k < 10; k++) cache[k] = 0;
129 
130  // Get the track number
131  Int_t iTr = hit->Track();
132  if (iTr < 0 || iTr >= fNTrack)
133  AliWarningF("Track # %d seems out of bounds [0,%d]",
134  iTr, fNTrack-1);
135 
136  // Get the track
137  AliMCParticle* particle =
138  static_cast<AliMCParticle*>(mcInput.GetTrack(iTr));
139 
140  // Particle type 0: unknown, 1: primary, 2: secondary
141  // Int_t flag = 0; - not used at the moment
142 
143  // Check if this charged and a primary
144  if (particle) {
145  // Only charged tracks
146  if (!particle->Charge() != 0) continue;
147  }
148  Bool_t isPrimary = stack->IsPhysicalPrimary(iTr) && iTr < fNPrimary;
149  Double_t eloss = hit->Edep();
150  Double_t length = hit->Length();
151  // Double_t eta = particle->Eta();
152  Double_t x = hit->X() - fIp[0];
153  Double_t y = hit->Y() - fIp[1];
154  Double_t z = hit->Z() - fIp[2];
155  Double_t r = TMath::Sqrt(x*x + y*y);
156  Double_t phi = TMath::ATan2(y, x);
157  Double_t theta = TMath::ATan2(r, z);
158  Double_t eta = -TMath::Log(TMath::Tan(theta/2));
159  UShort_t det = hit->Detector();
160  Char_t rng = hit->Ring();
161  UShort_t sec = hit->Sector();
162  UShort_t str = hit->Strip();
163 
164  Bool_t isDistinct = (iTr != oldTr) || (det != oldD) || (rng != oldR);
165  if (isDistinct) {
166  RingHistos* h = static_cast<RingHistos*>(GetRingHistos(det, rng));
167  h->fKind->Fill(eta, isPrimary ? .5 : 1.5);
168  }
169  oldTr = iTr;
170  oldD = det;
171  oldR = rng;
172 
173  AliFMDFloatMap& sum = (isPrimary ? fSumPrimary : fSumSecondary);
174  // Float_t old = sum(det, rng, sec, str);
175  sum(det, rng, sec, str) += eloss;
176 #if 0
177  Info("", "FMD%d%c[%02d,%03d]-%s: Adding %10.7f to %10.7f -> %10.7f",
178  det, rng, sec, str, (isPrimary ? "1st" : "2nd"),
179  eloss, old, sum(det, rng, sec, str));
180 #endif
181  if (!fTuple) continue;
182 
183  // Leaves primary:eta:phi:edep:length:detector:ring:sector:strip:ipz
184  cache[0] = isPrimary ? 1 : 0;
185  cache[1] = eta;
186  cache[2] = phi;
187  cache[3] = eloss;
188  cache[4] = length;
189  cache[5] = hit->Detector();
190  cache[6] = Int_t(hit->Ring());
191  cache[7] = hit->Sector();
192  cache[8] = hit->Strip();
193  cache[9] = fIp[2];
194  fTuple->Fill(cache);
195  }
196  return true;
197 }
198 
199 //____________________________________________________________________
201 {
202  DGUARD(fDebug, 5, "Convert MC hit energy loss (%p)", &input);
203 
204  AliESDFMD* esdFMD = input.GetFMDData();
205  for (UShort_t d = 1; d <= 3; d++) {
206  UShort_t nQ = (d == 1 ? 1 : 2);
207  for (UShort_t q = 0; q < nQ; q++) {
208  UShort_t nS = (q == 0 ? 20 : 40);
209  UShort_t nT = (q == 0 ? 512 : 256);
210  Char_t r = (q == 0 ? 'I' : 'O');
211  RingHistos* h =
213  if (!h) continue;
214 
215  for (UShort_t s = 0; s < nS; s++) {
216  for (UShort_t t = 0; t < nT; t++) {
217  Float_t totPrim = fSumPrimary(d, r, s, t);
218  Float_t totSec = fSumSecondary(d, r, s, t);
219  Float_t totAll = totPrim + totSec;
220  Double_t esdEta = esdFMD->Eta(d, r, s, t);
221  if (totAll > 0) h->FillMC(0, esdEta, totAll);
222  if (totPrim > 0) h->FillMC(1, esdEta, totPrim);
223  if (totSec > 0) h->FillMC(2, esdEta, totSec);
224  }
225  }
226  }
227  }
228  return true;
229 }
230 
231 //____________________________________________________________________
234 {
235  // DGUARD(fDebug, 1, "Create Ring cache of MC hit energy loss (FMD%d%c)",d,r);
236  return new AliFMDMCHitEnergyFitter::RingHistos(d,r);
237 }
238 
239 //====================================================================
242  fPrimary(0),
243  fSecondary(0),
244  fKind(0)
245 {}
246 
247 //____________________________________________________________________
250  fPrimary(0),
251  fSecondary(0),
252  fKind(0)
253 {}
254 //____________________________________________________________________
255 TArrayD
257  Double_t,
258  Double_t) const
259 {
260  // Make an increasing axis for ELoss distributions.
261  //
262  // We use the service function of AliFMDEncodedEdx to do this
263  TArrayD ret;
265  s.FillBinArray(ret);
266 
267  return ret;
268 }
269 //____________________________________________________________________
270 void
272  const TAxis& /*cAxis*/,
273  Double_t maxDE,
274  Int_t nDEbins,
275  Bool_t useIncrBin)
276 {
277  // AliFMDEnergyFitter::RingHistos::SetupForData(eAxis, cAxis, maxDE, nDEbins,
278  // useIncrBin);
279 
280  fHist = Make("eloss", "#sum#Delta_{true} of all",
281  eAxis, maxDE, nDEbins, useIncrBin);
282  fPrimary = Make("primary", "#sum#Delta_{true} of primaries",
283  eAxis, maxDE, nDEbins, useIncrBin);
284  fSecondary = Make("secondary","#sum#Delta_{true} of secondaries",
285  eAxis, maxDE, nDEbins, useIncrBin);
286  fHist->SetXTitle("#Delta_{true}");
287  fSecondary->SetXTitle("#Delta_{true}");
288  fPrimary->SetXTitle("#Delta_{true}");
289 
290  fKind = 0;
291  if (eAxis.GetXbins()->GetArray())
292  fKind = new TH2I("kind", "Particle types",
293  eAxis.GetNbins(), eAxis.GetXbins()->GetArray(),
294  2, 0, 2);
295  else
296  fKind = new TH2I("kind", "Particle types",
297  eAxis.GetNbins(), eAxis.GetXmin(), eAxis.GetXmax(),
298  2, 0, 2);
299  fKind->SetXTitle("#eta");
300  fKind->GetYaxis()->SetBinLabel(1, "Primary");
301  fKind->GetYaxis()->SetBinLabel(2, "Secondary");
302  fKind->SetDirectory(0);
303  fKind->Sumw2();
304 
305  fList->Add(fHist);
306  fList->Add(fPrimary);
307  fList->Add(fSecondary);
308  fList->Add(fKind);
309 }
310 
311 //____________________________________________________________________
312 void
314  Double_t mult)
315 {
316  switch (flag) {
317  // case 0: AliFMDEnergyFitter::RingHistos::Fill(false, eta, 0, mult); break;
318  case 0: fHist->Fill(eta, mult); break;
319  case 1: fPrimary->Fill(eta, mult); break;
320  case 2: fSecondary->Fill(eta, mult); break;
321  }
322 }
323 
324 //____________________________________________________________________
325 TObjArray*
327  Double_t lowCut,
328  UShort_t nParticles,
329  UShort_t minEntries,
330  UShort_t minusBins,
331  Double_t relErrorCut,
332  Double_t chi2nuCut,
333  Double_t minWeight,
334  Double_t regCut,
335  EResidualMethod residuals) const
336 {
337  TObjArray* all = FitSlices(dir, "eloss", lowCut, nParticles, minEntries,
338  minusBins, relErrorCut, chi2nuCut, minWeight,
339  regCut, residuals, false, 0);
340  TObjArray* prim = FitSlices(dir, "primary", lowCut, nParticles, minEntries,
341  minusBins, relErrorCut, chi2nuCut, minWeight,
342  regCut, residuals, false, 0);
343  TObjArray* sec = FitSlices(dir, "secondary", lowCut, nParticles,
344  minEntries, minusBins, relErrorCut, chi2nuCut,
345  minWeight, regCut, residuals, false, 0);
346  if (!all || !prim || !sec) {
347  AliWarningF("Failed to get results for all(%p), primary(%p), and/or "
348  "secondary(%p)", all, prim, sec);
349  return 0;
350  }
351  // Already added to the sub-folders
352  // dir->Add(all);
353  // dir->Add(prim);
354  // dir->Add(sec);
355 
356  Int_t nPrim = prim->GetEntriesFast();
357  TObjArray* out = new TObjArray;
358  for (Int_t i = 0; i < nPrim; i++) {
359  TH1* h = static_cast<TH1*>(prim->At(i));
360  if (!h) continue;
361 
362  TAxis* xAxis = h->GetXaxis();
363  TH2* hh = 0;
364  if (xAxis->GetXbins()->GetArray())
365  hh = new TH2D(h->GetName(), h->GetTitle(),
366  xAxis->GetNbins(), xAxis->GetXbins()->GetArray(),
367  3, 0, 3);
368  else
369  hh = new TH2D(h->GetName(), h->GetTitle(),
370  xAxis->GetNbins(), xAxis->GetXmin(), xAxis->GetXmax(),
371  3, 0, 3);
372  hh->GetYaxis()->SetBinLabel(1, "Primary");
373  hh->GetYaxis()->SetBinLabel(2, "Secondary");
374  hh->GetYaxis()->SetBinLabel(3, "All");
375  hh->GetXaxis()->SetTitle("#eta");
376  out->Add(hh);
377  }
378  for (Int_t i = 0; i < nPrim; i++) {
379  TH2* h = static_cast<TH2*>(out->At(i));
380  if (!h) continue;
381 
382  TH1* hp = static_cast<TH1*>(prim->At(i));
383  TH1* hs = static_cast<TH1*>(sec->At(i));
384  TH1* ha = static_cast<TH1*>(all->At(i));
385  TH1* hh[] = { hp, hs, ha };
386  for (Int_t j = 0; j < 3; j++) {
387  TH1* ph = hh[j];
388  if (!ph) continue;
389 
390  for (Int_t k = 1; k <= ph->GetNbinsX(); k++) {
391  Double_t c = ph->GetBinContent(k);
392  Double_t e = ph->GetBinError(k);
393  h->SetBinContent(k, j+1, c);
394  h->SetBinError(k, j+1, e);
395  }
396  }
397  }
398  TList* l = GetOutputList(dir);
399  if (!l) return 0;
400 
401  out->SetName("compartive");
402  out->SetOwner();
403  l->Add(out);
404  return out;
405 }
virtual Bool_t AccumulateHits(const AliMCEvent &mcInput, const TClonesArray &hits)
EResidualMethod fResidualMethod
double Double_t
Definition: External.C:58
virtual void SetupForData(const TAxis &eAxis, const TAxis &cAxis, Double_t maxDE=10, Int_t nDEbins=300, Bool_t useIncrBin=true)
const char * title
Definition: MakeQAPdf.C:27
virtual Bool_t PostEvent(const AliESDEvent &esdInput)
virtual void CreateOutputObjects(TList *dir)
char Char_t
Definition: External.C:18
static const Spec & GetdEAxis()
#define DMSG(L, N, F,...)
TCanvas * c
Definition: TestFitELoss.C:172
virtual void FillMC(UShort_t flag, Double_t eta, Double_t mult)
AliStack * stack
AliFMDEnergyFitter::RingHistos * CreateRingHistos(UShort_t d, Char_t r) const
TArrayD MakeIncreasingAxis(Int_t nBins, Double_t low, Double_t high) const
virtual void CreateOutputObjects(TList *dir)
int Int_t
Definition: External.C:63
float Float_t
Definition: External.C:68
virtual Bool_t PreEvent(const AliMCEvent &mcInput)
Definition: External.C:228
void FillBinArray(TArrayD &a) const
Class to encode the energy loss and path length of a particle into track reference bits...
#define DGUARD(L, N, F,...)
TH1 * Make(UShort_t d, Char_t r)
Definition: TestEtaPhi.C:3
Int_t GetNEntry() const
Definition: External.C:220
unsigned short UShort_t
Definition: External.C:28
RingHistos * GetRingHistos(UShort_t d, Char_t r) const
TObjArray * Fit(TList *dir, Double_t lowCut, UShort_t nParticles, UShort_t minEntries, UShort_t minusBins, Double_t relErrorCut, Double_t chi2nuCut, Double_t minWeight, Double_t regCut, EResidualMethod residuals) const
bool Bool_t
Definition: External.C:53
Definition: External.C:196
TDirectoryFile * dir
virtual Bool_t Event(const AliESDEvent &esdInput, const AliMCEvent &mcInput, AliMCAuxHandler &handler)
TClonesArray * GetEntryArray(Int_t entry)