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

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