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

Last change on this file since 1449 was 1449, checked in by tbretz, 22 years ago
*** empty log message ***
File size: 15.5 KB
Line 
1//Double_t kRad2Deg = 180./TMath::Pi();
2
3#include <float.h>
4#include <iomanip.h>
5
6#include <TH1.h>
7#include <TH2.h>
8#include <TLine.h>
9#include <TLeaf.h>
10#include <TChain.h>
11#include <TGaxis.h>
12#include <TStyle.h>
13#include <TCanvas.h>
14#include <TPolyLine.h>
15#include <TPolyMarker.h>
16
17#include "MH.h"
18#include "MBinning.h"
19
20#include "MFit.h"
21#include "MPhoton.h"
22
23void GetRange(TChain *chain, const char *name, Int_t &nbins, Double_t &min, Double_t &max, Double_t conv=1, Bool_t zero=kTRUE)
24{
25 TString str("MPhoton.MParticle.");
26 str += name;
27
28 MPhoton *p = new MPhoton;
29
30 chain->SetBranchAddress("MPhoton.", &p);
31 chain->SetBranchStatus("*", 0);
32 chain->SetBranchStatus(str, 1);
33
34 min = FLT_MAX;
35 max = -FLT_MAX;
36
37 Int_t numt = chain->GetTreeNumber();
38 TLeaf *leaf = chain->GetLeaf(str);
39
40 if (!leaf && numt>=0)
41 return;
42
43 Stat_t n = chain->GetEntries();
44 for (int i=0; i<n; i++)
45 {
46 chain->GetEntry(i);
47
48 if (numt != chain->GetTreeNumber())
49 {
50 if (!(leaf = chain->GetLeaf(str)))
51 return;
52
53 numt = chain->GetTreeNumber();
54 }
55
56 Double_t val = leaf->GetValue();
57
58 if (p->IsPrimary())
59 continue;
60
61 if (val < min && (zero || val!=0))
62 min = val;
63 if (val > max)
64 max = val;
65 }
66
67 min *= conv;
68 max *= conv;
69
70 Int_t newn=0;
71 MH::FindGoodLimits(nbins, newn, min, max, kFALSE);
72 nbins = newn;
73}
74
75void draw1h1426(Double_t E, Double_t plot=1)
76{
77 plot += 1;
78
79 const Double_t level = E / (5.73e-0*pow( 0.39e3, plot));
80
81 TPolyLine *l = new TPolyLine(6);
82 TPolyMarker *m = new TPolyMarker(6);
83
84 // m->SetPoint(0, log10( 0.28e3), log10(1.86e-0*pow( 0.28e3, plot)*level));
85 m->SetPoint(0, log10( 0.39e3)-0.3, log10(5.73e-0*pow( 0.39e3, plot)*level));
86 m->SetPoint(1, log10( 0.80e3)-0.3, log10(1.22e-0*pow( 0.80e3, plot)*level));
87 m->SetPoint(2, log10( 1.55e3)-0.3, log10(5.10e-2*pow( 1.55e3, plot)*level));
88 m->SetPoint(3, log10( 2.82e3)-0.3, log10(1.50e-2*pow( 2.82e3, plot)*level));
89 m->SetPoint(4, log10( 5.33e3)-0.3, log10(7.70e-3*pow( 5.33e3, plot)*level));
90 m->SetPoint(5, log10(10.00e3)-0.3, log10(1.20e-4*pow(10.00e3, plot)*level));
91 m->SetMarkerStyle(kOpenStar);
92 m->SetMarkerColor(kMagenta);
93 m->SetBit(kCanDelete);
94 m->Draw();
95
96 // l->SetPoint(0, log10( 0.28e3), log10(1.86e-0*pow( 0.28e3, plot)*level));
97 l->SetPoint(0, log10( 0.39e3)-0.3, log10(5.73e-0*pow( 0.39e3, plot)*level));
98 l->SetPoint(1, log10( 0.80e3)-0.3, log10(1.22e-0*pow( 0.80e3, plot)*level));
99 l->SetPoint(2, log10( 1.55e3)-0.3, log10(5.10e-2*pow( 1.55e3, plot)*level));
100 l->SetPoint(3, log10( 2.82e3)-0.3, log10(1.50e-2*pow( 2.82e3, plot)*level));
101 l->SetPoint(4, log10( 5.33e3)-0.3, log10(7.70e-3*pow( 5.33e3, plot)*level));
102 l->SetPoint(5, log10(10.00e3)-0.3, log10(1.20e-4*pow(10.00e3, plot)*level));
103 l->SetLineColor(kMagenta);
104 l->SetBit(kCanDelete);
105 l->Draw();
106}
107
108void analp()
109{
110 // -------------------------------------------------------------------
111
112 TChain chain("Photons");
113 chain.Add("cascade_0.1_18_1e2_1e5_B0_256_01.root");
114 chain.Add("cascade_0.1_18_1e2_1e5_B0_256_02.root");
115 chain.Add("cascade_0.1_18_1e2_1e5_B0_256_03.root");
116 chain.Add("cascade_0.1_18_1e2_1e5_B0_256_04.root");
117 chain.Add("cascade_0.1_18_1e2_1e5_B0_512_01.root");
118 chain.Add("cascade_0.1_18_1e2_1e5_B0_512_03.root");
119 chain.Add("cascade_0.1_18_1e2_1e5_B0_512_02.root");
120 /*
121 chain.Add("cascade_0.03_18_1e2_1e5_B0_256_3.root");
122 chain.Add("cascade_0.03_18_1e2_1e5_B0_256_4.root");
123 chain.Add("cascade_0.03_18_1e2_1e5_B0_256_5.root");
124 chain.Add("cascade_0.03_18_1e2_1e5_B0_256_6.root");
125 chain.Add("cascade_0.03_18_1e2_1e5_B0_256_7.root");
126 chain.Add("cascade_0.03_18_1e2_1e5_B0_256_8.root");
127 chain.Add("cascade_0.03_18_1e2_1e5_B0_256_9.root");
128
129 chain.Add("cascade_0.03_18_1e2_1e5_B0_512_01.root");
130 chain.Add("cascade_0.03_18_1e2_1e5_B0_512_02.root");
131 */
132
133 Int_t nbinse = 18; // number of bins in your histogram
134 Double_t lo = 1e2/1e3; // lower limit of the energy histogram
135 Double_t hi = 1e5; // upper limit of the energy histogram
136
137 Double_t minE = 1e2; // minimum value produced in the simulation
138
139 Double_t cutoff = 1e12; // arbitrary energy cutoff (throw away the events)
140
141 Double_t alpha = -1; // alpha of the energy spectrum (-1 is E^-2)
142 Double_t plot = 1; // plot of the enrgy spectrum (1 is E^2)
143
144 Int_t nbins = 16; // number of bins you want to have in psi/phi
145
146 // -------------------------------------------------------------------
147
148 MBinning binspolx;
149 binspolx.SetEdges(16, -180, 180);
150
151 // -------------------------------------------------------------------
152
153 Stat_t n = chain.GetEntries();
154
155 cout << endl << "Found " << n << " entries." << endl;
156
157 if (n==0)
158 return;
159
160 // -----------------------------------
161
162 MBinning bins;
163
164 cout << "Determin min/max..." << endl;
165
166 Double_t min, max;
167 GetRange(&chain, "fTheta", nbins, min, max, kRad2Deg*60);
168 cout << "fTheta, " << nbins << ": " << min << " ... " << max << " [']" << endl;
169
170 bins.SetEdges(nbins, min, max*1.1);
171
172 TH2D a, a2, a3;
173 MH::SetBinning(&a, &binspolx, &bins);
174 MH::SetBinning(&a2, &binspolx, &bins);
175 MH::SetBinning(&a3, &binspolx, &bins);
176
177 // -----------------------------------
178
179 GetRange(&chain, "fR", nbins, min, max);
180 cout << "fR, " << nbins << ": " << min << " ... " << max << " [kpc]" << endl;
181
182 bins.SetEdges(nbins, min, max*1.1);
183
184 TH2D r, r2, r3;
185 MH::SetBinning(&r, &binspolx, &bins);
186 MH::SetBinning(&r2, &binspolx, &bins);
187 MH::SetBinning(&r3, &binspolx, &bins);
188
189 // -----------------------------------
190
191 MBinning binsx;
192 binsx.SetEdgesLog(nbinse, lo, hi);
193
194 TH1D prim, h, h2, h3;
195 MH::SetBinning(&h, &binsx);
196 MH::SetBinning(&h2, &binsx);
197 MH::SetBinning(&h3, &binsx);
198 MH::SetBinning(&prim, &binsx);
199// prim.Sumw2();
200 h.Sumw2();
201 h2.Sumw2();
202 h3.Sumw2();
203
204 // -----------------------------------
205
206 MPhoton *p = new MPhoton;
207 chain.SetBranchAddress("MPhoton.", &p);
208 chain.SetBranchStatus("*", 1);
209
210 chain.GetEntry(0);
211 if (!p->IsPrimary())
212 return;
213
214 Double_t z = p->GetZ();
215 Double_t R = MParticle::RofZ(&z);
216
217 cout << "Z = " << z << endl;
218 cout << "R = " << R << "kpc" << endl;
219
220 // -----------------------------------
221
222 const Double_t skpc = 3600.*24.*365.*3.258*1e3; // [s/kpc]
223 GetRange(&chain, "fX", nbins, min, max, 1, kFALSE);
224 min *= skpc;
225 max *= skpc;
226 if (min<1e-2)
227 min=1e-2;
228 bins.SetEdgesLog(nbins, min, max);
229 cout << "dX, " << nbins << ": " << min << " ... " << max << " [s]" << endl;
230
231 TH1D t;
232 MH::SetBinning(&t, &bins);
233 t.Sumw2();
234
235 // ----------------------------------------------------------------------
236
237 chain.SetBranchAddress("MPhoton.", &p);
238 chain.SetBranchStatus("*", 1);
239
240 cout << "Filling: " << nbinse << " bins @ " << minE << " - " << hi << "... " << flush;
241
242 Double_t weight = 0;
243
244 Bool_t isprim = kFALSE;
245 Bool_t ignore = kFALSE;
246 for (int i=0; i<n; i++)
247 {
248 chain.GetEntry(i);
249
250 Double_t Ep = p->GetEnergy();
251
252 if (p->IsPrimary())
253 {
254 if (Ep>cutoff)
255 {
256 ignore = kTRUE;
257 continue;
258 }
259 ignore = kFALSE;
260 weight = pow(Ep, alpha);
261 prim.Fill(Ep, pow(Ep, plot+alpha));
262 isprim = kTRUE;
263 continue;
264 }
265
266 if (ignore)
267 continue;
268
269 Double_t phi = fmod(p->GetPhi()*kRad2Deg+180, 360)-180;
270 Double_t psi = fmod(p->GetPsi()*kRad2Deg+180, 360)-180;
271
272 if (isprim)
273 {
274 h3.Fill(Ep, pow(Ep,plot) * weight);
275 r3.Fill(phi, 0, weight);
276 a3.Fill(psi, 0, weight);
277 isprim = kFALSE;
278 continue;
279 }
280
281 h2.Fill(Ep, pow(Ep,plot) * weight);
282 r2.Fill(phi, p->GetR(), weight);
283 a2.Fill(psi, p->GetTheta()*kRad2Deg*60, weight);
284
285 Double_t tm1 = sqrt(R*R+p->GetR()*p->GetR());
286 Double_t tm2 = (R+p->GetX())/tm1 - 1;
287 t.Fill(p->GetX()<=0 ? 0 : R*tm2*skpc, weight);
288 }
289
290 // ----------------------------------------------------------------------
291
292 delete p;
293
294 cout << "Done." << endl;
295
296 Double_t N = prim.GetBinContent(nbinse);
297 Double_t Nerr = prim.GetBinError(nbinse);
298 for (int i=1; i<=nbinse; i++)
299 {
300 Double_t E = prim.GetXaxis()->GetBinCenter(i);
301 if (E>minE)
302 continue;
303
304 Double_t E1 = prim.GetXaxis()->GetBinLowEdge(i+1);
305 Double_t E0 = prim.GetXaxis()->GetBinLowEdge(i);
306
307 E = sqrt(E1*E0);
308
309 prim.SetBinContent(i, N);
310 h3.SetBinContent(i, N);
311 prim.SetBinError(i, Nerr);
312 h3.SetBinError(i, Nerr);
313 }
314
315 h.Add(&h2, &h3);
316 r.Add(&r2, &r3);
317 a.Add(&a2, &a3);
318
319 if (alpha==-plot)
320 {
321 cout << "Observed Energy: " << h.Integral() << "GeV" << endl;
322 cout << "Emitted Energy: " << prim.Integral() << "GeV" << endl;
323 cout << "Ratio: " << h.Integral()/prim.Integral() << endl;
324 }
325
326 cout << "Mean fPhi: " << r.GetMean() << " +- " << r.GetRMS() << endl;
327 cout << "Mean fPsi: " << a.GetMean() << " +- " << a.GetRMS() << endl;
328
329 gStyle->SetOptStat(10);
330 gStyle->SetPalette(1, 0);
331
332 TLine line;
333
334 // ----------------------------------------------------------------------
335
336// delete gROOT->FindObject("Analysis Arrival");
337 TCanvas *c = MH::MakeDefCanvas("Analysis Arrival", "");
338 c->Divide(2,2);
339
340 c->cd(1);
341 r.SetTitle(" Distance from Observer ");
342 r.GetXaxis()->SetLabelOffset(-0.015);
343 r.GetXaxis()->SetTitleOffset(1.1);
344 r.GetXaxis()->SetRangeUser(1e4, 1e9);
345 r.SetXTitle("\\Phi [\\circ]");
346 r.SetYTitle("R [kpc]");
347 TVirtualPad &pad = *gPad;
348 pad.Divide(2,2);
349 pad.cd(1);
350 gPad->SetLogy();
351 gPad->SetLogz();
352 gPad->SetTheta(0);
353 gPad->SetPhi(90);
354 TH1 *g = r.DrawCopy("surf1pol");
355 pad.cd(2);
356 gPad->SetLogy();
357 gPad->SetLogz();
358 gPad->SetTheta(70);
359 gPad->SetPhi(90);
360 g->Draw("surf1pol");
361 pad.cd(3);
362 gPad->SetLogy();
363 gPad->SetLogz();
364 gPad->SetTheta(90);
365 gPad->SetPhi(90);
366 g->Draw("surf1pol");
367 pad.cd(4);
368 gPad->SetLogy();
369 gPad->SetLogz();
370 gPad->SetTheta(20);
371 gPad->SetPhi(90);
372 g->Draw("surf1pol");
373
374 c->cd(2);
375 a.SetTitle(" Hit Angle for Observer ");
376 a.GetXaxis()->SetLabelOffset(-0.015);
377 a.GetXaxis()->SetTitleOffset(1.1);
378 a.GetXaxis()->SetRangeUser(1e4, 1e9);
379 a.SetXTitle("\\Psi [\\circ]");
380 a.SetYTitle("\\Theta [']");
381 TVirtualPad &pad2 = *gPad;
382 pad2.Divide(2,2);
383 pad2.cd(1);
384 gPad->SetLogy();
385 gPad->SetLogz();
386 gPad->SetTheta(0);
387 gPad->SetPhi(90);
388 g = a.DrawCopy("surf1pol");
389 pad2.cd(2);
390 gPad->SetLogy();
391 gPad->SetLogz();
392 gPad->SetTheta(70);
393 gPad->SetPhi(90);
394 g->Draw("surf1pol");
395 pad2.cd(3);
396 gPad->SetLogy();
397 gPad->SetLogz();
398 gPad->SetTheta(90);
399 gPad->SetPhi(90);
400 g->Draw("surf1pol");
401 pad2.cd(4);
402 gPad->SetLogy();
403 gPad->SetLogz();
404 gPad->SetTheta(20);
405 gPad->SetPhi(90);
406 g->Draw("surf1pol");
407
408 c->cd(3);
409 gPad->SetLogy();
410 g=r.ProjectionY();
411 g->SetDirectory(NULL);
412 TH1 *g2=r2.ProjectionY();
413 g->SetMinimum(MH::GetMinimumGT(*g2)/2);
414 g->SetBit(kCanDelete);
415 g->SetTitle(" Hit Observers Plain ");
416 g->GetXaxis()->SetLabelOffset(0);
417 g->GetXaxis()->SetTitleOffset(1.1);
418 g->GetYaxis()->SetTitleOffset(1.3);
419 g->SetXTitle("R [kpc]");
420 g->SetYTitle("Counts");
421 g->Draw();
422 g2->SetDirectory(NULL);
423 g2->SetBit(kCanDelete);
424 g2->SetLineColor(kBlue);
425 g2->Draw("same");
426 g=r3.ProjectionY();
427 g->SetDirectory(NULL);
428 g->SetBit(kCanDelete);
429 g->SetLineColor(kGreen);
430 g->Draw("same");
431
432 c->cd(4);
433 gPad->SetLogy();
434 g=a.ProjectionY();
435 g->SetMinimum(MH::GetMinimumGT(*g)/2);
436 g->SetDirectory(NULL);
437 g->SetBit(kCanDelete);
438 g->SetTitle("Hit Angle");
439 g->GetXaxis()->SetLabelOffset(0);
440 g->GetXaxis()->SetTitleOffset(1.1);
441 g->GetYaxis()->SetTitleOffset(1.3);
442 g->SetXTitle("\\Phi [']");
443 g->SetYTitle("Counts");
444 g->Draw();
445 g=a2.ProjectionY();
446 g->SetDirectory(NULL);
447 g->SetBit(kCanDelete);
448 g->SetLineColor(kBlue);
449 g->Draw("same");
450 g=a3.ProjectionY();
451 g->SetDirectory(NULL);
452 g->SetBit(kCanDelete);
453 g->SetLineColor(kGreen);
454 g->Draw("same");
455
456 // ----------------------------------------------------------------------
457
458 // delete gROOT->FindObject("Analysis Photons");
459 c = MH::MakeDefCanvas("Analysis Photons", "", 580, 870);
460 c->Divide(1,2);
461
462 c->cd(1);
463 gPad->SetLogx();
464 gPad->SetLogy();
465 h.SetFillStyle(0);
466 h.SetName("Photons");
467 h.SetTitle(Form(" E^{%.1f} Spectra ", alpha));
468 h.SetXTitle("E\\gamma [GeV]");
469 h.SetYTitle(Form("E^{%.1f} * Counts", plot));
470 h.GetXaxis()->SetLabelOffset(-0.015);
471 h.GetXaxis()->SetTitleOffset(1.1);
472 h.SetMarkerStyle(kPlus);
473 prim.SetFillStyle(0);
474 prim.SetMarkerStyle(kPlus);
475 prim.SetMarkerColor(kRed);
476 prim.SetLineColor(kRed);
477 TH1 &g1=*h.DrawCopy("P E4");
478 g2=prim.DrawCopy("Psame E0");
479 g2->Draw("Csame");
480 g1.Draw("Csame");
481 h2.SetFillStyle(0);
482 h2.SetName("SecPhotons");
483 h2.SetMarkerStyle(kMultiply);
484 h2.SetMarkerColor(kBlue);
485 h2.SetLineColor(kBlue);
486 TH1 &g3=*h2.DrawCopy("Psame E2");
487 g3.Draw("Csame");
488 h3.SetFillStyle(0);
489 h3.SetName("PrimPhotons");
490 h3.SetMarkerStyle(kMultiply);
491 h3.SetMarkerColor(kGreen);
492 h3.SetLineColor(kGreen);
493 TH1 &g4=*h3.DrawCopy("Psame E1");
494 g4.Draw("Csame");
495
496 draw1h1426(prim.GetBinContent(2),plot);
497
498 c->cd(2);
499 gPad->SetLogx();
500 TH1D div;
501 MH::SetBinning(&div, &binsx);
502 div.Sumw2();
503 div.Divide(&h, &prim);
504 div.SetTitle(" Spectrum / Primary Spectrum ");
505 div.GetXaxis()->SetLabelOffset(-0.015);
506 div.GetXaxis()->SetTitleOffset(1.1);
507 div.SetXTitle("E\\gamma [GeV]");
508 div.SetYTitle("Ratio");
509 TH1 *gHist = div.DrawCopy("E4");
510
511 line.SetLineWidth(1);
512 line.SetLineColor(kGreen);
513 line.DrawLine(log10(lo), exp(-1), log10(hi), exp(-1));
514 line.SetLineColor(kBlue);
515 line.DrawLine(log10(lo), 1, log10(hi), 1);
516
517 MFit fit("exp(-1)-[1]*log10(x/[0])");
518 for (int i=0; i<gHist->GetEntries(); i++)
519 {
520 if (gHist->GetBinContent(i+1)<exp(-1))
521 {
522 fit.SetRange(gHist->GetBinCenter(i-2), gHist->GetBinCenter(i+2));
523 cout << "Fitting from " << gHist->GetBinLowEdge(i-1) << " to ";
524 cout << gHist->GetBinLowEdge(i+1) << endl;
525 break;
526 }
527 }
528 fit.SetParameter(0, "t", 750, 10, 1e5);
529 fit.SetParameter(1, "m", 0.5, 0.1, 10);
530 fit.FitLog(gHist);
531 fit.Print();
532 fit.DrawCopy("same");
533
534 cout << "Cutoff: " << setprecision(2) << fit[0]/1e3 << "TeV +- ";
535 cout << fit.GetParError(0) << endl;
536
537 // ----------------------------------------------------------------------
538
539 c = MH::MakeDefCanvas("Time Analysis", "", 580, 435);
540 gPad->SetLogx();
541 // gPad->SetLogy();
542 t.SetName("Times");
543 t.SetTitle(" Arrival time distribution ");
544 t.SetXTitle("\\Delta t [s]");
545 t.GetXaxis()->SetLabelOffset(-0.015);
546 t.GetXaxis()->SetTitleOffset(1.1);
547 t.DrawCopy("E4");
548 line.DrawLine(log10(1), 0, log10(1), t.GetMaximum()*1.05);
549 line.DrawLine(log10(3600), 0, log10(3600), t.GetMaximum()*1.05);
550 line.DrawLine(log10(3600*24), 0, log10(3600*24), t.GetMaximum()*1.05);
551 line.DrawLine(log10(3600*24*365), 0, log10(3600*24*365), t.GetMaximum()*1.05);
552}
Note: See TracBrowser for help on using the repository browser.