98e3a39afb19e8b831cb0d6f2d74abb5fdc193da
[u/mrichter/AliRoot.git] / FMD / scripts / DrawHits.C
1 //____________________________________________________________________
2 //
3 // $Id$
4 //
5 // Script that contains a class to draw hits, using the
6 // 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 <AliFMDInput.h>
17 #include <iostream>
18 #include <TStyle.h>
19 #include <TArrayF.h>
20 #include <TParticle.h>
21 #include <TCanvas.h>
22 #include <TGraphErrors.h>
23
24 /** @class DrawHits
25     @brief Draw hit energy loss
26     @code 
27     Root> .L Compile.C
28     Root> Compile("DrawHits.C")
29     Root> DrawHits c
30     Root> c.Run();
31     @endcode
32     @ingroup FMD_script
33  */
34 class DrawHits : public AliFMDInput
35 {
36 private:
37   TH2D* fElossVsPMQ; // Histogram 
38   Int_t fPdg;
39 public:
40   //__________________________________________________________________
41   TArrayF MakeLogScale(Int_t n, Double_t min, Double_t max) 
42   {
43     TArrayF bins(n+1);
44     bins[0]      = min;
45     if (n <= 20) {
46       for (Int_t i = 1; i < n+1; i++) bins[i] = bins[i-1] + (max-min)/n;
47       return bins;
48     }
49     Float_t dp   = n / TMath::Log10(max / min);
50     Float_t pmin = TMath::Log10(min);
51     for (Int_t i = 1; i < n+1; i++) {
52       Float_t p = pmin + i / dp;
53       bins[i]   = TMath::Power(10, p);
54     }
55     return bins;
56   }
57   //__________________________________________________________________
58   DrawHits(Int_t m=1000, Double_t emin=1, Double_t emax=1000, 
59            Int_t n=900, Double_t tmin=1e-2, Double_t tmax=1e3, 
60            Int_t pdg=211) 
61     : fPdg(pdg)
62   { 
63     AddLoad(kKinematics);
64     AddLoad(kHits);
65     TArrayF tkine(MakeLogScale(n, tmin, tmax));
66     TArrayF eloss(MakeLogScale(m, emin, emax));
67     TString name(Form("elossVsP%s", (fPdg == 0 ? "MQ" : "")));
68     TString title(Form("#Delta E/#Delta x vs. p%s", 
69                        (fPdg == 0 ? "/(mq^{2})" : "")));
70     fElossVsPMQ = new TH2D(name.Data(), title.Data(), 
71                            tkine.fN-1, tkine.fArray, 
72                            eloss.fN-1, eloss.fArray);
73     fElossVsPMQ->SetXTitle(Form("p%s", (fPdg == 0 ? "/(mq^{2})" : " [GeV]")));
74     fElossVsPMQ->SetYTitle("#Delta E/#Delta x [MeV/cm]");
75   }
76   //__________________________________________________________________
77   Bool_t ProcessHit(AliFMDHit* hit, TParticle* p) 
78   {
79     if (!hit) {
80       std::cout << "No hit" << std::endl;
81       return kFALSE;
82     }
83
84     if (!p) {
85       std::cout << "No track" << std::endl;
86       return kFALSE;
87     }
88     if (!p->IsPrimary()) return kTRUE;
89     if (hit->IsStop()) return kTRUE;
90     Float_t x = hit->P();
91     Float_t y = hit->Edep()/hit->Length();
92     if (fPdg != 0) {
93       if (TMath::Abs(hit->Pdg()) != fPdg) return kTRUE;
94     }
95     else {
96       Float_t q = hit->Q() / 3.;
97       if (hit->M() == 0 || q == 0) return kTRUE;
98       x /= hit->M();
99       y /= q * q;
100     }
101     fElossVsPMQ->Fill(x, y);
102     return kTRUE;
103   }
104   //__________________________________________________________________
105   Bool_t Finish()
106   {
107     TCanvas* c = new TCanvas("c", "C");
108     c->SetLogy();
109     c->SetLogx();
110     gStyle->SetPalette(1);
111     gStyle->SetOptTitle(0);
112     gStyle->SetCanvasColor(0);
113     gStyle->SetCanvasBorderSize(0);
114     gStyle->SetPadColor(0);
115     gStyle->SetPadBorderSize(0);
116     fElossVsPMQ->SetStats(kFALSE);
117     fElossVsPMQ->Draw("COLZ box");
118     c->Modified();
119     c->Update();
120     c->cd();
121     return kTRUE;
122   }
123   void SuperImposeBetheBloc() 
124   {
125     // This is for pi+, made with MakeXsection.C and DrawXsection.C
126     TGraphErrors *gre = new TGraphErrors(91);
127     gre->SetName("BetheBlocPi");
128     gre->SetTitle("BetheBlocPi");
129     gre->SetFillColor(6);
130     gre->SetFillStyle(3001);
131     gre->SetLineWidth(2);
132     gre->SetPoint(0,7.16486e-05,1218.84);
133     gre->SetPointError(0,0,609.418);
134     gre->SetPoint(1,9.25378e-05,1221.38);
135     gre->SetPointError(1,0,610.689);
136     gre->SetPoint(2,0.000119517,1180.12);
137     gre->SetPointError(2,0,590.058);
138     gre->SetPoint(3,0.000154362,1100.31);
139     gre->SetPointError(3,0,550.156);
140     gre->SetPoint(4,0.000199367,996.621);
141     gre->SetPointError(4,0,498.31);
142     gre->SetPoint(5,0.000257492,886.005);
143     gre->SetPointError(5,0,443.003);
144     gre->SetPoint(6,0.000332563,780.483);
145     gre->SetPointError(6,0,390.241);
146     gre->SetPoint(7,0.000429522,684.927);
147     gre->SetPointError(7,0,342.463);
148     gre->SetPoint(8,0.000554749,599.407);
149     gre->SetPointError(8,0,299.703);
150     gre->SetPoint(9,0.000716486,522.375);
151     gre->SetPointError(9,0,261.187);
152     gre->SetPoint(10,0.000925378,452.497);
153     gre->SetPointError(10,0,226.249);
154     gre->SetPoint(11,0.00119517,389.101);
155     gre->SetPointError(11,0,194.551);
156     gre->SetPoint(12,0.00154362,331.974);
157     gre->SetPointError(12,0,165.987);
158     gre->SetPoint(13,0.00199367,280.969);
159     gre->SetPointError(13,0,140.485);
160     gre->SetPoint(14,0.00257492,235.689);
161     gre->SetPointError(14,0,117.844);
162     gre->SetPoint(15,0.00332564,196.156);
163     gre->SetPointError(15,0,98.078);
164     gre->SetPoint(16,0.00429522,162.402);
165     gre->SetPointError(16,0,81.2012);
166     gre->SetPoint(17,0.00554749,133.87);
167     gre->SetPointError(17,0,66.9351);
168     gre->SetPoint(18,0.00716486,109.959);
169     gre->SetPointError(18,0,54.9797);
170     gre->SetPoint(19,0.00925378,90.2035);
171     gre->SetPointError(19,0,45.1017);
172     gre->SetPoint(20,0.0119517,74.1317);
173     gre->SetPointError(20,0,37.0658);
174     gre->SetPoint(21,0.0154362,60.8988);
175     gre->SetPointError(21,0,30.4494);
176     gre->SetPoint(22,0.0199367,49.9915);
177     gre->SetPointError(22,0,24.9957);
178     gre->SetPoint(23,0.0257492,40.9812);
179     gre->SetPointError(23,0,20.4906);
180     gre->SetPoint(24,0.0332564,33.5739);
181     gre->SetPointError(24,0,16.7869);
182     gre->SetPoint(25,0.0429522,27.5127);
183     gre->SetPointError(25,0,13.7563);
184     gre->SetPoint(26,0.0554749,22.5744);
185     gre->SetPointError(26,0,11.2872);
186     gre->SetPoint(27,0.0716486,18.5674);
187     gre->SetPointError(27,0,9.28372);
188     gre->SetPoint(28,0.0925378,15.3292);
189     gre->SetPointError(28,0,7.66462);
190     gre->SetPoint(29,0.119517,12.7231);
191     gre->SetPointError(29,0,6.36156);
192     gre->SetPoint(30,0.154362,10.6352);
193     gre->SetPointError(30,0,5.31759);
194     gre->SetPoint(31,0.199367,8.97115);
195     gre->SetPointError(31,0,4.48558);
196     gre->SetPoint(32,0.257492,7.65358);
197     gre->SetPointError(32,0,3.82679);
198     gre->SetPoint(33,0.332564,6.61909);
199     gre->SetPointError(33,0,3.30955);
200     gre->SetPoint(34,0.429522,5.81614);
201     gre->SetPointError(34,0,2.90807);
202     gre->SetPoint(35,0.554749,5.20286);
203     gre->SetPointError(35,0,2.60143);
204     gre->SetPoint(36,0.716486,4.74533);
205     gre->SetPointError(36,0,2.37267);
206     gre->SetPoint(37,0.925378,4.40987);
207     gre->SetPointError(37,0,2.20494);
208     gre->SetPoint(38,1.19517,4.17077);
209     gre->SetPointError(38,0,2.08538);
210     gre->SetPoint(39,1.54362,4.014);
211     gre->SetPointError(39,0,2.007);
212     gre->SetPoint(40,1.99367,3.92577);
213     gre->SetPointError(40,0,1.96288);
214     gre->SetPoint(41,2.57492,3.89199);
215     gre->SetPointError(41,0,1.946);
216     gre->SetPoint(42,3.32564,3.90063);
217     gre->SetPointError(42,0,1.95032);
218     gre->SetPoint(43,4.29522,3.94146);
219     gre->SetPointError(43,0,1.97073);
220     gre->SetPoint(44,5.54749,4.00597);
221     gre->SetPointError(44,0,2.00299);
222     gre->SetPoint(45,7.16486,4.08725);
223     gre->SetPointError(45,0,2.04362);
224     gre->SetPoint(46,9.25378,4.17985);
225     gre->SetPointError(46,0,2.08993);
226     gre->SetPoint(47,11.9517,4.27962);
227     gre->SetPointError(47,0,2.13981);
228     gre->SetPoint(48,15.4362,4.38347);
229     gre->SetPointError(48,0,2.19174);
230     gre->SetPoint(49,19.9367,4.48919);
231     gre->SetPointError(49,0,2.2446);
232     gre->SetPoint(50,25.7492,4.59523);
233     gre->SetPointError(50,0,2.29762);
234     gre->SetPoint(51,33.2564,4.70052);
235     gre->SetPointError(51,0,2.35026);
236     gre->SetPoint(52,42.9522,4.80435);
237     gre->SetPointError(52,0,2.40218);
238     gre->SetPoint(53,55.4749,4.90625);
239     gre->SetPointError(53,0,2.45312);
240     gre->SetPoint(54,71.6486,5.00589);
241     gre->SetPointError(54,0,2.50295);
242     gre->SetPoint(55,92.5378,5.10279);
243     gre->SetPointError(55,0,2.55139);
244     gre->SetPoint(56,119.517,5.19654);
245     gre->SetPointError(56,0,2.59827);
246     gre->SetPoint(57,154.362,5.28758);
247     gre->SetPointError(57,0,2.64379);
248     gre->SetPoint(58,199.367,5.37581);
249     gre->SetPointError(58,0,2.6879);
250     gre->SetPoint(59,257.492,5.46109);
251     gre->SetPointError(59,0,2.73055);
252     gre->SetPoint(60,332.564,5.54335);
253     gre->SetPointError(60,0,2.77167);
254     gre->SetPoint(61,429.522,5.62248);
255     gre->SetPointError(61,0,2.81124);
256     gre->SetPoint(62,554.749,5.69843);
257     gre->SetPointError(62,0,2.84922);
258     gre->SetPoint(63,716.486,5.77122);
259     gre->SetPointError(63,0,2.88561);
260     gre->SetPoint(64,925.378,5.84093);
261     gre->SetPointError(64,0,2.92046);
262     gre->SetPoint(65,1195.17,5.9077);
263     gre->SetPointError(65,0,2.95385);
264     gre->SetPoint(66,1543.62,5.97165);
265     gre->SetPointError(66,0,2.98582);
266     gre->SetPoint(67,1993.67,6.03292);
267     gre->SetPointError(67,0,3.01646);
268     gre->SetPoint(68,2574.92,6.09171);
269     gre->SetPointError(68,0,3.04586);
270     gre->SetPoint(69,3325.64,6.14827);
271     gre->SetPointError(69,0,3.07413);
272     gre->SetPoint(70,4295.22,6.20286);
273     gre->SetPointError(70,0,3.10143);
274     gre->SetPoint(71,5547.49,6.25577);
275     gre->SetPointError(71,0,3.12788);
276     gre->SetPoint(72,7164.86,6.30725);
277     gre->SetPointError(72,0,3.15363);
278     gre->SetPoint(73,9253.78,6.35757);
279     gre->SetPointError(73,0,3.17878);
280     gre->SetPoint(74,11951.7,6.38446);
281     gre->SetPointError(74,0,3.19223);
282     gre->SetPoint(75,15436.2,6.38446);
283     gre->SetPointError(75,0,3.19223);
284     gre->SetPoint(76,19936.7,6.38446);
285     gre->SetPointError(76,0,3.19223);
286     gre->SetPoint(77,25749.2,6.38446);
287     gre->SetPointError(77,0,3.19223);
288     gre->SetPoint(78,33256.4,6.38446);
289     gre->SetPointError(78,0,3.19223);
290     gre->SetPoint(79,42952.2,6.38446);
291     gre->SetPointError(79,0,3.19223);
292     gre->SetPoint(80,55474.9,6.38446);
293     gre->SetPointError(80,0,3.19223);
294     gre->SetPoint(81,71648.6,6.38446);
295     gre->SetPointError(81,0,3.19223);
296     gre->SetPoint(82,92537.8,6.38446);
297     gre->SetPointError(82,0,3.19223);
298     gre->SetPoint(83,119517,6.38446);
299     gre->SetPointError(83,0,3.19223);
300     gre->SetPoint(84,154362,6.38446);
301     gre->SetPointError(84,0,3.19223);
302     gre->SetPoint(85,199367,6.38446);
303     gre->SetPointError(85,0,3.19223);
304     gre->SetPoint(86,257492,6.38446);
305     gre->SetPointError(86,0,3.19223);
306     gre->SetPoint(87,332563,6.38446);
307     gre->SetPointError(87,0,3.19223);
308     gre->SetPoint(88,429522,6.38446);
309     gre->SetPointError(88,0,3.19223);
310     gre->SetPoint(89,554749,6.38446);
311     gre->SetPointError(89,0,3.19223);
312     gre->SetPoint(90,716486,6.38446);
313     gre->SetPointError(90,0,3.19223);
314     gre->Draw("l same");
315     gre->DrawClone("l3 same");
316     gPad->Modified();
317     gPad->Update();
318     gPad->cd();
319   }
320   
321   ClassDef(DrawHits,0);
322 };
323
324 //____________________________________________________________________
325 //
326 // EOF
327 //