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

Last change on this file since 11787 was 1449, checked in by tbretz, 23 years ago
*** empty log message ***
File size: 7.5 KB
Line 
1///////////////////////////////////////////////////////////////////////
2//
3// MParticle
4//
5///////////////////////////////////////////////////////////////////////
6
7#include "MParticle.h"
8
9#include <TRandom.h>
10#include <TMatrixD.h>
11#include <TRotation.h>
12
13ClassImp(MParticle);
14
15/**************************************************
16 *
17 * H0 = 50./3.0857e19; // [km / Mpc s] --> [s^-1]
18 *
19 **************************************************/
20
21Double_t MParticle::ZofR(Double_t *x, Double_t *k)
22{
23 /*
24 const Double_t c = 299792458; // [m/s]
25 const Double_t H0 = 50./3.0857e19; // [km / Mpc s] --> [s^-1]
26
27 const Double_t ly = 3600.*24.*365.*c; // [m/ly]
28 const Double_t pc = 1./3.258; // [pc/ly]
29 const Double_t r = x[0] /pc*ly*1e3; // [m]
30
31 const Double_t R = r*H0/c; // [1]
32
33 return (R+1+sqrt(R*2+1))/2 - 1;
34 */
35 const Double_t c = 299792458; // [m/s]
36 const Double_t H0 = 50./3.0857e19; // [km / Mpc s] --> [s^-1]
37
38 const Double_t ly = 3600.*24.*365.*c; // [m/ly]
39 const Double_t pc = 1./3.258; // [pc/ly]
40 const Double_t r = x[0] /pc*ly*1e3; // [m]
41
42 const Double_t R = 1./(1.-r*H0/c/2); // [1]
43
44 return R*R - 1;
45}
46
47Double_t MParticle::RofZ(Double_t *x, Double_t *k)
48{
49 /*
50 Double_t z1 = x[0] + 1;
51
52 const Double_t c = 299792458; // [m/s]
53 const Double_t H0 = 50./3.0857e19; // [km / Mpc s] --> [s^-1]
54
55 const Double_t ly = 3600.*24.*365.*c; // [m/ly]
56 const Double_t pc = 1./3.258; // [pc/ly]
57
58 const Double_t R = c/H0 * 2 * (z1 - sqrt(z1)); // [m]
59
60 return R * pc/ly/1e3; // [kpc]
61 */
62 Double_t z1 = x[0] + 1;
63
64 const Double_t c = 299792458; // [m/s]
65 const Double_t H0 = 50./3.0857e19; // [km / Mpc s] --> [s^-1]
66
67 const Double_t ly = 3600.*24.*365.*c; // [m/ly]
68 const Double_t pc = 1./3.258; // [pc/ly]
69
70 const Double_t R = c/H0 * 2 * (1. - 1./sqrt(z1)); // [m]
71
72 return R * pc/ly/1e3; // [kpc]
73}
74
75#include <fstream.h>
76#include <TH2.h>
77#include "../mhist/MBinning.h"
78#include "../mhist/MH.h"
79
80TH2D *hist2;
81
82Double_t MParticle::Planck(Double_t *x, Double_t *k)
83{
84 static Bool_t isloaded = kFALSE;
85
86 if (!isloaded)
87 {
88 Double_t c = 299792458; // [m/s]
89 Double_t e = 1.602176462e-19; // [C]
90 Double_t h = 1e-9/e*6.62606876e-34; // [GeVs]
91 Double_t hc = h*c; // [GeVm]
92 Double_t konst = 4.*TMath::Pi() * 2. / (hc*hc*hc);
93
94 ifstream fin("mphys/background.txt");
95
96 hist2 = new TH2D;
97
98 MBinning binsz;
99 MBinning binse;
100 binsz.SetEdgesLog(100, 1e-6, 1); // --> 101 Edges / 100 bins
101 binse.SetEdgesLog(100, 7e-15, 3e-8); // --> 101 Edges / 100 bins
102
103 MH::SetBinning(hist2, &binsz, &binse);
104
105 for (int y=0; y<101; y++)
106 {
107 Double_t val;
108 fin >> val;
109 for (int x=0; x<101; x++)
110 {
111 fin >> val;
112
113 val += 9;
114
115 Double_t z = binsz.GetEdges()[x];
116 Double_t f = z+1;
117
118 Double_t newval = pow(10, val)/konst;
119 hist2->SetBinContent(x, y, newval*f*f*f);
120
121 }
122 }
123 isloaded = kTRUE;
124 }
125
126 static TAxis &axez = *hist2->GetXaxis();
127 static TAxis &axee = *hist2->GetYaxis();
128
129 //
130 // y = (y1-y0)/(x1-x0) * (x-x0) + y0
131 //
132 Double_t z = k ? k[0] : 0;
133 Double_t E = x[0];
134
135 Int_t i = axez.FindFixBin(z);
136 Int_t j = axee.FindFixBin(E);
137
138 Double_t z1 = axez.GetBinLowEdge(i+1);
139 Double_t z0 = axez.GetBinLowEdge(i);
140
141 Double_t E1 = axee.GetBinLowEdge(j+1);
142 Double_t E0 = axee.GetBinLowEdge(j);
143
144 Double_t n00 = hist2->GetBinContent(i, j);
145 Double_t n01 = hist2->GetBinContent(i+1, j);
146 Double_t n10 = hist2->GetBinContent(i, j+1);
147 Double_t n11 = hist2->GetBinContent(i+1, j+1);
148
149 Double_t dz = (z-z0)/(z1-z0);
150 Double_t dE = (E-E0)/(E1-E0);
151
152 Double_t n0 = dz*(n01-n00)+n00;
153 Double_t n1 = dz*(n11-n10)+n10;
154
155 Double_t n = dE*(n1-n0)+n0;
156
157 return n;
158 /*
159 //
160 // TANJA2
161 //
162 Double_t E1 = x[0];
163 Double_t E2 = x[0]/8;
164 return (MParticle::Planck0(&E1, k)+MParticle::Planck0(&E2, k)/40e3)*0.7/0.4;
165 */
166 /*
167 //
168 // TANJA
169 //
170 Double_t E1 = x[0];
171 Double_t E2 = x[0]/8;
172 return Planck0(&E1, k)+Planck0(&E2, k)/5e3;
173 */
174}
175
176Double_t MParticle::Planck0(Double_t *x, Double_t *k)
177{
178 //
179 // Planck, per unit volume, per unit energy
180 //
181 // constants (see below) moved out of function
182 //
183 const Double_t E = x[0]; // [GeV]
184 const Double_t z = k ? k[0] : 0;
185
186 const Double_t T = 2.96*(z+1); // [K]
187 const Double_t e = 1.602176462e-19; // [C]
188 const Double_t kB = 1e-9/e*1.3806503e-23; // [GeV/K]
189
190 const Double_t EkT = E/kB/T;
191
192 /*
193 //Double_t c = 299792458; // [m/s]
194 //Double_t h = 1e-9/e*6.62606876e-34; // [GeVs]
195 //Double_t hc = h*c; // [GeVm]
196 Double_t konst = 4.*TMath::Pi() * 2. / (hc*hc*hc);
197 return konst * E*E / (exp(EkT)-1.); // [1 / GeV / m^3 ]
198 */
199
200 return E*E / (exp(EkT)-1.); // [GeV^2]
201}
202
203MParticle::MParticle(ParticleType_t t, const char *name, const char *title)
204 : fPType(t), fZ(0), fR(0), fPhi(0), fTheta(0), fPsi(0), fX(0)
205{
206 //
207 // default constructor
208 //
209}
210
211void MParticle::InitRandom()
212{
213 fPhi = gRandom->Uniform(TMath::Pi()*2);
214 fPsi = gRandom->Uniform(TMath::Pi()*2);
215}
216
217void MParticle::SetNewDirection(Double_t theta, Double_t phi)
218{
219 static TMatrixD B(3, 3);
220
221 const Double_t ct = cos(fTheta);
222 const Double_t st = sin(fTheta);
223 const Double_t cp = cos(fPsi);
224 const Double_t sp = sin(fPsi);
225
226 /*
227 TRotation B( ct*cp, ct*sp, -st,
228 -sp, cp, 0,
229 st*cp, st*sp, ct);
230 */
231
232 // first row
233 B(0, 0) = ct*cp;
234 B(1, 0) = ct*sp;
235 B(2, 0) = -st;
236
237 // second row
238 B(0, 1) = -sp;
239 B(1, 1) = cp;
240 B(2, 1) = 0;
241
242 // third row
243 B(0, 2) = st*cp;
244 B(1, 2) = st*sp;
245 B(2, 2) = ct;
246
247 // ------------------------------
248
249 static TVectorD r(3);
250
251 const Double_t sint = sin(theta);
252
253 /*
254 TVector3 r(sint*cos(phi), sint*sin(phi), cos(theta));
255 */
256
257 r(0) = sint*cos(phi);
258 r(1) = sint*sin(phi);
259 r(2) = cos(theta);
260
261 // ------------------------------
262
263 r *= B;
264
265 fTheta = asin(sqrt(r(0)*r(0)+r(1)*r(1))); // Numerically bad: acos(r(2));
266 fPsi = atan2(r(1), r(0));
267
268 /*
269 if (fTheta*2 > TMath::Pi())
270 fTheta = fabs(fTheta-TMath::Pi());
271 */
272}
273
274Bool_t MParticle::SetNewPosition(Double_t dr)
275{
276 Bool_t rc=kTRUE;
277
278 const Double_t st = sin(fTheta);
279
280 TVector3 d(st*cos(fPsi), st*sin(fPsi), cos(fTheta));
281
282 // ------------------------------
283
284 const Double_t R = RofZ(&fZ);
285 const Double_t dx = R - dr*d.z();
286 if (dx < 0)
287 {
288 dr = R/d.z(); // R>0 --> x(2)>0
289 rc = kFALSE;
290 }
291
292 fX += dr*(1.-d(2));
293 d *= dr;
294
295 // ------------------------------
296
297 TVector3 r(fR*cos(fPhi), fR*sin(fPhi), R);
298
299 r -= d;
300
301 // ------------------------------
302
303 if (fR!=0)
304 fPhi = atan2(r.y(), r.x());
305 fR = sqrt(r.x()*r.x()+r.y()*r.y());
306 fZ = ZofR(&r(2));
307
308 return rc;
309}
310
311Bool_t MParticle::SetNewPosition()
312{
313 Double_t r = gRandom->Exp(GetInteractionLength());
314 return SetNewPosition(r);
315}
Note: See TracBrowser for help on using the repository browser.