source: trunk/WuerzburgSoft/Thomas/mphys/MElectron.cc@ 11787

Last change on this file since 11787 was 1510, checked in by tbretz, 23 years ago
*** empty log message ***
File size: 15.0 KB
Line 
1/* ======================================================================== *\
2!
3! *
4! * This file is part of MARS, the MAGIC Analysis and Reconstruction
5! * Software. It is distributed to you in the hope that it can be a useful
6! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
7! * It is distributed WITHOUT ANY WARRANTY.
8! *
9! * Permission to use, copy, modify and distribute this software and its
10! * documentation for any purpose is hereby granted without fee,
11! * provided that the above copyright notice appear in all copies and
12! * that both that copyright notice and this permission notice appear
13! * in supporting documentation. It is provided "as is" without express
14! * or implied warranty.
15! *
16!
17!
18! Author(s): Harald Kornmayer 1/2001 (harald@mppmu.mpg.de)
19! Author(s): Thomas Bretz 12/2000 (tbretz@uni-sw.gwdg.de)
20!
21! Copyright: MAGIC Software Development, 2000-2001
22!
23!
24\* ======================================================================== */
25
26//////////////////////////////////////////////////////////////////////////////
27// //
28// //
29//////////////////////////////////////////////////////////////////////////////
30#include "MElectron.h"
31
32#include <math.h> // for aqlphas
33#include <iostream.h>
34
35#include <TF1.h>
36#include <TH1.h>
37#include <TMatrixD.h>
38#include <TVectorD.h>
39
40#include <TRandom.h>
41
42#include <TPad.h>
43#include <TCanvas.h>
44
45#include "MPhoton.h"
46
47ClassImp(MElectron);
48
49Double_t MElectron::Sigma_ge(Double_t *x, Double_t *k)
50{
51 const Double_t E0 = 511e-6; // [GeV]
52 const Double_t c = 299792458; // [m/s]
53 const Double_t e = 1.602176462e-19; // [C]
54 const Double_t h = 1e-9/e*6.62606876e-34; // [GeVs]
55 const Double_t a = 1./137; // [1]
56
57 const Double_t o = x[0];
58
59 const Double_t s1 = TMath::Pi()*2;
60 const Double_t s2 = a*h*c/E0; // [m]
61
62 if (o<1e-4)
63 return s2*s2/s1 * 4/3;
64
65 const Double_t o1 = o+1;
66 const Double_t o21 = o*2+1;
67
68 const Double_t s3 = o1/(o*o*o);
69 const Double_t s4 = o*2*o1/o21;
70 const Double_t s5 = log(o21);
71 const Double_t s6 = s5/o/2;
72 const Double_t s7 = (o*3+1)/(o21*o21);
73
74 return s2*s2/s1*(s3*(s4-s5)+s6-s7); // [m^2]
75}
76
77Double_t MElectron::Li(Double_t *x, Double_t *k)
78{
79 const Double_t t = x[0];
80 return log(1.-t)/t;
81}
82
83Double_t MElectron::DiSum(Double_t *x, Double_t *k)
84{
85 Double_t t = x[0];
86
87 const Double_t eps = fabs(t*1e-2);
88
89 Double_t disum = t;
90 Double_t add = 0;
91
92 Int_t n = 2;
93 Double_t pow = t*t; // t^2
94
95 do
96 {
97 add = pow/n/n;
98
99 pow *= t; // pow = t^n
100 n++;
101
102 disum += add;
103
104 } while (fabs(add)>eps);
105
106 return disum;
107}
108
109Double_t MElectron::Li2(Double_t *x, Double_t *k)
110{
111 //
112 // Dilog, Li2
113 // ----------
114 //
115 // Integral(0, 1) = konst;
116 // Double_t konst = 1./6*TMath::Pi()*TMath::Pi();
117 //
118 // x[0]: z
119 //
120 const Double_t z = x[0];
121
122 if (fabs(z)<1)
123 return DiSum(x);
124
125 // TF1 IntLi("Li", Li, 0, z, 0);
126 static TF1 IntLi("Li", Li, 0, 0, 0);
127 const Double_t integ = IntLi.Integral(0, z, (Double_t*)NULL, 1e-2);
128 return -integ;
129}
130
131Double_t MElectron::Flim(Double_t *x, Double_t *k) // F(omegap)-F(omegam) mit b-->1 (Maple)
132{
133 const Double_t w = x[0];
134
135 const Double_t w4 = w*4;
136 const Double_t wsqr = w*w;
137
138 const Double_t u1 = (w*wsqr*16 + wsqr*40 + w*17 + 2)*log(w4 + 1);
139 const Double_t u2 = -w4*(wsqr*2 + w*9 + 2);
140 const Double_t d = w4*(w4 + 1);
141
142 Double_t s = -w*2*(1+1); // -2*omega*(1+beta)
143 const Double_t li2 = Li2(&s);
144
145 const Double_t res = (u1+u2)/d + li2;
146
147 return res; //<1e-10? 0 : res;
148}
149
150Double_t MElectron::Compton(Double_t *x, Double_t *k)
151{
152 const Double_t E0 = 511e-6; //[GeV]
153
154 Double_t epsilon = x[0];
155 Double_t z = k[1];
156
157 const Double_t E = k[0];
158
159 Double_t flim;
160 if (epsilon<1e-14)
161 {
162 const Double_t d = E/(E0*E0);
163
164 Double_t omega1 = 1e-13*d;
165 Double_t omega2 = 1e-12*d;
166
167 const Double_t f1 = Flim(&omega1);
168 const Double_t f2 = Flim(&omega2);
169
170 const Double_t m = log10(f2/f1);
171 const Double_t t = pow(f2, 13)/pow(f1, 12);
172
173 flim = pow(epsilon, m) * t;
174 }
175 else
176 {
177 Double_t omega = epsilon*E/(E0*E0);
178 flim = Flim(&omega);
179 }
180
181 const Double_t n = MParticle::Planck(&epsilon, &z)/epsilon/epsilon; // [1]
182 return flim*n;
183}
184
185Double_t MElectron::InteractionLength(Double_t *E, Double_t *k)
186{
187 // E = electron energy, ~ TeV(?) 1e12
188 // e = photon energy, ~ meV(?) 1e-3
189 // mc^2 = electron rest mass energy ~.5keV(?) .5e3
190 //
191 // x^-1 = int( n(epsilon)/2beta * ((mc^2)^2/eE)^2 * int ( omega*sigma(omega), omega=o-..o+), epsilon=0..inf)
192 //
193 // o+/- = omage_0 (1 +- beta)
194 //
195 // omega_0 = eE/(mc^2)^2 ~1e12*1e-3/.25e6=4e3
196 //
197 // --> x^-1 = (alpha*hc)^2/4pibetaE^2 * int(n(epsilon)/epsilon^2 *( F(o+)-F(o-)), epsilon=0..inf)
198 //
199 // F(o) = -o/4 + (9/4 + 1/o + o/2) * ln(1+2o) + 1/8(1+2o) - 3/8 + Li2(-2o)
200 //
201 // Li2(x) = int(ln(1-t)/t, t=0..x)
202 //
203 // F(o+)~F(2o) = -o/2 + (9/4 + 1/2o + o) * ln(1+4o) + 1/8(1+4o) - 3/8 + Li2(-4o)
204 // F(o-)~F(0) = 14/8 = 1.75
205
206 const Double_t E0 = 511e-6; // [GeV]
207 const Double_t E02 = E0*E0; // [GeV^2]
208 const Double_t c = 299792458; // [m/s]
209 const Double_t e = 1.602176462e-19; // [C]
210 const Double_t h = 1e-9/e*6.62606876e-34; // [GeVs]
211 const Double_t hc = h*c; // [GeVm]
212 const Double_t alpha = 1./137.; // [1]
213
214 const Double_t z = k ? k[0] : 0;
215
216 /* -------------- old ----------------
217 Double_t from = 1e-15;
218 Double_t to = 1e-11;
219 eps = [default];
220 -----------------------------------
221 */
222 static TF1 func("Compton", Compton, 0, 0, 2); // [0, inf]
223
224 const Double_t from = 1e-17;
225 const Double_t to = 2e-11;
226
227 Double_t val[3] = { E[0], z }; // E[GeV]
228
229 Double_t integ = func.Integral(from, to, val, 1e-2); // [Gev] [0, inf]
230
231 const Double_t aE = alpha/E[0]; // [1/GeV]
232
233 const Double_t beta = 1;
234
235 const Double_t konst = 2.*E02/hc/beta; // [1 / GeV m]
236 const Double_t ret = konst * (aE*aE) * integ; // [1 / m]
237
238 const Double_t ly = 3600.*24.*365.*c; // [m/ly]
239 const Double_t pc = 1./3.258; // [pc/ly]
240
241 return (1./ret)/ly*pc/1000; // [kpc]
242}
243
244Double_t MElectron::GetInteractionLength(Double_t energy, Double_t z)
245{
246 return InteractionLength(&energy, &z);
247}
248
249Double_t MElectron::GetInteractionLength() const
250{
251 return InteractionLength((Double_t*)&fEnergy, (Double_t*)&fZ);
252}
253
254// --------------------------------------------------------------------------
255
256/*inline*/ Double_t MElectron::p_e(Double_t *x, Double_t *k)
257{
258 Double_t e = pow(10, x[0]);
259 return Compton(&e, k);
260 /*
261 Double_t z = k[1];
262
263 const Double_t E = k[0];
264
265 const Double_t E0 = 511e-6; //[GeV]
266 const Double_t E02 = E0*E0;
267
268 Double_t omega = e*E/E02;
269
270 const Double_t n = Planck(&e, &z);
271
272 const Double_t F = Flim(&omega)/omega/omega;
273
274 return n*F*1e26;
275 */
276}
277
278Double_t MElectron::G_q(Double_t *x, Double_t *k)
279{
280 const Double_t q = x[0];
281 const Double_t Gamma = k[0];
282
283 const Double_t Gq = Gamma*q;
284
285 const Double_t s1 = 2.*q*log(q);
286 const Double_t s2 = (1.+2.*q);
287 const Double_t s3 = (Gq*Gq)/(1.+Gq)/2.;
288
289 return s1+(s2+s3)*(1.-q);
290}
291
292
293Double_t MElectron::EnergyLoss(Double_t *x, Double_t *k, Double_t *ep)
294{
295 const Double_t E = x[0];
296 const Double_t z = k ? k[0] : 0;
297
298 const Double_t E0 = 511e-6; //[GeV]
299
300 const Double_t lolim = -log10(E)/7*4-13.5;
301
302 static TF1 fP("p_e", p_e, lolim, -10.6, 2);
303 static TF1 fQ("G", G_q, 0, 1., 1);
304
305 fP.SetNpx(50);
306 fQ.SetNpx(50);
307
308 fP.SetRange(lolim, -10.6);
309 fP.SetParameter(0, E);
310 fP.SetParameter(1, z);
311
312 const Double_t e = pow(10, fP.GetRandom());
313
314 if (ep)
315 *ep = e;
316
317 const Double_t omega = (e/E0)*(E/E0);
318 const Double_t Gamma = 4.*omega;
319
320 fQ.SetParameter(0, Gamma);
321
322 const Double_t q = fQ.GetRandom();
323 const Double_t Gq = Gamma*q;
324
325 const Double_t e1 = Gq*E/(1.+Gq);
326
327 return e1;
328}
329
330Double_t MElectron::GetEnergyLoss(Double_t E, Double_t z, Double_t *ep)
331{
332 return EnergyLoss(&E, &z);
333}
334
335Double_t MElectron::GetEnergyLoss(Double_t *ep) const
336{
337 return EnergyLoss((Double_t*)&fEnergy, (Double_t*)&fZ, ep);
338}
339
340
341Double_t Omega_sigmae(Double_t *x, Double_t *k)
342{
343 Double_t sbar = pow(10,x[0]);
344
345 Double_t omega = (sbar-1)/2;
346
347 Double_t sigma = MElectron::Sigma_ge(&omega);
348
349 return (sbar-1)*sigma*1e28;
350}
351
352Double_t RandomThetaE(Double_t Ee, Double_t Ep)
353{
354 Double_t E0 = 511e-6; // [GeV]
355
356 Double_t f = 2*Ee/E0*Ep/E0;
357
358 static TF1 func("RndThetaE", Omega_sigmae, 0, 0, 0);
359
360 Double_t beta = sqrt(1-E0/Ee*E0/Ee);
361
362 //func.SetRange(0, log10(1+f*(1+beta)));
363 func.SetRange(log10(1+f*(1-beta)), log10(1+f*(1+beta)));
364 func.SetNpx(50);
365
366 Double_t sbar = pow(10, func.GetRandom());
367
368 Double_t bcost = 1 - (sbar-1)/f;
369 return bcost;
370
371 /*
372 Double_t theta = acos(bcost/beta);
373 return theta;
374 */
375}
376
377MPhoton *MElectron::DoInvCompton()
378{
379 const Double_t E0 = 511e-6; //[GeV]
380
381 Double_t epsilon;
382 const Double_t e = GetEnergyLoss(&epsilon);
383
384 // epsilon: photon energy befor interaction, lab
385 // e: photon energy after interaction, lab
386
387 const Double_t gamma = fEnergy/E0;
388 const Double_t gammabeta = sqrt(gamma*gamma-1);
389 const Double_t beta = gammabeta/gamma;
390
391 const Double_t bcost = RandomThetaE(fEnergy, epsilon);
392 const Double_t t = acos(bcost/beta);
393
394 const Double_t f = epsilon/fEnergy;
395 const Double_t r = gamma*(1-bcost);
396
397 Double_t arg = (1 - 1/(gamma*r) - f)/(beta-f);
398
399 if (arg<-1 || arg>1)
400 cout << "<" << (int)(t*180/TMath::Pi()) << "°>" << flush;
401
402 //
403 // Due to numerical uncertanties arg can be something like:
404 // 1+1e-15 which we cannot allow
405 //
406 if (arg<-1)
407 arg = -1;
408 if (arg>1)
409 arg = 1;
410
411 const Double_t theta1s = acos(arg);
412 const Double_t thetas = atan(-sin(t)/(beta*r)/*(1-bcost)/gammabeta*/);
413
414 const Double_t thetastar = thetas-theta1s;
415
416 // const Double_t theta1 = atan(sin(thetastar)/(gamma*cos(thetastar)+gammabeta));
417 const Double_t theta1 = atan(sin(thetastar)/(gamma*(cos(thetastar)+beta)));
418
419 fEnergy -= e;
420
421 const Double_t phi = gRandom->Uniform(TMath::Pi()*2);
422
423 MPhoton &p = *new MPhoton(*this, e);
424 p.SetNewDirection(theta1, phi);
425
426 /*
427 const Double_t beta2 = sqrt(1.-E0/fEnergy*E0/fEnergy);
428 const Double_t theta2 = asin((epsilon*sin(t)-e*sin(theta1))/fEnergy/beta2);
429 */
430
431 const Double_t div = sqrt(fEnergy*fEnergy-E0*E0);
432 const Double_t theta2 = asin((epsilon*sin(t)-e*sin(theta1))/div);
433
434 SetNewDirection(theta2, phi);
435
436 return &p;
437}
438
439void MElectron::DrawInteractionLength(Double_t z)
440{
441 if (!gPad)
442 new TCanvas("IL_Electron", "Mean Interaction Length Electron");
443 else
444 gPad->GetVirtCanvas()->cd(3);
445
446 TF1 f2("length", MElectron::InteractionLength, 1e3, 1e10, 0);
447 f2.SetParameter(0, z);
448
449 gPad->SetLogx();
450 gPad->SetLogy();
451 gPad->SetGrid();
452 f2.SetLineWidth(1);
453
454 TH1 &h=*f2.DrawCopy()->GetHistogram();
455
456 h.SetTitle("Mean Interaction Length (Electron)");
457 h.SetXTitle("E [GeV]");
458 h.SetYTitle("x [kpc]");
459
460 gPad->Modified();
461 gPad->Update();
462}
463
464void MElectron::DrawInteractionLength() const
465{
466 DrawInteractionLength(fZ);
467}
468
469Bool_t MElectron::SetNewPositionB(Double_t B)
470{
471 if (B==0)
472 return SetNewPosition();
473
474 Double_t x = gRandom->Exp(GetInteractionLength());
475
476 // -----------------------------
477
478 const Double_t E0 = 511e-6; //[GeV]
479
480 Double_t B_theta = gRandom->Uniform(TMath::Pi());
481 Double_t B_phi = gRandom->Uniform(TMath::Pi()*2);
482
483 static TMatrixD M(3,3);
484
485 M(0, 0) = sin(B_theta)*cos(B_phi);
486 M(1, 0) = cos(B_theta)*cos(B_phi);
487 M(2, 0) = -sin(B_phi);
488
489 M(0, 1) = sin(B_theta)*sin(B_phi);
490 M(1, 1) = cos(B_theta)*sin(B_phi);
491 M(2, 1) = cos(B_phi);
492
493 M(0, 2) = cos(B_theta);
494 M(1, 2) = -sin(B_theta);
495 M(2, 2) = 0;
496
497 const Double_t beta = sqrt(1.-E0/fEnergy*E0/fEnergy);
498
499 //
500 // rotate vector of velocity into a system in
501 // which the x-axis is identical with the B-Field
502 //
503 static TVectorD v(3);
504 v(0) = beta*sin(fTheta)*cos(fPsi);
505 v(1) = beta*sin(fTheta)*sin(fPsi);
506 v(2) = beta*cos(fTheta);
507 v *= M;
508
509 //
510 // Now rotate the system so, that the velocity vector
511 // lays in the y-z-plain
512 //
513 Double_t chi = atan(-v(2)/v(1));
514
515 // -----------------------------
516
517 static TMatrixD N(3,3);
518 N(0, 0) = 1;
519 N(1, 0) = 0;
520 N(2, 0) = 0;
521
522 N(0, 1) = 0;
523 N(1, 1) = -cos(chi);
524 N(2, 1) = -sin(chi);
525
526 N(0, 2) = 0;
527 N(1, 2) = sin(chi);
528 N(2, 2) = -cos(chi);
529 v *= N;
530
531 Double_t beta_p = v(0); // beta parallel
532 Double_t beta_o = v(1); // beta orthogonal
533 // v(2) = 0
534 // -----------------------------
535 static TVectorD p(3);
536
537 Double_t rho = 0;
538 if (B>0)
539 {
540 const Double_t c = 299792458; // [m/s]
541 const Double_t ly = 3600.*24.*365.*c; // [m/ly]
542 const Double_t pc = 1./3.258; // [pc/ly]
543
544 Double_t r = (fEnergy*1e9)*beta_o/(c*B)*(pc/ly/1e3); // [kpc]
545
546 rho = beta_o/beta*x/r; // [2pi]
547
548 // -----------------------------
549
550 if (fabs(rho*3437)>5*60) // > 5*60min
551 {
552 cout << "r" << flush;
553 return kFALSE;
554 }
555
556 p(0) = beta_p/beta*x;
557 p(1) = GetPType()==kEElectron?r*sin(rho):-r*sin(rho);
558 p(2) = r*(1.-cos(rho));
559 }
560 else
561 {
562 p(0) = beta_p/beta*x;
563 p(1) = beta_o/beta*x;
564 p(2) = 0;
565 cout << "------------- HEY! --------------" << endl;
566 }
567
568 static TVectorD s(3);
569 s(0) = fR*cos(fPhi);
570 s(1) = fR*sin(fPhi);
571 s(2) = RofZ(&fZ);
572
573 TMatrixD N2(TMatrixD::kTransposed, N);
574 TMatrixD M2(TMatrixD::kTransposed, M);
575 p *= N2;
576 p *= M2;
577
578 if (p(2)<x) // happens sometimes in case B==0
579 {
580 cout << "----- HA: " << B << " " << x << " " << p(2) << " " << x-p(2) << endl;
581 p(2)=x;
582 }
583
584 s -= p;
585
586 fR = sqrt(s(0)*s(0)+s(1)*s(1));
587 fPhi = atan2(s(1), s(0));
588 fZ = ZofR(&s(2));
589 fX += x-p(2);
590
591 // -----------------------------
592
593 static TVectorD w(3);
594 w(0) = beta_p/beta;
595 w(1) = beta_o/beta*cos(rho);
596 w(2) = beta_o/beta*sin(rho);
597
598 w *= N2;
599 w *= M2;
600
601 fPsi = atan2(w(1), w(0));
602 fTheta = asin(sqrt(w(0)*w(0)+w(1)*w(1)));
603
604 // -----------------------------
605
606 return fZ<0 ? kFALSE : kTRUE;
607}
Note: See TracBrowser for help on using the repository browser.