1 | //Double_t kRad2Deg = 180./TMath::Pi();
|
---|
2 |
|
---|
3 | #include <float.h>
|
---|
4 |
|
---|
5 | Double_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 |
|
---|
16 | void 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 |
|
---|
58 | void 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 |
|
---|
88 | void 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 | }
|
---|