source: trunk/WuerzburgSoft/Thomas/mphys/MCascade.cc@ 4939

Last change on this file since 4939 was 1507, checked in by tbretz, 23 years ago
*** empty log message ***
File size: 10.8 KB
Line 
1#include "MCascade.h"
2
3#include <math.h> // fabs, for alpha
4#include <iostream.h>
5
6#include <TF1.h>
7#include <TH2.h>
8#include <TList.h>
9#include <TFile.h>
10#include <TTree.h>
11#include <TTimer.h>
12#include <TStyle.h>
13#include <TBranch.h>
14#include <TCanvas.h>
15#include <TRandom3.h>
16#include <TStopwatch.h>
17
18#include "MPhoton.h"
19#include "MElectron.h"
20
21#include "MH.h"
22#include "MBinning.h"
23
24ClassImp(MCascade);
25
26Double_t PrimSpect(Double_t *x, Double_t *k)
27{
28 return pow(pow(10, x[0]), k[0]);
29}
30
31Double_t PhotonSpect(Double_t *x, Double_t *k=NULL)
32{
33 Double_t Ep = pow(10, x[0]);
34
35 Double_t res = MPhoton::Int2(&Ep, k);
36 return res*1e55; //65/k[0];
37
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 RandomThetaG(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 static TF1 func("RndThetaG", Sbar_sigmas, 0, 0, 0);
62
63 func.SetRange(0, log10(f));
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 MCascade::GetEnergy()
73{
74 static int bin=0;
75 Double_t w = log10(fEHi/fELo)/fNumBins;
76
77 Double_t E = fELo*pow(10, gRandom->Uniform(w) + w*(fNumBins-bin-1));
78
79 if (++bin==fNumBins)
80 bin=0;
81
82 return E;
83}
84
85TFile *MCascade::OpenFile(TString fFilename)
86{
87 TFile *fFile = new TFile(fFilename, "CREATE", "Intergalactic cascade", 9);
88
89 if (fFile->IsZombie())
90 {
91 delete fFile;
92 return NULL;
93 }
94
95 Write("Setup");
96
97 cout << "Trees... " << flush;
98
99 TTree *T1 = new TTree ("Photons", "Photons from Cascade");
100 TTree *T2 = new TTree ("Electrons", "Electrons in the Cascade");
101
102 cout << "Branches... " << flush;
103
104 MPhoton dummyp;
105 void *ptr = &dummyp;
106 fBranchGammas = T1->Branch("MPhoton.", "MPhoton", &ptr);
107
108 MElectron dummye;
109 ptr = &dummye;
110 fBranchElectrons = T2->Branch("MElectron.", "MElectron", &ptr);
111
112 return fFile;
113}
114
115void MCascade::CloseFile(TFile *fFile)
116{
117 fFile->Write();
118 cout << "Wrote: " << fFile->GetName() << endl;
119 delete fFile;
120}
121
122void MCascade::ProcessElectron(MElectron &e, TList &fListGammas)
123{
124 Double_t Ee = e.GetEnergy();
125
126 cout << ":" << flush;
127
128 int test = fNumMaxInvCompton<0 ? -1 : 0;
129 int n=0;
130 while ((test<0 ? true : (test++<fNumMaxInvCompton)) &&
131 (e.GetEnergy() > fRatioInvCompton*Ee))
132 {
133 n++;
134
135 if (!e.SetNewPositionB(e.GetZ()>fBubbleZ ? fB : 0))
136 {
137 cout << "!" << flush;
138 return;
139 }
140
141 MPhoton *p = e.DoInvCompton();
142
143 fBranchElectrons->GetTree()->Fill();
144
145 //cout << "." << flush;
146 fListGammas.Add(p);
147
148 if (fabs(e.GetTheta()*3437)>60) // < 60min
149 {
150 cout << "T" << flush;
151 return;
152 }
153
154 if (e.GetEnergy()<Ee*1e-3) // <2e3
155 {
156 cout << "E" << flush;
157 return;
158 }
159
160 if (e.GetEnergy()<1e2)
161 {
162 cout << "x" << flush;
163 return;
164 }
165 }
166 cout << n << flush;
167}
168
169Bool_t MCascade::ProcessGamma(MPhoton &p, Double_t weight, TList &fListElectrons)
170{
171 Double_t Eg = p.GetEnergy();
172
173 Double_t E0 = 511e-6;
174 Double_t z = p.GetZ()+1;
175 Double_t lolim = E0*E0/Eg;
176 Double_t inf = (Eg<1e6 ? 3e-11*z : 3e-12*z);
177 if (Eg<5e4)
178 inf = 3e-11*z*pow(10, 9.4-log10(Eg)*2);
179
180 TF1 phot("PhotonSpectrum", PhotonSpect, 0, 0, 2);
181 phot.SetRange(log10(lolim), log10(inf));
182 phot.SetNpx(50);
183 phot.SetParameter(0, Eg);
184 while (1)
185 {
186 if (!p.SetNewPosition())
187 return kTRUE;
188
189 //
190 // Sample phtoton from background and interaction angle
191 //
192 phot.SetParameter(1, p.GetZ());
193 Double_t pe = phot.GetRandom();
194 if (pe==0)
195 {
196 cout << "z" << flush;
197 continue;
198 }
199
200 Double_t Ep = pow(10, pe);
201 Double_t theta = RandomThetaG(Eg, Ep);
202 if (theta==0)
203 {
204 cout << "t" << flush;
205 continue;
206 }
207
208 if (!fPair.Process(&p, Ep, theta, &fListElectrons))
209 {
210 // should never happen
211 cout << "0" << flush;
212 continue;
213 }
214
215 return kFALSE;
216 }
217}
218
219Double_t MCascade::ProcessGammas(TList &fListGammas, TList &fListElectrons, Double_t weight)
220{
221 Double_t Esum = 0;
222
223 MPhoton *p = NULL;
224 fBranchGammas->SetAddress(&p);
225
226 cout << ":" << fListGammas.GetSize() << ":" << flush;
227
228 TIter NextP(&fListGammas);
229 while ((p=(MPhoton*)NextP()))
230 {
231 if (!ProcessGamma(*p, weight, fListElectrons))
232 {
233 delete fListGammas.Remove(p);
234 cout << "." << flush;
235 continue;
236 }
237
238 fBranchGammas->GetTree()->Fill();
239 Esum += p->GetEnergy();
240 //cout << "!" << flush;
241 }
242
243 return Esum;
244}
245
246Double_t MCascade::ProcessElectrons(TList &fListElectrons, TList &fListGammas)
247{
248 Double_t E = 0;
249
250 cout << ":" << fListElectrons.GetSize() << flush;
251
252 TIter Next(&fListElectrons);
253 MElectron *e = NULL;
254 fBranchElectrons->SetAddress(&e);
255 while ((e=(MElectron*)Next()))
256 {
257 e->SetIsPrimary(kTRUE);
258 fBranchElectrons->GetTree()->Fill();
259 e->SetIsPrimary(kFALSE);
260
261 ProcessElectron(*e, fListGammas);
262
263 E += e->GetEnergy();
264 }
265 fListElectrons.Delete();
266
267 return E;
268}
269
270void MCascade::ProcessPrimaryGamma(Double_t E, Double_t weight)
271{
272 TList fListGammas;
273 TList fListElectrons;
274
275 fListGammas.SetOwner();
276 fListElectrons.SetOwner();
277
278 MPhoton *gamma=new MPhoton(E, fSrcZ);
279 gamma->SetSrcR(fSrcR);
280 gamma->InitRandom();
281 fListGammas.Add(gamma);
282
283 gamma->SetIsPrimary(kTRUE);
284 fBranchGammas->SetAddress(&gamma);
285 fBranchGammas->GetTree()->Fill();
286 gamma->SetIsPrimary(kFALSE);
287
288 Double_t Esum=0; // sum of all energies
289 Double_t Emis=0; // sum of the energies thrown away
290 while (1)
291 {
292 if (fListGammas.GetSize()==0)
293 break;
294
295 cout << " |P" << flush;
296
297 Esum += ProcessGammas(fListGammas, fListElectrons, weight);
298
299 if (!fIsBatch)
300 fListGammas.ForEach(MPhoton, Fill)(fHist, fDisplayIndex, weight);
301 fListGammas.Delete();
302
303 if (fListElectrons.GetSize()==0)
304 break;
305
306 cout << " |E" << flush;
307
308 Emis += ProcessElectrons(fListElectrons, fListGammas);
309 }
310 Esum += Emis;
311
312 cout << " ----> " << Form("%3.1f %3.1e / %3.1f %3.1e", Emis/E, Emis, Esum/E, Esum) << endl;
313}
314
315MCascade::MCascade()
316{
317 if (gRandom)
318 delete gRandom;
319
320 TRandom r(0);
321 gRandom = new TRandom3(r.GetSeed());
322
323 fHist.SetName("Spectrum");
324 fHist.SetXTitle("E [GeV]");
325 fHist.SetYTitle(Form("E^{%.1f} Counts", fDisplayIndex));
326 fHist.GetXaxis()->SetLabelOffset(-0.015);
327 fHist.GetXaxis()->SetTitleOffset(1.1);
328 fHist.SetFillStyle(0);
329 fHist.SetMarkerStyle(kPlus);
330}
331
332MCascade::~MCascade()
333{
334 delete gRandom;
335 gRandom = 0;
336}
337
338void MCascade::SetSourceZ(Double_t z)
339{
340 fSrcZ = z;
341 fSrcR = MParticle::RofZ(&z);
342}
343
344void MCascade::SetSourceRZ(Double_t r) // [kpc]
345{
346 fSrcZ = MParticle::ZofR(&r);
347 fSrcR = r;
348}
349
350void MCascade::SetEnergyBins(Int_t n, Double_t lo, Double_t hi)
351{
352 fNumBins = n; // number of bins produced in energy spectrum
353
354 fELo = lo; // lower limit of displayed spectrum
355 fEHi = hi; // upper limit of spectrum (cutoff)
356}
357
358void MCascade::SetBradius(Double_t r) // [Mpc]
359{
360 Double_t bubbler = fSrcR-1e3*r;
361 fBubbleZ = MParticle::ZofR(&bubbler);
362}
363
364void MCascade::Run(TString filename, Bool_t draw)
365{
366 fIsBatch = gROOT->IsBatch() ? kFALSE : draw;
367
368 // ------------------------------
369
370 cout << "Output File '" << filename << "'... " << flush;
371
372 TFile *file=OpenFile(filename);
373 if (!file)
374 return;
375
376 // ------------------------------
377
378 cout << endl;
379
380 cout << "R = " << fSrcR << "kpc" << endl;
381 cout << "Z = " << fSrcZ << endl;
382
383 cout << "Setting up: Histograms... " << flush;
384
385 fHist.Reset();
386 TH1D histsrc;
387
388 MBinning bins;
389 bins.SetEdgesLog(fNumBins, fELo, fEHi);
390
391 MH::SetBinning(&fHist, &bins);
392 MH::SetBinning(&histsrc, &bins);
393
394 TCanvas *c=NULL;
395
396 if (!fIsBatch)
397 {
398 fHist.SetMinimum(pow(fELo, fSpectralIndex+fDisplayIndex)/100);
399 histsrc.SetMinimum(pow(fELo, fSpectralIndex+fDisplayIndex)/100);
400
401 gStyle->SetOptStat(10);
402
403 //
404 // Don't change the order!!!
405 //
406 histsrc.SetFillStyle(0);
407 histsrc.SetMarkerStyle(kMultiply);
408 histsrc.SetMarkerColor(kRed);
409 histsrc.SetLineColor(kRed);
410
411 c=MH::MakeDefCanvas("Cascade", "Cascade");
412
413 c->SetGrid();
414 c->SetLogx();
415 c->SetLogy();
416
417 fHist.Draw("P");
418 histsrc.Draw("Psame");
419 histsrc.Draw("same");
420 fHist.Draw("same");
421 }
422
423 // ------------------------------
424
425 cout << "Timers... " << flush;
426
427 TTimer timer("gSystem->ProcessEvents();", 333, kFALSE);
428 if (!fIsBatch)
429 timer.TurnOn();
430
431 TStopwatch clock;
432 clock.Start();
433
434 cout << "Done. " << endl;
435
436 Int_t n=0;
437 Double_t starttime = TStopwatch::GetRealTime(); // s
438 while (TStopwatch::GetRealTime()<starttime+fRuntime)
439 for (int i=0; i<fNumBins; i++)
440 {
441 n++;
442
443 Double_t E = GetEnergy();
444 Double_t weight = pow(E, fSpectralIndex);
445
446 if (!fIsBatch)
447 histsrc.Fill(E, pow(E, fDisplayIndex) * weight);
448
449 cout << "--> " << n << ". " << Form("%d: %3.1e", (int)(starttime+fRuntime-TStopwatch::GetRealTime()), E) << flush;
450
451 ProcessPrimaryGamma(E, weight);
452
453 if (fIsBatch)
454 continue;
455
456 fHist.SetTitle(Form("E^{%.1f} z=%f T=%d'%d\" N=%d", fSpectralIndex, fSrcZ,
457 (int)fRuntime/60, (int)fRuntime%60, n));
458
459 timer.Stop();
460 c->Update();
461 timer.Start(250);
462 }
463
464 cout << endl;
465
466 clock.Stop();
467 clock.Print();
468
469 timer.Stop();
470
471 CloseFile(file);
472
473 cout << "Created " << n << " gammas (" << n/fNumBins << "rows) from source with E^" << fSpectralIndex << endl;
474 cout << "Processing time: " << Form("%.1f", (TStopwatch::GetRealTime()-starttime)/n) << " sec/gamma (";
475 cout << Form("%.1f", (TStopwatch::GetRealTime()-starttime)/n*fNumBins/60) << " min/row)" << endl;
476
477 // ------------------------------
478
479 if (!fIsBatch)
480 {
481 c->Clear();
482
483 fHist.SetTitle(Form("E^{%.1f} z=%f T=%d'%d\" N=%d", fSpectralIndex, fSrcZ, (int)fRuntime/60, (int)fRuntime%60, n));
484
485 TH1 &h1 = *fHist.DrawCopy("P");
486 TH1 &h2 = *histsrc.DrawCopy("Psame");
487 h2.Draw("Csame");
488 h1.Draw("Csame");
489 }
490}
491
Note: See TracBrowser for help on using the repository browser.