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