]>
Commit | Line | Data |
---|---|---|
bf000c32 | 1 | //____________________________________________________________________ |
8f6ee336 | 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> | |
8f6ee336 | 16 | #include <AliFMDDigit.h> |
15b17c89 | 17 | #include <AliFMDRecPoint.h> |
8f6ee336 | 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> | |
2b893216 | 30 | #include <TF1.h> |
8f6ee336 | 31 | |
bf000c32 | 32 | //____________________________________________________________________ |
9b48326f | 33 | /** @class DrawHitsRecs |
34 | @brief Draw hit energy loss versus rec point mult | |
35 | @code | |
36 | Root> .L Compile.C | |
37 | Root> Compile("DrawHitsRecs.C") | |
38 | Root> DrawHitsRecs c | |
39 | Root> c.Run(); | |
40 | @endcode | |
41 | @ingroup FMD_script | |
42 | */ | |
15b17c89 | 43 | class DrawHitsRecs : public AliFMDInput |
8f6ee336 | 44 | { |
45 | private: | |
15b17c89 | 46 | TH2D* fHitEvsAdc; |
47 | TH2D* fHitEvsRecM; // Histogram | |
48 | TH2D* fHitEvsRecE; // Histogram | |
49 | TH1D* fDiffE; // Histogram | |
50 | TH2D* fHitsVsRecM; // Histogram | |
f95a63c4 | 51 | TH2D* fDiffM; // Histogram |
2b893216 | 52 | TH1* fHitEloss; |
53 | TH1* fRecEloss; | |
15b17c89 | 54 | AliFMDEdepMap fMap; |
8f6ee336 | 55 | AliFMDFloatMap fEta; |
56 | AliFMDFloatMap fPhi; | |
57 | AliFMDFloatMap fMult; | |
58 | Bool_t fPrimary; | |
59 | public: | |
bf000c32 | 60 | //__________________________________________________________________ |
8f6ee336 | 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; | |
66 | AddLoad(kRecPoints); | |
8ec606c2 | 67 | AddLoad(kHits); |
15b17c89 | 68 | AddLoad(kDigits); |
8f6ee336 | 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; | |
15b17c89 | 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"); | |
f48d9b11 | 79 | fHitEvsAdc->Sumw2(); |
15b17c89 | 80 | |
81 | fHitEvsRecM = new TH2D("hitEvsRecM", "#Delta E_{sim} vs. M_{rec}", | |
8f6ee336 | 82 | eloss.fN-1, eloss.fArray, mults.fN-1, mults.fArray); |
15b17c89 | 83 | fHitEvsRecM->SetXTitle("#Delta E_{sim} [MeV]"); |
84 | fHitEvsRecM->SetYTitle("M_{rec}"); | |
f48d9b11 | 85 | fHitEvsRecM->Sumw2(); |
15b17c89 | 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]"); | |
f48d9b11 | 91 | fHitEvsRecE->Sumw2(); |
8f6ee336 | 92 | |
15b17c89 | 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}}"); | |
f48d9b11 | 98 | fDiffE->Sumw2(); |
15b17c89 | 99 | |
69893a66 | 100 | Double_t omin = mmin; // -.5; |
101 | Double_t omax = mmax; // 7.5; | |
102 | Int_t o = m; // 8; | |
15b17c89 | 103 | fHitsVsRecM = new TH2D("hitsVsStrip", "# of Hits vs. M_{rec}", |
8f6ee336 | 104 | o, omin, omax, m, mmin, mmax); |
15b17c89 | 105 | fHitsVsRecM->SetXTitle("# of Hits"); |
106 | fHitsVsRecM->SetYTitle("M_{rec}"); | |
f48d9b11 | 107 | fHitsVsRecM->Sumw2(); |
f95a63c4 | 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|"); | |
f48d9b11 | 114 | fDiffM->Sumw2(); |
f95a63c4 | 115 | // fDiffM->SetYTitle("Detector"); |
2b893216 | 116 | |
69893a66 | 117 | fHitEloss = new TH1D("hitEloss","#frac{#Delta E_{sim}}{#Delta x} (MeV/cm)", |
2b893216 | 118 | 100, 0, 10); |
119 | fHitEloss->SetFillColor(2); | |
120 | fHitEloss->SetFillStyle(3001); | |
f48d9b11 | 121 | fHitEloss->Sumw2(); |
122 | ||
69893a66 | 123 | fRecEloss = new TH1D("recEloss","#frac{#Delta E_{rec}}{#Delta x} (MeV/cm)", |
2b893216 | 124 | 100, 0, 10); |
125 | fRecEloss->SetFillColor(4); | |
126 | fRecEloss->SetFillStyle(3001); | |
f48d9b11 | 127 | fRecEloss->Sumw2(); |
8f6ee336 | 128 | } |
bf000c32 | 129 | //__________________________________________________________________ |
9b48326f | 130 | /** Begining of event |
131 | @param ev Event number | |
132 | @return @c false on error */ | |
8f6ee336 | 133 | Bool_t Begin(Int_t ev) |
134 | { | |
135 | fMap.Reset(); | |
15b17c89 | 136 | return AliFMDInput::Begin(ev); |
8f6ee336 | 137 | } |
bf000c32 | 138 | //__________________________________________________________________ |
8f6ee336 | 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 | ||
69893a66 | 157 | if (hit->Edep()/hit->Length() > 0.1) |
158 | fHitEloss->Fill(hit->Edep() / hit->Length()); | |
8f6ee336 | 159 | fMap(det, rng, sec, str).fEdep += hit->Edep(); |
160 | fMap(det, rng, sec, str).fN++; | |
161 | return kTRUE; | |
162 | } | |
15b17c89 | 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 | ||
bf000c32 | 183 | //__________________________________________________________________ |
184 | Bool_t ProcessRecPoint(AliFMDRecPoint* single) | |
8f6ee336 | 185 | { |
15b17c89 | 186 | if (!single) return kTRUE; |
bf000c32 | 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)); | |
15b17c89 | 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); | |
8f6ee336 | 201 | } |
15b17c89 | 202 | if (nhit > 0) fHitsVsRecM->Fill(nhit, single->Particles()); |
f95a63c4 | 203 | fDiffM->Fill(nhit - single->Particles(), TMath::Abs(single->Eta())); |
2b893216 | 204 | if (single->Edep()/.03 > 0.1) fRecEloss->Fill(single->Edep() / 0.0300); |
8f6ee336 | 205 | return kTRUE; |
206 | } | |
bf000c32 | 207 | //__________________________________________________________________ |
8f6ee336 | 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); | |
f95a63c4 | 216 | TCanvas* c = 0; |
8f6ee336 | 217 | |
f95a63c4 | 218 | c = new TCanvas("c0", fHitEvsAdc->GetTitle()); |
15b17c89 | 219 | fHitEvsAdc->SetStats(kFALSE); |
220 | fHitEvsAdc->Draw("COLZ"); | |
221 | ||
f95a63c4 | 222 | c = new TCanvas("c1", fHitEvsRecM->GetTitle()); |
15b17c89 | 223 | fHitEvsRecM->SetStats(kFALSE); |
224 | fHitEvsRecM->Draw("COLZ"); | |
225 | ||
f95a63c4 | 226 | c = new TCanvas("c2", fHitEvsRecE->GetTitle()); |
15b17c89 | 227 | fHitEvsRecE->SetStats(kFALSE); |
228 | fHitEvsRecE->Draw("COLZ"); | |
229 | ||
f95a63c4 | 230 | c = new TCanvas("c3", fDiffE->GetTitle()); |
231 | c->SetLogz(); | |
15b17c89 | 232 | fDiffE->Draw(); |
8f6ee336 | 233 | |
f95a63c4 | 234 | c = new TCanvas("c4", fHitsVsRecM->GetTitle()); |
235 | c->SetLogz(); | |
15b17c89 | 236 | fHitsVsRecM->SetStats(kFALSE); |
237 | fHitsVsRecM->Draw("COLZ"); | |
8f6ee336 | 238 | |
f95a63c4 | 239 | c = new TCanvas("c5", fDiffM->GetTitle()); |
240 | fDiffM->SetFillColor(2); | |
241 | fDiffM->SetFillStyle(3001); | |
242 | c->SetLogz(); | |
243 | fDiffM->Draw("colz"); | |
244 | ||
2b893216 | 245 | c = new TCanvas("c6", "Hit Eloss, Reco Eloss"); |
f48d9b11 | 246 | fRecEloss->Scale(1./fRecEloss->GetEntries()); |
247 | fRecEloss->Draw("hist e"); | |
2b893216 | 248 | fRecEloss->Fit("landau", "", "SAME", 2, 4); |
249 | TF1* recResp = new TF1(*fRecEloss->GetFunction("landau")); | |
f48d9b11 | 250 | fHitEloss->Scale(1./fHitEloss->GetEntries()); |
251 | fHitEloss->Draw("same hist e"); | |
2b893216 | 252 | fHitEloss->Fit("landau", "", "SAME", 2, 10); |
253 | TF1* hitResp = new TF1(*fHitEloss->GetFunction("landau")); | |
f48d9b11 | 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; | |
2b893216 | 265 | c->SetLogy(); |
f95a63c4 | 266 | |
8f6ee336 | 267 | return kTRUE; |
268 | } | |
269 | ||
270 | ClassDef(DrawHitsRecs,0); | |
271 | }; | |
272 | ||
273 | //____________________________________________________________________ | |
274 | // | |
275 | // EOF | |
276 | // |