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

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