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

Last change on this file since 1360 was 1358, checked in by tbretz, 22 years ago
*** empty log message ***
File size: 9.6 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 <iostream.h>
33
34#include <TF1.h>
35#include <TH1.h>
36#include <TPad.h>
37#include <TCanvas.h>
38#include <TRandom.h>
39
40#include "MPhoton.h"
41
42ClassImp(MElectron);
43
44Double_t MElectron::Planck(Double_t *x, Double_t *k=NULL)
45{
46 //
47 // Planck, per unit volume, per unit energy
48 //
49 // constants moved out of function
50 //
51 Double_t E = x[0]; // [GeV]
52 Double_t z = k ? k[0] : 0;
53
54 Double_t T = 2.96*(z+1); // [K]
55 Double_t e = 1.602176462e-19; // [C]
56 Double_t kB = 1e-9/e*1.3806503e-23; // [GeV/K]
57
58 Double_t EkT = E/kB/T;
59
60 /*
61 //Double_t c = 299792458; // [m/s]
62 //Double_t h = 1e-9/e*6.62606876e-34; // [GeVs]
63 //Double_t hc = h*c; // [GeVm]
64 Double_t konst = 4.*TMath::Pi() * 2. / (hc*hc*hc);
65 return konst * E*E / (exp(EkT)-1.); // [1 / GeV / m^3 ]
66 */
67
68 return E*E / (exp(EkT)-1.); // [GeV^2]
69}
70
71Double_t MElectron::Li(Double_t *x, Double_t *k)
72{
73 Double_t t = x[0];
74
75 return log(1.-t)/t;
76}
77
78Double_t MElectron::Li2(Double_t *x, Double_t *k=NULL)
79{
80 //
81 // Dilog, Li2
82 //
83 Double_t z = x[0];
84
85 TF1 IntLi("Li", Li, 0, z, 0);
86 Double_t integ = IntLi.Integral(0, z);
87
88 /*
89 if (fabs(z)<1)
90 {
91 Double_t disum = DiSum(&z);
92 cout << "Disum (" << z << ") " << disum << "=" << -integ << "\t" << disum-integ << endl;
93 return disum;
94 }
95 */
96
97 /*
98 Integral(0, 1) = konst;
99 Double_t konst = 1./6*TMath::Pi()*TMath::Pi();
100 */
101
102 return -integ;
103}
104
105Double_t MElectron::Flim(Double_t *x, Double_t *k=NULL) // F(omegap)-F(omegam) mit b-->1 (Maple)
106{
107 Double_t w = x[0];
108
109 Double_t w4 = w*4;
110 Double_t wsqr = w*w;
111
112 Double_t u1 = (w*wsqr*16 + wsqr*40 + w*17 + 2)*log(w4 + 1);
113 Double_t u2 = -w4*(wsqr*2 + w*9 + 2);
114 Double_t d = w4*(w4 + 1);
115
116 Double_t s = -w*2*(1+1); // -2*omega*(1+beta)
117 Double_t li2 = Li2(&s);
118
119 Double_t res = (u1+u2)/d + li2;
120
121 return res; //<1e-10? 0 : res;
122}
123
124Double_t MElectron::Compton(Double_t *x, Double_t *k)
125{
126 Double_t E0 = 511e-6; //[GeV]
127 Double_t E02 = E0*E0;
128
129 Double_t epsilon = x[0];
130 Double_t E = k[0];
131 // Double_t beta = k[1];
132 Double_t z = k[2];
133
134 Double_t omega = epsilon*E/E02;
135
136 Double_t n = Planck(&epsilon, &z)/epsilon/epsilon; // [1]
137 return Flim(&omega)*n;
138}
139
140
141Double_t MElectron::InteractionLength(Double_t *E, Double_t *k=NULL)
142{
143 // E = electron energy, ~ TeV(?) 1e12
144 // e = photon energy, ~ meV(?) 1e-3
145 // mc^2 = electron rest mass energy ~.5keV(?) .5e3
146 //
147 // x^-1 = int( n(epsilon)/2beta * ((mc^2)^2/eE)^2 * int ( omega*sigma(omega), omega=o-..o+), epsilon=0..inf)
148 //
149 // o+/- = omage_0 (1 +- beta)
150 //
151 // omega_0 = eE/(mc^2)^2 ~1e12*1e-3/.25e6=4e3
152 //
153 // --> x^-1 = (alpha*hc)^2/4pibetaE^2 * int(n(epsilon)/epsilon^2 *( F(o+)-F(o-)), epsilon=0..inf)
154 //
155 // F(o) = -o/4 + (9/4 + 1/o + o/2) * ln(1+2o) + 1/8(1+2o) - 3/8 + Li2(-2o)
156 //
157 // Li2(x) = int(ln(1-t)/t, t=0..x)
158 //
159 // F(o+)~F(2o) = -o/2 + (9/4 + 1/2o + o) * ln(1+4o) + 1/8(1+4o) - 3/8 + Li2(-4o)
160 // F(o-)~F(0) = 14/8 = 1.75
161
162 Double_t E0 = 511e-6; // [GeV]
163 Double_t E02 = E0*E0; // [GeV^2]
164 Double_t c = 299792458; // [m/s]
165 Double_t e = 1.602176462e-19; // [C]
166 Double_t h = 1e-9/e*6.62606876e-34; // [GeVs]
167 Double_t hc = h*c; // [GeVm]
168 Double_t alpha = 1./137.; // [1]
169
170 Double_t z = k ? k[0] : 0;
171
172 // Double_t beta = sqrt(1-E0/E*E0/E);
173 Double_t beta = 1; //0.999999999999999999999999999;
174
175 Double_t val[3] = { E[0], beta, z }; // E[GeV]
176
177 Double_t from = 1e-17;
178 Double_t to = 1e-11;
179
180 /* -------------- old ----------------
181 Double_t from = 1e-15;
182 Double_t to = 1e-11;
183 eps = [default];
184 -----------------------------------
185 */
186 TF1 func("Compton", Compton, from, to, 3); // [0, inf]
187
188 Double_t integ = func.Integral(from, to, val, 1e-15); // [Gev] [0, inf]
189
190 Double_t aE = alpha/E[0]; // [1/GeV]
191
192 Double_t konst = 2.*E02/hc/beta; // [1 / GeV m]
193 Double_t ret = konst * (aE*aE) * integ; // [1 / m]
194
195 Double_t ly = 3600.*24.*365.*c; // [m/ly]
196 Double_t pc = 1./3.258; // [pc/ly]
197
198 return (1./ret)/ly*pc/1000; // [kpc]
199}
200
201Double_t MElectron::GetInteractionLength(Double_t energy, Double_t z)
202{
203 return InteractionLength(&energy, &z);
204}
205
206Double_t MElectron::GetInteractionLength() const
207{
208 return InteractionLength((Double_t*)&fEnergy, (Double_t*)&fZ);
209}
210
211// --------------------------------------------------------------------------
212
213Double_t MElectron::p_e(Double_t *x, Double_t *k)
214{
215 Double_t e = pow(10, x[0]);
216 Double_t E = k[0];
217 Double_t z = k[1];
218
219 Double_t E0 = 511e-6; //[GeV]
220 Double_t E02 = E0*E0;
221
222 Double_t omega = e*E/E02;
223
224 Double_t n = Planck(&e, &z);
225
226 Double_t F = Flim(&omega)/omega/omega;
227
228 return n*F*1e26;
229}
230
231Double_t MElectron::G_q(Double_t *x, Double_t *k)
232{
233 Double_t q = x[0];
234 Double_t Gamma = k[0];
235
236 Double_t Gq = Gamma*q;
237
238 Double_t s1 = 2.*q*log(q);
239 Double_t s2 = (1.+2.*q);
240 Double_t s3 = (Gq*Gq)/(1.+Gq)/2.;
241
242 return s1+(s2+s3)*(1.-q);
243}
244
245
246Double_t MElectron::EnergyLoss(Double_t *x, Double_t *k, Double_t *ep)
247{
248 Double_t E = x[0];
249 Double_t z = k ? k[0] : 0;
250
251 Double_t E0 = 511e-6; //[GeV]
252
253 Double_t lolim = -log10(E)/7*4-13;
254
255 TF1 fP("p_e", p_e, lolim, -10.8, 2);
256 TF1 fQ("G", G_q, 0, 1., 1);
257
258 fP.SetParameter(0, E);
259 fP.SetParameter(1, z);
260
261 Double_t e = pow(10, fP.GetRandom());
262
263 if (ep)
264 *ep = e;
265
266 Double_t omega = e*E/E0/E0;
267 Double_t Gamma = 4.*omega;
268
269 // --old-- fQ.SetRange(1e-6, 1./(1.+ 2.*Gamma));
270 fQ.SetParameter(0, Gamma);
271
272 Double_t q = fQ.GetRandom();
273 Double_t Gq = Gamma*q;
274
275 Double_t e1 = Gq*E/(1.+Gq);
276
277 return e1;
278}
279
280Double_t MElectron::GetEnergyLoss(Double_t E, Double_t z, Double_t *ep)
281{
282 return EnergyLoss(&E, &z);
283}
284
285Double_t MElectron::GetEnergyLoss(Double_t *ep) const
286{
287 return EnergyLoss((Double_t*)&fEnergy, (Double_t*)&fZ, ep);
288}
289
290MPhoton *MElectron::DoInvCompton(Double_t theta)
291{
292 static TRandom rand(0);
293
294 Double_t E0 = 511e-6; //[GeV]
295
296 Double_t epsilon;
297 Double_t e = GetEnergyLoss(&epsilon);
298
299 // er: photon energy before interaction, rest frame
300 // e: photon energy after interaction, lab
301
302 Double_t gamma = fEnergy/E0;
303 Double_t beta = sqrt(1.-1./(gamma*gamma));
304 //Double_t gammabeta = sqrt(gamma*gamma-1);
305
306 Double_t f = fEnergy/e;
307
308 Double_t t;
309 Double_t arg;
310 do
311 {
312 t = rand.Uniform(TMath::Pi()*2);
313 Double_t er = gamma*epsilon*(1.-beta*cos(t)); // photon energy rest frame
314 arg = (f - E0/er - 1)/(f*beta+1);
315 cout << "~" << flush;
316
317 } while (arg<-1 || arg>1);
318
319 Double_t theta1s = acos(arg);
320 Double_t thetas = atan(sin(t)/(gamma*(cos(t)-beta)));
321
322 Double_t thetastar = thetas-theta1s;
323
324 Double_t theta1 = atan(sin(thetastar)/(gamma*(cos(thetastar)+beta)));
325
326 /*
327 cout << "(" << theta1/TMath::Pi()*180 << ",";
328 cout << theta1s/TMath::Pi()*180<< ",";
329 cout << arg << ")" << flush;
330 */
331
332 fEnergy -= e;
333
334 MPhoton &p = *new MPhoton(e, fZ);
335 p = *this;
336 p.SetNewDirection(theta1, rand.Uniform(TMath::Pi()*2));
337
338 // MISSING: Electron angle
339 //
340 // E1 = fEnergy (after!)
341 //
342 // sin(t) = (epsilon sin(theta) - e sin(atan(tg)))/sqrt(E1*E1-E0*E0)
343
344 return &p;
345}
346
347void MElectron::DrawInteractionLength(Double_t z)
348{
349 if (!gPad)
350 new TCanvas("IL_Electron", "Mean Interaction Length Electron");
351 else
352 gPad->GetVirtCanvas()->cd(3);
353
354 TF1 f2("length", MElectron::InteractionLength, 1e3, 1e10, 0);
355 f2.SetParameter(0, z);
356
357 gPad->SetLogx();
358 gPad->SetLogy();
359 gPad->SetGrid();
360 f2.SetLineWidth(1);
361
362 TH1 &h=*f2.DrawCopy()->GetHistogram();
363
364 h.SetTitle("Mean Interaction Length (Electron)");
365 h.SetXTitle("E [GeV]");
366 h.SetYTitle("x [kpc]");
367
368 gPad->Modified();
369 gPad->Update();
370}
371
372void MElectron::DrawInteractionLength() const
373{
374 DrawInteractionLength(fZ);
375}
Note: See TracBrowser for help on using the repository browser.