source: trunk/WuerzburgSoft/Thomas/mphys/analp.C@ 1428

Last change on this file since 1428 was 1428, checked in by tbretz, 22 years ago
*** empty log message ***
File size: 10.3 KB
Line 
1//Double_t kRad2Deg = 180./TMath::Pi();
2
3#include <float.h>
4
5Double_t GetMin(TH1 *h)
6{
7 Double_t min = FLT_MAX;
8 for (int i=1; i<=h->GetXaxis()->GetNbins(); i++)
9 {
10 if (h->GetBinContent(i)<min && h->GetBinContent(i)!=0)
11 min = h->GetBinContent(i);
12 }
13 return min;
14}
15
16void GetRange(TChain *chain, const char *name, Int_t &nbins, Double_t &min, Double_t &max, Double_t conv=1)
17{
18 TString str("MPhoton.MParticle.");
19
20 MPhoton *p = new MPhoton;
21
22 chain->SetBranchAddress("MPhoton.", &p);
23
24 min = FLT_MAX;
25 max = -FLT_MAX;
26
27 Int_t n = chain->GetEntries();
28 for (int i=0; i<n; i++)
29 {
30 chain->GetEntry(i);
31
32 TLeaf *leaf = chain->GetLeaf(str+name);
33 if (!leaf)
34 return 0;
35
36 Double_t val = leaf->GetValue();
37
38 if (p->IsPrimary())
39 continue;
40
41 if (val < min)
42 min = val;
43 if (val > max)
44 max = val;
45 }
46
47 min *= conv;
48 max *= conv;
49
50 chain.GetPlayer();
51 TTreePlayer &play = *chain.GetPlayer();
52
53 Int_t n;
54 play.FindGoodLimits(nbins, n, min, max, kFALSE);
55 nbins = n;
56}
57
58void draw1h1426(Double_t E, Double_t plot=1)
59{
60 plot=2;
61 Double_t level = log10(E);
62 level -= log10(5.73e-0*pow( 0.39e3, plot));
63
64 cout << "Level: " << level << endl;
65
66 TPolyMarker *m = new TPolyMarker(7);
67 /*
68 m->SetPoint(0, 0.28e3, 1.86e-0*pow( 0.28e3, plot)*level);
69 m->SetPoint(1, 0.39e3, 5.73e-0*pow( 0.39e3, plot)*level);
70 m->SetPoint(2, 0.80e3, 1.22e-0*pow( 0.80e3, plot)*level);
71 m->SetPoint(3, 1.55e3, 5.10e-2*pow( 1.55e3, plot)*level);
72 m->SetPoint(4, 2.82e3, 1.50e-2*pow( 2.82e3, plot)*level);
73 m->SetPoint(5, 5.33e3, 7.70e-3*pow( 5.33e3, plot)*level);
74 m->SetPoint(6, 10.00e3, 1.20e-4*pow(10.00e3, plot)*level);
75 */
76 m->SetPoint(0, log10( 0.28e3), log10(1.86e-0*pow( 0.28e3, plot))+level);
77 m->SetPoint(1, log10( 0.39e3), log10(5.73e-0*pow( 0.39e3, plot))+level);
78 m->SetPoint(2, log10( 0.80e3), log10(1.22e-0*pow( 0.80e3, plot))+level);
79 m->SetPoint(3, log10( 1.55e3), log10(5.10e-2*pow( 1.55e3, plot))+level);
80 m->SetPoint(4, log10( 2.82e3), log10(1.50e-2*pow( 2.82e3, plot))+level);
81 m->SetPoint(5, log10( 5.33e3), log10(7.70e-3*pow( 5.33e3, plot))+level);
82 m->SetPoint(6, log10(10.00e3), log10(1.20e-4*pow(10.00e3, plot))+level);
83 m->SetMarkerStyle(kOpenStar);
84// m->SetMarkerSize(1);
85 m->Draw();
86}
87
88void analp()
89{
90 TChain chain("Photons");
91 /*
92 chain.Add("delme3.root");
93 chain.Add("delme3B.root");
94 chain.Add("delme3C.root");
95 chain.Add("delme3D.root");
96 chain.Add("delme3E.root");
97 */
98 chain.Add("delme3H.root");
99
100 Int_t nbins = 24;
101 Double_t lo = 1e2;
102 Double_t hi = 1e6;
103
104 MBinning binsx;
105 binsx.SetEdgesLog(nbins, lo, hi);
106
107 // ------------------------------
108
109 Double_t cutoff = 1e12;
110
111 Double_t alpha = -1;
112 Double_t plot = 1;
113
114 Int_t nbins = 16;
115
116 // ------------------------------
117
118 Int_t n = chain.GetEntries();
119
120 cout << endl << "Found " << n << " entries." << endl;
121
122 if (n==0)
123 return;
124
125 MBinning binspolx;
126 MBinning binspola;
127 MBinning binspolr;
128 binspolx.SetEdges(16, -180, 180);
129
130 Double_t min, max;
131 GetRange(&chain, "fTheta", nbins, min, max, kRad2Deg*60);
132 binspola.SetEdges(nbins, min, max*1.1);
133 cout << "fTheta, " << nbins << ": " << min << " ... " << max << " [']" << endl;
134
135 GetRange(&chain, "fR", nbins, min, max);
136 binspolr.SetEdges(nbins, min, max*1.1);
137 cout << "fR, " << nbins << ": " << min << " ... " << max << " [kpc]" << endl;
138
139 TH1D h, h2, h3;
140 TH1D prim;
141 MH::SetBinning(&h, &binsx);
142 MH::SetBinning(&h2, &binsx);
143 MH::SetBinning(&h3, &binsx);
144 MH::SetBinning(&prim, &binsx);
145
146 TH2D r, r2, r3;
147 TH2D a, a2, a3;
148 MH::SetBinning(&a, &binspolx, &binspola);
149 MH::SetBinning(&r, &binspolx, &binspolr);
150 MH::SetBinning(&a2, &binspolx, &binspola);
151 MH::SetBinning(&r2, &binspolx, &binspolr);
152 MH::SetBinning(&a3, &binspolx, &binspola);
153 MH::SetBinning(&r3, &binspolx, &binspolr);
154
155 MPhoton *p = new MPhoton;
156 chain.SetBranchAddress("MPhoton.", &p);
157
158 chain.GetEntry(0);
159 if (!p->IsPrimary())
160 return;
161
162 Double_t z = p->GetZ();
163 Double_t R = MParticle::RofZ(&z);
164
165 cout << "Z = " << z << endl;
166 cout << "R = " << R << "kpc" << endl;
167
168 Double_t weight = 0;
169
170 Bool_t isprim = kFALSE;
171 Bool_t ignore = kFALSE;
172 for (int i=0; i<n; i++)
173 {
174 chain.GetEntry(i);
175
176 Double_t Ep = p->GetEnergy();
177
178 if (p->IsPrimary())
179 {
180 if (Ep>cutoff)
181 {
182 ignore = kTRUE;
183 continue;
184 }
185 ignore = kFALSE;
186 weight = pow(Ep, alpha);
187 prim.Fill(Ep, pow(Ep, plot+alpha));
188 isprim = kTRUE;
189 continue;
190 }
191
192 if (ignore)
193 continue;
194
195 h.Fill(Ep, pow(Ep,plot) * weight);
196
197 Double_t phi = fmod(p->GetPhi()*kRad2Deg+180, 360)-180;
198 Double_t psi = fmod(p->GetPsi()*kRad2Deg+180, 360)-180;
199
200 r.Fill(phi, isprim?0:p->GetR(), weight);
201 a.Fill(psi, isprim?0:p->GetTheta()*kRad2Deg*60, weight);
202
203 if (isprim)
204 {
205 h3.Fill(Ep, pow(Ep,plot) * weight);
206 r3.Fill(phi, 0, weight);
207 a3.Fill(psi, 0, weight);
208 isprim = kFALSE;
209 continue;
210 }
211
212 h2.Fill(Ep, pow(Ep,plot) * weight);
213 r2.Fill(phi, p->GetR(), weight);
214 a2.Fill(psi, p->GetTheta()*kRad2Deg*60, weight);
215 }
216
217 delete p;
218
219 TH1 *g=r.ProjectionX();
220 cout << "Mean fPhi: " << g->GetMean() << " +- " << g->GetRMS() << endl;
221 delete g;
222 g=a.ProjectionX();
223 cout << "Mean fPsi: " << g->GetMean() << " +- " << g->GetRMS() << endl;
224 delete g;
225
226 gStyle->SetOptStat(10);
227
228 TLine line;
229
230// delete gROOT->FindObject("Analysis Arrival");
231 c = MH::MakeDefCanvas("Analysis Arrival", "");
232 c->Divide(2,2);
233
234 gStyle->SetPalette(1, 0);
235
236 c->cd(1);
237 r.SetTitle(" Distance from Observer ");
238 r.GetXaxis()->SetLabelOffset(-0.015);
239 r.GetXaxis()->SetTitleOffset(1.1);
240 r.GetXaxis()->SetRangeUser(1e4, 1e9);
241 r.SetXTitle("\\Phi [\\circ]");
242 r.SetYTitle("R [kpc]");
243 TPad &pad = *gPad;
244 pad.Divide(2,2);
245 pad.cd(1);
246 gPad->SetLogy();
247 gPad->SetLogz();
248 gPad->SetTheta(0);
249 gPad->SetPhi(90);
250 g = r.DrawCopy("surf1pol");
251 pad.cd(2);
252 gPad->SetLogy();
253 gPad->SetLogz();
254 gPad->SetTheta(70);
255 gPad->SetPhi(90);
256 g->Draw("surf1pol");
257 pad.cd(3);
258 gPad->SetLogy();
259 gPad->SetLogz();
260 gPad->SetTheta(90);
261 gPad->SetPhi(90);
262 g->Draw("surf1pol");
263 pad.cd(4);
264 gPad->SetLogy();
265 gPad->SetLogz();
266 gPad->SetTheta(20);
267 gPad->SetPhi(90);
268 g->Draw("surf1pol");
269
270 c->cd(2);
271 a.SetTitle(" Hit Angle for Observer ");
272 a.GetXaxis()->SetLabelOffset(-0.015);
273 a.GetXaxis()->SetTitleOffset(1.1);
274 a.GetXaxis()->SetRangeUser(1e4, 1e9);
275 a.SetXTitle("\\Psi [\\circ]");
276 a.SetYTitle("\\Theta [']");
277 TPad &pad = *gPad;
278 pad.Divide(2,2);
279 pad.cd(1);
280 gPad->SetLogy();
281 gPad->SetLogz();
282 gPad->SetTheta(0);
283 gPad->SetPhi(90);
284 g = a.DrawCopy("surf1pol");
285 pad.cd(2);
286 gPad->SetLogy();
287 gPad->SetLogz();
288 gPad->SetTheta(70);
289 gPad->SetPhi(90);
290 g->Draw("surf1pol");
291 pad.cd(3);
292 gPad->SetLogy();
293 gPad->SetLogz();
294 gPad->SetTheta(90);
295 gPad->SetPhi(90);
296 g->Draw("surf1pol");
297 pad.cd(4);
298 gPad->SetLogy();
299 gPad->SetLogz();
300 gPad->SetTheta(20);
301 gPad->SetPhi(90);
302 g->Draw("surf1pol");
303
304 c->cd(3);
305 gPad->SetLogy();
306 g=r.ProjectionY();
307 g->SetDirectory(NULL);
308 TH1 *g2=r2.ProjectionY();
309 g->SetMinimum(GetMin(g2)/2);
310 g->SetBit(kCanDelete);
311 g->SetTitle(" Hit Observers Plain ");
312 g->GetXaxis()->SetLabelOffset(0);
313 g->GetXaxis()->SetTitleOffset(1.1);
314 g->GetYaxis()->SetTitleOffset(1.3);
315 g->SetXTitle("R [kpc]");
316 g->SetYTitle("Counts");
317 g->Draw();
318 g2->SetDirectory(NULL);
319 g2->SetBit(kCanDelete);
320 g2->SetLineColor(kBlue);
321 g2->Draw("same");
322 g=r3.ProjectionY();
323 g->SetDirectory(NULL);
324 g->SetBit(kCanDelete);
325 g->SetLineColor(kGreen);
326 g->Draw("same");
327
328 c->cd(4);
329 gPad->SetLogy();
330 g=a.ProjectionY();
331 g->SetMinimum(GetMin(g)/2);
332 g->SetDirectory(NULL);
333 g->SetBit(kCanDelete);
334 g->SetTitle("Hit Angle");
335 g->GetXaxis()->SetLabelOffset(0);
336 g->GetXaxis()->SetTitleOffset(1.1);
337 g->GetYaxis()->SetTitleOffset(1.3);
338 g->SetXTitle("\\Phi [']");
339 g->SetYTitle("Counts");
340 g->Draw();
341 g=a2.ProjectionY();
342 g->SetDirectory(NULL);
343 g->SetBit(kCanDelete);
344 g->SetLineColor(kBlue);
345 g->Draw("same");
346 g=a3.ProjectionY();
347 g->SetDirectory(NULL);
348 g->SetBit(kCanDelete);
349 g->SetLineColor(kGreen);
350 g->Draw("same");
351
352 // delete gROOT->FindObject("Analysis Photons");
353 TCanvas *c = MH::MakeDefCanvas("Analysis Photons", "", 580, 870);
354 c->Divide(1,2);
355
356 c->cd(1);
357 gPad->SetLogx();
358 gPad->SetLogy();
359 h.SetFillStyle(0);
360 h.SetName("Photons");
361 h.SetTitle(Form(" E^{%.1f} Spectra ", alpha));
362 h.SetXTitle("E\\gamma [GeV]");
363 h.SetYTitle(Form("E^{%.1f} * Counts", plot));
364 h.GetXaxis()->SetLabelOffset(-0.015);
365 h.GetXaxis()->SetTitleOffset(1.1);
366 h.SetMarkerStyle(kMultiply);
367 prim.SetFillStyle(0);
368 prim.SetMarkerStyle(kPlus);
369 prim.SetMarkerColor(kRed);
370 prim.SetLineColor(kRed);
371 TH1 &g1=*h.DrawCopy("P");
372 g2=*prim.DrawCopy("Psame");
373 g2->Draw("Csame");
374 g1.Draw("Csame");
375 h2.SetFillStyle(0);
376 h2.SetName("SecPhotons");
377 h2.SetMarkerStyle(kMultiply);
378 h2.SetMarkerColor(kBlue);
379 h2.SetLineColor(kBlue);
380 TH1 &g3=*h2.DrawCopy("Psame");
381 g3.Draw("Csame");
382 h3.SetFillStyle(0);
383 h3.SetName("PrimPhotons");
384 h3.SetMarkerStyle(kMultiply);
385 h3.SetMarkerColor(kGreen);
386 h3.SetLineColor(kGreen);
387 TH1 &g4=*h3.DrawCopy("Psame");
388 g4.Draw("Csame");
389
390 draw1h1426(prim.GetBinContent(2),plot);
391
392 c->cd(2);
393 gPad->SetLogx();
394 TH1D div;
395 MH::SetBinning(&div, &binsx);
396 div.Divide(&h, &prim);
397 div.SetTitle("Spectrum / Primary Spectrum");
398 div.GetXaxis()->SetLabelOffset(-0.015);
399 div.GetXaxis()->SetTitleOffset(1.1);
400 div.SetXTitle("E\\gamma [GeV]");
401 div.SetYTitle("Ratio");
402 line.SetLineColor(kBlue);
403 line.SetLineWidth(1);
404 div.DrawCopy();
405 line.DrawLine(log10(lo), 1, log10(hi), 1);
406}
Note: See TracBrowser for help on using the repository browser.