]>
Commit | Line | Data |
---|---|---|
0ca57c2f | 1 | #include "TauolaParticle.h" |
e1938fed | 2 | #include "TauolaLog.h" |
0ca57c2f | 3 | #include <stdexcept> |
4 | ||
5 | namespace Tauolapp | |
6 | { | |
7 | ||
8 | double TauolaParticle::getPolarimetricX(){ | |
9 | return m_pol_x; | |
10 | } | |
11 | ||
12 | double TauolaParticle::getPolarimetricY(){ | |
13 | return m_pol_y; | |
14 | } | |
15 | ||
16 | double TauolaParticle::getPolarimetricZ(){ | |
17 | return m_pol_z; | |
18 | } | |
19 | ||
20 | ||
21 | TauolaParticle * TauolaParticle::clone(){ | |
22 | ||
23 | return createNewParticle(getPdgID(),getStatus(),getMass(), | |
24 | getPx(),getPy(),getPz(),getE()); | |
25 | ||
26 | } | |
27 | ||
28 | // NOTE: Not executed by release examples | |
29 | double TauolaParticle::getAngle(TauolaParticle * other_particle){ | |
30 | ||
31 | //use the dot product | |
32 | double x1 = getPx(); | |
33 | double y1 = getPy(); | |
34 | double z1 = getPz(); | |
35 | double x2 = other_particle->getPx(); | |
36 | double y2 = other_particle->getPy(); | |
37 | double z2 = other_particle->getPz(); | |
38 | ||
39 | return acos( (x1*x2+y1*y2+z1*z2) / sqrt((x1*x1+y1*y1+z1*z1)*(x2*x2+y2*y2+z2*z2)) ); | |
40 | ||
41 | } | |
42 | ||
43 | // NOTE: Not executed by release examples | |
44 | void TauolaParticle::add(TauolaParticle * other_particle){ | |
45 | ||
46 | setPx(getPx() + other_particle->getPx()); | |
47 | setPy(getPy() + other_particle->getPy()); | |
48 | setPz(getPz() + other_particle->getPz()); | |
49 | setE(getE() + other_particle->getE()); | |
50 | setMass( sqrt( getE()*getE()-getPx()*getPx()-getPy()*getPy()-getPz()*getPz() )); | |
51 | } | |
52 | ||
53 | void TauolaParticle::subtract(TauolaParticle * other_particle){ | |
54 | ||
55 | setPx(getPx() - other_particle->getPx()); | |
56 | setPy(getPy() - other_particle->getPy()); | |
57 | setPz(getPz() - other_particle->getPz()); | |
58 | setE(getE() - other_particle->getE()); | |
59 | setMass( sqrt( getE()*getE()-getPx()*getPx()-getPy()*getPy()-getPz()*getPz() )); | |
60 | } | |
61 | ||
62 | int TauolaParticle::getSign(){ | |
63 | if(getPdgID()== Tauola::getDecayingParticle()) | |
64 | return SAME_SIGN; | |
65 | else if(getPdgID()== -1 * Tauola::getDecayingParticle()) | |
66 | return OPPOSITE_SIGN; | |
67 | else | |
68 | return NA_SIGN; | |
69 | } | |
70 | ||
71 | bool TauolaParticle::hasDaughters(){ | |
72 | if(getDaughters().size()==0) | |
73 | return 0; | |
74 | else | |
75 | return 1; | |
76 | } | |
77 | ||
78 | TauolaParticle * TauolaParticle::findLastSelf(){ | |
79 | vector<TauolaParticle*> daughters = getDaughters(); | |
80 | vector<TauolaParticle*>::iterator pcl_itr = daughters.begin(); | |
81 | ||
82 | //get all daughters and look for stable with same pgd id | |
83 | for(;pcl_itr != daughters.end();pcl_itr++){ | |
84 | if((*pcl_itr)->getPdgID()==this->getPdgID()) | |
85 | return (*pcl_itr)->findLastSelf(); | |
86 | } | |
87 | ||
88 | return this; | |
89 | } | |
90 | ||
91 | std::vector<TauolaParticle*> TauolaParticle::findProductionMothers(){ | |
92 | vector<TauolaParticle*> mothers = getMothers(); | |
93 | vector<TauolaParticle*>::iterator pcl_itr = mothers.begin(); | |
94 | ||
95 | //get all mothers and check none have pdg id of this one | |
96 | for(;pcl_itr != mothers.end();pcl_itr++){ | |
97 | if((*pcl_itr)->getPdgID()==this->getPdgID()) | |
98 | return (*pcl_itr)->findProductionMothers(); | |
99 | } | |
100 | return mothers; | |
101 | } | |
102 | ||
103 | void TauolaParticle::decay(){ | |
104 | ||
105 | //Do the decay and set the polarimetric vectors | |
106 | TauolaDecay(getSign(),&m_pol_x, &m_pol_y, &m_pol_z, &m_pol_n); | |
107 | } | |
108 | ||
109 | void TauolaParticle::addDecayToEventRecord(){ | |
110 | ||
111 | //Add to decay list used by f_filhep.c | |
112 | DecayList::addToEnd(this); | |
113 | TauolaWriteDecayToEventRecord(getSign()); | |
114 | ||
115 | double xmom[4]={0}; | |
116 | double *pp=xmom; | |
117 | ||
118 | if (Tauola::ion[2]) | |
119 | for(int i=1;1;i++) | |
120 | { | |
121 | TauolaParticle *x; | |
122 | try{ x=DecayList::getParticle(i); } | |
123 | catch(std::out_of_range d) {break;} | |
124 | if(x->getPdgID()==221){ | |
125 | ||
126 | // Fix 28.04.2011 The pp vector must have boost for eta undone | |
127 | TauolaParticle *x_copy = x->clone(); | |
128 | if(getP(TauolaParticle::Z_AXIS)>0) | |
129 | x_copy->boostAlongZ(-getP(),getE()); | |
130 | else | |
131 | x_copy->boostAlongZ(getP(),getE()); | |
132 | ||
133 | pp[3]=x_copy->getE(); | |
134 | pp[0]=x_copy->getPx(); | |
135 | pp[1]=x_copy->getPy(); | |
136 | pp[2]=x_copy->getPz(); | |
137 | taueta_(pp,&i); | |
138 | } | |
139 | } | |
140 | ||
141 | if (Tauola::ion[1]) | |
142 | for(int i=1;1;i++) | |
143 | { | |
144 | TauolaParticle *x; | |
145 | try{ x=DecayList::getParticle(i); } | |
146 | catch(std::out_of_range d) {break;} | |
147 | if(x->getPdgID()==310){ | |
148 | ||
149 | // Fix 28.04.2011 The pp vector must have boost for k0 undone | |
150 | TauolaParticle *x_copy = x->clone(); | |
151 | if(getP(TauolaParticle::Z_AXIS)>0) | |
152 | x_copy->boostAlongZ(-getP(),getE()); | |
153 | else | |
154 | x_copy->boostAlongZ(getP(),getE()); | |
155 | ||
156 | pp[3]=x_copy->getE(); | |
157 | pp[0]=x_copy->getPx(); | |
158 | pp[1]=x_copy->getPy(); | |
159 | pp[2]=x_copy->getPz(); | |
160 | tauk0s_(pp,&i); | |
161 | } | |
162 | } | |
163 | ||
164 | if (Tauola::ion[0]) | |
165 | for(int i=1;1;i++) | |
166 | { | |
167 | TauolaParticle *x; | |
168 | try{ x=DecayList::getParticle(i); } | |
169 | catch(std::out_of_range d) {break;} | |
170 | if(x->getPdgID()==111){ | |
171 | ||
172 | // Fix 28.04.2011 The pp vector must have boost for pi0 undone | |
173 | TauolaParticle *x_copy = x->clone(); | |
174 | if(getP(TauolaParticle::Z_AXIS)>0) | |
175 | x_copy->boostAlongZ(-getP(),getE()); | |
176 | else | |
177 | x_copy->boostAlongZ(getP(),getE()); | |
178 | ||
179 | pp[3]=x_copy->getE(); | |
180 | pp[0]=x_copy->getPx(); | |
181 | pp[1]=x_copy->getPy(); | |
182 | pp[2]=x_copy->getPz(); | |
183 | taupi0_(pp,&i); | |
184 | } | |
185 | } | |
186 | DecayList::clear(); | |
187 | ||
188 | if(!hasDaughters()) | |
189 | Log::Fatal("TAUOLA failed. No decay was created",5); | |
190 | // checkMomentumConservation(); | |
191 | // decayEndgame(); // vertex shift was wrongly calculated, | |
192 | // used 4-momenta should be in lab frame, | |
193 | // thanks to Sho Iwamoto for debug info. | |
194 | ||
195 | } | |
196 | ||
197 | ||
198 | void TauolaParticle::boostDaughtersFromRestFrame(TauolaParticle * tau_momentum){ | |
199 | ||
200 | if(!hasDaughters()) //if there are no daughters | |
201 | return; | |
202 | ||
203 | vector<TauolaParticle*> daughters = getDaughters(); | |
204 | vector<TauolaParticle*>::iterator pcl_itr = daughters.begin(); | |
205 | ||
206 | //get all daughters then rotate and boost them. | |
207 | for(;pcl_itr != daughters.end();pcl_itr++){ | |
208 | ||
209 | (*pcl_itr)->boostFromRestFrame(tau_momentum); | |
210 | (*pcl_itr)->boostDaughtersFromRestFrame(tau_momentum); | |
211 | } | |
212 | //checkMomentumConservation(); | |
213 | } | |
214 | ||
215 | void TauolaParticle::boostDaughtersToRestFrame(TauolaParticle * tau_momentum){ | |
216 | ||
217 | if(!hasDaughters()) //if there are no daughters | |
218 | return; | |
219 | // NOTE: Not executed by release examples | |
220 | // because !hasDaughters() is always true | |
221 | vector<TauolaParticle*> daughters = getDaughters(); | |
222 | vector<TauolaParticle*>::iterator pcl_itr = daughters.begin(); | |
223 | ||
224 | //get all daughters then rotate and boost them. | |
225 | for(;pcl_itr != daughters.end();pcl_itr++){ | |
226 | ||
227 | (*pcl_itr)->boostToRestFrame(tau_momentum); | |
228 | (*pcl_itr)->boostDaughtersToRestFrame(tau_momentum); | |
229 | } | |
230 | //checkMomentumConservation(); | |
231 | } | |
232 | ||
233 | ||
234 | void TauolaParticle::boostToRestFrame(TauolaParticle * tau_momentum){ | |
235 | ||
236 | double theta = tau_momentum->getRotationAngle(Y_AXIS); | |
237 | tau_momentum->rotate(Y_AXIS,theta); | |
238 | double phi = tau_momentum->getRotationAngle(X_AXIS); | |
239 | tau_momentum->rotate(Y_AXIS,-theta); | |
240 | ||
241 | //Now rotate coordinates to get boost in Z direction. | |
242 | rotate(Y_AXIS,theta); | |
243 | rotate(X_AXIS,phi); | |
244 | boostAlongZ(-1*tau_momentum->getP(),tau_momentum->getE()); | |
245 | rotate(X_AXIS,-phi); | |
246 | rotate(Y_AXIS,-theta); | |
247 | ||
248 | } | |
249 | ||
250 | void TauolaParticle::boostFromRestFrame(TauolaParticle * tau_momentum){ | |
251 | //get the rotation angles | |
252 | //and boost z | |
253 | ||
254 | double theta = tau_momentum->getRotationAngle(Y_AXIS); | |
255 | tau_momentum->rotate(Y_AXIS,theta); | |
256 | double phi = tau_momentum->getRotationAngle(X_AXIS); | |
257 | tau_momentum->rotate(Y_AXIS,-theta); | |
258 | ||
259 | //Now rotate coordinates to get boost in Z direction. | |
260 | rotate(Y_AXIS,theta); | |
261 | rotate(X_AXIS,phi); | |
262 | boostAlongZ(tau_momentum->getP(),tau_momentum->getE()); | |
263 | rotate(X_AXIS,-phi); | |
264 | rotate(Y_AXIS,-theta); | |
265 | } | |
266 | ||
267 | /** Get the angle needed to rotate the 4 momentum vector so that | |
268 | the x (y) component disapears. (and the Z component is > 0) */ | |
269 | double TauolaParticle::getRotationAngle(int axis, int second_axis){ | |
270 | ||
271 | /**if(getP(axis)==0){ | |
272 | if(getPz()>0) | |
273 | return 0; //no rotaion required | |
274 | else | |
275 | return M_PI; | |
276 | }**/ | |
277 | if(getP(second_axis)==0){ | |
278 | if(getP(axis)>0) | |
279 | return -M_PI/2.0; | |
280 | else | |
281 | return M_PI/2.0; | |
282 | } | |
283 | if(getP(second_axis)>0) | |
284 | return -atan(getP(axis)/getP(second_axis)); | |
285 | else | |
286 | return M_PI-atan(getP(axis)/getP(second_axis)); | |
287 | ||
288 | } | |
289 | ||
290 | /** Boost this vector along the Z direction. | |
291 | Assume no momentum components in the X or Y directions. */ | |
292 | void TauolaParticle::boostAlongZ(double boost_pz, double boost_e){ | |
293 | ||
294 | // Boost along the Z axis | |
295 | double m=sqrt(boost_e*boost_e-boost_pz*boost_pz); | |
296 | ||
297 | double p=getPz(); | |
298 | double e=getE(); | |
299 | ||
300 | setPz((boost_e*p + boost_pz*e)/m); | |
301 | setE((boost_pz*p + boost_e*e )/m); | |
302 | } | |
303 | ||
304 | /** Rotation around an axis X or Y */ | |
305 | void TauolaParticle::rotate(int axis,double theta, int second_axis){ | |
306 | ||
307 | double temp_px=getP(axis); | |
308 | double temp_pz=getP(second_axis); | |
309 | setP(axis,cos(theta)*temp_px + sin(theta)*temp_pz); | |
310 | setP(second_axis,-sin(theta)*temp_px + cos(theta)*temp_pz); | |
311 | } | |
312 | ||
313 | void TauolaParticle::rotateDaughters(int axis,double theta, int second_axis){ | |
314 | if(!hasDaughters()) //if there are no daughters | |
315 | return; | |
316 | ||
317 | vector<TauolaParticle*> daughters = getDaughters(); | |
318 | vector<TauolaParticle*>::iterator pcl_itr = daughters.begin(); | |
319 | ||
320 | //get all daughters then rotate and boost them. | |
321 | for(;pcl_itr != daughters.end();pcl_itr++){ | |
322 | ||
323 | (*pcl_itr)->rotate(axis,theta,second_axis); | |
324 | (*pcl_itr)->rotateDaughters(axis,theta,second_axis); | |
325 | } | |
326 | //checkMomentumConservation(); | |
327 | } | |
328 | ||
329 | double TauolaParticle::getMass(){ | |
330 | double e_sq=getE()*getE(); | |
331 | double p_sq=getP()*getP(); | |
332 | ||
333 | if(e_sq>p_sq) | |
334 | return sqrt(e_sq-p_sq); | |
335 | else | |
336 | return -1*sqrt(p_sq-e_sq); //if it's negative | |
337 | } | |
338 | ||
339 | double TauolaParticle::getP(){ | |
340 | return sqrt(getPx()*getPx()+getPy()*getPy()+getPz()*getPz()); | |
341 | } | |
342 | ||
343 | double TauolaParticle::getP(int axis){ | |
344 | if(axis==X_AXIS) | |
345 | return getPx(); | |
346 | ||
347 | if(axis==Y_AXIS) | |
348 | return getPy(); | |
349 | ||
350 | if(axis==Z_AXIS) | |
351 | return getPz(); | |
352 | ||
353 | return 0; | |
354 | } | |
355 | ||
356 | void TauolaParticle::setP(int axis, double p_component){ | |
357 | if(axis==X_AXIS) | |
358 | setPx(p_component); | |
359 | if(axis==Y_AXIS) | |
360 | setPy(p_component); | |
361 | if(axis==Z_AXIS) | |
362 | setPz(p_component); | |
363 | } | |
364 | ||
365 | } // namespace Tauolapp |