AliRoot Core  3dc7879 (3dc7879)
DrawHitsRecs.C
Go to the documentation of this file.
1 //____________________________________________________________________
2 //
3 // $Id$
4 //
5 // Script that contains a class to draw eloss from hits, versus ADC
6 // counts from digits, using the AliFMDInputHits class in the util library.
7 //
8 // It draws the energy loss versus the p/(mq^2). It can be overlayed
9 // with the Bethe-Bloc curve to show how the simulation behaves
10 // relative to the expected.
11 //
12 // Use the script `Compile.C' to compile this class using ACLic.
13 //
14 #include <TH2D.h>
15 #include <AliFMDHit.h>
16 #include <AliFMDDigit.h>
17 #include <AliFMDRecPoint.h>
18 #include <AliFMDInput.h>
19 #include <AliFMDEdepMap.h>
20 #include <AliFMDFloatMap.h>
21 #include <iostream>
22 #include <TStyle.h>
23 #include <TArrayF.h>
24 #include <TCanvas.h>
25 #include <TParticle.h>
26 #include <TClonesArray.h>
27 #include <TTree.h>
28 #include <AliStack.h>
29 #include <AliLog.h>
30 #include <TF1.h>
31 
32 //____________________________________________________________________
43 class DrawHitsRecs : public AliFMDInput
44 {
45 private:
46  TH2D* fHitEvsAdc;
47  TH2D* fHitEvsRecM; // Histogram
48  TH2D* fHitEvsRecE; // Histogram
49  TH1D* fDiffE; // Histogram
50  TH2D* fHitsVsRecM; // Histogram
51  TH2D* fDiffM; // Histogram
52  TH1* fHitEloss;
53  TH1* fRecEloss;
58  Bool_t fPrimary;
59 public:
60  //__________________________________________________________________
61  DrawHitsRecs(Bool_t primary=kFALSE,
62  Int_t n=900, Double_t emin=1e-3, Double_t emax=10,
63  Int_t m=21, Double_t mmin=-0.5, Double_t mmax=20.5)
64  {
65  fPrimary = primary;
67  AddLoad(kHits);
69  if (fPrimary) AddLoad(kKinematics);
70  TArrayF eloss(MakeLogScale(n, emin, emax));
71  TArrayF mults(m+1);
72  mults[0] = mmin;
73  for (Int_t i = 1; i < m+1; i++) mults[i] = mults[i-1] + (mmax-mmin)/m;
74 
75  fHitEvsAdc = new TH2D("hitEvsAdc", "#Delta E_{sim} vs. ADC",
76  n, emin, emax, 1025, -.5, 1024.5);
77  fHitEvsAdc->SetXTitle("#Delta E_{sim} [MeV]");
78  fHitEvsAdc->SetYTitle("ADC");
79  fHitEvsAdc->Sumw2();
80 
81  fHitEvsRecM = new TH2D("hitEvsRecM", "#Delta E_{sim} vs. M_{rec}",
82  eloss.fN-1, eloss.fArray, mults.fN-1, mults.fArray);
83  fHitEvsRecM->SetXTitle("#Delta E_{sim} [MeV]");
84  fHitEvsRecM->SetYTitle("M_{rec}");
85  fHitEvsRecM->Sumw2();
86 
87  fHitEvsRecE = new TH2D("hitEvsRecE", "#Delta E_{sim} vs. #Delta E_{rec}",
88  n, emin, emax, n, emin, emax);
89  fHitEvsRecE->SetXTitle("#Delta E_{sim} [MeV]");
90  fHitEvsRecE->SetYTitle("#Delta E_{rec} [MeV]");
91  fHitEvsRecE->Sumw2();
92 
93 
94  fDiffE = new TH1D("diffE",
95  "#frac{#Delta E_{sim}-#Delta E_{rec}}{#Delta E_{sim}}",
96  1100, -1, 1.1);
97  fDiffE->SetXTitle("#frac{#Delta E_{sim}-#Delta E_{rec}}{#Delta E_{sim}}");
98  fDiffE->Sumw2();
99 
100  Double_t omin = mmin; // -.5;
101  Double_t omax = mmax; // 7.5;
102  Int_t o = m; // 8;
103  fHitsVsRecM = new TH2D("hitsVsStrip", "# of Hits vs. M_{rec}",
104  o, omin, omax, m, mmin, mmax);
105  fHitsVsRecM->SetXTitle("# of Hits");
106  fHitsVsRecM->SetYTitle("M_{rec}");
107  fHitsVsRecM->Sumw2();
108 
109  fDiffM = new TH2D("diffM", "M_{sim} - M_{rec}",
110  41, -20.5, 20.5, 70, 1.5, 5);
111  // 36, -TMath::Pi(),TMath::Pi());
112  fDiffM->SetXTitle("M_{sim} - M_{rec}");
113  fDiffM->SetYTitle("|#eta|");
114  fDiffM->Sumw2();
115  // fDiffM->SetYTitle("Detector");
116 
117  fHitEloss = new TH1D("hitEloss","#frac{#Delta E_{sim}}{#Delta x} (MeV/cm)",
118  100, 0, 10);
119  fHitEloss->SetFillColor(2);
120  fHitEloss->SetFillStyle(3001);
121  fHitEloss->Sumw2();
122 
123  fRecEloss = new TH1D("recEloss","#frac{#Delta E_{rec}}{#Delta x} (MeV/cm)",
124  100, 0, 10);
125  fRecEloss->SetFillColor(4);
126  fRecEloss->SetFillStyle(3001);
127  fRecEloss->Sumw2();
128  }
129  //__________________________________________________________________
133  Bool_t Begin(Int_t ev)
134  {
135  fMap.Reset();
136  return AliFMDInput::Begin(ev);
137  }
138  //__________________________________________________________________
139  Bool_t ProcessHit(AliFMDHit* hit, TParticle*)
140  {
141  // Cache the energy loss
142  if (!hit) return kFALSE;
143  UShort_t det = hit->Detector();
144  Char_t rng = hit->Ring();
145  UShort_t sec = hit->Sector();
146  UShort_t str = hit->Strip();
147  if (str > 511) {
148  AliWarning(Form("Bad strip number %d in hit", str));
149  return kTRUE;
150  }
151  if (fPrimary) {
152  TParticle* kine = fStack->Particle(hit->Track());
153  if (!kine) return kTRUE;
154  if (!kine->IsPrimary()) return kTRUE;
155  }
156 
157  if (hit->Edep()/hit->Length() > 0.1)
158  fHitEloss->Fill(hit->Edep() / hit->Length());
159  fMap(det, rng, sec, str).fEdep += hit->Edep();
160  fMap(det, rng, sec, str).fN++;
161  return kTRUE;
162  }
163 
164  //__________________________________________________________________
165  Bool_t ProcessDigit(AliFMDDigit* digit)
166  {
167  if (!digit) return kTRUE;
168 
169  UShort_t det = digit->Detector();
170  Char_t rng = digit->Ring();
171  UShort_t sec = digit->Sector();
172  UShort_t str = digit->Strip();
173  if (str > 511) {
174  AliWarning(Form("Bad strip number %d in digit", str));
175  return kTRUE;
176  }
177  Double_t edep = fMap(det, rng, sec, str).fEdep;
178  if (edep > 0) fHitEvsAdc->Fill(edep, digit->Counts());
179 
180  return kTRUE;
181  }
182 
183  //__________________________________________________________________
185  {
186  if (!single) return kTRUE;
187  UShort_t det = single->Detector();
188  Char_t rng = single->Ring();
189  UShort_t sec = single->Sector();
190  UShort_t str = single->Strip();
191  if (str > 511) {
192  AliWarning(Form("Bad strip number %d in single", str));
193  return kTRUE;
194  }
195  Double_t edep = fMap(det, rng, sec, str).fEdep;
196  Int_t nhit = fMap(det, rng, sec, str).fN;
197  if (edep > 0) {
198  fHitEvsRecM->Fill(edep, single->Particles());
199  fHitEvsRecE->Fill(edep, single->Edep());
200  fDiffE->Fill((single->Edep() - edep) / edep);
201  }
202  if (nhit > 0) fHitsVsRecM->Fill(nhit, single->Particles());
203  fDiffM->Fill(nhit - single->Particles(), TMath::Abs(single->Eta()));
204  if (single->Edep()/.03 > 0.1) fRecEloss->Fill(single->Edep() / 0.0300);
205  return kTRUE;
206  }
207  //__________________________________________________________________
208  Bool_t Finish()
209  {
210  gStyle->SetPalette(1);
211  gStyle->SetOptTitle(0);
212  gStyle->SetCanvasColor(0);
213  gStyle->SetCanvasBorderSize(0);
214  gStyle->SetPadColor(0);
215  gStyle->SetPadBorderSize(0);
216  TCanvas* c = 0;
217 
218  c = new TCanvas("c0", fHitEvsAdc->GetTitle());
219  fHitEvsAdc->SetStats(kFALSE);
220  fHitEvsAdc->Draw("COLZ");
221 
222  c = new TCanvas("c1", fHitEvsRecM->GetTitle());
223  fHitEvsRecM->SetStats(kFALSE);
224  fHitEvsRecM->Draw("COLZ");
225 
226  c = new TCanvas("c2", fHitEvsRecE->GetTitle());
227  fHitEvsRecE->SetStats(kFALSE);
228  fHitEvsRecE->Draw("COLZ");
229 
230  c = new TCanvas("c3", fDiffE->GetTitle());
231  c->SetLogz();
232  fDiffE->Draw();
233 
234  c = new TCanvas("c4", fHitsVsRecM->GetTitle());
235  c->SetLogz();
236  fHitsVsRecM->SetStats(kFALSE);
237  fHitsVsRecM->Draw("COLZ");
238 
239  c = new TCanvas("c5", fDiffM->GetTitle());
240  fDiffM->SetFillColor(2);
241  fDiffM->SetFillStyle(3001);
242  c->SetLogz();
243  fDiffM->Draw("colz");
244 
245  c = new TCanvas("c6", "Hit Eloss, Reco Eloss");
246  fRecEloss->Scale(1./fRecEloss->GetEntries());
247  fRecEloss->Draw("hist e");
248  fRecEloss->Fit("landau", "", "SAME", 2, 4);
249  TF1* recResp = new TF1(*fRecEloss->GetFunction("landau"));
250  fHitEloss->Scale(1./fHitEloss->GetEntries());
251  fHitEloss->Draw("same hist e");
252  fHitEloss->Fit("landau", "", "SAME", 2, 10);
253  TF1* hitResp = new TF1(*fHitEloss->GetFunction("landau"));
254  std::cout << "From hits: MPV="
255  << hitResp->GetParameter(1) << "+/-"
256  << hitResp->GetParError(1) << " Width="
257  << hitResp->GetParameter(2) << "+/-"
258  << hitResp->GetParError(2) << "\nFrom recs: MPV="
259  << recResp->GetParameter(1) << "+/-"
260  << recResp->GetParError(1) << " Width="
261  << recResp->GetParameter(2) << "+/-"
262  << recResp->GetParError(2) << "\nRatio: MPV(hit/rec)="
263  << hitResp->GetParameter(1) / recResp->GetParameter(1)
264  << std::endl;
265  c->SetLogy();
266 
267  return kTRUE;
268  }
269 
270  ClassDef(DrawHitsRecs,0);
271 };
272 
273 //____________________________________________________________________
274 //
275 // EOF
276 //
class for digits
Definition: AliFMDDigit.h:28
Bool_t Finish()
Definition: DrawHitsRecs.C:208
Reconstructed FMD points. It contains the pseudo-inclusive multiplicity.
Map of Energy deposited, hit information per strip. Contains a pair of energy deposited fEdep and num...
Definition: AliFMDEdepMap.h:33
TStyle * gStyle
AliFMDFloatMap fMult
Definition: DrawHitsRecs.C:57
Hit in the FMD.
Digits for the FMD.
FMD utility classes for reading FMD data.
UShort_t Sector() const
Char_t Ring() const
Per strip map of energy deposited and number of hits.
UShort_t Sector() const
TH2D * fHitEvsRecE
Definition: DrawHitsRecs.C:48
AliFMDFloatMap fPhi
Definition: DrawHitsRecs.C:56
Bool_t Begin(Int_t ev)
Definition: DrawHitsRecs.C:133
virtual Bool_t Begin(Int_t event)
UShort_t Detector() const
TH1D * fDiffE
Definition: DrawHitsRecs.C:49
Pseudo reconstructed charged particle multiplicity.
TH2D * fDiffM
Definition: DrawHitsRecs.C:51
#define AliWarning(message)
Definition: AliLog.h:541
Int_t Track() const
Definition: AliHit.h:24
Bool_t ProcessDigit(AliFMDDigit *digit)
Definition: DrawHitsRecs.C:165
Float_t Edep() const
DrawHitsRecs(Bool_t primary=kFALSE, Int_t n=900, Double_t emin=1e-3, Double_t emax=10, Int_t m=21, Double_t mmin=-0.5, Double_t mmax=20.5)
Definition: DrawHitsRecs.C:61
UShort_t Strip() const
Float_t Edep() const
Definition: AliFMDHit.h:82
Char_t Ring() const
Definition: AliFMDHit.h:76
TH1 * fHitEloss
Definition: DrawHitsRecs.C:52
Float_t Eta() const
Bool_t fPrimary
Definition: DrawHitsRecs.C:58
static TArrayF MakeLogScale(Int_t n, Double_t min, Double_t max)
UShort_t Detector() const
Definition: AliFMDHit.h:74
UShort_t Strip() const
Draw hit energy loss versus rec point mult.
Definition: DrawHitsRecs.C:43
virtual void AddLoad(ETrees tree)
Definition: AliFMDInput.h:134
UShort_t Sector() const
Definition: AliFMDHit.h:78
TH2D * fHitsVsRecM
Definition: DrawHitsRecs.C:50
Bool_t ProcessHit(AliFMDHit *hit, TParticle *)
Definition: DrawHitsRecs.C:139
Char_t Ring() const
AliStack * fStack
Definition: AliFMDInput.h:416
TH2D * fHitEvsAdc
Definition: DrawHitsRecs.C:46
virtual void Reset()
Base class for reading in various FMD data. The class loops over all found events. For each event the specified data is read in. The class then loops over all elements of the read data, and process these with user defined code.
Definition: AliFMDInput.h:106
TParticle * Particle(Int_t id, Bool_t useInEmbedding=kFALSE)
Definition: AliStack.cxx:689
UShort_t Strip() const
Definition: AliFMDHit.h:80
Float_t Length() const
Definition: AliFMDHit.h:100
AliFMDEdepMap fMap
Definition: DrawHitsRecs.C:54
AliFMDhit is the hit class for the FMD. Hits are the information that comes from a Monte Carlo at eac...
Definition: AliFMDHit.h:30
AliFMDFloatMap fEta
Definition: DrawHitsRecs.C:55
UShort_t Detector() const
UShort_t Counts() const
Definition: AliFMDDigit.h:147
Bool_t ProcessRecPoint(AliFMDRecPoint *single)
Definition: DrawHitsRecs.C:184
TH2D * fHitEvsRecM
Definition: DrawHitsRecs.C:47
Float_t Particles() const
TH1 * fRecEloss
Definition: DrawHitsRecs.C:53