/**************************************************************************/ /**************************************************************************/ /* This program calculates the kaon photo/electro production */ /* differential cross sections (e,e' K+ Lambda), (e,e' K+ Sigma0), */ /* on a nucleon */ /* at high s and low t (t-channel Regge exchanges) */ /* or high s and low u (u-channel Regge exchanges). */ /* */ /* All dimensional quantities in this program are in units GeV. */ /* */ /**************************************************************************/ /**************************************************************************/ /*************************************************************/ /* */ /* Authors : Marc Vanderhaeghen & Michel Guidal */ /* */ /* First run : 24/04/1997 */ /* */ /* Last update : 18/10/1999 */ /* */ /*************************************************************/ #include #include #include #include #include #include #include #include #include #include /**************************************************************************/ /**************************************************************************/ double dsigma_dt(int reaction, int mech, int deg_pi, int deg_kstar, double el_hel, double phi_min, double phi_max, double epsilon, double Q_sqr, double s, double t); double int_dsigma_dt(double phi); double dsigma_dtdphi(int reaction, int mech, int deg_pi, int deg_kstar, double el_hel, double phi, double epsilon, double Q_sqr, double s, double t); double dsigma_dt_anal(int reaction, int mech, int deg_pi, int deg_kstar, double el_hel, double epsilon, double Q_sqr, double s, double t); double dsigmal_dt(int reaction, int mech, int deg_pi, int deg_kstar, double el_hel, double phi_min, double phi_max, double epsilon, double Q_sqr, double s, double t); double dsigmatt_dt(int reaction, int mech, int deg_pi, int deg_kstar, double el_hel, double phi_min, double phi_max, double epsilon, double Q_sqr, double s, double t); double dsigmatl_dt(int reaction, int mech, int deg_pi, int deg_kstar, double el_hel, double phi_min, double phi_max, double epsilon, double Q_sqr, double s, double t); double int_dsigmal_dt(double phi); double int_dsigmatt_dt(double phi); double int_dsigmatl_dt(double phi); double dsigmal_dtdphi(int reaction, int mech, int deg_pi, int deg_kstar, double el_hel, double phi, double epsilon, double Q_sqr, double s, double t); double dsigmatt_dtdphi(int reaction, int mech, int deg_pi, int deg_kstar, double el_hel, double phi, double epsilon, double Q_sqr, double s, double t); double dsigmatl_dtdphi(int reaction, int mech, int deg_pi, int deg_kstar, double el_hel, double phi, double epsilon, double Q_sqr, double s, double t); double sigma_trans_tot(int reaction, int mech, int deg_k, int deg_kstar, double Q_sqr, double s, double tmin, double tmax); double int_sigma_trans_tot(double t); double pol_tot(int reaction, int mech, int deg_k, int deg_kstar, double Q_sqr, double s, double epsilon, double phi, int iintegr, int polbeam, int xyz); double int_pol(double t); /**************************************************************************/ /**************************************************************************/ double klam_t_matrix(int mech, int deg_k, int deg_kstar, double el_hel, double sp_in, double sp_out, double phi, vec k_in, vec k_out, vec p_in, vec p_out, vec kaon_out); double klam_hel_spin_diff_cross_lab(int mech,int deg_k, int deg_kstar, double el_hel, double sp_out, double k_in_mom_lab, double k_out_mom_lab, double k_out_th_lab, double kaon_th, double kaon_phi); double v_klam_5fold_diff_cross_test(int reaction, int mech, int deg_k, int deg_kstar, double el_hel, double k_in_mom_lab, double k_out_mom_lab, double k_out_th_lab, double kaon_th, double kaon_phi); /**************************************************************************/ /**************************************************************************/ double response_R_T(int reaction, int mech, int deg_pi, int deg_kstar, double Q_sqr, double s, double t); double response_R_L(int reaction, int mech, int deg_pi, int deg_kstar, double Q_sqr, double s, double t); double response_R_TT(int reaction, int mech, int deg_pi, int deg_kstar, double Q_sqr, double s, double t); double response_R_TL(int reaction, int mech, int deg_pi, int deg_kstar, double Q_sqr, double s, double t); double response_R_TLp(int reaction, int mech, int deg_pi, int deg_kstar, double Q_sqr, double s, double t); double response_R_T_ortho(int reaction, int mech, int deg_pi, int deg_kstar, double Q_sqr, double s, double t); double response_R_L_ortho(int reaction, int mech, int deg_pi, int deg_kstar, double Q_sqr, double s, double t); double response_R_TL_ortho(int reaction, int mech, int deg_pi, int deg_kstar, double Q_sqr, double s, double t); double response_R_TT_ortho(int reaction, int mech, int deg_pi, int deg_kstar, double Q_sqr, double s, double t); double response_R_TLp_ortho(int reaction, int mech, int deg_pi, int deg_kstar, double Q_sqr, double s, double t); double response_R_TT_par(int reaction, int mech, int deg_pi, int deg_kstar, double Q_sqr, double s, double t); double response_R_TL_par(int reaction, int mech, int deg_pi, int deg_kstar, double Q_sqr, double s, double t); double response_R_TLp_par(int reaction, int mech, int deg_pi, int deg_kstar, double Q_sqr, double s, double t); double response_R_Tp_par(int reaction, int mech, int deg_pi, int deg_kstar, double Q_sqr, double s, double t); double response_R_TT_long(int reaction, int mech, int deg_pi, int deg_kstar, double Q_sqr, double s, double t); double response_R_TL_long(int reaction, int mech, int deg_pi, int deg_kstar, double Q_sqr, double s, double t); double response_R_TLp_long(int reaction, int mech, int deg_pi, int deg_kstar, double Q_sqr, double s, double t); double response_R_Tp_long(int reaction, int mech, int deg_pi, int deg_kstar, double Q_sqr, double s, double t); double recoil_asymm(int reaction, int mech, int deg_pi, int deg_kstar, double Q_sqr, double s, double t); dcomplex J_hadron(int reaction, int mech, int deg_pi, int deg_kstar, double Q_sqr, double s, double t, double sp_out, int pol_ph, double sp_in); /**************************************************************************/ /**************************************************************************/ dcomplex kaon_regge_trajectory(int deg_pi, double s, double t); dcomplex kstar_regge_trajectory(int deg_kstar, double s, double t); /*dcomplex omega_regge_trajectory(double s, double t);*/ /* ********************************************************************* */ double kaon_emformfactor(double Q_sqr, double t); double kstarkgamma_formfactor(double Q_sqr, double t); /* ********************************************************************* */ dcomplex t_vgap_klam_regge_phi(int mech, int deg_pi, int deg_kstar, int pol_ph, double sp_out, double sp_in, double Q_sqr, double s, double t, double phi); dcomplex t_vgap_klam_regge_forward(int mech, int deg_k, int deg_kstar, int pol_ph, double sp_out, double sp_in, double Q_sqr, double s, double t); /* ********************************************************************* */ dcomplex t_vgap_ksigma_regge_forward(int mech, int deg_pi, int deg_kstar, int pol_ph, double sp_out, double sp_in, double Q_sqr, double s, double t); /* ********************************************************************* */ /* ********************************************************************* */ FILE *fp, *fp1, *fp2, *fp3, *fp4, *fp5; static int reaction_s, mech_s, deg_k_s, deg_kstar_s, ichoicek, ichoicekst, isys, targrec, iintegr_s, polbeam_s, xyz_s, rep_s, k0sig_s; static double el_hel_s, epsilon_s, Q_sqr_s, s_s, t_s, mscalk, mscalkst, c_offshell, epsilon_s, phi_s; main() { int reaction; double Q_sqr, phi_min, phi_max, epsilon, w, s, t, u, mt, phi, diff, ga_mom_cm, ga_en_cm, k_mom_cm, k_en, k_th, difft, diffl, difftt, difftl; char fileb[80], file2[80]; /* ********************************************************************* */ /* ********************************************************************* */ /* */ /* Choose the reaction you want to calculate : */ /* 1 = (e, e' K Lambda) */ /* 2 = (e, e' K Sigma) */ /* */ /* ********************************************************************* */ /* ********************************************************************* */ printf("\n\n Choose the reaction you want to calculate :"); printf("\n 1 = (e, e' K+ Lambda)"); printf("\n 2 = (e, e' K+ Sigma0)"); printf("\n 3 = (e, e' K0 Sigma+)"); printf("\n"); scanf("%d", &reaction); k0sig_s=0; switch(reaction) { /***********************************************************************/ /* */ /* (e, e' K LAMBDA) reaction */ /* */ /********************************************************************* */ case 1: /* for (e, e' K+ Lambda) reaction */ { /* 1. (e, e' K+ Lambda) differential cross section : forward */ int pol = 1; switch(pol) { case 1: /* (e, e' K Lambda) reaction differential cross section : forward */ { int mech, deg_k, deg_kstar, separated; double res_trans, res_long, res_tt, res_tl, res_tlp; printf("\n\n Give choice for FORWARD scattering mechanism"); printf("\n 1 = t-channel kaon (no form factor)"); printf("\n 2 = kstar"); printf("\n 3 = kaon gaugeinvariant"); printf("\n 4 = kaon gaugeinvariant + kstar"); printf("\n"); scanf("%d", &mech); printf("\n\n Give choice for degeneracy of KAON trajectory"); printf("\n 1 = non-degenerate trajectory"); printf("\n 2 = degenerate trajectory with rotating phase (default)"); printf("\n 3 = degenerate trajectory with NON-rotating phase"); printf("\n 4 = single pole"); printf("\n"); scanf("%d", °_k); printf("\n\n Give choice for degeneracy of KSTAR trajectory"); printf("\n 1 = non-degenerate trajectory"); printf("\n 2 = degenerate trajectory with rotating phase (default)"); printf("\n 3 = degenerate trajectory with NON-rotating phase"); printf("\n 4 = single pole"); printf("\n"); scanf("%d", °_kstar); switch(mech) { case 1 : { break; } case 2 : { if (deg_kstar == 2) strcpy(fileb,"_kstar_regge.dat"); if (deg_kstar == 4) strcpy(fileb,"_kstar_born.dat"); break; } case 3 : { if (deg_k == 2) strcpy(fileb,"_k_regge.dat"); if (deg_k == 4) strcpy(fileb,"_k_born.dat"); break; } case 4 : { if ((deg_k == 2) && (deg_kstar == 2)) strcpy(fileb,"_kkstar_regge.dat"); if ((deg_k == 4) && (deg_kstar == 4)) strcpy(fileb,"_kkstar_born.dat"); break; } } printf("\n\n Give choice for electromagnetic form factor of K"); printf("\n 1 = Monopole parametrization"); printf("\n 2 = Monopole at large Q_sqr, correct charge radius at low Q_sqr"); printf("\n 3 = VMD form (rho, omega, phi poles)"); printf("\n 4 = FF with off-shell extrapolation"); printf("\n"); scanf("%d", &ichoicek); printf("\n\n Give choice for electromagnetic form factor of K*"); printf("\n 1 = Monopole parametrization"); printf("\n 2 = Monopole at large Q_sqr, correct charge radius at low Q_sqr"); printf("\n 3 = VMD form (rho, omega, phi poles)"); printf("\n"); scanf("%d", &ichoicekst); if((ichoicek == 1)||(ichoicek == 2)) { printf("\n\n K form factor mass scale/cut-off"); printf("\n"); scanf("%lf", &mscalk); } if(ichoicek == 4) { printf("\n Give the value of the dimensionless parameter c in the kaon off-shell FF (e.g. 1) \n"); scanf("%lf",&c_offshell); } if((ichoicekst == 1)||(ichoicekst == 2)) { printf("\n\n K* form factor mass scale/cut-off"); printf("\n"); scanf("%lf", &mscalkst); } printf("\n\n Give choice for OBSERVABLE to calculate"); printf("\n 1 = photo- and electro-production : RESPONSE functions"); printf("\n 2 = theta dependence (Cornell experiments : Bebek et al.)"); printf("\n 3 = q2 dependence"); printf("\n 4 = recoil polarization"); printf("\n 5 = double polarization (method 1 ; response function formalism)"); printf("\n 6 = double polarization (method 2 ; full matrix element calculated for checking-)"); printf("\n 7 = cross sections (T,L,TT,LT) : user-defined kinematics"); printf("\n"); scanf("%d", &separated); switch(separated) { case 1 : /* separated cross section */ { char file1t[80]="ep_ekphyp_w2p14_qsqr1p19"; char file2t[80]; char trans[]="_trans_kkstar.dat"; char longi[]="_long_kkstar.dat"; char tt[]="_tt_pi_gaugeinv.dat"; char tl[]="_tl_pi_gaugeinv.dat"; char tlp[]="_tlp_k_gaugeinv.dat"; char file1l[80],file1tt[80],file1tl[80],file1tlp[80]; int photoelectro; strcpy(file1l,file1t); strcpy(file1tt,file1t); strcpy(file1tl,file1t); strcpy(file1tlp,file1t); printf("\n photo or electro ?"); printf("\n 1 = photo"); printf("\n 2 = electro"); printf("\n"); scanf("%d", &photoelectro); switch(photoelectro) { case 1 : /*photoproduction */ { int totdiff; printf("\n total or differential cross section ?"); printf("\n 1 = total"); printf("\n 2 = differential"); printf("\n"); scanf("%d", &totdiff); switch(totdiff) { case 1 : /*photoproduction : total cross section -as function of photon energy-*/ { char file1[80]="phototot.dat"; double ga_en_lab, result1, tmin, tmax; Q_sqr=0.; if ( (fp1 = fopen(file1,"w")) == NULL) fprintf(stderr,"\nUnable to open file"); printf("OUTPUT : E_gam(GeV), sig_tot(mubarn) \n"); for( ga_en_lab= 1.; ga_en_lab <= 4.; ga_en_lab += .25) { s=pow(NU_M,2.)+2.*ga_en_lab*NU_M; tmin=pow((-Q_sqr-pow(K_M,2.)-pow(NU_M,2.)+pow(LAM_M,2.))/2./sqrt(s),2.) -pow(sqrt(pow((s-Q_sqr-pow(NU_M,2.))/2./sqrt(s),2.)+Q_sqr) -sqrt(pow((s+pow(K_M,2.)-pow(LAM_M,2.))/2./sqrt(s),2.)-pow(K_M,2.)),2.); tmax=pow((-Q_sqr-pow(K_M,2.)-pow(NU_M,2.)+pow(LAM_M,2.))/2./sqrt(s),2.) -pow(sqrt(pow((s-Q_sqr-pow(NU_M,2.))/2./sqrt(s),2.)+Q_sqr) +sqrt(pow((s+pow(K_M,2.)-pow(LAM_M,2.))/2./sqrt(s),2.)-pow(K_M,2.)),2.); result1 = sigma_trans_tot(reaction, mech, deg_k, deg_kstar, Q_sqr, s, tmin, tmax); /* printf("\n%f \t %f \t %f \t %f \t %e", Q_sqr, s, tmin, tmax, result1);*/ fprintf(fp1,"%f \t %e \n", ga_en_lab, result1); } fclose(fp1); break; } case 2 : /*photoproduction : differential cross section*/ { double ga_en_lab, bidon1, cos_k_th; int wort; char file1[80]; char file2[80]; printf("\n As a function of t or W ?"); printf("\n 1: As a function of t"); printf("\n 2: As a function of W"); printf("\n"); scanf("%d", &wort); if (wort==1) { printf("\n photon energy (in GeV) ?"); printf("\n"); scanf("%lf", &ga_en_lab); bidon1=(ga_en_lab-(int) ga_en_lab +.01)*10.; sprintf(file1,"gap_kl_diff_%dp%dgev.dat", (int) ga_en_lab, (int) bidon1); sprintf(file2,"gap_kl_ass_%dp%dgev.dat", (int) ga_en_lab, (int) bidon1); s = 2. * NU_M * ga_en_lab + pow(NU_M,2.); Q_sqr = 0.; ga_mom_cm = sqrt( ( pow(s - pow(NU_M,2.) + Q_sqr, 2.) + pow(2. * NU_M, 2.) * Q_sqr ) / (4. * s) ); ga_en_cm = sqrt(s) - sqrt(pow(ga_mom_cm,2.) + pow(NU_M,2.)); k_mom_cm = sqrt( ( pow(s - pow(LAM_M,2.) - pow(K_M,2.), 2.) - pow(2. * LAM_M * K_M, 2.) ) / (4. * s) ); k_en = sqrt( pow(k_mom_cm, 2.) + pow(K_M, 2.) ); if ( (fp1 = fopen(file1,"w")) == NULL) fprintf(stderr,"\nUnable to open file"); if ( (fp2 = fopen(file2,"w")) == NULL) fprintf(stderr,"\nUnable to open file"); printf("OUTPUT1 : -t(GeV2), dsig_T/dt(mubarn/GeV2) \n"); printf("OUTPUT2 : -t(GeV2), -dsig_TT/dt/dsig_T/dt \n"); for(mt = 0.001; mt <= 1.5; mt += 0.02) { cos_k_th = (- mt - pow(K_M,2.) + 2. * ga_en_cm * k_en) / (2. * ga_mom_cm * k_mom_cm) ; res_trans = response_R_T(reaction, mech, deg_k, deg_kstar, Q_sqr, s, - mt); res_long = response_R_TT(reaction, mech, deg_k, deg_kstar, Q_sqr, s, - mt); res_long = res_long / res_trans; fprintf(fp1, "%f \t %f \t %e \n", mt, cos_k_th, res_trans); fprintf(fp2, "%f \t %f \t %e \n", mt, cos_k_th, res_long); } fclose(fp2); } if (wort==2) { double angle, conv; printf("\n Angle (in degrees) ?"); printf("\n"); scanf("%lf", &angle); bidon1=(angle-(int) angle +.01)*10.; sprintf(file1,"gap_kl_diff_%dp%ddeg.dat", (int) angle, (int) bidon1); Q_sqr = 0.; if ( (fp1 = fopen(file1,"w")) == NULL) fprintf(stderr,"\nUnable to open file"); printf("OUTPUT1 : -t(GeV2), dsig_T/dt(mubarn/GeV2), dsig_T/dcos(thcm) (mubarn/sr) \n"); for(w = 1.61; w <= 2.4; w += 0.1) { s=pow(w,2.); ga_mom_cm = sqrt( ( pow(s - pow(NU_M,2.) + Q_sqr, 2.) + pow(2. * NU_M, 2.) * Q_sqr ) / (4. * s) ); ga_en_cm = sqrt(s) - sqrt(pow(ga_mom_cm,2.) + pow(NU_M,2.)); k_mom_cm = sqrt( ( pow(s - pow(LAM_M,2.) - pow(K_M,2.), 2.) - pow(2. * LAM_M * K_M, 2.) ) / (4. * s) ); k_en = sqrt( pow(k_mom_cm, 2.) + pow(K_M, 2.) ); mt = - 2. * ga_mom_cm * k_mom_cm *cos(angle*PI/180.) - pow(K_M,2.) + 2. * ga_en_cm * k_en ; conv = 2. * k_mom_cm * ga_mom_cm ; res_trans = response_R_T(reaction, mech, deg_k, deg_kstar, Q_sqr, s, - mt); res_long = response_R_TT(reaction, mech, deg_k, deg_kstar, Q_sqr, s, - mt); res_long = res_long / res_trans; fprintf(fp1, "%f \t %f \t %e \n", w, res_trans, conv * res_trans); } } fclose(fp1); break; } } break; } case 2 : /*electroproduction*/ { int sep; printf("\n Which kinematics ?"); printf("\n Ratio R_L/R_T : "); printf("\n 1 : W = 2.14 GeV, Q^2 = 1.19 GeV^2 (Bebek et al.)"); printf("\n 2 : W = 2.65 GeV, Q^2 = 1.99 GeV^2 (Bebek et al.)"); printf("\n 3 : W = 2.56 GeV, Q^2 = 3.38 GeV^2 (Bebek et al.)"); printf("\n R_T+epsilon R_L, R_TT, R_TL : "); printf("\n 4 : W = 2.21 GeV, Q^2 = .7 GeV^2, epsilon =.85 (Brauel et al.)"); printf("\n 5 : W = 2.21 GeV, Q^2 = 1.35 GeV^2, epsilon =.82 (Brauel et al.)"); printf("\n 6 : W = 2.21 GeV, Q^2 = 0.06 GeV^2, epsilon =.41 (Brauel et al.)"); printf("\n"); scanf("%d", &sep); switch(sep) { case 1 : /*ratio rlrt*/ { char file1[80]="ep_ekl_w2p14_qsqr1p19_rap"; s = pow(2.14, 2.); Q_sqr = 1.19; if ( (fp1 = fopen(strcat(file1,fileb),"w")) == NULL) fprintf(stderr,"\nUnable to open file"); printf("OUTPUT : theta_cm(deg), -t(GeV2), L/T \n"); ga_mom_cm = sqrt( ( pow(s - pow(NU_M,2.) + Q_sqr, 2.) + pow(2. * NU_M, 2.) * Q_sqr ) / (4. * s) ); ga_en_cm = sqrt(s) - sqrt(pow(ga_mom_cm,2.) + pow(NU_M,2.)); k_mom_cm = sqrt( ( pow(s - pow(LAM_M,2.) - pow(K_M,2.), 2.) - pow(2. * LAM_M * K_M, 2.) ) / (4. * s) ); k_en = sqrt( pow(k_mom_cm, 2.) + pow(K_M, 2.) ); for(k_th = 0.001; k_th <= 21.; k_th += 2.) { t = - Q_sqr + pow(K_M,2.) - 2. * ga_en_cm * k_en + 2. * ga_mom_cm * k_mom_cm * cos(k_th * PI / 180.); res_trans = response_R_T(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t); res_long = response_R_L(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t); res_long = res_long / res_trans; fprintf(fp1, "%f \t %f \t %e \n", k_th, - t, res_long); } fclose(fp1); break; } case 2 : /*ratio rlrt*/ { char file1[80]="ep_ekl_w2p65_qsqr1p99_rap"; s = pow(2.65, 2.); Q_sqr = 1.99; if ( (fp1 = fopen(strcat(file1,fileb),"w")) == NULL) fprintf(stderr,"\nUnable to open file"); printf("OUTPUT : theta_cm(deg), -t(GeV2), L/T \n"); ga_mom_cm = sqrt( ( pow(s - pow(NU_M,2.) + Q_sqr, 2.) + pow(2. * NU_M, 2.) * Q_sqr ) / (4. * s) ); ga_en_cm = sqrt(s) - sqrt(pow(ga_mom_cm,2.) + pow(NU_M,2.)); k_mom_cm = sqrt( ( pow(s - pow(LAM_M,2.) - pow(K_M,2.), 2.) - pow(2. * LAM_M * K_M, 2.) ) / (4. * s) ); k_en = sqrt( pow(k_mom_cm, 2.) + pow(K_M, 2.) ); for(k_th = 0.001; k_th <= 21.; k_th += 2.) { t = - Q_sqr + pow(K_M,2.) - 2. * ga_en_cm * k_en + 2. * ga_mom_cm * k_mom_cm * cos(k_th * PI / 180.); res_trans = response_R_T(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t); res_long = response_R_L(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t); res_long = res_long / res_trans; fprintf(fp1, "%f \t %f \t %e \n", k_th, - t, res_long); } fclose(fp1); break; } case 3 : /*ratio rlrt*/ { char file1[80]="ep_ekl_w2p56_qsqr3p38_rap"; s = pow(2.56, 2.); Q_sqr = 3.38; if ( (fp1 = fopen(strcat(file1,fileb),"w")) == NULL) fprintf(stderr,"\nUnable to open file"); printf("OUTPUT : theta_cm(deg), -t(GeV2), L/T \n"); ga_mom_cm = sqrt( ( pow(s - pow(NU_M,2.) + Q_sqr, 2.) + pow(2. * NU_M, 2.) * Q_sqr ) / (4. * s) ); ga_en_cm = sqrt(s) - sqrt(pow(ga_mom_cm,2.) + pow(NU_M,2.)); k_mom_cm = sqrt( ( pow(s - pow(LAM_M,2.) - pow(K_M,2.), 2.) - pow(2. * LAM_M * K_M, 2.) ) / (4. * s) ); k_en = sqrt( pow(k_mom_cm, 2.) + pow(K_M, 2.) ); for(k_th = 0.001; k_th <= 21.; k_th += 2.) { t = - Q_sqr + pow(K_M,2.) - 2. * ga_en_cm * k_en + 2. * ga_mom_cm * k_mom_cm * cos(k_th * PI / 180.); res_trans = response_R_T(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t); res_long = response_R_L(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t); res_long = res_long / res_trans; fprintf(fp1, "%f \t %f \t %e \n", k_th, - t, res_long); } fclose(fp1); break; } case 4 : /*separation rt+eps*rl, rtt, rtl*/ { char file1[80]="ep_ekl_w2p21_qsqr0p7_tpepsl"; char file2[80]="ep_ekl_w2p21_qsqr0p7_tt"; char file3[80]="ep_ekl_w2p21_qsqr0p7_tl"; s = pow(2.21, 2.); Q_sqr = .7; epsilon = .85; if ( (fp1 = fopen(strcat(file1,fileb),"w")) == NULL) fprintf(stderr,"\nUnable to open file"); if ( (fp2 = fopen(strcat(file2,fileb),"w")) == NULL) fprintf(stderr,"\nUnable to open file"); if ( (fp3 = fopen(strcat(file3,fileb),"w")) == NULL) fprintf(stderr,"\nUnable to open file"); printf("OUTPUT1 : -t(GeV2), T+eps*L(mubarn/GeV2) \n"); printf("OUTPUT2 : -t(GeV2), -TT(mubarn/GeV2) \n"); printf("OUTPUT3 : -t(GeV2), TL(mubarn/GeV2) \n"); ga_mom_cm = sqrt( ( pow(s - pow(NU_M,2.) + Q_sqr, 2.) + pow(2. * NU_M, 2.) * Q_sqr ) / (4. * s) ); ga_en_cm = sqrt(s) - sqrt(pow(ga_mom_cm,2.) + pow(NU_M,2.)); k_mom_cm = sqrt( ( pow(s - pow(LAM_M,2.) - pow(K_M,2.), 2.) - pow(2. * LAM_M * K_M, 2.) ) / (4. * s) ); k_en = sqrt( pow(k_mom_cm, 2.) + pow(K_M, 2.) ); for(mt = 0.01; mt <= .7; mt += 0.1) { res_trans = response_R_T(reaction, mech, deg_k, deg_kstar, Q_sqr, s, -mt); res_long = response_R_L(reaction, mech, deg_k, deg_kstar, Q_sqr, s, -mt); res_trans = res_trans+epsilon*res_long; res_tt = response_R_TT(reaction, mech, deg_k, deg_kstar, Q_sqr, s, - mt); res_tl = response_R_TL(reaction, mech, deg_k, deg_kstar, Q_sqr, s, - mt); /* printf("%f \t %e \n", mt, res_long);*/ fprintf(fp1, "%f \t %e \n", mt, res_trans); fprintf(fp2, "%f \t %e \n", mt, res_tt); fprintf(fp3, "%f \t %e \n", mt, res_tl); } fclose(fp1); fclose(fp2); fclose(fp3); break; } case 5 : /*separation rt+eps*rl, rtt, rtl*/ { char file1[80]="ep_ekl_w2p21_qsqr1p35_tpepsl"; char file2[80]="ep_ekl_w2p21_qsqr1p35_tt"; char file3[80]="ep_ekl_w2p21_qsqr1p35_tl"; s = pow(2.21, 2.); Q_sqr = 1.35; epsilon = .82; if ( (fp1 = fopen(strcat(file1,fileb),"w")) == NULL) fprintf(stderr,"\nUnable to open file"); if ( (fp2 = fopen(strcat(file2,fileb),"w")) == NULL) fprintf(stderr,"\nUnable to open file"); if ( (fp3 = fopen(strcat(file3,fileb),"w")) == NULL) fprintf(stderr,"\nUnable to open file"); printf("OUTPUT1 : -t(GeV2), T+eps*L(mubarn/GeV2) \n"); printf("OUTPUT2 : -t(GeV2), -TT(mubarn/GeV2) \n"); printf("OUTPUT3 : -t(GeV2), TL(mubarn/GeV2) \n"); ga_mom_cm = sqrt( ( pow(s - pow(NU_M,2.) + Q_sqr, 2.) + pow(2. * NU_M, 2.) * Q_sqr ) / (4. * s) ); ga_en_cm = sqrt(s) - sqrt(pow(ga_mom_cm,2.) + pow(NU_M,2.)); k_mom_cm = sqrt( ( pow(s - pow(LAM_M,2.) - pow(K_M,2.), 2.) - pow(2. * LAM_M * K_M, 2.) ) / (4. * s) ); k_en = sqrt( pow(k_mom_cm, 2.) + pow(K_M, 2.) ); for(mt = 0.01; mt <= .7; mt += 0.1) { res_trans = response_R_T(reaction, mech, deg_k, deg_kstar, Q_sqr, s, -mt); res_long = response_R_L(reaction, mech, deg_k, deg_kstar, Q_sqr, s, -mt); res_trans = res_trans+epsilon*res_long; res_tt = response_R_TT(reaction, mech, deg_k, deg_kstar, Q_sqr, s, - mt); res_tl = response_R_TL(reaction, mech, deg_k, deg_kstar, Q_sqr, s, - mt); fprintf(fp1, "%f \t %e \n", mt, res_trans); fprintf(fp2, "%f \t %e \n", mt, res_tt); fprintf(fp3, "%f \t %e \n", mt, res_tl); } fclose(fp1); fclose(fp2); fclose(fp3); break; } case 6 : /*separation rt+eps*rl, rtt, rtl*/ { char file1[80]="ep_ekl_w2p21_qsqr0p06_tpepsl"; char file2[80]="ep_ekl_w2p21_qsqr0p06_tt"; char file3[80]="ep_ekl_w2p21_qsqr0p06_tl"; s = pow(2.21, 2.); Q_sqr = .06; epsilon = .41; if ( (fp1 = fopen(strcat(file1,fileb),"w")) == NULL) fprintf(stderr,"\nUnable to open file"); if ( (fp2 = fopen(strcat(file2,fileb),"w")) == NULL) fprintf(stderr,"\nUnable to open file"); if ( (fp3 = fopen(strcat(file3,fileb),"w")) == NULL) fprintf(stderr,"\nUnable to open file"); printf("OUTPUT1 : -t(GeV2), T+eps*L(mubarn/GeV2) \n"); printf("OUTPUT2 : -t(GeV2), -TT(mubarn/GeV2) \n"); printf("OUTPUT3 : -t(GeV2), TL(mubarn/GeV2) \n"); ga_mom_cm = sqrt( ( pow(s - pow(NU_M,2.) + Q_sqr, 2.) + pow(2. * NU_M, 2.) * Q_sqr ) / (4. * s) ); ga_en_cm = sqrt(s) - sqrt(pow(ga_mom_cm,2.) + pow(NU_M,2.)); k_mom_cm = sqrt( ( pow(s - pow(LAM_M,2.) - pow(K_M,2.), 2.) - pow(2. * LAM_M * K_M, 2.) ) / (4. * s) ); k_en = sqrt( pow(k_mom_cm, 2.) + pow(K_M, 2.) ); for(mt = 0.01; mt <= .7; mt += 0.1) { res_trans = response_R_T(reaction, mech, deg_k, deg_kstar, Q_sqr, s, -mt); res_long = response_R_L(reaction, mech, deg_k, deg_kstar, Q_sqr, s, -mt); /* res_trans = res_trans+epsilon*res_long;*/ res_tt = response_R_TT(reaction, mech, deg_k, deg_kstar, Q_sqr, s, - mt); res_tl = response_R_TL(reaction, mech, deg_k, deg_kstar, Q_sqr, s, - mt); fprintf(fp1, "%f \t %e \n", mt, res_trans); fprintf(fp2, "%f \t %e \n", mt, res_tt); fprintf(fp3, "%f \t %e \n", mt, res_tl); } fclose(fp1); fclose(fp2); fclose(fp3); break; } } } } break; } case 2 : /* theta dependence (Cornell and DESY experiments : Bebek et al., Brauel et al.) */ { int i, kin, choiceintegr; double mt_min,cos_min,lbin; printf("\n\n Choose the kinematics of the Cornell experiments (Bebek et al.)"); printf("\n 1 : W = 2.2 GeV, Q^2 = 1.18 GeV^2, epsilon = 0.94 (Bebek et al.)"); printf("\n 2 : W = 2.17 GeV, Q^2 = 1.97 GeV^2, epsilon = 0.94 (Bebek et al.)"); printf("\n 3 : W = 2.19 GeV, Q^2 = 3.95 GeV^2, epsilon = 0.88 (Bebek et al.)"); printf("\n 4 : W = 2.78 GeV, Q^2 = 1.38 GeV^2, epsilon = 0.88 (Bebek et al.)"); printf("\n 5 : W = 2.47 GeV, Q^2 = 3.52 GeV^2, epsilon = 0.86 (Bebek et al.)"); printf("\n 6 : W = 2.66 GeV, Q^2 = 2.04 GeV^2, epsilon = 0.88 (Bebek et al.)"); printf("\n 7 : W = 2.21 GeV, Q^2 = .06 GeV^2, epsilon = 0.41 (Brauel et al.)"); printf("\n 8 : W = 2.21 GeV, Q^2 = .28 GeV^2, epsilon = 0.74 (Brauel et al.)"); printf("\n 9 : W = 2.21 GeV, Q^2 = .7 GeV^2, epsilon = 0.85 (Brauel et al.)"); printf("\n 10 : W = 2.21 GeV, Q^2 = 1.35 GeV^2, epsilon = 0.82 (Brauel et al.)"); printf("\n"); scanf("%d", &kin); /* printf("\n\n Integration on Phi (120=>240) or not?"); printf("\n 0 : no"); printf("\n 1 : yes"); printf("\n"); scanf("%d", &choiceintegr);*/ choiceintegr=1; switch(kin) { case 1 : { strcpy(file2,"ep_ekl_w2p2_qsqr1p18_eps0p94"); w=2.2; Q_sqr=1.18; epsilon=.94; break; } case 2 : { strcpy(file2,"ep_ekl_w2p17_qsqr1p97_eps0p94"); w=2.17; Q_sqr=1.97; epsilon=.94; break; } case 3 : { strcpy(file2,"ep_ekl_w2p19_qsqr3p95_eps0p88"); w=2.19; Q_sqr=3.95; epsilon=.88; break; } case 4 : { strcpy(file2,"ep_ekl_w2p78_qsqr1p38_eps0p88"); w=2.78; Q_sqr=1.38; epsilon=.88; break; } case 5 : { strcpy(file2,"ep_ekl_w2p47_qsqr3p52_eps0p86"); w=2.47; Q_sqr=3.52; epsilon=.86; break; } case 6 : { strcpy(file2,"ep_ekl_w2p66_qsqr2p04_eps0p88"); w=2.66; Q_sqr=2.04; epsilon=.88; break; } case 7 : { strcpy(file2,"ep_ekl_w2p21_qsqr0p06_eps0p41"); w=2.21; Q_sqr=.06; /* w=1.9; Q_sqr=.02;*/ epsilon=.41; break; } case 8 : { strcpy(file2,"ep_ekl_w2p21_qsqr0p28_eps0p74"); w=2.21; Q_sqr=.28; /* w=1.9; Q_sqr=.2;*/ epsilon=.74; break; } case 9 : { strcpy(file2,"ep_ekl_w2p21_qsqr0p70_eps0p85"); w=2.21; Q_sqr=.7; /* w=1.9; Q_sqr=.55;*/ epsilon=.85; break; } case 10 : { strcpy(file2,"ep_ekl_w2p21_qsqr1p35_eps0p82"); w=2.21; Q_sqr=1.35; /* w=1.9; Q_sqr=1.;*/ epsilon=.82; break; } } if ( (fp2 = fopen(strcat(file2,fileb),"w")) == NULL) fprintf(stderr,"\nUnable to open file"); s = pow(w, 2.); ga_mom_cm = sqrt( ( pow(s - pow(NU_M,2.) + Q_sqr, 2.) + pow(2. * NU_M, 2.) * Q_sqr ) / (4. * s) ); ga_en_cm = sqrt(s) - sqrt(pow(ga_mom_cm,2.) + pow(NU_M,2.)); k_mom_cm = sqrt( ( pow(s - pow(LAM_M,2.) - pow(K_M,2.), 2.) - pow(2. * LAM_M * K_M, 2.) ) / (4. * s) ); k_en = sqrt( pow(k_mom_cm, 2.) + pow(K_M, 2.) ); if (kin<7){ printf("OUTPUT : thetaK_cm(deg), -t(GeV2), T+eps*L(mubarn/sr) \n"); for(k_th = 0.001; k_th <= 21.; k_th += 2.) { double conv; t = - Q_sqr + pow(K_M,2.) - 2. * ga_en_cm * k_en + 2. * ga_mom_cm * k_mom_cm * cos(k_th * PI / 180.); conv = 2. * k_mom_cm * ga_mom_cm / (2.*PI) ; diff = conv * dsigma_dt_anal(reaction, mech, deg_k, deg_kstar, 0., epsilon, Q_sqr, s, t); fprintf(fp2, "%f \t %f \t %e \n", k_th, -t, diff); } fclose(fp2); } if ((kin>6)&&(kin<11)){ /* double conv; conv = 2. * k_mom_cm * ga_mom_cm / (2.*PI) ; printf("\n %f",conv);*/ printf("Below, the dsig/dtdphi's have been Phi-integrated from 120 to 240 deg. Therefore, outputs have been divided by Phi-bin width"); if (choiceintegr==0) printf("OUTPUT : -t(GeV2), 2PI*(T+eps*L+eps*cos(2*Phi)*TT+sqrt(2*eps*(1+eps))*cos(Phi)*TL)(mubarn/GeV2), 2PI*T(mubarn/GeV2), ignore, ignore, ignore \n"); if (choiceintegr==1) printf("OUTPUT : -t(GeV2), 2PI*(T+eps*L+eps*cos(2*Phi)*TT+sqrt(2*eps*(1+eps))*cos(Phi)*TL)(mubarn/GeV2), 2PI*T(mubarn/GeV2), 2PI*eps*L(mubarn/GeV2), 2PI*eps*TT(mubarn/GeV2), 2PI*sqrt(2*eps*(1+eps))*TL(mubarn/GeV2) \n"); printf("OUTPUT : -t(GeV2), 2PI*dsig_tot/dt/dphi, 2PI*dsig_T/dt/dphi, 2PI*dsig_L/dt/dphi, 2PI*dsig_TT/dt/dphi, 2PI*dsig_TL/dt/dphi(mubarn/sr) \n"); mt_min=0.; do{ mt_min=mt_min+.01; cos_min = (- Q_sqr + pow(K_M,2.) - 2. * ga_en_cm * k_en + mt_min)/ (-2. * ga_mom_cm * k_mom_cm) ; }while(fabs(cos_min)>1.); phi_min = 120. * PI / 180.; phi_max = 240. * PI / 180.; lbin=phi_max-phi_min; for(mt = mt_min; mt <= 1.7; mt += 0.2) { if (choiceintegr==0) { diff = 2.*PI* dsigma_dtdphi(reaction, mech, deg_k, deg_kstar, 0., PI ,epsilon, Q_sqr, s, -mt); difft = 2.*PI* dsigma_dtdphi(reaction, mech, deg_k, deg_kstar, 0., PI ,0., Q_sqr, s, -mt); } else { diff = 2. * PI * dsigma_dt(reaction, mech, deg_k, deg_kstar, 0., phi_min, phi_max, epsilon, Q_sqr, s, - mt); difft = 2. * PI * dsigma_dt(reaction, mech, deg_k, deg_kstar, 0., phi_min, phi_max, 0., Q_sqr, s, - mt); diffl = 2. * PI * dsigmal_dt(reaction, mech, deg_k, deg_kstar, 0., phi_min, phi_max, epsilon, Q_sqr, s, - mt); difftt = 2. * PI * dsigmatt_dt(reaction, mech, deg_k, deg_kstar, 0., phi_min, phi_max, epsilon, Q_sqr, s, - mt); difftl = 2. * PI * dsigmatl_dt(reaction, mech, deg_k, deg_kstar, 0., phi_min, phi_max, epsilon, Q_sqr, s, - mt); diff=diff/lbin; difft=difft/lbin; diffl=diffl/lbin; difftt=difftt/lbin; difftl=difftl/lbin; } if (diff!=0) fprintf(fp2, "%f \t %e \t %e \t %e \t %e \t %e \n", mt, diff, difft, diffl, difftt, difftl); } fclose(fp2); } break; } case 3 : /* q2 dependence */ { int obsq2; printf("\n\n Which observable ?"); printf("\n 1 = dsig/dt (Bebek et al.)"); printf("\n 2 = dsigT/dt, dsigL/dt, sigL/sigT (Baker et al.)"); printf("\n 3 = dsigT/dt, dsigL/dt, sigL/sigT (Bydzo)"); printf("\n"); scanf("%d", &obsq2); switch(obsq2) { case 1: /* dsigdt */ { int kin; double diff_test; char filea[80]; printf("\n\n Choose the kinematics"); printf("\n 1 : K Lambda electroproduction : Q2 dependence"); printf("\n : W = 2.17 GeV, theta = 7 deg, epsilon = 0.85"); printf("\n 2 : K Lambda electroproduction : Q2 dependance"); printf("\n : W = 2.17 GeV, theta = 1 deg, epsilon = 0.85"); printf("\n 3 : K Lambda electroproduction : Q2 dependance"); printf("\n : W = 2.17 GeV, theta = 15 deg, epsilon = 0.85"); printf("\n 4 : K Lambda electroproduction : Q2 dependence"); printf("\n : W = 2.66 GeV, theta = 8 deg, epsilon = 0.85"); printf("\n"); scanf("%d", &kin); if (kin==1) { strcpy(filea,"ep_ekl_w2p17_the7p0_eps0p85_q2"); k_th=7.; } if (kin==2) { strcpy(filea,"ep_ekl_w2p17_the1p0_eps0p85_q2"); k_th=1.; } if (kin==3) { strcpy(filea,"ep_ekl_w2p17_the15p0_eps0p85_q2"); k_th=15.; } if (kin==4) { strcpy(filea,"ep_ekl_w2p66_the8p0_eps0p85_q2"); k_th=8.; } if (kin<=3) { s = pow(2.17, 2.); } if (kin==4) { s = pow(2.66, 2.); } epsilon = 0.85; if ( (fp1 = fopen(strcat(filea,fileb),"w")) == NULL) fprintf(stderr,"\nUnable to open file"); printf("OUTPUT : Q2(GeV2), T+eps*L(at tmin)(mubarn/sr), L(at tmin)(mubarn/sr) \n"); for(Q_sqr = 0.001; Q_sqr <= 5.; Q_sqr += 0.1) { double conv; ga_mom_cm = sqrt( ( pow(s - pow(NU_M,2.) + Q_sqr, 2.) + pow(2. * NU_M, 2.) * Q_sqr ) / (4. * s) ); ga_en_cm = sqrt(s) - sqrt(pow(ga_mom_cm,2.) + pow(NU_M,2.)); k_mom_cm = sqrt( ( pow(s - pow(LAM_M,2.) - pow(K_M,2.), 2.) - pow(2. * LAM_M * K_M, 2.) ) / (4. * s) ); k_en = sqrt( pow(k_mom_cm, 2.) + pow(K_M, 2.) ); t = - Q_sqr + pow(K_M,2.) - 2. * ga_en_cm * k_en + 2. * ga_mom_cm * k_mom_cm * cos(k_th * PI/180.); conv = 2. * k_mom_cm * ga_mom_cm / (2. * PI); diff = conv * dsigma_dt_anal(reaction, mech, deg_k, deg_kstar, 0.,epsilon,Q_sqr,s,t); diffl = conv * response_R_L(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t); fprintf(fp1, "%f \t %e \t %e \n", Q_sqr, diff, diffl); } fclose(fp1); break; } case 2 : /*ratio siglsigt*/ { char file1[80]="ep_ekl_cebaf_rap.dat"; char file2[80]="ep_ekl_cebaf_seplt.dat"; double rap, conv, cos_min; int iii; /* if ( (fp1 = fopen(file1,"w")) == NULL) fprintf(stderr,"\nUnable to open file");*/ if ( (fp2 = fopen(file2,"w")) == NULL) fprintf(stderr,"\nUnable to open file"); printf("OUTPUT : t(GeV2), Q2(GeV2), T(mubarn/sr), L(mubarn/sr), L/T \n"); for(iii = 1; iii <= 6; iii += 1) { if (iii==1) { s = pow(1.84, 2.); Q_sqr = .1; /* t = -.15;*/ } if (iii==2) { s = pow(1.84, 2.); Q_sqr = .52; /* t = -.319;*/ } if (iii==3) { s = pow(1.83, 2.); Q_sqr = .75; /* t = -.42;*/ } if (iii==4) { s = pow(1.81, 2.); Q_sqr = 1.; /* t = -.557;*/ } if (iii==5) { s = pow(1.84, 2.); Q_sqr = 2.; /* t = -.961;*/ } if (iii==6) { s = pow(1.84, 2.); Q_sqr = 3.9; /* t = -1.84;*/ } ga_mom_cm = sqrt( ( pow(s - pow(NU_M,2.) + Q_sqr, 2.) + pow(2. * NU_M, 2.) * Q_sqr ) / (4. * s) ); k_mom_cm = sqrt( (pow(s - pow(LAM_M,2.) - pow(K_M,2.), 2.) - pow(2. * LAM_M * K_M,2.) ) /(4. * s)); ga_en_cm = sqrt(s) - sqrt(pow(ga_mom_cm,2.)+pow(NU_M,2.)); k_en = sqrt( pow(k_mom_cm, 2.) + pow(K_M, 2.) ); conv = 2. * k_mom_cm * ga_mom_cm / (2.*PI) ; t=pow((-Q_sqr-pow(K_M,2.)-pow(NU_M,2.)+pow(LAM_M,2.))/2./sqrt(s),2.) -pow(sqrt(pow((s-Q_sqr-pow(NU_M,2.))/2./sqrt(s),2.)+Q_sqr) -sqrt(pow((s+pow(K_M,2.)-pow(LAM_M,2.))/2./sqrt(s),2.)-pow(K_M,2.)),2.); res_trans = response_R_T(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t); res_long = response_R_L(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t); /* printf("\n %f %f %f", conv, res_trans,res_long);*/ rap = res_long / res_trans; res_trans = res_trans*conv; res_long = res_long*conv; /* fprintf(fp1, "%f \t %f \t %f \n", t, Q_sqr, rap);*/ fprintf(fp2, "%f \t %f \t %f \t %f \t %f \n", t, Q_sqr, res_trans, res_long, rap); } /* fclose(fp1);*/ fclose(fp2); break; } case 3: /* dsigdt */ { int kin; double diff_test; char filea[80]; printf("\n\n Choose the kinematics"); printf("\n 1 : K Lambda electroproduction : Q2 dependence"); printf("\n : W = 2.24 GeV, t = -.15, epsilon = 0.72"); printf("\n"); scanf("%d", &kin); if (kin==1) { strcpy(filea,"ep_ekl_w2p24_t0p15_eps0p72_q2"); t=-.15; s = pow(2.24, 2.); epsilon = 0.72; } if ( (fp1 = fopen(strcat(filea,fileb),"w")) == NULL) fprintf(stderr,"\nUnable to open file"); printf("OUTPUT : Q2(GeV2), T+eps*L(mubarn/GeV2), T(mubarn/GeV2), L(mubarn/GeV2) \n"); for(Q_sqr = 0.001; Q_sqr <= 1.; Q_sqr += 0.05) { double conv; ga_mom_cm = sqrt( ( pow(s - pow(NU_M,2.) + Q_sqr, 2.) + pow(2. * NU_M, 2.) * Q_sqr ) / (4. * s) ); ga_en_cm = sqrt(s) - sqrt(pow(ga_mom_cm,2.) + pow(NU_M,2.)); k_mom_cm = sqrt( ( pow(s - pow(LAM_M,2.) - pow(K_M,2.), 2.) - pow(2. * LAM_M * K_M, 2.) ) / (4. * s) ); k_en = sqrt( pow(k_mom_cm, 2.) + pow(K_M, 2.) ); /* conv = 2. * k_mom_cm * ga_mom_cm / (2. * PI);*/ conv = 1.; diff = conv * dsigma_dt_anal(reaction, mech, deg_k, deg_kstar, 0.,epsilon,Q_sqr,s,t); difft = conv * response_R_T(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t); diffl = conv * response_R_L(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t); fprintf(fp1, "%f \t %e \t %e \t %e \n", Q_sqr, diff, difft,diffl); } fclose(fp1); break; } } break; } case 4 : /* recoil polarization */ { double asymm_lam, asymm_lam2, tortho, lortho, res_trans, res_long, W, epsilon ; printf("\n Which W ?"); printf("\n"); scanf("%lf", &W); printf("\n Which Q2 ?"); printf("\n"); scanf("%lf", &Q_sqr); printf("\n Which epsilon ?"); printf("\n"); scanf("%lf", &epsilon); if ( (fp1 = fopen("recoillam.dat", "w")) == NULL) fprintf(stderr,"\nUnable to open file"); s=pow(W,2.); ga_mom_cm = sqrt( ( pow(s - pow(NU_M,2.) + Q_sqr, 2.) + pow(2. * NU_M, 2.) * Q_sqr ) / (4. * s) ); ga_en_cm = sqrt(s) - sqrt(pow(ga_mom_cm,2.) + pow(NU_M,2.)); k_mom_cm = sqrt( ( pow(s - pow(LAM_M,2.) - pow(K_M,2.), 2.) - pow(2. * LAM_M * K_M, 2.) ) / (4. * s) ); k_en = sqrt( pow(k_mom_cm, 2.) + pow(K_M, 2.) ); printf("OUTPUT : -t(GeV2), (Ty+eps*Ly)/(T+eps*L),ignore, Ty, Ly, T, L (mubarn/GeV2), thetaK_cm(deg) \n"); for(k_th = 0.; k_th <= 180.; k_th += 5.) { t = - Q_sqr + pow(K_M,2.) - 2. * ga_en_cm * k_en + 2. * ga_mom_cm * k_mom_cm * cos(k_th * PI / 180.); asymm_lam2 = recoil_asymm(reaction,mech,deg_k,deg_kstar, Q_sqr, s, t); tortho = response_R_T_ortho(reaction,mech,deg_k,deg_kstar, Q_sqr, s, t); res_trans = response_R_T(reaction,mech,deg_k,deg_kstar, Q_sqr, s, t); lortho = response_R_L_ortho(reaction,mech,deg_k,deg_kstar, Q_sqr, s, t); res_long = response_R_L(reaction,mech,deg_k,deg_kstar, Q_sqr, s, t); if (Q_sqr==0.) { asymm_lam = -tortho/res_trans; lortho=0.; res_long=0.; } else { asymm_lam = (tortho+epsilon*lortho) / (res_trans+epsilon*res_long); } fprintf(fp1,"%f \t %f \t %f \t %f \t %f \t %f \t %f \t %f\n", -t, asymm_lam, asymm_lam2, tortho, lortho, res_trans, res_long, k_th ); } fclose(fp); break; } case 5 : /* double polarization */ { double asymm_perp, asymm_par, asymm_long, asymm_perp0, tortho, lortho, res_trans, res_long, W, epsilon, ebeam, el_hel, x_b, ga_mom_in_lab, k_out_mom_lab, k_out_th_lab, phi, den, num_perp, num_par, num_long, num_perp0, conv, asymm_perp_prime, asymm_par_prime, asymm_long_prime, asymm_perp_prime2, asymm_par_prime2, asymm_long_prime2, virtual_photon_flux, asymm_perp2, asymm_par2, asymm_long2; int iq2t, iintegr, polbeam, iidum; printf("\n As a function of : "); printf("\n 1 : t"); printf("\n 2 : W (integrated over t)"); printf("\n"); scanf("%d", &iq2t); printf("\n Which Q2 ?"); printf("\n"); scanf("%lf", &Q_sqr); printf("\n Which beam energy ?"); printf("\n"); scanf("%lf", &ebeam); printf("\n Integration over Phi_K or specific Phi_K ?"); printf("\n 1 : integration"); printf("\n 2 : specific Phi_K"); printf("\n"); scanf("%d", &iintegr); if(iintegr==2){ printf("\n Which Phi_K ?"); printf("\n"); scanf("%lf", &phi); phi = phi * PI /180.; /*pour etre en accord avec la methode 2 qui a une definition de Phi differente : */ phi=-phi; printf("\n Double asymmetries : "); printf("\n 1 : P' (only transferred)"); printf("\n 2 : P + P'(induced + transferred) -can be compared with method 2-"); printf("\n"); scanf("%d", &polbeam); } else polbeam=0; el_hel=1.; if(iq2t==1){ if ( (fp1 = fopen("double.dat", "w")) == NULL) fprintf(stderr,"\nUnable to open file"); if ( (fp2 = fopen("double_prime.dat", "w")) == NULL) fprintf(stderr,"\nUnable to open file"); printf("\n Which W ?"); printf("\n"); scanf("%lf", &W); s=pow(W,2.); printf("In the following : \n"); if(iintegr==2){ if(polbeam==2) { printf("'den=T+eps*L+eps*cos(2Phi)*TT+sqrt(2*eps(1+eps))*cos(Phi)*TL \n"); printf("'polx=sqrt(2*eps(1+eps))*sin(Phi)*TLx+eps*sin(2Phi)*TTx+sqrt(2*eps*(1-eps))*cos(Phi)*TL'x+sqrt(1-eps*eps)*T'x \n"); printf("'polz=sqrt(2*eps(1+eps))*sin(Phi)*TLz+eps*sin(2Phi)*TTz+sqrt(2*eps*(1-eps))*cos(Phi)*TL'z+sqrt(1-eps*eps)*T'z \n"); printf("'poly=sqrt(2*eps* (1-eps))*sin(Phi)*TL'y \n"); /* printf("'pol0= Ty+eps*Ly+sqrt(2*eps*(1+eps))*cos(Phi)*TLy+eps*cos(2*Phi)*TTy \n");*/ } if(polbeam==1) { printf("'den=T+eps*L+eps*cos(2Phi)*TT+sqrt(2*eps(1+eps))*cos(Phi)*TL \n"); printf("'polx=sqrt(2*eps*(1-eps))*cos(Phi)*TL'x+sqrt(1-eps*eps)*T'x \n"); printf("'polz=sqrt(2*eps*(1-eps))*cos(Phi)*TL'z+sqrt(1-eps*eps)*T'z \n"); printf("'poly=sqrt(2*eps* (1-eps))*sin(Phi)*TL'y \n"); /* printf("'pol0= Ty+eps*Ly+sqrt(2*eps*(1+eps))*cos(Phi)*TLy+eps*cos(2*Phi)*TTy \n");*/ } } if(iintegr==1) { printf("'den'=T+eps*L \n"); printf("'polx'=sqrt(1-eps*eps)*T'x \n"); printf("'polz'=sqrt(1-eps*eps)*T'z \n"); /*printf("'poly'=Ty+eps*Ly \n");*/ printf("'poly'=0. \n"); printf("'pol0'= Ty+eps*Ly \n"); } x_b = Q_sqr / (s + Q_sqr - pow(NU_M,2.)); k_out_mom_lab = ebeam - Q_sqr / (2. * NU_M * x_b); k_out_th_lab = 2. * asin(sqrt(Q_sqr / (4. * ebeam * k_out_mom_lab))); ga_mom_in_lab = sqrt(pow(ebeam, 2.) + pow(k_out_mom_lab, 2.) - 2. * ebeam * k_out_mom_lab * cos(k_out_th_lab) ); epsilon = 1. / (1. + 2. * pow(ga_mom_in_lab, 2.) / Q_sqr * pow(tan(0.5 * k_out_th_lab), 2.) ); printf("hello %lf ",epsilon); ga_mom_cm = sqrt( ( pow(s - pow(NU_M,2.) + Q_sqr, 2.) + pow(2. * NU_M, 2.) * Q_sqr ) / (4. * s) ); ga_en_cm = sqrt(s) - sqrt(pow(ga_mom_cm,2.) + pow(NU_M,2.)); k_mom_cm = sqrt( ( pow(s - pow(LAM_M,2.) - pow(K_M,2.), 2.) - pow(2. * LAM_M * K_M, 2.) ) / (4. * s) ); k_en = sqrt( pow(k_mom_cm, 2.) + pow(K_M, 2.) ); printf("OUTPUT : cos(thetaK_cm),-t(GeV2),poly/den,polx/den,polz/den\n"); for(k_th = 0.; k_th <= 180.; k_th += 10.) { t = - Q_sqr + pow(K_M,2.) - 2. * ga_en_cm * k_en + 2. * ga_mom_cm * k_mom_cm * cos(k_th * PI / 180.); /*Dans la derivation des fonctions de reponse en terme de courants J, il y a un signe moins problematique pour toutes les fonctions de reponse de la forme TL ; plutot que de changer dans la definition des fonctions de reponse, on met un signe moins devant chaque TL dans ce qui suit SAUF pour R_TL non-pol pour lequel le changement a ete fait au niveau de la definition */ if(iintegr==2){ den = response_R_T(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t) + epsilon * response_R_L(reaction,mech, deg_k, deg_kstar, Q_sqr, s, t) - epsilon * cos(2. * phi) * response_R_TT(reaction,mech, deg_k, deg_kstar, Q_sqr, s, t) + sqrt(2. * epsilon * (1. + epsilon)) * cos(phi) * response_R_TL(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t); if(polbeam==2) num_par = -sqrt(2. * epsilon * (1. + epsilon)) * sin(phi) * response_R_TL_par(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t) - epsilon * sin(2. * phi) * response_R_TT_par(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t) + el_hel * ( -sqrt(2. * epsilon * (1. - epsilon)) * cos(phi) * response_R_TLp_par(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t) + sqrt(1.-pow(epsilon,2.)) * response_R_Tp_par(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t)); if(polbeam==1) num_par = el_hel * ( -sqrt(2. * epsilon * (1. - epsilon)) * cos(phi) * response_R_TLp_par(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t) + sqrt(1.-pow(epsilon,2.)) * response_R_Tp_par(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t)); if(polbeam==2) num_long = -sqrt(2. * epsilon * (1. + epsilon)) * sin(phi) * response_R_TL_long(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t) - epsilon * sin(2. * phi) * response_R_TT_long(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t) + el_hel * ( -sqrt(2. * epsilon * (1. - epsilon)) * cos(phi) * response_R_TLp_long(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t) + sqrt(1.-pow(epsilon,2.)) * response_R_Tp_long(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t)); if(polbeam==1) num_long = el_hel * ( -sqrt(2. * epsilon * (1. - epsilon)) * cos(phi) * response_R_TLp_long(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t) + sqrt(1.-pow(epsilon,2.)) * response_R_Tp_long(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t)); if(polbeam==2) num_perp = response_R_T_ortho(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t) +epsilon * response_R_L_ortho(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t) -sqrt(2. * epsilon * (1. + epsilon)) * cos(phi) * response_R_TL_ortho(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t) - epsilon * cos(2. * phi) * response_R_TT_ortho(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t) + el_hel * (-sqrt(2. * epsilon * (1. - epsilon))) * sin(phi) * response_R_TLp_ortho(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t); if(polbeam==1) num_perp = el_hel *( -sqrt(2. * epsilon * (1. - epsilon)) * sin(phi) * response_R_TLp_ortho(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t)); num_perp0 = response_R_T_ortho(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t) +epsilon * response_R_L_ortho(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t) -sqrt(2. * epsilon * (1. + epsilon)) * cos(phi) * response_R_TL_ortho(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t) - epsilon * cos(2. * phi) * response_R_TT_ortho(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t); } if(iintegr==1){ den = response_R_T(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t) + epsilon * response_R_L(reaction,mech, deg_k, deg_kstar, Q_sqr, s, t); num_par = el_hel * sqrt(1.-pow(epsilon,2.)) * response_R_Tp_par(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t); num_long = el_hel * sqrt(1.-pow(epsilon,2.)) * response_R_Tp_long(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t); num_perp = 0.; num_perp0 = response_R_T_ortho(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t) +epsilon * response_R_L_ortho(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t); phi=0.; printf("\n TOT,TTs,TTl %lf %lf %lf %lf %lf", den,response_R_Tp_par(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t), response_R_Tp_long(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t),num_par / den,num_long / den); } /* printf("\n T,L,TT,TL %lf %lf %lf %lf %lf \n T,L,TT,TLnorm %lf %lf %lf %lf \n TT,TLs %lf %lf TT,TLl %lf %lf", k_th,response_R_T(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t), response_R_TT(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t), response_R_L(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t), response_R_TL(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t)*ga_en_cm/sqrt(Q_sqr)*2., response_R_T_ortho(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t), response_R_L_ortho(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t), response_R_TT_ortho(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t), response_R_TL_ortho(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t)*ga_en_cm/sqrt(Q_sqr)*2., response_R_Tp_par(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t), response_R_TLp_par(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t)*ga_en_cm/sqrt(Q_sqr)*2., response_R_Tp_long(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t), response_R_TLp_long(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t)*ga_en_cm/sqrt(Q_sqr)*2.); */ asymm_perp = num_perp / den; asymm_par = num_par / den; asymm_long = num_long / den; asymm_perp0 = num_perp0 / den; /* Faire une rotation en Phi car, dans cette methode, le repere est attache au plan hadronique : il faut repasser dans le plan leptonique */ asymm_par2 = asymm_par*cos(phi)-asymm_perp*sin(phi); asymm_perp2 = asymm_par*sin(phi)+asymm_perp*cos(phi); asymm_long2 = asymm_long; /* signe - en dessous pour la composante perp car dans le formalisme des fonctions de reponse, le Phi a ete inverse tout le long (phi->-phi)*/ asymm_perp = -asymm_perp2; asymm_par = asymm_par2; asymm_long = asymm_long2; /*On repasse dans le plan hadronique...*/ if(iintegr==2) { asymm_par_prime = asymm_par*cos(-phi)-asymm_perp*sin(-phi); asymm_perp_prime = asymm_par*sin(-phi)+asymm_perp*cos(-phi); asymm_long_prime = asymm_long; } /*S'il y a integration, phi avait ete defini egal a 0 quand on etait passe dans le plan leptonique, la rotation est triviale...*/ if(iintegr==1) { asymm_par_prime = asymm_par; asymm_perp_prime = asymm_perp; asymm_long_prime = asymm_long; } /* si on veut definir l'axe z selon la direction du Lambda de recul (utile pour comparer la composante "long" avec la methode 2 et les spineurs d'helicite) k_th=k_th+180.; */ asymm_par_prime2 = asymm_par_prime*cos(k_th* PI / 180.)-asymm_long_prime*sin(k_th* PI / 180.); asymm_perp_prime2 = asymm_perp_prime; asymm_long_prime2 = asymm_par_prime*sin(k_th* PI / 180.)+asymm_long_prime*cos(k_th* PI / 180.); /* k_th=k_th-180.; */ virtual_photon_flux = pow(ELEC,2.)/(4. * PI) /(2. * pow(PI,2.)) / (2. * pow(ebeam,2.)) * (s-pow(NU_M,2.))/(2.*NU_M) /Q_sqr /(1. - epsilon); fprintf(fp1,"%f \t %f \t %f \t %f \t %f \t %f \n", cos(k_th* PI / 180.), -t, asymm_perp, asymm_par, asymm_long, den); fprintf(fp2,"%f \t %f \t %f \t %f \t %f \n", cos(k_th* PI / 180.), -t, asymm_perp_prime2, asymm_par_prime2, asymm_long_prime2,den); printf("%f \t %f \t %f \t %f \t %f \n", cos(k_th* PI / 180.), -t, asymm_perp, asymm_par, asymm_long); } } if(iq2t==2){ printf("\n Which frame ?"); printf("\n 1 : xyz ?"); printf("\n 2 : x'y'z' ?"); printf("\n"); scanf("%d", &rep_s); if(rep_s==1) { if ( (fp1 = fopen("double.dat", "w")) == NULL) fprintf(stderr,"\nUnable to open file"); } if(rep_s==2) { if ( (fp2 = fopen("double_prime.dat", "w")) == NULL) fprintf(stderr,"\nUnable to open file"); } printf("In the following : \n"); if(iintegr==2){ if(polbeam==2) { printf("'den=T+eps*L+eps*cos(2Phi)*TT+sqrt(2*eps(1+eps))*cos(Phi)*TL \n"); printf("'polx=sqrt(2*eps(1+eps))*sin(Phi)*TLx+eps*sin(2Phi)*TTx+sqrt(2*eps*(1-eps))*cos(Phi)*TL'x+sqrt(1-eps*eps)*T'x \n"); printf("'polz=sqrt(2*eps(1+eps))*sin(Phi)*TLz+eps*sin(2Phi)*TTz+sqrt(2*eps*(1-eps))*cos(Phi)*TL'z+sqrt(1-eps*eps)*T'z \n"); printf("'poly=sqrt(2*eps* (1-eps))*sin(Phi)*TL'y \n"); /* printf("'pol0= Ty+eps*Ly+sqrt(2*eps*(1+eps))*cos(Phi)*TLy+eps*cos(2*Phi)*TTy \n");*/ } if(polbeam==1) { printf("'den=T+eps*L+eps*cos(2Phi)*TT+sqrt(2*eps(1+eps))*cos(Phi)*TL \n"); printf("'polx=sqrt(2*eps*(1-eps))*cos(Phi)*TL'x+sqrt(1-eps*eps)*T'x \n"); printf("'polz=sqrt(2*eps*(1-eps))*cos(Phi)*TL'z+sqrt(1-eps*eps)*T'z \n"); printf("'poly=sqrt(2*eps* (1-eps))*sin(Phi)*TL'y \n"); /* printf("'pol0= Ty+eps*Ly+sqrt(2*eps*(1+eps))*cos(Phi)*TLy+eps*cos(2*Phi)*TTy \n");*/ } } if(iintegr==1) { printf("'den'=T+eps*L \n"); printf("'polx'=sqrt(1-eps*eps)*T'x \n"); printf("'polz'=sqrt(1-eps*eps)*T'z \n"); /*printf("'poly'=Ty+eps*Ly \n");*/ printf("'poly'=0. \n"); printf("'pol0'= Ty+eps*Ly \n"); } printf("OUTPUT : W(GeV),poly/den,polx/den,polz/den\n"); for(W = 1.6; W <= 2.3; W += .1) { s=pow(W,2.); x_b = Q_sqr / (s + Q_sqr - pow(NU_M,2.)); k_out_mom_lab = ebeam - Q_sqr / (2. * NU_M * x_b); k_out_th_lab = 2. * asin(sqrt(Q_sqr / (4. * ebeam * k_out_mom_lab))); ga_mom_in_lab = sqrt(pow(ebeam, 2.) + pow(k_out_mom_lab, 2.) - 2. * ebeam * k_out_mom_lab * cos(k_out_th_lab) ); epsilon = 1. / (1. + 2. * pow(ga_mom_in_lab, 2.) / Q_sqr * pow(tan(0.5 * k_out_th_lab), 2.) ); iidum=1; asymm_par = pol_tot(reaction, mech, deg_k, deg_kstar, Q_sqr, s, epsilon, phi, iintegr, polbeam, iidum); iidum=2; asymm_perp = pol_tot(reaction, mech, deg_k, deg_kstar, Q_sqr, s, epsilon, phi, iintegr, polbeam, iidum); iidum=3; asymm_long = pol_tot(reaction, mech, deg_k, deg_kstar, Q_sqr, s, epsilon, phi, iintegr, polbeam, iidum); if(rep_s==1) fprintf(fp1,"%f \t %f \t %f \t %f \t %f \n", W, asymm_perp, asymm_par, asymm_long, W); if(rep_s==2) fprintf(fp2,"%f \t %f \t %f \t %f \t %f \n", W, asymm_perp, asymm_par, asymm_long, W); printf("%f \t %f \t %f \t %f \n", W, asymm_perp, asymm_par, asymm_long); } } fclose(fp1); fclose(fp2); break; } case 6 : /* double polarization : method 2 */ { double k_in_mom_lab, k_out_mom_lab, k_out_th_lab, kaon_th, ga_en_in, ga_en_in_lab, ga_mom_in_lab, ga_mom_in, mes_mom_lab, Q_sqr, s, t, W, phi, el_hel, result1, result2, result3, result4, result5, result6, result7, result8, result9, result10, asymm_x, asymm_y, asymm_z, asymm_0, cross, k_en, k_mom_cm, asymm_x_prime, asymm_y_prime, asymm_z_prime, asymm_x_prime2, asymm_y_prime2, asymm_z_prime2 ; double result15,result16,result17,result18,result47,result48,result57,result58; int iintegr, isys2; printf("\n\n Give the value of beam energy in GeV (e.g. 6.)"); printf("\n"); scanf("%lf", &k_in_mom_lab); printf("\n\n Give the value of Q^2 in GeV^2 (e.g. 2.0)"); printf("\n"); scanf("%lf", &Q_sqr); printf("\n\n Give the value of W (e.g. 2.)"); printf("\n"); scanf("%lf", &W); s=pow(W,2.); /* printf("\n Integration over Phi_K or specific Phi_K ?"); printf("\n 1 : integration"); printf("\n 2 : specific Phi_K"); printf("\n"); scanf("%d", &iintegr); */ iintegr=2; if(iintegr==2){ printf("\n Which Phi_K ?"); printf("\n"); scanf("%lf", &phi); phi = phi * PI /180.; } /* printf("\n Which reference system for polarization ?"); printf("\n 1 : axis z along incoming virtual photon"); printf("\n 2 : axis z along recoil nucleon"); printf("\n"); scanf("%d", &isys2); */ isys2=1; /* printf("\n Which nucleon polarization ?"); printf("\n 1 : target"); printf("\n 2 : recoil"); printf("\n"); scanf("%d", &targrec); */ targrec=2; ga_en_in_lab= (s - pow(NU_M, 2.) + Q_sqr)/ 2. / NU_M; k_out_mom_lab = k_in_mom_lab - ga_en_in_lab; k_out_th_lab = 2. * asin(sqrt(Q_sqr / (4. * k_in_mom_lab * k_out_mom_lab))); ga_mom_cm = sqrt( ( pow(s - pow(NU_M,2.) + Q_sqr, 2.) + pow(2. * NU_M, 2.) * Q_sqr ) / (4. * s) ); ga_en_cm = sqrt(s) - sqrt(pow(ga_mom_cm,2.) + pow(NU_M,2.)); k_mom_cm = sqrt( ( pow(s - pow(LAM_M,2.) - pow(K_M,2.), 2.) - pow(2. * LAM_M * K_M, 2.) ) / (4. * s) ); k_en = sqrt( pow(k_mom_cm, 2.) + pow(K_M, 2.) ); if ( (fp1 = fopen("double_phi.dat", "w")) == NULL) fprintf(stderr,"\nUnable to open file"); if ( (fp2 = fopen("double_phi_prime.dat", "w")) == NULL) fprintf(stderr,"\nUnable to open file"); if (iintegr==2) { for( kaon_th = 0.; kaon_th <= 180.; kaon_th += 10.) { t = -Q_sqr + pow(K_M,2.) - 2. * ga_en_cm * k_en + 2. * ga_mom_cm * k_mom_cm * cos(kaon_th * PI /180.); isys=isys2; /* nucleon polarization along z-axis : isys=1 or 2 */ result1 = klam_hel_spin_diff_cross_lab(mech,deg_k,deg_kstar,1.,+.5, k_in_mom_lab, k_out_mom_lab, k_out_th_lab, kaon_th * PI/180., phi); result2 = klam_hel_spin_diff_cross_lab(mech,deg_k,deg_kstar,1.,-.5, k_in_mom_lab, k_out_mom_lab, k_out_th_lab, kaon_th * PI/180., phi); /* needed for nucleon polarization along z-axis but without electron polarization */ result3 = klam_hel_spin_diff_cross_lab(mech,deg_k,deg_kstar,-1.,+.5, k_in_mom_lab, k_out_mom_lab, k_out_th_lab, kaon_th * PI/180., phi); result4 = klam_hel_spin_diff_cross_lab(mech,deg_k,deg_kstar,-1.,-.5, k_in_mom_lab, k_out_mom_lab, k_out_th_lab, kaon_th * PI/180., phi); isys=3; /* nucleon polarization along x-axis : isys=3 */ result5 = klam_hel_spin_diff_cross_lab(mech,deg_k,deg_kstar,1.,+.5, k_in_mom_lab, k_out_mom_lab, k_out_th_lab, kaon_th * PI/180., phi); result6 = klam_hel_spin_diff_cross_lab(mech,deg_k,deg_kstar,1.,-.5, k_in_mom_lab, k_out_mom_lab, k_out_th_lab, kaon_th * PI/180., phi); /* result15 = klam_hel_spin_diff_cross_lab(mech,deg_k,deg_kstar,-1.,+.5, k_in_mom_lab, k_out_mom_lab, k_out_th_lab, kaon_th * PI/180., phi); result16 = klam_hel_spin_diff_cross_lab(mech,deg_k,deg_kstar,-1.,-.5, k_in_mom_lab, k_out_mom_lab, k_out_th_lab, kaon_th * PI/180., phi); */ isys=4; /* nucleon polarization along y-axis : isys=4 */ result7 = klam_hel_spin_diff_cross_lab(mech,deg_k,deg_kstar,1.,+.5, k_in_mom_lab, k_out_mom_lab, k_out_th_lab, kaon_th * PI/180., phi); result8 = klam_hel_spin_diff_cross_lab(mech,deg_k,deg_kstar,1.,-.5, k_in_mom_lab, k_out_mom_lab, k_out_th_lab, kaon_th * PI/180., phi); result17 = klam_hel_spin_diff_cross_lab(mech,deg_k,deg_kstar,-1.,+.5, k_in_mom_lab, k_out_mom_lab, k_out_th_lab, kaon_th * PI/180., phi); result18 = klam_hel_spin_diff_cross_lab(mech,deg_k,deg_kstar,-1.,-.5, k_in_mom_lab, k_out_mom_lab, k_out_th_lab, kaon_th * PI/180., phi); /* result47 = klam_hel_spin_diff_cross_lab(mech,deg_k,deg_kstar,1.,+.5, k_in_mom_lab, k_out_mom_lab, k_out_th_lab, kaon_th * PI/180., phi+PI); result48 = klam_hel_spin_diff_cross_lab(mech,deg_k,deg_kstar,1.,-.5, k_in_mom_lab, k_out_mom_lab, k_out_th_lab, kaon_th * PI/180., phi+PI); result57 = klam_hel_spin_diff_cross_lab(mech,deg_k,deg_kstar,-1.,+.5, k_in_mom_lab, k_out_mom_lab, k_out_th_lab, kaon_th * PI/180., phi+PI); result58 = klam_hel_spin_diff_cross_lab(mech,deg_k,deg_kstar,-1.,-.5, k_in_mom_lab, k_out_mom_lab, k_out_th_lab, kaon_th * PI/180., phi+PI); */ asymm_z = (result1 - result2)/(result1 + result2); asymm_x = (result5 - result6)/(result5 + result6); asymm_y = (result7 - result8)/(result7 + result8); asymm_0 = (result7 + result17 - result8 -result18)/ (result7 + result17 + result8 + result18); /* asymm_y = asymm_y - asymm_0;*/ asymm_x_prime = asymm_x*cos(phi)-asymm_y*sin(phi); asymm_y_prime = asymm_x*sin(phi)+asymm_y*cos(phi); asymm_z_prime = asymm_z; asymm_x_prime2 = asymm_x_prime*cos(kaon_th* PI / 180.)-asymm_z_prime*sin(kaon_th* PI / 180.); asymm_y_prime2 = asymm_y_prime; asymm_z_prime2 = asymm_x_prime*sin(kaon_th* PI / 180.)+asymm_z_prime*cos(kaon_th* PI / 180.); fprintf(fp1,"\n %f \t %e \t %e \t %e \t %e \t %e \t %e", cos(kaon_th* PI / 180.), -t, asymm_y, asymm_x, asymm_z, asymm_0, .5*(result1+result2+result3+result4)); fprintf(fp2,"\n %f \t %e \t %e \t %e \t %e", cos(kaon_th* PI / 180.), -t, asymm_y_prime2, asymm_x_prime2, asymm_z_prime2); printf("%f \t %f \t %f \t %f \t %f \n", cos(kaon_th* PI / 180.), -t, asymm_y, asymm_x, asymm_z); } } fclose(fp1); break; } case 7 : /* User-defined */ { double a, a_min, a_max, a_step, W_min, W_max, W_step, t_min, t_max, t_step, Q2_min, Q2_max, Q2_step, bidon1, bidon2, cos_th_cm, cth_min, cth_max, cth_step; double ga_en_in_lab, ga_mom_in_lab, pi_mom_cm, W, tmin_test, mes_mom_lab_test, mes_en_lab, mes_mom_lab; double x_b,k_out_mom_lab,k_out_th_lab,ebeam,s,epsilon; double res_u; int rep, rep2; char file1t[80],file1l[80],file1tt[80],file1tl[80],file1tlp[80],file1u[80]; printf("\n 1 : Q2 dependence"); printf("\n 2 : t dependence"); printf("\n 3 : W dependence"); printf("\n 4 : cthkcm dependence"); printf("\n"); scanf("%d", &rep); if (rep==1) { printf("\n Which W"); printf("\n"); scanf("%lf", &W); s = pow(W,2.); printf("\n Which cos_th_cm ?"); printf("\n"); scanf("%lf", &cos_th_cm); /*printf("\n Q2_min ?"); printf("\n"); scanf("%lf", &Q2_min); printf("\n Which Q2_max ?"); printf("\n"); scanf("%lf", &Q2_max); printf("\n Which Q2_step ?"); printf("\n"); scanf("%lf", &Q2_step);*/ Q2_min=0.2; Q2_max=5.0; Q2_step=0.1; a_min=Q2_min; a_max=Q2_max; a_step=Q2_step; bidon1=(W-(int) W)*10.; bidon2=(t-(int) t)*10.; /*sprintf(file1t,"ep_klambda_t%dp%d_w%dp%d_trans", (int) t, (int) bidon2, (int) W, (int) bidon1); sprintf(file1l,"ep_klambda_t%dp%d_w%dp%d_long", (int) t, (int) bidon2, (int) W, (int) bidon1); sprintf(file1tt,"ep_klambda_t%dp%d_w%dp%d_tt", (int) t, (int) bidon2, (int) W, (int) bidon1); sprintf(file1tl,"ep_klambda_t%dp%d_w%dp%d_tl", (int) t, (int) bidon2, (int) W, (int) bidon1); printf("OUTPUT1 : Q2 (GeV2), dsig_T/dt (mubarn/GeV2), dsig_T/dOmega_cm (mubarn/sr) \n"); printf("OUTPUT2 : Q2 (GeV2), dsig_L/dt (mubarn/GeV2), dsig_L/dOmega_cm (mubarn/sr) \n"); printf("OUTPUT3 : Q2 (GeV2), dsig_TT/dt (mubarn/GeV2), dsig_TT/dOmega_cm (mubarn/sr) \n"); printf("OUTPUT4 : Q2 (GeV2), dsig_TL/dt (mubarn/GeV2), dsig_TL/dOmega_cm (mubarn/sr) \n");*/ sprintf(file1u,"ep_klambda_q2", (int) Q_sqr, (int) bidon1, (int) t, (int) bidon2); printf("OUTPUT : Q2 (GeV2), sigu, siglt, sigtt, sigltp (nb/sr) \n"); } if (rep==2) { printf("\n Which W"); printf("\n"); scanf("%lf", &W); printf("\n Which Q2"); printf("\n"); scanf("%lf", &Q_sqr); s = pow(W,2.); printf("\n t_min ? (if 0 => kinematical t_min)"); printf("\n"); scanf("%lf", &t_min); printf("\n t_max ?"); printf("\n"); scanf("%lf", &t_max); printf("\n t_step ?"); printf("\n"); scanf("%lf", &t_step); a_min=t_min; a_max=t_max; a_step=t_step; if(t_min==0.) { a_min=pow((-Q_sqr-pow(K_M,2.)-pow(NU_M,2.)+pow(LAM_M,2.))/2./sqrt(s),2.) -pow(sqrt(pow((s-Q_sqr-pow(NU_M,2.))/2./sqrt(s),2.)+Q_sqr) -sqrt(pow((s+pow(K_M,2.)-pow(LAM_M,2.))/2./sqrt(s),2.)-pow(K_M,2.)),2.); a_min=-a_min; } bidon1=(Q_sqr-(int) Q_sqr)*10.; bidon2=(W-(int) W)*10.; sprintf(file1t,"ep_klambda_q2%dp%d_w%dp%d_trans", (int) Q_sqr, (int) bidon1, (int) W, (int) bidon2); sprintf(file1l,"ep_klambda_q2%dp%d_w%dp%d_long", (int) Q_sqr, (int) bidon1, (int) W, (int) bidon2); sprintf(file1tt,"ep_klambda_q2%dp%d_w%dp%d_tt", (int) Q_sqr, (int) bidon1, (int) W, (int) bidon2); sprintf(file1tl,"ep_klambda_q2%dp%d_w%dp%d_tl", (int) Q_sqr, (int) bidon1, (int) W, (int) bidon2); printf("OUTPUT1 : -t (GeV2), dsig_T/dt (mubarn/GeV2), dsig_T/dOmega_cm (mubarn/sr) \n"); printf("OUTPUT2 : -t (GeV2), dsig_L/dt (mubarn/GeV2), dsig_L/dOmega_cm (mubarn/sr) \n"); printf("OUTPUT3 : -t (GeV2), dsig_TT/dt (mubarn/GeV2), dsig_TT/dOmega_cm (mubarn/sr) \n"); printf("OUTPUT4 : -t (GeV2), dsig_TL/dt, (mubarn/GeV2),dsig_TL/dOmega_cm (mubarn/sr) \n"); } if (rep==3) { printf("\n Which Q2"); printf("\n"); scanf("%lf", &Q_sqr); printf("\n Fixed t or fixed cos_th_cm"); printf("\n 1 : Fixed t"); printf("\n 2 : Fixed cos_th_cm"); printf("\n"); scanf("%d", &rep2); if (rep2==1) { printf("\n Which t (t>0) -if t=0, t=t_min-"); printf("\n"); scanf("%lf", &t); mt=-t; } if (rep2==2) { printf("\n Which cos_th_cm ?"); printf("\n"); scanf("%lf", &cos_th_cm); } /*printf("\n W_min ?"); printf("\n"); scanf("%lf", &W_min); printf("\n W_max ?"); printf("\n"); scanf("%lf", &W_max); printf("\n W_step ?"); printf("\n"); scanf("%lf", &W_step);*/ W_min=1.61; W_max=2.80; W_step=0.005; a_min=W_min; a_max=W_max; a_step=W_step; bidon1=(Q_sqr-(int) Q_sqr)*10.; bidon2=(t-(int) t)*10.; /*sprintf(file1t,"ep_klambda_q2%dp%d_t%dp%d_trans", (int) Q_sqr, (int) bidon1, (int) t, (int) bidon2); sprintf(file1l,"ep_klambda_q2%dp%d_t%dp%d_long", (int) Q_sqr, (int) bidon1, (int) t, (int) bidon2); sprintf(file1tt,"ep_klambda_q2%dp%d_t%dp%d_tt", (int) Q_sqr, (int) bidon1, (int) t, (int) bidon2); sprintf(file1tl,"ep_klambda_q2%dp%d_t%dp%d_tl", (int) Q_sqr, (int) bidon1, (int) t, (int) bidon2); sprintf(file1tlp,"ep_klambda_q2%dp%d_t%dp%d_tlp", (int) Q_sqr, (int) bidon1, (int) t, (int) bidon2); printf("OUTPUT1 : W (GeV), dsig_T/dt (nb/GeV2), dsig_T/dOmega_cm (nb/sr) \n"); printf("OUTPUT2 : W (GeV), dsig_L/dt (nb/GeV2), dsig_L/dOmega_cm (nb/sr) \n"); printf("OUTPUT3 : W (GeV), dsig_TT/dt (nb/GeV2), dsig_TT/dOmega_cm (nb/sr) \n"); printf("OUTPUT4 : W (GeV), dsig_TL/dt (nb/GeV2), dsig_TL/dOmega_cm (nb/sr) \n"); printf("OUTPUT5 : W (GeV), dsig_TLp/dt (nb/GeV2), dsig_TLp/dOmega_cm (nb/sr) \n");*/ sprintf(file1u,"ep_klambda_w", (int) Q_sqr, (int) bidon1, (int) t, (int) bidon2); printf("OUTPUT : W (GeV), sigu, siglt, sigtt, sigltp (nb/sr) \n"); } if (rep==4) { printf("\n Which Q2"); printf("\n"); scanf("%lf", &Q_sqr); printf("\n Which W"); printf("\n"); scanf("%lf", &W); s = pow(W,2.); /*printf("\n cth_min ?"); printf("\n"); scanf("%lf", &cth_min); printf("\n cth_max ?"); printf("\n"); scanf("%lf", &cth_max); printf("\n cth_step ?"); printf("\n"); scanf("%lf", &cth_step);*/ cth_min=-1.0; cth_max=1.0; cth_step=0.02; a_min=cth_min; a_max=cth_max; a_step=cth_step; bidon1=(Q_sqr-(int) Q_sqr)*10.; bidon2=(W-(int) W)*10.; /*sprintf(file1t,"ep_klambda_q2%dp%d_t%dp%d_trans", (int) Q_sqr, (int) bidon1, (int) W, (int) bidon2); sprintf(file1l,"ep_klambda_q2%dp%d_t%dp%d_long", (int) Q_sqr, (int) bidon1, (int) W, (int) bidon2); sprintf(file1tt,"ep_klambda_q2%dp%d_t%dp%d_tt", (int) Q_sqr, (int) bidon1, (int) W, (int) bidon2); sprintf(file1tl,"ep_klambda_q2%dp%d_t%dp%d_tl", (int) Q_sqr, (int) bidon1, (int) W, (int) bidon2); sprintf(file1tlp,"ep_klambda_q2%dp%d_t%dp%d_tlp", (int) Q_sqr, (int) bidon1, (int) W, (int) bidon2); printf("OUTPUT1 : W (GeV), dsig_T/dt (nb/GeV2), dsig_T/dOmega_cm (nb/sr) \n"); printf("OUTPUT2 : W (GeV), dsig_L/dt (nb/GeV2), dsig_L/dOmega_cm (nb/sr) \n"); printf("OUTPUT3 : W (GeV), dsig_TT/dt (nb/GeV2), dsig_TT/dOmega_cm (nb/sr) \n"); printf("OUTPUT4 : W (GeV), dsig_TL/dt (nb/GeV2), dsig_TL/dOmega_cm (nb/sr) \n"); printf("OUTPUT5 : W (GeV), dsig_TLp/dt (nb/GeV2), dsig_TLp/dOmega_cm (nb/sr) \n");*/ sprintf(file1u,"ep_klambda_a", (int) Q_sqr, (int) bidon1, (int) t, (int) bidon2); printf("OUTPUT : cthkcm, sigu, siglt, sigtt, sigltp (nb/sr) \n"); } /*if ( (fp1 = fopen(strcat(file1t,fileb),"w")) == NULL) fprintf(stderr,"\nUnable to open file"); if ( (fp2 = fopen(strcat(file1l,fileb),"w")) == NULL) fprintf(stderr,"\nUnable to open file"); if ( (fp3 = fopen(strcat(file1tt,fileb),"w")) == NULL) fprintf(stderr,"\nUnable to open file"); if ( (fp4 = fopen(strcat(file1tl,fileb),"w")) == NULL) fprintf(stderr,"\nUnable to open file"); if ( (fp5 = fopen(strcat(file1tlp,fileb),"w")) == NULL) fprintf(stderr,"\nUnable to open file");*/ if ( (fp1 = fopen(strcat(file1u,fileb),"w")) == NULL) fprintf(stderr,"\nUnable to open file"); for(a = a_min; a <= a_max; a += a_step) { double conv, cos_k, ga_en_cm2; if (rep==1) Q_sqr=a; if (rep==2) mt=-a; if (rep==3) s=pow(a,2.); if (rep==4) cos_th_cm=a; if (rep==1) { ga_en_cm = sqrt(s) - sqrt(pow(ga_mom_cm,2.) + pow(NU_M,2.)); k_en = sqrt( pow(pi_mom_cm, 2.) + pow(K_M, 2.) ); mt = -Q_sqr + pow(K_M,2.) - 2. * ga_en_cm * k_en + 2. * ga_mom_cm * pi_mom_cm * cos_th_cm; } if((t==0)&&(rep!=2)) { mt=pow((-Q_sqr-pow(K_M,2.)-pow(NU_M,2.)+pow(LAM_M,2.))/2./sqrt(s),2.) -pow(sqrt(pow((s-Q_sqr-pow(NU_M,2.))/2./sqrt(s),2.)+Q_sqr) -sqrt(pow((s+pow(K_M,2.)-pow(LAM_M,2.))/2./sqrt(s),2.)-pow(K_M,2.)),2.); } ga_mom_cm = sqrt( ( pow(s - pow(NU_M,2.) + Q_sqr, 2.) + pow(2. * NU_M, 2.) * Q_sqr ) / (4. * s) ); pi_mom_cm = sqrt( ( pow(s - pow(LAM_M,2.) - pow(K_M,2.), 2.) - pow(2. * LAM_M * K_M, 2.) ) / (4. * s) ); if (rep2==2) { ga_en_cm = sqrt(s) - sqrt(pow(ga_mom_cm,2.) + pow(NU_M,2.)); k_en = sqrt( pow(pi_mom_cm, 2.) + pow(K_M, 2.) ); mt = -Q_sqr + pow(K_M,2.) - 2. * ga_en_cm * k_en + 2. * ga_mom_cm * pi_mom_cm * cos_th_cm; } if (rep==4) { ga_en_cm = sqrt(s) - sqrt(pow(ga_mom_cm,2.) + pow(NU_M,2.)); k_en = sqrt( pow(pi_mom_cm, 2.) + pow(K_M, 2.) ); mt = -Q_sqr + pow(K_M,2.) - 2. * ga_en_cm * k_en + 2. * ga_mom_cm * pi_mom_cm * cos_th_cm; } conv = 1. / (2. * PI) * 2. * pi_mom_cm * ga_mom_cm; res_trans = response_R_T(reaction, mech, deg_k, deg_kstar, Q_sqr, s, mt); res_tt = response_R_TT(reaction, mech, deg_k, deg_kstar, Q_sqr, s, mt); if (Q_sqr==0.) { res_long=0.; res_tl=0.; res_tlp=0.; } else { res_long = response_R_L(reaction, mech, deg_k, deg_kstar, Q_sqr, s, mt); res_tl = response_R_TL(reaction, mech, deg_k, deg_kstar, Q_sqr, s, mt); res_tlp = response_R_TLp(reaction, mech, deg_k, deg_kstar, Q_sqr, s, mt); } /*in order to compare to Feuerbach's results who has a different convention*/ /* ga_en_in_lab= (s - pow(NU_M, 2.) + Q_sqr)/ 2. / NU_M; res_trans=res_trans; res_long=res_long/Q_sqr*pow(ga_en_in_lab,2.); res_tt=res_tt; res_tl=res_tl/sqrt(Q_sqr)*ga_en_in_lab; */ /* Changes based on e1f paper conventions. Also use nb/sr for output */ res_trans = res_trans*1000.; res_long = res_long*1000.; res_tt = res_tt*1000.; res_tl = sqrt(2.)*res_tl*1000.; res_tlp = sqrt(2.)*res_tlp*1000.; ebeam = 5.499; x_b = Q_sqr / (s + Q_sqr - pow(NU_M,2.)); k_out_mom_lab = ebeam - Q_sqr / (2. * NU_M * x_b); k_out_th_lab = 2. * asin(sqrt(Q_sqr / (4. * ebeam * k_out_mom_lab))); ga_mom_in_lab = sqrt(pow(ebeam, 2.) + pow(k_out_mom_lab, 2.) - 2. * ebeam * k_out_mom_lab * cos(k_out_th_lab) ); epsilon = 1. / (1. + 2. * pow(ga_mom_in_lab, 2.) / Q_sqr * pow(tan(0.5 * k_out_th_lab), 2.) ); res_u = res_trans + epsilon*res_long; /*fprintf(fp1, "%f \t %e \t %e \n", a, res_trans, res_trans*conv); fprintf(fp2, "%f \t %e \t %e \n", a, res_long, res_long*conv); fprintf(fp3, "%f \t %e \t %e \n", a, res_tt, res_tt*conv); fprintf(fp4, "%f \t %e \t %e \n", a, res_tl, res_tl*conv); fprintf(fp5, "%f \t %e \t %e \n", a, res_tlp, res_tlp*conv);*/ fprintf(fp1, "%f \t %e \t %e \t %e \t %e \n", a, res_u*conv, res_tl*conv, res_tt*conv, res_tlp*conv); } /*fclose(fp1); fclose(fp2); fclose(fp3); fclose(fp4); fclose(fp5);*/ fclose(fp1); break; } } break; } case 2: /* */ { break; } } break; } /***********************************************************************/ /* */ /* (e, e' K+ SIGMA0) reaction */ /* */ /********************************************************************* */ case 2: /* for (e, e' K+ Sigma0) reaction */ { /* 1. unpolarized forward differential cross section */ /* 2. unpolarized backward differential cross section */ int pol = 1; switch(pol) { case 1: /* unpolarized (e, e' K Sigma) differential cross section : forward */ { int mech, deg_k, deg_kstar, separated; double res_trans, res_long, res_tt, res_tl, res_tlp; printf("\n\n Give choice for FORWARD scattering mechanism"); printf("\n 1 = t-channel kaon"); printf("\n 2 = kstar"); printf("\n 3 = kaon gaugeinvariant"); printf("\n 4 = kaon gaugeinvariant + kstar"); printf("\n"); scanf("%d", &mech); printf("\n\n Give choice for degeneracy of KAON trajectory"); printf("\n 1 = non-degenerate trajectory"); printf("\n 2 = degenerate trajectory with rotating phase (default)"); printf("\n 3 = degenerate trajectory with NON-rotating phase"); printf("\n 4 = single pole"); printf("\n"); scanf("%d", °_k); printf("\n\n Give choice for degeneracy of KSTAR trajectory"); printf("\n 1 = non-degenerate trajectory"); printf("\n 2 = degenerate trajectory with rotating phase (default)"); printf("\n 3 = degenerate trajectory with NON-rotating phase"); printf("\n 4 = single pole"); printf("\n"); scanf("%d", °_kstar); deg_kstar=2; switch(mech) { case 1 : { break; } case 2 : { if (deg_kstar == 2) strcpy(fileb,"_kstar_regge.dat"); if (deg_kstar == 4) strcpy(fileb,"_kstar_born.dat"); break; } case 3 : { if (deg_k == 2) strcpy(fileb,"_k_regge.dat"); if (deg_k == 4) strcpy(fileb,"_k_born.dat"); break; } case 4 : { if ((deg_k == 2) && (deg_kstar == 2)) strcpy(fileb,"_kkstar_regge.dat"); if ((deg_k == 4) && (deg_kstar == 4)) strcpy(fileb,"_kkstar_born.dat"); break; } } printf("\n\n Give choice for electromagnetic form factor of K"); printf("\n 1 = Monopole parametrization"); printf("\n 2 = Monopole at large Q_sqr, correct charge radius at low Q_sqr"); printf("\n 3 = VMD form (rho, omega, phi mesons)"); printf("\n 4 = FF with off-shell extrapolation"); printf("\n"); scanf("%d", &ichoicek); printf("\n\n Give choice for electromagnetic form factor of K*"); printf("\n 1 = Monopole parametrization"); printf("\n 2 = Monopole at large Q_sqr, correct charge radius at low Q_sqr"); printf("\n 3 = VMD form (rho, omega, phi mesons)"); printf("\n"); scanf("%d", &ichoicekst); if((ichoicek == 1)||(ichoicek == 2)) { printf("\n\n K form factor mass scale/cut-off"); printf("\n"); scanf("%lf", &mscalk); } if(ichoicek == 4) { printf("\n Give the value of the dimensionless parameter c in the kaon off-shell FF (e.g. 1) \n"); scanf("%lf",&c_offshell); } if((ichoicekst == 1)||(ichoicekst == 2)) { printf("\n\n K* form factor mass scale/cut-off"); printf("\n"); scanf("%lf", &mscalkst); } printf("\n\n Give choice for OBSERVABLE to calculate"); printf("\n 1 = separated cross section : RESPONSE functions"); printf("\n 2 = theta dependence (Cornell experiments : Bebek et al.)"); printf("\n 3 = q2 dependence"); printf("\n 4 = recoil polarization"); printf("\n 7 = cross sections (T,L,TT,LT) : user-defined kinematics"); printf("\n"); scanf("%d", &separated); switch(separated) { case 1 : /* separated cross section */ { char file1t[80]="ep_ekphyp_w2p15_qsqr1p19"; char file2t[80]; char trans[]="_trans_kkstar.dat"; char longi[]="_long_kkstar.dat"; char tt[]="_tt_pi_gaugeinv.dat"; char tl[]="_tl_pi_gaugeinv.dat"; char tlp[]="_tlp_k_gaugeinv.dat"; char file1l[80],file1tt[80],file1tl[80],file1tlp[80]; int photoelectro; strcpy(file1l,file1t); strcpy(file1tt,file1t); strcpy(file1tl,file1t); strcpy(file1tlp,file1t); printf("\n photo or electro ?"); printf("\n 1 = photo"); printf("\n 2 = electro"); printf("\n"); scanf("%d", &photoelectro); switch(photoelectro) { case 1 : /*photoproduction*/ { double ga_en_lab, bidon1, cos_kth, conv; char file1[80]; char file2[80]; int totdiff; printf("\n total or differential cross section ?"); printf("\n 1 = total"); printf("\n 2 = differential"); printf("\n"); scanf("%d", &totdiff); switch(totdiff) { case 1 : /*photoproduction : total cross section*/ { double tmin, tmax, result1; Q_sqr=0.; sprintf(file1,"gap_kps0_tot.dat"); if ( (fp1 = fopen(file1,"w")) == NULL) fprintf(stderr,"\nUnable to open file"); printf("OUTPUT : Egamma (GeV), W(GeV), sig_tot(mubarn) \n"); for( ga_en_lab= 1.; ga_en_lab <= 4.; ga_en_lab += .25) { s=pow(NU_M,2.)+2.*ga_en_lab*NU_M; w=sqrt(s); tmin=pow((-Q_sqr-pow(K_M,2.)-pow(NU_M,2.)+pow(SIG_M,2.))/2./sqrt(s),2.) -pow(sqrt(pow((s-Q_sqr-pow(NU_M,2.))/2./sqrt(s),2.)+Q_sqr) -sqrt(pow((s+pow(K_M,2.)-pow(SIG_M,2.))/2./sqrt(s),2.)-pow(K_M,2.)),2.); tmax=pow((-Q_sqr-pow(K_M,2.)-pow(NU_M,2.)+pow(SIG_M,2.))/2./sqrt(s),2.) -pow(sqrt(pow((s-Q_sqr-pow(NU_M,2.))/2./sqrt(s),2.)+Q_sqr) +sqrt(pow((s+pow(K_M,2.)-pow(SIG_M,2.))/2./sqrt(s),2.)-pow(K_M,2.)),2.); result1 = sigma_trans_tot(reaction, mech, deg_k, deg_kstar, Q_sqr, s, tmin, tmax); fprintf(fp1,"%f \t %e \t %e \n", ga_en_lab, w, result1); } fclose(fp1); break; } case 2 : /*photoproduction : differential cross section*/ { double ga_en_lab, bidon1, cos_k_th; int wort; char file1[80]; char file2[80]; printf("\n As a function of t or W ?"); printf("\n 1: As a function of t"); printf("\n 2: As a function of W"); printf("\n"); scanf("%d", &wort); if (wort==1) { printf("\n photon energy (in GeV) ?"); printf("\n"); scanf("%lf", &ga_en_lab); bidon1=(ga_en_lab-(int) ga_en_lab +.01)*10.; sprintf(file1,"gap_ks_diff_%dp%dgev.dat", (int) ga_en_lab, (int) bidon1); sprintf(file2,"gap_ks_ass_%dp%dgev.dat", (int) ga_en_lab, (int) bidon1); s = 2. * NU_M * ga_en_lab + pow(NU_M,2.); Q_sqr = 0.; ga_mom_cm = sqrt( ( pow(s - pow(NU_M,2.) + Q_sqr, 2.) + pow(2. * NU_M, 2.) * Q_sqr ) / (4. * s) ); ga_en_cm = sqrt(s) - sqrt(pow(ga_mom_cm,2.) + pow(NU_M,2.)); k_mom_cm = sqrt( ( pow(s - pow(SIG_M,2.) - pow(K_M,2.), 2.) - pow(2. * SIG_M * K_M, 2.) ) / (4. * s) ); k_en = sqrt( pow(k_mom_cm, 2.) + pow(K_M, 2.) ); if ( (fp1 = fopen(file1,"w")) == NULL) fprintf(stderr,"\nUnable to open file"); if ( (fp2 = fopen(file2,"w")) == NULL) fprintf(stderr,"\nUnable to open file"); printf("OUTPUT1 : -t(GeV2), cos(thcm), dsig_T/dt(mubarn/GeV2) \n"); printf("OUTPUT2 : -t(GeV2), cos(thcm), -dsig_TT/dsig_T_dt \n"); for(mt = 0.001; mt <= 1.5; mt += 0.02) { cos_k_th = (- mt - pow(K_M,2.) + 2. * ga_en_cm * k_en) / (2. * ga_mom_cm * k_mom_cm) ; res_trans = response_R_T(reaction, mech, deg_k, deg_kstar, Q_sqr, s, - mt); res_long = response_R_TT(reaction, mech, deg_k, deg_kstar, Q_sqr, s, - mt); res_long = res_long / res_trans; fprintf(fp1, "%f \t %f \t %e \n", mt, cos_k_th, res_trans); fprintf(fp2, "%f \t %f \t %e \n", mt, cos_k_th, res_long); } fclose(fp2); } if (wort==2) { double angle, conv; printf("\n Angle (in degrees) ?"); printf("\n"); scanf("%lf", &angle); bidon1=(angle-(int) angle +.01)*10.; sprintf(file1,"gap_ks_diff_%dp%ddeg.dat", (int) angle, (int) bidon1); Q_sqr = 0.; if ( (fp1 = fopen(file1,"w")) == NULL) fprintf(stderr,"\nUnable to open file"); printf("OUTPUT1 : W(GeV), dsig_T/dt(mubarn/GeV2), dsig_T/dcos(thcm) (mubarn)\n"); for(w = 1.69; w <= 2.4; w += 0.1) { s=pow(w,2.); ga_mom_cm = sqrt( ( pow(s - pow(NU_M,2.) + Q_sqr, 2.) + pow(2. * NU_M, 2.) * Q_sqr ) / (4. * s) ); ga_en_cm = sqrt(s) - sqrt(pow(ga_mom_cm,2.) + pow(NU_M,2.)); k_mom_cm = sqrt( ( pow(s - pow(SIG_M,2.) - pow(K_M,2.), 2.) - pow(2. * SIG_M * K_M, 2.) ) / (4. * s) ); k_en = sqrt( pow(k_mom_cm, 2.) + pow(K_M, 2.) ); mt = - 2. * ga_mom_cm * k_mom_cm *cos(angle*PI/180.) - pow(K_M,2.) + 2. * ga_en_cm * k_en ; conv = 2. * k_mom_cm * ga_mom_cm ; res_trans = response_R_T(reaction, mech, deg_k, deg_kstar, Q_sqr, s, - mt); res_long = response_R_TT(reaction, mech, deg_k, deg_kstar, Q_sqr, s, - mt); res_long = res_long / res_trans; fprintf(fp1, "%f \t %f \t %e \n", w, res_trans, conv * res_trans); } } fclose(fp1); break; } } break; } case 2 : /*electroproduction*/ { int sep; printf("\n Which kinematics ?"); printf("\n Ratio R_L/R_T : "); printf("\n 1 : W = 2.14 GeV, Q^2 = 1.19 GeV^2"); printf("\n 2 : W = 2.66 GeV, Q^2 = 1.99 GeV^2"); printf("\n 3 : W = 2.56 GeV, Q^2 = 3.38 GeV^2"); printf("\n R_T+epsilon R_L, R_TT, R_TL : "); printf("\n 4 : W = 2.21 GeV, Q^2 = .7 GeV^2, epsilon =.85 (Brauel et al.)"); printf("\n 5 : W = 2.21 GeV, Q^2 = 1.35 GeV^2, epsilon =.82 (Brauel et al.)"); printf("\n 6 : W = 2.21 GeV, Q^2 = .06 GeV^2, epsilon =.41 (Brauel et al.)"); printf("\n"); scanf("%d", &sep); switch(sep) { case 1 : /*ratio rlrt*/ { char file1[80]="ep_eks_w2p14_qsqr1p19_rap"; s = pow(2.14, 2.); Q_sqr = 1.19; if ( (fp1 = fopen(strcat(file1,fileb),"w")) == NULL) fprintf(stderr,"\nUnable to open file"); ga_mom_cm = sqrt( ( pow(s - pow(NU_M,2.) + Q_sqr, 2.) + pow(2. * NU_M, 2.) * Q_sqr ) / (4. * s) ); ga_en_cm = sqrt(s) - sqrt(pow(ga_mom_cm,2.) + pow(NU_M,2.)); k_mom_cm = sqrt( ( pow(s - pow(SIG_M,2.) - pow(K_M,2.), 2.) - pow(2. * SIG_M * K_M, 2.) ) / (4. * s) ); k_en = sqrt( pow(k_mom_cm, 2.) + pow(K_M, 2.) ); for(k_th = 0.001; k_th <= 21.; k_th += 2.) { t = - Q_sqr + pow(K_M,2.) - 2. * ga_en_cm * k_en + 2. * ga_mom_cm * k_mom_cm * cos(k_th * PI / 180.); res_trans = response_R_T(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t); res_long = response_R_L(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t); res_long = res_long / res_trans; fprintf(fp1, "%f \t %f \t %e \n", k_th, - t, res_long); } fclose(fp1); break; } case 2 : /*ratio rlrt*/ { char file1[80]="ep_eks_w2p66_qsqr1p99_rap"; s = pow(2.66, 2.); Q_sqr = 1.99; if ( (fp1 = fopen(strcat(file1,fileb),"w")) == NULL) fprintf(stderr,"\nUnable to open file"); ga_mom_cm = sqrt( ( pow(s - pow(NU_M,2.) + Q_sqr, 2.) + pow(2. * NU_M, 2.) * Q_sqr ) / (4. * s) ); ga_en_cm = sqrt(s) - sqrt(pow(ga_mom_cm,2.) + pow(NU_M,2.)); k_mom_cm = sqrt( ( pow(s - pow(SIG_M,2.) - pow(K_M,2.), 2.) - pow(2. * SIG_M * K_M, 2.) ) / (4. * s) ); k_en = sqrt( pow(k_mom_cm, 2.) + pow(K_M, 2.) ); for(k_th = 0.001; k_th <= 21.; k_th += 2.) { t = - Q_sqr + pow(K_M,2.) - 2. * ga_en_cm * k_en + 2. * ga_mom_cm * k_mom_cm * cos(k_th * PI / 180.); res_trans = response_R_T(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t); res_long = response_R_L(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t); res_long = res_long / res_trans; fprintf(fp1, "%f \t %f \t %e \n", k_th, - t, res_long); } fclose(fp1); break; } case 3 : /*ratio rlrt*/ { char file1[80]="ep_eks_w2p56_qsqr3p38_rap"; s = pow(2.56, 2.); Q_sqr = 3.38; if ( (fp1 = fopen(strcat(file1,fileb),"w")) == NULL) fprintf(stderr,"\nUnable to open file"); ga_mom_cm = sqrt( ( pow(s - pow(NU_M,2.) + Q_sqr, 2.) + pow(2. * NU_M, 2.) * Q_sqr ) / (4. * s) ); ga_en_cm = sqrt(s) - sqrt(pow(ga_mom_cm,2.) + pow(NU_M,2.)); k_mom_cm = sqrt( ( pow(s - pow(SIG_M,2.) - pow(K_M,2.), 2.) - pow(2. * SIG_M * K_M, 2.) ) / (4. * s) ); k_en = sqrt( pow(k_mom_cm, 2.) + pow(K_M, 2.) ); for(k_th = 0.001; k_th <= 21.; k_th += 2.) { t = - Q_sqr + pow(K_M,2.) - 2. * ga_en_cm * k_en + 2. * ga_mom_cm * k_mom_cm * cos(k_th * PI / 180.); res_trans = response_R_T(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t); res_long = response_R_L(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t); res_long = res_long / res_trans; fprintf(fp1, "%f \t %f \t %e \n", k_th, - t, res_long); } fclose(fp1); break; } case 4 : /*separation rt+eps*rl, rtt, rtl*/ { char file1[80]="ep_eks_w2p21_qsqr0p7_tpepsl"; char file2[80]="ep_eks_w2p21_qsqr0p7_tt"; char file3[80]="ep_eks_w2p21_qsqr0p7_tl"; s = pow(2.21, 2.); Q_sqr = .7; epsilon = .85; if ( (fp1 = fopen(strcat(file1,fileb),"w")) == NULL) fprintf(stderr,"\nUnable to open file"); if ( (fp2 = fopen(strcat(file2,fileb),"w")) == NULL) fprintf(stderr,"\nUnable to open file"); if ( (fp3 = fopen(strcat(file3,fileb),"w")) == NULL) fprintf(stderr,"\nUnable to open file"); ga_mom_cm = sqrt( ( pow(s - pow(NU_M,2.) + Q_sqr, 2.) + pow(2. * NU_M, 2.) * Q_sqr ) / (4. * s) ); ga_en_cm = sqrt(s) - sqrt(pow(ga_mom_cm,2.) + pow(NU_M,2.)); k_mom_cm = sqrt( ( pow(s - pow(SIG_M,2.) - pow(K_M,2.), 2.) - pow(2. * SIG_M * K_M, 2.) ) / (4. * s) ); k_en = sqrt( pow(k_mom_cm, 2.) + pow(K_M, 2.) ); for(mt = 0.01; mt <= .7; mt += 0.1) { res_trans = response_R_T(reaction, mech, deg_k, deg_kstar, Q_sqr, s, -mt); res_long = response_R_L(reaction, mech, deg_k, deg_kstar, Q_sqr, s, -mt); res_trans = res_trans+epsilon*res_long; res_tt = response_R_TT(reaction, mech, deg_k, deg_kstar, Q_sqr, s, - mt); res_tl = response_R_TL(reaction, mech, deg_k, deg_kstar, Q_sqr, s, - mt); /* printf("%f \t %e \n", mt, res_long);*/ fprintf(fp1, "%f \t %e \n", mt, res_trans); fprintf(fp2, "%f \t %e \n", mt, res_tt); fprintf(fp3, "%f \t %e \n", mt, res_tl); } fclose(fp1); fclose(fp2); fclose(fp3); break; } case 5 : /*separation rt+eps*rl, rtt, rtl*/ { char file1[80]="ep_eks_w2p21_qsqr1p35_tpepsl"; char file2[80]="ep_eks_w2p21_qsqr1p35_tt"; char file3[80]="ep_eks_w2p21_qsqr1p35_tl"; s = pow(2.21, 2.); Q_sqr = 1.35; epsilon = .82; if ( (fp1 = fopen(strcat(file1,fileb),"w")) == NULL) fprintf(stderr,"\nUnable to open file"); if ( (fp2 = fopen(strcat(file2,fileb),"w")) == NULL) fprintf(stderr,"\nUnable to open file"); if ( (fp3 = fopen(strcat(file3,fileb),"w")) == NULL) fprintf(stderr,"\nUnable to open file"); ga_mom_cm = sqrt( ( pow(s - pow(NU_M,2.) + Q_sqr, 2.) + pow(2. * NU_M, 2.) * Q_sqr ) / (4. * s) ); ga_en_cm = sqrt(s) - sqrt(pow(ga_mom_cm,2.) + pow(NU_M,2.)); k_mom_cm = sqrt( ( pow(s - pow(SIG_M,2.) - pow(K_M,2.), 2.) - pow(2. * SIG_M * K_M, 2.) ) / (4. * s) ); k_en = sqrt( pow(k_mom_cm, 2.) + pow(K_M, 2.) ); for(mt = 0.01; mt <= .7; mt += 0.1) { res_trans = response_R_T(reaction, mech, deg_k, deg_kstar, Q_sqr, s, -mt); res_long = response_R_L(reaction, mech, deg_k, deg_kstar, Q_sqr, s, -mt); res_trans = res_trans+epsilon*res_long; res_tt = response_R_TT(reaction, mech, deg_k, deg_kstar, Q_sqr, s, - mt); res_tl = response_R_TL(reaction, mech, deg_k, deg_kstar, Q_sqr, s, - mt); fprintf(fp1, "%f \t %e \n", mt, res_trans); fprintf(fp2, "%f \t %e \n", mt, res_tt); fprintf(fp3, "%f \t %e \n", mt, res_tl); } fclose(fp1); fclose(fp2); fclose(fp3); break; } case 6 : /*separation rt+eps*rl, rtt, rtl*/ { char file1[80]="ep_eks_w2p21_qsqr0p06_tpepsl"; char file2[80]="ep_eks_w2p21_qsqr0p06_tt"; char file3[80]="ep_eks_w2p21_qsqr0p06_tl"; s = pow(2.21, 2.); Q_sqr = .06; epsilon = .41; if ( (fp1 = fopen(strcat(file1,fileb),"w")) == NULL) fprintf(stderr,"\nUnable to open file"); if ( (fp2 = fopen(strcat(file2,fileb),"w")) == NULL) fprintf(stderr,"\nUnable to open file"); if ( (fp3 = fopen(strcat(file3,fileb),"w")) == NULL) fprintf(stderr,"\nUnable to open file"); ga_mom_cm = sqrt( ( pow(s - pow(NU_M,2.) + Q_sqr, 2.) + pow(2. * NU_M, 2.) * Q_sqr ) / (4. * s) ); ga_en_cm = sqrt(s) - sqrt(pow(ga_mom_cm,2.) + pow(NU_M,2.)); k_mom_cm = sqrt( ( pow(s - pow(SIG_M,2.) - pow(K_M,2.), 2.) - pow(2. * SIG_M * K_M, 2.) ) / (4. * s) ); k_en = sqrt( pow(k_mom_cm, 2.) + pow(K_M, 2.) ); for(mt = 0.01; mt <= .7; mt += 0.1) { res_trans = response_R_T(reaction, mech, deg_k, deg_kstar, Q_sqr, s, -mt); res_long = response_R_L(reaction, mech, deg_k, deg_kstar, Q_sqr, s, -mt); /* res_trans = res_trans+epsilon*res_long;*/ res_tt = response_R_TT(reaction, mech, deg_k, deg_kstar, Q_sqr, s, - mt); res_tl = response_R_TL(reaction, mech, deg_k, deg_kstar, Q_sqr, s, - mt); fprintf(fp1, "%f \t %e \n", mt, res_trans); fprintf(fp2, "%f \t %e \n", mt, res_tt); fprintf(fp3, "%f \t %e \n", mt, res_tl); } fclose(fp1); fclose(fp2); fclose(fp3); break; } } } } break; } case 2 : /* theta dependence (Cornell experiments : Bebek et al.) */ { int i, kin, choiceintegr; double cos_min,mt_min,lbin; printf("\n\n Choose the kinematics of the Cornell experiments (Bebek et al.)"); printf("\n 1 : W = 2.2 GeV, Q^2 = 1.18 GeV^2, epsilon = 0.94"); printf("\n 2 : W = 2.17 GeV, Q^2 = 1.97 GeV^2, epsilon = 0.94"); printf("\n 3 : W = 2.19 GeV, Q^2 = 3.95 GeV^2, epsilon = 0.88"); printf("\n 4 : W = 2.78 GeV, Q^2 = 1.38 GeV^2, epsilon = 0.88"); printf("\n 5 : W = 2.47 GeV, Q^2 = 3.52 GeV^2, epsilon = 0.86"); printf("\n 6 : W = 2.66 GeV, Q^2 = 2.04 GeV^2, epsilon = 0.88"); printf("\n 7 : W = 2.21 GeV, Q^2 = .06 GeV^2, epsilon = 0.41 (Brauel et al.)"); printf("\n 8 : W = 2.21 GeV, Q^2 = .28 GeV^2, epsilon = 0.74 (Brauel et al.)"); printf("\n 9 : W = 2.21 GeV, Q^2 = .7 GeV^2, epsilon = 0.85 (Brauel et al.)"); printf("\n 10 : W = 2.21 GeV, Q^2 = 1.35 GeV^2, epsilon = 0.82 (Brauel et al.)"); printf("\n"); printf("\n"); scanf("%d", &kin); /* printf("\n\n Integration on Phi (120=>240) or not?"); printf("\n 0 : no"); printf("\n 1 : yes"); printf("\n"); scanf("%d", &choiceintegr);*/ choiceintegr=1; switch(kin) { case 1 : { strcpy(file2,"ep_eks_w2p2_qsqr1p18_eps0p94"); w=2.2; Q_sqr=1.18; epsilon=.94; break; } case 2 : { strcpy(file2,"ep_eks_w2p17_qsqr1p97_eps0p94"); w=2.17; Q_sqr=1.97; epsilon=.94; break; } case 3 : { strcpy(file2,"ep_eks_w2p19_qsqr3p95_eps0p88"); w=2.19; Q_sqr=3.95; epsilon=.88; break; } case 4 : { strcpy(file2,"ep_eks_w2p78_qsqr1p38_eps0p88"); w=2.78; Q_sqr=1.38; epsilon=.88; break; } case 5 : { strcpy(file2,"ep_eks_w2p47_qsqr3p52_eps0p86"); w=2.47; Q_sqr=3.52; epsilon=.86; break; } case 6 : { strcpy(file2,"ep_eks_w2p66_qsqr2p04_eps0p88"); w=2.66; Q_sqr=2.04; epsilon=.88; /* w=2.21; Q_sqr=.06; epsilon=.41;*/ break; } case 7 : { strcpy(file2,"ep_eks_w2p21_qsqr0p06_eps0p41"); w=2.21; Q_sqr=.06; epsilon=.41; break; } case 8 : { strcpy(file2,"ep_eks_w2p21_qsqr0p28_eps0p74"); w=2.21; Q_sqr=.28; epsilon=.74; break; } case 9 : { strcpy(file2,"ep_eks_w2p21_qsqr0p70_eps0p85"); w=2.21; Q_sqr=.7; epsilon=.85; break; } case 10 : { strcpy(file2,"ep_eks_w2p21_qsqr1p35_eps0p82"); w=2.21; Q_sqr=1.35; epsilon=.82; break; } } if ( (fp2 = fopen(strcat(file2,fileb),"w")) == NULL) fprintf(stderr,"\nUnable to open file"); s = pow(w, 2.); ga_mom_cm = sqrt( ( pow(s - pow(NU_M,2.) + Q_sqr, 2.) + pow(2. * NU_M, 2.) * Q_sqr ) / (4. * s) ); ga_en_cm = sqrt(s) - sqrt(pow(ga_mom_cm,2.) + pow(NU_M,2.)); k_mom_cm = sqrt( ( pow(s - pow(SIG_M,2.) - pow(K_M,2.), 2.) - pow(2. * SIG_M * K_M, 2.) ) / (4. * s) ); k_en = sqrt( pow(k_mom_cm, 2.) + pow(K_M, 2.) ); if (kin<7){ for(k_th = 0.001; k_th <= 21.; k_th += 2.) /* for(k_th = 0.001; k_th <= 180.; k_th += 2.)*/ { double conv; t = - Q_sqr + pow(K_M,2.) - 2. * ga_en_cm * k_en + 2. * ga_mom_cm * k_mom_cm * cos(k_th * PI / 180.); conv = 2. * k_mom_cm * ga_mom_cm / (2.*PI) ; diff = conv * dsigma_dt_anal(reaction, mech, deg_k, deg_kstar, 0., epsilon, Q_sqr, s, t); /* fprintf(fp2, "%f \t %f \t %e \n", cos(k_th*PI/180.), -t, diff); */ fprintf(fp2, "%f \t %f \t %e \n", k_th, -t, diff); } fclose(fp2); } mt_min=0.; do{ mt_min=mt_min+.01; cos_min = (- Q_sqr + pow(K_M,2.) - 2. * ga_en_cm * k_en + mt_min)/ (-2. * ga_mom_cm * k_mom_cm) ; }while(fabs(cos_min)>1.); if ((kin>6)&&(kin<11)){ phi_min = 120. * PI / 180.; phi_max = 240. * PI / 180.; /* phi_min = 0. * PI / 180.; phi_max = 360. * PI / 180.;*/ lbin=phi_max-phi_min; for(mt = mt_min; mt <= 1.7; mt += 0.2) { if (choiceintegr==0) { diff = 2.*PI* dsigma_dtdphi(reaction, mech, deg_k, deg_kstar, 0., PI ,epsilon, Q_sqr, s, -mt); difft = 2.*PI* dsigma_dtdphi(reaction, mech, deg_k, deg_kstar, 0., PI ,0., Q_sqr, s, -mt); diffl = 2.*PI* dsigmal_dtdphi(reaction, mech, deg_k, deg_kstar, 0., PI ,epsilon, Q_sqr, s, -mt); } else { diff = 2. * PI * dsigma_dt(reaction, mech, deg_k, deg_kstar, 0., phi_min, phi_max, epsilon, Q_sqr, s, - mt); difft = 2. * PI * dsigma_dt(reaction, mech, deg_k, deg_kstar, 0., phi_min, phi_max, 0. , Q_sqr, s, - mt); diffl = 2. * PI * dsigmal_dt(reaction, mech, deg_k, deg_kstar, 0., phi_min, phi_max, epsilon , Q_sqr, s, - mt); difftt = 2. * PI * dsigmatt_dt(reaction, mech, deg_k, deg_kstar, 0., phi_min, phi_max, epsilon , Q_sqr, s, - mt); difftl = 2. * PI * dsigmatl_dt(reaction, mech, deg_k, deg_kstar, 0., phi_min, phi_max, epsilon , Q_sqr, s, - mt); diff=diff/lbin; difft=difft/lbin; diffl=diffl/lbin; difftt=difftt/lbin; difftl=difftl/lbin; } if (diff!=0) fprintf(fp2, "%f \t %e \t %e \t %e \t %e \t %e \n", mt, diff, difft, diffl, difftt, difftl); } fclose(fp2); } break; } case 3 : /* q2 dependence */ { int obsq2; printf("\n\n Which observable ?"); printf("\n 1 = dsig/dt (Bebek et al.)"); printf("\n 2 = dsigT/dt, dsigL/dt, sigL/sigT (CEBAF)"); printf("\n 3 = dsigT/dt, dsigL/dt, sigL/sigT (Bydzo.)"); printf("\n"); scanf("%d", &obsq2); switch(obsq2) { case 1: /* dsigdt */ { int kin; double diff_test; char filea[80]; printf("\n\n Choose the angle"); printf("\n 1 : K Sigma electroproduction : Q2 dependence"); printf("\n : W = 2.17 GeV, theta = 7 deg, epsilon = 0.85"); printf("\n 2 : K Sigma electroproduction : Q2 dependance"); printf("\n : W = 2.17 GeV, theta = 1 deg, epsilon = 0.85"); printf("\n 3 : K Sigma electroproduction : Q2 dependance"); printf("\n : W = 2.17 GeV, theta = 15 deg, epsilon = 0.85"); printf("\n 4 : K Sigma electroproduction : Q2 dependance"); printf("\n : W = 2.66 GeV, theta = 7 deg, epsilon = 0.85"); printf("\n"); scanf("%d", &kin); if (kin==1) { strcpy(filea,"ep_eks_w2p17_the7p0_eps0p85_q2"); k_th=7.; } if (kin==2) { strcpy(filea,"ep_eks_w2p17_the1p0_eps0p85_q2"); k_th=1.; } if (kin==3) { strcpy(filea,"ep_eks_w2p17_the15p0_eps0p85_q2"); k_th=15.; } if (kin==4) { strcpy(filea,"ep_eks_w2p66_the8p0_eps0p85_q2"); k_th=8.; } if (kin<=3) { s = pow(2.17, 2.); } if (kin==4) { s = pow(2.66, 2.); } epsilon = 0.85; if ( (fp1 = fopen(strcat(filea,fileb),"w")) == NULL) fprintf(stderr,"\nUnable to open file"); for(Q_sqr = 0.001; Q_sqr <= 5.; Q_sqr += 0.1) { double conv; ga_mom_cm = sqrt( ( pow(s - pow(NU_M,2.) + Q_sqr, 2.) + pow(2. * NU_M, 2.) * Q_sqr ) / (4. * s) ); ga_en_cm = sqrt(s) - sqrt(pow(ga_mom_cm,2.) + pow(NU_M,2.)); k_mom_cm = sqrt( ( pow(s - pow(SIG_M,2.) - pow(K_M,2.), 2.) - pow(2. * SIG_M * K_M, 2.) ) / (4. * s) ); k_en = sqrt( pow(k_mom_cm, 2.) + pow(K_M, 2.) ); t = - Q_sqr + pow(K_M,2.) - 2. * ga_en_cm * k_en + 2. * ga_mom_cm * k_mom_cm * cos(k_th * PI / 180.); conv = 2. * k_mom_cm * ga_mom_cm / (2. * PI); diff = conv * dsigma_dt_anal(reaction, mech, deg_k, deg_kstar, 0., epsilon, Q_sqr, s, t); diffl = conv * response_R_L(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t); fprintf(fp1, "%f \t %e \t %e \n", Q_sqr, diff, diffl); } fclose(fp1); break; } case 2 : /*ratio siglsigt*/ { char file1[80]="ep_eks_cebaf_rap.dat"; char file2[80]="ep_eks_cebaf_seplt.dat"; double rap, conv, cos_min; int iii; if ( (fp2 = fopen(file2,"w")) == NULL) fprintf(stderr,"\nUnable to open file"); for(iii = 1; iii <= 6; iii += 1) { if (iii==1) { s = pow(1.84, 2.); Q_sqr = .1; t = -.15; } if (iii==2) { s = pow(1.84, 2.); Q_sqr = .52; t = -.219; } if (iii==3) { s = pow(1.83, 2.); Q_sqr = .75; t = -.3; } if (iii==4) { s = pow(1.81, 2.); Q_sqr = 1.; t = -.407; } if (iii==5) { s = pow(1.84, 2.); Q_sqr = 2.; t = -.741; } if (iii==6) { s = pow(1.84, 2.); Q_sqr = 3.9; t = -1.4; } ga_mom_cm = sqrt( ( pow(s - pow(NU_M,2.) + Q_sqr, 2.) + pow(2. * NU_M, 2.) * Q_sqr ) / (4. * s) ); k_mom_cm = sqrt( (pow(s - pow(SIG_M,2.) - pow(K_M,2.), 2.) - pow(2. * SIG_M * K_M,2.) ) /(4. * s)); ga_en_cm = sqrt(s) - sqrt(pow(ga_mom_cm,2.)+pow(NU_M,2.)); k_en = sqrt( pow(k_mom_cm, 2.) + pow(K_M, 2.) ); conv = 2. * k_mom_cm * ga_mom_cm / (2.*PI) ; t=t+.01; do{ t=t-.01; cos_min = (- Q_sqr + pow(K_M,2.) - 2. * ga_en_cm * k_en -t)/(-2. * ga_mom_cm * k_mom_cm) ; }while(fabs(cos_min)>1.); res_trans = response_R_T(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t); res_long = response_R_L(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t); rap = res_long / res_trans; res_trans = res_trans*conv; res_long = res_long*conv; fprintf(fp2, "%f \t %f \t %f \t %f \t %f \n", t, Q_sqr, res_trans, res_long, rap); } fclose(fp2); break; } case 3: /* dsigdt */ { int kin; double diff_test; char filea[80]; printf("\n\n Choose the kinematics"); printf("\n 1 : K Sigma electroproduction : Q2 dependence"); printf("\n : W = 2.24 GeV, t = -.15, epsilon = 0.72"); printf("\n"); scanf("%d", &kin); if (kin==1) { strcpy(filea,"ep_eks_w2p24_t0p15_eps0p72_q2"); t=-.15; s = pow(2.24, 2.); epsilon = 0.72; } if ( (fp1 = fopen(strcat(filea,fileb),"w")) == NULL) fprintf(stderr,"\nUnable to open file"); for(Q_sqr = 0.001; Q_sqr <= 1.; Q_sqr += 0.05) { double conv; ga_mom_cm = sqrt( ( pow(s - pow(NU_M,2.) + Q_sqr, 2.) + pow(2. * NU_M, 2.) * Q_sqr ) / (4. * s) ); ga_en_cm = sqrt(s) - sqrt(pow(ga_mom_cm,2.) + pow(NU_M,2.)); k_mom_cm = sqrt( ( pow(s - pow(SIG_M,2.) - pow(K_M,2.), 2.) - pow(2. * SIG_M * K_M, 2.) ) / (4. * s) ); k_en = sqrt( pow(k_mom_cm, 2.) + pow(K_M, 2.) ); /* conv = 2. * k_mom_cm * ga_mom_cm / (2. * PI);*/ conv = 1.; diff = conv * dsigma_dt_anal(reaction, mech, deg_k, deg_kstar, 0.,epsilon,Q_sqr,s,t); difft = conv * response_R_T(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t); diffl = conv * response_R_L(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t); fprintf(fp1, "%f \t %e \t %e \t %e \n", Q_sqr, diff, difft,diffl); } fclose(fp1); break; } } break; } case 4 : /* recoil polarization */ { double asymm_lam, asymm_lam2, tortho, lortho, res_trans, res_long, W, epsilon ; printf("\n Which W ?"); printf("\n"); scanf("%lf", &W); printf("\n Which Q2 ?"); printf("\n"); scanf("%lf", &Q_sqr); printf("\n Which epsilon ?"); printf("\n"); scanf("%lf", &epsilon); if ( (fp1 = fopen("recoilsig.dat", "w")) == NULL) fprintf(stderr,"\nUnable to open file"); s=pow(W,2.); ga_mom_cm = sqrt( ( pow(s - pow(NU_M,2.) + Q_sqr, 2.) + pow(2. * NU_M, 2.) * Q_sqr ) / (4. * s) ); ga_en_cm = sqrt(s) - sqrt(pow(ga_mom_cm,2.) + pow(NU_M,2.)); k_mom_cm = sqrt( ( pow(s - pow(SIG_M,2.) - pow(K_M,2.), 2.) - pow(2. * SIG_M * K_M, 2.) ) / (4. * s) ); k_en = sqrt( pow(k_mom_cm, 2.) + pow(K_M, 2.) ); for(k_th = 0.; k_th <= 180.; k_th += 5.) { t = - Q_sqr + pow(K_M,2.) - 2. * ga_en_cm * k_en + 2. * ga_mom_cm * k_mom_cm * cos(k_th * PI / 180.); asymm_lam2 = recoil_asymm(reaction,mech,deg_k,deg_kstar, Q_sqr, s, t); tortho = response_R_T_ortho(reaction,mech,deg_k,deg_kstar, Q_sqr, s, t); res_trans = response_R_T(reaction,mech,deg_k,deg_kstar, Q_sqr, s, t); lortho = response_R_L_ortho(reaction,mech,deg_k,deg_kstar, Q_sqr, s, t); res_long = response_R_L(reaction,mech,deg_k,deg_kstar, Q_sqr, s, t); if (Q_sqr==0.) { asymm_lam = -tortho/res_trans; lortho=0.; res_long=0.; } else { asymm_lam = (tortho+epsilon*lortho) / (res_trans+epsilon*res_long); } fprintf(fp1,"%f \t %f \t %f \t %f \t %f \t %f \t %f \t %f\n", -t, asymm_lam, asymm_lam2, tortho, lortho, res_trans, res_long, k_th ); } fclose(fp); break; } case 7 : /* User-defined */ { double a, a_min, a_max, a_step, W_min, W_max, W_step, t_min, t_max, t_step, Q2_min, Q2_max, Q2_step, bidon1, bidon2, cos_th_cm, cth_min, cth_max, cth_step; double ga_en_in_lab, ga_mom_in_lab, pi_mom_cm, W, tmin_test, mes_mom_lab_test, mes_en_lab, mes_mom_lab; double x_b,k_out_mom_lab,k_out_th_lab,ebeam,s,epsilon; double res_u; int rep, rep2; char file1t[80],file1l[80],file1tt[80],file1tl[80],file1tlp[80],file1u[80]; printf("\n 1 : Q2 dependence"); printf("\n 2 : t dependence"); printf("\n 3 : W dependence"); printf("\n 4 : cthkcm dependence"); printf("\n"); scanf("%d", &rep); if (rep==1) { printf("\n Which W"); printf("\n"); scanf("%lf", &W); s = pow(W,2.); printf("\n Which cos_th_cm ?"); printf("\n"); scanf("%lf", &cos_th_cm); /*printf("\n Q2_min ?"); printf("\n"); scanf("%lf", &Q2_min); printf("\n Which Q2_max ?"); printf("\n"); scanf("%lf", &Q2_max); printf("\n Which Q2_step ?"); printf("\n"); scanf("%lf", &Q2_step);*/ Q2_min=0.2; Q2_max=5.0; Q2_step=0.1; a_min=Q2_min; a_max=Q2_max; a_step=Q2_step; bidon1=(W-(int) W)*10.; bidon2=(t-(int) t)*10.; /*sprintf(file1t,"ep_ksigma_t%dp%d_w%dp%d_trans", (int) t, (int) bidon2, (int) W, (int) bidon1); sprintf(file1l,"ep_ksigma_t%dp%d_w%dp%d_long", (int) t, (int) bidon2, (int) W, (int) bidon1); sprintf(file1tt,"ep_ksigma_t%dp%d_w%dp%d_tt", (int) t, (int) bidon2, (int) W, (int) bidon1); sprintf(file1tl,"ep_ksigma_t%dp%d_w%dp%d_tl", (int) t, (int) bidon2, (int) W, (int) bidon1); printf("OUTPUT1 : Q2 (GeV2), dsig_T/dt (mubarn/GeV2), dsig_T/dOmega_cm (mubarn/sr) \n"); printf("OUTPUT2 : Q2 (GeV2), dsig_L/dt (mubarn/GeV2), dsig_L/dOmega_cm (mubarn/sr) \n"); printf("OUTPUT3 : Q2 (GeV2), dsig_TT/dt (mubarn/GeV2), dsig_TT/dOmega_cm (mubarn/sr) \n"); printf("OUTPUT4 : Q2 (GeV2), dsig_TL/dt (mubarn/GeV2), dsig_TL/dOmega_cm (mubarn/sr) \n");*/ sprintf(file1u,"ep_ksigma_q2", (int) Q_sqr, (int) bidon1, (int) t, (int) bidon2); printf("OUTPUT : Q2 (GeV2), sigu, siglt, sigtt, sigltp (nb/sr) \n"); } if (rep==2) { printf("\n Which W"); printf("\n"); scanf("%lf", &W); printf("\n Which Q2"); printf("\n"); scanf("%lf", &Q_sqr); s = pow(W,2.); printf("\n t_min ? (if 0 => kinematical t_min)"); printf("\n"); scanf("%lf", &t_min); printf("\n t_max ?"); printf("\n"); scanf("%lf", &t_max); printf("\n t_step ?"); printf("\n"); scanf("%lf", &t_step); a_min=t_min; a_max=t_max; a_step=t_step; if(t_min==0.) { a_min=pow((-Q_sqr-pow(K_M,2.)-pow(NU_M,2.)+pow(SIG_M,2.))/2./sqrt(s),2.) -pow(sqrt(pow((s-Q_sqr-pow(NU_M,2.))/2./sqrt(s),2.)+Q_sqr) -sqrt(pow((s+pow(K_M,2.)-pow(SIG_M,2.))/2./sqrt(s),2.)-pow(K_M,2.)),2.); a_min=-a_min; } bidon1=(Q_sqr-(int) Q_sqr)*10.; bidon2=(W-(int) W)*10.; sprintf(file1t,"ep_ksigma_q2%dp%d_w%dp%d_trans", (int) Q_sqr, (int) bidon1, (int) W, (int) bidon2); sprintf(file1l,"ep_ksigma_q2%dp%d_w%dp%d_long", (int) Q_sqr, (int) bidon1, (int) W, (int) bidon2); sprintf(file1tt,"ep_ksigma_q2%dp%d_w%dp%d_tt", (int) Q_sqr, (int) bidon1, (int) W, (int) bidon2); sprintf(file1tl,"ep_ksigma_q2%dp%d_w%dp%d_tl", (int) Q_sqr, (int) bidon1, (int) W, (int) bidon2); printf("OUTPUT1 : -t (GeV2), dsig_T/dt (mubarn/GeV2), dsig_T/dOmega_cm (mubarn/sr) \n"); printf("OUTPUT2 : -t (GeV2), dsig_L/dt (mubarn/GeV2), dsig_L/dOmega_cm (mubarn/sr) \n"); printf("OUTPUT3 : -t (GeV2), dsig_TT/dt (mubarn/GeV2), dsig_TT/dOmega_cm (mubarn/sr) \n"); printf("OUTPUT4 : -t (GeV2), dsig_TL/dt, (mubarn/GeV2),dsig_TL/dOmega_cm (mubarn/sr) \n"); } if (rep==3) { printf("\n Which Q2"); printf("\n"); scanf("%lf", &Q_sqr); printf("\n Fixed t or fixed cos_th_cm"); printf("\n 1 : Fixed t"); printf("\n 2 : Fixed cos_th_cm"); printf("\n"); scanf("%d", &rep2); if (rep2==1) { printf("\n Which t (t>0) -if t=0, t=t_min-"); printf("\n"); scanf("%lf", &t); mt=-t; } if (rep2==2) { printf("\n Which cos_th_cm ?"); printf("\n"); scanf("%lf", &cos_th_cm); } /*printf("\n W_min ?"); printf("\n"); scanf("%lf", &W_min); printf("\n W_max ?"); printf("\n"); scanf("%lf", &W_max); printf("\n W_step ?"); printf("\n"); scanf("%lf", &W_step);*/ W_min=1.69; W_max=2.80; W_step=0.005; a_min=W_min; a_max=W_max; a_step=W_step; bidon1=(Q_sqr-(int) Q_sqr)*10.; bidon2=(t-(int) t)*10.; /*sprintf(file1t,"ep_ksigma_q2%dp%d_t%dp%d_trans", (int) Q_sqr, (int) bidon1, (int) t, (int) bidon2); sprintf(file1l,"ep_ksigma_q2%dp%d_t%dp%d_long", (int) Q_sqr, (int) bidon1, (int) t, (int) bidon2); sprintf(file1tt,"ep_ksigma_q2%dp%d_t%dp%d_tt", (int) Q_sqr, (int) bidon1, (int) t, (int) bidon2); sprintf(file1tl,"ep_ksigma_q2%dp%d_t%dp%d_tl", (int) Q_sqr, (int) bidon1, (int) t, (int) bidon2); sprintf(file1tlp,"ep_ksigma_q2%dp%d_t%dp%d_tlp", (int) Q_sqr, (int) bidon1, (int) t, (int) bidon2); printf("OUTPUT1 : W (GeV), dsig_T/dt (nb/GeV2), dsig_T/dOmega_cm (nb/sr) \n"); printf("OUTPUT2 : W (GeV), dsig_L/dt (nb/GeV2), dsig_L/dOmega_cm (nb/sr) \n"); printf("OUTPUT3 : W (GeV), dsig_TT/dt (nb/GeV2), dsig_TT/dOmega_cm (nb/sr) \n"); printf("OUTPUT4 : W (GeV), dsig_TL/dt (nb/GeV2), dsig_TL/dOmega_cm (nb/sr) \n"); printf("OUTPUT5 : W (GeV), dsig_TLp/dt (nb/GeV2), dsig_TLp/dOmega_cm (nb/sr) \n");*/ sprintf(file1u,"ep_ksigma_w", (int) Q_sqr, (int) bidon1, (int) t, (int) bidon2); printf("OUTPUT : W (GeV), sigu, siglt, sigtt, sigltp (nb/sr) \n"); } if (rep==4) { printf("\n Which Q2"); printf("\n"); scanf("%lf", &Q_sqr); printf("\n Which W"); printf("\n"); scanf("%lf", &W); s = pow(W,2.); /*printf("\n cth_min ?"); printf("\n"); scanf("%lf", &cth_min); printf("\n cth_max ?"); printf("\n"); scanf("%lf", &cth_max); printf("\n cth_step ?"); printf("\n"); scanf("%lf", &cth_step);*/ cth_min=-1.0; cth_max=1.0; cth_step=0.02; a_min=cth_min; a_max=cth_max; a_step=cth_step; bidon1=(Q_sqr-(int) Q_sqr)*10.; bidon2=(W-(int) W)*10.; /*sprintf(file1t,"ep_ksigma_q2%dp%d_t%dp%d_trans", (int) Q_sqr, (int) bidon1, (int) W, (int) bidon2); sprintf(file1l,"ep_ksigma_q2%dp%d_t%dp%d_long", (int) Q_sqr, (int) bidon1, (int) W, (int) bidon2); sprintf(file1tt,"ep_ksigma_q2%dp%d_t%dp%d_tt", (int) Q_sqr, (int) bidon1, (int) W, (int) bidon2); sprintf(file1tl,"ep_ksigma_q2%dp%d_t%dp%d_tl", (int) Q_sqr, (int) bidon1, (int) W, (int) bidon2); sprintf(file1tlp,"ep_ksigma_q2%dp%d_t%dp%d_tlp", (int) Q_sqr, (int) bidon1, (int) W, (int) bidon2); printf("OUTPUT1 : W (GeV), dsig_T/dt (nb/GeV2), dsig_T/dOmega_cm (nb/sr) \n"); printf("OUTPUT2 : W (GeV), dsig_L/dt (nb/GeV2), dsig_L/dOmega_cm (nb/sr) \n"); printf("OUTPUT3 : W (GeV), dsig_TT/dt (nb/GeV2), dsig_TT/dOmega_cm (nb/sr) \n"); printf("OUTPUT4 : W (GeV), dsig_TL/dt (nb/GeV2), dsig_TL/dOmega_cm (nb/sr) \n"); printf("OUTPUT5 : W (GeV), dsig_TLp/dt (nb/GeV2), dsig_TLp/dOmega_cm (nb/sr) \n");*/ sprintf(file1u,"ep_ksigma_a", (int) Q_sqr, (int) bidon1, (int) t, (int) bidon2); printf("OUTPUT : cthkcm, sigu, siglt, sigtt, sigltp (nb/sr) \n"); } /*if ( (fp1 = fopen(strcat(file1t,fileb),"w")) == NULL) fprintf(stderr,"\nUnable to open file"); if ( (fp2 = fopen(strcat(file1l,fileb),"w")) == NULL) fprintf(stderr,"\nUnable to open file"); if ( (fp3 = fopen(strcat(file1tt,fileb),"w")) == NULL) fprintf(stderr,"\nUnable to open file"); if ( (fp4 = fopen(strcat(file1tl,fileb),"w")) == NULL) fprintf(stderr,"\nUnable to open file"); if ( (fp5 = fopen(strcat(file1tlp,fileb),"w")) == NULL) fprintf(stderr,"\nUnable to open file");*/ if ( (fp1 = fopen(strcat(file1u,fileb),"w")) == NULL) fprintf(stderr,"\nUnable to open file"); for(a = a_min; a <= a_max; a += a_step) { double conv, cos_k, ga_en_cm2; if (rep==1) Q_sqr=a; if (rep==2) mt=-a; if (rep==3) s=pow(a,2.); if (rep==4) cos_th_cm=a; if (rep==1) { ga_en_cm = sqrt(s) - sqrt(pow(ga_mom_cm,2.) + pow(NU_M,2.)); k_en = sqrt( pow(pi_mom_cm, 2.) + pow(K_M, 2.) ); mt = -Q_sqr + pow(K_M,2.) - 2. * ga_en_cm * k_en + 2. * ga_mom_cm * pi_mom_cm * cos_th_cm; } if((t==0)&&(rep!=2)) { mt=pow((-Q_sqr-pow(K_M,2.)-pow(NU_M,2.)+pow(SIG_M,2.))/2./sqrt(s),2.) -pow(sqrt(pow((s-Q_sqr-pow(NU_M,2.))/2./sqrt(s),2.)+Q_sqr) -sqrt(pow((s+pow(K_M,2.)-pow(SIG_M,2.))/2./sqrt(s),2.)-pow(K_M,2.)),2.); } ga_mom_cm = sqrt( ( pow(s - pow(NU_M,2.) + Q_sqr, 2.) + pow(2. * NU_M, 2.) * Q_sqr ) / (4. * s) ); pi_mom_cm = sqrt( ( pow(s - pow(SIG_M,2.) - pow(K_M,2.), 2.) - pow(2. * SIG_M * K_M, 2.) ) / (4. * s) ); if (rep2==2) { ga_en_cm = sqrt(s) - sqrt(pow(ga_mom_cm,2.) + pow(NU_M,2.)); k_en = sqrt( pow(pi_mom_cm, 2.) + pow(K_M, 2.) ); mt = -Q_sqr + pow(K_M,2.) - 2. * ga_en_cm * k_en + 2. * ga_mom_cm * pi_mom_cm * cos_th_cm; } if (rep==4) { ga_en_cm = sqrt(s) - sqrt(pow(ga_mom_cm,2.) + pow(NU_M,2.)); k_en = sqrt( pow(pi_mom_cm, 2.) + pow(K_M, 2.) ); mt = -Q_sqr + pow(K_M,2.) - 2. * ga_en_cm * k_en + 2. * ga_mom_cm * pi_mom_cm * cos_th_cm; } conv = 1. / (2. * PI) * 2. * pi_mom_cm * ga_mom_cm; res_trans = response_R_T(reaction, mech, deg_k, deg_kstar, Q_sqr, s, mt); res_tt = response_R_TT(reaction, mech, deg_k, deg_kstar, Q_sqr, s, mt); if (Q_sqr==0.) { res_long=0.; res_tl=0.; res_tlp=0.; } else { res_long = response_R_L(reaction, mech, deg_k, deg_kstar, Q_sqr, s, mt); res_tl = response_R_TL(reaction, mech, deg_k, deg_kstar, Q_sqr, s, mt); res_tlp = response_R_TLp(reaction, mech, deg_k, deg_kstar, Q_sqr, s, mt); } /* Changes based on e1f paper conventions. Also use nb/sr for output */ res_trans = res_trans*1000.; res_long = res_long*1000.; res_tt = res_tt*1000.; res_tl = sqrt(2.)*res_tl*1000.; res_tlp = sqrt(2.)*res_tlp*1000.; ebeam = 5.499; x_b = Q_sqr / (s + Q_sqr - pow(NU_M,2.)); k_out_mom_lab = ebeam - Q_sqr / (2. * NU_M * x_b); k_out_th_lab = 2. * asin(sqrt(Q_sqr / (4. * ebeam * k_out_mom_lab))); ga_mom_in_lab = sqrt(pow(ebeam, 2.) + pow(k_out_mom_lab, 2.) - 2. * ebeam * k_out_mom_lab * cos(k_out_th_lab) ); epsilon = 1. / (1. + 2. * pow(ga_mom_in_lab, 2.) / Q_sqr * pow(tan(0.5 * k_out_th_lab), 2.) ); res_u = res_trans + epsilon*res_long; /*fprintf(fp1, "%f \t %e \t %e \n", a, res_trans, res_trans*conv); fprintf(fp2, "%f \t %e \t %e \n", a, res_long, res_long*conv); fprintf(fp3, "%f \t %e \t %e \n", a, res_tt, res_tt*conv); fprintf(fp4, "%f \t %e \t %e \n", a, res_tl, res_tl*conv); fprintf(fp5, "%f \t %e \t %e \n", a, res_tlp, res_tlp*conv);*/ fprintf(fp1, "%f \t %e \t %e \t %e \t %e \n", a, res_u*conv, res_tl*conv, res_tt*conv, res_tlp*conv); } /*fclose(fp1); fclose(fp2); fclose(fp3); fclose(fp4); fclose(fp5);*/ fclose(fp1); break; } } break; } } break; } /***********************************************************************/ /* */ /* (e, e' K0 SIGMA+) reaction */ /* */ /********************************************************************* */ case 3: /* for (e, e' K0 Sigma+) reaction */ { /* 1. unpolarized forward differential cross section */ /* 2. unpolarized backward differential cross section */ int pol = 1; k0sig_s=1; switch(pol) { case 1: /* unpolarized (e, e' K0 Sigma+) differential cross section : forward */ { int mech, deg_k, deg_kstar, separated; double res_trans, res_long, res_tt, res_tl, res_tlp; /* printf("\n\n Give choice for FORWARD scattering mechanism"); printf("\n 1 = t-channel kaon"); printf("\n 2 = kstar"); printf("\n 3 = kaon gaugeinvariant"); printf("\n 4 = kaon gaugeinvariant + kstar"); printf("\n"); scanf("%d", &mech); printf("\n\n Give choice for degeneracy of KAON trajectory"); printf("\n 1 = non-degenerate trajectory"); printf("\n 2 = degenerate trajectory with rotating phase (default)"); printf("\n 3 = degenerate trajectory with NON-rotating phase"); printf("\n 4 = single pole"); printf("\n"); scanf("%d", °_k); */ printf("\n\n Give choice for degeneracy of KSTAR trajectory"); printf("\n 1 = non-degenerate trajectory"); printf("\n 2 = degenerate trajectory with rotating phase (default)"); printf("\n 3 = degenerate trajectory with NON-rotating phase"); printf("\n 4 = single pole"); printf("\n"); scanf("%d", °_kstar); mech=2; deg_k=0; switch(mech) { case 1 : { break; } case 2 : { if (deg_kstar == 2) strcpy(fileb,"_kstar_regge.dat"); if (deg_kstar == 4) strcpy(fileb,"_kstar_born.dat"); break; } case 3 : { if (deg_k == 2) strcpy(fileb,"_k_regge.dat"); if (deg_k == 4) strcpy(fileb,"_k_born.dat"); break; } case 4 : { if ((deg_k == 2) && (deg_kstar == 2)) strcpy(fileb,"_kkstar_regge.dat"); if ((deg_k == 4) && (deg_kstar == 4)) strcpy(fileb,"_kkstar_born.dat"); break; } } /* printf("\n\n Give choice for electromagnetic form factor of K"); printf("\n 1 = Monopole parametrization"); printf("\n 2 = Monopole at large Q_sqr, correct charge radius at low Q_sqr"); printf("\n 3 = VMD form (rho, omega, phi mesons)"); printf("\n 4 = FF with off-shell extrapolation"); printf("\n"); scanf("%d", &ichoicek); printf("\n\n Give choice for electromagnetic form factor of K*"); printf("\n 1 = Monopole parametrization"); printf("\n 2 = Monopole at large Q_sqr, correct charge radius at low Q_sqr"); printf("\n 3 = VMD form (rho, omega, phi mesons)"); printf("\n"); scanf("%d", &ichoicekst); if((ichoicek == 1)||(ichoicek == 2)) { printf("\n\n K form factor mass scale/cut-off"); printf("\n"); scanf("%lf", &mscalk); } if(ichoicek == 4) { printf("\n Give the value of the dimensionless parameter c in the kaon off-shell FF (e.g. 1) \n"); scanf("%lf",&c_offshell); } if((ichoicekst == 1)||(ichoicekst == 2)) { printf("\n\n K* form factor mass scale/cut-off"); printf("\n"); scanf("%lf", &mscalkst); } */ ichoicek=1; ichoicekst=1; /* printf("\n\n Give choice for OBSERVABLE to calculate"); printf("\n 1 = separated cross section : RESPONSE functions"); printf("\n 2 = theta dependence (Cornell experiments : Bebek et al.)"); printf("\n 3 = q2 dependence"); printf("\n 4 = recoil polarization"); printf("\n 7 = cross sections (T,L,TT,LT) : user-defined kinematics"); printf("\n"); scanf("%d", &separated); */ separated=1; switch(separated) { case 1 : /* separated cross section */ { char file1t[80]="ep_ekphyp_w2p15_qsqr1p19"; char file2t[80]; char trans[]="_trans_kkstar.dat"; char longi[]="_long_kkstar.dat"; char tt[]="_tt_pi_gaugeinv.dat"; char tl[]="_tl_pi_gaugeinv.dat"; char tlp[]="_tlp_k_gaugeinv.dat"; char file1l[80],file1tt[80],file1tl[80],file1tlp[80]; int photoelectro; strcpy(file1l,file1t); strcpy(file1tt,file1t); strcpy(file1tl,file1t); strcpy(file1tlp,file1t); /* printf("\n photo or electro ?"); printf("\n 1 = photo"); printf("\n 2 = electro"); printf("\n"); scanf("%d", &photoelectro); */ photoelectro=1; switch(photoelectro) { case 1 : /*photoproduction*/ { double ga_en_lab, bidon1, cos_kth, conv; char file1[80]; char file2[80]; int totdiff; printf("\n total or differential cross section ?"); printf("\n 1 = total"); printf("\n 2 = differential"); printf("\n"); scanf("%d", &totdiff); switch(totdiff) { case 1 : /*photoproduction : total cross section*/ { double tmin, tmax, result1; Q_sqr=0.; sprintf(file1,"gap_k0sp_tot.dat"); if ( (fp1 = fopen(file1,"w")) == NULL) fprintf(stderr,"\nUnable to open file"); printf("OUTPUT : Egamma (GeV), W(GeV), sig_tot(mubarn) \n"); for( ga_en_lab= 1.; ga_en_lab <= 4.; ga_en_lab += .25) { s=pow(NU_M,2.)+2.*ga_en_lab*NU_M; w=sqrt(s); tmin=pow((-Q_sqr-pow(K0_M,2.)-pow(NU_M,2.)+pow(SIG_M,2.))/2./sqrt(s),2.) -pow(sqrt(pow((s-Q_sqr-pow(NU_M,2.))/2./sqrt(s),2.)+Q_sqr) -sqrt(pow((s+pow(K0_M,2.)-pow(SIG_M,2.))/2./sqrt(s),2.)-pow(K0_M,2.)),2.); tmax=pow((-Q_sqr-pow(K0_M,2.)-pow(NU_M,2.)+pow(SIG_M,2.))/2./sqrt(s),2.) -pow(sqrt(pow((s-Q_sqr-pow(NU_M,2.))/2./sqrt(s),2.)+Q_sqr) +sqrt(pow((s+pow(K0_M,2.)-pow(SIG_M,2.))/2./sqrt(s),2.)-pow(K0_M,2.)),2.); result1 = sigma_trans_tot(reaction, mech, deg_k, deg_kstar, Q_sqr, s, tmin, tmax); fprintf(fp1,"%f \t %e \t %e \n", ga_en_lab, w, result1); } fclose(fp1); break; } case 2 : /*photoproduction : differential cross section*/ { double ga_en_lab, bidon1, cos_k_th; int wort; char file1[80]; char file2[80]; printf("\n As a function of t or W ?"); printf("\n 1: As a function of t"); printf("\n 2: As a function of W"); printf("\n"); scanf("%d", &wort); if (wort==1) { printf("\n photon energy (in GeV) ?"); printf("\n"); scanf("%lf", &ga_en_lab); bidon1=(ga_en_lab-(int) ga_en_lab +.01)*10.; sprintf(file1,"gap_k0sp_diff_%dp%dgev.dat", (int) ga_en_lab, (int) bidon1); sprintf(file2,"gap_k0sp_ass_%dp%dgev.dat", (int) ga_en_lab, (int) bidon1); s = 2. * NU_M * ga_en_lab + pow(NU_M,2.); Q_sqr = 0.; ga_mom_cm = sqrt( ( pow(s - pow(NU_M,2.) + Q_sqr, 2.) + pow(2. * NU_M, 2.) * Q_sqr ) / (4. * s) ); ga_en_cm = sqrt(s) - sqrt(pow(ga_mom_cm,2.) + pow(NU_M,2.)); k_mom_cm = sqrt( ( pow(s - pow(SIG_M,2.) - pow(K0_M,2.), 2.) - pow(2. * SIG_M * K0_M, 2.) ) / (4. * s) ); k_en = sqrt( pow(k_mom_cm, 2.) + pow(K0_M, 2.) ); if ( (fp1 = fopen(file1,"w")) == NULL) fprintf(stderr,"\nUnable to open file"); if ( (fp2 = fopen(file2,"w")) == NULL) fprintf(stderr,"\nUnable to open file"); printf("OUTPUT1 : -t(GeV2), cos(thcm), dsig_T/dt(mubarn/GeV2) \n"); printf("OUTPUT2 : -t(GeV2), cos(thcm), -dsig_TT/dsig_T_dt \n"); for(mt = 0.001; mt <= 1.5; mt += 0.02) { cos_k_th = (- mt - pow(K_M,2.) + 2. * ga_en_cm * k_en) / (2. * ga_mom_cm * k_mom_cm) ; res_trans = response_R_T(reaction, mech, deg_k, deg_kstar, Q_sqr, s, - mt); res_long = response_R_TT(reaction, mech, deg_k, deg_kstar, Q_sqr, s, - mt); res_long = res_long / res_trans; fprintf(fp1, "%f \t %f \t %e \n", mt, cos_k_th, res_trans); fprintf(fp2, "%f \t %f \t %e \n", mt, cos_k_th, res_long); } fclose(fp2); } if (wort==2) { double angle, conv; printf("\n Angle (in degrees) ?"); printf("\n"); scanf("%lf", &angle); bidon1=(angle-(int) angle +.01)*10.; sprintf(file1,"gap_k0sp_diff_%dp%ddeg.dat", (int) angle, (int) bidon1); Q_sqr = 0.; if ( (fp1 = fopen(file1,"w")) == NULL) fprintf(stderr,"\nUnable to open file"); printf("OUTPUT1 : W(GeV), dsig_T/dt(mubarn/GeV2), dsig_T/dcos(thcm) (mubarn)\n"); for(w = 1.69; w <= 2.4; w += 0.1) { s=pow(w,2.); ga_mom_cm = sqrt( ( pow(s - pow(NU_M,2.) + Q_sqr, 2.) + pow(2. * NU_M, 2.) * Q_sqr ) / (4. * s) ); ga_en_cm = sqrt(s) - sqrt(pow(ga_mom_cm,2.) + pow(NU_M,2.)); k_mom_cm = sqrt( ( pow(s - pow(SIG_M,2.) - pow(K_M,2.), 2.) - pow(2. * SIG_M * K_M, 2.) ) / (4. * s) ); k_en = sqrt( pow(k_mom_cm, 2.) + pow(K_M, 2.) ); mt = - 2. * ga_mom_cm * k_mom_cm *cos(angle*PI/180.) - pow(K_M,2.) + 2. * ga_en_cm * k_en ; conv = 2. * k_mom_cm * ga_mom_cm ; res_trans = response_R_T(reaction, mech, deg_k, deg_kstar, Q_sqr, s, - mt); res_long = response_R_TT(reaction, mech, deg_k, deg_kstar, Q_sqr, s, - mt); res_long = res_long / res_trans; fprintf(fp1, "%f \t %f \t %e \n", w, res_trans, conv * res_trans); } } fclose(fp1); break; } } break; } case 2 : /*electroproduction*/ { int sep; printf("\n Which kinematics ?"); printf("\n Ratio R_L/R_T : "); printf("\n 1 : W = 2.14 GeV, Q^2 = 1.19 GeV^2"); printf("\n 2 : W = 2.66 GeV, Q^2 = 1.99 GeV^2"); printf("\n 3 : W = 2.56 GeV, Q^2 = 3.38 GeV^2"); printf("\n R_T+epsilon R_L, R_TT, R_TL : "); printf("\n 4 : W = 2.21 GeV, Q^2 = .7 GeV^2, epsilon =.85 (Brauel et al.)"); printf("\n 5 : W = 2.21 GeV, Q^2 = 1.35 GeV^2, epsilon =.82 (Brauel et al.)"); printf("\n 6 : W = 2.21 GeV, Q^2 = .06 GeV^2, epsilon =.41 (Brauel et al.)"); printf("\n"); scanf("%d", &sep); switch(sep) { case 1 : /*ratio rlrt*/ { char file1[80]="ep_eks_w2p14_qsqr1p19_rap"; s = pow(2.14, 2.); Q_sqr = 1.19; if ( (fp1 = fopen(strcat(file1,fileb),"w")) == NULL) fprintf(stderr,"\nUnable to open file"); ga_mom_cm = sqrt( ( pow(s - pow(NU_M,2.) + Q_sqr, 2.) + pow(2. * NU_M, 2.) * Q_sqr ) / (4. * s) ); ga_en_cm = sqrt(s) - sqrt(pow(ga_mom_cm,2.) + pow(NU_M,2.)); k_mom_cm = sqrt( ( pow(s - pow(SIG_M,2.) - pow(K_M,2.), 2.) - pow(2. * SIG_M * K_M, 2.) ) / (4. * s) ); k_en = sqrt( pow(k_mom_cm, 2.) + pow(K_M, 2.) ); for(k_th = 0.001; k_th <= 21.; k_th += 2.) { t = - Q_sqr + pow(K_M,2.) - 2. * ga_en_cm * k_en + 2. * ga_mom_cm * k_mom_cm * cos(k_th * PI / 180.); res_trans = response_R_T(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t); res_long = response_R_L(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t); res_long = res_long / res_trans; fprintf(fp1, "%f \t %f \t %e \n", k_th, - t, res_long); } fclose(fp1); break; } case 2 : /*ratio rlrt*/ { char file1[80]="ep_eks_w2p66_qsqr1p99_rap"; s = pow(2.66, 2.); Q_sqr = 1.99; if ( (fp1 = fopen(strcat(file1,fileb),"w")) == NULL) fprintf(stderr,"\nUnable to open file"); ga_mom_cm = sqrt( ( pow(s - pow(NU_M,2.) + Q_sqr, 2.) + pow(2. * NU_M, 2.) * Q_sqr ) / (4. * s) ); ga_en_cm = sqrt(s) - sqrt(pow(ga_mom_cm,2.) + pow(NU_M,2.)); k_mom_cm = sqrt( ( pow(s - pow(SIG_M,2.) - pow(K_M,2.), 2.) - pow(2. * SIG_M * K_M, 2.) ) / (4. * s) ); k_en = sqrt( pow(k_mom_cm, 2.) + pow(K_M, 2.) ); for(k_th = 0.001; k_th <= 21.; k_th += 2.) { t = - Q_sqr + pow(K_M,2.) - 2. * ga_en_cm * k_en + 2. * ga_mom_cm * k_mom_cm * cos(k_th * PI / 180.); res_trans = response_R_T(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t); res_long = response_R_L(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t); res_long = res_long / res_trans; fprintf(fp1, "%f \t %f \t %e \n", k_th, - t, res_long); } fclose(fp1); break; } case 3 : /*ratio rlrt*/ { char file1[80]="ep_eks_w2p56_qsqr3p38_rap"; s = pow(2.56, 2.); Q_sqr = 3.38; if ( (fp1 = fopen(strcat(file1,fileb),"w")) == NULL) fprintf(stderr,"\nUnable to open file"); ga_mom_cm = sqrt( ( pow(s - pow(NU_M,2.) + Q_sqr, 2.) + pow(2. * NU_M, 2.) * Q_sqr ) / (4. * s) ); ga_en_cm = sqrt(s) - sqrt(pow(ga_mom_cm,2.) + pow(NU_M,2.)); k_mom_cm = sqrt( ( pow(s - pow(SIG_M,2.) - pow(K_M,2.), 2.) - pow(2. * SIG_M * K_M, 2.) ) / (4. * s) ); k_en = sqrt( pow(k_mom_cm, 2.) + pow(K_M, 2.) ); for(k_th = 0.001; k_th <= 21.; k_th += 2.) { t = - Q_sqr + pow(K_M,2.) - 2. * ga_en_cm * k_en + 2. * ga_mom_cm * k_mom_cm * cos(k_th * PI / 180.); res_trans = response_R_T(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t); res_long = response_R_L(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t); res_long = res_long / res_trans; fprintf(fp1, "%f \t %f \t %e \n", k_th, - t, res_long); } fclose(fp1); break; } case 4 : /*separation rt+eps*rl, rtt, rtl*/ { char file1[80]="ep_eks_w2p21_qsqr0p7_tpepsl"; char file2[80]="ep_eks_w2p21_qsqr0p7_tt"; char file3[80]="ep_eks_w2p21_qsqr0p7_tl"; s = pow(2.21, 2.); Q_sqr = .7; epsilon = .85; if ( (fp1 = fopen(strcat(file1,fileb),"w")) == NULL) fprintf(stderr,"\nUnable to open file"); if ( (fp2 = fopen(strcat(file2,fileb),"w")) == NULL) fprintf(stderr,"\nUnable to open file"); if ( (fp3 = fopen(strcat(file3,fileb),"w")) == NULL) fprintf(stderr,"\nUnable to open file"); ga_mom_cm = sqrt( ( pow(s - pow(NU_M,2.) + Q_sqr, 2.) + pow(2. * NU_M, 2.) * Q_sqr ) / (4. * s) ); ga_en_cm = sqrt(s) - sqrt(pow(ga_mom_cm,2.) + pow(NU_M,2.)); k_mom_cm = sqrt( ( pow(s - pow(SIG_M,2.) - pow(K_M,2.), 2.) - pow(2. * SIG_M * K_M, 2.) ) / (4. * s) ); k_en = sqrt( pow(k_mom_cm, 2.) + pow(K_M, 2.) ); for(mt = 0.01; mt <= .7; mt += 0.1) { res_trans = response_R_T(reaction, mech, deg_k, deg_kstar, Q_sqr, s, -mt); res_long = response_R_L(reaction, mech, deg_k, deg_kstar, Q_sqr, s, -mt); res_trans = res_trans+epsilon*res_long; res_tt = response_R_TT(reaction, mech, deg_k, deg_kstar, Q_sqr, s, - mt); res_tl = response_R_TL(reaction, mech, deg_k, deg_kstar, Q_sqr, s, - mt); /* printf("%f \t %e \n", mt, res_long);*/ fprintf(fp1, "%f \t %e \n", mt, res_trans); fprintf(fp2, "%f \t %e \n", mt, res_tt); fprintf(fp3, "%f \t %e \n", mt, res_tl); } fclose(fp1); fclose(fp2); fclose(fp3); break; } case 5 : /*separation rt+eps*rl, rtt, rtl*/ { char file1[80]="ep_eks_w2p21_qsqr1p35_tpepsl"; char file2[80]="ep_eks_w2p21_qsqr1p35_tt"; char file3[80]="ep_eks_w2p21_qsqr1p35_tl"; s = pow(2.21, 2.); Q_sqr = 1.35; epsilon = .82; if ( (fp1 = fopen(strcat(file1,fileb),"w")) == NULL) fprintf(stderr,"\nUnable to open file"); if ( (fp2 = fopen(strcat(file2,fileb),"w")) == NULL) fprintf(stderr,"\nUnable to open file"); if ( (fp3 = fopen(strcat(file3,fileb),"w")) == NULL) fprintf(stderr,"\nUnable to open file"); ga_mom_cm = sqrt( ( pow(s - pow(NU_M,2.) + Q_sqr, 2.) + pow(2. * NU_M, 2.) * Q_sqr ) / (4. * s) ); ga_en_cm = sqrt(s) - sqrt(pow(ga_mom_cm,2.) + pow(NU_M,2.)); k_mom_cm = sqrt( ( pow(s - pow(SIG_M,2.) - pow(K_M,2.), 2.) - pow(2. * SIG_M * K_M, 2.) ) / (4. * s) ); k_en = sqrt( pow(k_mom_cm, 2.) + pow(K_M, 2.) ); for(mt = 0.01; mt <= .7; mt += 0.1) { res_trans = response_R_T(reaction, mech, deg_k, deg_kstar, Q_sqr, s, -mt); res_long = response_R_L(reaction, mech, deg_k, deg_kstar, Q_sqr, s, -mt); res_trans = res_trans+epsilon*res_long; res_tt = response_R_TT(reaction, mech, deg_k, deg_kstar, Q_sqr, s, - mt); res_tl = response_R_TL(reaction, mech, deg_k, deg_kstar, Q_sqr, s, - mt); fprintf(fp1, "%f \t %e \n", mt, res_trans); fprintf(fp2, "%f \t %e \n", mt, res_tt); fprintf(fp3, "%f \t %e \n", mt, res_tl); } fclose(fp1); fclose(fp2); fclose(fp3); break; } case 6 : /*separation rt+eps*rl, rtt, rtl*/ { char file1[80]="ep_eks_w2p21_qsqr0p06_tpepsl"; char file2[80]="ep_eks_w2p21_qsqr0p06_tt"; char file3[80]="ep_eks_w2p21_qsqr0p06_tl"; s = pow(2.21, 2.); Q_sqr = .06; epsilon = .41; if ( (fp1 = fopen(strcat(file1,fileb),"w")) == NULL) fprintf(stderr,"\nUnable to open file"); if ( (fp2 = fopen(strcat(file2,fileb),"w")) == NULL) fprintf(stderr,"\nUnable to open file"); if ( (fp3 = fopen(strcat(file3,fileb),"w")) == NULL) fprintf(stderr,"\nUnable to open file"); ga_mom_cm = sqrt( ( pow(s - pow(NU_M,2.) + Q_sqr, 2.) + pow(2. * NU_M, 2.) * Q_sqr ) / (4. * s) ); ga_en_cm = sqrt(s) - sqrt(pow(ga_mom_cm,2.) + pow(NU_M,2.)); k_mom_cm = sqrt( ( pow(s - pow(SIG_M,2.) - pow(K_M,2.), 2.) - pow(2. * SIG_M * K_M, 2.) ) / (4. * s) ); k_en = sqrt( pow(k_mom_cm, 2.) + pow(K_M, 2.) ); for(mt = 0.01; mt <= .7; mt += 0.1) { res_trans = response_R_T(reaction, mech, deg_k, deg_kstar, Q_sqr, s, -mt); res_long = response_R_L(reaction, mech, deg_k, deg_kstar, Q_sqr, s, -mt); /* res_trans = res_trans+epsilon*res_long;*/ res_tt = response_R_TT(reaction, mech, deg_k, deg_kstar, Q_sqr, s, - mt); res_tl = response_R_TL(reaction, mech, deg_k, deg_kstar, Q_sqr, s, - mt); fprintf(fp1, "%f \t %e \n", mt, res_trans); fprintf(fp2, "%f \t %e \n", mt, res_tt); fprintf(fp3, "%f \t %e \n", mt, res_tl); } fclose(fp1); fclose(fp2); fclose(fp3); break; } } } } break; } case 2 : /* theta dependence (Cornell experiments : Bebek et al.) */ { int i, kin, choiceintegr; double cos_min,mt_min,lbin; printf("\n\n Choose the kinematics of the Cornell experiments (Bebek et al.)"); printf("\n 1 : W = 2.2 GeV, Q^2 = 1.18 GeV^2, epsilon = 0.94"); printf("\n 2 : W = 2.17 GeV, Q^2 = 1.97 GeV^2, epsilon = 0.94"); printf("\n 3 : W = 2.19 GeV, Q^2 = 3.95 GeV^2, epsilon = 0.88"); printf("\n 4 : W = 2.78 GeV, Q^2 = 1.38 GeV^2, epsilon = 0.88"); printf("\n 5 : W = 2.47 GeV, Q^2 = 3.52 GeV^2, epsilon = 0.86"); printf("\n 6 : W = 2.66 GeV, Q^2 = 2.04 GeV^2, epsilon = 0.88"); printf("\n 7 : W = 2.21 GeV, Q^2 = .06 GeV^2, epsilon = 0.41 (Brauel et al.)"); printf("\n 8 : W = 2.21 GeV, Q^2 = .28 GeV^2, epsilon = 0.74 (Brauel et al.)"); printf("\n 9 : W = 2.21 GeV, Q^2 = .7 GeV^2, epsilon = 0.85 (Brauel et al.)"); printf("\n 10 : W = 2.21 GeV, Q^2 = 1.35 GeV^2, epsilon = 0.82 (Brauel et al.)"); printf("\n"); printf("\n"); scanf("%d", &kin); /* printf("\n\n Integration on Phi (120=>240) or not?"); printf("\n 0 : no"); printf("\n 1 : yes"); printf("\n"); scanf("%d", &choiceintegr);*/ choiceintegr=1; switch(kin) { case 1 : { strcpy(file2,"ep_eks_w2p2_qsqr1p18_eps0p94"); w=2.2; Q_sqr=1.18; epsilon=.94; break; } case 2 : { strcpy(file2,"ep_eks_w2p17_qsqr1p97_eps0p94"); w=2.17; Q_sqr=1.97; epsilon=.94; break; } case 3 : { strcpy(file2,"ep_eks_w2p19_qsqr3p95_eps0p88"); w=2.19; Q_sqr=3.95; epsilon=.88; break; } case 4 : { strcpy(file2,"ep_eks_w2p78_qsqr1p38_eps0p88"); w=2.78; Q_sqr=1.38; epsilon=.88; break; } case 5 : { strcpy(file2,"ep_eks_w2p47_qsqr3p52_eps0p86"); w=2.47; Q_sqr=3.52; epsilon=.86; break; } case 6 : { strcpy(file2,"ep_eks_w2p66_qsqr2p04_eps0p88"); w=2.66; Q_sqr=2.04; epsilon=.88; /* w=2.21; Q_sqr=.06; epsilon=.41;*/ break; } case 7 : { strcpy(file2,"ep_eks_w2p21_qsqr0p06_eps0p41"); w=2.21; Q_sqr=.06; epsilon=.41; break; } case 8 : { strcpy(file2,"ep_eks_w2p21_qsqr0p28_eps0p74"); w=2.21; Q_sqr=.28; epsilon=.74; break; } case 9 : { strcpy(file2,"ep_eks_w2p21_qsqr0p70_eps0p85"); w=2.21; Q_sqr=.7; epsilon=.85; break; } case 10 : { strcpy(file2,"ep_eks_w2p21_qsqr1p35_eps0p82"); w=2.21; Q_sqr=1.35; epsilon=.82; break; } } if ( (fp2 = fopen(strcat(file2,fileb),"w")) == NULL) fprintf(stderr,"\nUnable to open file"); s = pow(w, 2.); ga_mom_cm = sqrt( ( pow(s - pow(NU_M,2.) + Q_sqr, 2.) + pow(2. * NU_M, 2.) * Q_sqr ) / (4. * s) ); ga_en_cm = sqrt(s) - sqrt(pow(ga_mom_cm,2.) + pow(NU_M,2.)); k_mom_cm = sqrt( ( pow(s - pow(SIG_M,2.) - pow(K_M,2.), 2.) - pow(2. * SIG_M * K_M, 2.) ) / (4. * s) ); k_en = sqrt( pow(k_mom_cm, 2.) + pow(K_M, 2.) ); if (kin<7){ for(k_th = 0.001; k_th <= 21.; k_th += 2.) /* for(k_th = 0.001; k_th <= 180.; k_th += 2.)*/ { double conv; t = - Q_sqr + pow(K_M,2.) - 2. * ga_en_cm * k_en + 2. * ga_mom_cm * k_mom_cm * cos(k_th * PI / 180.); conv = 2. * k_mom_cm * ga_mom_cm / (2.*PI) ; diff = conv * dsigma_dt_anal(reaction, mech, deg_k, deg_kstar, 0., epsilon, Q_sqr, s, t); /* fprintf(fp2, "%f \t %f \t %e \n", cos(k_th*PI/180.), -t, diff); */ fprintf(fp2, "%f \t %f \t %e \n", k_th, -t, diff); } fclose(fp2); } mt_min=0.; do{ mt_min=mt_min+.01; cos_min = (- Q_sqr + pow(K_M,2.) - 2. * ga_en_cm * k_en + mt_min)/ (-2. * ga_mom_cm * k_mom_cm) ; }while(fabs(cos_min)>1.); if ((kin>6)&&(kin<11)){ phi_min = 120. * PI / 180.; phi_max = 240. * PI / 180.; /* phi_min = 0. * PI / 180.; phi_max = 360. * PI / 180.;*/ lbin=phi_max-phi_min; for(mt = mt_min; mt <= 1.7; mt += 0.2) { if (choiceintegr==0) { diff = 2.*PI* dsigma_dtdphi(reaction, mech, deg_k, deg_kstar, 0., PI ,epsilon, Q_sqr, s, -mt); difft = 2.*PI* dsigma_dtdphi(reaction, mech, deg_k, deg_kstar, 0., PI ,0., Q_sqr, s, -mt); diffl = 2.*PI* dsigmal_dtdphi(reaction, mech, deg_k, deg_kstar, 0., PI ,epsilon, Q_sqr, s, -mt); } else { diff = 2. * PI * dsigma_dt(reaction, mech, deg_k, deg_kstar, 0., phi_min, phi_max, epsilon, Q_sqr, s, - mt); difft = 2. * PI * dsigma_dt(reaction, mech, deg_k, deg_kstar, 0., phi_min, phi_max, 0. , Q_sqr, s, - mt); diffl = 2. * PI * dsigmal_dt(reaction, mech, deg_k, deg_kstar, 0., phi_min, phi_max, epsilon , Q_sqr, s, - mt); difftt = 2. * PI * dsigmatt_dt(reaction, mech, deg_k, deg_kstar, 0., phi_min, phi_max, epsilon , Q_sqr, s, - mt); difftl = 2. * PI * dsigmatl_dt(reaction, mech, deg_k, deg_kstar, 0., phi_min, phi_max, epsilon , Q_sqr, s, - mt); diff=diff/lbin; difft=difft/lbin; diffl=diffl/lbin; difftt=difftt/lbin; difftl=difftl/lbin; } if (diff!=0) fprintf(fp2, "%f \t %e \t %e \t %e \t %e \t %e \n", mt, diff, difft, diffl, difftt, difftl); } fclose(fp2); } break; } case 3 : /* q2 dependence */ { int obsq2; printf("\n\n Which observable ?"); printf("\n 1 = dsig/dt (Bebek et al.)"); printf("\n 2 = dsigT/dt, dsigL/dt, sigL/sigT (CEBAF)"); printf("\n 3 = dsigT/dt, dsigL/dt, sigL/sigT (Bydzo.)"); printf("\n"); scanf("%d", &obsq2); switch(obsq2) { case 1: /* dsigdt */ { int kin; double diff_test; char filea[80]; printf("\n\n Choose the angle"); printf("\n 1 : K Sigma electroproduction : Q2 dependence"); printf("\n : W = 2.17 GeV, theta = 7 deg, epsilon = 0.85"); printf("\n 2 : K Sigma electroproduction : Q2 dependance"); printf("\n : W = 2.17 GeV, theta = 1 deg, epsilon = 0.85"); printf("\n 3 : K Sigma electroproduction : Q2 dependance"); printf("\n : W = 2.17 GeV, theta = 15 deg, epsilon = 0.85"); printf("\n 4 : K Sigma electroproduction : Q2 dependance"); printf("\n : W = 2.66 GeV, theta = 7 deg, epsilon = 0.85"); printf("\n"); scanf("%d", &kin); if (kin==1) { strcpy(filea,"ep_eks_w2p17_the7p0_eps0p85_q2"); k_th=7.; } if (kin==2) { strcpy(filea,"ep_eks_w2p17_the1p0_eps0p85_q2"); k_th=1.; } if (kin==3) { strcpy(filea,"ep_eks_w2p17_the15p0_eps0p85_q2"); k_th=15.; } if (kin==4) { strcpy(filea,"ep_eks_w2p66_the8p0_eps0p85_q2"); k_th=8.; } if (kin<=3) { s = pow(2.17, 2.); } if (kin==4) { s = pow(2.66, 2.); } epsilon = 0.85; if ( (fp1 = fopen(strcat(filea,fileb),"w")) == NULL) fprintf(stderr,"\nUnable to open file"); for(Q_sqr = 0.001; Q_sqr <= 5.; Q_sqr += 0.1) { double conv; ga_mom_cm = sqrt( ( pow(s - pow(NU_M,2.) + Q_sqr, 2.) + pow(2. * NU_M, 2.) * Q_sqr ) / (4. * s) ); ga_en_cm = sqrt(s) - sqrt(pow(ga_mom_cm,2.) + pow(NU_M,2.)); k_mom_cm = sqrt( ( pow(s - pow(SIG_M,2.) - pow(K_M,2.), 2.) - pow(2. * SIG_M * K_M, 2.) ) / (4. * s) ); k_en = sqrt( pow(k_mom_cm, 2.) + pow(K_M, 2.) ); t = - Q_sqr + pow(K_M,2.) - 2. * ga_en_cm * k_en + 2. * ga_mom_cm * k_mom_cm * cos(k_th * PI / 180.); conv = 2. * k_mom_cm * ga_mom_cm / (2. * PI); diff = conv * dsigma_dt_anal(reaction, mech, deg_k, deg_kstar, 0., epsilon, Q_sqr, s, t); diffl = conv * response_R_L(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t); fprintf(fp1, "%f \t %e \t %e \n", Q_sqr, diff, diffl); } fclose(fp1); break; } case 2 : /*ratio siglsigt*/ { char file1[80]="ep_eks_cebaf_rap.dat"; char file2[80]="ep_eks_cebaf_seplt.dat"; double rap, conv, cos_min; int iii; if ( (fp2 = fopen(file2,"w")) == NULL) fprintf(stderr,"\nUnable to open file"); for(iii = 1; iii <= 6; iii += 1) { if (iii==1) { s = pow(1.84, 2.); Q_sqr = .1; t = -.15; } if (iii==2) { s = pow(1.84, 2.); Q_sqr = .52; t = -.219; } if (iii==3) { s = pow(1.83, 2.); Q_sqr = .75; t = -.3; } if (iii==4) { s = pow(1.81, 2.); Q_sqr = 1.; t = -.407; } if (iii==5) { s = pow(1.84, 2.); Q_sqr = 2.; t = -.741; } if (iii==6) { s = pow(1.84, 2.); Q_sqr = 3.9; t = -1.4; } ga_mom_cm = sqrt( ( pow(s - pow(NU_M,2.) + Q_sqr, 2.) + pow(2. * NU_M, 2.) * Q_sqr ) / (4. * s) ); k_mom_cm = sqrt( (pow(s - pow(SIG_M,2.) - pow(K_M,2.), 2.) - pow(2. * SIG_M * K_M,2.) ) /(4. * s)); ga_en_cm = sqrt(s) - sqrt(pow(ga_mom_cm,2.)+pow(NU_M,2.)); k_en = sqrt( pow(k_mom_cm, 2.) + pow(K_M, 2.) ); conv = 2. * k_mom_cm * ga_mom_cm / (2.*PI) ; t=t+.01; do{ t=t-.01; cos_min = (- Q_sqr + pow(K_M,2.) - 2. * ga_en_cm * k_en -t)/(-2. * ga_mom_cm * k_mom_cm) ; }while(fabs(cos_min)>1.); res_trans = response_R_T(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t); res_long = response_R_L(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t); rap = res_long / res_trans; res_trans = res_trans*conv; res_long = res_long*conv; fprintf(fp2, "%f \t %f \t %f \t %f \t %f \n", t, Q_sqr, res_trans, res_long, rap); } fclose(fp2); break; } case 3: /* dsigdt */ { int kin; double diff_test; char filea[80]; printf("\n\n Choose the kinematics"); printf("\n 1 : K Sigma electroproduction : Q2 dependence"); printf("\n : W = 2.24 GeV, t = -.15, epsilon = 0.72"); printf("\n"); scanf("%d", &kin); if (kin==1) { strcpy(filea,"ep_eks_w2p24_t0p15_eps0p72_q2"); t=-.15; s = pow(2.24, 2.); epsilon = 0.72; } if ( (fp1 = fopen(strcat(filea,fileb),"w")) == NULL) fprintf(stderr,"\nUnable to open file"); for(Q_sqr = 0.001; Q_sqr <= 1.; Q_sqr += 0.05) { double conv; ga_mom_cm = sqrt( ( pow(s - pow(NU_M,2.) + Q_sqr, 2.) + pow(2. * NU_M, 2.) * Q_sqr ) / (4. * s) ); ga_en_cm = sqrt(s) - sqrt(pow(ga_mom_cm,2.) + pow(NU_M,2.)); k_mom_cm = sqrt( ( pow(s - pow(SIG_M,2.) - pow(K_M,2.), 2.) - pow(2. * SIG_M * K_M, 2.) ) / (4. * s) ); k_en = sqrt( pow(k_mom_cm, 2.) + pow(K_M, 2.) ); /* conv = 2. * k_mom_cm * ga_mom_cm / (2. * PI);*/ conv = 1.; diff = conv * dsigma_dt_anal(reaction, mech, deg_k, deg_kstar, 0.,epsilon,Q_sqr,s,t); difft = conv * response_R_T(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t); diffl = conv * response_R_L(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t); fprintf(fp1, "%f \t %e \t %e \t %e \n", Q_sqr, diff, difft,diffl); } fclose(fp1); break; } } break; } case 4 : /* recoil polarization */ { double asymm_lam, asymm_lam2, tortho, lortho, res_trans, res_long, W, epsilon ; printf("\n Which W ?"); printf("\n"); scanf("%lf", &W); printf("\n Which Q2 ?"); printf("\n"); scanf("%lf", &Q_sqr); printf("\n Which epsilon ?"); printf("\n"); scanf("%lf", &epsilon); if ( (fp1 = fopen("recoilsig.dat", "w")) == NULL) fprintf(stderr,"\nUnable to open file"); s=pow(W,2.); ga_mom_cm = sqrt( ( pow(s - pow(NU_M,2.) + Q_sqr, 2.) + pow(2. * NU_M, 2.) * Q_sqr ) / (4. * s) ); ga_en_cm = sqrt(s) - sqrt(pow(ga_mom_cm,2.) + pow(NU_M,2.)); k_mom_cm = sqrt( ( pow(s - pow(SIG_M,2.) - pow(K_M,2.), 2.) - pow(2. * SIG_M * K_M, 2.) ) / (4. * s) ); k_en = sqrt( pow(k_mom_cm, 2.) + pow(K_M, 2.) ); for(k_th = 0.; k_th <= 180.; k_th += 5.) { t = - Q_sqr + pow(K_M,2.) - 2. * ga_en_cm * k_en + 2. * ga_mom_cm * k_mom_cm * cos(k_th * PI / 180.); asymm_lam2 = recoil_asymm(reaction,mech,deg_k,deg_kstar, Q_sqr, s, t); tortho = response_R_T_ortho(reaction,mech,deg_k,deg_kstar, Q_sqr, s, t); res_trans = response_R_T(reaction,mech,deg_k,deg_kstar, Q_sqr, s, t); lortho = response_R_L_ortho(reaction,mech,deg_k,deg_kstar, Q_sqr, s, t); res_long = response_R_L(reaction,mech,deg_k,deg_kstar, Q_sqr, s, t); if (Q_sqr==0.) { asymm_lam = -tortho/res_trans; lortho=0.; res_long=0.; } else { asymm_lam = (tortho+epsilon*lortho) / (res_trans+epsilon*res_long); } fprintf(fp1,"%f \t %f \t %f \t %f \t %f \t %f \t %f \t %f\n", -t, asymm_lam, asymm_lam2, tortho, lortho, res_trans, res_long, k_th ); } fclose(fp); break; } case 7 : /* User-defined */ { double a, a_min, a_max, a_step, W_min, W_max, W_step, t_min, t_max, t_step, Q2_min, Q2_max, Q2_step, bidon1, bidon2; double ga_en_in_lab, ga_mom_in_lab, pi_mom_cm, W, tmin_test, mes_mom_lab_test, mes_en_lab, mes_mom_lab; int rep; char file1t[80],file1l[80],file1tt[80],file1tl[80]; printf("\n 1 : Q2 dependence"); printf("\n 2 : t dependence"); printf("\n 3 : W dependence"); printf("\n"); scanf("%d", &rep); if (rep==1) { printf("\n Which W"); printf("\n"); scanf("%lf", &W); printf("\n Which t (t>0) -if t=0, t=t_min-"); printf("\n"); scanf("%lf", &t); s = pow(W,2.); mt=-t; printf("\n Q2_min ?"); printf("\n"); scanf("%lf", &Q2_min); printf("\n Which Q2_max ?"); printf("\n"); scanf("%lf", &Q2_max); printf("\n Which Q2_step ?"); printf("\n"); scanf("%lf", &Q2_step); a_min=Q2_min; a_max=Q2_max; a_step=Q2_step; bidon1=(W-(int) W)*10.; bidon2=(t-(int) t)*10.; sprintf(file1t,"ep_ksigma_t%dp%d_w%dp%d_trans", (int) t, (int) bidon2, (int) W, (int) bidon1); sprintf(file1l,"ep_ksigma_t%dp%d_w%dp%d_long", (int) t, (int) bidon2, (int) W, (int) bidon1); sprintf(file1tt,"ep_ksigma_t%dp%d_w%dp%d_tt", (int) t, (int) bidon2, (int) W, (int) bidon1); sprintf(file1tl,"ep_ksigma_t%dp%d_w%dp%d_tl", (int) t, (int) bidon2, (int) W, (int) bidon1); printf("OUTPUT1 : Q2 (GeV2), dsig_T/dt (mubarn/GeV2), dsig_T/dOmega_cm (mubarn/sr) \n"); printf("OUTPUT2 : Q2 (GeV2), dsig_L/dt (mubarn/GeV2), dsig_L/dOmega_cm (mubarn/sr) \n"); printf("OUTPUT3 : Q2 (GeV2), dsig_TT/dt (mubarn/GeV2), dsig_TT/dOmega_cm (mubarn/sr) \n"); printf("OUTPUT4 : Q2 (GeV2), dsig_TL/dt (mubarn/GeV2), dsig_TL/dOmega_cm (mubarn/sr) \n"); } if (rep==2) { printf("\n Which W"); printf("\n"); scanf("%lf", &W); printf("\n Which Q2"); printf("\n"); scanf("%lf", &Q_sqr); s = pow(W,2.); printf("\n t_min ? (if 0 => kinematical t_min)"); printf("\n"); scanf("%lf", &t_min); printf("\n t_max ?"); printf("\n"); scanf("%lf", &t_max); printf("\n t_step ?"); printf("\n"); scanf("%lf", &t_step); a_min=t_min; a_max=t_max; a_step=t_step; if(t_min==0.) { a_min=pow((-Q_sqr-pow(K_M,2.)-pow(NU_M,2.)+pow(SIG_M,2.))/2./sqrt(s),2.) -pow(sqrt(pow((s-Q_sqr-pow(NU_M,2.))/2./sqrt(s),2.)+Q_sqr) -sqrt(pow((s+pow(K_M,2.)-pow(SIG_M,2.))/2./sqrt(s),2.)-pow(K_M,2.)),2.); a_min=-a_min; } bidon1=(Q_sqr-(int) Q_sqr)*10.; bidon2=(W-(int) W)*10.; sprintf(file1t,"ep_ksigma_q2%dp%d_w%dp%d_trans", (int) Q_sqr, (int) bidon1, (int) W, (int) bidon2); sprintf(file1l,"ep_ksigma_q2%dp%d_w%dp%d_long", (int) Q_sqr, (int) bidon1, (int) W, (int) bidon2); sprintf(file1tt,"ep_ksigma_q2%dp%d_w%dp%d_tt", (int) Q_sqr, (int) bidon1, (int) W, (int) bidon2); sprintf(file1tl,"ep_ksigma_q2%dp%d_w%dp%d_tl", (int) Q_sqr, (int) bidon1, (int) W, (int) bidon2); printf("OUTPUT1 : -t (GeV2), dsig_T/dt (mubarn/GeV2), dsig_T/dOmega_cm (mubarn/sr) \n"); printf("OUTPUT2 : -t (GeV2), dsig_L/dt (mubarn/GeV2), dsig_L/dOmega_cm (mubarn/sr) \n"); printf("OUTPUT3 : -t (GeV2), dsig_TT/dt (mubarn/GeV2), dsig_TT/dOmega_cm (mubarn/sr) \n"); printf("OUTPUT4 : -t (GeV2), dsig_TL/dt, (mubarn/GeV2),dsig_TL/dOmega_cm (mubarn/sr) \n"); } if (rep==3) { printf("\n Which Q2"); printf("\n"); scanf("%lf", &Q_sqr); printf("\n Which t (t>0) -if t=0, t=t_min-"); printf("\n"); scanf("%lf", &t); mt=-t; printf("\n W_min ?"); printf("\n"); scanf("%lf", &W_min); printf("\n W_max ?"); printf("\n"); scanf("%lf", &W_max); printf("\n W_step ?"); printf("\n"); scanf("%lf", &W_step); a_min=W_min; a_max=W_max; a_step=W_step; bidon1=(Q_sqr-(int) Q_sqr)*10.; bidon2=(t-(int) t)*10.; sprintf(file1t,"ep_ksigma_q2%dp%d_t%dp%d_trans", (int) Q_sqr, (int) bidon1, (int) t, (int) bidon2); sprintf(file1l,"ep_ksigma_q2%dp%d_t%dp%d_long", (int) Q_sqr, (int) bidon1, (int) t, (int) bidon2); sprintf(file1tt,"ep_ksigma_q2%dp%d_t%dp%d_tt", (int) Q_sqr, (int) bidon1, (int) t, (int) bidon2); sprintf(file1tl,"ep_ksigma_q2%dp%d_t%dp%d_tl", (int) Q_sqr, (int) bidon1, (int) t, (int) bidon2); printf("OUTPUT1 : W (GeV), dsig_T/dt (mubarn/GeV2), dsig_T/dOmega_cm (mubarn/sr) \n"); printf("OUTPUT2 : W (GeV), dsig_L/dt (mubarn/GeV2)dsig_L/dOmega_cm (mubarn/sr) \n"); printf("OUTPUT3 : W (GeV), dsig_TT/dt (mubarn/GeV2)dsig_TT/dOmega_cm (mubarn/sr) \n"); printf("OUTPUT4 : W (GeV), dsig_TL/dt (mubarn/GeV2)dsig_TL/dOmega_cm (mubarn/sr) \n"); } if ( (fp1 = fopen(strcat(file1t,fileb),"w")) == NULL) fprintf(stderr,"\nUnable to open file"); if ( (fp2 = fopen(strcat(file1l,fileb),"w")) == NULL) fprintf(stderr,"\nUnable to open file"); if ( (fp3 = fopen(strcat(file1tt,fileb),"w")) == NULL) fprintf(stderr,"\nUnable to open file"); if ( (fp4 = fopen(strcat(file1tl,fileb),"w")) == NULL) fprintf(stderr,"\nUnable to open file"); for(a = a_min; a <= a_max; a += a_step) { double conv; if (rep==1) Q_sqr=a; if (rep==2) mt=-a; if (rep==3) s=pow(a,2.); if((t==0)&&(rep!=2)) { mt=pow((-Q_sqr-pow(K_M,2.)-pow(NU_M,2.)+pow(SIG_M,2.))/2./sqrt(s),2.) -pow(sqrt(pow((s-Q_sqr-pow(NU_M,2.))/2./sqrt(s),2.)+Q_sqr) -sqrt(pow((s+pow(K_M,2.)-pow(SIG_M,2.))/2./sqrt(s),2.)-pow(K_M,2.)),2.); } ga_mom_cm = sqrt( ( pow(s - pow(NU_M,2.) + Q_sqr, 2.) + pow(2. * NU_M, 2.) * Q_sqr ) / (4. * s) ); pi_mom_cm = sqrt( ( pow(s - pow(SIG_M,2.) - pow(K_M,2.), 2.) - pow(2. * SIG_M * K_M, 2.) ) / (4. * s) ); conv = 1. / (2. * PI) * 2. * pi_mom_cm * ga_mom_cm; res_trans = response_R_T(reaction, mech, deg_k, deg_kstar, Q_sqr, s, mt); res_tt = response_R_TT(reaction, mech, deg_k, deg_kstar, Q_sqr, s, mt); if (Q_sqr==0.) { res_long=0.; res_tl=0.; } else { res_long = response_R_L(reaction, mech, deg_k, deg_kstar, Q_sqr, s, mt); res_tl = response_R_TL(reaction, mech, deg_k, deg_kstar, Q_sqr, s, mt); } fprintf(fp1, "%f \t %e \t %e \n", a, res_trans, res_trans*conv); fprintf(fp2, "%f \t %e \t %e \n", a, res_long, res_long*conv); fprintf(fp3, "%f \t %e \t %e \n", a, res_tt, res_tt*conv); fprintf(fp4, "%f \t %e \t %e \n", a, res_tl, res_tl*conv); } fclose(fp1); fclose(fp2); fclose(fp3); fclose(fp4); break; } } break; } } break; } /***********************************************************************/ /* */ /* STRANGENESS FORM FACTORS */ /* */ /***********************************************************************/ case 4: { double ff1, ff2; for(Q_sqr = 0.; Q_sqr <= 2.; Q_sqr += 0.01) { ichoicek = 1; mscalk = 1.5; ff1 = kaon_emformfactor(Q_sqr, pow(K_M,2.)); ichoicek = 3; ff2 = kaon_emformfactor(Q_sqr, pow(K_M,2.)); printf("\n%e \t %e \t %e", Q_sqr, ff1, ff2); } break; } } } /*************************************************************************/ /*************************************************************************/ /* */ /* e + N -> e + K + HYP : CROSS SECTIONS and ASYMMETRIES */ /* */ /*************************************************************************/ /*************************************************************************/ /* The function dsigma_dt() gives the e + N -> e + K + HYP differential cross section dsigma/dt integrated over Phi as function of electron helicity : el_hel, epsilon, Q_sqr, s and t. The units of the cross section are (microbarns / GeV^2). */ double dsigma_dt(int reaction, int mech, int deg_k, int deg_kstar, double el_hel, double phi_min, double phi_max, double epsilon, double Q_sqr, double s, double t) { reaction_s = reaction; mech_s = mech; deg_k_s = deg_k; deg_kstar_s = deg_kstar; el_hel_s = el_hel; epsilon_s = epsilon; Q_sqr_s = Q_sqr; s_s = s; t_s = t; return int1d(int_dsigma_dt, phi_min, phi_max, 5); } double dsigmal_dt(int reaction, int mech, int deg_k, int deg_kstar, double el_hel, double phi_min, double phi_max, double epsilon, double Q_sqr, double s, double t) { reaction_s = reaction; mech_s = mech; deg_k_s = deg_k; deg_kstar_s = deg_kstar; el_hel_s = el_hel; epsilon_s = epsilon; Q_sqr_s = Q_sqr; s_s = s; t_s = t; return int1d(int_dsigmal_dt, phi_min, phi_max, 5); } double dsigmatt_dt(int reaction, int mech, int deg_k, int deg_kstar, double el_hel, double phi_min, double phi_max, double epsilon, double Q_sqr, double s, double t) { reaction_s = reaction; mech_s = mech; deg_k_s = deg_k; deg_kstar_s = deg_kstar; el_hel_s = el_hel; epsilon_s = epsilon; Q_sqr_s = Q_sqr; s_s = s; t_s = t; return int1d(int_dsigmatt_dt, phi_min, phi_max, 5); } double dsigmatl_dt(int reaction, int mech, int deg_k, int deg_kstar, double el_hel, double phi_min, double phi_max, double epsilon, double Q_sqr, double s, double t) { reaction_s = reaction; mech_s = mech; deg_k_s = deg_k; deg_kstar_s = deg_kstar; el_hel_s = el_hel; epsilon_s = epsilon; Q_sqr_s = Q_sqr; s_s = s; t_s = t; return int1d(int_dsigmatl_dt, phi_min, phi_max, 5); } double int_dsigma_dt(double phi) { return dsigma_dtdphi(reaction_s, mech_s, deg_k_s, deg_kstar_s, el_hel_s, phi, epsilon_s, Q_sqr_s, s_s, t_s); } double int_dsigmal_dt(double phi) { return dsigmal_dtdphi(reaction_s, mech_s, deg_k_s, deg_kstar_s, el_hel_s, phi, epsilon_s, Q_sqr_s, s_s, t_s); } double int_dsigmatt_dt(double phi) { return dsigmatt_dtdphi(reaction_s, mech_s, deg_k_s, deg_kstar_s, el_hel_s, phi, epsilon_s, Q_sqr_s, s_s, t_s); } double int_dsigmatl_dt(double phi) { return dsigmatl_dtdphi(reaction_s, mech_s, deg_k_s, deg_kstar_s, el_hel_s, phi, epsilon_s, Q_sqr_s, s_s, t_s); } /* The function dsigma_dtdphi() gives the e + N -> e + K + HYP differential cross section dsigma/dt dPhi as function of electron helicity : el_hel, out-of-plane angle Phi, epsilon, Q_sqr, s and t. The units of the cross section are (microbarns / GeV^2). */ double dsigma_dtdphi(int reaction, int mech, int deg_k, int deg_kstar, double el_hel, double phi, double epsilon, double Q_sqr, double s, double t) { double unit_conv; unit_conv = 1. / (2. * PI); return unit_conv * ( response_R_T(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t) + epsilon * response_R_L(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t) - epsilon * cos(2. * phi) * response_R_TT(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t) + sqrt(2. * epsilon * (1. + epsilon)) * cos(phi) * response_R_TL(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t) + el_hel * sqrt(2. * epsilon * (1. - epsilon)) * sin(phi) * response_R_TLp(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t) ); } double dsigmal_dtdphi(int reaction, int mech, int deg_k, int deg_kstar, double el_hel, double phi, double epsilon, double Q_sqr, double s, double t) { double unit_conv; unit_conv = 1. / (2. * PI); return unit_conv * epsilon * response_R_L(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t); } double dsigmatt_dtdphi(int reaction, int mech, int deg_k, int deg_kstar, double el_hel, double phi, double epsilon, double Q_sqr, double s, double t) { double unit_conv; unit_conv = 1. / (2. * PI); return unit_conv * (- epsilon * cos(2. * phi) * response_R_TT(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t) ); } double dsigmatl_dtdphi(int reaction, int mech, int deg_k, int deg_kstar, double el_hel, double phi, double epsilon, double Q_sqr, double s, double t) { double unit_conv; unit_conv = 1. / (2. * PI); return unit_conv * (sqrt(2. * epsilon * (1. + epsilon)) * cos(phi) * response_R_TL(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t) ); } double dsigma_dt_anal(int reaction, int mech, int deg_k, int deg_kstar, double el_hel, double epsilon, double Q_sqr, double s, double t) { double unit_conv; unit_conv = 1. ; return unit_conv * ( response_R_T(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t) + epsilon * response_R_L(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t) + el_hel * response_R_TLp(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t) ); } double sigma_trans_tot(int reaction, int mech, int deg_k, int deg_kstar, double Q_sqr, double s, double tmin, double tmax) { double result; reaction_s = reaction; mech_s = mech; deg_k_s = deg_k; deg_kstar_s = deg_kstar; Q_sqr_s = Q_sqr; s_s = s; result = int_gauss(int_sigma_trans_tot, tmax, tmin, 50); return result; } double int_sigma_trans_tot(double t) { return response_R_T(reaction_s, mech_s, deg_k_s, deg_kstar_s, Q_sqr_s, s_s, t); } double pol_tot(int reaction, int mech, int deg_k, int deg_kstar, double Q_sqr, double s, double epsilon, double phi, int iintegr, int polbeam, int xyz) { double result,tmin,tmax; reaction_s = reaction; mech_s = mech; deg_k_s = deg_k; deg_kstar_s = deg_kstar; Q_sqr_s = Q_sqr; s_s = s; epsilon_s = epsilon; phi_s = phi; iintegr_s = iintegr; polbeam_s = polbeam; xyz_s = xyz; if(((pow((s-Q_sqr-pow(NU_M,2.))/2./sqrt(s),2.)+Q_sqr)>0.) &&((pow((s+pow(K_M,2.)-pow(LAM_M,2.))/2./sqrt(s),2.)-pow(K_M,2.))>0.)) { tmin=pow((-Q_sqr-pow(K_M,2.)-pow(NU_M,2.)+pow(LAM_M,2.))/2./sqrt(s),2.) -pow(sqrt(pow((s-Q_sqr-pow(NU_M,2.))/2./sqrt(s),2.)+Q_sqr) -sqrt(pow((s+pow(K_M,2.)-pow(LAM_M,2.))/2./sqrt(s),2.)-pow(K_M,2.)),2.); tmax=pow((-Q_sqr-pow(K_M,2.)-pow(NU_M,2.)+pow(LAM_M,2.))/2./sqrt(s),2.) -pow(sqrt(pow((s-Q_sqr-pow(NU_M,2.))/2./sqrt(s),2.)+Q_sqr) +sqrt(pow((s+pow(K_M,2.)-pow(LAM_M,2.))/2./sqrt(s),2.)-pow(K_M,2.)),2.); result = int_gauss(int_pol, tmax, tmin, 30); } else { result=0.; } return result; } double int_pol(double t) { double result,el_hel=1,den,num_par,num_long,num_perp, asymm_par,asymm_long,asymm_perp,asymm_par2,asymm_long2,asymm_perp2, asymm_par_prime,asymm_long_prime,asymm_perp_prime, asymm_par_prime2,asymm_long_prime2,asymm_perp_prime2, ga_mom_cm, ga_en_cm, k_mom_cm, k_en_cm, k_th, cos_k_th; if(iintegr_s==2){ den = response_R_T(reaction_s, mech_s, deg_k_s, deg_kstar_s, Q_sqr_s, s_s, t) + epsilon_s * response_R_L(reaction_s,mech_s, deg_k_s, deg_kstar_s, Q_sqr_s, s_s, t) - epsilon_s * cos(2. * phi_s) * response_R_TT(reaction_s,mech_s, deg_k_s, deg_kstar_s, Q_sqr_s, s_s, t) + sqrt(2. * epsilon_s * (1. + epsilon_s)) * cos(phi_s) * response_R_TL(reaction_s, mech_s, deg_k_s, deg_kstar_s, Q_sqr_s, s_s, t); if(polbeam_s==2) num_par = -sqrt(2. * epsilon_s * (1. + epsilon_s)) * sin(phi_s) * response_R_TL_par(reaction_s, mech_s, deg_k_s, deg_kstar_s, Q_sqr_s, s_s, t) - epsilon_s * sin(2. * phi_s) * response_R_TT_par(reaction_s, mech_s, deg_k_s, deg_kstar_s, Q_sqr_s, s_s, t) + el_hel * ( -sqrt(2. * epsilon_s * (1. - epsilon_s)) * cos(phi_s) * response_R_TLp_par(reaction_s, mech_s, deg_k_s, deg_kstar_s, Q_sqr_s, s_s, t) + sqrt(1.-pow(epsilon_s,2.)) * response_R_Tp_par(reaction_s, mech_s, deg_k_s, deg_kstar_s, Q_sqr_s, s_s, t)); if(polbeam_s==1) num_par = el_hel * ( -sqrt(2. * epsilon_s * (1. - epsilon_s)) * cos(phi_s) * response_R_TLp_par(reaction_s, mech_s, deg_k_s, deg_kstar_s, Q_sqr_s, s_s, t) + sqrt(1.-pow(epsilon_s,2.)) * response_R_Tp_par(reaction_s, mech_s, deg_k_s, deg_kstar_s, Q_sqr_s, s_s, t)); if(polbeam_s==2) num_long = -sqrt(2. * epsilon_s * (1. + epsilon_s)) * sin(phi_s) * response_R_TL_long(reaction_s, mech_s, deg_k_s, deg_kstar_s, Q_sqr_s, s_s, t) - epsilon_s * sin(2. * phi_s) * response_R_TT_long(reaction_s, mech_s, deg_k_s, deg_kstar_s, Q_sqr_s, s_s, t) + el_hel * ( -sqrt(2. * epsilon_s * (1. - epsilon_s)) * cos(phi_s) * response_R_TLp_long(reaction_s, mech_s, deg_k_s, deg_kstar_s, Q_sqr_s, s_s, t) + sqrt(1.-pow(epsilon_s,2.)) * response_R_Tp_long(reaction_s, mech_s, deg_k_s, deg_kstar_s, Q_sqr_s, s_s, t)); if(polbeam_s==1) num_long = el_hel * ( -sqrt(2. * epsilon_s * (1. - epsilon_s)) * cos(phi_s) * response_R_TLp_long(reaction_s, mech_s, deg_k_s, deg_kstar_s, Q_sqr_s, s_s, t) + sqrt(1.-pow(epsilon_s,2.)) * response_R_Tp_long(reaction_s, mech_s, deg_k_s, deg_kstar_s, Q_sqr_s, s_s, t)); if(polbeam_s==2) num_perp = response_R_T_ortho(reaction_s, mech_s, deg_k_s, deg_kstar_s, Q_sqr_s, s_s, t) +epsilon_s * response_R_L_ortho(reaction_s, mech_s, deg_k_s, deg_kstar_s, Q_sqr_s, s_s, t) -sqrt(2. * epsilon_s * (1. + epsilon_s)) * cos(phi_s) * response_R_TL_ortho(reaction_s, mech_s, deg_k_s, deg_kstar_s, Q_sqr_s, s_s, t) - epsilon_s * cos(2. * phi_s) * response_R_TT_ortho(reaction_s, mech_s, deg_k_s, deg_kstar_s, Q_sqr_s, s_s, t) + el_hel * (-sqrt(2. * epsilon_s * (1. - epsilon_s))) * sin(phi_s) * response_R_TLp_ortho(reaction_s, mech_s, deg_k_s, deg_kstar_s, Q_sqr_s, s_s, t); if(polbeam_s==1) num_perp = el_hel *( -sqrt(2. * epsilon_s * (1. - epsilon_s)) * sin(phi_s) * response_R_TLp_ortho(reaction_s, mech_s, deg_k_s, deg_kstar_s, Q_sqr_s, s_s, t)); } if(iintegr_s==1){ den = response_R_T(reaction_s, mech_s, deg_k_s, deg_kstar_s, Q_sqr_s, s_s, t) + epsilon_s * response_R_L(reaction_s,mech_s, deg_k_s, deg_kstar_s, Q_sqr_s, s_s, t); num_par = el_hel * sqrt(1.-pow(epsilon_s,2.)) * response_R_Tp_par(reaction_s, mech_s, deg_k_s, deg_kstar_s, Q_sqr_s, s_s, t); num_long = el_hel * sqrt(1.-pow(epsilon_s,2.)) * response_R_Tp_long(reaction_s, mech_s, deg_k_s, deg_kstar_s, Q_sqr_s, s_s, t); num_perp = 0.; } asymm_perp2 = num_perp / den; asymm_par2 = num_par / den; asymm_long2 = num_long / den; if (iintegr_s==2) { /* Faire une rotation en Phi car, dans cette methode, le repere est attache au plan hadronique : il faut repasser dans le plan leptonique */ asymm_par = asymm_par2*cos(phi_s)-asymm_perp2*sin(phi_s); asymm_perp = asymm_par2*sin(phi_s)+asymm_perp2*cos(phi_s); asymm_long = asymm_long2; } /* signe - en dessous pour la composante perp car dans le formalisme des fonctions de reponse, le Phi a ete inverse tout le long (phi->-phi)*/ asymm_perp = -asymm_perp; asymm_par = asymm_par; asymm_long = asymm_long; if (iintegr_s==2) { asymm_par_prime = asymm_par*cos(-phi_s)-asymm_perp*sin(-phi_s); asymm_perp_prime = asymm_par*sin(-phi_s)+asymm_perp*cos(-phi_s); asymm_long_prime = asymm_long; } if (iintegr_s==1) { asymm_par_prime = 0.; asymm_perp_prime = 0.; asymm_long_prime = asymm_long; } ga_mom_cm = sqrt( ( pow(s_s - pow(NU_M,2.) + Q_sqr_s, 2.) + pow(2. * NU_M, 2.) * Q_sqr_s ) / (4. * s_s) ); ga_en_cm = sqrt(s_s) - sqrt(pow(ga_mom_cm,2.) + pow(NU_M,2.)); k_mom_cm = sqrt( ( pow(s_s - pow(LAM_M,2.) - pow(K_M,2.), 2.) - pow(2. * LAM_M * K_M, 2.) ) / (4. * s_s) ); k_en_cm = sqrt( pow(k_mom_cm, 2.) + pow(K_M, 2.) ); cos_k_th = (t + Q_sqr_s - pow(K_M,2.) + 2. * ga_en_cm * k_en_cm) / (2. * ga_mom_cm * k_mom_cm) ; k_th=acos(cos_k_th); asymm_par_prime2 = asymm_par_prime*cos(k_th)-asymm_long_prime*sin(k_th); asymm_perp_prime2 = asymm_perp_prime; asymm_long_prime2 = asymm_par_prime*sin(k_th)+asymm_long_prime*cos(k_th); if(rep_s==1) { if (xyz_s==1) return asymm_par; if (xyz_s==2) return asymm_perp; if (xyz_s==3) return asymm_long; } if(rep_s==2) { if (xyz_s==1) return asymm_par_prime2; if (xyz_s==2) return asymm_perp_prime2; if (xyz_s==3) return asymm_long_prime2; } } /* The function v_compton_diff_cross_lab_test() returns the unpolarized differential cross section in the Lab system for the DVCS calculated from the response functions. The cross section dsig/dQ2dnudt is expressed in nb / ( GeV^5 ). */ double v_klam_5fold_diff_cross_test(int reaction, int mech, int deg_k, int deg_kstar, double el_hel, double k_in_mom_lab, double k_out_mom_lab, double k_out_th_lab, double kaon_th, double kaon_phi) { double ga_en_in_lab, ga_mom_in_lab, ga_en_cm, ga_mom_cm, k_mom_cm, k_en, epsilon, s, Q_sqr, t, phi; ga_en_in_lab = k_in_mom_lab - k_out_mom_lab; ga_mom_in_lab = sqrt(pow(k_in_mom_lab, 2.) + pow(k_out_mom_lab, 2.) - 2.*k_in_mom_lab*k_out_mom_lab * cos(k_out_th_lab) ); Q_sqr = 4. * k_in_mom_lab*k_out_mom_lab * pow(sin(0.5 * k_out_th_lab),2.); s = pow(NU_M, 2.) - Q_sqr + 2. * NU_M * ga_en_in_lab; epsilon = 1./(1. + 2. * 1./Q_sqr * pow(ga_mom_in_lab * tan(k_out_th_lab/2.), 2.) ); ga_mom_cm = sqrt( ( pow(s - pow(NU_M,2.) + Q_sqr, 2.) + pow(2. * NU_M, 2.) * Q_sqr ) / (4. * s) ); ga_en_cm = sqrt(s) - sqrt(pow(ga_mom_cm,2.) + pow(NU_M,2.)); k_mom_cm = sqrt( ( pow(s - pow(LAM_M,2.) - pow(K_M,2.), 2.) - pow(2. * LAM_M * K_M, 2.) ) / (4. * s) ); k_en = sqrt( pow(k_mom_cm, 2.) + pow(K_M, 2.) ); t = -Q_sqr + pow(K_M,2.) - 2. * ga_en_cm * k_en + 2. * ga_mom_cm * k_mom_cm * cos(kaon_th); /* phi=PI-kaon_phi;*/ phi=kaon_phi; /* printf("\n %lf %lf %lf %lf %lf", response_R_T(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t), response_R_L(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t), response_R_TT(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t), response_R_TL(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t), response_R_TLp(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t)); */ return pow(ELEC,2.) / (4. * PI) / (2. * pow( PI, 2.)) /( 2. * pow(k_in_mom_lab,2.)) * ga_mom_in_lab /Q_sqr /(1. - epsilon) * ( response_R_T(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t) + epsilon * response_R_L(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t) - epsilon * cos(2. * phi) * response_R_TT(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t) + sqrt(2. * epsilon * (1. + epsilon)) * cos(phi) * response_R_TL(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t) + el_hel * sqrt(2. * epsilon * (1. - epsilon)) * sin(phi) * response_R_TLp(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t) ) *(s-pow(NU_M,2.))/2./NU_M/ga_mom_in_lab; } /*************************************************************************/ /*************************************************************************/ /* */ /* e + N -> e + K + HYP : RESPONSE FUNCTIONS */ /* */ /*************************************************************************/ /*************************************************************************/ /* The function response_R_T() gives the response function R_T for the virtual photon + N -> pi + N process for kaon (1 = pi^+ n, 2 = pi^- p, 3 = pi^0 p) and the kinematic variables Q_sqr, s, and t. */ double response_R_T(int reaction, int mech, int deg_k, int deg_kstar, double Q_sqr, double s, double t) { int choice; double sp_out, sp_in, result_plus = 0., result_minus = 0., unit_conv, cross, Lambda; dcomplex contr_plus, contr_minus; for(sp_out = -0.5; sp_out <= 0.5; sp_out += 1.0) { for(sp_in = -0.5; sp_in <= 0.5; sp_in += 1.0) { contr_plus = J_hadron(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t, sp_out, +1, sp_in); contr_minus = J_hadron(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t, sp_out, -1, sp_in); result_plus += pow(Cabs(contr_plus),2.); result_minus += pow(Cabs(contr_minus),2.); } } choice = 2; /* default = 2 */ switch(choice) { case 1 : /* for the T cross section in units (microbarns / GeV^2) */ { unit_conv = 1. / (16. * PI * pow(s - pow(NU_M, 2.),2.)) * pow(0.197, 2.) * pow(10., 4.); return unit_conv * 0.25 * (result_plus + result_minus); break; } case 2 : /* for the T cross section in units (microbarns / GeV^2) (remark choice of virtual photon flux factor to compare to experimentally extracted reduced cross sections) */ { Lambda = sqrt( pow(s,2.) + pow(Q_sqr,2.) + pow(NU_M,4.) + 2. * s * Q_sqr - 2. * s * pow(NU_M,2.) + 2. * pow(NU_M,2.) * Q_sqr ); unit_conv = 1. / (16. * PI) * 1./(s - pow(NU_M, 2.)) * 1./Lambda * pow(0.197, 2.) * pow(10., 4.); return unit_conv * 0.25 * (result_plus + result_minus); break; } case 3 : /* for the R_T response function */ { return 0.25 * (result_plus + result_minus); break; } } return 0.; } /* The function response_R_L() gives the response function R_L for the virtual photon + N -> pi + N process for kaon (1 = pi^+ n, 2 = pi^- p, 3 = pi^0 p) and the kinematic variables Q_sqr, s, and t. */ double response_R_L(int reaction, int mech, int deg_k, int deg_kstar, double Q_sqr, double s, double t) { int choice; double sp_out, sp_in, result_zero = 0., unit_conv, cross, Lambda; dcomplex contr_zero; for(sp_out = -0.5; sp_out <= 0.5; sp_out += 1.0) { for(sp_in = -0.5; sp_in <= 0.5; sp_in += 1.0) { contr_zero = J_hadron(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t, sp_out, 0, sp_in); result_zero += pow(Cabs(contr_zero),2.); } } choice = 2; /* default = 2 */ switch(choice) { case 1 : /* for the L cross section in units (microbarns / GeV^2) */ { unit_conv = 1. / (16. * PI * pow(s - pow(NU_M, 2.),2.)) * pow(0.197, 2.) * pow(10., 4.); return unit_conv * 0.5 * result_zero; break; } case 2 : /* for the L cross section in units (microbarns / GeV^2) (remark choice of virtual photon flux factor to compare to experimentally extracted reduced cross sections) */ { Lambda = sqrt( pow(s,2.) + pow(Q_sqr,2.) + pow(NU_M,4.) + 2. * s * Q_sqr - 2. * s * pow(NU_M,2.) + 2. * pow(NU_M,2.) * Q_sqr ); unit_conv = 1. / (16. * PI) * 1./(s - pow(NU_M, 2.)) * 1./Lambda * pow(0.197, 2.) * pow(10., 4.); return unit_conv * 0.5 * result_zero; break; } case 3 : /* for the R_L response function */ { return 0.5 * result_zero; break; } } return 0.; } /* The function response_R_TT() gives the response function R_TT for the virtual photon + N -> pi + N process for kaon (1 = pi^+ n, 2 = pi^- p, 3 = pi^0 p) and the kinematic variables Q_sqr, s, and t. */ double response_R_TT(int reaction, int mech, int deg_k, int deg_kstar, double Q_sqr, double s, double t) { int choice; double sp_out, sp_in, result = 0., unit_conv, cross, Lambda; dcomplex contr; for(sp_in = -0.5; sp_in <= 0.5; sp_in += 1.0) { contr = Cmul(J_hadron(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t, .5, +1, sp_in), Conjg(J_hadron(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t, .5, -1, sp_in) ) ); result += 2.*contr.re; } choice = 2; /* default = 2 */ switch(choice) { case 1 : /* for the TT cross section in units (microbarns / GeV^2). */ { unit_conv = 1. / (16. * PI * pow(s - pow(NU_M, 2.),2.)) * pow(0.197, 2.) * pow(10., 4.); return unit_conv * 0.5 * result; break; } case 2 : /* for the TT cross section in units (microbarns / GeV^2) (remark choice of virtual photon flux factor to compare to experimentally extracted reduced cross sections) */ { Lambda = sqrt( pow(s,2.) + pow(Q_sqr,2.) + pow(NU_M,4.) + 2. * s * Q_sqr - 2. * s * pow(NU_M,2.) + 2. * pow(NU_M,2.) * Q_sqr ); unit_conv = 1. / (16. * PI) * 1./(s - pow(NU_M, 2.)) * 1./Lambda * pow(0.197, 2.) * pow(10., 4.); return unit_conv * 0.5 * result; break; } case 3 : /* for the R_TT response function */ { return 0.5 * result; break; } } return 0.; } /* The function response_R_TL() gives the response function R_TL for the virtual photon + N -> pi + N process for kaon (1 = pi^+ n, 2 = pi^- p, 3 = pi^0 p) and the kinematic variables Q_sqr, s, and t. */ double response_R_TL(int reaction, int mech, int deg_k, int deg_kstar, double Q_sqr, double s, double t) { int choice; double sp_out, sp_in, result = 0., unit_conv, cross, Lambda; dcomplex contr_plus, contr_minus, contr_zero; /* for(sp_out = -0.5; sp_out <= 0.5; sp_out += 1.0) { contr_plus = J_hadron(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t, sp_out, +1, .5); contr_minus = J_hadron(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t, sp_out, -1, .5); contr_zero = J_hadron(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t, sp_out, 0, .5); result += 2.*Cmul(Csub(contr_plus, contr_minus), Conjg(contr_zero) ).re; } */ for(sp_in = -0.5; sp_in <= 0.5; sp_in += 1.0) { contr_plus = J_hadron(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t, .5, +1, sp_in); contr_minus = J_hadron(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t, .5, -1, sp_in); contr_zero = J_hadron(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t, .5, 0, sp_in); contr_zero=Cmul(Complex(-1.,0.),contr_zero); result += 2.*Cmul(Csub(contr_plus, contr_minus), Conjg(contr_zero) ).re; } /* for(sp_out = -0.5; sp_out <= 0.5; sp_out += 1.0) { for(sp_in = -0.5; sp_in <= 0.5; sp_in += 1.0) { contr_plus = J_hadron(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t, sp_out, +1, sp_in); contr_minus = J_hadron(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t, sp_out, -1, sp_in); contr_zero = J_hadron(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t, sp_out, 0, sp_in); result += Cmul(Csub(contr_plus, contr_minus), Conjg(contr_zero) ).re; } } */ choice = 2; /* default = 2 */ switch(choice) { case 1 : /* for the TL cross section in units (microbarns / GeV^2). */ { unit_conv = 1. / (16. * PI * pow(s - pow(NU_M, 2.),2.)) * pow(0.197, 2.) * pow(10., 4.); return unit_conv * 0.5 / sqrt(2.) * result; break; } case 2 : /* for the TL cross section in units (microbarns / GeV^2) (remark choice of virtual photon flux factor to compare to experimentally extracted reduced cross sections) */ { Lambda = sqrt( pow(s,2.) + pow(Q_sqr,2.) + pow(NU_M,4.) + 2. * s * Q_sqr - 2. * s * pow(NU_M,2.) + 2. * pow(NU_M,2.) * Q_sqr ); unit_conv = 1. / (16. * PI) * 1./(s - pow(NU_M, 2.)) * 1./Lambda * pow(0.197, 2.) * pow(10., 4.); return unit_conv * 0.5 / sqrt(2.) * result ; break; } case 3 : /* for the R_TL response function */ { return 0.5 / sqrt(2.) * result; break; } } return 0.; } /* The function response_R_TLp() gives the response function R_TL' for the virtual photon + N -> pi + N process for kaon (1 = pi^+ n, 2 = pi^- p, 3 = pi^0 p) and the kinematic variables Q_sqr, s, and t. */ double response_R_TLp(int reaction, int mech, int deg_k, int deg_kstar, double Q_sqr, double s, double t) { int choice; double sp_out, sp_in, result = 0., unit_conv, cross, Lambda; dcomplex contr_plus, contr_minus, contr_zero; for(sp_in = -0.5; sp_in <= 0.5; sp_in += 1.0) { contr_plus = J_hadron(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t, .5, +1, sp_in); contr_minus = J_hadron(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t, .5, -1, sp_in); contr_zero = J_hadron(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t, .5, 0, sp_in); /* contr_zero=Cmul(Complex(-1.,0.),contr_zero);*/ result += 2.*Cmul(Csub(contr_plus, contr_minus), Conjg(contr_zero) ).im; } /* for(sp_out = -0.5; sp_out <= 0.5; sp_out += 1.0) { for(sp_in = -0.5; sp_in <= 0.5; sp_in += 1.0) { contr_plus = J_hadron(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t, sp_out, +1, sp_in); contr_minus = J_hadron(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t, sp_out, -1, sp_in); contr_zero = J_hadron(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t, sp_out, 0, sp_in); result += Cmul(Csub(contr_plus, contr_minus), Conjg(contr_zero) ).im; } } */ choice = 2; /* default = 2 */ switch(choice) { case 1 : /* for the TL' cross section in units (microbarns / GeV^2). */ { unit_conv = 1. / (16. * PI * pow(s - pow(NU_M, 2.),2.)) * pow(0.197, 2.) * pow(10., 4.); return unit_conv * 0.5 / sqrt(2.) * result; break; } case 2 : /* for the TL' cross section in units (microbarns / GeV^2) (remark choice of virtual photon flux factor to compare to experimentally extracted reduced cross sections) */ { Lambda = sqrt( pow(s,2.) + pow(Q_sqr,2.) + pow(NU_M,4.) + 2. * s * Q_sqr - 2. * s * pow(NU_M,2.) + 2. * pow(NU_M,2.) * Q_sqr ); unit_conv = 1. / (16. * PI) * 1./(s - pow(NU_M, 2.)) * 1./Lambda * pow(0.197, 2.) * pow(10., 4.); return unit_conv * 0.5 / sqrt(2.) * result ; break; } case 3 : /* for the R_TL' response function */ { return 0.5 / sqrt(2.) * result; break; } } return 0.; } /* The function response_R_T_ortho() gives the response function R_T_ortho for the virtual photon + N -> pi + N process for kaon (1 = pi^+ n, 2 = pi^- p, 3 = pi^0 p) and the kinematic variables Q_sqr, s, and t. */ double response_R_T_ortho(int reaction, int mech, int deg_k, int deg_kstar, double Q_sqr, double s, double t) { int choice; double sp_out, sp_in, result = 0., unit_conv, cross, Lambda; dcomplex contr; for(sp_in = -0.5; sp_in <= 0.5; sp_in += 1.0) { contr = Cmul(J_hadron(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t, .5, +1, sp_in), Conjg(J_hadron(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t, -.5, 1, sp_in) ) ); /* for(sp_out = -0.5; sp_out <= 0.5; sp_out += 1.0) { contr = Cmul(J_hadron(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t, sp_out, +1, .5), Conjg(J_hadron(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t, sp_out, 1, -.5) ) ); */ result += 2.*contr.im; } choice = 2; /* default = 2 */ switch(choice) { case 1 : /* for the T_ortho cross section in units (microbarns / GeV^2). */ { unit_conv = 1. / (16. * PI * pow(s - pow(NU_M, 2.),2.)) * pow(0.197, 2.) * pow(10., 4.); return unit_conv * 0.5 * result; break; } case 2 : /* for the T_ortho cross section in units (microbarns / GeV^2) (remark choice of virtual photon flux factor to compare to experimentally extracted reduced cross sections) */ { Lambda = sqrt( pow(s,2.) + pow(Q_sqr,2.) + pow(NU_M,4.) + 2. * s * Q_sqr - 2. * s * pow(NU_M,2.) + 2. * pow(NU_M,2.) * Q_sqr ); unit_conv = 1. / (16. * PI) * 1./(s - pow(NU_M, 2.)) * 1./Lambda * pow(0.197, 2.) * pow(10., 4.); return unit_conv * 0.5 * result; break; } case 3 : /* for the R_T_ortho response function */ { return 0.5 * result; break; } } return 0.; } /* The function response_R_L_ortho() gives the response function R_L_ortho for the virtual photon + N -> pi + N process for kaon (1 = pi^+ n, 2 = pi^- p, 3 = pi^0 p) and the kinematic variables Q_sqr, s, and t. */ double response_R_L_ortho(int reaction, int mech, int deg_k, int deg_kstar, double Q_sqr, double s, double t) { int choice; double sp_out, sp_in, result = 0., unit_conv, cross, Lambda; dcomplex contr; for(sp_in = -0.5; sp_in <= 0.5; sp_in += 1.0) { contr = Cmul(J_hadron(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t, .5, 0, sp_in), Conjg(J_hadron(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t, -.5, 0, sp_in) ) ); /* for(sp_out = -0.5; sp_out <= 0.5; sp_out += 1.0) { contr = Cmul(J_hadron(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t, sp_out, 0, .5), Conjg(J_hadron(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t, sp_out, 0, -.5) ) ); */ result += 2.*contr.im; } choice = 2; /* default = 2 */ switch(choice) { case 1 : /* for the L_ortho cross section in units (microbarns / GeV^2). */ { unit_conv = 1. / (16. * PI * pow(s - pow(NU_M, 2.),2.)) * pow(0.197, 2.) * pow(10., 4.); return unit_conv * 0.5 * result; break; } case 2 : /* for the L_ortho cross section in units (microbarns / GeV^2) (remark choice of virtual photon flux factor to compare to experimentally extracted reduced cross sections) */ { Lambda = sqrt( pow(s,2.) + pow(Q_sqr,2.) + pow(NU_M,4.) + 2. * s * Q_sqr - 2. * s * pow(NU_M,2.) + 2. * pow(NU_M,2.) * Q_sqr ); unit_conv = 1. / (16. * PI) * 1./(s - pow(NU_M, 2.)) * 1./Lambda * pow(0.197, 2.) * pow(10., 4.); return unit_conv * 0.5 * result; break; } case 3 : /* for the R_L_ortho response function */ { return 0.5 * result; break; } } return 0.; } /* The function response_R_TT_ortho() gives the response function R_TT_ortho for the virtual photon + N -> pi + N process for kaon (1 = pi^+ n, 2 = pi^- p, 3 = pi^0 p) and the kinematic variables Q_sqr, s, and t. */ double response_R_TT_ortho(int reaction, int mech, int deg_k, int deg_kstar, double Q_sqr, double s, double t) { int choice; double sp_out, sp_in, result = 0., unit_conv, cross, Lambda; dcomplex contr_a, contr_b; for(sp_in = -0.5; sp_in <= 0.5; sp_in += 1.0) { contr_a = Cmul(J_hadron(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t, .5, +1, sp_in), Conjg(J_hadron(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t, -.5, -1, sp_in) ) ); contr_b = Cmul(J_hadron(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t, .5, -1, sp_in), Conjg(J_hadron(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t, -.5, 1, sp_in) ) ); /* for(sp_out = -0.5; sp_out <= 0.5; sp_out += 1.0) { contr_a = Cmul(J_hadron(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t, sp_out, +1, .5), Conjg(J_hadron(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t, sp_out, -1, -.5) ) ); contr_b = Cmul(J_hadron(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t, sp_out, -1, .5), Conjg(J_hadron(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t, sp_out, 1, -.5) ) ); */ result += contr_a.im + contr_b.im; } choice = 2; /* default = 2 */ switch(choice) { case 1 : /* for the TT_ortho cross section in units (microbarns / GeV^2). */ { unit_conv = 1. / (16. * PI * pow(s - pow(NU_M, 2.),2.)) * pow(0.197, 2.) * pow(10., 4.); return unit_conv * 0.5 * result; break; } case 2 : /* for the TT_ortho cross section in units (microbarns / GeV^2) (remark choice of virtual photon flux factor to compare to experimentally extracted reduced cross sections) */ { Lambda = sqrt( pow(s,2.) + pow(Q_sqr,2.) + pow(NU_M,4.) + 2. * s * Q_sqr - 2. * s * pow(NU_M,2.) + 2. * pow(NU_M,2.) * Q_sqr ); unit_conv = 1. / (16. * PI) * 1./(s - pow(NU_M, 2.)) * 1./Lambda * pow(0.197, 2.) * pow(10., 4.); return unit_conv * 0.5 * result; break; } case 3 : /* for the R_TT_ortho response function */ { return 0.5 * result; break; } } return 0.; } /* The function response_R_TL_ortho() gives the response function R_TL_ortho for the virtual photon + N -> pi + N process for kaon (1 = pi^+ n, 2 = pi^- p, 3 = pi^0 p) and the kinematic variables Q_sqr, s, and t. */ double response_R_TL_ortho(int reaction, int mech, int deg_k, int deg_kstar, double Q_sqr, double s, double t) { int choice; double sp_out, sp_in, result = 0., unit_conv, cross, Lambda; dcomplex contr_plus, contr_minus, contr_zero; for(sp_in = -0.5; sp_in <= 0.5; sp_in += 1.0) { contr_plus = J_hadron(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t, .5, +1, sp_in); contr_minus = J_hadron(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t, .5, -1, sp_in); contr_zero = J_hadron(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t, -.5, 0, sp_in); /* for(sp_out = -0.5; sp_out <= 0.5; sp_out += 1.0) { contr_plus = J_hadron(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t, sp_out, +1, .5); contr_minus = J_hadron(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t, sp_out, -1, .5); contr_zero = J_hadron(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t, sp_out, 0, -.5); */ result += 2. * Cmul(Csub(contr_plus, contr_minus), Conjg(contr_zero) ).im; } choice = 2; /* default = 2 */ switch(choice) { case 1 : /* for the TL_ortho cross section in units (microbarns / GeV^2). */ { unit_conv = 1. / (16. * PI * pow(s - pow(NU_M, 2.),2.)) * pow(0.197, 2.) * pow(10., 4.); return unit_conv * 0.5 / sqrt(2.) * result; break; } case 2 : /* for the TL_ortho cross section in units (microbarns / GeV^2) (remark choice of virtual photon flux factor to compare to experimentally extracted reduced cross sections) */ { Lambda = sqrt( pow(s,2.) + pow(Q_sqr,2.) + pow(NU_M,4.) + 2. * s * Q_sqr - 2. * s * pow(NU_M,2.) + 2. * pow(NU_M,2.) * Q_sqr ); unit_conv = 1. / (16. * PI) * 1./(s - pow(NU_M, 2.)) * 1./Lambda * pow(0.197, 2.) * pow(10., 4.); return unit_conv * 0.5 / sqrt(2.) * result ; break; } case 3 : /* for the R_TL_ortho response function */ { return 0.5 / sqrt(2.) * result; break; } } return 0.; } /* The function response_R_TLp_ortho() gives the response function R_TL'_ortho for the virtual photon + N -> pi + N process for kaon (1 = pi^+ n, 2 = pi^- p, 3 = pi^0 p) and the kinematic variables Q_sqr, s, and t. */ double response_R_TLp_ortho(int reaction, int mech, int deg_k, int deg_kstar, double Q_sqr, double s, double t) { int choice; double sp_out, sp_in, result = 0., unit_conv, cross, Lambda; dcomplex contr_plus, contr_minus, contr_zero; for(sp_in = -0.5; sp_in <= 0.5; sp_in += 1.0) { contr_plus = J_hadron(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t, .5, +1, sp_in); contr_minus = J_hadron(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t, .5, -1, sp_in); contr_zero = J_hadron(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t, -.5, 0, sp_in); /* contr_zero=Cmul(Complex(-1.,0.),contr_zero);*/ result += -2. * Cmul(Csub(contr_plus, contr_minus), Conjg(contr_zero) ).re; } choice = 2; /* default = 2 */ switch(choice) { case 1 : /* for the TL' cross section in units (microbarns / GeV^2). */ { unit_conv = 1. / (16. * PI * pow(s - pow(NU_M, 2.),2.)) * pow(0.197, 2.) * pow(10., 4.); return unit_conv * 0.5 / sqrt(2.) * result; break; } case 2 : /* for the TL' cross section in units (microbarns / GeV^2) (remark choice of virtual photon flux factor to compare to experimentally extracted reduced cross sections) */ { Lambda = sqrt( pow(s,2.) + pow(Q_sqr,2.) + pow(NU_M,4.) + 2. * s * Q_sqr - 2. * s * pow(NU_M,2.) + 2. * pow(NU_M,2.) * Q_sqr ); unit_conv = 1. / (16. * PI) * 1./(s - pow(NU_M, 2.)) * 1./Lambda * pow(0.197, 2.) * pow(10., 4.); return unit_conv * 0.5 / sqrt(2.) * result ; break; } case 3 : /* for the R_TL' response function */ { return 0.5 / sqrt(2.) * result; break; } } return 0.; } /* The function response_R_TL_par() gives the response function R_TL_par for the virtual photon + N -> pi + N process for kaon (1 = pi^+ n, 2 = pi^- p, 3 = pi^0 p) and the kinematic variables Q_sqr, s, and t. */ double response_R_TL_par(int reaction, int mech, int deg_k, int deg_kstar, double Q_sqr, double s, double t) { int choice; double sp_out, sp_in, result = 0., unit_conv, cross, Lambda; dcomplex contr_plus, contr_minus, contr_zero; for(sp_in = -0.5; sp_in <= 0.5; sp_in += 1.0) { contr_plus = J_hadron(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t, .5, +1, sp_in); contr_minus = J_hadron(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t, .5, -1, sp_in); contr_zero = J_hadron(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t, -.5, 0, sp_in); /* contr_zero=Cmul(Complex(-1.,0.),contr_zero);*/ result += 2.*Cmul(Cadd(contr_plus, contr_minus), Conjg(contr_zero) ).im; } choice = 2; /* default = 2 */ switch(choice) { case 1 : /* for the TL_par cross section in units (microbarns / GeV^2). */ { unit_conv = 1. / (16. * PI * pow(s - pow(NU_M, 2.),2.)) * pow(0.197, 2.) * pow(10., 4.); return unit_conv * 0.5 / sqrt(2.) * result; break; } case 2 : /* for the TL_par cross section in units (microbarns / GeV^2) (remark choice of virtual photon flux factor to compare to experimentally extracted reduced cross sections) */ { Lambda = sqrt( pow(s,2.) + pow(Q_sqr,2.) + pow(NU_M,4.) + 2. * s * Q_sqr - 2. * s * pow(NU_M,2.) + 2. * pow(NU_M,2.) * Q_sqr ); unit_conv = 1. / (16. * PI) * 1./(s - pow(NU_M, 2.)) * 1./Lambda * pow(0.197, 2.) * pow(10., 4.); return unit_conv * 0.5 / sqrt(2.) * result ; break; } case 3 : /* for the R_TL_par response function */ { return 0.5 / sqrt(2.) * result; break; } } return 0.; } /* The function response_R_TT_par() gives the response function R_TT_par for the virtual photon + N -> pi + N process for kaon (1 = pi^+ n, 2 = pi^- p, 3 = pi^0 p) and the kinematic variables Q_sqr, s, and t. */ double response_R_TT_par(int reaction, int mech, int deg_k, int deg_kstar, double Q_sqr, double s, double t) { int choice; double sp_out, sp_in, result = 0., unit_conv, cross, Lambda; dcomplex contr_a, contr_b; for(sp_in = -0.5; sp_in <= 0.5; sp_in += 1.0) { contr_a = Cmul(J_hadron(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t, .5, +1, sp_in), Conjg(J_hadron(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t, -.5, -1, sp_in) ) ); contr_b = Cmul(J_hadron(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t, .5, -1, sp_in), Conjg(J_hadron(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t, -.5, 1, sp_in) ) ); result += contr_a.im - contr_b.im; } choice = 2; /* default = 2 */ switch(choice) { case 1 : /* for the TT_par cross section in units (microbarns / GeV^2). */ { unit_conv = 1. / (16. * PI * pow(s - pow(NU_M, 2.),2.)) * pow(0.197, 2.) * pow(10., 4.); return unit_conv * 0.5 * result; break; } case 2 : /* for the TT_par cross section in units (microbarns / GeV^2) (remark choice of virtual photon flux factor to compare to experimentally extracted reduced cross sections) */ { Lambda = sqrt( pow(s,2.) + pow(Q_sqr,2.) + pow(NU_M,4.) + 2. * s * Q_sqr - 2. * s * pow(NU_M,2.) + 2. * pow(NU_M,2.) * Q_sqr ); unit_conv = 1. / (16. * PI) * 1./(s - pow(NU_M, 2.)) * 1./Lambda * pow(0.197, 2.) * pow(10., 4.); return unit_conv * 0.5 * result; break; } case 3 : /* for the R_TT_par response function */ { return 0.5 * result; break; } } return 0.; } /* The function response_R_TLp_par() gives the response function R_TLp_par for the virtual photon + N -> pi + N process for kaon (1 = pi^+ n, 2 = pi^- p, 3 = pi^0 p) and the kinematic variables Q_sqr, s, and t. */ double response_R_TLp_par(int reaction, int mech, int deg_k, int deg_kstar, double Q_sqr, double s, double t) { int choice; double sp_out, sp_in, result = 0., unit_conv, cross, Lambda; dcomplex contr_plus, contr_minus, contr_zero; for(sp_in = -0.5; sp_in <= 0.5; sp_in += 1.0) { contr_plus = J_hadron(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t, .5, +1, sp_in); contr_minus = J_hadron(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t, .5, -1, sp_in); contr_zero = J_hadron(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t, -.5, 0, sp_in); /* contr_zero=Cmul(Complex(-1.,0.),contr_zero);*/ result += 2.*Cmul(Cadd(contr_plus, contr_minus), Conjg(contr_zero) ).re; } choice = 2; /* default = 2 */ switch(choice) { case 1 : /* for the TLp_par cross section in units (microbarns / GeV^2). */ { unit_conv = 1. / (16. * PI * pow(s - pow(NU_M, 2.),2.)) * pow(0.197, 2.) * pow(10., 4.); return unit_conv * 0.5 / sqrt(2.) * result; break; } case 2 : /* for the TLp_par cross section in units (microbarns / GeV^2) (remark choice of virtual photon flux factor to compare to experimentally extracted reduced cross sections) */ { Lambda = sqrt( pow(s,2.) + pow(Q_sqr,2.) + pow(NU_M,4.) + 2. * s * Q_sqr - 2. * s * pow(NU_M,2.) + 2. * pow(NU_M,2.) * Q_sqr ); unit_conv = 1. / (16. * PI) * 1./(s - pow(NU_M, 2.)) * 1./Lambda * pow(0.197, 2.) * pow(10., 4.); return unit_conv * 0.5 / sqrt(2.) * result ; break; } case 3 : /* for the R_TLp_par response function */ { return 0.5 / sqrt(2.) * result; break; } } return 0.; } /* The function response_R_Tp_par() gives the response function R_Tp_par for the virtual photon + N -> pi + N process for kaon (1 = pi^+ n, 2 = pi^- p, 3 = pi^0 p) and the kinematic variables Q_sqr, s, and t. */ double response_R_Tp_par(int reaction, int mech, int deg_k, int deg_kstar, double Q_sqr, double s, double t) { int choice; double sp_out, sp_in, result = 0., unit_conv, cross, Lambda; dcomplex contr; for(sp_in = -0.5; sp_in <= 0.5; sp_in += 1.0) { contr = Cmul(J_hadron(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t, .5, +1, sp_in), Conjg(J_hadron(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t, -.5, 1, sp_in) ) ); result += 2.*contr.re; } choice = 2; /* default = 2 */ switch(choice) { case 1 : /* for the Tp_par cross section in units (microbarns / GeV^2). */ { unit_conv = 1. / (16. * PI * pow(s - pow(NU_M, 2.),2.)) * pow(0.197, 2.) * pow(10., 4.); return unit_conv * 0.5 * result; break; } case 2 : /* for the Tp_par cross section in units (microbarns / GeV^2) (remark choice of virtual photon flux factor to compare to experimentally extracted reduced cross sections) */ { Lambda = sqrt( pow(s,2.) + pow(Q_sqr,2.) + pow(NU_M,4.) + 2. * s * Q_sqr - 2. * s * pow(NU_M,2.) + 2. * pow(NU_M,2.) * Q_sqr ); unit_conv = 1. / (16. * PI) * 1./(s - pow(NU_M, 2.)) * 1./Lambda * pow(0.197, 2.) * pow(10., 4.); return unit_conv * 0.5 * result; break; } case 3 : /* for the R_Tp_par response function */ { return 0.5 * result; break; } } return 0.; } /* The function response_R_TL_long() gives the response function R_TL_long for the virtual photon + N -> pi + N process for kaon (1 = pi^+ n, 2 = pi^- p, 3 = pi^0 p) and the kinematic variables Q_sqr, s, and t. */ double response_R_TL_long(int reaction, int mech, int deg_k, int deg_kstar, double Q_sqr, double s, double t) { int choice; double sp_out, sp_in, result = 0., unit_conv, cross, Lambda; dcomplex contr_plus, contr_minus, contr_zero; for(sp_in = -0.5; sp_in <= 0.5; sp_in += 1.0) { contr_plus = J_hadron(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t, .5, +1, sp_in); contr_minus = J_hadron(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t, .5, -1, sp_in); contr_zero = J_hadron(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t, .5, 0, sp_in); /* contr_zero=Cmul(Complex(-1.,0.),contr_zero);*/ result += 2.*Cmul(Cadd(contr_plus, contr_minus), Conjg(contr_zero) ).im; } choice = 2; /* default = 2 */ switch(choice) { case 1 : /* for the TL_long cross section in units (microbarns / GeV^2). */ { unit_conv = 1. / (16. * PI * pow(s - pow(NU_M, 2.),2.)) * pow(0.197, 2.) * pow(10., 4.); return unit_conv * 0.5 / sqrt(2.) * result; break; } case 2 : /* for the TL_long cross section in units (microbarns / GeV^2) (remark choice of virtual photon flux factor to compare to experimentally extracted reduced cross sections) */ { Lambda = sqrt( pow(s,2.) + pow(Q_sqr,2.) + pow(NU_M,4.) + 2. * s * Q_sqr - 2. * s * pow(NU_M,2.) + 2. * pow(NU_M,2.) * Q_sqr ); unit_conv = 1. / (16. * PI) * 1./(s - pow(NU_M, 2.)) * 1./Lambda * pow(0.197, 2.) * pow(10., 4.); return unit_conv * 0.5 / sqrt(2.) * result ; break; } case 3 : /* for the R_TL_long response function */ { return 0.5 / sqrt(2.) * result; break; } } return 0.; } /* The function response_R_TT_long() gives the response function R_TT_long for the virtual photon + N -> pi + N process for kaon (1 = pi^+ n, 2 = pi^- p, 3 = pi^0 p) and the kinematic variables Q_sqr, s, and t. */ double response_R_TT_long(int reaction, int mech, int deg_k, int deg_kstar, double Q_sqr, double s, double t) { int choice; double sp_out, sp_in, result = 0., unit_conv, cross, Lambda; dcomplex contr; for(sp_in = -0.5; sp_in <= 0.5; sp_in += 1.0) { contr = Cmul(J_hadron(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t, .5, +1, sp_in), Conjg(J_hadron(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t, .5, -1, sp_in) ) ); result += 2. * contr.im; } choice = 2; /* default = 2 */ switch(choice) { case 1 : /* for the TT_long cross section in units (microbarns / GeV^2). */ { unit_conv = 1. / (16. * PI * pow(s - pow(NU_M, 2.),2.)) * pow(0.197, 2.) * pow(10., 4.); return unit_conv * 0.5 * result; break; } case 2 : /* for the TT_long cross section in units (microbarns / GeV^2) (remark choice of virtual photon flux factor to compare to experimentally extracted reduced cross sections) */ { Lambda = sqrt( pow(s,2.) + pow(Q_sqr,2.) + pow(NU_M,4.) + 2. * s * Q_sqr - 2. * s * pow(NU_M,2.) + 2. * pow(NU_M,2.) * Q_sqr ); unit_conv = 1. / (16. * PI) * 1./(s - pow(NU_M, 2.)) * 1./Lambda * pow(0.197, 2.) * pow(10., 4.); return unit_conv * 0.5 * result; break; } case 3 : /* for the R_TT_long response function */ { return 0.5 * result; break; } } return 0.; } /* The function response_R_TLp_long() gives the response function R_TLp_long for the virtual photon + N -> pi + N process for kaon (1 = pi^+ n, 2 = pi^- p, 3 = pi^0 p) and the kinematic variables Q_sqr, s, and t. */ double response_R_TLp_long(int reaction, int mech, int deg_k, int deg_kstar, double Q_sqr, double s, double t) { int choice; double sp_out, sp_in, result = 0., unit_conv, cross, Lambda; dcomplex contr_plus, contr_minus, contr_zero; for(sp_in = -0.5; sp_in <= 0.5; sp_in += 1.0) { contr_plus = J_hadron(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t, .5, +1, sp_in); contr_minus = J_hadron(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t, .5, -1, sp_in); contr_zero = J_hadron(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t, .5, 0, sp_in); /* contr_zero=Cmul(Complex(-1.,0.),contr_zero);*/ result += 2.*Cmul(Cadd(contr_plus, contr_minus), Conjg(contr_zero) ).re; } choice = 2; /* default = 2 */ switch(choice) { case 1 : /* for the TLp_long cross section in units (microbarns / GeV^2). */ { unit_conv = 1. / (16. * PI * pow(s - pow(NU_M, 2.),2.)) * pow(0.197, 2.) * pow(10., 4.); return unit_conv * 0.5 / sqrt(2.) * result; break; } case 2 : /* for the TLp_long cross section in units (microbarns / GeV^2) (remark choice of virtual photon flux factor to compare to experimentally extracted reduced cross sections) */ { Lambda = sqrt( pow(s,2.) + pow(Q_sqr,2.) + pow(NU_M,4.) + 2. * s * Q_sqr - 2. * s * pow(NU_M,2.) + 2. * pow(NU_M,2.) * Q_sqr ); unit_conv = 1. / (16. * PI) * 1./(s - pow(NU_M, 2.)) * 1./Lambda * pow(0.197, 2.) * pow(10., 4.); return unit_conv * 0.5 / sqrt(2.) * result; break; } case 3 : /* for the R_TLp_long response function */ { return 0.5 / sqrt(2.) * result ; break; } } return 0.; } /* The function response_R_Tp_long() gives the response function R_Tp_long for the virtual photon + N -> pi + N process for kaon (1 = pi^+ n, 2 = pi^- p, 3 = pi^0 p) and the kinematic variables Q_sqr, s, and t. */ double response_R_Tp_long(int reaction, int mech, int deg_k, int deg_kstar, double Q_sqr, double s, double t) { int choice; double sp_out, sp_in, result_plus= 0., result_minus= 0., unit_conv, cross, Lambda; dcomplex contr_plus, contr_minus; for(sp_in = -0.5; sp_in <= 0.5; sp_in += 1.0) { contr_plus = J_hadron(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t, .5, +1, sp_in); contr_minus = J_hadron(reaction, mech, deg_k, deg_kstar, Q_sqr, s, t, .5, -1, sp_in); result_plus += pow(Cabs(contr_plus),2.); result_minus += pow(Cabs(contr_minus),2.); } choice = 2; /* default = 2 */ switch(choice) { case 1 : /* for the Tp_par cross section in units (microbarns / GeV^2). */ { unit_conv = 1. / (16. * PI * pow(s - pow(NU_M, 2.),2.)) * pow(0.197, 2.) * pow(10., 4.); return unit_conv * 0.5 * (result_plus - result_minus); break; } case 2 : /* for the Tp_par cross section in units (microbarns / GeV^2) (remark choice of virtual photon flux factor to compare to experimentally extracted reduced cross sections) */ { Lambda = sqrt( pow(s,2.) + pow(Q_sqr,2.) + pow(NU_M,4.) + 2. * s * Q_sqr - 2. * s * pow(NU_M,2.) + 2. * pow(NU_M,2.) * Q_sqr ); unit_conv = 1. / (16. * PI) * 1./(s - pow(NU_M, 2.)) * 1./Lambda * pow(0.197, 2.) * pow(10., 4.); return unit_conv * 0.5 * (result_plus - result_minus); break; } case 3 : /* for the R_Tp_par response function */ { return 0.5 * (result_plus - result_minus); break; } } return 0.; } /* The function recoil_asymm() returns the recoil nucleon asymmetry for the kaon photoproduction process on a nucleon : ga + N to pi + N. */ double recoil_asymm(int reaction, int mech, int deg_k, int deg_kstar, double Q_sqr, double s, double t) { int lam; double sp_in; double asymm_up = 0., asymm_down = 0.; dcomplex t_up, ampl_up, ampl1_up, ampl2_up, t_down, ampl_down, ampl1_down, ampl2_down; for( lam = -1; lam <= 1; lam += 2) { for( sp_in = -.5; sp_in <= .5; sp_in += 1.) { ampl_up= t_vgap_klam_regge_phi(mech, deg_k, deg_kstar, lam, .5, sp_in, Q_sqr, s, t, 0.); ampl_down= t_vgap_klam_regge_phi(mech, deg_k, deg_kstar, lam, -.5, sp_in, Q_sqr, s, t, 0.); t_up = Cadd(Cmul(Complex(1./ sqrt(2.), 0.), ampl_up), Cmul(Complex(0., - 1./ sqrt(2.)), ampl_down) ); t_down = Cadd(Cmul(Complex(- 1./ sqrt(2.), 0.), ampl_up), Cmul(Complex(0., - 1./ sqrt(2.)), ampl_down) ); asymm_up += pow(Cabs(t_up), 2.); asymm_down += pow(Cabs(t_down), 2.); } } return ( asymm_up - asymm_down ) / (asymm_up + asymm_down); } /*************************************************************************/ /*************************************************************************/ /* The relativistic hadronic current operator considered is given by J_hadron(). */ dcomplex J_hadron(int reaction, int mech, int deg_k, int deg_kstar, double Q_sqr, double s, double t, double sp_out, int pol_ph, double sp_in) { switch(reaction) { case 1: { return t_vgap_klam_regge_forward(mech, deg_k, deg_kstar, pol_ph, sp_out, sp_in, Q_sqr, s, t); break; } case 2: { return t_vgap_ksigma_regge_forward(mech, deg_k, deg_kstar, pol_ph, sp_out, sp_in, Q_sqr, s, t); break; } case 3: { return t_vgap_ksigma_regge_forward(mech, deg_k, deg_kstar, pol_ph, sp_out, sp_in, Q_sqr, s, t); break; } } return Complex(0.,0.); } /*************************************************************************/ /*************************************************************************/ /* The function klam_hel_spin_diff_cross_lab() returns the differential cross sections in the Lab system for the p(e, e'Lambda)K reaction for an electron of definite helicity and a Lambda of definite spin. The cross section dsig/dQ2dnu/dt is expressed in nanobarns / ( GeV^5). */ double klam_hel_spin_diff_cross_lab(int mech,int deg_k, int deg_kstar, double el_hel, double sp_out, double k_in_mom_lab, double k_out_mom_lab, double k_out_th_lab, double kaon_th, double kaon_phi) { vec k_in, k_out, p_in, p_out, q_in, kaon_out; double ga_en_in_lab, ga_mom_in_lab, sp_in, v, gamma, unit_conv, p_out_mom, p_out_en, s, Q_sqr, ga_en_in, ga_mom_in, k_in_x_lab, k_in_z_lab, k_out_x_lab, k_out_z_lab, k_in_x, k_in_z, k_out_x, k_out_z, kaon_mom_cm, sp_out2; double tmat = 0., t_matrix; ga_en_in_lab = k_in_mom_lab - k_out_mom_lab; ga_mom_in_lab = sqrt(pow(k_in_mom_lab, 2.) + pow(k_out_mom_lab, 2.) - 2. * k_in_mom_lab * k_out_mom_lab * cos(k_out_th_lab) ); Q_sqr = 4.*k_in_mom_lab*k_out_mom_lab * pow(sin(0.5 * k_out_th_lab),2.); s = pow(NU_M, 2.) - Q_sqr + 2. * NU_M * ga_en_in_lab; ga_en_in = (NU_M * ga_en_in_lab - Q_sqr) / sqrt(s); ga_mom_in = NU_M / sqrt(s) * ga_mom_in_lab; kaon_mom_cm = sqrt( ( pow(s - pow(LAM_M,2.) - pow(K_M,2.), 2.) - pow(2. * LAM_M * K_M, 2.) ) / (4. * s) ); k_in_x_lab = k_in_mom_lab*k_out_mom_lab/ga_mom_in_lab * sin(k_out_th_lab); k_out_x_lab = k_in_x_lab; k_in_z_lab = sqrt(pow(k_in_mom_lab, 2.) - pow(k_in_x_lab, 2.)); k_out_z_lab = k_in_z_lab - ga_mom_in_lab; v = ga_mom_in_lab / (NU_M + k_in_mom_lab - k_out_mom_lab); gamma = 1. / sqrt(1. - pow(v, 2.)); k_in_x = k_in_x_lab; k_out_x = k_out_x_lab; k_in_z = gamma * (k_in_z_lab - v * k_in_mom_lab); k_out_z = gamma * (k_out_z_lab - v * k_out_mom_lab); k_in = Vini(Complex(k_in_x,0.),Complex(0.,0.),Complex(k_in_z,0.)); k_out = Vini(Complex(k_out_x,0.),Complex(0.,0.),Complex(k_out_z,0.)); q_in = Vsub(k_in,k_out); p_in = Vsub(k_out, k_in); kaon_out = Vpolar(kaon_mom_cm, kaon_th, kaon_phi); p_out = VRmul(kaon_out, - 1.); if (targrec==1) { sp_in=sp_out; for( sp_out2 = -.5; sp_out2 <= .5; sp_out2 += 1. ) { t_matrix = klam_t_matrix(mech, deg_k, deg_kstar, el_hel, sp_in, sp_out2, kaon_phi, k_in, k_out, p_in, p_out, kaon_out); tmat += t_matrix; } } if (targrec==2) { for( sp_in = -.5; sp_in <= .5; sp_in += 1. ) { t_matrix = klam_t_matrix(mech, deg_k, deg_kstar, el_hel, sp_in, sp_out, kaon_phi, k_in, k_out, p_in, p_out, kaon_out); tmat += t_matrix; } } tmat=.5 * tmat; unit_conv = pow(.197, 2.) * pow(10., 4.); return unit_conv / pow(2. * PI, 4.) / (64. * NU_M) / ga_mom_in / sqrt(s) * (k_out_mom_lab / k_in_mom_lab ) / 2. / k_in_mom_lab / k_out_mom_lab * tmat; } double klam_t_matrix(int mech, int deg_k, int deg_kstar, double el_hel, double sp_in, double sp_out, double phi, vec k_in, vec k_out, vec p_in, vec p_out, vec kaon_out) { int nu; dcomplex nucleon_part, electron_part, e_times_p, e_times_p_incr; double k_in_mom, k_out_mom, p_in_en, p_out_en, kaon_out_en, ga_en_in, s, Q_sqr, t, e_times_p_mod2, result; fvec p4_in, p4_out, k4_in, k4_out, kaon4_out, q4_in; spinor spinor_e_in, spinor_e_out; e_times_p = Complex(0., 0.); k_in_mom = Vnorm(k_in); k_out_mom = Vnorm(k_out); p_in_en = sqrt(pow(Vnorm(p_in), 2.) + pow(NU_M, 2.)); p_out_en = sqrt(pow(Vnorm(p_out), 2.) + pow(LAM_M, 2.)); kaon_out_en = sqrt(pow(Vnorm(kaon_out), 2.) + pow(K_M, 2.)); k4_in = V4ini(Complex(k_in_mom, 0.0), k_in); k4_out = V4ini(Complex(k_out_mom, 0.0), k_out); p4_in = V4ini(Complex(p_in_en, 0.0), p_in); p4_out = V4ini(Complex(p_out_en, 0.0), p_out); q4_in = V4sub(k4_in, k4_out); kaon4_out = V4ini(Complex(kaon_out_en, 0.0), kaon_out); t = V4mul(V4sub(p4_out, p4_in),V4sub(p4_out, p4_in)).re; Q_sqr = CRmul(-1., V4mul(q4_in, q4_in)).re; s = V4mul(V4add(p4_in,q4_in),V4add(p4_in,q4_in)).re; spinor_e_in = spinor_hel(k_in, el_hel); spinor_e_out = spinor_hel(k_out, el_hel); for(nu = 0; nu <= 3; nu++) { electron_part = spinleft_mat_spinright(spinor_e_out, dirac_gamma(nu), spinor_e_in); nucleon_part = t_vgap_klam_regge_phi(mech, deg_k, deg_kstar, nu, sp_out, sp_in, Q_sqr, s, t, phi); e_times_p_incr = CRmul(metric(nu, nu), Cmul(electron_part, nucleon_part) ); /* e_times_p_incr = Cmul(electron_part, nucleon_part);*/ e_times_p = Cadd(e_times_p, e_times_p_incr); } e_times_p_mod2 = pow(Cabs(e_times_p),2.); result = pow(sqrt(k_in_mom * k_out_mom) * ELEC / Q_sqr,2.) * e_times_p_mod2 ; return result; } /*************************************************************************/ /*************************************************************************/ /* */ /* The following functions define the hadron REGGE TRAJECTORIES. */ /* */ /*************************************************************************/ /*************************************************************************/ dcomplex kaon_regge_trajectory(int deg_k, double s, double t) { int choice_alpha; double alpha_lin, alpha, alpha_sat; dcomplex traj; choice_alpha = 1; switch(choice_alpha) { case 1: /* linear trajectory */ { alpha = ALPHAp_K * (t - pow(KM_M, 2.)); break; } } switch(deg_k) /* 2 = default */ { case 1: /* non-degenerate trajectory */ { traj = CRmul(pow(s,alpha) / gamma_function(alpha + 1.) * PI * ALPHAp_K / sin(PI * alpha), CRmul(0.5, Complex(1. + cos(PI * alpha),- sin(PI * alpha)) ) ); break; } case 2: /* degenerate trajectory : rotating phase */ { traj = CRmul(pow(s,alpha) / gamma_function(alpha + 1.) * PI * ALPHAp_K / sin(PI * alpha), Complex(cos(PI * alpha), - sin(PI * alpha)) ); break; } case 3: /* degenerate trajectory : non-rotating phase */ { traj = CRmul(pow(s,alpha) / gamma_function(alpha + 1.) * PI * ALPHAp_K / sin(PI * alpha), Complex(1., 0.) ); break; } case 4: /* single pole propagator */ { traj = Complex(1./(t - pow(K_M, 2.) ), 0.); break; } } return traj; } dcomplex kstar_regge_trajectory(int deg_kstar, double s, double t) { int choice_alpha; double alpha_lin, alpha_sat, alpha; dcomplex traj; choice_alpha = 2; switch(choice_alpha) { case 1: /* linear trajectory which passes exactly through kstar pole */ { alpha = 1. + ALPHAp_KSTAR * (t - pow(KSTAR_M, 2.)); break; } case 2: /* linear trajectory */ { alpha = 0.25 + ALPHAp_KSTAR * t; break; } } switch(deg_kstar) /* 2 = default */ { case 1: /* non-degenerate trajectory */ { traj = CRmul(pow(s,alpha - 1.) / gamma_function(alpha) * 0.5 * PI * ALPHAp_KSTAR / sin(PI * alpha), Complex(-1. + cos(PI * alpha), - sin(PI * alpha)) ); break; } case 2: /* degenerate trajectory : rotating phase */ { traj = CRmul(pow(s,alpha - 1.) / gamma_function(alpha) * PI * ALPHAp_KSTAR / sin(PI * alpha), Complex(cos(PI * alpha), - sin(PI * alpha)) ); break; } case 3: /* degenerate trajectory : non-rotating phase */ { traj = CRmul(pow(s,alpha - 1.) / gamma_function(alpha) * PI * ALPHAp_KSTAR / sin(PI * alpha), Complex(- 1., 0.) ); break; } case 4: /* single pole propagator */ { double gamma_K=.0508; traj = Complex(1./(t - pow(KSTAR_M, 2.) ), 0.); /* traj=Complex(t-pow(KSTAR_M,2.),-KSTAR_M*gamma_K); traj=CRmul(1./(pow(t-pow(KSTAR_M,2.),2.)+pow(KSTAR_M*gamma_K,2.)),traj); */ break; } } return traj; } /*************************************************************************/ /*************************************************************************/ /* */ /* The following functions define the FORM FACTORS that are needed. */ /* */ /*************************************************************************/ /*************************************************************************/ double kaon_emformfactor(double Q_sqr, double t) { switch(ichoicek) { case 1 : /* monopole form with mass scale as fit parameter */ { return 1./(1. + Q_sqr/mscalk); break; } case 2 : /* form with mass scale as fit parameter */ { double a; a=2./mscalk-1.46; return (1. + a*Q_sqr) / pow((1.+ Q_sqr/mscalk),2.); break; } case 3 : /* VMD (rho, omega, phi) form */ { double rho_mass_sqr, omega_mass_sqr, phi_mass_sqr; rho_mass_sqr = pow(0.770, 2.); omega_mass_sqr = pow(0.782, 2.); phi_mass_sqr = pow(1.019, 2.); return 1./2. * 1./(1. + Q_sqr/rho_mass_sqr) + 1./6. * 1./(1. + Q_sqr/omega_mass_sqr) + 1./3. * 1./(1. + Q_sqr/phi_mass_sqr); break; } case 4 : /* FF with off-shell parametrization */ { return 1./(1. + Q_sqr/(.7 + c_offshell * (-t + pow(K_M,2.))) ); break; } } } double kstarkgamma_formfactor(double Q_sqr, double t) { switch(ichoicekst) { case 1 : /* monopole form with mass scale as fit parameter */ { return 1./(1. + Q_sqr/mscalkst); break; } case 2 : /* monopole form with mass scale as fit parameter */ { double a; a=2./mscalkst-1.46; return (1.+a*Q_sqr) / pow((1.+ Q_sqr/mscalkst),2.); break; } case 3 : /* VMD (rho, omega, phi) form */ { double rho_mass_sqr, omega_mass_sqr, phi_mass_sqr; rho_mass_sqr = pow(0.770, 2.); omega_mass_sqr = pow(0.782, 2.); phi_mass_sqr = pow(1.019, 2.); return 1./2. * 1./(1. + Q_sqr/rho_mass_sqr) + 1./6. * 1./(1. + Q_sqr/omega_mass_sqr) + 1./3. * 1./(1. + Q_sqr/phi_mass_sqr); break; } } } /*************************************************************************/ /*************************************************************************/ /* */ /* The following amplitudes are for the (e, e' K Lambda) reaction. */ /* */ /*************************************************************************/ /*************************************************************************/ dcomplex t_vgap_klam_regge_phi(int mech, int deg_k, int deg_kstar, int nu, double sp_out, double sp_in, double Q_sqr, double s, double t, double phi) { int mu, lambda, alpha, beta, off_shell; double ga_mom_cm, ga_en_cm, k_mom_cm, k_en, k_th, cos_th, po_s, E_in, E_out, norm_in, norm_out; dcomplex contact, kaon, gauge_corr, kaon_gaugeinv, temp_result, kstarkga; vec pol_ph_left, pol_ph_right, pol_ph_long, k_ga, p_k, p_in, p_out, pi_s; fvec k4_ga, p4_k, p_s, diff, mom, pol4_ph; spinor spinor_p_in, spinor_p_out; spinor_mat contact_vertex_times_eps4, knn_vertex, dirac_sN, pauli_sN, sNU_vertex, kr_vertex, pauli, pauli_incr, kstarNN_vertex, kstarNN_matrix; ga_mom_cm = sqrt( ( pow(s - pow(NU_M,2.) + Q_sqr, 2.) + pow(2. * NU_M, 2.) * Q_sqr ) / (4. * s) ); ga_en_cm = sqrt(s) - sqrt(pow(ga_mom_cm,2.) + pow(NU_M,2.)); k_ga = Vpolar(ga_mom_cm, 0., 0.); k4_ga = V4ini(Complex(ga_en_cm,0.), k_ga); k_mom_cm = sqrt( ( pow(s - pow(LAM_M,2.) - pow(K_M,2.), 2.) - pow(2. * LAM_M * K_M, 2.) ) / (4. * s) ); k_en = sqrt( pow(k_mom_cm, 2.) + pow(K_M, 2.) ); cos_th = (t + Q_sqr - pow(K_M,2.) + 2. * ga_en_cm * k_en) / (2. * ga_mom_cm * k_mom_cm); if (cos_th > 1.) cos_th = 1.; if (cos_th < -1.) cos_th = -1.; if ((cos_th < -1.) || (cos_th > 1.)) return Complex(0., 0.); k_th = acos(cos_th); p_k = Vpolar(k_mom_cm, k_th, phi); p4_k = V4ini(Complex(k_en,0.), p_k); p_in = VRmul(k_ga, - 1.0); p_out = VRmul(p_k,- 1.0); E_in = sqrt( pow( Vnorm(k_ga), 2.) + pow(NU_M,2.) ); E_out = sqrt( pow( Vnorm(p_k), 2.) + pow(LAM_M,2.) ); pi_s = Vadd(p_out, p_k); po_s = sqrt(s); p_s = V4ini(Complex(po_s, 0.), pi_s); if (isys==1) { spinor_p_out = spinor_pos(p_out, LAM_M, sp_out); spinor_p_in = spinor_pos(p_in, NU_M, sp_in); } if (isys==2) { spinor_p_out = spinor_part_hel(p_out, LAM_M, sp_out); spinor_p_in = spinor_part_hel(p_in, NU_M, sp_in); } if (isys==3) { spinor_p_out = spinor_pos_xaxis(p_out, LAM_M, sp_out); spinor_p_in = spinor_pos_xaxis(p_in, NU_M, sp_in); } if (isys==4) { spinor_p_out = spinor_pos_yaxis(p_out, LAM_M, sp_out); spinor_p_in = spinor_pos_yaxis(p_in, NU_M, sp_in); } norm_in = sqrt(E_in + NU_M); norm_out = sqrt(E_out + LAM_M); diff = V4sub(k4_ga, p4_k); mom = V4sub(V4Rmul(p4_k, 2.), k4_ga); /**************************************/ /* */ /* KAON exchange contribution */ /* */ /**************************************/ knn_vertex = SMCmul(dirac_gamma5(), Complex(G_KNLAM, 0.) ); kaon = Cmul(Cmul(Complex(0., norm_in * norm_out * ELEC), Cmul(kaon_regge_trajectory(deg_k, s, t), mom.fv[nu]) ), spinleft_mat_spinright(spinor_p_out, knn_vertex, spinor_p_in) ); off_shell = 1; switch(off_shell) { case 1 : /* PS coupling for the gauge correction terms */ { dirac_sN = zero_spinor_mat(); pauli_sN = zero_spinor_mat(); dirac_sN = dirac_gamma(nu); sNU_vertex = SMCmul(SMmul(dirac_gamma5(), SMmul(SMadd(fvec_slash(p_s), SMCmul(unit_spinor_mat(), Complex(NU_M, 0.) ) ), SMadd(dirac_sN, pauli_sN) ) ), Complex(G_KNLAM * (t - pow(KM_M,2.))/(s - pow(NU_M,2.)) , 0.) ); gauge_corr = Cmul(Cmul(Complex(0., norm_in * norm_out * ELEC), kaon_regge_trajectory(deg_k, s, t) ), spinleft_mat_spinright(spinor_p_out, sNU_vertex, spinor_p_in) ); break; } } kaon_gaugeinv = Cadd(kaon, gauge_corr); /**************************************/ /* */ /* KSTAR exchange contribution */ /* */ /**************************************/ if ((mech != 1) && (mech != 3)) { kstarNN_matrix = zero_spinor_mat(); for(beta = 0; beta <= 3; beta++) { pauli = zero_spinor_mat(); for(lambda = 0; lambda <= 3; lambda++) { pauli_incr = SMCmul(dirac_sigma(beta, lambda), Cmul(diff.fv[lambda], Complex(0.0,metric(lambda,lambda) * KAPPA_KSTARLAM /(NU_M + LAM_M)) ) ); pauli = SMadd(pauli, pauli_incr); } kstarNN_vertex = SMCmul(SMadd(dirac_gamma(beta), pauli), Complex(G_KSTARNLAM, 0.) ); temp_result = Complex(0.,0.); for(mu = 0; mu <= 3; mu++) { for(alpha = 0; alpha <= 3; alpha++) { temp_result = Cadd(temp_result, CRmul(antisymm_epsilon(mu,nu,alpha,beta), CRmul(metric(nu,nu), Cmul(k4_ga.fv[mu], diff.fv[alpha]))) ); } } kstarNN_matrix = SMadd(kstarNN_matrix, SMCmul(kstarNN_vertex, temp_result) ); } kstarkga = Cmul(CRmul( norm_in * norm_out * G_KSTARKGA * ELEC , kstar_regge_trajectory(deg_kstar, s, t) ), spinleft_mat_spinright(spinor_p_out, kstarNN_matrix, spinor_p_in) ); } switch(mech) { case 1 : /* kaon exchange */ { double ff_kaon; ff_kaon = kaon_emformfactor(Q_sqr, t); return CRmul(ff_kaon, kaon); break; } case 2 : /* kstar exchange */ { double ff_kstar; ff_kstar = kstarkgamma_formfactor(Q_sqr, t); return CRmul(ff_kstar, kstarkga); break; } case 3 : /* gauge invariant kaon exchange */ { double ff_kaon; ff_kaon = kaon_emformfactor(Q_sqr, t); return CRmul(ff_kaon, kaon_gaugeinv); break; } case 4 : /* gauge invariant kaon exchange + kstar exchange */ { double ff_kaon, ff_kstar; ff_kaon = kaon_emformfactor(Q_sqr, t); ff_kstar = kstarkgamma_formfactor(Q_sqr, t); return Cadd(CRmul(ff_kaon, kaon_gaugeinv), CRmul(ff_kstar, kstarkga) ); break; } } return Complex(0.,0.); } /*************************************************************************/ /*************************************************************************/ /* */ /* The following amplitudes are for the (e, e' K Lambda) reaction. */ /* */ /*************************************************************************/ /*************************************************************************/ dcomplex t_vgap_klam_regge_forward(int mech, int deg_k, int deg_kstar, int pol_ph, double sp_out, double sp_in, double Q_sqr, double s, double t) { int mu, lambda, nu, alpha, beta, off_shell; double ga_mom_cm, ga_en_cm, k_mom_cm, k_en, k_th, cos_th, po_s, E_in, E_out, norm_in, norm_out; dcomplex contact, kaon, gauge_corr, kaon_gaugeinv, temp_result, kstarkga, mom_times_eps4; vec pol_ph_left, pol_ph_right, pol_ph_long, k_ga, p_k, p_in, p_out, pi_s; fvec k4_ga, p4_k, p_s, diff, mom, pol4_ph; spinor spinor_p_in, spinor_p_out; spinor_mat contact_vertex_times_eps4, knn_vertex, dirac_sN, pauli_sN, sNU_vertex, kr_vertex, pauli, pauli_incr, kstarNN_vertex, kstarNN_matrix; ga_mom_cm = sqrt( ( pow(s - pow(NU_M,2.) + Q_sqr, 2.) + pow(2. * NU_M, 2.) * Q_sqr ) / (4. * s) ); ga_en_cm = sqrt(s) - sqrt(pow(ga_mom_cm,2.) + pow(NU_M,2.)); k_ga = Vpolar(ga_mom_cm, 0., 0.); k4_ga = V4ini(Complex(ga_en_cm,0.), k_ga); k_mom_cm = sqrt( ( pow(s - pow(LAM_M,2.) - pow(K_M,2.), 2.) - pow(2. * LAM_M * K_M, 2.) ) / (4. * s) ); k_en = sqrt( pow(k_mom_cm, 2.) + pow(K_M, 2.) ); cos_th = (t + Q_sqr - pow(K_M,2.) + 2. * ga_en_cm * k_en) / (2. * ga_mom_cm * k_mom_cm); if (cos_th > 1.) cos_th = 1.; if (cos_th < -1.) cos_th = -1.; /* if ((cos_th < -1.) || (cos_th > 1.)) return Complex(0., 0.); */ k_th = acos(cos_th); p_k = Vpolar(k_mom_cm, k_th, 0.); p4_k = V4ini(Complex(k_en,0.), p_k); p_in = VRmul(k_ga, - 1.0); p_out = VRmul(p_k,- 1.0); E_in = sqrt( pow( Vnorm(k_ga), 2.) + pow(NU_M,2.) ); E_out = sqrt( pow( Vnorm(p_k), 2.) + pow(LAM_M,2.) ); pi_s = Vadd(p_out, p_k); po_s = sqrt(s); p_s = V4ini(Complex(po_s, 0.), pi_s); /* spinor_p_out = spinor_part_hel(p_out, LAM_M, sp_out); spinor_p_in = spinor_part_hel(p_in, NU_M, sp_in); */ spinor_p_out = spinor_pos(p_out, LAM_M, sp_out); spinor_p_in = spinor_pos(p_in, NU_M, sp_in); norm_in = sqrt(E_in + NU_M); norm_out = sqrt(E_out + LAM_M); diff = V4sub(k4_ga, p4_k); mom = V4sub(V4Rmul(p4_k, 2.), k4_ga); switch(pol_ph) { case -1 : { pol_ph_left = Vini( Complex(1./ sqrt(2.), 0.), Complex(0., - 1./sqrt(2.)), Complex(0., 0.) ); pol4_ph = V4ini(Complex(0.,0.), pol_ph_left); break; } case 0 : { pol_ph_long = Vini( Complex(0., 0.), Complex(0., 0.), Complex(ga_en_cm/sqrt(Q_sqr), 0.) ); pol4_ph = V4ini(Complex(ga_mom_cm/sqrt(Q_sqr), 0.), pol_ph_long); /* pol4_ph = V4ini(Complex(0., 0.), Vini(Complex(0.,0.),Complex(0.,0.), Complex(sqrt(Q_sqr)/ga_en_cm,0.)) ); */ break; } case 1 : { pol_ph_right = Vini( Complex(- 1./ sqrt(2.), 0.), Complex(0., - 1./sqrt(2.)), Complex(0., 0.) ); pol4_ph = V4ini(Complex(0., 0.), pol_ph_right); break; } } /* use the following to check the gauge invariance */ /* pol4_ph = k4_ga; */ /**************************************/ /* */ /* KAON exchange contribution */ /* */ /**************************************/ mom_times_eps4 = V4mul(mom, pol4_ph); knn_vertex = SMCmul(dirac_gamma5(), Complex(G_KNLAM, 0.) ); kaon = Cmul(Cmul(Complex(0., norm_in * norm_out * ELEC), Cmul(kaon_regge_trajectory(deg_k, s, t), mom_times_eps4) ), spinleft_mat_spinright(spinor_p_out, knn_vertex, spinor_p_in) ); off_shell = 1; switch(off_shell) { case 1 : /* PS coupling for the gauge correction terms */ { dirac_sN = zero_spinor_mat(); pauli_sN = zero_spinor_mat(); for(mu = 0; mu <= 3; mu++) dirac_sN = SMadd(dirac_sN, SMCmul(dirac_gamma(mu), CRmul(metric(mu,mu), pol4_ph.fv[mu]) ) ); sNU_vertex = SMCmul(SMmul(dirac_gamma5(), SMmul(SMadd(fvec_slash(p_s), SMCmul(unit_spinor_mat(), Complex(NU_M, 0.) ) ), SMadd(dirac_sN, pauli_sN) ) ), Complex(G_KNLAM * (t - pow(KM_M,2.))/(s - pow(NU_M,2.)) , 0.) ); gauge_corr = Cmul(Cmul(Complex(0., norm_in * norm_out * ELEC), kaon_regge_trajectory(deg_k, s, t) ), spinleft_mat_spinright(spinor_p_out, sNU_vertex, spinor_p_in) ); break; } } kaon_gaugeinv = Cadd(kaon, gauge_corr); /**************************************/ /* */ /* KSTAR exchange contribution */ /* */ /**************************************/ if ((mech != 1) && (mech != 3)) { kstarNN_matrix = zero_spinor_mat(); for(beta = 0; beta <= 3; beta++) { pauli = zero_spinor_mat(); for(lambda = 0; lambda <= 3; lambda++) { pauli_incr = SMCmul(dirac_sigma(beta, lambda), Cmul(diff.fv[lambda], Complex(0.0,metric(lambda,lambda) * KAPPA_KSTARLAM /(NU_M + LAM_M)) ) ); /* pauli_incr = SMCmul(dirac_sigma(beta, lambda), Cmul(diff.fv[lambda], Complex(0.0,metric(lambda,lambda) * KAPPA_KSTARLAM /(NU_M + NU_M)) ) ); */ pauli = SMadd(pauli, pauli_incr); } kstarNN_vertex = SMCmul(SMadd(dirac_gamma(beta), pauli), Complex(G_KSTARNLAM, 0.) ); temp_result = Complex(0.,0.); for(mu = 0; mu <= 3; mu++) { for(nu = 0; nu <= 3; nu++) { for(alpha = 0; alpha <= 3; alpha++) { temp_result = Cadd(temp_result, CRmul(antisymm_epsilon(nu,mu,alpha,beta), Cmul(Cmul(k4_ga.fv[nu], diff.fv[alpha]), pol4_ph.fv[mu]) ) ); } } } /* for(mu = 0; mu <= 3; mu++) {if (pol_ph==0) printf("\n %f %f %f",k4_ga.fv[mu].re,diff.fv[mu].re, pol4_ph.fv[mu].re);}*/ kstarNN_matrix = SMadd(kstarNN_matrix, SMCmul(kstarNN_vertex, temp_result) ); } kstarkga = Cmul(CRmul( norm_in * norm_out * G_KSTARKGA * ELEC , kstar_regge_trajectory(deg_kstar, s, t) ), spinleft_mat_spinright(spinor_p_out, kstarNN_matrix, spinor_p_in) ); } /* printf("\n %e \t", Q_sqr); Cecho(kstarkga); */ switch(mech) { case 1 : /* kaon exchange */ { double ff_kaon; ff_kaon = kaon_emformfactor(Q_sqr, t); return CRmul(ff_kaon, kaon); break; } case 2 : /* kstar exchange */ { double ff_kstar; ff_kstar = kstarkgamma_formfactor(Q_sqr, t); return CRmul(ff_kstar, kstarkga); break; } case 3 : /* gauge invariant kaon exchange */ { double ff_kaon; ff_kaon = kaon_emformfactor(Q_sqr, t); return CRmul(ff_kaon, kaon_gaugeinv); break; } case 4 : /* gauge invariant kaon exchange + kstar exchange */ { double ff_kaon, ff_kstar; ff_kaon = kaon_emformfactor(Q_sqr, t); ff_kstar = kstarkgamma_formfactor(Q_sqr, t); return Cadd(CRmul(ff_kaon, kaon_gaugeinv), CRmul(ff_kstar, kstarkga) ); break; } } return Complex(0.,0.); } /*************************************************************************/ /*************************************************************************/ /* */ /* The following t matrices are for the (e, e' K Sigma) reactions. */ /* */ /*************************************************************************/ /*************************************************************************/ dcomplex t_vgap_ksigma_regge_forward(int mech, int deg_k, int deg_kstar, int pol_ph, double sp_out, double sp_in, double Q_sqr, double s, double t) { int mu, lambda, nu, alpha, beta, off_shell; double ga_mom_cm, ga_en_cm, k_mom_cm, k_en, k_th, cos_th, po_s, E_in, E_out, norm_in, norm_out; dcomplex contact, kaon, gauge_corr, kaon_gaugeinv, temp_result, kstarkga, mom_times_eps4; vec pol_ph_left, pol_ph_right, pol_ph_long, k_ga, p_k, p_in, p_out, pi_s; fvec k4_ga, p4_k, p_s, diff, mom, pol4_ph; spinor spinor_p_in, spinor_p_out; spinor_mat contact_vertex_times_eps4, knn_vertex, dirac_sN, pauli_sN, sNU_vertex, kr_vertex, pauli, pauli_incr, kstarNN_vertex, kstarNN_matrix; ga_mom_cm = sqrt( ( pow(s - pow(NU_M,2.) + Q_sqr, 2.) + pow(2. * NU_M, 2.) * Q_sqr ) / (4. * s) ); ga_en_cm = sqrt(s) - sqrt(pow(ga_mom_cm,2.) + pow(NU_M,2.)); k_ga = Vpolar(ga_mom_cm, 0., 0.); k4_ga = V4ini(Complex(ga_en_cm,0.), k_ga); k_mom_cm = sqrt( ( pow(s - pow(SIG_M,2.) - pow(K_M,2.), 2.) - pow(2. * SIG_M * K_M, 2.) ) / (4. * s) ); k_en = sqrt( pow(k_mom_cm, 2.) + pow(K_M, 2.) ); cos_th = (t + Q_sqr - pow(K_M,2.) + 2. * ga_en_cm * k_en) / (2. * ga_mom_cm * k_mom_cm); if (cos_th > 1.) cos_th = 1.; if (cos_th < -1.) cos_th = -1.; if ((cos_th < -1.) || (cos_th > 1.)) return Complex(0., 0.); k_th = acos(cos_th); p_k = Vpolar(k_mom_cm, k_th, 0.); p4_k = V4ini(Complex(k_en,0.), p_k); p_in = VRmul(k_ga, - 1.0); p_out = VRmul(p_k,- 1.0); E_in = sqrt( pow( Vnorm(k_ga), 2.) + pow(NU_M,2.) ); E_out = sqrt( pow( Vnorm(p_k), 2.) + pow(SIG_M,2.) ); pi_s = Vadd(p_out, p_k); po_s = sqrt(s); p_s = V4ini(Complex(po_s, 0.), pi_s); spinor_p_out = spinor_pos(p_out, SIG_M, sp_out); spinor_p_in = spinor_pos(p_in, NU_M, sp_in); norm_in = sqrt(E_in + NU_M); norm_out = sqrt(E_out + SIG_M); diff = V4sub(k4_ga, p4_k); mom = V4sub(V4Rmul(p4_k, 2.), k4_ga); switch(pol_ph) { case -1 : { pol_ph_left = Vini( Complex(1./ sqrt(2.), 0.), Complex(0., - 1./sqrt(2.)), Complex(0., 0.) ); pol4_ph = V4ini(Complex(0.,0.), pol_ph_left); break; } case 0 : { pol_ph_long = Vini( Complex(0., 0.), Complex(0., 0.), Complex(ga_en_cm/sqrt(Q_sqr), 0.) ); pol4_ph = V4ini(Complex(ga_mom_cm/sqrt(Q_sqr), 0.), pol_ph_long); /* pol4_ph = V4ini(Complex(0., 0.), Vini(Complex(0.,0.),Complex(0.,0.), Complex(sqrt(Q_sqr)/ga_en_cm,0.)) ); */ break; } case 1 : { pol_ph_right = Vini( Complex(- 1./ sqrt(2.), 0.), Complex(0., - 1./sqrt(2.)), Complex(0., 0.) ); pol4_ph = V4ini(Complex(0., 0.), pol_ph_right); break; } } /**************************************/ /* */ /* KAON exchange contribution */ /* */ /**************************************/ dirac_sN = zero_spinor_mat(); pauli_sN = zero_spinor_mat(); mom_times_eps4 = V4mul(mom, pol4_ph); knn_vertex = SMCmul(dirac_gamma5(), Complex(G_KNSIG, 0.) ); kaon = Cmul(Cmul(Complex(0., norm_in * norm_out * ELEC), Cmul(kaon_regge_trajectory(deg_k, s, t), mom_times_eps4) ), spinleft_mat_spinright(spinor_p_out, knn_vertex, spinor_p_in) ); off_shell = 1; switch(off_shell) { case 1 : /* PS coupling for the gauge correction terms */ { dirac_sN = zero_spinor_mat(); pauli_sN = zero_spinor_mat(); for(mu = 0; mu <= 3; mu++) dirac_sN = SMadd(dirac_sN, SMCmul(dirac_gamma(mu), CRmul(metric(mu,mu), pol4_ph.fv[mu]) ) ); sNU_vertex = SMCmul(SMmul(dirac_gamma5(), SMmul(SMadd(fvec_slash(p_s), SMCmul(unit_spinor_mat(), Complex(NU_M, 0.) ) ), SMadd(dirac_sN, pauli_sN) ) ), Complex(G_KNSIG * (t - pow(KM_M,2.))/(s - pow(NU_M,2.)) , 0.) ); gauge_corr = Cmul(Cmul(Complex(0., norm_in * norm_out * ELEC), kaon_regge_trajectory(deg_k, s, t) ), spinleft_mat_spinright(spinor_p_out, sNU_vertex, spinor_p_in) ); break; } } kaon_gaugeinv = Cadd(kaon, gauge_corr); /**************************************/ /* */ /* KSTAR exchange contribution */ /* */ /**************************************/ if ((mech != 1) && (mech != 3)) { kstarNN_matrix = zero_spinor_mat(); for(beta = 0; beta <= 3; beta++) { pauli = zero_spinor_mat(); for(lambda = 0; lambda <= 3; lambda++) { pauli_incr = SMCmul(dirac_sigma(beta, lambda), Cmul(diff.fv[lambda], Complex(0.0,metric(lambda,lambda) * KAPPA_KSTARSIG /(NU_M + SIG_M)) ) ); pauli = SMadd(pauli, pauli_incr); } kstarNN_vertex = SMCmul(SMadd(dirac_gamma(beta), pauli), Complex(G_KSTARNSIG, 0.) ); temp_result = Complex(0.,0.); for(mu = 0; mu <= 3; mu++) { for(nu = 0; nu <= 3; nu++) { for(alpha = 0; alpha <= 3; alpha++) { temp_result = Cadd(temp_result, CRmul(antisymm_epsilon(nu,mu,alpha,beta), Cmul(Cmul(k4_ga.fv[nu], diff.fv[alpha]), pol4_ph.fv[mu] ) ) ); } } } /* for(mu = 0; mu <= 3; mu++) {if (pol_ph==0) printf("\n %f %f %f",k4_ga.fv[mu].re,diff.fv[mu].re, pol4_ph.fv[mu].re);} */ kstarNN_matrix = SMadd(kstarNN_matrix, SMCmul(kstarNN_vertex, temp_result) ); } kstarkga = Cmul(CRmul( norm_in * norm_out * G_KSTARKGA * ELEC , kstar_regge_trajectory(deg_kstar, s, t) ), spinleft_mat_spinright(spinor_p_out, kstarNN_matrix, spinor_p_in) ); if (k0sig_s==1) kstarkga=CRmul(sqrt(2.3/.99),kstarkga); } switch(mech) { case 1 : /* kaon exchange */ { double ff_kaon; ff_kaon = kaon_emformfactor(Q_sqr, t); return CRmul(ff_kaon, kaon); break; } case 2 : /* kstar exchange */ { double ff_kstar; if (k0sig_s==1) ff_kstar=1. ; else ff_kstar = kstarkgamma_formfactor(Q_sqr, t); return CRmul(ff_kstar,kstarkga); break; } case 3 : /* gauge invariant kaon exchange */ { double ff_kaon; ff_kaon = kaon_emformfactor(Q_sqr, t); return CRmul(ff_kaon, kaon_gaugeinv); break; } case 4 : /* gauge invariant kaon exchange + kstar exchange */ { double ff_kaon, ff_kstar; ff_kaon = kaon_emformfactor(Q_sqr, t); ff_kstar = kstarkgamma_formfactor(Q_sqr, t); return Cadd(CRmul(ff_kaon, kaon_gaugeinv), CRmul(ff_kstar, kstarkga) ); break; } } return Complex(0.,0.); } /*************************************************************************/