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

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