src/eic_evgen/KPlus_sig_Scaling.cc

0001 // 31/01/23 - SJDK
0002 // New file to get the KPlus cross section via the scaling method that Ali Usman utilised in preliminary kaon studies
0003 // This method relies upon getting the pion cross section before subsequently scaling it
0004
0005 #include "KPlus_sig_Scaling.h"
0006 #include "PiPlus_sig.h"
0007 #include "PiPlus_sig_Param.h"
0008 #include "reaction_routine.h"
0009 #include "TF1.h"
0010 #include "TString.h"
0011
0012 double GetKPlus_CrossSection_Scaling(double ft, double fw, double fqsq, double feps, double meson_mass, TString recoil_hadron){
0013
0014   double sig_total;
0015
0016   // --------------------------------------------------------------------------------------------------
0017   // CKY sigma L and T starts
0018   // --------------------------------------------------------------------------------------------------
0019   double lpar0 = 0., lpar1 = 0., lpar2 = 0., lpar3 = 0., lpar4 = 0., lpar5 = 0., lpar6 = 0.;
0020   double tpar0 = 0., tpar1 = 0., tpar2 = 0., tpar3 = 0., tpar4 = 0.;
0021
0022   lpar0 = 0.;    lpar1 = 0.;    lpar2 = 0.;    lpar3 = 0.;    lpar4 = 0.;    lpar5 = 0.;    lpar6 = 0.;
0023   tpar0 = 0.;    tpar1 = 0.;    tpar2 = 0.;    tpar3 = 0.;    tpar4 = 0.;
0024
0025   fSig_L = 0;
0026   fSig_T = 0;
0027
0028   if ( ( ft > 0. ) && ( ft < 0.15 ) ) {
0029     PiPlus_sigmaL_Param( fw,  fqsq, lpar0, lpar1, lpar2 , lpar3 , lpar4 , lpar5 , lpar6 );
0030     TF1 *fitCKYLonglandau = new TF1("sigmaL","landau", 0.0 , 0.15 );
0031     fitCKYLonglandau->FixParameter( 0 , lpar0 );
0032     fitCKYLonglandau->FixParameter( 1 , lpar1 );
0033     fitCKYLonglandau->FixParameter( 2 , lpar2 );
0034     fSig_L = fitCKYLonglandau->Eval(ft);
0035     if ( lpar0 == 0 || lpar1 == 0 || lpar2 == 0 )
0036       fSig_L = 0;
0037     fitCKYLonglandau = NULL;
0038     delete fitCKYLonglandau;
0039   }
0040   else if ( ( ft > 0.15 ) && ( ft < 0.5 ) ) {
0041     PiPlus_sigmaL_Param( fw,  fqsq, lpar0, lpar1, lpar2 , lpar3 , lpar4 , lpar5 , lpar6 );
0042     TF1 *fitCKYLongexpo1 = new TF1("sigmaL","expo", 0.15 , 0.5 );
0043     fitCKYLongexpo1->FixParameter( 0 , lpar3 );
0044     fitCKYLongexpo1->FixParameter( 1 , lpar4 );
0045     fSig_L = fitCKYLongexpo1->Eval(ft);
0046     if ( lpar3 == 0 || lpar4 == 0 )
0047       fSig_L = 0;
0048     fitCKYLongexpo1 = NULL;
0049     delete fitCKYLongexpo1;
0050   }
0051   else if ( ( ft > 0.5 ) && ( ft < 1.3 ) ) {
0052     PiPlus_sigmaL_Param( fw,  fqsq, lpar0, lpar1, lpar2 , lpar3 , lpar4 , lpar5 , lpar6 );
0053     TF1 *fitCKYLongexpo2 = new TF1("sigmaL","expo", 0.5 , 1.3 );
0054     fitCKYLongexpo2->FixParameter( 0 , lpar5 );
0055     fitCKYLongexpo2->FixParameter( 1 , lpar6 );
0056     fSig_L = fitCKYLongexpo2->Eval(ft);
0057     if ( lpar5 == 0 || lpar6 == 0 )
0058       fSig_L = 0;
0059     fitCKYLongexpo2 = NULL;
0060     delete fitCKYLongexpo2;
0061   }
0062   else {
0063     fSig_L = 0;
0064   }
0065    // SJDK - 02/06/22 - The validity range here was inconsistent, this only went from 0.0 to 0.15, leaving a gap between 0.15 to 0.2
0066   // I changed the range to remove this gap.
0067   if ( ( ft > 0.0 ) && ( ft < 0.2 ) ) {
0068     PiPlus_sigmaT_Param( fw,  fqsq, tpar0, tpar1, tpar2 , tpar3 , tpar4 );
0069     TF1 *fitCKYTranspol2 = new TF1("sigmaL","pol2", 0.0 , 0.2 );
0070     fitCKYTranspol2->FixParameter( 0 , tpar0 );
0071     fitCKYTranspol2->FixParameter( 1 , tpar1 );
0072     fitCKYTranspol2->FixParameter( 2 , tpar2 );
0073     fSig_T = fitCKYTranspol2->Eval(ft);
0074     if ( tpar0 == 0 || tpar1 == 0 || tpar2 == 0 )
0075       fSig_T = 0;
0076     fitCKYTranspol2 = NULL;
0077     delete fitCKYTranspol2;
0078   }
0079   else if ( ( ft > 0.2 ) && ( ft < 1.3 ) ) {
0080     PiPlus_sigmaT_Param( fw,  fqsq, tpar0, tpar1, tpar2 , tpar3 , tpar4 );
0081     TF1 *fitCKYTransexpo = new TF1("sigmaL","expo", 0.2 , 1.3 );
0082     fitCKYTransexpo->FixParameter( 0 , tpar3 );
0083     fitCKYTransexpo->FixParameter( 1 , tpar4 );
0084     fSig_T = fitCKYTransexpo->Eval(ft);
0085     if ( tpar3 == 0 || tpar4 == 0 )
0086       fSig_T = 0;
0087     fitCKYTransexpo = NULL;
0088     delete fitCKYTransexpo;
0089   }
0090
0091   // ------------------------------------------------------------------------------------------------
0092   // Improving the Kaon Sigma_L following GH's fortran code (pole_ration.f) located in main driectory
0093   // Added by AU on July 10, 2021
0094   // ------------------------------------------------------------------------------------------------
0095
0096   double mkg = 0, mpig = 0, hbarc = 0, gpoleKL = 0, gpoleKS = 0, gpolepi = 0, gpoleKhyp = 0, r2_dip = 0, r2_mono = 0, Fpi_mono = 0, Fpi_dip = 0, pmono = 0, Fpi_fit = 0, fpisq = 0, Fkk = 0, fksq = 0, lNpi = 0, lNk = 0, gkhypn = 0, gpinn = 0, dl_poleKhyp = 0, dl_polepi = 0, ratio_Khyp_pi = 0, FT_GeV_neg = 0;
0097
0098   gpoleKL = -13.3;
0099   gpoleKS = -3.5;
0100   gpolepi = 13.1;
0101
0102   r2_dip = 0.411;
0103   r2_mono = 0.431;
0104
0105   hbarc = 0.197;
0106   mkg = meson_mass;
0107   mpig = 0.13957;
0108
0109   FT_GeV_neg = -1.0 * ft;
0110
0111   if ( recoil_hadron == "Lambda" ) {
0112     gpoleKhyp = gpoleKL;
0113   }
0114   else if ( recoil_hadron == "Sigma0" ) {
0115     gpoleKhyp = gpoleKS;
0116   }
0117
0118   Fpi_mono = 1.0 / ( 1.0 + ( r2_mono * fqsq) / (6 * (hbarc * hbarc)));
0119   Fpi_dip = 1.0 / ( 1.0 + ( r2_dip * fqsq) / (12 * (hbarc * hbarc))) * ( 1.0 + ( r2_dip * fqsq) / (12 * (hbarc * hbarc)));
0120   pmono = 0.85;
0121
0122   Fpi_fit = pmono * Fpi_mono + (1 - pmono) * Fpi_dip;
0123   fpisq = Fpi_fit * Fpi_fit;
0124
0125   Fkk = 0.9 / (1.0 + fqsq / 0.462 );
0126   fksq = Fkk * Fkk;
0127
0128   lNpi = 0.44;
0129   lNk = (0.44+0.80)/2;
0130
0131   gkhypn = gpoleKhyp * ((lNk * lNk) - (mkg * mkg)) / ((lNk * lNk) - FT_GeV_neg);
0132   gpinn = gpolepi * ((lNpi * lNpi) - (mpig * mpig)) / ((lNpi * lNpi) - FT_GeV_neg);
0133
0134   dl_poleKhyp = ((gkhypn * gkhypn) * fksq) / ((FT_GeV_neg - (mkg * mkg))*(FT_GeV_neg - (mkg * mkg)));
0135   dl_polepi = ((gpinn * gpinn) * fpisq) / ((FT_GeV_neg - (mpig * mpig))*(FT_GeV_neg - (mpig * mpig)));
0136
0137   ratio_Khyp_pi = dl_poleKhyp / dl_polepi;
0138
0139   // --------------------------------------------------------------------------------
0140
0141   fSig_VR = (0.1* fSig_T) + feps * (ratio_Khyp_pi* fSig_L);
0142   sig_total = fSig_VR;
0143
0144   return sig_total;
0145 }