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