00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021 #include "EvtGenBase/EvtPatches.hh"
00022 #include <math.h>
00023 #include <stdlib.h>
00024 #include "EvtGenBase/EvtParticle.hh"
00025 #include "EvtGenBase/EvtKine.hh"
00026
00027 #include "EvtGenBase/EvtVector4C.hh"
00028 #include "EvtGenBase/EvtVector4R.hh"
00029 #include "EvtGenBase/EvtTensor4C.hh"
00030 #include "EvtGenBase/EvtHelSys.hh"
00031 #include "EvtGenModels/EvtDIY.hh"
00032 #include <string>
00033 using namespace std;
00034
00035
00036
00037
00038 class AngularSam{
00039
00040 public:
00041 AngularSam(double alpha, EvtVector4R p1){
00042 _a=alpha;
00043 pd1=p1;
00044 }
00045
00046 double amps();
00047
00048 private:
00049 EvtVector4R pd1;
00050 double _a;
00051 };
00052
00053 double AngularSam::amps(){
00054 double costheta=pd1.get(3)/pd1.d3mag();
00055 return 1+_a*pow(costheta,2.);
00056 }
00057
00058
00059
00060 class rhopi{
00061 public:
00062 rhopi(EvtVector4R pd1, EvtVector4R pd2,EvtVector4R pd3){
00063 _pd[0]=pd1;
00064 _pd[1]=pd2;
00065 _pd[2]=pd3;
00066 }
00067 double F00(double s);
00068 double F10(double s);
00069 double amps1(double s,int i, int j);
00070 double amps( );
00071
00072 private:
00073 EvtVector4R _pd[3];
00074 };
00075
00076 double rhopi::F00(double s){
00077 double mpi=0.1395;
00078 return sqrt(s-4*mpi*mpi)/sqrt(s);
00079 }
00080
00081 double rhopi::F10(double s){
00082 double mpi=0.1395,mpsi=3.096916;
00083 double tep=sqrt((mpsi*mpsi-pow(mpi+sqrt(s),2.))*(mpsi*mpsi-pow(mpi-sqrt(s),2.)));
00084 return tep;
00085 }
00086
00087 double rhopi::amps1(double s, int i, int j){
00088 double mrho=0.771,wrho=0.1492,dpro;
00089 EvtComplex img(0.0,1.0);
00090 dpro=pow(abs(s-mrho*mrho+img*sqrt(s)*wrho),2.);
00091 EvtVector4R prho;
00092 prho=_pd[i]+_pd[j];
00093 EvtHelSys angles(prho,_pd[i]),labAngles;
00094 double theta,phi,ct1,st1,phi1,st,ct,temp;
00095 theta=angles.getHelAng(1);
00096 phi =angles.getHelAng(2);
00097 ct1 =labAngles.Angles(prho,1);
00098 phi1=labAngles.Angles(prho,2);
00099 st=sin(theta);ct=cos(theta);
00100 temp=pow(F00(s),2.)*pow(F10(s),2.)*pow(st,2.)/dpro;
00101 return temp;
00102 }
00103
00104 double rhopi::amps(){
00105 double temp,s12,s13,s23;
00106 s12=(_pd[0]+_pd[1]).mass2();
00107 s13=(_pd[0]+_pd[2]).mass2();
00108 s23=(_pd[1]+_pd[2]).mass2();
00109 temp=amps1(s12,0,1)+amps1(s13,0,2)+amps1(s23,1,2);
00110 return temp;
00111 }
00112
00113
00114
00115
00116 class rhopifull{
00117 public:
00118 rhopifull(EvtVector4R pd1, EvtVector4R pd2,EvtVector4R pd3){
00119 _pd[0]=pd1;
00120 _pd[1]=pd2;
00121 _pd[2]=pd3;
00122 }
00123 double Fij(int i, int j, double r);
00124 double R00(double r);
00125 EvtComplex amps1(int m, int i, int j);
00126 double amps( );
00127
00128 private:
00129 EvtVector4R _pd[3];
00130 };
00131
00132 double rhopifull::R00(double r){
00133 double mpi=0.1395;
00134 return r;
00135 }
00136
00137 double rhopifull::Fij(int i, int j, double r){
00138 double mpi=0.1395,mpsi=3.096916;
00139 double temp=mpsi*r;
00140 if(i==0&& j==0) return 0;
00141 if(i==1&& j==0) return temp;
00142 if(i==-1&& j==0) return -temp;
00143
00144 }
00145
00146 EvtComplex rhopifull::amps1(int m, int i, int j){
00147 double mrho=0.771,wrho=0.1492,s;
00148 EvtComplex img(0.0,1.0),dpro;
00149 EvtVector4R prho;
00150 prho=_pd[i]+_pd[j];
00151 s=prho.mass2();
00152 dpro=s-mrho*mrho+img*sqrt(s)*wrho;
00153 EvtHelSys angles(prho,_pd[i]),labAngles;
00154 double theta,phi,ct1,st1,phi1,st,ct;
00155 double rpp=angles.getHelAng(0);
00156 theta=angles.getHelAng(1);
00157 phi =angles.getHelAng(2);
00158 ct1 =labAngles.Angles(prho,1);
00159 phi1=labAngles.Angles(prho,2);
00160 int lamb;
00161 EvtComplex temp(0.0,0.0);
00162 for(lamb=-1;lamb<=1;lamb++) temp=temp+Fij(lamb,0,prho.d3mag())*Djmn(1,m,lamb,phi1,ct1,0.0)/dpro*R00(rpp)*Djmn(1,lamb,0,phi,theta,0.0);
00163 return temp;
00164 }
00165
00166 double rhopifull::amps(){
00167 double temp=0.0;
00168 int m;
00169 for(m=-1;m<=1;m+=2) temp=temp+pow(abs(amps1(m,0,1)+amps1(m,0,2)+amps1(m,1,2)),2.);
00170 return temp;
00171 }
00172
00173
00174
00176 double EvtDIY::AmplitudeSquare(){
00177 EvtVector4R dp1=GetDaugMomLab(0);
00178 EvtVector4R dp2=GetDaugMomLab(1),dp3=GetDaugMomLab(2);
00179
00180
00181
00182
00183
00184
00185
00186 rhopifull Jpsi2rhopi(dp1,dp2,dp3);
00187 return Jpsi2rhopi.amps();
00188 }
00189