1 | #include <iostream.h>
|
---|
2 |
|
---|
3 | #include <fstream.h>
|
---|
4 |
|
---|
5 | #include <TStopwatch.h>
|
---|
6 | #include <TF1.h>
|
---|
7 | #include <TH1.h>
|
---|
8 | #include <TH2.h>
|
---|
9 | #include <TList.h>
|
---|
10 | #include <TFile.h>
|
---|
11 | #include <TTree.h>
|
---|
12 | #include <TTimer.h>
|
---|
13 | #include <TStyle.h>
|
---|
14 | #include <TBranch.h>
|
---|
15 | #include <TRandom.h>
|
---|
16 | #include <TCanvas.h>
|
---|
17 |
|
---|
18 | #include "mbase/MParContainer.h"
|
---|
19 | #include "mphys/MPhoton.h"
|
---|
20 | #include "mphys/MElectron.h"
|
---|
21 | #include "mphys/MPairProduction.h"
|
---|
22 | #include "mhist/MBinning.h"
|
---|
23 | #include "mhist/MH.h"
|
---|
24 |
|
---|
25 | // 2.96
|
---|
26 | // 2.87
|
---|
27 | // 2.73
|
---|
28 | Double_t PrimSpect(Double_t *x, Double_t *k)
|
---|
29 | {
|
---|
30 | return pow(pow(10, x[0]), k[0]);
|
---|
31 | }
|
---|
32 |
|
---|
33 | Double_t PhotonSpect(Double_t *x, Double_t *k=NULL)
|
---|
34 | {
|
---|
35 | Double_t Ep = pow(10, x[0]);
|
---|
36 |
|
---|
37 | Double_t res = MPhoton::Int2(&Ep, k);
|
---|
38 | return res*1e55; //65/k[0];
|
---|
39 |
|
---|
40 | //return MPhoton::Planck(&Ep, &k[1]);
|
---|
41 | }
|
---|
42 |
|
---|
43 | Double_t Sbar_sigmas(Double_t *x, Double_t *k)
|
---|
44 | {
|
---|
45 | Double_t sbar = pow(10, x[0]);
|
---|
46 |
|
---|
47 | Double_t s = 1./(sbar*4);
|
---|
48 |
|
---|
49 | Double_t sigma = MPhoton::Sigma_gg(&s);
|
---|
50 |
|
---|
51 | return sigma*sbar*1e28;
|
---|
52 | }
|
---|
53 |
|
---|
54 | Double_t RandomTheta(Double_t Eg, Double_t Ep)
|
---|
55 | {
|
---|
56 | Double_t E0 = 511e-6; // [GeV]
|
---|
57 |
|
---|
58 | Double_t f = Eg/E0*Ep/E0;
|
---|
59 |
|
---|
60 | if (f<1)
|
---|
61 | return 0;
|
---|
62 |
|
---|
63 | TF1 func("RndTheta", Sbar_sigmas, 0, log10(f), 0);
|
---|
64 |
|
---|
65 | func.SetNpx(50);
|
---|
66 |
|
---|
67 | Double_t sbar = pow(10, func.GetRandom());
|
---|
68 | Double_t theta = acos(1.-sbar*2/f);
|
---|
69 |
|
---|
70 | return theta;
|
---|
71 | }
|
---|
72 |
|
---|
73 | Double_t GetEnergy(Int_t n, Double_t lo, Double_t hi)
|
---|
74 | {
|
---|
75 | static int bin=0;
|
---|
76 | static TRandom rnd(0);
|
---|
77 |
|
---|
78 | Double_t w = log10(hi/lo)/n;
|
---|
79 |
|
---|
80 | Double_t E = lo*pow(10, rnd.Uniform(w) + w*(n-bin-1));
|
---|
81 |
|
---|
82 | //cout << endl << w << " " << n << " " << w*bin << " " << E << endl;
|
---|
83 |
|
---|
84 | if (++bin==n)
|
---|
85 | bin=0;
|
---|
86 |
|
---|
87 | return E;
|
---|
88 |
|
---|
89 | /*
|
---|
90 | // TF1 src ("PrimarySpectrum", PrimSpect, log10(lo), log10(hi), 1); // 10GeV-1000PeV
|
---|
91 | // src.SetParameter(0, 0);
|
---|
92 | // Double_t E = pow(10, src.GetRandom());
|
---|
93 |
|
---|
94 | //
|
---|
95 | // 0: 11111111111111|22222222222222 2 2^0 + 1 1
|
---|
96 | // 1: 11111111|22222222222|33333333 3 2^1 + 1 2
|
---|
97 | // 2: 11111|22222|33333|44444|55555 5 2^2 + 1 7
|
---|
98 | // 3: | | | | | | | | 15
|
---|
99 | //
|
---|
100 |
|
---|
101 | static int run=0;
|
---|
102 | static int num=1;
|
---|
103 |
|
---|
104 | Double_t w = log10(hi/lo)/((1<<run) + 1);
|
---|
105 |
|
---|
106 | Double_t E = lo*pow(10, num*w);
|
---|
107 |
|
---|
108 | if (num == (1<<run))
|
---|
109 | {
|
---|
110 | run++;
|
---|
111 | num=1;
|
---|
112 | }
|
---|
113 | else
|
---|
114 | num++;
|
---|
115 |
|
---|
116 | return E;
|
---|
117 | */
|
---|
118 | }
|
---|
119 |
|
---|
120 | void DoIt()
|
---|
121 | {
|
---|
122 | // a --> 5
|
---|
123 | // b --> 10
|
---|
124 | // c --> 2
|
---|
125 |
|
---|
126 | // delme, delmeB starting integration at lolim
|
---|
127 | // delme2, delme2B starting integration at -log10(Eg)-8
|
---|
128 | // delme3, lolim, 24, 1e2-1e6, 2x inv.Compton
|
---|
129 | // delme3B, lolim, 24, 1e2-1e6, 4x inv.Compton
|
---|
130 | // delme3C, lolim, 24, 1e2-1e6, 8x inv.Compton
|
---|
131 | // delme3D, lolim, 24, 1e2-1e6, 16x inv.Compton
|
---|
132 | // delme3E, lolim, 24, 1e2-1e6, 32x inv.Compton
|
---|
133 |
|
---|
134 | // delme3H, lolim, 24, 1e2-1e6, 256x inv.Compton 8*60
|
---|
135 |
|
---|
136 | Double_t startz = 0.03; //MParticle::ZofR(&R);
|
---|
137 | Double_t R = MParticle::RofZ(&startz); // [kpc]
|
---|
138 |
|
---|
139 | const char *filename = "cascade_0.03_18_1e2_1e5_B0_256_1.root";
|
---|
140 | //const char *filename = "test.root";
|
---|
141 |
|
---|
142 | const Double_t B = 0;//1e-6*1e-4; // [T=1e4G] mean magnetic field
|
---|
143 |
|
---|
144 | Double_t runtime = 45*60; // [s] maximum time to run the simulation
|
---|
145 |
|
---|
146 | Int_t nbins = 18; // number of bins produced in energy spectrum
|
---|
147 |
|
---|
148 | Double_t lo = 1e2; // lower limit of displayed spectrum
|
---|
149 | Double_t hi = 1e5; // upper limit of spectrum (cutoff)
|
---|
150 |
|
---|
151 | Int_t counter = 256; // maximum number of inv. Compton (-1 means infinite)
|
---|
152 |
|
---|
153 | Double_t alpha = -1; // -1 means a E^2 plot
|
---|
154 | Double_t plot = 1; // 1 means a E^-2 spectrum
|
---|
155 |
|
---|
156 | Double_t bubbler = R-1e3*10; // [kpc]
|
---|
157 | Double_t bubblez = MParticle::ZofR(&bubbler);
|
---|
158 |
|
---|
159 | // --------------------------------------------------------------------
|
---|
160 |
|
---|
161 | cout << "R = " << R << "kpc" << endl;
|
---|
162 | cout << "Z = " << startz << endl;
|
---|
163 |
|
---|
164 | cout << "Setting up: Histograms... " << flush;
|
---|
165 |
|
---|
166 | MPairProduction pair;
|
---|
167 |
|
---|
168 | TH1D hist;
|
---|
169 | TH1D histsrc;
|
---|
170 |
|
---|
171 | hist.SetMinimum(pow(lo, plot+alpha)/100);
|
---|
172 | histsrc.SetMinimum(pow(lo, plot+alpha)/100);
|
---|
173 |
|
---|
174 | TList listg;
|
---|
175 | TList liste;
|
---|
176 |
|
---|
177 | listg.SetOwner();
|
---|
178 | liste.SetOwner();
|
---|
179 |
|
---|
180 | gStyle->SetOptStat(10);
|
---|
181 |
|
---|
182 | //
|
---|
183 | // Don't change the order!!!
|
---|
184 | //
|
---|
185 | hist.SetFillStyle(0);
|
---|
186 | hist.SetMarkerStyle(kPlus);
|
---|
187 | histsrc.SetFillStyle(0);
|
---|
188 | histsrc.SetMarkerStyle(kMultiply);
|
---|
189 | histsrc.SetMarkerColor(kRed);
|
---|
190 | histsrc.SetLineColor(kRed);
|
---|
191 |
|
---|
192 | MBinning bins;
|
---|
193 | bins.SetEdgesLog(nbins, lo, hi);
|
---|
194 |
|
---|
195 | MH::SetBinning(&hist, &bins);
|
---|
196 | MH::SetBinning(&histsrc, &bins);
|
---|
197 |
|
---|
198 | TCanvas *c=MH::MakeDefCanvas();
|
---|
199 |
|
---|
200 | gPad->SetGrid();
|
---|
201 | gPad->SetLogx();
|
---|
202 | gPad->SetLogy();
|
---|
203 |
|
---|
204 | hist.SetName("Spectrum");
|
---|
205 | hist.SetXTitle("E [GeV]");
|
---|
206 | hist.SetYTitle(Form("E^{%.1f} Counts", plot));
|
---|
207 | hist.GetXaxis()->SetLabelOffset(-0.015);
|
---|
208 | hist.GetXaxis()->SetTitleOffset(1.1);
|
---|
209 |
|
---|
210 | hist.Draw("P");
|
---|
211 | histsrc.Draw("Psame");
|
---|
212 | histsrc.Draw("same");
|
---|
213 | hist.Draw("same");
|
---|
214 |
|
---|
215 | // ------------------------------
|
---|
216 |
|
---|
217 | cout << "Output File... " << flush;
|
---|
218 |
|
---|
219 | TFile file(filename, "RECREATE", "Intergalactic cascade", 9);
|
---|
220 | cout << "Trees... " << flush;
|
---|
221 | TTree *T1 = new TTree ("Photons", "Photons from Cascade");
|
---|
222 | TTree *T2 = new TTree ("Electrons", "Electrons in the Cascade");
|
---|
223 | cout << "Branches... " << flush;
|
---|
224 | MPhoton dummyp;
|
---|
225 | void *ptr = &dummyp;
|
---|
226 | TBranch *B1 = T1->Branch("MPhoton.", "MPhoton", &ptr);
|
---|
227 | MPhoton dummye;
|
---|
228 | ptr = &dummye;
|
---|
229 | TBranch *B2 = T2->Branch("MElectron.", "MElectron", &ptr);
|
---|
230 |
|
---|
231 | // ------------------------------
|
---|
232 |
|
---|
233 | cout << "Timers... " << flush;
|
---|
234 |
|
---|
235 | TTimer timer("gSystem->ProcessEvents();", 250, kFALSE);
|
---|
236 | timer.TurnOn();
|
---|
237 |
|
---|
238 | TStopwatch clock;
|
---|
239 | clock.Start();
|
---|
240 |
|
---|
241 | cout << "Done. " << endl;
|
---|
242 |
|
---|
243 | Int_t n=0;
|
---|
244 | Double_t starttime = TStopwatch::GetRealTime(); // s
|
---|
245 | while (TStopwatch::GetRealTime()<starttime+runtime)
|
---|
246 | for (int i=0; i<nbins; i++)
|
---|
247 | {
|
---|
248 | n++;
|
---|
249 |
|
---|
250 | Double_t E = GetEnergy(nbins, lo, hi);
|
---|
251 | Double_t weight = pow(E, alpha);
|
---|
252 | histsrc.Fill(E, pow(E, plot) * weight);
|
---|
253 |
|
---|
254 | MPhoton *gamma=new MPhoton(E, startz);
|
---|
255 | gamma->SetIsPrimary();
|
---|
256 | gamma->SetSrcR(R);
|
---|
257 | gamma->InitRandom();
|
---|
258 | listg.Add(gamma);
|
---|
259 |
|
---|
260 | B1->SetAddress(&gamma);
|
---|
261 | T1->Fill();
|
---|
262 |
|
---|
263 | cout << "--> " << n << ". " << Form("%d: %3.1e", (int)(starttime+runtime-TStopwatch::GetRealTime()), E) << flush;
|
---|
264 |
|
---|
265 | Double_t Esum=0;
|
---|
266 | Double_t Emis=0;
|
---|
267 | while (1)
|
---|
268 | {
|
---|
269 | if (listg.GetSize()!=0)
|
---|
270 | cout << " |P" << flush;
|
---|
271 |
|
---|
272 | TIter NextP(&listg);
|
---|
273 | MPhoton *p = NULL;
|
---|
274 | B1->SetAddress(&p);
|
---|
275 | while ((p=(MPhoton*)NextP()))
|
---|
276 | {
|
---|
277 | cout << ":" << flush;
|
---|
278 | Double_t Eg = p->GetEnergy();
|
---|
279 |
|
---|
280 | Double_t E0 = 511e-6;
|
---|
281 | Double_t z = p->GetZ();
|
---|
282 | Double_t lolim = E0*E0/Eg;
|
---|
283 | Double_t inf = (Eg<1e6 ? 3e-11*(z+1) : 3e-12*(z+1));
|
---|
284 | if (Eg<5e4)
|
---|
285 | inf = 3e-11*(z+1)*pow(10, 9.4-log10(Eg)*2);
|
---|
286 |
|
---|
287 | //TF1 phot("PhotonSpectrum", PhotonSpect, -log10(Eg)-8, -10.5, 2);
|
---|
288 | //Double_t E0 = 511e-6; // [GeV]
|
---|
289 | TF1 phot("PhotonSpectrum", PhotonSpect, log10(lolim), log10(inf), 2);
|
---|
290 | phot.SetParameter(0, Eg);
|
---|
291 | while (1)
|
---|
292 | {
|
---|
293 | if (!p->SetNewPosition())
|
---|
294 | {
|
---|
295 | cout << "!" << flush;
|
---|
296 |
|
---|
297 | hist.Fill(Eg, pow(Eg, plot)*weight);
|
---|
298 |
|
---|
299 | p->SetIsPrimary(kFALSE);
|
---|
300 | T1->Fill();
|
---|
301 |
|
---|
302 | Esum += Eg;
|
---|
303 |
|
---|
304 | break;
|
---|
305 | }
|
---|
306 |
|
---|
307 | //
|
---|
308 | // Sample phtoton from background and interaction angle
|
---|
309 | //
|
---|
310 | phot.SetParameter(1, p->GetZ());
|
---|
311 | Double_t pe = phot.GetRandom();
|
---|
312 | if (pe==0)
|
---|
313 | {
|
---|
314 | cout << "z" << flush;
|
---|
315 | continue;
|
---|
316 | }
|
---|
317 |
|
---|
318 | Double_t Ep = pow(10, pe);
|
---|
319 | Double_t theta = RandomTheta(Eg, Ep);
|
---|
320 | if (theta==0)
|
---|
321 | {
|
---|
322 | cout << "t" << flush;
|
---|
323 | continue;
|
---|
324 | }
|
---|
325 |
|
---|
326 | if (!pair.Process(p, Ep, theta, &liste))
|
---|
327 | {
|
---|
328 | cout << "0" << flush;
|
---|
329 | continue;
|
---|
330 | }
|
---|
331 |
|
---|
332 | cout << "." << flush;
|
---|
333 | break;
|
---|
334 | }
|
---|
335 | delete listg.Remove(p);
|
---|
336 | }
|
---|
337 |
|
---|
338 | if (liste.GetSize()==0 && listg.GetSize()==0)
|
---|
339 | break;
|
---|
340 |
|
---|
341 | if (liste.GetSize()!=0)
|
---|
342 | cout << " |E" << flush;
|
---|
343 |
|
---|
344 | TIter Next(&liste);
|
---|
345 | MElectron *e = NULL;
|
---|
346 | B2->SetAddress(&e);
|
---|
347 | while ((e=(MElectron*)Next()))
|
---|
348 | {
|
---|
349 | e->SetIsPrimary(kTRUE);
|
---|
350 | T2->Fill();
|
---|
351 | e->SetIsPrimary(kFALSE);
|
---|
352 |
|
---|
353 | Double_t Ee = e->GetEnergy();
|
---|
354 |
|
---|
355 | if (Ee<lo)
|
---|
356 | continue;
|
---|
357 |
|
---|
358 | cout << ":" << flush;
|
---|
359 | int test = counter<0 ? -1 : 0;
|
---|
360 | while (test<0 ? true : (test++<counter))
|
---|
361 | {
|
---|
362 | if (!e->SetNewPositionB(e->GetZ()>bubblez ? B : 0))
|
---|
363 | {
|
---|
364 | cout << "!" << flush;
|
---|
365 | break;
|
---|
366 | }
|
---|
367 |
|
---|
368 | // WRONG!
|
---|
369 | MPhoton *p = e->DoInvCompton(0);
|
---|
370 |
|
---|
371 | T2->Fill();
|
---|
372 |
|
---|
373 | cout << "." << flush;
|
---|
374 | listg.Add(p);
|
---|
375 | /*
|
---|
376 | if (fabs(p->GetTheta()*3437)<60 && // < 60min
|
---|
377 | p->GetEnergy()>lo)
|
---|
378 | {
|
---|
379 | cout << "." << flush;
|
---|
380 | listg.Add(p);
|
---|
381 | }
|
---|
382 | else
|
---|
383 | {
|
---|
384 | if (p->GetEnergy()<E*pow(10, alpha)/5 || p->GetEnergy()<=lo)
|
---|
385 | cout << "e" << flush;
|
---|
386 | else
|
---|
387 | cout << "t" << flush; // ignored
|
---|
388 | delete p;
|
---|
389 | Esum += p->GetEnergy();
|
---|
390 | Emis += p->GetEnergy();
|
---|
391 | }
|
---|
392 | */
|
---|
393 | if (fabs(e->GetTheta()*3437)>60) // < 60min
|
---|
394 | {
|
---|
395 | cout << "T" << flush;
|
---|
396 | break;
|
---|
397 | }
|
---|
398 |
|
---|
399 | if (e->GetEnergy()<Ee*1e-3) // <2e3
|
---|
400 | {
|
---|
401 | cout << "E" << flush;
|
---|
402 | break;
|
---|
403 | }
|
---|
404 | }
|
---|
405 | Esum += e->GetEnergy();
|
---|
406 | Emis += e->GetEnergy();
|
---|
407 | }
|
---|
408 | liste.Delete();
|
---|
409 | }
|
---|
410 | cout << " ----> " << Form("%3.1f %3.1e / %3.1f %3.1e", Emis/E, Emis, Esum/E, Esum) << endl;
|
---|
411 |
|
---|
412 | hist.SetTitle(Form("E^{%.1f} z=%f T=%d'%d\" N=%d", alpha, startz,
|
---|
413 | (int)runtime/60, (int)runtime%60, n));
|
---|
414 |
|
---|
415 | timer.Stop();
|
---|
416 | c->Update();
|
---|
417 | timer.Start(250);
|
---|
418 |
|
---|
419 | }
|
---|
420 | cout << endl;
|
---|
421 |
|
---|
422 | clock.Stop();
|
---|
423 | clock.Print();
|
---|
424 |
|
---|
425 | timer.Stop();
|
---|
426 |
|
---|
427 | file.Write();
|
---|
428 |
|
---|
429 | cout << "Wrote: " << filename << endl;
|
---|
430 | cout << "Created " << n << " gammas from source with E^" << alpha << endl;
|
---|
431 | cout << "Processing time: " << Form("%.1f", (TStopwatch::GetRealTime()-starttime)/n) << " sec/gamma." << endl;
|
---|
432 |
|
---|
433 | // ------------------------------
|
---|
434 |
|
---|
435 | gPad->Clear();
|
---|
436 |
|
---|
437 | hist.SetTitle(Form("E^{%.1f} z=%f T=%d'%d\" N=%d", alpha, startz, (int)runtime/60, (int)runtime%60, n));
|
---|
438 |
|
---|
439 | TH1 &h1 = *hist.DrawCopy("P");
|
---|
440 | TH1 &h2 = *histsrc.DrawCopy("Psame");
|
---|
441 | h2.Draw("Csame");
|
---|
442 | h1.Draw("Csame");
|
---|
443 | }
|
---|
444 |
|
---|
445 | // -------------------------------------------------------------------
|
---|
446 | void phys()
|
---|
447 | {
|
---|
448 | /*
|
---|
449 | Double_t Eg = 10;
|
---|
450 | Double_t E0 = 511e-6;
|
---|
451 | Double_t z = 0.03;
|
---|
452 | Double_t lolim = E0*E0/Eg;
|
---|
453 | Double_t inf = (Eg<1e6 ? 3e-11*(z+1) : 3e-12*(z+1));
|
---|
454 | if (Eg<5e4)
|
---|
455 | inf = 3e-11*(z+1)*pow(10, 9.4-log10(Eg)*2);
|
---|
456 | TF1 phot("PhotonSpectrum", PhotonSpect, log10(lolim), log10(inf), 2);
|
---|
457 | phot.SetParameter(0, Eg);
|
---|
458 | phot.SetParameter(1, z);
|
---|
459 |
|
---|
460 | Double_t val[2] = {Eg,z};
|
---|
461 | cout << phot.Integral(log10(lolim), log10(inf), val, 1e-2) << endl;
|
---|
462 |
|
---|
463 | cout << lolim << " " << inf << endl;
|
---|
464 |
|
---|
465 | phot.DrawCopy();
|
---|
466 | return;
|
---|
467 | */
|
---|
468 |
|
---|
469 | DoIt();
|
---|
470 |
|
---|
471 | /*
|
---|
472 | Double_t E = 1e10;
|
---|
473 | TF1 phot("PhotonSpectrum", PhotonSpect, -log10(E)-8, -10.5, 2);
|
---|
474 | phot.SetParameter(0, E);
|
---|
475 | phot.SetParameter(1, 5);
|
---|
476 | phot.DrawCopy();
|
---|
477 | return;
|
---|
478 | */
|
---|
479 |
|
---|
480 | /*
|
---|
481 | Double_t Eg = 1e5;
|
---|
482 |
|
---|
483 | Double_t E0 = 511e-6; // [GeV]
|
---|
484 | Double_t lolim = E0*E0/Eg;
|
---|
485 | Double_t inf = 4e-12; //E0*E0/Eg * sqrt(Eg);
|
---|
486 |
|
---|
487 | // 1e5: 5e-13, 4e-12
|
---|
488 | // 1e6: 5e-14, 2e-12
|
---|
489 | // 1e8: 5e-16, 1e-12
|
---|
490 | // 1e10: 5e-18, 1e-12
|
---|
491 |
|
---|
492 | // 1e5: 2.6e-12, 4e-12
|
---|
493 | // 1e6: 2.6e-13, 2e-12
|
---|
494 | // 1e8: 2.6e-15, 1e-12
|
---|
495 | // 1e10: 1.6e-17, 1e-12
|
---|
496 |
|
---|
497 | TF1 f("int2", MPhoton::Int2, lolim, inf, 2);
|
---|
498 | f.SetParameter(0, Eg);
|
---|
499 | f.SetParameter(1, 0);
|
---|
500 |
|
---|
501 | MH::MakeDefCanvas();
|
---|
502 | gPad->SetLogx();
|
---|
503 | f.DrawCopy();
|
---|
504 | */
|
---|
505 | }
|
---|
506 |
|
---|