1 | ///////////////////////////////////////////////////////////////////////
|
---|
2 | //
|
---|
3 | // MParticle
|
---|
4 | //
|
---|
5 | ///////////////////////////////////////////////////////////////////////
|
---|
6 |
|
---|
7 | #include "MParticle.h"
|
---|
8 |
|
---|
9 | #include <TRandom.h>
|
---|
10 | #include <TMatrixD.h>
|
---|
11 |
|
---|
12 | ClassImp(MParticle);
|
---|
13 |
|
---|
14 | /***********************
|
---|
15 | *
|
---|
16 | * calc r from z:
|
---|
17 | *
|
---|
18 | * R = 2*c/H0 *(1+z - sqrt(1+z))
|
---|
19 | *
|
---|
20 | * z = -0.5 * (3 + R' +/- sqrt(R'+1)) R' = R*H0/c
|
---|
21 | *
|
---|
22 | * H0 = 50./3.0857e19; // [km / Mpc s] --> [s^-1]
|
---|
23 | *
|
---|
24 | ************************
|
---|
25 | */
|
---|
26 |
|
---|
27 | Double_t ZofR(Double_t *x, Double_t *k=NULL)
|
---|
28 | {
|
---|
29 | const Double_t c = 299792458; // [m/s]
|
---|
30 | const Double_t H0 = 50./3.0857e19; // [km / Mpc s] --> [s^-1]
|
---|
31 |
|
---|
32 | const Double_t ly = 3600.*24.*365.*c; // [m/ly]
|
---|
33 | const Double_t pc = 1./3.258; // [pc/ly]
|
---|
34 | const Double_t r = x[0] /pc*ly*1e3; // [m]
|
---|
35 |
|
---|
36 | const Double_t R = r*H0/c; // [1]
|
---|
37 |
|
---|
38 | return (R+1+sqrt(R*2+1))/2 - 1;
|
---|
39 | }
|
---|
40 |
|
---|
41 | Double_t RofZ(Double_t *x, Double_t *k=NULL)
|
---|
42 | {
|
---|
43 | Double_t z1 = x[0] + 1;
|
---|
44 |
|
---|
45 | const Double_t c = 299792458; // [m/s]
|
---|
46 | const Double_t H0 = 50./3.0857e19; // [km / Mpc s] --> [s^-1]
|
---|
47 |
|
---|
48 | const Double_t ly = 3600.*24.*365.*c; // [m/ly]
|
---|
49 | const Double_t pc = 1./3.258; // [pc/ly]
|
---|
50 |
|
---|
51 | const Double_t R = c/H0 * 2 * (z1 - sqrt(z1)); // [m]
|
---|
52 |
|
---|
53 | return R * pc/ly/1e3; // [kpc]
|
---|
54 | }
|
---|
55 |
|
---|
56 | MParticle::MParticle(ParticleType_t t, const char *name, const char *title)
|
---|
57 | : fPType(t), fZ(0), fR(0), fPhi(0), fTheta(0), fPsi(0)
|
---|
58 | {
|
---|
59 | //
|
---|
60 | // default constructor
|
---|
61 | // creates an a list of histograms for all pixels and both gain channels
|
---|
62 | //
|
---|
63 |
|
---|
64 | //
|
---|
65 | // set the name and title of this object
|
---|
66 | //
|
---|
67 |
|
---|
68 | fName = name ? name : "MParticle";
|
---|
69 | fTitle = title ? title : "Container for a particle";
|
---|
70 | }
|
---|
71 |
|
---|
72 | void MParticle::SetNewDirection(Double_t theta, Double_t phi)
|
---|
73 | {
|
---|
74 | TMatrixD B(3, 3);
|
---|
75 |
|
---|
76 | B(0, 0) = cos(fTheta)*cos(fPsi);
|
---|
77 | B(1, 0) = cos(fTheta)*sin(fPsi);
|
---|
78 | B(2, 0) = -sin(fTheta);
|
---|
79 |
|
---|
80 | B(0, 1) = -sin(fPsi);
|
---|
81 | B(1, 1) = cos(fPsi);
|
---|
82 | B(2, 1) = 0;
|
---|
83 |
|
---|
84 | B(0, 2) = sin(fTheta)*cos(fPsi);
|
---|
85 | B(1, 2) = sin(fTheta)*sin(fPsi);
|
---|
86 | B(2, 2) = cos(fTheta);
|
---|
87 |
|
---|
88 | // ------------------------------
|
---|
89 |
|
---|
90 | TVectorD r(3);
|
---|
91 |
|
---|
92 | r(0) = sin(theta)*cos(phi);
|
---|
93 | r(1) = sin(theta)*sin(phi);
|
---|
94 | r(2) = cos(theta);
|
---|
95 |
|
---|
96 | // ------------------------------
|
---|
97 |
|
---|
98 | r *= B;
|
---|
99 |
|
---|
100 | fTheta = sqrt(r(0)*r(0)+r(1)*r(1)); // Numerically bad: acos(r(2));
|
---|
101 | fPsi = atan2(r(1), r(0));
|
---|
102 |
|
---|
103 | if (fTheta*2 > TMath::Pi())
|
---|
104 | fTheta = fabs(fTheta-TMath::Pi());
|
---|
105 | }
|
---|
106 |
|
---|
107 | #include <iostream.h>
|
---|
108 | #include <iomanip.h>
|
---|
109 |
|
---|
110 | Bool_t MParticle::SetNewPosition(Double_t dr)
|
---|
111 | {
|
---|
112 | Bool_t rc=kTRUE;
|
---|
113 |
|
---|
114 | TVectorD x(3);
|
---|
115 |
|
---|
116 | x(0) = sin(fTheta)*cos(fPsi);
|
---|
117 | x(1) = sin(fTheta)*sin(fPsi);
|
---|
118 | x(2) = cos(fTheta);
|
---|
119 |
|
---|
120 | x *= dr;
|
---|
121 |
|
---|
122 | // ------------------------------
|
---|
123 |
|
---|
124 | const Double_t R = RofZ(&fZ);
|
---|
125 |
|
---|
126 | if (x(2) > R*cos(fTheta))
|
---|
127 | {
|
---|
128 | x *= R/dr;
|
---|
129 | rc = kFALSE;
|
---|
130 | }
|
---|
131 |
|
---|
132 | // ------------------------------
|
---|
133 |
|
---|
134 | TVectorD r(3);
|
---|
135 |
|
---|
136 | r(0) = fR*cos(fPhi);
|
---|
137 | r(1) = fR*sin(fPhi);
|
---|
138 | r(2) = R;
|
---|
139 |
|
---|
140 | // ------------------------------
|
---|
141 |
|
---|
142 | r -= x;
|
---|
143 |
|
---|
144 | fR = sqrt(r(0)*r(0)+r(1)*r(1));
|
---|
145 | fPhi = atan2(r(1), r(0));
|
---|
146 | fZ = ZofR(&r(2));
|
---|
147 |
|
---|
148 | return rc;
|
---|
149 | }
|
---|
150 |
|
---|
151 | Bool_t MParticle::SetNewPosition()
|
---|
152 | {
|
---|
153 | static TRandom rand(0);
|
---|
154 | Double_t r = rand.Exp(GetInteractionLength());
|
---|
155 | return SetNewPosition(r);
|
---|
156 | }
|
---|