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

Last change on this file since 1373 was 1373, checked in by tbretz, 22 years ago
*** empty log message ***
File size: 12.4 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 = 100; //MParticle::RofZ(&startz); // [kpc]
121 Double_t startz = MParticle::ZofR(&R);
122
123 const char *filename = "data/cascade_delme.root";
124
125 const Double_t B = 0;
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 = 60; // [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 Double_t Esum=0;
249 Double_t Emis=0;
250
251 while (1)
252 {
253 if (listg.GetSize()!=0)
254 cout << " |P:" << flush;
255
256 TIter NextP(&listg);
257 MPhoton *p = NULL;
258 B1->SetAddress(&p);
259 while ((p=(MPhoton*)NextP()))
260 {
261 Double_t Eg = p->GetEnergy();
262
263 if (!p->SetNewPosition())
264 {
265 cout << "!" << flush;
266
267 hist.Fill(Eg, Eg*Eg*weight);
268 position.Fill(p->GetPhi()*kRad2Deg, p->GetR(), weight);
269 angle.Fill(p->GetPsi()*kRad2Deg, p->GetTheta()*kRad2Deg, weight);
270 histpos.Fill(p->GetR(), weight);
271 hista.Fill(p->GetTheta()*kRad2Deg, weight);
272
273 p->SetIsPrimary(kFALSE);
274 T1->Fill();
275
276 Esum += Eg;
277
278 delete listg.Remove(p);
279 continue;
280 }
281
282 //
283 // Sample phtoton from background and interaction angle
284 //
285 TF1 phot("PhotonSpectrum", PhotonSpect, -log10(Eg)-8, -10.5, 2);
286 phot.SetParameter(0, Eg);
287 phot.SetParameter(1, p->GetZ());
288 /*
289 if (phot.Integral(-log10(Eg)-8, -10.5, (Double_t*)NULL, 1e-2)==0)
290 {
291 cout << "z" << flush;
292 continue;
293 }
294 */
295 Double_t Ep = pow(10, phot.GetRandom());
296 if (Ep==pow(10,0))
297 {
298 cout << "z" << flush;
299 continue;
300 }
301 Double_t theta = RandomTheta(Eg, Ep);
302 if (theta==0)
303 {
304 cout << "t" << flush;
305 continue;
306 }
307
308 if (!pair.Process(p, Ep, theta, &liste))
309 {
310 cout << "0" << flush;
311 continue;
312 }
313
314 delete listg.Remove(p);
315 cout << "." << flush;
316 }
317
318 if (liste.GetSize()==0 && listg.GetSize()==0)
319 break;
320
321 if (liste.GetSize()!=0)
322 cout << " |E" << flush;
323
324 TIter Next(&liste);
325 MElectron *e = NULL;
326 B2->SetAddress(&e);
327 while ((e=(MElectron*)Next()))
328 {
329 e->SetIsPrimary();
330 T2->Fill();
331 e->SetIsPrimary(kFALSE);
332
333 if (e->GetEnergy()<lo)
334 continue;
335
336 cout << ":" << flush;
337 while(1)
338 {
339 if (!e->SetNewPositionB(B))
340 {
341 cout << "!" << flush;
342 Esum += e->GetEnergy();
343 Emis += e->GetEnergy();
344 break;
345 }
346
347 // WRONG!
348 Double_t theta = 0; //rand.Uniform(TMath::Pi()/2)+TMath::Pi()*3/4;
349 MPhoton *p = e->DoInvCompton(theta);
350
351 T2->Fill();
352
353 if (fabs(p->GetTheta()*3437)<60 && // < 60min
354 p->GetEnergy()>lo)
355 {
356 cout << "." << flush;
357 listg.Add(p);
358 }
359 else
360 {
361 if (p->GetEnergy()<=lo)
362 cout << "e" << flush;
363 else
364 cout << "t" << flush; // ignored
365 delete p;
366 Esum += p->GetEnergy();
367 Emis += p->GetEnergy();
368 }
369
370 if (fabs(e->GetTheta()*3437)>60) // < 60min
371 {
372 cout << "T" << flush;
373 Esum += e->GetEnergy();
374 Emis += e->GetEnergy();
375 break;
376 }
377
378 if (e->GetEnergy()<lo)//*5)
379 {
380 cout << "E" << flush;
381 Esum += e->GetEnergy();
382 Emis += e->GetEnergy();
383 break;
384 }
385 }
386 }
387 liste.Delete();
388 }
389 cout << " ----> " << Form("%3.1f %3.1e / %3.1f %3.1e", Emis/E, Emis, Esum/E, Esum) << endl;
390
391 Int_t now = (int)(TStopwatch::GetRealTime()- starttime)/5;
392 if (now!=timer || n<10)
393 {
394 histsrc.SetTitle(Form("E^%.1f, z=%f, T=%d'%d\", N=%d", alpha, startz,
395 (int)runtime/60, (int)runtime%60, hist.GetEntries()));
396 gPad->Modified();
397 gPad->Update();
398 timer = now;
399 }
400
401 }
402 cout << endl;
403
404 clock.Stop();
405 clock.Print();
406
407 file.Write();
408
409 cout << "Created " << n << " gammas from source with E^" << alpha << endl;
410 cout << "Processing time: " << Form("%.1f", (TStopwatch::GetRealTime()-starttime)/n) << " sec/gamma." << endl;
411
412 // ------------------------------
413
414 gPad->Clear();
415
416 hist.SetTitle(Form("E^%.1f, z=%f, T=%d'%d\", N=%d", alpha, startz, (int)runtime/60, (int)runtime%60, n));
417 gPad->Modified();
418 gPad->Update();
419
420 hist.DrawCopy("P");
421 histsrc.DrawCopy("Psame");
422 hist.DrawCopy("Csame");
423 histsrc.DrawCopy("Csame");
424
425 gPad->SetGrid();
426 gPad->SetLogx();
427 gPad->SetLogy();
428
429 MH::MakeDefCanvas();
430 position.SetXTitle("Pos: \\Phi [\\circ]");
431 position.SetYTitle("Pos: R [kpc]");
432 position.DrawCopy("surf1pol");
433 //gPad->SetLogy();
434
435 MH::MakeDefCanvas();
436 angle.SetXTitle("Angle: \\Psi [\\circ]");
437 angle.SetYTitle("Angle: \\Theta [\\circ]");
438 angle.DrawCopy("surf1pol");
439 //gPad->SetLogy();
440
441 MH::MakeDefCanvas();
442 histpos.SetXTitle("R [kpc]");
443 histpos.SetYTitle("Counts");
444 histpos.DrawCopy();
445 gPad->SetLogx();
446
447 MH::MakeDefCanvas();
448 hista.SetXTitle("\\Theta [\\circ]");
449 hista.SetYTitle("Counts");
450 hista.DrawCopy();
451 gPad->SetLogx();
452}
453
454// -------------------------------------------------------------------
455void phys()
456{
457 DoIt();
458
459 /*
460 Double_t E = 1e10;
461 TF1 phot("PhotonSpectrum", PhotonSpect, -log10(E)-8, -10.5, 2);
462 phot.SetParameter(0, E);
463 phot.SetParameter(1, 5);
464 phot.DrawCopy();
465 return;
466 */
467
468 /*
469 Double_t Eg = 1e5;
470
471 Double_t E0 = 511e-6; // [GeV]
472 Double_t lolim = E0*E0/Eg;
473 Double_t inf = 4e-12; //E0*E0/Eg * sqrt(Eg);
474
475 // 1e5: 5e-13, 4e-12
476 // 1e6: 5e-14, 2e-12
477 // 1e8: 5e-16, 1e-12
478 // 1e10: 5e-18, 1e-12
479
480 // 1e5: 2.6e-12, 4e-12
481 // 1e6: 2.6e-13, 2e-12
482 // 1e8: 2.6e-15, 1e-12
483 // 1e10: 1.6e-17, 1e-12
484
485 TF1 f("int2", MPhoton::Int2, lolim, inf, 2);
486 f.SetParameter(0, Eg);
487 f.SetParameter(1, 0);
488
489 MH::MakeDefCanvas();
490 gPad->SetLogx();
491 f.DrawCopy();
492 */
493}
494
Note: See TracBrowser for help on using the repository browser.