src/eic_evgen/PiPlus_sig.cc

0001 // 21/12/22 - SJDK - The PiPlus cross section routine stripped out of the DEMP_Reaction.cc code in Process routine
0002 // Inputs are -t, W, Qsq and epsilon (in GeV where appropriate)
0003
0004 // Note, there isn't really any reason to split out the parameter determination from this, they could be done in the same file. The PiPlus parameters are a mess anyway.
0005 // In future, these files should also probably be moved elsewhere. The includes/dependencies probably need cleaning up too.
0006
0007 #include "PiPlus_sig.h"
0008 #include "PiPlus_sig_Param.h"
0009 #include "reaction_routine.h"
0010 #include "TF1.h"
0011
0012 double GetPiPlus_CrossSection(double ft, double fw, double fqsq, double feps ){
0013
0014   double_t 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     delete fitCKYLonglandau;
0038     fitCKYLonglandau = NULL;
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     delete fitCKYLongexpo1;
0049     fitCKYLongexpo1 = NULL;
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     delete fitCKYLongexpo2;
0060     fitCKYLongexpo2 = NULL;
0061   }
0062   else {
0063     fSig_L = 0;
0064   }
0065
0066   // -------------------------------------------------------------------------------------------
0067   // 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
0068   // I changed the range to remove this gap.
0069   if ( ( ft > 0.0 ) && ( ft < 0.2 ) ) {
0070     PiPlus_sigmaT_Param( fw,  fqsq, tpar0, tpar1, tpar2 , tpar3 , tpar4 );
0071     TF1 *fitCKYTranspol2 = new TF1("sigmaL","pol2", 0.0 , 0.2 );
0072     fitCKYTranspol2->FixParameter( 0 , tpar0 );
0073     fitCKYTranspol2->FixParameter( 1 , tpar1 );
0074     fitCKYTranspol2->FixParameter( 2 , tpar2 );
0075     fSig_T = fitCKYTranspol2->Eval(ft);
0076     if ( tpar0 == 0 || tpar1 == 0 || tpar2 == 0 )
0077       fSig_T = 0;
0078     delete fitCKYTranspol2;
0079     fitCKYTranspol2 = NULL;
0080   }
0081   else if ( ( ft > 0.2 ) && ( ft < 1.3 ) ) {
0082     PiPlus_sigmaT_Param( fw,  fqsq, tpar0, tpar1, tpar2 , tpar3 , tpar4 );
0083     TF1 *fitCKYTransexpo = new TF1("sigmaL","expo", 0.2 , 1.3 );
0084     fitCKYTransexpo->FixParameter( 0 , tpar3 );
0085     fitCKYTransexpo->FixParameter( 1 , tpar4 );
0086     fSig_T = fitCKYTransexpo->Eval(ft);
0087     if ( tpar3 == 0 || tpar4 == 0 )
0088       fSig_T = 0;
0089     delete fitCKYTransexpo;
0090     fitCKYTransexpo = NULL;
0091   }
0092
0093   // -------------------------------------------------------------------------------------------
0094
0095   fSig_VR = fSig_T + feps * fSig_L;
0096
0097   sig_total = fSig_VR;
0098
0099   return sig_total;
0100 }