source: trunk/WuerzburgSoft/Thomas/mphys/energyloss.C@ 2699

Last change on this file since 2699 was 1363, checked in by tbretz, 23 years ago
*** empty log message ***
File size: 12.7 KB
Line 
1#include <math.h>
2#include <iostream.h>
3#include <fstream.h>
4
5#include <TStopwatch.h>
6#include <TF1.h>
7#include <TH1.h>
8#include <TRandom.h>
9#include <TGraph.h>
10#include <TCanvas.h>
11#include <TStyle.h>
12
13#include "mbase/MParContainer.h"
14#include "mphys/MPhoton.h"
15#include "mphys/MElectron.h"
16#include "mphys/MPairProduction.h"
17#include "mhist/MBinning.h"
18#include "mhist/MH.h"
19
20// 2.96
21// 2.87
22// 2.73
23
24// ================================================================
25Double_t DiSum(Double_t *x, Double_t *k=NULL)
26{
27 Double_t t = x[0];
28
29 Double_t disum = t;
30 Double_t add = 0;
31
32 Double_t eps = fabs(t*1e-1);
33
34 Int_t n = 2;
35 Double_t pow = t*t; // t^2
36
37 do
38 {
39 add = pow/n/n;
40
41 pow *= t; // pow = t^n
42 n++;
43
44 disum += add;
45
46 } while (fabs(add)>eps);
47
48 return disum;
49}
50
51Double_t F(Double_t *x, Double_t *k=NULL)
52{
53 Double_t o = x[0];
54 Double_t s = -2.*o;
55
56// if (o<1e-10)
57// return 2.125; //-3./8.;
58
59 return -o/4. + (9./4. + 1./o + o/2.) * log(1.+2.*o) + 1./8.*(1.+2.*o) + MElectron::Li2(&s); //- 3./8.
60}
61
62void rkck(Double_t y[], Double_t dydx[], int n, Double_t x, Double_t h, Double_t yout[],
63 Double_t yerr[], void (*derivs)(Double_t, Double_t[], Double_t[]))
64{
65 /*
66 * ---------------------------------------------------------
67 * Numerical recipes for C, Chapter 16.1, Runge-Kutta Method
68 * ---------------------------------------------------------
69 */
70 const Double_t a2 = 0.2;
71 const Double_t a3 = 0.3;
72 const Double_t a4 = 0.6;
73 const Double_t a5 = 1.0;
74 const Double_t a6 = 0.875;
75 const Double_t b21 = 0.2;
76 const Double_t b31 = 3.0/40.0;
77 const Double_t b32 = 9.0/40.0;
78 const Double_t b41 = 0.3;
79 const Double_t b42 = -0.9;
80 const Double_t b43 = 1.2;
81 const Double_t b51 = -11.0/54.0;
82 const Double_t b52 = 2.5;
83 const Double_t b53 = -70.0/27.0;
84 const Double_t b54 = 35.0/27.0;
85 const Double_t b61 = 1631.0/55296.0;
86 const Double_t b62 = 175.0/512.0;
87 const Double_t b63 = 575.0/13824.0;
88 const Double_t b64 = 44275.0/110592.0;
89 const Double_t b65 = 253.0/4096.0;
90 const Double_t c1 = 37.0/378.0;
91 const Double_t c3 = 250.0/621.0;
92 const Double_t c4 = 125.0/594.0;
93 const Double_t c6 = 512.0/1771.0;
94 const Double_t dc5 = -277.00/14336.0;
95
96 const Double_t dc1 = c1-2825.0/27648.0;
97 const Double_t dc3 = c3-18575.0/48384.0;
98 const Double_t dc4 = c4-13525.0/55296.0;
99 const Double_t dc6 = c6-0.25;
100
101 Double_t ak2[n];
102 Double_t ak3[n];
103 Double_t ak4[n];
104 Double_t ak5[n];
105 Double_t ak6[n];
106 Double_t ytemp[n];
107
108 for (int i=0; i<n; i++)
109 ytemp[i] = y[i]+b21*h*dydx[i];
110
111 (*derivs)(x+a2*h,ytemp,ak2);
112
113 for (int i=0; i<n; i++)
114 ytemp[i] = y[i]+h*(b31*dydx[i]+b32*ak2[i]);
115
116 (*derivs)(x+a3*h,ytemp,ak3);
117
118 for (int i=0; i<n; i++)
119 ytemp[i] = y[i]+h*(b41*dydx[i]+b42*ak2[i]+b43*ak3[i]);
120
121 (*derivs)(x+a4*h,ytemp,ak4);
122
123 for (int i=0; i<n; i++)
124 ytemp[i] = y[i]+h*(b51*dydx[i]+b52*ak2[i]+b53*ak3[i]+b54*ak4[i]);
125
126 (*derivs)(x+a5*h,ytemp,ak5);
127
128 for (int i=0; i<n; i++)
129 ytemp[i] = y[i]+h*(b61*dydx[i]+b62*ak2[i]+b63*ak3[i]+b64*ak4[i]+b65*ak5[i]);
130
131 (*derivs)(x+a6*h,ytemp,ak6);
132
133 for (int i=0; i<n; i++)
134 yout[i]=y[i]+h*(c1*dydx[i]+c3*ak3[i]+c4*ak4[i]+c6*ak6[i]);
135
136 for (int i=0; i<n; i++)
137 yerr[i]=h*(dc1*dydx[i]+dc3*ak3[i]+dc4*ak4[i]+dc5*ak5[i]+dc6*ak6[i]);
138}
139
140Double_t FMIN(Double_t a, Double_t b)
141{
142 return a<b ? a : b;
143}
144
145Double_t FMAX(Double_t a, Double_t b)
146{
147 return a>b ? a : b;
148}
149
150void SolvEq(Double_t y[], Double_t dydx[], int n, Double_t *x, Double_t htry, Double_t eps,
151 Double_t yscal[], Double_t *hdid, Double_t *hnext,
152 void (*derivs)(Double_t, Double_t[], Double_t[])) // rkqs
153{
154 /*
155 * ---------------------------------------------------------
156 * Numerical recipes for C, Chapter 16.1, Runge-Kutta Method
157 * ---------------------------------------------------------
158 */
159 const Double_t SAFETY = 0.9;
160 const Double_t PGROW = -0.2;
161 const Double_t PSHRNK = -0.25;
162 const Double_t ERRCON = 1.89e-4;
163
164 Double_t yerr[n];
165 Double_t ytemp[n];
166
167 Double_t h = htry;
168 Double_t errmax;
169 while (1)
170 {
171 rkck(y, dydx, n, *x, h, ytemp, yerr, derivs);
172
173 errmax=0.0;
174
175 for (int i=0; i<n; i++)
176 errmax = FMAX(errmax, fabs(yerr[i]/yscal[i]) );
177
178 errmax /= eps;
179
180 if (errmax <= 1.0)
181 break;
182
183 Double_t htemp = SAFETY*h*pow(errmax,PSHRNK);
184
185 h = (h >= 0.0 ? FMAX(htemp,0.1*h) : FMIN(htemp,0.1*h));
186
187 Double_t xnew= (*x) + h;
188
189 if (xnew != *x)
190 continue;
191
192 cout << "stepsize underflow in rkqs" << endl;
193 break;
194 }
195
196 *hnext = errmax>ERRCON ? SAFETY*h*pow(errmax,PGROW) : 5.0*h;
197
198 *x += (*hdid=h);
199
200 for (int i=0; i<n; i++)
201 y[i] = ytemp[i];
202}
203
204Double_t dEdt(Double_t E, Double_t t, Double_t z0=-1)
205{
206 /* ------------------------------------
207 * Lower limit as differential Equation
208 * ------------------------------------
209
210 Double_t T = 2.96; // [K]
211 Double_t sigma = 6.653e-29; // [m^2]
212 Double_t E0 = 511e-6; // [GeV]
213 Double_t e = 1.602176462e-19; // [C]
214 Double_t kB = 1e-9/e*1.3806503e-23; // [GeV/K]
215 Double_t h = 1e-9/e*6.62606876e-34; // [GeVs]
216 Double_t hc = h*c; // [GeVm]
217 Double_t pi = TMath::Pi(); // [1]
218 Double_t khc = pi*kB/hc; // [1 / K m]
219 Double_t a = 8./15 * pi * khc*khc*khc*khc * hc; // [Gev / K^4 / m^3]
220
221 Double_t konst = 4./3. * sigma * a *T*T*T*T *c /E0 /E0;
222
223 return -konst *E*E;
224 */
225 Double_t pc = 1./3.258; // [pc/ly]
226 Double_t c = 299792458; // [m/s]
227 Double_t ly = 3600.*24.*365.*c; // [m/ly]
228 Double_t r = t * c / ly * pc / 1000; // [kpc]
229 Double_t R = MParticle::RofZ(&z0) - r;
230 Double_t z = z0>=0 ? MParticle::ZofR(&R) : 0;
231
232 Double_t e = 1.602176462e-19; // [C]
233 Double_t T = 2.96*(z+1); // [K]
234 Double_t kB = 1e-9/e*1.3806503e-23; // [GeV/K]
235 Double_t kBT = kB*T; // [GeV]
236 Double_t alpha = 1./137; // [|]
237 Double_t h = 1e-9/e*6.62606876e-34; // [GeVs]
238 Double_t E0 = 511e-6; // [GeV]
239
240
241 Double_t k = TMath::Pi() * alpha * kBT;
242
243 Double_t ln = 4.*kBT/E0/E0;
244
245 return -1./3/h* k*k * (log(ln*E)-1.9805);
246}
247
248void dEdt(Double_t t, Double_t E[], Double_t dedt[])
249{
250 dedt[0] = dEdt(E[0], t);
251 //cout << t << "\t" << E[0] << "\t" << dedt[0] << endl;
252}
253
254void DrawDevelopmentHiLim(Double_t E0, Double_t z, Option_t *opt="")
255{
256 Double_t t = 0;
257 Double_t E[1] = { E0 };
258 Double_t yscal[1] = { 1 };
259
260 Double_t dedt[1];
261
262 Double_t eps;// = 1e5;
263 Double_t step = 5e6;
264
265 Double_t hdid;
266 Double_t hnext;
267
268 cout << "Start: " << dedt[0] << endl;
269
270 Int_t n = 15000;
271 Double_t tres[n];
272 Double_t Eres[n];
273 int i;
274 for (i=0; i<n; i++)
275 {
276 tres[i] = t;
277
278 eps = E[0]/1e9; //9;
279
280 dedt[0] = dEdt(E[0], t);
281
282 SolvEq(E, dedt, 1, &t, step, eps, yscal, &hdid, &hnext, dEdt);
283
284 step = hnext;
285
286 Eres[i] = E[0];
287
288 if (i==0) cout << "Did: " << hdid << endl;
289
290 if (t>1.5e14 || E[0]<5e6)
291 break;
292// cout << tres[i] << "\t" << Eres[i] << "\t(" << step << ")" << endl;
293 }
294
295 cout << i << endl;
296
297 TGraph grp(i<n?i:n, tres, Eres);
298 grp.Clone()->Draw(opt);
299
300}
301
302Double_t EnergyLossRateLoLim(Double_t *x, Double_t *k=NULL)
303{
304 Double_t t = x[0];
305 Double_t E = k[0];
306 Double_t t0 = k[1];
307
308 Double_t c = 299792458; // [m/s]
309 Double_t ly = 3600.*24.*365.*c; // [m/ly]
310 Double_t pc = 1./3.258; // [pc/ly]
311
312 Double_t r = t * c / ly * pc / 1000; // [kpc]
313 Double_t R = MParticle::RofZ(&k[2]) - r;
314 Double_t z = k[2]>=0 ? MParticle::ZofR(&R) : 0;
315
316 Double_t T = 2.96*(z+1); // [K]
317 // Double_t alpha = 1./137; // [1]
318 Double_t sigma = 6.653e-29; // [m^2]
319 Double_t E0 = 511e-6; // [GeV]
320 Double_t e = 1.602176462e-19; // [C]
321 Double_t kB = 1e-9/e*1.3806503e-23; // [GeV/K]
322 Double_t h = 1e-9/e*6.62606876e-34; // [GeVs]
323 Double_t hc = h*c; // [GeVm]
324 Double_t pi = TMath::Pi(); // [1]
325 Double_t khc = pi*kB/hc; // [1 / K m]
326 Double_t a = 8./15 * pi * khc*khc*khc*khc * hc; // [Gev / K^4 / m^3]
327 Double_t konst = 4./3 * sigma * a * T*T*T*T * c / (E0* E0); // [1 / GeV s]
328
329 Double_t ret = 1./(konst*(t-t0) + 1./E);
330
331 return ret;
332}
333
334void DrawDevelopmentLoLim(Double_t t0, Double_t E0, Double_t z=-1, Option_t *opt="")
335{
336 // 8.7
337
338 Double_t val[] = { E0, t0, z };
339 Double_t t = 1.5e14;
340 while (EnergyLossRateLoLim(&t, val)<1e4)
341 t -= 0.01e14;
342
343 TF1 *f1=new TF1("LoLim", EnergyLossRateLoLim, t, 1.5e14, 2);
344 f1->SetParameter(0, E0);
345 f1->SetParameter(1, t0);
346
347 f1->Draw(opt);
348 f1->SetBit(kCanDelete);
349}
350
351//
352// (3) Energy loss rate of electrons and 'high energy particle'
353//
354Double_t DrawDevelopment(Double_t E, Double_t z, Option_t *opt="", TH1 *hist=NULL)
355{
356 Double_t ly = 3600.*24.*365.; // [s/ly]
357 Double_t pc = 1./3.258; // [pc/ly]
358
359 TGraph *graph = new TGraph;
360 graph->SetPoint(0, 0, E);
361 graph->SetMaximum(E*3); // *MENU*
362
363 cout << "------ " << endl;
364
365 static TRandom rand;
366
367 Double_t x = 0;
368 for (int i=1; i<10; i++)
369 {
370 Double_t l = rand.Exp(MElectron::InteractionLength(&E, &z));
371
372 if (z>=0)
373 {
374 Double_t r = MParticle::RofZ(&z) - l;
375
376 cout << " " << i << ". R=" << MParticle::RofZ(&z) << " l=" << l << " z=" << MParticle::ZofR(&r) << endl;
377
378 z = r>=0 ? MParticle::ZofR(&r) : 0;
379
380 if (z==0)
381 cout << "z<0" << endl;
382 }
383
384 x += l;
385
386 Double_t t = x/pc*ly*1000;
387 graph->SetPoint(i*2-1, t, E);
388
389 Double_t e1 = MElectron::GetEnergyLoss(E, z<0?0:z);
390
391 E -= e1;
392
393 if (hist)
394 hist->Fill(e1);
395
396 cout << " Ep=" << e1 << flush;
397
398 graph->SetPoint(i*2, t, E);
399 }
400
401 graph->SetMinimum(E/3); // *MENU*
402 graph->Draw(opt);
403 graph->SetBit(kCanDelete);
404
405 //if (E<31500)
406 cout << "t=" << x*ly/pc*1000 << "\tE=" << E << "\tz=" << z << endl;
407
408 return E;
409}
410
411void EnergyLossRate()
412{
413 if (gPad)
414 gPad->Clear();
415
416 Double_t E = 1.5e9; // [GeV]
417 Double_t z = 0.03;
418
419 MBinning bins;
420 bins.SetEdgesLog(18, 0.1, 1e9);
421
422 TH1D hist;
423 hist.SetName("Phot");
424 hist.SetTitle("Photons from inverse Compton");
425 MH::SetBinning(&hist, &bins);
426
427 cout << "Working..." << flush;
428
429 for (int i=0; i<50; i++)
430 {
431 cout << i << "." << flush;
432 DrawDevelopment(E, z, i?"L":"AL", &hist);
433 }
434
435 //DrawDevelopmentLoLim(2e14, 1.64e2, "Lsame"); // seen
436 DrawDevelopmentLoLim(1.78e14, 280, z, "Lsame"); // seen
437 DrawDevelopmentHiLim(E, z, "L");
438 gPad->SetLogy();
439
440 new TCanvas("Photons", "Photons created in inverse Compton");
441 hist.DrawCopy();
442 gPad->SetLogx();
443
444 cout << "...done." << endl;
445}
446
447void DrawRZ()
448{
449 new TCanvas("RZ", "r and z");
450
451 TF1 f1("ZofR", MParticle::ZofR, 0, 7.1e6, 0);
452 TF1 f2("RofZ", MParticle::RofZ, 0, 5, 0);
453
454 gPad->Divide(2,2);
455
456 gPad->GetVirtCanvas()->cd(1);
457 TH1 *h = f1.DrawCopy()->GetHistogram();
458
459 h->SetTitle("z(r)");
460 h->SetXTitle("r [kpc]");
461 h->SetYTitle("z");
462
463 gPad->Modified();
464 gPad->Update();
465
466 gPad->GetVirtCanvas()->cd(2);
467 h = f2.DrawCopy()->GetHistogram();
468
469 h->SetTitle("r(z)");
470 h->SetXTitle("z");
471 h->SetYTitle("r [kpc]");
472
473 gPad->Modified();
474 gPad->Update();
475}
476
477// -------------------------------------------------------------------
478
479Double_t func(Double_t *x, Double_t *k)
480{
481 return MPhoton::Int2(x, k)*1e68;
482}
483
484void energyloss()
485{
486/* Double_t E0 = 511e-6; // [GeV]
487
488 Double_t Eg = 1e4; //3.6e4;
489 Double_t z = 5;
490
491 Double_t val[2] = { Eg, z };
492
493 Double_t lolim = E0*E0/Eg;
494 Double_t inf = Eg<1e6 ? 3e-11*z : 3e-12*z;
495
496 cout << Eg << " " << z << " " << lolim << " " << inf << endl;
497
498 TF1 f("int2", func, lolim, inf, 2);
499
500 Double_t int2 = f.Integral(lolim, inf, val); //[GeV^3 m^2]
501
502 f.SetParameter(0, Eg);
503 f.SetParameter(1, z);
504
505 cout << int2 << endl;
506
507 new TCanvas("ILPhoton", "Mean Interaction Length Photon");
508
509 gPad->SetLogx();
510// gPad->SetLogy();
511 gPad->SetGrid();
512
513 f.SetLineWidth(1);
514 f.DrawCopy();
515
516 gPad->Modified();
517 gPad->Update();
518
519 // return;
520 */ MPhoton p;
521 p.DrawInteractionLength();
522
523 return;
524// EnergyLossRate();
525
526 DrawRZ();
527
528 return;
529
530}
Note: See TracBrowser for help on using the repository browser.