00001 #include "T3piCrossPart.h"
00002 #include "TConstant.h"
00003
00004 #include <complex>
00005 #include <iostream>
00006
00007 using namespace rb;
00008 typedef std::complex<double> complex_t;
00009
00010 inline complex_t operator+(const double &x, const complex_t &y) {
00011 return y + x;
00012 }
00013
00014 inline complex_t operator-(const double &x, const complex_t &y) {
00015 return -(y - x);
00016 }
00017
00018 inline double Gamma(const double Q2, const double M2, const double G){
00019 double beta2 = (Q2-4*mpi2)/(M2-4*mpi2);
00020 double beta = sqrt(beta2);
00021 return G*M2/Q2*beta2*beta;
00022 }
00023
00024
00025 inline complex_t Rrho(double Q2, bool IsIsospinBreak){
00026 if(Q2<4*mpi2) return 0;
00027 const complex_t I(0,1);
00028 double Q = sqrt(Q2);
00029 complex_t prop_rho = Mrho2/(Q2 - Mrho2 + I*(Gamma(Q2,Mrho2,Grho)*Q));
00030 return prop_rho;
00031 complex_t prop_rhop = Mrhop2/(Q2 - Mrhop2 + I*(Gamma(Q2,Mrhop2,Grhop)*Q));
00032 complex_t prop_rhopp = Mrhopp2/(Q2 - Mrhopp2 + I*(Gamma(Q2,Mrhopp2,Grhopp)*Q));
00033 complex_t prop_rhoppp = Mrhoppp2/(Q2 - Mrhoppp2 + I*(Gamma(Q2,Mrhoppp2,Grhoppp)*Q));
00034
00035
00036
00037
00038 const complex_t beta(-0.178042,-0.369325);
00039
00040
00041 const complex_t gamma(-0.068488,0.0280249);
00042
00043 const complex_t xi(0.05,0.0);
00044 if(IsIsospinBreak){
00045
00046 const complex_t delta(0.002,0);
00047 complex_t prop_omega = Momega2/(Q2 - Momega2 + I*(Gomega*Momega));
00048
00049 return (prop_rho*(1-delta*Q2/Momega2*prop_omega) + beta*prop_rhop + gamma*prop_rhopp + xi*prop_rhoppp)/(1+beta+gamma+xi);
00050 } else {
00051 return (prop_rho + beta*prop_rhop + gamma*prop_rhopp + xi*prop_rhoppp)/(1+beta+gamma+xi);
00052 }
00053 }
00054
00055 inline complex_t Amp(double s){
00056 const complex_t I(0,1);
00057 const double Apsi = 4*alpha*Gpsi/(3*Geepsi);
00058 const double phipsi = -M_PI/2;
00059 const double Apsip = alpha*Gpsip/(3*Geepsip);
00060 const double phipsip = -M_PI/2;
00061 const double Apsipp = alpha*Gpsipp/(3*Geepsipp);
00062 const double phipsipp = -M_PI/2;
00063 complex_t res = 3*sqrt(s)/alpha*
00064 (
00065 Geepsi*(1 + Apsi*exp(I*phipsi))/(s - Mpsi2 + I*Mpsi*Gpsi) +
00066 Geepsip*(1 + Apsip*exp(I*phipsip))/(s - Mpsip2 + I*Mpsip*Gpsip) +
00067 Geepsipp*(1 + Apsipp*exp(I*phipsipp))/(s - Mpsipp2 + I*Mpsipp*Gpsipp));
00068 return res;
00069 }
00070
00071 inline complex_t AmpOmegaPhi(double s){
00072
00073
00074
00075
00076 const double comega = 6.23657803940740774e-01;
00077 const double cphi = -4.63075347511148622e-02;
00078 complex_t Romega = Momega2/complex_t(s-Momega2,Gomega*Momega);
00079 complex_t Rphi = Mphi2/complex_t(s-Mphi2,Gphi*Mphi);
00080 return (comega*Romega + cphi*Rphi);
00081 }
00082
00083 complex_t F3pi(double s, double Q02, double Qp2, double Qm2){
00084
00085 complex_t H = Rrho(Q02,true) + Rrho(Qm2,false) + Rrho(Qp2,false);
00086 return (5e-3*(1/s))*H*(1+Amp(s));
00087
00088
00089
00090
00091
00092
00093
00094 }
00095
00096 T3piCrossPart::T3piCrossPart(double e, double de, double nth0):
00097 TCrossPart(e,de,nth0)
00098 {
00099 const double m[] = { mpi, mpi, mpi0};
00100 const int pid[] = { 211, -211, 111};
00101 double emax = e*(1 - pow((m[0]+m[1]+m[2])/(2*e),2));
00102
00103 fphot->SetEnergyRange(de,emax);
00104 std::cout<<de<<"MeV <E_gamma< "<<emax<<" MeV"<<std::endl;
00105
00106 SetFinalParticles(3, m, pid);
00107 std::cout<<"Process -- e^+e^- -> rho pi -> pi^+pi^-pi^0"<<std::endl;
00108 }
00109
00110 void T3piCrossPart::SetJ(){
00111
00112 TLorentzVector &qp = *fres[0];
00113 TLorentzVector &qm = *fres[1];
00114 TLorentzVector &q0 = *fres[2];
00115
00116 double qpqm = qp*qm;
00117 double qpq0 = qp*q0;
00118 double qmq0 = qm*q0;
00119
00120
00121 double Q02 = 2*(mpi2 + qpqm);
00122 double Qp2 = mpi2 + mpi02 + 2*qmq0;
00123 double Qm2 = mpi2 + mpi02 + 2*qpq0;
00124
00125 complex_t F = F3pi(fq2,Q02,Qp2,Qm2);
00126
00127 J3PseudoScalars();
00128
00129 fJc *= F;
00130 }
00131
00132 bool T3piCrossPart::Accepted(){
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155 return true;
00156 }
00157