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

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