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

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