PhotosParticle.cxx
1#include <vector>
2#include <math.h>
3#include "PhotosParticle.h"
4#include "Log.h"
5using std::vector;
6
7namespace Photospp
8{
9
11{
12 if(getDaughters().size()==0) return false;
13 else return true;
14}
15
17{
18 vector<PhotosParticle*> daughters = getDaughters();
19 vector<PhotosParticle*>::iterator pcl_itr = daughters.begin();
20
21 //get all daughters and look for stable with same pgd id
22 for(;pcl_itr != daughters.end();pcl_itr++)
23 {
24 if((*pcl_itr)->getPdgID()==this->getPdgID())
25 return (*pcl_itr)->findLastSelf();
26 }
27
28 return this;
29}
30
32{
33 vector<PhotosParticle*> mothers = getMothers();
34 vector<PhotosParticle*>::iterator pcl_itr = mothers.begin();
35
36 //get all mothers and check none have pdg id of this one
37 for(;pcl_itr != mothers.end();pcl_itr++)
38 {
39 if((*pcl_itr)->getPdgID()==this->getPdgID())
40 return (*pcl_itr)->findProductionMothers();
41 }
42 return mothers;
43}
44
45vector<PhotosParticle *> PhotosParticle::getDecayTree()
46{
47 vector<PhotosParticle *> particles;
48 particles.push_back(this);
49 vector<PhotosParticle *> daughters = getDaughters();
50 for(int i=0;i<(int)daughters.size();i++)
51 {
52 // Check if we are the first mother of each daughters
53 // If not - skip this daughter
54 PhotosParticle *p = daughters.at(i);
55 vector<PhotosParticle *> mothers = p->getMothers();
56 if(mothers.size()>1 && mothers.at(0)->getBarcode()!=getBarcode()) continue;
57 vector<PhotosParticle *> tree = p->getDecayTree();
58 particles.insert(particles.end(),tree.begin(),tree.end());
59 }
60 return particles;
61}
62
64{
65 if(!hasDaughters()) //if there are no daughters
66 return;
67
68 // get all daughters, granddaughters, etc. then rotate and boost them
69 vector<PhotosParticle*> list = getAllDecayProducts();
71
72 for(;pcl_itr != list.end();pcl_itr++)
73 {
74 (*pcl_itr)->boostFromRestFrame(tau_momentum);
75 }
76
77 //checkMomentumConservation();
78}
79
81{
82 if(!hasDaughters()) //if there are no daughters
83 return;
84
85 // get all daughters, granddaughters, etc. then rotate and boost them
86 vector<PhotosParticle*> list = getAllDecayProducts();
88
89 for(;pcl_itr != list.end();pcl_itr++)
90 {
91 (*pcl_itr)->boostToRestFrame(tau_momentum);
92 }
93
94 //checkMomentumConservation();
95}
96
97
99{
100 double theta = tau_momentum->getRotationAngle(Y_AXIS);
101 tau_momentum->rotate(Y_AXIS,theta);
102
103 double phi = tau_momentum->getRotationAngle(X_AXIS);
104 tau_momentum->rotate(Y_AXIS,-theta);
105
106 //Now rotate coordinates to get boost in Z direction.
107 rotate(Y_AXIS,theta);
108 rotate(X_AXIS,phi);
109 boostAlongZ(-1*tau_momentum->getP(),tau_momentum->getE());
110 rotate(X_AXIS,-phi);
111 rotate(Y_AXIS,-theta);
112}
113
115{
116 //get the rotation angles
117 //and boost z
118
119 double theta = tau_momentum->getRotationAngle(Y_AXIS);
120 tau_momentum->rotate(Y_AXIS,theta);
121
122 double phi = tau_momentum->getRotationAngle(X_AXIS);
123 tau_momentum->rotate(Y_AXIS,-theta);
124
125 //Now rotate coordinates to get boost in Z direction.
126 rotate(Y_AXIS,theta);
127 rotate(X_AXIS,phi);
128 boostAlongZ(tau_momentum->getP(),tau_momentum->getE());
129 rotate(X_AXIS,-phi);
130 rotate(Y_AXIS,-theta);
131}
132
133/** Get the angle needed to rotate the 4 momentum vector so that
134 the x (y) component disapears. (and the Z component is > 0) */
135double PhotosParticle::getRotationAngle(int axis, int second_axis)
136{
137 /**if(getP(axis)==0){
138 if(getPz()>0)
139 return 0; //no rotaion required
140 else
141 return M_PI;
142 }**/
143 if(getP(second_axis)==0)
144 {
145 if(getP(axis)>0) return -M_PI/2.0;
146 else return M_PI/2.0;
147 }
148 if(getP(second_axis)>0) return -atan(getP(axis)/getP(second_axis));
149 else return M_PI-atan(getP(axis)/getP(second_axis));
150
151}
152
153/** Boost this vector along the Z direction.
154 Assume no momentum components in the X or Y directions. */
155void PhotosParticle::boostAlongZ(double boost_pz, double boost_e)
156{
157 // Boost along the Z axis
158 double m_tau=sqrt(boost_e*boost_e-boost_pz*boost_pz);
159
160 double p=getPz();
161 double e=getE();
162
163 setPz((boost_e*p + boost_pz*e)/m_tau);
164 setE((boost_pz*p + boost_e*e )/m_tau);
165}
166
167/** Rotation around an axis X or Y */
168void PhotosParticle::rotate(int axis,double theta, int second_axis)
169{
170 double temp_px=getP(axis);
171 double temp_pz=getP(second_axis);
172 setP(axis,cos(theta)*temp_px + sin(theta)*temp_pz);
173 setP(second_axis,-sin(theta)*temp_px + cos(theta)*temp_pz);
174}
175
176void PhotosParticle::rotateDaughters(int axis,double theta, int second_axis)
177{
178 if(!hasDaughters()) //if there are no daughters
179 return;
180
181 vector<PhotosParticle*> daughters = getDaughters();
182 vector<PhotosParticle*>::iterator pcl_itr = daughters.begin();
183
184 //get all daughters then rotate and boost them.
185 for(;pcl_itr != daughters.end();pcl_itr++)
186 {
187 (*pcl_itr)->rotate(axis,theta,second_axis);
188 (*pcl_itr)->rotateDaughters(axis,theta,second_axis);
189 }
190 //checkMomentumConservation();
191}
192
194{
195 double e_sq=getE()*getE();
196 double p_sq=getP()*getP();
197
198 if(e_sq>p_sq) return sqrt(e_sq-p_sq);
199 else return -1*sqrt(p_sq-e_sq); //if it's negative
200}
201
203{
204 return sqrt(getPx()*getPx()+getPy()*getPy()+getPz()*getPz());
205}
206
207double PhotosParticle::getP(int axis)
208{
209 if(axis==X_AXIS) return getPx();
210 if(axis==Y_AXIS) return getPy();
211 if(axis==Z_AXIS) return getPz();
212 return 0;
213}
214
215void PhotosParticle::setP(int axis, double p_component)
216{
217 if(axis==X_AXIS) setPx(p_component);
218 if(axis==Y_AXIS) setPy(p_component);
219 if(axis==Z_AXIS) setPz(p_component);
220}
221
222} // namespace Photospp
virtual void setE(double e)=0
void rotate(int axis, double phi, int second_axis=Z_AXIS)
virtual std::vector< PhotosParticle * > getDaughters()=0
virtual std::vector< PhotosParticle * > getMothers()=0
PhotosParticle * findLastSelf()
void boostToRestFrame(PhotosParticle *boost)
virtual int getBarcode()=0
virtual void setPy(double py)=0
void rotateDaughters(int axis, double phi, int second_axis=Z_AXIS)
virtual double getPz()=0
virtual double getE()=0
virtual double getPy()=0
virtual double getPx()=0
void setP(int axis, double p_component)
std::vector< PhotosParticle * > getDecayTree()
virtual double getVirtuality()
std::vector< PhotosParticle * > findProductionMothers()
virtual void setPx(double px)=0
void boostAlongZ(double pz, double e)
virtual void setPz(double pz)=0
double getRotationAngle(int axis, int second_axis=Z_AXIS)
void boostDaughtersFromRestFrame(PhotosParticle *boost)
virtual std::vector< PhotosParticle * > getAllDecayProducts()=0
void boostFromRestFrame(PhotosParticle *boost)
void boostDaughtersToRestFrame(PhotosParticle *boost)
STL class.