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

Last change on this file since 1364 was 1364, checked in by tbretz, 22 years ago
*** empty log message ***
File size: 11.2 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 <TBranch.h>
13#include <TStyle.h>
14#include <TRandom.h>
15#include <TCanvas.h>
16
17#include "mbase/MParContainer.h"
18#include "mphys/MPhoton.h"
19#include "mphys/MElectron.h"
20#include "mphys/MPairProduction.h"
21#include "mhist/MBinning.h"
22#include "mhist/MH.h"
23
24// 2.96
25// 2.87
26// 2.73
27Double_t PrimSpect(Double_t *x, Double_t *k)
28{
29 return pow(pow(10, x[0]), k[0]);
30}
31
32Double_t PhotonSpect(Double_t *x, Double_t *k=NULL)
33{
34 Double_t Ep = pow(10, x[0]);
35 Double_t res = MPhoton::Int2(&Ep, k);
36
37 return res*1e55;
38}
39
40Double_t Sbar_sigmas(Double_t *x, Double_t *k)
41{
42 Double_t sbar = pow(10, x[0]);
43
44 Double_t s = 1./(sbar*4);
45
46 Double_t sigma = MPhoton::Sigma_gg(&s);
47
48 return sigma*sbar*1e28;
49}
50
51Double_t RandomTheta(Double_t Eg, Double_t Ep)
52{
53 Double_t E0 = 511e-6; // [GeV]
54
55 Double_t f = Eg/E0*Ep/E0;
56
57 if (f<1)
58 return 0;
59
60 TF1 func("RndTheta", Sbar_sigmas, 0, log10(f), 0);
61
62 Double_t sbar = pow(10, func.GetRandom());
63 Double_t theta = acos(1.-sbar*2/f);
64
65 return theta;
66}
67
68Double_t GetEnergy(Int_t n, Double_t lo, Double_t hi)
69{
70 static int bin=0;
71 static TRandom rnd(0);
72
73 Double_t w = log10(hi/lo)/n;
74
75 Double_t E = lo*pow(10, rnd.Uniform(w) + w*bin);
76
77 //cout << endl << w << " " << n << " " << w*bin << " " << E << endl;
78
79 if (++bin==n)
80 bin=0;
81
82 return E;
83
84 /*
85 // TF1 src ("PrimarySpectrum", PrimSpect, log10(lo), log10(hi), 1); // 10GeV-1000PeV
86 // src.SetParameter(0, 0);
87 // Double_t E = pow(10, src.GetRandom());
88
89 //
90 // 0: 11111111111111|22222222222222 2 2^0 + 1 1
91 // 1: 11111111|22222222222|33333333 3 2^1 + 1 2
92 // 2: 11111|22222|33333|44444|55555 5 2^2 + 1 7
93 // 3: | | | | | | | | 15
94 //
95
96 static int run=0;
97 static int num=1;
98
99 Double_t w = log10(hi/lo)/((1<<run) + 1);
100
101 Double_t E = lo*pow(10, num*w);
102
103 if (num == (1<<run))
104 {
105 run++;
106 num=1;
107 }
108 else
109 num++;
110
111 return E;
112 */
113}
114
115void DoIt()
116{
117 Double_t R = 500; // [kpc]
118 Double_t startz = MParticle::ZofR(&R);
119
120 cout << "R = " << R << "kpc" << endl;
121 cout << "Z = " << startz << endl;
122
123 static TRandom rand(0);
124 MPairProduction pair;
125
126 Double_t runtime = 7*60*60; // [s]
127
128 Double_t lo = 1e4;
129 Double_t hi = 1e11;
130 Double_t alpha = -2;
131
132 Int_t nbins = 21; //4*log10(hi/lo);
133
134 TH1D hist;
135 TH1D histsrc;
136
137 hist.SetMinimum(lo*lo * pow(lo, alpha)/10);
138 histsrc.SetMinimum(lo*lo * pow(lo, alpha)/10);
139
140 TList listg;
141 TList liste;
142
143 listg.SetOwner();
144 liste.SetOwner();
145
146 gStyle->SetOptStat(10);
147
148 TH2D position;
149 TH2D angle;
150 TH1D histpos;
151 TH1D hista;
152
153 MBinning binspolx;
154 MBinning binspoly1;
155 MBinning binspoly2;
156 binspolx.SetEdges(16, -180, 180);
157 binspoly1.SetEdges(20, 0, 2e-6);
158 binspoly2.SetEdges(20, 0, 1e-5);
159 MH::SetBinning(&angle, &binspolx, &binspoly1);
160 MH::SetBinning(&position, &binspolx, &binspoly2);
161 MH::SetBinning(&hista, &binspoly1);
162 MH::SetBinning(&histpos, &binspoly2);
163
164 hista.SetTitle("Particle Angle");
165 angle.SetTitle("Particle Angle");
166 histpos.SetTitle("Particle Position");
167 position.SetTitle("Particle Position");
168
169 //
170 // Don't change the order!!!
171 //
172 hist.SetFillStyle(0);
173 hist.SetMarkerStyle(kPlus);
174 histsrc.SetFillStyle(0);
175 histsrc.SetMarkerStyle(kMultiply);
176 histsrc.SetMarkerColor(kRed);
177 histsrc.SetLineColor(kRed);
178
179 MBinning bins;
180 bins.SetEdgesLog(nbins, lo, hi);
181
182 MH::SetBinning(&hist, &bins);
183 MH::SetBinning(&histsrc, &bins);
184
185 MH::MakeDefCanvas();
186
187 gPad->SetGrid();
188 gPad->SetLogx();
189 gPad->SetLogy();
190
191 histsrc.SetName("Spectrum");
192 histsrc.SetXTitle("E [GeV]");
193 histsrc.SetYTitle("E^{2}\\dotCounts");
194 histsrc.GetXaxis()->SetLabelOffset(-0.015);
195 histsrc.GetXaxis()->SetTitleOffset(1.1);
196
197 histsrc.Draw("P");
198 hist.Draw("Psame");
199 histsrc.Draw("Csame");
200 hist.Draw("Csame");
201
202 // ------------------------------
203
204 MPhoton *p4file = NULL;
205 MElectron *e4file = NULL;
206
207 TFile file("cascade.root", "RECREATE", "Intergalactic cascade", 9);
208 TTree *T = new TTree("Photons", "Photons from Cascade");
209 TTree *T2 = new TTree("Electrons", "Electrons in the Cascade");
210 TBranch *B =T ->Branch("MPhoton.", "MPhoton", &p4file, 32000);
211 TBranch *B2=T2->Branch("MElectron.", "MElectron", &e4file, 32000);
212
213 // ------------------------------
214
215 TStopwatch clock;
216 clock.Start();
217
218 Int_t timer = 0;
219
220 Int_t n=0;
221 Double_t starttime = TStopwatch::GetRealTime(); // s
222 while (TStopwatch::GetRealTime()<starttime+runtime)
223 for (int i=0; i<nbins; i++)
224 {
225 n++;
226
227 Double_t E = GetEnergy(nbins, lo, hi);
228 Double_t weight = pow(E, alpha);
229 histsrc.Fill(E, E*E * weight);
230
231 MPhoton *gamma=new MPhoton(E, startz);
232 gamma->SetIsPrimary();
233 listg.Add(gamma);
234
235 B->SetAddress(&gamma);
236 T->Fill();
237
238 cout << "--> " << n << ". " << Form("%d: %3.1e", (int)(starttime+runtime-TStopwatch::GetRealTime()), E) << flush;
239
240 while (1)
241 {
242 if (listg.GetSize()!=0)
243 cout << " |P:" << flush;
244
245 TIter NextP(&listg);
246 MPhoton *p = NULL;
247 B->SetAddress(&p);
248 while ((p=(MPhoton*)NextP()))
249 {
250 Double_t Eg = p->GetEnergy();
251
252 if (!p->SetNewPosition())
253 {
254 cout << "!" << flush;
255
256 hist.Fill(Eg, Eg*Eg*weight);
257 position.Fill(p->GetPhi()*kRad2Deg, p->GetR(), weight);
258 angle.Fill(p->GetPsi()*kRad2Deg, p->GetTheta()*kRad2Deg, weight);
259 histpos.Fill(p->GetR(), weight);
260 hista.Fill(p->GetTheta()*kRad2Deg, weight);
261
262 p->SetIsPrimary(kFALSE);
263 T->Fill();
264
265 delete listg.Remove(p);
266 continue;
267 }
268
269 //
270 // Sample phtoton from background and interaction angle
271 //
272 TF1 phot("PhotonSpectrum", PhotonSpect, -log10(Eg)-8, -10.5, 2);
273 phot.SetParameter(0, Eg);
274 phot.SetParameter(1, p->GetZ());
275 if (phot.Integral(-log10(Eg)-8, -10.5, (Double_t*)NULL, 1e-2)==0)
276 {
277 cout << "z" << flush;
278 continue;
279 }
280
281 Double_t Ep = pow(10, phot.GetRandom());
282 Double_t theta = RandomTheta(Eg, Ep);
283 if (theta==0)
284 {
285 cout << "t" << flush;
286 continue;
287 }
288
289 if (!pair.Process(p, Ep, theta, &liste))
290 {
291 cout << "0" << flush;
292 continue;
293 }
294
295 delete listg.Remove(p);
296 cout << "." << flush;
297 }
298
299 if (liste.GetSize()==0 && listg.GetSize()==0)
300 break;
301
302 if (liste.GetSize()!=0)
303 cout << " |E" << flush;
304
305 TIter Next(&liste);
306 MElectron *e = NULL;
307 B2->SetAddress(&e);
308 while ((e=(MElectron*)Next()))
309 {
310 e->SetIsPrimary();
311 T2->Fill();
312 e->SetIsPrimary(kFALSE);
313
314 if (e->GetEnergy()<lo)
315 continue;
316
317 cout << ":" << flush;
318 while(1)
319 {
320 if (!e->SetNewPosition())
321 {
322 cout << "!" << flush;
323 break;
324 }
325
326 // WRONG!
327 Double_t theta = rand.Uniform(TMath::Pi()/2)+TMath::Pi()*3/4;
328 MPhoton *p = e->DoInvCompton(theta);
329
330 T2->Fill();
331
332 if (fabs(p->GetTheta()*3437)<1 && // < 1min
333 p->GetEnergy()>lo)
334 {
335 cout << "." << flush;
336 listg.Add(p);
337 }
338 else
339 {
340 cout << "i" << flush; // ignored
341 delete p;
342 }
343
344 if (fabs(e->GetTheta()*3437)<1 && // < 1min
345 e->GetEnergy()>lo &&
346 e->GetEnergy()>3e4)
347 continue;
348
349 cout << "0" << flush;
350 break;
351 }
352 }
353 liste.Delete();
354 }
355 cout << endl;
356
357 Int_t now = (int)(TStopwatch::GetRealTime()- starttime)/5;
358 if (now!=timer || n<10)
359 {
360 histsrc.SetTitle(Form("E^%.1f, z=%f, T=%d'%d\", N=%d", alpha, startz,
361 (int)runtime/60, (int)runtime%60, hist.GetEntries()));
362 gPad->Modified();
363 gPad->Update();
364 timer = now;
365 }
366
367 }
368 cout << endl;
369
370 clock.Stop();
371 clock.Print();
372
373 file.Write();
374
375 cout << "Created " << n << " gammas from source with E^" << alpha << endl;
376 cout << "Processed " << Form("%.1f", n/(TStopwatch::GetRealTime()-starttime)) << " gammas/sec." << endl;
377
378 // ------------------------------
379
380 gPad->Clear();
381
382 hist.SetTitle(Form("E^%.1f, z=%f, T=%d'%d\", N=%d", alpha, startz, (int)runtime/60, (int)runtime%60, n));
383 gPad->Modified();
384 gPad->Update();
385
386 hist.DrawCopy("P");
387 histsrc.DrawCopy("Psame");
388 hist.DrawCopy("Csame");
389 histsrc.DrawCopy("Csame");
390
391 gPad->SetGrid();
392 gPad->SetLogx();
393 gPad->SetLogy();
394
395 MH::MakeDefCanvas();
396 position.SetXTitle("Pos: \\Phi [\\circ]");
397 position.SetYTitle("Pos: R [kpc]");
398 position.DrawCopy("surf1pol");
399 //gPad->SetLogy();
400
401 MH::MakeDefCanvas();
402 angle.SetXTitle("Angle: \\Psi [\\circ]");
403 angle.SetYTitle("Angle: \\Theta [\\circ]");
404 angle.DrawCopy("surf1pol");
405 //gPad->SetLogy();
406
407 MH::MakeDefCanvas();
408 histpos.SetXTitle("R [kpc]");
409 histpos.SetYTitle("Counts");
410 histpos.DrawCopy();
411 gPad->SetLogx();
412
413 MH::MakeDefCanvas();
414 hista.SetXTitle("\\Theta [\\circ]");
415 hista.SetYTitle("Counts");
416 hista.DrawCopy();
417 gPad->SetLogx();
418}
419
420// -------------------------------------------------------------------
421void phys()
422{
423 DoIt();
424
425 /*
426 Double_t E = 1e10;
427 TF1 phot("PhotonSpectrum", PhotonSpect, -log10(E)-8, -10.5, 2);
428 phot.SetParameter(0, E);
429 phot.SetParameter(1, 5);
430 phot.DrawCopy();
431 return;
432 */
433
434 /*
435 Double_t Eg = 1e5;
436
437 Double_t E0 = 511e-6; // [GeV]
438 Double_t lolim = E0*E0/Eg;
439 Double_t inf = 4e-12; //E0*E0/Eg * sqrt(Eg);
440
441 // 1e5: 5e-13, 4e-12
442 // 1e6: 5e-14, 2e-12
443 // 1e8: 5e-16, 1e-12
444 // 1e10: 5e-18, 1e-12
445
446 // 1e5: 2.6e-12, 4e-12
447 // 1e6: 2.6e-13, 2e-12
448 // 1e8: 2.6e-15, 1e-12
449 // 1e10: 1.6e-17, 1e-12
450
451 TF1 f("int2", MPhoton::Int2, lolim, inf, 2);
452 f.SetParameter(0, Eg);
453 f.SetParameter(1, 0);
454
455 MH::MakeDefCanvas();
456 gPad->SetLogx();
457 f.DrawCopy();
458 */
459}
460
Note: See TracBrowser for help on using the repository browser.