source: trunk/WuerzburgSoft/Thomas/mphys/phys.C@ 1448

Last change on this file since 1448 was 1448, checked in by tbretz, 22 years ago
*** empty log message ***
File size: 13.0 KB
Line 
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
28Double_t PrimSpect(Double_t *x, Double_t *k)
29{
30 return pow(pow(10, x[0]), k[0]);
31}
32
33Double_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
43Double_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
54Double_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
73Double_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
120void 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// -------------------------------------------------------------------
446void 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
Note: See TracBrowser for help on using the repository browser.