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

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