flavour
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
widthcalc.h
Go to the documentation of this file.
1 #ifndef WIDTHCALC_H
2 #define WIDTHCALC_H
3 
4 #define _USE_MATH_DEFINES
5 #include <fstream>
6 #include <iostream>
7 #include <set>
8 #include <ginac/ginac.h>
9 #include <cmath>
10 #include <complex>
11 #include "multivector.h"
12 
13 
14 //g++ teste.cpp -o teste -lcln -lginac
15 using namespace std;
16 using namespace GiNaC;
17 
21 class widthcalc{
22 
23 public:
24 
25 widthcalc(): M2(0,2,2,2,2,2,2,2), M22(0,2,2,2,2,2,2), s2("s2"), s3("s3"), mq1("mq1"), mq2("mq2"), mq3("mq3"), mq4("mq4"){
26 
27  integral::max_integration_level=100;
28  integral::relative_integration_error=1e-3;
29 
30  genM2();
31  genM22();
32 }
33 
34 void genM22(){
35  cout<<"Generating M22.dat"<<endl;
36 
37  realsymbol mu("mu"), nu("nu"), alpha("alpha"), beta("beta"), rho("rho"), zeta("zeta");
38  realsymbol q1("q1"), q2("q2"), q3("q3"), q4("q4");
39  realsymbol h1("h1"), h2("h2"), h3("h3"), h4("h4");
40 
41  varidx imu(mu,4,0), inu(nu,4,0), irho(rho,4,0);
42  varidx ialpha(alpha,4,0), ibeta(beta,4,0), izeta(zeta,4,0);
43 
44  varidx jmu(mu,4,1), jnu(nu,4,1), jrho(rho,4,1);
45 
46  ex m2q1=mq1*mq1, m2q2=mq2*mq2, m2q3=mq3*mq3, m2q4=mq4*mq4;
47 
48  ex q2mu=indexed(q2,imu);
49  ex q3mu=indexed(q3,imu);
50  ex q4mu=indexed(q4,imu);
51  ex q1mu=q2mu+q3mu+q4mu;
52 
53  ex s4=m2q1+m2q2+m2q3+m2q4-s2-s3;
54 
55  ex vq1=dirac_slash(q2,4)+dirac_slash(q3,4)+dirac_slash(q4,4)+mq1*dirac_ONE(), vq2=dirac_slash(q2,4)+mq2*dirac_ONE();
56  ex vq3=dirac_slash(q3,4)+mq3*dirac_ONE(), vq4=dirac_slash(q4,4)-mq4*dirac_ONE();
57 
58  scalar_products sp;
59  //sp.add(q2, q3, (s4-m2q2-m2q3)/2);
60  sp.add(q4, q3, (s2-m2q4-m2q3)/2);
61  //sp.add(q2, q4, (s3-m2q2-m2q4)/2);
62 
63  //sp.add(q2, q2, m2q2);
64  sp.add(q3, q3, m2q3);
65  sp.add(q4, q4, m2q4);
66 
67  //sp.add(h1,h1,-1);
68  //sp.add(h2,h2,-1);
69  //sp.add(h3,h3,-1);
70  //sp.add(h4,h4,-1);
71  //sp.add(h1,h1,-1);
72 
73  multivector<ex,3> v(0,2,2,2);
74  v[0][0][0]=dirac_gammaL(); v[0][0][1]=dirac_gammaR();
75  v[0][1][0]=dirac_gammaR(); v[0][1][1]=dirac_gammaL();
76  v[1][0][0]=dirac_gammaL(); v[1][0][1]=dirac_gammaL();
77  v[1][1][0]=dirac_gammaR(); v[1][1][1]=dirac_gammaR();
78 
79  for(uint i=0;i<2;i++)
80  for(uint j=0;j<2;j++){
81  v[0][i][j]=-v[0][i][j]*s2/(mq1+mq2);
82  v[1][i][j]=(dirac_slash(q3,4)+dirac_slash(q4,4))*v[1][i][j];
83  }
84 
85  //vector<ex> prop(2,0);
86  //prop[0]=-s2/(mq1+mq2);
87  //prop[1]=indexed(q3,imu.toggle_variance())+indexed(q4,imu.toggle_variance());
88  /*
89  multivector<ex,2> prop2(0,2,2);
90  for(uint i=0;i<2;i++)
91  for(uint j=0;j<2;j++){
92  prop2[i][j]=prop[i]*prop[j].subs(mu==nu);
93  }
94  */
95 
96  //ofstream f("M22.dat");
97  for(uint i=0;i<2;i++)
98  for(uint j=0;j<2;j++)
99  for(uint k=0;k<2;k++)
100  for(uint l=0;l<2;l++)
101  for(uint m=0;m<2;m++)
102  for(uint n=0;n<2;n++){
103  //cout<<i<<" "<<j<<" "<<k<<" "<<l<<" "<<m<<" "<<n<<endl;
104  //cout<<dirac_trace(vq3*v[i][m][1]*vq4*v[j][n][0].subs(mu==nu))<<endl;
105 
106  ex tmp=dirac_trace(vq3*v[i][m][0]*vq4*v[j][n][1])*int(pow(-1.0,double(k+1+l)));
107  M22[i][j][k][l][m][n]=tmp.simplify_indexed(sp);
108  //cout<<M22[i][j][k][l][m][n]<<endl<<endl;
109  //f<<M22[i][j][k][l][m][n]<<endl;
110  }
111 }
112 
113 void genM2(){
114  cout<<"Generating M2.dat"<<endl;
115 
116  realsymbol mu("mu"), nu("nu"), alpha("alpha"), beta("beta"), rho("rho"), zeta("zeta");
117  realsymbol q1("q1"), q2("q2"), q3("q3"), q4("q4");
118  realsymbol h1("h1"), h2("h2"), h3("h3"), h4("h4");
119 
120  varidx imu(mu,4,0), inu(nu,4,0), irho(rho,4,0);
121  varidx ialpha(alpha,4,0), ibeta(beta,4,0), izeta(zeta,4,0);
122 
123  varidx jmu(mu,4,1), jnu(nu,4,1), jrho(rho,4,1);
124 
125  ex m2q1=mq1*mq1, m2q2=mq2*mq2, m2q3=mq3*mq3, m2q4=mq4*mq4;
126 
127  ex q2mu=indexed(q2,imu);
128  ex q3mu=indexed(q3,imu);
129  ex q4mu=indexed(q4,imu);
130  ex q1mu=q2mu+q3mu+q4mu;
131 
132 
133 
134  ex vq1=dirac_slash(q2,4)+dirac_slash(q3,4)+dirac_slash(q4,4)+mq1*dirac_ONE(), vq2=dirac_slash(q2,4)+mq2*dirac_ONE();
135  ex vq3=dirac_slash(q3,4)+mq3*dirac_ONE(), vq4=dirac_slash(q4,4)-mq4*dirac_ONE();
136 
137  ex s4=m2q1+m2q2+m2q3+m2q4-s2-s3;
138  scalar_products sp;
139  sp.add(q2, q3, (s4-m2q2-m2q3)/2);
140  sp.add(q4, q3, (s2-m2q4-m2q3)/2);
141  sp.add(q2, q4, (s3-m2q2-m2q4)/2);
142 
143  sp.add(q2, q2, m2q2);
144  sp.add(q3, q3, m2q3);
145  sp.add(q4, q4, m2q4);
146 
147  sp.add(h1,h1,-1);
148  sp.add(h2,h2,-1);
149  sp.add(h3,h3,-1);
150  sp.add(h4,h4,-1);
151 
152  sp.add(h2,q2,0);
153  sp.add(h3,q3,0);
154  sp.add(h4,q4,0);
155 
156  multivector<ex,3> v(0,2,2,2);
157  v[0][0][0]=dirac_gammaL(); v[0][0][1]=dirac_gammaR();
158  v[0][1][0]=dirac_gammaR(); v[0][1][1]=dirac_gammaL();
159  v[1][0][0]=dirac_gamma(jmu)*dirac_gammaL(); v[1][0][1]=dirac_gamma(jmu)*dirac_gammaL();
160  v[1][1][0]=dirac_gamma(jmu)*dirac_gammaR(); v[1][1][1]=dirac_gamma(jmu)*dirac_gammaR();
161 
162  multivector<ex,7> traces(0,2,2,2,2,2,2,2);
163  for(uint i=0;i<2;i++)
164  for(uint j=0;j<2;j++)
165  for(uint k=0;k<2;k++)
166  for(uint l=0;l<2;l++)
167  for(uint m=0;m<2;m++)
168  for(uint n=0;n<2;n++){
169  ex vik=v[i][k][0];
170  ex vim=v[i][m][0].subs(mu==nu);
171  ex vjl=v[j][l][1].subs(mu==alpha);
172  ex vjn=v[j][n][1].subs(mu==beta);
173 
174  traces[i][j][k][l][m][n][0]=dirac_trace(vq2*vik*vq1*vjl)*dirac_trace(vq3*vim*vq4*vjn);
175  traces[i][j][k][l][m][n][1]=-dirac_trace(vq2*vik*vq1*vjl*vq3*vim*vq4*vjn);
176  }
177 
178  vector<ex> prop(2,0);
179  prop[0]=1;
180  prop[1]=lorentz_g(imu,inu);
181 
182  multivector<ex,2> prop2(0,2,2);
183  for(uint i=0;i<2;i++)
184  for(uint j=0;j<2;j++){
185  prop2[i][j]=prop[i]*prop[j].subs(lst(mu==alpha,nu==beta));
186  }
187 
188  //ofstream f("M2.dat");
189  for(uint i=0;i<2;i++)
190  for(uint j=0;j<2;j++)
191  for(uint k=0;k<2;k++)
192  for(uint l=0;l<2;l++)
193  for(uint m=0;m<2;m++)
194  for(uint n=0;n<2;n++)
195  for(uint o=0;o<2;o++)
196  {
197  //cout<<i<<" "<<j<<" "<<k<<" "<<l<<" "<<m<<endl;
198  M2[i][j][k][l][m][n][o]=(traces[i][j][k][l][m][n][o]*prop2[i][j]).simplify_indexed(sp);
199  //cout<<M2[i][j][k][l][m][n][o]<<endl<<endl;
200  //f<<M2[i][j][k][l][m][n][o]<<endl;
201  }
202 }
203 
204 
205 ex get_integral(const multivector<ex,4>& a, const vector<ex>& mass, const vector<int>& op, double m1, double m2, double m3, double m4) const{
206 
207  ex q10=(s2+mq1*mq1-mq2*mq2)/(2*sqrt(s2)), lq1l=sqrt(q10*q10-mq1*mq1);
208  ex q30=(s2+mq3*mq3-mq4*mq4)/(2*sqrt(s2)), lq3l=sqrt(q30*q30-mq3*mq3);
209  ex q20=(mq1*mq1+mq2*mq2-s2)/(2*mq1), lq2l=sqrt(q20*q20-mq2*mq2);
210 
211  ex total=0;
212  for(uint i=0;i<a.size();i++) if(!mass[i].is_zero())
213  for(uint j=0;j<a.size();j++)
214  for(uint k=0;k<2;k++)
215  for(uint l=0;l<2;l++)
216  for(uint m=0;m<2;m++)
217  for(uint n=0;n<2;n++)
218  for(uint r=0;r<2;r++)
219  for(uint s=0;s<2;s++){
220  ex coup=a[i][r][k][m]*a[j][s][l][n].conjugate();
221  if(!coup.is_zero()){
222  //cout<<i<<" "<<j<<" "<<k<<" "<<l<<endl;
223  ex integrand=M2[op[i]][op[j]][k][l][m][n][(r+s)%2];
224  integrand=expand(integral(s3, mq1*mq1+mq3*mq3-2*q10*q30-2*lq1l*lq3l, mq1*mq1+mq3*mq3-2*q10*q30+2*lq1l*lq3l, integrand)\
225  .eval_integ()/lq1l/sqrt(s2)*lq2l/mq1/mq1);
226  double mm2=m2, mm3=m3;
227  if(l) {mm2=m3; mm3=m2;}
228  double result=ex_to<numeric>(integral(s2,std::pow(mm3+m4,2),std::pow(m1-mm2,2),integrand.subs(lst(mq1 == m1, mq2 == mm2, mq3 == mm3, mq4 == m4))).evalf()).to_double()/std::pow(M_PI,3)/512;
229  ex partial=result*coup/(pow(mass[i],2)*pow(mass[j],2));
230  //cout<<partial<<endl;
231  total=total+partial;
232  }
233  }
234 
235  return total;
236 }
237 ex get_integral_symb(const multivector<ex,4>& a, const vector<ex>& mass, const vector<int>& op, ex m1) const{
238 
239  ex q10=(s2+m1*m1)/(2*sqrt(s2)), lq1l=(m1*m1-s2)/(2*sqrt(s2));
240  ex q30=sqrt(s2)/2, lq3l=q30;
241  ex q20=(m1*m1-s2)/(2*m1), lq2l=q20;
242 
243  ex total=0;
244  for(uint i=0;i<a.size();i++) if(!mass[i].is_zero())
245  for(uint j=0;j<a.size();j++)
246  for(uint k=0;k<2;k++)
247  for(uint l=0;l<2;l++)
248  for(uint m=0;m<2;m++)
249  for(uint n=0;n<2;n++)
250  for(uint r=0;r<2;r++)
251  for(uint s=0;s<2;s++){
252  ex coup=a[i][r][k][m]*a[j][s][l][n].conjugate();
253  if(!coup.is_zero()){
254  //cout<<i<<" "<<j<<" "<<k<<" "<<l<<endl;
255  ex integrand=M2[op[i]][op[j]][k][l][m][n][(r+s)%2].subs(lst(mq1 == m1, mq2 == 0, mq3 == 0, mq4 == 0));
256  integrand=expand(integral(s3, m1*m1-2*q10*q30-2*lq1l*lq3l, m1*m1-2*q10*q30+2*lq1l*lq3l, integrand)\
257  .eval_integ()/lq1l/sqrt(s2)*lq2l/m1/m1);
258  //integrand=integrand.subs(lst(mq1 == m1, mq2 == 0, mq3 == 0, mq4 == 0));
259  //integrand=integrand.subs(pow(m1,2)/4-s2/2+pow(s2/m1,2)/4==pow((m1-s2/m1)/2,2));
260 
261 
262  double mm2=0, mm3=0, m4=0;
263  //if(l) {mm2=m3; mm3=m2;}
264  ex result=integral(s2,std::pow(mm3+m4,2),pow(m1-mm2,2),integrand).eval_integ()/pow(Pi,3)/512;
265 
266  ex partial=result*coup/(pow(mass[i],2)*pow(mass[j],2));
267  total=total+partial;
268  }
269  }
270  return total;
271 }
272 ex get_integral_meson(const multivector<ex,4>& a, const vector<ex>& mass, const vector<int>& op, ex mm, ex m1, ex m2, ex m3, ex m4) const{
273 
274  ex q10=(s2+mq1*mq1-mq2*mq2)/(2*sqrt(s2)), lq1l=sqrt(q10*q10-mq1*mq1);
275  ex q30=(s2+mq3*mq3-mq4*mq4)/(2*sqrt(s2)), lq3l=sqrt(q30*q30-mq3*mq3);
276  ex total=0;
277  for(uint i=0;i<a.size();i++)
278  for(uint j=0;j<a.size();j++)
279  for(uint k=0;k<2;k++)
280  for(uint l=0;l<2;l++)
281  for(uint m=0;m<2;m++)
282  for(uint n=0;n<2;n++){
283  ex coup=a[i][0][k][m]*a[j][0][l][n].conjugate();
284  if(!coup.is_zero()){
285  //cout<<i<<" "<<j<<" "<<k<<" "<<l<<" "<<m<<" "<<n<<" "<<endl;
286  ex integrand=M22[op[i]][op[j]][k][l][m][n];
287  //cout<<collect_common_factors(expand(a[i][0][k][m]))<<endl;
288  //cout<<collect_common_factors(expand(a[j][0][l][n].conjugate()))<<endl;
289  integrand=expand(integral(s3, mq1*mq1+mq3*mq3-2*q10*q30-2*lq1l*lq3l, mq1*mq1+mq3*mq3-2*q10*q30+2*lq1l*lq3l, integrand)/lq1l/s2);
290  ex result=integrand.subs(lst(sqrt(s2) == mm, s2==mm*mm, mq1 == m1, mq2 == m2, mq3 == m3, mq4 == m4))/Pi/128;
291  ex mi=mass[i];
292  if(mi.is_zero()) mi=mm;
293  ex mj=mass[j];
294  if(mj.is_zero()) mj=mm;
295  ex partial=result*coup/(pow(mi,2)*pow(mj,2));
296  //cout<<i<<"-"<<op[i]<<" "<<j<<"-"<<op[j]<<" "<<a[i]*a[j].conjugate()/(pow(mass[i],2)*pow(mass[j],2))<<endl<<endl;
297 
298  total=total+partial;
299  }
300  }
301 
302  return total;
303 }
304 
305 ex get_integral_meson2(const multivector<ex,4>& a, const vector<ex>& mass, const vector<int>& op, ex mm, ex m1, ex m2, ex m3, ex m4) const{
306 
307  ex q10=(s2+mq1*mq1-mq2*mq2)/(2*sqrt(s2)), lq1l=sqrt(q10*q10-mq1*mq1);
308  ex q30=(s2+mq3*mq3-mq4*mq4)/(2*sqrt(s2)), lq3l=sqrt(q30*q30-mq3*mq3);
309  ex total=0;
310  for(uint i=0;i<a.size();i++)
311  for(uint j=0;j<a.size();j++)
312  for(uint k=0;k<2;k++)
313  for(uint l=0;l<2;l++)
314  for(uint m=0;m<2;m++)
315  for(uint n=0;n<2;n++){
316  ex coup=a[i][0][k][m]*a[j][0][l][n].conjugate();
317  if(!coup.is_zero()){
318  //cout<<i<<" "<<j<<" "<<k<<" "<<l<<" "<<m<<" "<<n<<" "<<endl;
319  ex integrand=M22[op[i]][op[j]][k][l][m][n];
320  //cout<<collect_common_factors(expand(a[i][0][k][m]))<<endl;
321  //cout<<collect_common_factors(expand(a[j][0][l][n].conjugate()))<<endl;
322  integrand=expand(integral(s3, mq1*mq1+mq3*mq3-2*q10*q30-2*lq1l*lq3l, mq1*mq1+mq3*mq3-2*q10*q30+2*lq1l*lq3l, integrand)/lq1l/s2);
323  ex result=integrand.subs(lst(sqrt(s2) == mm, s2==mm*mm, mq1 == m1, mq2 == m2, mq3 == m3, mq4 == m4))/Pi/128;
324  ex mi=mass[i];
325  if(mi.is_zero()) mi=mm;
326  ex mj=mass[j];
327  if(mj.is_zero()) mj=mm;
328  ex partial=result*coup/(pow(mi,2)*pow(mj,2));
329  //cout<<i<<"-"<<op[i]<<" "<<j<<"-"<<op[j]<<" "<<a[i]*a[j].conjugate()/(pow(mass[i],2)*pow(mass[j],2))<<endl<<endl;
330 
331  total=total+partial;
332  }
333  }
334 
335  return total;
336 }
337 /*
338 ex get_integral_meson(const multivector<ex,2>& a, const vector<ex>& mass, const vector<int>& op, double meson_mass, double m3, double m4) const{
339 
340  ex q10=(meson_mass)/(2), lq1l=sqrt(q10*q10-mq1*mq1);
341  ex q30=(s2+mq3*mq3-mq4*mq4)/(2*sqrt(s2)), lq3l=sqrt(q30*q30-mq3*mq3);
342  ex q20=(mq1*mq1+mq2*mq2-s2)/(2*mq1), lq2l=sqrt(q20*q20-mq2*mq2);
343 
344 
345  ex total=0;
346  for(uint i=0;i<a.size();i++)
347  for(uint j=0;j<a.size();j++)
348  for(uint k=0;k<2;k++)
349  for(uint l=0;l<2;l++) if(!(a[i][k]*a[j][l]).is_zero()){
350  //cout<<i<<" "<<j<<" "<<k<<" "<<l<<endl;
351  ex integrand=M2[op[i]/2][op[j]/2][op[i]%2][op[j]%2][(k+l)%2];
352  integrand=expand(integral(s3, mq1*mq1+mq3*mq3-2*q10*q30-2*lq1l*lq3l, mq1*mq1+mq3*mq3-2*q10*q30+2*lq1l*lq3l, integrand)\
353  .eval_integ()/lq1l/sqrt(s2)*lq2l/mq1/mq1);
354  double mm2=m2, mm3=m3;
355  if(l) {mm2=m3; mm3=m2;}
356  double result=ex_to<numeric>(integral(s2,pow(mm3+m4,2),pow(m1-mm2,2),integrand.subs(lst(mq1 == m1, mq2 == mm2, mq3 == mm3, mq4 == m4))).evalf()).to_double()/pow(M_PI,3)/512;
357  ex partial=result*a[i][k]*a[j][l].conjugate()/(pow(mass[i],2)*pow(mass[j],2));
358  //cout<<partial<<endl;
359  total=total+partial;
360  }
361 
362  return total;
363  }
364 */
365 multivector<ex,7> M2;
366 multivector<ex,6> M22;
367 realsymbol s2, s3;
368 realsymbol mq1, mq2, mq3, mq4;
369 
370 };
371 
372 
373 
374 #endif
widthcalc
this class calculates decay widths of one lepton to 3 leptons
Definition: widthcalc.h:21
widthcalc::get_integral
ex get_integral(const multivector< ex, 4 > &a, const vector< ex > &mass, const vector< int > &op, double m1, double m2, double m3, double m4) const
Definition: widthcalc.h:205
std::multivector
A vector of vectors of vectors of... (N times) of class T objects.
Definition: multivector.h:8
std
Definition: multivector.h:4
widthcalc::get_integral_meson
ex get_integral_meson(const multivector< ex, 4 > &a, const vector< ex > &mass, const vector< int > &op, ex mm, ex m1, ex m2, ex m3, ex m4) const
Definition: widthcalc.h:272
widthcalc::get_integral_symb
ex get_integral_symb(const multivector< ex, 4 > &a, const vector< ex > &mass, const vector< int > &op, ex m1) const
Definition: widthcalc.h:237
widthcalc::widthcalc
widthcalc()
Definition: widthcalc.h:25
widthcalc::s3
realsymbol s3
Definition: widthcalc.h:367
widthcalc::get_integral_meson2
ex get_integral_meson2(const multivector< ex, 4 > &a, const vector< ex > &mass, const vector< int > &op, ex mm, ex m1, ex m2, ex m3, ex m4) const
Definition: widthcalc.h:305
widthcalc::M2
multivector< ex, 7 > M2
Definition: widthcalc.h:365
uint
unsigned int uint
Definition: script.cpp:4
widthcalc::genM22
void genM22()
Definition: widthcalc.h:34
widthcalc::M22
multivector< ex, 6 > M22
Definition: widthcalc.h:366
widthcalc::mq4
realsymbol mq4
Definition: widthcalc.h:368
widthcalc::genM2
void genM2()
Definition: widthcalc.h:113