AliPhysics  64f4410 (64f4410)
AliFMDMCTrackInspector.cxx
Go to the documentation of this file.
2 #include "AliMCEvent.h"
3 #include "AliESDEvent.h"
4 #include "AliStack.h"
5 #include "AliESDFMD.h"
6 #include "AliGenEventHeader.h"
7 #include "AliHeader.h"
8 #include "AliLog.h"
9 #include "AliFMDEncodedEdx.h"
10 #include <TH1.h>
11 #include <TH2.h>
12 #include <TMath.h>
13 #include <TClonesArray.h>
14 
15 //____________________________________________________________________
18  fTracker(),
19  fIp(),
20  fNTrack(0),
21  fNPrimary(0)
22 {
23 }
24 
25 //____________________________________________________________________
27  : AliFMDEnergyFitter(title),
28  fTracker("tracker"),
29  fIp(),
30  fNTrack(0),
31  fNPrimary(0)
32 {
33  // Some defaults
34  fUseIncreasingBins = true;
35  fDoFits = true;
36  fDoMakeObject = false;
38 }
39 //____________________________________________________________________
41 {}
42 
43 //____________________________________________________________________
44 void
46 {
47  DGUARD(fDebug, 2, "Create output objects w/MC hits (%p)", dir);
50 
51  TIter next(&fRingHistos);
52  RingHistos* o = 0;
53  while ((o = static_cast<RingHistos*>(next()))) {
54  o->fBetaGammadEdx = static_cast<TH2*>(fTracker.GetBetaGammadEdx()->Clone());
55  o->fBetaGammaEta = static_cast<TH2*>(fTracker.GetBetaGammaEta()->Clone());
56  o->fDedxEta = static_cast<TH2*>(fTracker.GetDEdxEta()->Clone());
57  o->fBetaGammadEdx->SetDirectory(0);
58  o->fBetaGammaEta ->SetDirectory(0);
59  o->fDedxEta ->SetDirectory(0);
60  o->fList->Add(o->fBetaGammadEdx);
61  o->fList->Add(o->fBetaGammaEta );
62  o->fList->Add(o->fDedxEta );
63  }
64 
65  // TList* d = static_cast<TList*>(l->FindObject(GetName()));
66 }
67 //____________________________________________________________________
68 Bool_t
69 AliFMDMCTrackInspector::PreEvent(const AliMCEvent& mcInput)
70 {
71  DGUARD(fDebug, 5, "Reset for event w/MC hits (%p)", &mcInput);
72  AliMCEvent& mc = const_cast<AliMCEvent&>(mcInput);
73  AliHeader* header = mc.Header();
74  AliStack* stack = mc.Stack();
75  AliGenEventHeader* genHeader = header->GenEventHeader();
76 
77  genHeader->PrimaryVertex(fIp);
78  fNTrack = stack->GetNtrack();
79  fNPrimary = stack->GetNprimary();
80 
81  return true;
82 }
83 
84 //____________________________________________________________________
85 Bool_t
87  const AliMCEvent& mcInput,
88  Double_t cent)
89 {
90  DGUARD(fDebug, 5, "Process an event for MC hits (%p,%p)",
91  &esdInput, &mcInput);
92  PreEvent(mcInput);
93 
94  AliESDFMD* esdFMD = esdInput.GetFMDData();
95  if (!esdFMD) return true;
96 
97  TVector3 ip(fIp[0], fIp[1], fIp[2]);
98  fTracker.Calculate(*esdFMD, mcInput, ip, cent);
99 
100  return PostEvent();
101 }
102 
103 
104 //____________________________________________________________________
106 {
107  DGUARD(fDebug, 5, "Fill MC hit energy loss");
108 
109  // AliESDFMD* esdFMD = input.GetFMDData();
110  for (UShort_t d = 1; d <= 3; d++) {
111  UShort_t nQ = (d == 1 ? 1 : 2);
112  for (UShort_t q = 0; q < nQ; q++) {
113  UShort_t nS = (q == 0 ? 20 : 40);
114  UShort_t nT = (q == 0 ? 512 : 256);
115  Char_t r = (q == 0 ? 'I' : 'O');
116  RingHistos* h =
118  if (!h) continue;
119 
120  for (UShort_t s = 0; s < nS; s++) {
121  for (UShort_t t = 0; t < nT; t++) {
122  Float_t totPrim = fTracker.GetPrimaries()(d, r, s, t);
123  Float_t totSec = fTracker.GetSecondaries()(d, r, s, t);
124  Float_t totAll = fTracker.GetAll()(d, r, s, t);
125  Double_t esdEta = fTracker.GetEta()(d, r, s, t);
126  if (totAll > 0) h->FillMC(0, esdEta, totAll);
127  if (totPrim > 0) h->FillMC(1, esdEta, totPrim);
128  if (totSec > 0) h->FillMC(2, esdEta, totSec);
129  }
130  }
131  }
132  }
133  TClonesArray* hits = fTracker.GetHits();
134  TIter next(hits);
135  AliFMDMCTrackELoss::Hit* hit = 0;
136  while ((hit = static_cast<AliFMDMCTrackELoss::Hit*>(next()))) {
137  RingHistos* h =
139  hit->R()));
140  if (!h) continue;
141  h->fBetaGammadEdx->Fill(hit->BetaGamma(), hit->DeDx());
142  h->fBetaGammaEta->Fill(hit->Eta(), hit->BetaGamma());
143  h->fDedxEta->Fill(hit->Eta(), hit->DeDx());
144 
145  }
146  return true;
147 }
148 
149 //____________________________________________________________________
152 {
153  // DGUARD(fDebug, 1, "Create Ring cache of MC hit energy loss (FMD%d%c)",d,r);
154  return new AliFMDMCTrackInspector::RingHistos(d,r);
155 }
156 
157 //====================================================================
160  fPrimary(0),
161  fSecondary(0),
162  fBetaGammadEdx(0),
163  fBetaGammaEta(0),
164  fDedxEta(0)
165 {}
166 
167 //____________________________________________________________________
170  fPrimary(0),
171  fSecondary(0),
172  fBetaGammadEdx(0),
173  fBetaGammaEta(0),
174  fDedxEta(0)
175 {}
176 //____________________________________________________________________
177 TArrayD
179  Double_t,
180  Double_t) const
181 {
182  // Make an increasing axis for ELoss distributions.
183  //
184  // We use the service function of AliFMDEncodedEdx to do this
185  TArrayD ret;
187  s.FillBinArray(ret);
188 
189  return ret;
190 }
191 //____________________________________________________________________
192 void
194  const TAxis& /*cAxis*/,
195  Double_t maxDE,
196  Int_t nDEbins,
197  Bool_t useIncrBin)
198 {
199  // AliFMDEnergyFitter::RingHistos::SetupForData(eAxis, cAxis, maxDE, nDEbins,
200  // useIncrBin);
201 
202  fHist = Make("eloss", "#sum#Delta_{true} of all",
203  eAxis, maxDE, nDEbins, useIncrBin);
204  fPrimary = Make("primary", "#sum#Delta_{true} of primaries",
205  eAxis, maxDE, nDEbins, useIncrBin);
206  fPrimary->SetMarkerStyle(24);
207  fPrimary->SetMarkerSize(fPrimary->GetMarkerSize()*1.2);
208  fSecondary = Make("secondary","#sum#Delta_{true} of secondaries",
209  eAxis, maxDE, nDEbins, useIncrBin);
210  fSecondary->SetMarkerStyle(25);
211  fSecondary->SetMarkerSize(fSecondary->GetMarkerSize()*1.2);
212  fHist->SetXTitle("#Delta_{true}");
213  fSecondary->SetXTitle("#Delta_{true}");
214  fPrimary->SetXTitle("#Delta_{true}");
215 
216  fList->Add(fHist);
217  fList->Add(fPrimary);
218  fList->Add(fSecondary);
219 }
220 
221 //____________________________________________________________________
222 void
224  Double_t mult)
225 {
226  switch (flag) {
227  // case 0: AliFMDEnergyFitter::RingHistos::Fill(false, eta, 0, mult); break;
228  case 0: fHist->Fill(eta, mult); break;
229  case 1: fPrimary->Fill(eta, mult); break;
230  case 2: fSecondary->Fill(eta, mult); break;
231  }
232 }
233 
234 //____________________________________________________________________
235 void
237 {
238  // First scale to bin width
240  // Then smoothen the histogram
241  dist->Smooth(2);
242 }
243 
244 //____________________________________________________________________
245 TObjArray*
247  Double_t lowCut,
248  UShort_t nParticles,
249  UShort_t minEntries,
250  UShort_t minusBins,
251  Double_t relErrorCut,
252  Double_t chi2nuCut,
253  Double_t minWeight,
254  Double_t regCut,
255  EResidualMethod residuals) const
256 {
257  TObjArray* all = FitSlices(dir, "eloss", lowCut, nParticles, minEntries,
258  minusBins, relErrorCut, chi2nuCut, minWeight,
259  regCut, residuals, false, 0);
260  TObjArray* prim = FitSlices(dir, "primary", lowCut, nParticles, minEntries,
261  minusBins, relErrorCut, chi2nuCut, minWeight,
262  regCut, residuals, false, 0);
263  TObjArray* sec = FitSlices(dir, "secondary", lowCut, nParticles,
264  minEntries, minusBins, relErrorCut, chi2nuCut,
265  minWeight, regCut, residuals, false, 0);
266  if (!all || !prim || !sec) {
267  AliWarningF("Failed to get results for all(%p), primary(%p), and/or "
268  "secondary(%p)", all, prim, sec);
269  return 0;
270  }
271  // Already added to the sub-folders
272  // dir->Add(all);
273  // dir->Add(prim);
274  // dir->Add(sec);
275 
276  Int_t nPrim = prim->GetEntriesFast();
277  TObjArray* out = new TObjArray;
278  for (Int_t i = 0; i < nPrim; i++) {
279  TH1* h = static_cast<TH1*>(prim->At(i));
280  if (!h) continue;
281 
282  TAxis* xAxis = h->GetXaxis();
283  TH2* hh = 0;
284  if (xAxis->GetXbins()->GetArray())
285  hh = new TH2D(h->GetName(), h->GetTitle(),
286  xAxis->GetNbins(), xAxis->GetXbins()->GetArray(),
287  3, 0, 3);
288  else
289  hh = new TH2D(h->GetName(), h->GetTitle(),
290  xAxis->GetNbins(), xAxis->GetXmin(), xAxis->GetXmax(),
291  3, 0, 3);
292  hh->GetYaxis()->SetBinLabel(1, "Primary");
293  hh->GetYaxis()->SetBinLabel(2, "Secondary");
294  hh->GetYaxis()->SetBinLabel(3, "All");
295  hh->GetXaxis()->SetTitle("#eta");
296  out->Add(hh);
297  }
298  for (Int_t i = 0; i < nPrim; i++) {
299  TH2* h = static_cast<TH2*>(out->At(i));
300  if (!h) continue;
301 
302  TH1* hp = static_cast<TH1*>(prim->At(i));
303  TH1* hs = static_cast<TH1*>(sec->At(i));
304  TH1* ha = static_cast<TH1*>(all->At(i));
305  TH1* hh[] = { hp, hs, ha };
306  for (Int_t j = 0; j < 3; j++) {
307  TH1* ph = hh[j];
308  if (!ph) continue;
309 
310  for (Int_t k = 1; k <= ph->GetNbinsX(); k++) {
311  Double_t c = ph->GetBinContent(k);
312  Double_t e = ph->GetBinError(k);
313  h->SetBinContent(k, j+1, c);
314  h->SetBinError(k, j+1, e);
315  }
316  }
317  }
318  TList* l = GetOutputList(dir);
319  if (!l) return 0;
320 
321  out->SetName("comparative");
322  out->SetOwner();
323  l->Add(out);
324  return out;
325 }
326 //
327 // EOF
328 //
EResidualMethod fResidualMethod
double Double_t
Definition: External.C:58
Double_t BetaGamma() const
const char * title
Definition: MakeQAPdf.C:27
TH2 * GetBetaGammadEdx() const
virtual Bool_t PreEvent(const AliMCEvent &mcInput)
virtual void SetupForData(const TAxis &eAxis, const TAxis &cAxis, Double_t maxDE=10, Int_t nDEbins=300, Bool_t useIncrBin=true)
char Char_t
Definition: External.C:18
virtual void Scale(TH1 *dist) const
virtual void Scale(TH1 *dist) const
static const Spec & GetdEAxis()
TCanvas * c
Definition: TestFitELoss.C:172
AliFMDFloatMap & GetPrimaries()
AliStack * stack
TH2 * GetBetaGammaEta() const
virtual Bool_t Event(const AliESDEvent &esdInput, const AliMCEvent &mcInput, Double_t cent=-1)
TH2 * Make(const char *name, const char *title, const TAxis &eAxis, Double_t deMax=12, Int_t nDeBins=300, Bool_t incr=true)
virtual void CreateOutputObjects(TList *dir)
AliFMDFloatMap & GetEta()
int Int_t
Definition: External.C:63
float Float_t
Definition: External.C:68
Definition: External.C:228
virtual void FillMC(UShort_t flag, Double_t eta, Double_t mult)
void FillBinArray(TArrayD &a) const
TClonesArray * GetHits() const
Class to encode the energy loss and path length of a particle into track reference bits...
#define DGUARD(L, N, F,...)
Definition: External.C:220
Bool_t Calculate(const AliESDFMD &esd, const AliMCEvent &event, const TVector3 &ip, Double_t cent)
TArrayD MakeIncreasingAxis(Int_t nBins, Double_t low, Double_t high) const
AliFMDFloatMap & GetSecondaries()
unsigned short UShort_t
Definition: External.C:28
RingHistos * GetRingHistos(UShort_t d, Char_t r) const
virtual TObjArray * FitSlices(TList *dir, const char *name, 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, Bool_t scaleToPeak=true, TObjArray *best=0) const
bool Bool_t
Definition: External.C:53
TH2 * GetDEdxEta() const
void CreateOutputObjects(TList *list)
AliFMDEnergyFitter::RingHistos * CreateRingHistos(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
Definition: External.C:196
TList * GetOutputList(const TList *d) const
virtual void CreateOutputObjects(TList *dir)
AliFMDFloatMap & GetAll()
TDirectoryFile * dir