eic_evgen/process_routine/DEMP_Reaction.cc

0001 #include "reaction_routine.h"
0002 #include "eic.h"
0003
0004 #include <sys/stat.h>
0005
0006 using namespace std;
0007
0008 DEMP_Reaction::DEMP_Reaction(){
0009
0010   cout << "Program Start" << endl;
0011
0012 }
0013
0014 //-------------------Love Preet - Added for phase space factor calculations----------------------------//
0015
0016 Double_t DEMP_Reaction::calculate_psf_t(const TLorentzVector& psf_photon, const TLorentzVector& psf_ejectile){ // Calculate the Lorentz variable (t) for Phase space factor (psf)/
0017   return -1.0*((psf_photon - psf_ejectile).Mag2());
0018 }
0019
0020 void DEMP_Reaction::calculate_psf_max_min(double value, double& maxValue, double& minValue){// Find the max and min values for scattered electron's energy, theta and ejetile's theta
0021   if (value > maxValue) {
0022     maxValue = value;
0023   }
0024   if (value < minValue) {
0025     minValue = value;
0026   }
0027 }
0028
0029 //-------------------Love Preet - Added for phase space factor (PSF) calculations----------------------------//
0030 Double_t DEMP_Reaction::psf(){
0031
0032   // Added a print out to clarify what is happening
0033   cout << "Beginning phase space calculation -" << endl;
0034
0035   //----Note that all the calculations have been done in GeV for the psf...............................................//
0036   //-----------------------------Calculate variables for scattered electron...........................................//
0037   psf_ScatElec_E_Stepsize = (fEBeam*(fScatElec_E_Hi - fScatElec_E_Lo))/psf_steps; //in GeV // Defined stepsize for scattered electron's energy, theta, and phi
0038   psf_ScatElec_Theta_Stepsize = (fScatElec_Theta_F - fScatElec_Theta_I)/psf_steps; // in radian
0039   psf_ScatElec_Phi_Stepsize = (2.0 * fPi)/4.0;
0040
0041   psf_Ejec_Theta_Stepsize = (fEjectileX_Theta_F - fEjectileX_Theta_I)/psf_steps; // in radian  // Defined stepsize for ejectile theta
0042
0043   for (int psf_E = 0; psf_E <= psf_steps; psf_E++){ // Loop over the scattered electron's energy
0044     // SJDK - 29/04/24 - Added progress report
0045     dFractTime = time(0);
0046
0047     if ( psf_E % ( psf_steps / 10 ) == 0 ){
0048       cout << ((1.0*psf_E)/(1.0*psf_steps))*100.0 << setw(4) << " % of phase space checked"
0049        << "   Day: " <<  dFractTime.GetDay()
0050        << "   Time:   " << dFractTime.GetHour()
0051        << ":" << dFractTime.GetMinute()
0052        << ":" << dFractTime.GetSecond()
0053        << endl;
0054     }
0055
0056     psf_ScatElec_E = (fEBeam * fScatElec_E_Lo + (psf_E * psf_ScatElec_E_Stepsize));
0057     psf_ScalElec_Mom = sqrt(pow(psf_ScatElec_E,2) - pow(fElectron_Mass_GeV,2));
0058
0059     for (int psf_Theta = 0; psf_Theta <= psf_steps; psf_Theta++){  // Loop over the scattered electron's theta
0060
0061       psf_ScatElec_Theta = (fScatElec_Theta_I + (psf_Theta * psf_ScatElec_Theta_Stepsize));
0062
0063       for (int psf_Phi = 0; psf_Phi <= 4; psf_Phi++){ // Loop over the scattered electron's Phi
0064
0065     psf_ScatElec_Phi = (0.0 + (psf_Phi * psf_ScatElec_Phi_Stepsize));
0066     psf_scatelec.SetPxPyPzE(psf_ScalElec_Mom * sin(psf_ScatElec_Theta) * cos( psf_ScatElec_Phi),   // Scattered's electron four momentum
0067                 psf_ScalElec_Mom * sin(psf_ScatElec_Theta) * sin( psf_ScatElec_Phi),
0068                 psf_ScalElec_Mom * cos(psf_ScatElec_Theta),
0069                 psf_ScatElec_E);
0070     psf_photon = r_lelectrong - psf_scatelec; // virtual photon four momentum
0071     psf_Q2 = -1.*(psf_photon.Mag2());         // Lorentz variable Q2
0072     psf_W  =  ((psf_photon + r_lprotong).Mag()); // Lorentz variable W
0073     psf_W2 =  ((psf_photon + r_lprotong).Mag2()); // Lorentz variable W2
0074
0075     if (( psf_Q2 >= fQsq_Min && psf_Q2 <= fQsq_Max) && psf_W2 >= 0.0 && ( psf_W >= fW_Min && psf_W <= fW_Max)){
0076       for (int psf_Ejec_Theta = 0; psf_Ejec_Theta <= psf_steps; psf_Ejec_Theta++){ // Loop over the scattered ejectile's theta
0077
0078         //-----------------------------Calculate variables for ejectile..........................................//
0079         psf_Ejectile_Theta = (fEjectileX_Theta_I + (psf_Ejec_Theta *psf_Ejec_Theta_Stepsize));
0080         psf_Ejectile_Phi = 0.0;
0081
0082         double fupx = sin( psf_Ejectile_Theta ) * cos( psf_Ejectile_Phi );
0083         double fupy = sin( psf_Ejectile_Theta ) * sin( psf_Ejectile_Phi );
0084         double fupz = cos( psf_Ejectile_Theta );
0085
0086         double fuqx = sin( psf_photon.Theta() ) * cos( psf_photon.Phi() );
0087         double fuqy = sin( psf_photon.Theta() ) * sin( psf_photon.Phi() );
0088         double fuqz = cos( psf_photon.Theta() );
0089
0090         double fa = -(psf_photon.Vect()).Mag() * ( fupx * fuqx +  fupy * fuqy +  fupz * fuqz);
0091         double fb = pow ( (psf_photon.Vect()).Mag() , 2 );
0092         double fc = psf_photon.E() +  fProton_Mass_GeV;
0093
0094         fa = ( fa - std::abs( (r_lprotong.Vect()).Mag() ) * ( ( ( r_lprotong.X() / (r_lprotong.Vect()).Mag() ) * fupx ) +
0095                                   ( ( r_lprotong.Y() / (r_lprotong.Vect()).Mag() ) * fupy ) +
0096                                   ( ( r_lprotong.Z() / (r_lprotong.Vect()).Mag() ) * fupz ) ) );
0097
0098         double factor = ( pow( (r_lprotong.Vect()).Mag() , 2 ) + 2.0 * (psf_photon.Vect()).Mag() * (r_lprotong.Vect()).Mag() *
0099                   ( ( ( r_lprotong.X() / (r_lprotong.Vect()).Mag() ) * fuqx ) +
0100                 ( ( r_lprotong.Y() / (r_lprotong.Vect()).Mag() ) * fuqy ) +
0101                 ( ( r_lprotong.Z() / (r_lprotong.Vect()).Mag() ) * fuqz ) ) );
0102
0103         fb =  fb + factor;
0104         fc = psf_photon.E() + r_lprotong.E();
0105
0106         double ft = fc * fc - fb + f_Ejectile_Mass_GeV * f_Ejectile_Mass_GeV -  f_Recoil_Mass_GeV *  f_Recoil_Mass_GeV;
0107
0108         double fQA = 4.0 * ( fa * fa - fc * fc );
0109         double fQB = 4.0 * fc * ft;
0110
0111         double fQC = -4.0 * fa * fa * f_Ejectile_Mass_GeV * f_Ejectile_Mass_GeV - ft * ft;
0112
0113         fradical = fQB * fQB - 4.0 * fQA * fQC;
0114
0115         fepi1 = ( -fQB - sqrt( fradical ) ) / ( 2.0 * fQA );
0116         fepi2 = ( -fQB + sqrt( fradical ) ) / ( 2.0 * fQA );
0117
0118         psf_ejectile.SetPxPyPzE( ( sqrt( pow( fepi1 , 2) - pow(f_Ejectile_Mass_GeV , 2) ) ) * sin( psf_Ejectile_Theta ) * cos(  psf_Ejectile_Phi ),
0119                      ( sqrt( pow( fepi1 , 2) - pow(f_Ejectile_Mass_GeV , 2) ) ) * sin( psf_Ejectile_Theta ) * sin(  psf_Ejectile_Phi ),
0120                      ( sqrt( pow( fepi1 , 2) - pow(f_Ejectile_Mass_GeV , 2) ) ) * cos( psf_Ejectile_Theta ),
0121                      fepi1 );
0122
0123         if (!TMath::IsNaN(psf_ejectile.E())) { // If the ejectile's energy is a number
0124
0125           psf_t = calculate_psf_t(psf_photon, psf_ejectile);
0126
0127           if (psf_t >= 0.0 && psf_t <= fT_Max) { // If the t value is in this range
0128
0129         calculate_psf_max_min(psf_ScatElec_E, psf_ScatElec_E_max, psf_ScatElec_E_min);
0130         calculate_psf_max_min(psf_ScatElec_Theta, psf_ScatElec_Theta_max, psf_ScatElec_Theta_min);
0131         calculate_psf_max_min(psf_Ejectile_Theta, psf_Ejectile_Theta_max, psf_Ejectile_Theta_min);
0132
0133           } // if statement over psf_t
0134         } // if statement over ejectile energy
0135       } // End of for loop over psf_Theta_Ejec
0136     } // If condition over psf_Q2,psf_W,psf_W2
0137       } // End of for loop over psf_Phi
0138     } // End of for loop over psf_Theta
0139   } // End of for loop over psf_Phi
0140
0141   if ((psf_ScatElec_E_max == psf_ScatElec_Theta_max) && (psf_ScatElec_Theta_max == psf_Ejectile_Theta_max) && (psf_Ejectile_Theta_max == psf_ScatElec_E_max)){
0142     cout << "!!!!! - ERROR - !!!!!" << endl;
0143     cout << "Specified scattered electron and ejectile energy/angle ranges in input .json file did not yield any valid allowed phase space." << endl << "Processing stopped, please re-specify values in input .json file." << endl;
0144     cout << "!!!!! - ERROR - !!!!!" << endl;
0145     fPSF = 0.0;
0146   }
0147   else if ((psf_ScatElec_E_min == psf_ScatElec_Theta_min) && (psf_ScatElec_Theta_min == psf_Ejectile_Theta_min) && (psf_Ejectile_Theta_min == psf_ScatElec_E_min)){
0148     cout << "!!!!! - ERROR - !!!!!" << endl;
0149     cout << "Specified scattered electron and ejectile energy/angle ranges in input .json file did not yield any valid allowed phase space." << endl << " Processing stopped, please re-specify values in input .json file." << endl;
0150     cout << "!!!!! - ERROR - !!!!!" << endl;
0151     fPSF = 0.0;
0152   }
0153   else if ((psf_ScatElec_E_max == psf_ScatElec_E_min) || (psf_ScatElec_Theta_max == psf_ScatElec_Theta_min) || (psf_Ejectile_Theta_max == psf_Ejectile_Theta_min)) {
0154     cout << "!!!!! - ERROR - !!!!!" << endl;
0155     cout << "Specified scattered electron and ejectile energy/angle ranges in input .json file include one (or more) values with equal min/max values in range." << endl << "Processing stopped, please re-specify values in input .json file." << endl;
0156     cout << "!!!!! - ERROR - !!!!!" << endl;
0157     fPSF = 0.0;
0158   }
0159   else {
0160     // Calculate the PSF by expanding the limits on both sides by a step size of 1.0.
0161     psf_ScatElec_E_min      = psf_ScatElec_E_min     - ( 1.0 * psf_ScatElec_E_Stepsize ); // Units are in GeV
0162     psf_ScatElec_E_max      = psf_ScatElec_E_max     + ( 1.0 * psf_ScatElec_E_Stepsize );
0163     psf_ScatElec_Theta_min  = psf_ScatElec_Theta_min - ( 1.0 * psf_ScatElec_Theta_Stepsize ); // Units are in Radian
0164     psf_ScatElec_Theta_max  = psf_ScatElec_Theta_max + ( 1.0 * psf_ScatElec_Theta_Stepsize );
0165     psf_Ejectile_Theta_min  = psf_Ejectile_Theta_min - ( 1.0 * psf_Ejec_Theta_Stepsize ); // Units are in Radian
0166     psf_Ejectile_Theta_max  = psf_Ejectile_Theta_max + ( 1.0 * psf_Ejec_Theta_Stepsize );
0167
0168     if ( psf_ScatElec_E_min < ( fEBeam * fScatElec_E_Lo ) || psf_ScatElec_E_min == ( fEBeam * fScatElec_E_Lo ) ) { psf_ScatElec_E_min = ( fEBeam * fScatElec_E_Lo ); }
0169     if ( psf_ScatElec_E_max > ( fEBeam * fScatElec_E_Hi ) || psf_ScatElec_E_max == ( fEBeam * fScatElec_E_Hi ) ) { psf_ScatElec_E_max = ( fEBeam * fScatElec_E_Hi ); }
0170     if ( psf_ScatElec_Theta_min < fScatElec_Theta_I || psf_ScatElec_Theta_min == fScatElec_Theta_I ) { psf_ScatElec_Theta_min = fScatElec_Theta_I; }
0171     if ( psf_ScatElec_Theta_max > fScatElec_Theta_F || psf_ScatElec_Theta_max == fScatElec_Theta_F ) { psf_ScatElec_Theta_max = fScatElec_Theta_F; }
0172     if ( psf_Ejectile_Theta_min < fEjectileX_Theta_I || psf_Ejectile_Theta_min == fEjectileX_Theta_I ) { psf_Ejectile_Theta_min = fEjectileX_Theta_I; }
0173     if ( psf_Ejectile_Theta_max > fEjectileX_Theta_F || psf_Ejectile_Theta_max == fEjectileX_Theta_F ) { psf_Ejectile_Theta_max = fEjectileX_Theta_F; }
0174
0175     fPSF = (( psf_ScatElec_E_max - psf_ScatElec_E_min ) *( cos( psf_ScatElec_Theta_max ) - cos( psf_ScatElec_Theta_min ) ) * 2 * fPI *( cos( psf_Ejectile_Theta_max ) - cos( psf_Ejectile_Theta_min ) ) * 2 * fPI );
0176   }
0177   return fPSF;
0178 }
0179
0180 /*--------------------------------------------------*/
0181 /// DEMP_Reaction
0182 DEMP_Reaction::DEMP_Reaction(TString particle_str, TString hadron_str){
0183
0184   rEjectile = particle_str;
0185   rRecoil = hadron_str;
0186
0187 }
0188
0189 DEMP_Reaction::~DEMP_Reaction(){
0190
0191   // Text output data file and generation information
0192   DEMPOut.close();
0193   DEMPDetails.close();
0194
0195   if (gROOTOut == true){
0196     // Diagnostic root file with plots
0197     dRootFile->Write(); // Write the contents of the ROOT file to disk.
0198     delete dRootTree;   // Delete the dynamically allocated memory for the ROOT tree.
0199     dRootFile->Close(); // Close the ROOT file.
0200     delete dRootFile;   // Delete the dynamically allocated memory for the ROOT file object.
0201   }
0202 }
0203
0204 void DEMP_Reaction::process_reaction(){
0205
0206   Init();
0207
0208   if (fPSF <= 0){// If phase space factor is zero or less than zero, stop further processing // Love Preet - Added for phase space factor calculations
0209     return;
0210   }
0211
0212   //-------------------Love Preet - Added to print out the considered ranges for PSF calculations----------------------------//
0213
0214   cout << "------------------------------------------------------------------------------------------------------------------------------------------------------" << endl;
0215   cout << "Specified scattered electron and ejectile energy/angle ranges in input .json file are ajusted to ensure that they fall within the allowed phase space." << endl;
0216   cout << "Ee_Low = " << psf_ScatElec_E_min << ", Ee_High = " << psf_ScatElec_E_max << ", e_Theta_Low = " << psf_ScatElec_Theta_min * TMath::RadToDeg() << ", e_Theta_High = " << psf_ScatElec_Theta_max * TMath::RadToDeg() << ", EjectileX_Theta_Low = " << psf_Ejectile_Theta_min * TMath::RadToDeg() << ", EjectileX_Theta_High = " << psf_Ejectile_Theta_max * TMath::RadToDeg() << endl;
0217   cout << "------------------------------------------------------------------------------------------------------------------------------------------------------" << endl;
0218
0219   //-------------------Love Preet - Added to print out the considered ranges for PSF calculations----------------------------//
0220
0221   if (gOutputType == "Pythia6"){
0222     DEMPReact_Pythia6_Out_Init();
0223   }
0224   else if (gOutputType == "HEPMC3"){
0225     DEMPReact_HEPMC3_Out_Init();
0226   }
0227
0228   for( long long int i = 0; i < rNEvents; i++ ){
0229
0230     rNEvent_itt = i;
0231     fNGenerated ++;
0232
0233     Progress_Report();  // This happens at each 10% of the total event is processed
0234     Processing_Event();
0235   }
0236
0237   //-------------------Love Preet - Added for actual phase space factor calculations----------------------------//
0238   fPSF_org = (( fScatElec_Energy_Col_max * fm - fScatElec_Energy_Col_min * fm ) *( cos( fScatElec_Theta_Col_max ) - cos( fScatElec_Theta_Col_min ) ) * 2 * fPI *( cos( f_Ejectile_Theta_Col_max ) - cos( f_Ejectile_Theta_Col_min ) ) * 2 * fPI ); // Calculate the actual phase space factor, this is determined from the actual scattered electron/ejectile values in the generated events that passed all cuts
0239
0240   //-------------------Love Preet - Added for actual phase space factor calculations----------------------------//
0241
0242   Detail_Output();
0243
0244 }
0245
0246 void DEMP_Reaction::Init(){
0247
0248   pim* myPim;
0249
0250   pd = dynamic_cast<pim*>(myPim);
0251
0252   rEjectile_charge = ExtractCharge(rEjectile);
0253
0254   const char* dir_name = "OutputFiles";
0255   struct stat sb;
0256
0257   if (stat(dir_name, &sb) == 0) {
0258     cout << "Output file directory found from DEMPgen directory - " << dir_name  << endl;
0259   }
0260   else {
0261     cout << "Output file directory not found from DEMPgen directory - " << dir_name  << endl;
0262     cout << "Making OutputFiles directory!" << endl;
0263     mkdir(dir_name,0777);
0264   }
0265
0266   sTFile = Form("./%s/eic_%s.txt", dir_name, gfile_name.Data());
0267   sLFile = Form("./%s/eic_input_%s.dat", dir_name, gfile_name.Data());
0268
0269   DEMPOut.open( sLFile.c_str() );
0270   DEMPDetails.open( sTFile.c_str() );
0271
0272   if (gROOTOut == true){ // Only initialise and open root file if output is enabled
0273     sDFile = Form("./%s/eic_%s.root", dir_name, gfile_name.Data()); // LovePreet changed to make the files name consistent
0274     dRootFile = new TFile(sDFile.c_str(),"RECREATE");
0275     dRootTree = new TTree("Events", "Description of a tree");  //Love Preet added all these new braches to be stored in a root ttree
0276     dRootTree->Branch("EventWeight", &fEventWeight, "EventWeight/D");
0277     dRootTree->Branch("scat_e_px", &scat_e_px, "scat_e_px/D");
0278     dRootTree->Branch("scat_e_py", &scat_e_py, "scat_e_py/D");
0279     dRootTree->Branch("scat_e_pz", &scat_e_pz, "scat_e_pz/D");
0280     dRootTree->Branch("scat_e_E",  &scat_e_E, "scat_e_E/D");
0281     dRootTree->Branch("ejec_px", &ejec_px, "ejec_px/D");
0282     dRootTree->Branch("ejec_py", &ejec_py, "ejec_py/D");
0283     dRootTree->Branch("ejec_pz", &ejec_pz, "ejec_pz/D");
0284     dRootTree->Branch("ejec_E",  &ejec_E, "ejec_E/D");
0285     dRootTree->Branch("rclH_px", &rclH_px, "rclH_px/D");
0286     dRootTree->Branch("rclH_py", &rclH_py, "rclH_py/D");
0287     dRootTree->Branch("rclH_pz", &rclH_pz, "rclH_pz/D");
0288     dRootTree->Branch("rclH_E",  &rclH_E, "rclH_E/D");
0289     dRootTree->Branch("Q2", &fQsq_GeV, "Q2/D");
0290     dRootTree->Branch("W",  &fW_GeV, "W/D");
0291     dRootTree->Branch("t", &fT_GeV, "t/D");
0292     dRootTree->Branch("x_b", &fx, "x_b/D");
0293     dRootTree->Branch("y_E", &fy, "y_E/D");
0294   }
0295   /*--------------------------------------------------*/
0296   qsq_ev = 0, t_ev = 0, w_neg_ev = 0, w_ev = 0;
0297   rNEvents = fNEvents;
0298   rNEvent_itt = 0;
0299
0300   // 02/06/21 SJDK
0301   // Set these values once the beam energies are read in
0302   fElectron_Kin_Col_GeV = fEBeam;
0303   fElectron_Kin_Col = fElectron_Kin_Col_GeV * 1000.0;
0304
0305   // cout << rNEvents << "    " << fNEvents << endl;
0306
0307   // 08/09/23 - SJDK - Fermi momentum commented out for now, this is not fully implemented yet
0308   // In future, this will be enabled/disabled automatically depending upon the specified hadron beam
0309   //rFermiMomentum = pd->fermiMomentum();
0310
0311   // ----------------------------------------------------
0312   // Proton in collider (lab) frame
0313   // 08/09/23 - SJDK - The naming needs to be adjusted to be the incoming hadron beam, not the proton. Make this generic.
0314   r_lproton = GetProtonVector_lab();
0315   r_lprotong = GetProtonVector_lab() * fm;
0316
0317   // Getting the mass of the hadron beam
0318   r_lhadron_beam_mass = ParticleMass(ParticleEnum(gBeamPart.c_str()))*1000; // in MeV
0319
0320   //  cout << gBeamPart << endl;
0321   //  cout << ParticleEnum(gBeamPart.c_str()) << endl;
0322   //  cout << ParticleMass(ParticleEnum(gBeamPart.c_str())) << endl;
0323   //  cout << r_lhadron_beam_mass << endl;
0324   //  exit(0);
0325
0326   // ----------------------------------------------------
0327   // Electron in collider (lab) frame
0328
0329   // 06/09/23 - SJDK - Commenting out for now, should be disabled for e/p collisions
0330   //cout << "Fermi momentum: " << rFermiMomentum << endl;
0331
0332   r_lelectron    = GetElectronVector_lab();
0333   r_lelectrong = r_lelectron * fm;
0334   cout << "Define: " << fElectron_MomZ_Col << "    "<< fElectron_Mom_Col << "  " << cos(fElectron_Theta_Col) << endl;
0335
0336   ///*--------------------------------------------------*/
0337   /// Getting the ejectile (produced meson) particle mass from the data base
0338
0339   produced_X = ParticleEnum(rEjectile);
0340   f_Ejectile_Mass = ParticleMass(produced_X)*1000; //MeV
0341   f_Ejectile_Mass_GeV = f_Ejectile_Mass/1000; //GeV
0342
0343   cout << rEjectile << "  " << produced_X << "  " << f_Ejectile_Mass_GeV <<  endl;
0344   cout << rEjectile_charge << endl;
0345
0346   if (rRecoil == "Neutron" ) {
0347     rEjectile_scat_hadron  = "Neutron";
0348     recoil_hadron  = Neutron;
0349     f_Recoil_Mass     = fNeutron_Mass;
0350     f_Recoil_Mass_GeV = f_Recoil_Mass/1000;
0351   }
0352   else if (rRecoil == "Proton" ) {
0353     rEjectile_scat_hadron  = "Proton";
0354     recoil_hadron  = Proton;
0355     f_Recoil_Mass = fProton_Mass;
0356     f_Recoil_Mass_GeV = f_Recoil_Mass/1000;
0357   }
0358   else if(rRecoil == "Lambda"){
0359     rEjectile_scat_hadron  = "Lambda";
0360     recoil_hadron  = Lambda;
0361     f_Recoil_Mass = fLambda_Mass;
0362     f_Recoil_Mass_GeV = f_Recoil_Mass/1000;
0363     cout<<"Particle = "<<rEjectile_scat_hadron<<" with mass = "<< f_Recoil_Mass << endl;
0364   }
0365   else if (rRecoil == "Sigma0"){
0366     rEjectile_scat_hadron  = "Sigma0";
0367     recoil_hadron  = Sigma0;
0368     f_Recoil_Mass  = fSigma_Mass;
0369     f_Recoil_Mass_GeV = f_Recoil_Mass/1000;
0370     cout<<"Particle = "<<rEjectile_scat_hadron<<" with mass = "<< f_Recoil_Mass << endl;
0371   }
0372
0373   rDEG2RAD   = fPI/180.0;
0374
0375   f_Ejectile_Theta_I = fEjectileX_Theta_I ;
0376   f_Ejectile_Theta_F = fEjectileX_Theta_F;
0377
0378   fPSF = psf(); // Love Preet - Added for phase space factor calculations
0379   if (fPSF <= 0){ // if phase space factor is zero or less than zero, stop further processing // Love Preet - Added for phase space factor calculations
0380     return;
0381   }
0382
0383   cout << "Produced particle in exclusive production: " << rEjectile << ";  with mass: " << f_Ejectile_Mass << " MeV "<< endl;
0384   cout << fEBeam << " GeV electrons on " << fHBeam << " GeV ions" << endl;
0385   if(UseSolve == true){
0386     cout << rEjectile << " and " << rEjectile_scat_hadron << " 4-vectors calculated using Solve function" << endl;
0387   }
0388   else if(UseSolve == false){
0389     cout << rEjectile << " and " << rEjectile_scat_hadron << " 4-vectors calculated using analytical solution" << endl;
0390   }
0391   // Set luminosity value based upon beam energy combination, note that if no case matches, a default of 1e33 is assumed. Cases are a set of nominal planned beam energy combinations for the EIC (and EICC)
0392   // See slide 11 in https://indico.cern.ch/event/1072579/contributions/4796856/attachments/2456676/4210776/CAP-EIC-June-7-2022-Seryi-r2.pdf
0393   // If available in the future, this could be replaced by some fixed function
0394   if ((fEBeam == 5.0 ) && (fHBeam == 41.0) ){
0395     fLumi = 0.44e33;
0396   }
0397   else if ((fEBeam == 5.0 ) && (fHBeam == 100.0) ){
0398     fLumi = 3.68e33;
0399   }
0400   else if ((fEBeam == 10.0 ) && (fHBeam == 100.0) ){
0401     fLumi = 4.48e33;
0402   }
0403   else if ((fEBeam == 18.0 ) && (fHBeam == 275.0) ){
0404     fLumi = 1.54e33;
0405   }
0406   else if ((fEBeam == 3.5 ) && (fHBeam == 20) ){ // EICC optimal beam energy combination
0407     fLumi = 2e33;
0408   }
0409   else if ((fEBeam == 2.8 ) && (fHBeam == 13) ){ // EICC lowest beam energy combination
0410     fLumi = 0.7e33;
0411   }
0412   else{
0413     cout << "!!! Notice !!! The beam energy combination simulated does not match an expected case, a default luminosity value of - " << fLumi << " cm^2s^-1 has been assumed. !!! Notice !!!" << endl;
0414   }
0415
0416   if(UseSolve == true){
0417     /*--------------------------------------------------*/
0418     // SJDK 03/04/22 -  New set of initialisation stuff for the solve function from Ishan and Bill
0419
0420     CoinToss = new TRandom3();
0421
0422     F = new TF1("F",
0423         "[6]-sqrt([7]**2+x**2)-sqrt([8]**2+([3]-[0]*x)**2+([4]-[1]*x)**2+([5]-[2]*x)**2)",
0424         0, r_lproton.E());
0425
0426     char AngleGenName[100] = "AngleGen";
0427     double dummy[2] = {0,1};
0428
0429     f_Ejectile_Theta_I = fEjectileX_Theta_I;
0430     f_Ejectile_Theta_F = fEjectileX_Theta_F;
0431
0432     double ThetaRange[2] = {f_Ejectile_Theta_I, f_Ejectile_Theta_F};
0433     double PhiRange[2] = {0, 360*TMath::DegToRad()};
0434
0435     AngleGen = new CustomRand(AngleGenName, dummy,
0436                   ThetaRange, PhiRange);
0437
0438     UnitVect = new TVector3(0,0,1);
0439
0440     //  ///*--------------------------------------------------*/
0441     //  // Produced hadron and recoilded hadron from the solve function
0442
0443     VertBeamElec = new TLorentzVector();
0444     VertScatElec = new TLorentzVector();
0445
0446     Initial      = new TLorentzVector();
0447     Target       = new TLorentzVector();
0448     Photon       = new TLorentzVector();
0449     Interaction  = new TLorentzVector();
0450     Final        = new TLorentzVector();
0451   }
0452 }
0453
0454 void DEMP_Reaction::Processing_Event(){
0455
0456   // ----------------------------------------------------
0457   // Considering Fermi momentum for the proton
0458   // ----------------------------------------------------
0459   // SJDK - 06/06/23 - Commenting out for now, this increases the number of recorded events - Should have no fermi momentum for e/p collisions
0460   // if( kCalcFermi ) {
0461   //   Consider_Proton_Fermi_Momentum();
0462   // }
0463   // ----------------------------------------------------
0464   // Boost vector from collider (lab) frame to protons rest frame (Fix target)
0465   // ----------------------------------------------------
0466   beta_col_rf = r_lproton.BoostVector();
0467   fGamma_Col_RF = 1.0/sqrt( 1 - pow( beta_col_rf.Mag() , 2 ) );
0468
0469   // ---------------------------------------------------------------------
0470   // Specify the energy and solid angle of scatterd electron in Collider (lab) frame
0471   // ---------------------------------------------------------------------
0472   //fScatElec_Theta_Col  = acos( fRandom->Uniform( cos( fScatElec_Theta_I ) , cos( fScatElec_Theta_F ) ) );
0473   fScatElec_Theta_Col  = acos( fRandom->Uniform( cos( psf_ScatElec_Theta_min ) , cos( psf_ScatElec_Theta_max ) ) ); // Love Preet changed to the allowed region calculated from psf
0474   fScatElec_Phi_Col    = fRandom->Uniform( 0 , 2.0 * fPi);
0475   //fScatElec_Energy_Col = fRandom->Uniform( fScatElec_E_Lo * fElectron_Energy_Col , fScatElec_E_Hi * fElectron_Energy_Col ); // in MeV
0476   fScatElec_Energy_Col = fRandom->Uniform( psf_ScatElec_E_min * fK, psf_ScatElec_E_max * fK ); // in MeV // Love Preet changed to the allowed region calculated from psf
0477
0478   // ----------------------------------------------------
0479   // Produced ejectile in Collider frame
0480   // ----------------------------------------------------
0481   //f_Ejectile_Theta_Col      = acos( fRandom->Uniform( cos(f_Ejectile_Theta_I), cos(f_Ejectile_Theta_F ) ) );
0482   f_Ejectile_Theta_Col      = acos( fRandom->Uniform( cos( psf_Ejectile_Theta_min ), cos( psf_Ejectile_Theta_max ) ) ); // Love Preet changed to the allowed region calculated from psf
0483   f_Ejectile_Phi_Col        = fRandom->Uniform( 0 , 2.0 * fPi );
0484
0485   //---------------------------------------------------------------------
0486   // Specify the energy and solid angle of scatterd electron in Collider (lab) frame
0487   // ---------------------------------------------------------------------
0488
0489   fScatElec_Mom_Col  = sqrt( pow( fScatElec_Energy_Col,2) - pow( fElectron_Mass , 2) );
0490   fScatElec_MomZ_Col = ( fScatElec_Mom_Col * cos(fScatElec_Theta_Col) );
0491   fScatElec_MomX_Col = ( fScatElec_Mom_Col * sin(fScatElec_Theta_Col) * cos(fScatElec_Phi_Col) );
0492   fScatElec_MomY_Col = ( fScatElec_Mom_Col * sin(fScatElec_Theta_Col) * sin(fScatElec_Phi_Col) );
0493
0494   r_lscatelec.SetPxPyPzE( fScatElec_MomX_Col, fScatElec_MomY_Col, fScatElec_MomZ_Col, fScatElec_Energy_Col);
0495
0496   r_lscatelecg = r_lscatelec * fm;
0497
0498   // ----------------------------------------------------
0499   // Photon in collider (lab) frame and Qsq
0500   // ----------------------------------------------------
0501
0502   r_lphoton  = r_lelectron - r_lscatelec;
0503   r_lphotong = r_lelectrong - r_lscatelecg;
0504
0505   fQsq_GeV = -1.* r_lphotong.Mag2();
0506
0507   // ----------------------------------------------------
0508   // W square, Invariant Mass (P_g + P_p)^2
0509   // ----------------------------------------------------
0510
0511   TLorentzVector lwg;
0512   lwg = r_lprotong + r_lphotong;
0513   fW_GeV    = lwg.Mag();
0514   fWSq_GeV  = lwg.Mag2();
0515
0516   // SJDK - 06/09/23 - Check UseSolved boolean, process through relevant loop
0517   if (UseSolve == false){
0518     // ---------------------------------------------------------
0519     // Pion momentum in collider frame, analytic solution starts
0520     // ---------------------------------------------------------
0521
0522     double fupx = sin( f_Ejectile_Theta_Col ) * cos( f_Ejectile_Phi_Col );
0523     double fupy = sin( f_Ejectile_Theta_Col ) * sin( f_Ejectile_Phi_Col );
0524     double fupz = cos( f_Ejectile_Theta_Col );
0525
0526     double fuqx = sin( r_lphoton.Theta() ) * cos( r_lphoton.Phi() );
0527     double fuqy = sin( r_lphoton.Theta() ) * sin( r_lphoton.Phi() );
0528     double fuqz = cos( r_lphoton.Theta() );
0529
0530     double fa = -(r_lphoton.Vect()).Mag() * ( fupx * fuqx +  fupy * fuqy +  fupz * fuqz );
0531     double fb = pow ( (r_lphoton.Vect()).Mag() , 2 );
0532     double fc = r_lphoton.E() + fProton_Mass;
0533
0534     fa = ( fa - std::abs( (r_lproton.Vect()).Mag() ) * ( ( ( r_lproton.X() / (r_lproton.Vect()).Mag() ) * fupx ) +
0535                              ( ( r_lproton.Y() / (r_lproton.Vect()).Mag() ) * fupy ) +
0536                              ( ( r_lproton.Z() / (r_lproton.Vect()).Mag() ) * fupz ) ) );
0537
0538     double factor = ( pow( (r_lproton.Vect()).Mag() , 2 ) + 2.0 * (r_lphoton.Vect()).Mag() * (r_lproton.Vect()).Mag() *
0539               ( ( ( r_lproton.X() / (r_lproton.Vect()).Mag() ) * fuqx ) +
0540             ( ( r_lproton.Y() / (r_lproton.Vect()).Mag() ) * fuqy ) +
0541             ( ( r_lproton.Z() / (r_lproton.Vect()).Mag() ) * fuqz ) ) );
0542
0543     fb =  fb + factor;
0544     fc = r_lphoton.E() + r_lproton.E();
0545
0546     double ft = fc * fc - fb + f_Ejectile_Mass * f_Ejectile_Mass - f_Recoil_Mass * f_Recoil_Mass;
0547
0548     double fQA = 4.0 * ( fa * fa - fc * fc );
0549     double fQB = 4.0 * fc * ft;
0550
0551     double fQC = -4.0 * fa * fa * f_Ejectile_Mass * f_Ejectile_Mass - ft * ft;
0552
0553     fradical = fQB * fQB - 4.0 * fQA * fQC;
0554
0555     fepi1 = ( -fQB - sqrt( fradical ) ) / ( 2.0 * fQA );
0556     fepi2 = ( -fQB + sqrt( fradical ) ) / ( 2.0 * fQA );
0557
0558     ///---------------------------------------------------------
0559     /// Particle X momentum in collider frame, analytic solution
0560     /// And obtain recoiled proton in collider (lab) frame
0561     ///---------------------------------------------------------
0562
0563     r_l_Ejectile.SetPxPyPzE( (sqrt( pow( fepi1 , 2) - pow(f_Ejectile_Mass , 2) ) ) * sin(f_Ejectile_Theta_Col) * cos(f_Ejectile_Phi_Col),
0564                  ( sqrt( pow( fepi1 , 2) - pow(f_Ejectile_Mass , 2) ) ) * sin(f_Ejectile_Theta_Col) * sin(f_Ejectile_Phi_Col),
0565                  ( sqrt( pow( fepi1 , 2) - pow(f_Ejectile_Mass , 2) ) ) * cos(f_Ejectile_Theta_Col),
0566                  fepi1 );
0567
0568     l_Recoil.SetPxPyPzE( ( r_lproton + r_lelectron - r_lscatelec - r_l_Ejectile).X(),
0569              ( r_lproton + r_lelectron - r_lscatelec - r_l_Ejectile ).Y(),
0570              ( r_lproton + r_lelectron - r_lscatelec - r_l_Ejectile ).Z(),
0571              sqrt( pow( ( ( ( r_lproton + r_lelectron - r_lscatelec - r_l_Ejectile ).Vect() ).Mag()),2) +
0572                    pow( f_Recoil_Mass , 2) ) );
0573   }
0574   else if (UseSolve == true){
0575     if(!Solve()){
0576       return;
0577     }
0578     r_l_Ejectile = r_l_Ejectile_solved;
0579     l_Recoil = r_l_Recoil_solved;
0580   }
0581
0582   ///--------------------------------------------------
0583
0584   r_l_Ejectile_g = r_l_Ejectile * fm;
0585   l_Recoil_g = l_Recoil * fm;
0586
0587   // ----------------------------------------------------------------------------------------------
0588   // Calculate w = (proton + photon)^2
0589   // ----------------------------------------------------------------------------------------------
0590
0591   r_lw = r_lproton + r_lphoton;
0592   fW = r_lw.Mag();
0593
0594   // SJDK 15/06/21 - Added integer counters for conservation law check and for NaN check
0595   if (r_l_Ejectile.E() != r_l_Ejectile.E()){ // SJDK 15/06/21 - If the energy of the produced meson is not a number, return and add to counter
0596     fNaN++;
0597     return;
0598   }
0599   //*--------------------------------------------------*/
0600   //-> 10/05/23 - Love added a slimmed down, simpler to read version of the CheckLaws fn
0601   //
0602   // To check the conservation of the energy and momentum, there two methods avalaible:
0603   // Method 1: Give the four-vectors of the initial and final states partciles,
0604   //           tolerance factor will be defaulted 1e-6 MeV
0605   //           CheckLaws(e_beam, h_beam, scatt_e, ejectile, recoil) <- input 4 vectors
0606   // Method 2: Give the four-vectors of the initial and final states partciles,
0607   //           and the prefered tolerance factor.
0608   //           CheckLaws(e_beam, h_beam, scatt_e, ejectile, recoil, tolerance) <- input 4 vectors and tolerance value in GeV
0609   // Both functions return 1 if conservations laws are satisified
0610
0611   if( pd->CheckLaws(r_lelectron, r_lproton, r_lscatelec, r_l_Ejectile, l_Recoil) !=1 ){ // Love Preet changed the order of the conditions so that can identify unphysical events based on NaN ejectile energy, conservation law, and negative W.
0612     fConserve++;
0613     return;
0614   }
0615
0616   if ( fWSq_GeV < 0 ) {
0617     w_neg_ev++;
0618     return;
0619   }
0620
0621   // SJDK 03/04/23 - Qsq an W ranges now variables set by particle type in .json read in. See eic.cc
0622   if ( fQsq_GeV < fQsq_Min || fQsq_GeV > fQsq_Max ) {
0623     qsq_ev++;
0624     return;
0625   }
0626
0627   if ( fW_GeV < fW_Min || fW_GeV > fW_Max ) { // SJDK 03/04/23 - Switched to the new variable, set by particle type
0628     w_ev++;
0629     return;
0630   }
0631
0632   ////////////////////////////////////////////////////////////////////////////////////////////
0633   //                                          Start                                         //
0634   // Transformation of e', pi- and recoil proton to target's rest frame without energy loss //
0635   ////////////////////////////////////////////////////////////////////////////////////////////
0636
0637   lproton_rf = r_lproton;
0638   lproton_rf.Boost(-beta_col_rf);
0639   lproton_rfg = lproton_rf * fm;
0640
0641   lelectron_rf = r_lelectron;
0642   lelectron_rf.Boost(-beta_col_rf);
0643   lelectron_rfg = lelectron_rf * fm;
0644
0645   lscatelec_rf = r_lscatelec;
0646   lscatelec_rf.Boost(-beta_col_rf);
0647   lscatelec_rfg = lscatelec_rf * fm;
0648
0649   lphoton_rf = r_lphoton;
0650   lphoton_rf.Boost(-beta_col_rf);
0651   lphoton_rfg = lphoton_rf * fm;
0652
0653   l_Ejectile_rf = r_l_Ejectile;
0654   l_Ejectile_rf.Boost(-beta_col_rf);
0655   l_Ejectile_rfg = l_Ejectile_rf * fm;
0656
0657   l_scat_hadron_rf = l_Recoil;
0658   l_scat_hadron_rf.Boost(-beta_col_rf);
0659   l_scat_hadron_rf_g = l_scat_hadron_rf * fm;
0660
0661   ////////////////////////////////////////////////////////////////////////////////////////////
0662   //                                          End                                           //
0663   // Transformation of e', pi- and recoil proton to target's rest frmae without energy loss //
0664   ////////////////////////////////////////////////////////////////////////////////////////////
0665
0666   // -----------------------------------------------------------------------------------------
0667   // Calculate -t
0668   // -----------------------------------------------------------------------------------------
0669   // 18/05/23 - SJDK - Should these be proton mass still or not?
0670
0671   //fBeta_CM_RF        = (lphoton_rf.Vect()).Mag() / (lphoton_rf.E() + r_lhadron_beam_mass);
0672   fBeta_CM_RF        = ( lphoton_rfg.Vect() ).Mag() / ( lphoton_rfg.E() + ( r_lhadron_beam_mass * fm ) ); // involved quantities are in GeV
0673
0674   //fGamma_CM_RF       = (lphoton_rf.E() + r_lhadron_beam_mass) / fW;
0675   fGamma_CM_RF       = ( lphoton_rfg.E() + ( r_lhadron_beam_mass * fm ) ) / ( fW * fm ); // involved quantities are in GeV
0676   f_Ejectile_Energy_CM       = (pow(fW, 2) + pow(f_Ejectile_Mass,2) - pow(f_Recoil_Mass,2) ) / (2.0* fW);
0677   f_Ejectile_Mom_CM          = sqrt(pow(f_Ejectile_Energy_CM,2) - pow(f_Ejectile_Mass,2));
0678   //f_Ejectile_Energy_CM_GeV   = f_Ejectile_Energy_CM / 1000.0;
0679   //f_Ejectile_Mom_CM_GeV      = f_Ejectile_Mom_CM / 1000.0;
0680   f_Ejectile_Energy_CM_GeV   = f_Ejectile_Energy_CM * fm;
0681   f_Ejectile_Mom_CM_GeV      = f_Ejectile_Mom_CM * fm;
0682
0683   // this equation is valid for parallel kinematics only!
0684   fT_Para = ( pow(((r_lphoton.Vect()).Mag() - (r_l_Ejectile.Vect()).Mag()),2) - pow((r_lphoton.E() - r_l_Ejectile.E()),2));
0685   fT_Para_GeV = fT_Para/1000000.0;
0686
0687   lt = r_lphoton - r_l_Ejectile;
0688   ltg = lt * fm;
0689
0690   fT = -1.*lt.Mag2();
0691   fT_GeV = -1.*ltg.Mag2();
0692
0693   // 31/01/23 - SJDK - Kinematics type 1 is FF and 2 is TSSA, this reaction class shouldn't care about this and only have a limit depending upon particle type?
0694   /*
0695     if ( gKinematics_type == 1 && fT_GeV > 0.5 ) {
0696     t_ev++;
0697     return;
0698     }
0699
0700     if ( gKinematics_type == 2 && fT_GeV > 1.3 ) {
0701     t_ev++;
0702     return;
0703     }
0704   */
0705
0706   // 31/01/23 SJDK - New limit on t, remove only events outside the parameterisation range
0707   // 06/09/23 SJDK - fT_Max set in eic.cc depending upon ejectile type
0708   if (fT_GeV > fT_Max ) {
0709     t_ev++;
0710     return;
0711   }
0712
0713   fx = fQsq_GeV / ( 2.0 * r_lprotong.Dot( r_lphotong ) );
0714   fy = r_lprotong.Dot( r_lphotong ) / r_lprotong.Dot( r_lelectrong );
0715   fz = r_l_Ejectile.E()/r_lphoton.E();
0716
0717   // -------------------------------------------------------------------------------------------------------
0718   // Calculation of Phi  ( azimuthal angle of pion momentum w.r.t lepton plane in target's rest frame)
0719   // Calculation of PhiS ( azimuthal angle of target polarization w.r.t lepton plane in target's rest frame)
0720   // -------------------------------------------------------------------------------------------------------
0721
0722   v3Photon.SetX( lphoton_rfg.X() );
0723   v3Photon.SetY( lphoton_rfg.Y() );
0724   v3Photon.SetZ( lphoton_rfg.Z() );
0725
0726   v3Electron.SetX( lelectron_rfg.X() );
0727   v3Electron.SetY( lelectron_rfg.Y() );
0728   v3Electron.SetZ( lelectron_rfg.Z() );
0729
0730   v3X.SetX( l_Ejectile_rfg.X() ) ;
0731   v3X.SetY( l_Ejectile_rfg.Y() ) ;
0732   v3X.SetZ( l_Ejectile_rfg.Z() );
0733
0734   v3S.SetX( -1 );
0735   v3S.SetY( 0 );
0736   v3S.SetZ( 0 );
0737
0738   v3PhotonUnit = v3Photon.Unit();
0739   v3QxL        = v3Photon.Cross(v3Electron);
0740   v3QxP        = v3Photon.Cross(v3X);
0741   v3QxS        = v3Photon.Cross(v3S);
0742   v3LxP        = v3Electron.Cross(v3X);
0743   v3LxS        = v3Electron.Cross(v3S);
0744   v3PxL        = v3X.Cross(v3Electron);
0745   v3QUnitxL    = v3PhotonUnit.Cross(v3Electron);
0746   v3QUnitxP    = v3PhotonUnit.Cross(v3X);
0747   v3QUnitxS    = v3PhotonUnit.Cross(v3S);
0748
0749   /*--------------------------------------------------*/
0750   // Get the Phi scattering angle with respect to the electron scattering plane
0751   fPhi   = Get_Phi_X_LeptonPlane_RF ();
0752
0753   /*--------------------------------------------------*/
0754   // Get the Phi scattering angle with respect to the electron scattering plane
0755   fPhiS  = Get_Phi_TargPol_LeptonPlane_RF();
0756
0757   fTheta_X_Photon_RF       = fRAD2DEG * acos( ( v3Photon.Dot( v3X     ) ) / ( v3Photon.Mag()  * v3X.Mag()    ) );
0758   if ( fTheta_X_Photon_RF < 0 ) { fTheta_X_Photon_RF = 180.0 + fTheta_X_Photon_RF; }
0759
0760   // -----------------------------------------------------------------------------------
0761   // If we have fermi momentum then epsilon should be in rest frame
0762   // The theta angle of scattered angle used in expression of epsilon is the angle
0763   // with respect to direction of incoming electron in the rest frame of target hadron
0764   // epsilon=1./(1.+ 2.*(pgam_restg**2)/q2g * *(tand(thscat_rest/2.))**2)
0765   // -----------------------------------------------------------------------------------
0766
0767   //double fTheta_EEp = (lelectron_rf.Vect()).Angle(lscatelec_rf.Vect());
0768   double fTheta_EEp = ( lelectron_rfg.Vect() ).Angle( lscatelec_rfg.Vect() ); // involved quantities are in GeV
0769
0770   fEpsilon = 1.0 / ( 1.0 + 2.0 * ( pow( (lphoton_rfg.Vect()).Mag(),2)/fQsq_GeV ) * pow( tan( fTheta_EEp / 2 ) , 2 ) );
0771
0772   // ----------------------------------------------------
0773   // Virtual Photon flux factor in units of 1/(GeV*Sr)
0774   // ----------------------------------------------------
0775   // fFlux_Factor_Col = (fAlpha/(2.0*pow(fPi,2))) * (r_lscatelecg.E() / r_lelectrong.E()) *
0776   //   ( pow(fW_GeV,2) - pow(fProton_Mass_GeV,2) ) / (2.0*fProton_Mass_GeV*fQsq_GeV*(1.0 - fEpsilon));
0777
0778   fFlux_Factor_Col = ( fAlpha / ( 2.0 * pow( fPi , 2 ) ) ) * ( r_lscatelecg.E() / r_lelectrong.E() ) *  // All are in GeV (more organised form)
0779     ( pow( fW_GeV , 2 ) - pow( fProton_Mass_GeV , 2 ) ) /
0780     ( 2.0 * fProton_Mass_GeV * fQsq_GeV * ( 1.0 - fEpsilon ) );
0781
0782   fFlux_Factor_RF = ( fAlpha / ( 2.0 * pow( fPi , 2 ) ) ) * ( lscatelec_rfg.E() / lelectron_rfg.E() ) *
0783     ( pow( fW_GeV , 2 ) - pow( fProton_Mass_GeV , 2 ) ) /
0784     ( 2.0 * fProton_Mass_GeV * fQsq_GeV * ( 1.0 - fEpsilon ) );
0785
0786   // cout<<"  "<<fFlux_Factor_Col<<"  "<<fFlux_Factor_RF<<endl;
0787
0788   // ----------------------------------------------------
0789   //  Jacobian  dt/dcos(theta*)dphi in units of GeV2/sr
0790   // ----------------------------------------------------
0791   // fJacobian_CM = ( (lphoton_rfg.Vect()).Mag() - fBeta_CM_RF * lphoton_rfg.E() ) / ( fGamma_CM_RF * ( 1.0 - pow(fBeta_CM_RF,2) ) ); // Eqn 22 in paper
0792   fJacobian_CM = fGamma_CM_RF * ( ( lphoton_rfg.Vect() ).Mag() - ( fBeta_CM_RF * lphoton_rfg.E() ) ); // All are in GeV // Love Preeet changed it
0793
0794   fA = fJacobian_CM * f_Ejectile_Mom_CM_GeV / fPi; // Eqn 21 in paper
0795
0796   // ----------------------------------------------------
0797   // Jacobian dOmega* / dOmega dimensionless
0798   // ----------------------------------------------------
0799   /*fJacobian_CM_RF  = ( pow((l_Ejectile_rf.Vect()).Mag(),2)*fW) /
0800     ( f_Ejectile_Mom_CM * std::abs( ( fProton_Mass + lphoton_rf.E()) * (l_Ejectile_rf.Vect()).Mag() -
0801     ( l_Ejectile_rf.E() * (lphoton_rf.Vect()).Mag() * cos( l_Ejectile_rf.Theta() ) ) ) ); // Differs from next line in photon vect -> lphoton_rf vs r_lphoton
0802
0803     fJacobian_CM_Col = ( ( pow((r_l_Ejectile.Vect()).Mag(),2) * fW ) / // This one is actually used subsequently, so this must be Eqn 20
0804     ( f_Ejectile_Mom_CM * std::abs( ( fProton_Mass + r_lphoton.E() ) * (r_l_Ejectile.Vect()).Mag() -
0805     ( r_l_Ejectile.E() * (r_lphoton.Vect()).Mag() * cos( r_l_Ejectile.Theta() ) ) ) ) ); */
0806
0807   fBeta_Col  = ( r_lphotong.Vect() + r_lprotong.Vect() ).Mag() / ( r_lphotong.E() + r_lprotong.E() ); // All are in GeV // Love Preeet changed it
0808   fGamma_Col = ( r_lphotong.E() + r_lprotong.E() ) / ( fW * fm );
0809   ftheta_Col = ( r_lphotong.Angle( r_l_Ejectile_g.Vect() ) ); // angle between the virtualphoton and the ejectile in the collider frame
0810
0811   fJacobian_CM_Col = ( pow( ( r_l_Ejectile_g.Vect()).Mag() , 2 ) /
0812                fGamma_Col * f_Ejectile_Mom_CM_GeV * ( ( r_l_Ejectile_g.Vect() ).Mag() - ( fBeta_Col * r_l_Ejectile_g.E() * cos( ftheta_Col ) ) ) );
0813
0814
0815   //     cout <<  l_Ejectile_rf.Vect().Mag() << "  " << << << << << << << << endl;
0816   //     cout << fJacobian_CM_RF << "    " << fJacobian_CM_Col << endl;
0817
0818   // -----------------------------------------------------------------------------------------------------------
0819   // CKY sigma L and T starts
0820   // -----------------------------------------------------------------------------------------------------------
0821   //     r_fSig_T = 1;
0822   //     r_fSig_L = 1;
0823   // -------------------------------------------------------------------------------------------
0824   r_fSig = Get_Total_Cross_Section();
0825
0826   // -----------------------------------------------------------------------------------------------------------
0827   // CKY sigma L and T ends
0828   // -----------------------------------------------------------------------------------------------------------
0829
0830   fSigma_Col = r_fSig * fFlux_Factor_Col * fA * fJacobian_CM_Col;
0831   // fSigma_Col = r_fSig * fFlux_Factor_RF * fA * fJacobian_CM_Col; // Love Preet changed flux factor from collider to rest frame
0832
0833   if ( ( fSigma_Col <= 0 ) || std::isnan( fSigma_Col ) ) {
0834     fNSigmaNeg ++;
0835     return;
0836   }
0837
0838   // -----------------------------------------------------------------------------------------------------------
0839   //             Lab cross section     Phase Space   Conversion     Luminosity                Total events tried
0840   // Hz        = ub / ( sr^2 * GeV ) * GeV * sr^2 * ( cm^2 / ub ) * ( # / ( cm^2 * sec ) ) / ( # )
0841
0842   // SJDK 24/06/21 - Explicitly taking the absolute value of the weight such that the value is positive! Shouldn't matter since any -ve cross section events should be dumped above
0843   //fEventWeight = abs(fSigma_Col * fPSF * fuBcm2 * fLumi / fNEvents);   // in Hz
0844   fEventWeight = fSigma_Col * fPSF * fuBcm2 * fLumi / fNEvents;   // in Hz // Love Preet removed the abs on the fEventWeight
0845
0846   if ( ( fEventWeight <= 0 ) || std::isnan( fEventWeight ) ) {  // Love Preet added new counter to track negative and NaN weights
0847     fNWeightNeg ++;
0848     return;
0849   }
0850
0851   fNRecorded++;
0852   fRatio = fNRecorded / fNGenerated;
0853   //Love Preet - Added for actual phase space factor calculations
0854   calculate_psf_max_min( fScatElec_Energy_Col, fScatElec_Energy_Col_max, fScatElec_Energy_Col_min ); // -> Love Preet added to find the max and min values to calculate the actual PSF
0855   calculate_psf_max_min( fScatElec_Theta_Col,  fScatElec_Theta_Col_max,  fScatElec_Theta_Col_min );
0856   calculate_psf_max_min( f_Ejectile_Theta_Col, f_Ejectile_Theta_Col_max, f_Ejectile_Theta_Col_min );
0857
0858   if (gOutputType == "Pythia6"){
0859     DEMPReact_Pythia6_Output();
0860   }
0861   else if (gOutputType == "LUND"){
0862     Lund_Output();
0863   }
0864   else if (gOutputType == "HEPMC3"){
0865     DEMPReact_HEPMC3_Output();
0866   }
0867
0868   if (gROOTOut == true){
0869     scat_e_px =  r_lscatelecg.X(); // Love Preet - Added to be stored in the root tree
0870     scat_e_py =  r_lscatelecg.Y();
0871     scat_e_pz =  r_lscatelecg.Z();
0872     scat_e_E  =  r_lscatelecg.E();
0873     ejec_px   =  r_l_Ejectile_g.X();
0874     ejec_py   =  r_l_Ejectile_g.Y();
0875     ejec_pz   =  r_l_Ejectile_g.Z();
0876     ejec_E    =  r_l_Ejectile_g.E();
0877     rclH_px   =  l_Recoil_g.X();
0878     rclH_py   =  l_Recoil_g.Y();
0879     rclH_pz   =  l_Recoil_g.Z();
0880     rclH_E    =  l_Recoil_g.E();
0881     dRootTree->Fill();
0882   }
0883 }
0884
0885 void DEMP_Reaction::Progress_Report(){
0886
0887   dFractTime = time(0);
0888
0889   if ( rNEvent_itt % ( rNEvents / 10 ) == 0 ) {
0890     cout << "Event: " << setw(8) << rNEvent_itt
0891      << "     % of events " << setw(4) << ((1.0*rNEvent_itt)/(1.0*rNEvents))*100.0
0892      << "   Day: " <<  dFractTime.GetDay()
0893      << "   Time:   " << dFractTime.GetHour()
0894      << ":" << dFractTime.GetMinute()
0895      << ":" << dFractTime.GetSecond()
0896      << endl;
0897   }
0898 }
0899
0900 TLorentzVector DEMP_Reaction::GetProtonVector_lab(){
0901
0902   // Crossing angle
0903   //     fProton_Theta_Col = 0.050;
0904   //     fProton_Theta_Col = 0.025;
0905   // SJDK - 12/01/22
0906   // Set crossing angle to 0 for fun4all, also required for ATHENA simulations
0907   fProton_Theta_Col = 0.0;
0908
0909   ///*--------------------------------------------------*/
0910   //     fProton_Phi_Col   = fPi;
0911   fProton_Phi_Col   = fProton_incidence_phi;
0912
0913   fProton_Mom_Col   = fHBeam * 1e3;
0914   fVertex_X         = 0.;
0915   fVertex_Y         = 0.;
0916   fVertex_Z         = 0.;
0917
0918   TLorentzVector lproton( fProton_Mom_Col * sin(fProton_Theta_Col) * cos(fProton_Phi_Col),
0919               fProton_Mom_Col * sin(fProton_Theta_Col) * sin(fProton_Phi_Col),
0920               fProton_Mom_Col * cos(fProton_Theta_Col),
0921               sqrt( pow( fProton_Mom_Col , 2 ) + pow( fProton_Mass , 2 ) ) );
0922
0923   return lproton;
0924
0925 }
0926
0927 //*--------------------------------------------------*/
0928 // Proton in collider (lab) frame
0929 // ----------------------------------------------------
0930
0931 void DEMP_Reaction::Consider_Proton_Fermi_Momentum(){
0932
0933   fProton_Mom_Col   = fProton_Mom_Col + rFermiMomentum;
0934   fProton_Theta_Col = acos( fRandom->Uniform( cos(0.0) , cos(fPi) ) );
0935   fProton_Phi_Col   = fRandom->Uniform( 0 , 360 );
0936
0937   double px, py, pz, e;
0938
0939   px = fProton_Mom_Col * sin(fProton_Theta_Col) * cos(fProton_Phi_Col);
0940   py = fProton_Mom_Col * sin(fProton_Theta_Col) * sin(fProton_Phi_Col);
0941   pz = fProton_Mom_Col * cos(fProton_Theta_Col);
0942   e  = sqrt( pow( fProton_Mom_Col , 2 ) + pow( fProton_Mass , 2 ) );
0943
0944   r_lproton.SetPxPyPzE(px,py,pz,e);
0945
0946   r_lprotong = r_lproton*fm;
0947
0948 }
0949
0950 // ----------------------------------------------------
0951 // Electron in collider (lab) frame
0952 // ----------------------------------------------------
0953
0954 TLorentzVector DEMP_Reaction::GetElectronVector_lab(){
0955
0956   fElectron_Energy_Col = fElectron_Kin_Col;
0957   fElectron_Mom_Col    = sqrt( pow(fElectron_Energy_Col , 2) - pow(fElectron_Mass , 2) );
0958   fElectron_Theta_Col  = fPi;
0959   fElectron_Phi_Col    = 0.0;
0960   fElectron_MomZ_Col   = fElectron_Mom_Col * cos(fElectron_Theta_Col);
0961   fElectron_MomX_Col   = fElectron_Mom_Col * sin(fElectron_Theta_Col) * cos(fElectron_Phi_Col);
0962   fElectron_MomY_Col   = fElectron_Mom_Col * sin(fElectron_Theta_Col) * sin(fElectron_Phi_Col);
0963
0964   //cout << "Define: " << fElectron_MomZ_Col << "    "<< fElectron_Mom_Col << "  " << cos(fElectron_Theta_Col) << endl;
0965
0966   TLorentzVector  lelectron( fElectron_MomX_Col, fElectron_MomY_Col, fElectron_MomZ_Col, fElectron_Energy_Col);
0967
0968   return lelectron;
0969
0970 }
0971
0972 Double_t DEMP_Reaction::Get_Phi_X_LeptonPlane_RF(){
0973
0974   fCos_Phi_X_LeptonPlane_RF = ( ( v3QUnitxL.Dot( v3QUnitxP ) ) / ( v3QUnitxL.Mag() * v3QUnitxP.Mag() ) ); // hep-ph/0410050v2
0975   fSin_Phi_X_LeptonPlane_RF = ( ( v3LxP.Dot( v3PhotonUnit  ) ) / ( v3QUnitxL.Mag() * v3QUnitxP.Mag() ) ); // hep-ph/0410050v2
0976   if ( fSin_Phi_X_LeptonPlane_RF >= 0 )
0977     fPhi_X_LeptonPlane_RF   = fRAD2DEG * acos( ( v3QUnitxL.Dot( v3QUnitxP ) ) / ( v3QUnitxL.Mag() * v3QUnitxP.Mag() ) );
0978   if ( fSin_Phi_X_LeptonPlane_RF < 0 )
0979     fPhi_X_LeptonPlane_RF   = 360.0 - std::abs( fRAD2DEG * acos( ( v3QUnitxL.Dot( v3QUnitxP ) ) / ( v3QUnitxL.Mag() * v3QUnitxP.Mag() ) ) );
0980
0981   return fPhi_X_LeptonPlane_RF;
0982
0983 }
0984
0985 Double_t DEMP_Reaction::Get_Phi_TargPol_LeptonPlane_RF(){
0986
0987   fCos_Phi_TargPol_LeptonPlane_RF = ( ( v3QUnitxL.Dot( v3QUnitxS ) ) / ( v3QUnitxL.Mag() * v3QUnitxS.Mag() ) ); // hep-ph/0410050v2
0988   fSin_Phi_TargPol_LeptonPlane_RF = ( ( v3LxS.Dot( v3PhotonUnit  ) ) / ( v3QUnitxL.Mag() * v3QUnitxS.Mag() ) ); // hep-ph/0410050v2
0989   if ( fSin_Phi_TargPol_LeptonPlane_RF >= 0 )
0990     fPhi_TargPol_LeptonPlane_RF = fRAD2DEG * acos( ( v3QUnitxL.Dot( v3QUnitxS ) ) / ( v3QUnitxL.Mag() * v3QUnitxS.Mag() ) );
0991   if ( fSin_Phi_TargPol_LeptonPlane_RF < 0 )
0992     fPhi_TargPol_LeptonPlane_RF = 360.0 - std::abs( fRAD2DEG * acos( ( v3QUnitxL.Dot( v3QUnitxS ) ) / ( v3QUnitxL.Mag() * v3QUnitxS.Mag() ) ) );
0993
0994   return fPhi_TargPol_LeptonPlane_RF;
0995
0996 }
0997
0998 Double_t DEMP_Reaction::Get_Total_Cross_Section(){
0999
1000   Double_t total_sig, total_sig2;
1001
1002   if (rEjectile == "Pi+"){
1003     total_sig = GetPiPlus_CrossSection(fT_GeV, fW_GeV, fQsq_GeV, fEpsilon);
1004   }
1005   else if (rEjectile == "Pi0"){
1006     total_sig = GetKPlus_CrossSection(fT_GeV, fW_GeV, fQsq_GeV, fEpsilon, rRecoil);
1007   }
1008   else if (rEjectile == "K+"){
1009     total_sig = GetKPlus_CrossSection(fT_GeV, fW_GeV, fQsq_GeV, fEpsilon, rRecoil);
1010     // SJDK - 31/01/23 - This second case gets the total cross section for the K+ as calculated from Ali's earlier scaling calculation, retain this for comparison.
1011     //total_sig2 = GetKPlus_CrossSection_Scaling(fT_GeV, fW_GeV, fQsq_GeV, fEpsilon, f_Ejectile_Mass_GeV, rRecoil);
1012   }
1013
1014   //cout << fT_GeV <<  "  " << fW_GeV << "   " << fQsq_GeV << "   " << "KPlus Paramterisation - " << total_sig << " !!! Old Scaling method - " << total_sig2 << endl;
1015
1016   return total_sig;
1017 }
1018
1019 //// SJDK 21/12/22 - This function needs updating!
1020 //Double_t DEMP_Reaction::GetPi0_CrossSection() {
1021 //
1022 //  double_t sig_total;
1023 //  return sig_total;
1024 //
1025 //}
1026
1027 /*--------------------------------------------------*/
1028 /// Output generator detail
1029 // 06/09/23 SJDK - Cuts are now ordered as they are applied in the generator
1030 void DEMP_Reaction::Detail_Output(){ // Love Preet changed the order of conditions and added new statements
1031
1032   DEMPDetails << left << setw(70) << "Seed used for the Random Number Generator" << right << setw(20) << fSeed << endl;
1033   DEMPDetails << endl;
1034   DEMPDetails << left << setw(70) << "Total events tried" << right << setw(20) << fNGenerated << endl;
1035   if(UseSolve == true){
1036     DEMPDetails << left << setw(70) << "Total events cut" << right << setw(20) << (fNaN + fConserve + w_neg_ev + qsq_ev + w_ev  + t_ev + fNSigmaNeg + fNWeightNeg + fSolveEvents_0Sol) << right << setw(20) << ((double) (fNaN + fConserve + w_neg_ev + qsq_ev + w_ev  + t_ev + fNSigmaNeg + fNWeightNeg + fSolveEvents_0Sol)/(double)fNGenerated)*100 << " %" << endl;
1037     DEMPDetails << left << setw(70) << "Total events recorded" << right << setw(20) << fNRecorded << right << setw(20) << ((double)fNRecorded/(double)fNGenerated)*100 << " %" << endl;
1038     if (fNGenerated != (fNaN + fConserve + w_neg_ev + qsq_ev + w_ev  + t_ev + fNSigmaNeg + fNWeightNeg + fSolveEvents_0Sol + fNRecorded)){
1039       DEMPDetails << left << setw(70) << "Total events cut + recorded = events tried?" << right << setw(20) << "NO! ERROR!" << endl;
1040     }
1041     else{
1042       DEMPDetails << left << setw(70) << "Total events cut + recorded = events tried?" << right << setw(20) << "Yes! :)" << endl;
1043     }
1044   }
1045   else{
1046     DEMPDetails << left << setw(70) << "Total events cut" << right << setw(20) << (fNaN + fConserve + w_neg_ev + qsq_ev + w_ev  + t_ev + fNSigmaNeg + fNWeightNeg) << right << setw(20) << ((double) (fNaN + fConserve + w_neg_ev + qsq_ev + w_ev  + t_ev + fNSigmaNeg + fNWeightNeg)/(double)fNGenerated)*100 << " %" << endl;
1047     DEMPDetails << left << setw(70) << "Total events recorded" << right << setw(20) << fNRecorded << right << setw(20) << ((double)fNRecorded/(double)fNGenerated)*100 << " %" << endl;
1048     if (fNGenerated != (fNaN + fConserve + w_neg_ev + qsq_ev + w_ev  + t_ev + fNSigmaNeg + fNWeightNeg + fNRecorded)){
1049       DEMPDetails << left << setw(70) << "Total events cut + recorded = events tried?" << right << setw(20) << "NO! ERROR!" << endl;
1050     }
1051     else{
1052       DEMPDetails << left << setw(70) << "Total events cut + recorded = events tried?" << right << setw(20) << "Yes! :)" << endl;
1053     }
1054   }
1055
1056   DEMPDetails << left << setw(70) << endl << "Cut details -" << endl;
1057   DEMPDetails << left << setw(70) << "Events cut due to ejectile (X) energy NaN" << right << setw(20) << fNaN << right << setw(20) << ((double)fNaN/(double)fNGenerated)*100 << " %" << endl;
1058   DEMPDetails << left << setw(70) << "Events cut due to conservation law check failure" << right << setw(20) << fConserve << right << setw(20) << ((double)fConserve/(double)fNGenerated)*100 << " %" << endl;
1059   DEMPDetails << left << setw(70) << "Events cut due to negative Wsq value " << right << setw(20) << w_neg_ev << right << setw(20) << ((double)w_neg_ev/(double)fNGenerated)*100 << " %" << endl;
1060   DEMPDetails << left << setw(70) << Form("Events cut due to qsq < %.1lf or qsq > %.1lf", fQsq_Min, fQsq_Max) << right << setw(20) << qsq_ev << right << setw(20) << ((double)qsq_ev/(double)fNGenerated)*100 << " %" << endl;
1061   DEMPDetails << left << setw(70) << Form("Events cut due to W < %.1lf or W > %.1lf", fW_Min, fW_Max) << right << setw(20) << w_ev << right << setw(20) << ((double)w_ev/(double)fNGenerated)*100 << " %" << endl;
1062   if(UseSolve == true){
1063     DEMPDetails << left << setw(70) << "Events cut due to solve function finding 0 solutions" << right << setw(20) << fSolveEvents_0Sol << right << setw(20) << ((double)fSolveEvents_0Sol/(double)fNGenerated)*100 << " %" << endl;
1064   }
1065
1066   DEMPDetails << left << setw(70) << Form("Events cut due to -t > %.1lf GeV", fT_Max) << right << setw(20) << t_ev << right << setw(20) << ((double)t_ev/(double)fNGenerated)*100 << " %" << endl;
1067   DEMPDetails << left << setw(70) << "Events cut due to -ve cross section value" << right << setw(20) << fNSigmaNeg << right << setw(20) << ((double)fNSigmaNeg/(double)fNGenerated)*100 << " %" << endl;
1068   DEMPDetails << left << setw(70) << "Events cut due to -ve weight value" << right << setw(20) << fNWeightNeg << right << setw(20) << ((double)fNWeightNeg/(double)fNGenerated)*100 << " %" << endl;
1069
1070   DEMPDetails << left << setw(70) << endl << "Conservation law checks details -" << endl;
1071   DEMPDetails << left << setw(70) << Form("Total events PASSING conservation law check with tolerance %.2e", fDiff) << right << setw(20) << conserve << endl;
1072   DEMPDetails << left << setw(70) << "Events cut due to energy conservation check ONLY" << right << setw(20) << ene << right << setw(20) << ((double)ene/(double)fNGenerated)*100 << " %" << endl;
1073   DEMPDetails << left << setw(70) << "Events cut due to momentum conservation check ONLY" << right << setw(20) << mom << right << setw(20) << ((double)mom/(double)fNGenerated)*100 << " %" << endl;
1074   DEMPDetails << left << setw(70) << "Events cut due to energy AND momentum conservation checks" << right << setw(20) << ene_mom << right << setw(20) << ((double)ene_mom/(double)fNGenerated)*100 << " %" << endl;
1075   DEMPDetails << left << setw(70) << "Events cut due to px conservation law check" << right << setw(20) << mom_px << right << setw(20) << ((double)mom_px/(double)fNGenerated)*100 << " %" << endl;
1076   DEMPDetails << left << setw(70) << "Events cut due to py conservation law check" << right << setw(20) << mom_py << right << setw(20) << ((double)mom_py/(double)fNGenerated)*100 << " %" << endl;
1077   DEMPDetails << left << setw(70) << "Events cut due to pz conservation law check" << right << setw(20) << mom_pz << right << setw(20) << ((double)mom_pz/(double)fNGenerated)*100 << " %" << endl;
1078   DEMPDetails << left << setw(70) << "Events cut due to px and py conservation law checks" << right << setw(20) << mom_pxpy << right << setw(20) << ((double)mom_pxpy/(double)fNGenerated)*100 << " %" << endl;
1079   DEMPDetails << left << setw(70) << "Events cut due to px and pz conservation law checks" << right << setw(20) << mom_pxpz << right << setw(20) << ((double)mom_pxpz/(double)fNGenerated)*100 << " %" << endl;
1080   DEMPDetails << left << setw(70) << "Events cut due to py and pz conservation law checks" << right << setw(20) << mom_pypz << right << setw(20) << ((double)mom_pypz/(double)fNGenerated)*100 << " %" << endl;
1081   DEMPDetails << left << setw(70) << "Events cut due to px, py and pz conservation law checks" << right << setw(20) << mom_pxpypz << right << setw(20) << ((double)mom_pxpypz/(double)fNGenerated)*100 << " %" << endl;
1082
1083   DEMPDetails << left << setw(70) << endl << "Weight correction factors -" << endl;
1084   DEMPDetails << left << setw(70) << "Ratio of phase space factors (fPSF_org/fPSF)" << right << setw(20) << ((double)fPSF_org/(double)fPSF) << endl;
1085   DEMPDetails << left << setw(70) << "Ratio of tried to physical events" << right << setw(20) << ((double)fNGenerated / (double) (fNGenerated - fNaN - fConserve - w_neg_ev))  << endl;
1086   DEMPDetails << left << setw(70) <<" 1. If the first ratio is not ~1.0 (+/- 0.05), check the number of recorded events. If the number of recorded events is < 50,000, increase the number of attempted events to generate more data. If the the ratio does not converge to ~1.0, multiply all individual event weights by this factor during analysis. Alternatively, this number can be treated as a tolerance on weights."<<endl;
1087   DEMPDetails << left << setw(70) <<" 2. If the second ratio is large (> XX), then again, multiply all invidiual event weights by this factor during analysis. Again, this could also be considered as a tolerance on weights."<<endl;
1088
1089   DEMPDetails << left << setw(70) << endl << "Energies, angles, and phase space factors -" << endl;
1090
1091   DEMPDetails << right << setw(90) << "User-defined" << right << setw(20) << "Calculated" << right <<  setw(20) << "Actual" << endl;
1092
1093   DEMPDetails << left << setw(70) << "Scattered electron energies (Low)" << right << setw(20) << fScatElec_E_Lo * fElectron_Energy_Col *fm << right << setw(20) << psf_ScatElec_E_min << right << setw(20) << fScatElec_Energy_Col_min * fm << endl;
1094
1095   DEMPDetails << left << setw(70) << "Scattered electron energies (High)" << right << setw(20) << fScatElec_E_Hi * fElectron_Energy_Col *fm << right << setw(20) << psf_ScatElec_E_max << right << setw(20) << fScatElec_Energy_Col_max * fm << endl;
1096
1097   DEMPDetails << left << setw(70) << "Scattered electron angles (Low)" << right << setw(20) << (fScatElec_Theta_I * TMath::RadToDeg()) << right << setw(20) << psf_ScatElec_Theta_min * TMath::RadToDeg() << right << setw(20) << fScatElec_Theta_Col_min * TMath::RadToDeg() << endl;
1098
1099   DEMPDetails << left << setw(70) << "Scattered electron angles (High)" << right << setw(20) << (fScatElec_Theta_F * TMath::RadToDeg()) << right << setw(20) << psf_ScatElec_Theta_max * TMath::RadToDeg() << right << setw(20) << fScatElec_Theta_Col_max * TMath::RadToDeg() << endl;
1100
1101   DEMPDetails << left << setw(70) << "Scattered ejectile angles (Low)" << right << setw(20) << (f_Ejectile_Theta_I * TMath::RadToDeg()) << right << setw(20) << psf_Ejectile_Theta_min * TMath::RadToDeg() << right << setw(20) << f_Ejectile_Theta_Col_min * TMath::RadToDeg() << endl;
1102
1103   DEMPDetails << left << setw(70) << "Scattered ejectile angles (High)" << right << setw(20) << (f_Ejectile_Theta_F * TMath::RadToDeg()) << right << setw(20) << psf_Ejectile_Theta_max * TMath::RadToDeg() << right << setw(20) << f_Ejectile_Theta_Col_max * TMath::RadToDeg() << endl;
1104
1105   DEMPDetails << left << setw(70) << "Phase space factors (fPSF fSF_org)" << right << setw(20) << fPSF << right << setw(20) << fPSF_org << endl;
1106
1107   if(UseSolve == true){
1108     DEMPDetails << left << setw(70) << endl << "Solve function, addtional info -" << endl;
1109     DEMPDetails << left << setw(70) << "Number of events with 0 Solution" << right << setw(20) << fSolveEvents_0Sol << endl;
1110     DEMPDetails << left << setw(70) << "Number of events with 1 Solution" << right << setw(20) << fSolveEvents_1Sol << endl;
1111     DEMPDetails << left << setw(70) << "Number of events with 2 Solution" << right << setw(20) << fSolveEvents_2Sol << endl;
1112   }
1113 }
1114
1115 ////*--------------------------------------------------
1116 /// Functions for different output formats follow
1117
1118 void DEMP_Reaction::Lund_Output(){
1119
1120   DEMPOut << "3"
1121       << " \t " << fPhi           // var 1
1122       << " \t " << fPhiS          // var 2
1123       << " \t " << fx             // var 3
1124       << " \t " << "1"
1125       << " \t " << fQsq_GeV       // var 4
1126       << " \t " << fT_GeV         // var 5
1127       << " \t " << fW_GeV          // var 6
1128       << " \t " << fEpsilon       // var 7
1129       << " \t " << fEventWeight   // var 8
1130       << endl;
1131
1132   // Produced Particle X
1133   DEMPOut << setw(10) << "1"
1134       << setw(10) << "1"
1135       << setw(10) << "1"
1136       << setw(10) << PDGtype(produced_X)
1137       << setw(10) << "0"
1138       << setw(10) << "0"
1139       << setw(16) << r_l_Ejectile_g.X()
1140       << setw(16) << r_l_Ejectile_g.Y()
1141       << setw(16) << r_l_Ejectile_g.Z()
1142       << setw(16) << r_l_Ejectile_g.E()
1143       << setw(16) << r_l_Ejectile_g.M() //15/05/23 - Love - Was fX_Mass_GeV
1144       << setw(16) << fVertex_X
1145       << setw(16) << fVertex_Y
1146       << setw(16) << fVertex_Z
1147       << endl;
1148
1149   // Scattered electron
1150   DEMPOut << setw(10) << "2"
1151       << setw(10) << "-1"
1152       << setw(10) << "1"
1153       << setw(10) << "11"
1154       << setw(10) << "0"
1155       << setw(10) << "0"
1156       << setw(16) << r_lscatelecg.X()
1157       << setw(16) << r_lscatelecg.Y()
1158       << setw(16) << r_lscatelecg.Z()
1159       << setw(16) << r_lscatelecg.E()
1160       << setw(16) << r_lscatelecg.M() //15/05/23 - Love- Was fElectron_Mass_GeV
1161       << setw(16) << fVertex_X
1162       << setw(16) << fVertex_Y
1163       << setw(16) << fVertex_Z
1164       << endl;
1165
1166   // Recoiled neutron
1167   DEMPOut << setw(10) << "3"
1168       << setw(10) << "1"
1169       << setw(10) << "1"
1170       << setw(10) << PDGtype(recoil_hadron)
1171       << setw(10) << "0"
1172       << setw(10) << "0"
1173       << setw(16) << l_Recoil_g.X()
1174       << setw(16) << l_Recoil_g.Y()
1175       << setw(16) << l_Recoil_g.Z()
1176       << setw(16) << l_Recoil_g.E()
1177       << setw(16) << l_Recoil_g.M() // 15/05/23 - Love - Was f_Scat_hadron_Mass_GeV
1178       << setw(16) << fVertex_X
1179       << setw(16) << fVertex_Y
1180       << setw(16) << fVertex_Z
1181       << endl;
1182 }
1183
1184 void DEMP_Reaction::DEMPReact_Pythia6_Out_Init(){
1185
1186   print_itt = 0;
1187
1188   DEMPOut << "DEMP Event FILE" << endl;
1189   DEMPOut << "============================================" << endl;
1190   DEMPOut << "I, ievent, nParticles, Weight" << endl;
1191   DEMPOut << "============================================" << endl;
1192   DEMPOut << "I  K(I,1)  K(I,2)  K(I,3)  K(I,4)  K(I,5)  P(I,1)  P(I,2)  P(I,3)  P(I,4)  P(I,5)  V(I,1)  V(I,2)  V(I,3)" << endl;
1193   DEMPOut << "============================================" << endl;
1194
1195 }
1196
1197 void DEMP_Reaction::DEMPReact_Pythia6_Output(){
1198
1199   DEMPOut << "0" << " \t\t\t "  << print_itt << " \t\t\t " << "1" << " \t\t\t " << fEventWeight << endl; // var 1
1200
1201   print_itt++;
1202
1203   DEMPOut << "============================================" << endl;
1204
1205   ///*--------------------------------------------------*/
1206   // Initial State
1207
1208   DEMPOut  << "1"
1209        << setw(6) << "21"
1210        << setw(6) << "11"
1211        << setw(6) << "0"
1212        << setw(6) << "3"
1213        << setw(6) << "4"
1214
1215        << setw(14) << r_lelectrong.X()
1216        << setw(14) << r_lelectrong.Y()
1217        << setw(14) << r_lelectrong.Z()
1218        << setw(14) << r_lelectrong.E()
1219        << setw(14) <<  r_lelectrong.M() // 15/05/23 - Love - Was fElectron_Mass_GeV
1220        << setw(6) << fVertex_X
1221        << setw(6) << fVertex_Y
1222        << setw(6) << fVertex_Z
1223        << endl;
1224
1225   DEMPOut << "2"
1226       << setw(6) << "21"
1227       << setw(6) << "2212"
1228       << setw(6) << "0"
1229       << setw(6) << "5"
1230       << setw(6) << "6"
1231
1232       << setw(14) << r_lprotong.X()
1233       << setw(14) << r_lprotong.Y()
1234       << setw(14) << r_lprotong.Z()
1235       << setw(14) << r_lprotong.E()
1236       << setw(14) << r_lprotong.M() // 15/05/23 - Love - Was fProton_Mass_GeV
1237       << setw(6) << fVertex_X
1238       << setw(6) << fVertex_Y
1239       << setw(6) << fVertex_Z
1240       << endl;
1241
1242   DEMPOut << "3"
1243       << setw(6) << "21"
1244       << setw(6) << "22"
1245       << setw(6) << "1"
1246       << setw(6) << "0"
1247       << setw(6) << "0"
1248
1249       << setw(14) << r_lphotong.X()
1250       << setw(14) << r_lphotong.Y()
1251       << setw(14) << r_lphotong.Z()
1252       << setw(14) << r_lphotong.E()
1253       << setw(14) << r_lphotong.M()
1254       << setw(6) << fVertex_X
1255       << setw(6) << fVertex_Y
1256       << setw(6) << fVertex_Z
1257       << endl;
1258
1259   ///*--------------------------------------------------*/
1260   // Final State
1261
1262   // Scattered electron
1263   DEMPOut << "4"
1264       << setw(6) << "1"
1265       << setw(6) << "11"
1266       << setw(6) << "1"
1267       << setw(6) << "0"
1268       << setw(6) << "0"
1269
1270       << setw(14) << r_lscatelecg.X()
1271       << setw(14) << r_lscatelecg.Y()
1272       << setw(14) << r_lscatelecg.Z()
1273       << setw(14) << r_lscatelecg.E()
1274       << setw(14) << r_lscatelecg.M() // 15/05/23 - Love - Was fElectron_Mass_GeV
1275       << setw(6) << fVertex_X
1276       << setw(6) << fVertex_Y
1277       << setw(6) << fVertex_Z
1278       << endl;
1279
1280   // Recoiled hadron
1281   DEMPOut << "5"
1282       << setw(6) << "1"
1283       << setw(6) << PDGtype(recoil_hadron)
1284       << setw(6) << "2"
1285       << setw(6) << "0"
1286       << setw(6) << "0"
1287
1288       << setw(14) << l_Recoil_g.X()
1289       << setw(14) << l_Recoil_g.Y()
1290       << setw(14) << l_Recoil_g.Z()
1291       << setw(14) << l_Recoil_g.E()
1292       << setw(14) <<  l_Recoil_g.M() // 15/05/23 - Love - Was f_Scat_hadron_Mass_GeV
1293       << setw(6) << fVertex_X
1294       << setw(6) << fVertex_Y
1295       << setw(6) << fVertex_Z
1296       << endl;
1297
1298   // Produced Particle X
1299   DEMPOut << "6"
1300       << setw(6) << "1"
1301       << setw(6) << PDGtype(produced_X)
1302       << setw(6) << "2"
1303       << setw(6) << "0"
1304       << setw(6) << "0"
1305
1306       << setw(14) << r_l_Ejectile_g.X()
1307       << setw(14) << r_l_Ejectile_g.Y()
1308       << setw(14) << r_l_Ejectile_g.Z()
1309       << setw(14) << r_l_Ejectile_g.E()
1310       << setw(14) <<  r_l_Ejectile_g.M() // 15/05/23 - Love - Was fX_Mass_GeV
1311       << setw(6) << fVertex_X
1312       << setw(6) << fVertex_Y
1313       << setw(6) << fVertex_Z
1314       << endl;
1315
1316   DEMPOut << "=============== Event finished ===============" << endl;
1317
1318 }
1319
1320 /*--------------------------------------------------*/
1321
1322 void DEMP_Reaction::DEMPReact_HEPMC3_Out_Init(){
1323
1324   print_itt = 0;
1325   DEMPOut << "HepMC::Version 3.02.02" << endl;
1326   DEMPOut << "HepMC::Asciiv3-START_EVENT_LISTING" << endl;
1327   // Book 1 weight
1328   DEMPOut << "W" << " " << "1" << endl;
1329   // Name it "weight"
1330   DEMPOut << "N" << " " << "weight" << endl;
1331
1332 }
1333
1334 /*--------------------------------------------------*/
1335
1336 void DEMP_Reaction::DEMPReact_HEPMC3_Output(){
1337
1338   // HEPMC3 output for Athena/ePIC simulations
1339   // First line - E - Event# - #Vertices - #Particles
1340   DEMPOut << std::scientific << std::setprecision(15) << "E" << " "  << print_itt <<  " " << "1" << " " << 5 << endl;
1341   print_itt++;
1342   // Second line, Units - U - ENERGY UNIT - DISTANCE UNIT
1343   DEMPOut << "U" << " " << "GEV" << " " << "MM" << endl;
1344   // Third line, optional attributes, the weight (LEGACY)
1345   DEMPOut << "A" << " " << "0" << " " << "weight" << " " <<  fEventWeight << endl;
1346   // Weight value
1347   DEMPOut << "W" << " " << fEventWeight << endl;
1348   // Beam particles, particle line - P - Particle ID - Parent Vertex ID - PDG id - px - py - pz - energy - particle mass - status (4, incoming beam particle)
1349   DEMPOut << "P" << " " << "1" << " " << "0" << " " << "11" << " " << r_lelectrong.X() << " " << r_lelectrong.Y() << " " << r_lelectrong.Z() << " " << r_lelectrong.E() << " " << r_lelectrong.M() << " " << "4" << endl;
1350   DEMPOut << "P" << " " << "2" << " " << "0" << " " << "2212" << " " << r_lprotong.X() << " " << r_lprotong.Y() << " " << r_lprotong.Z() << " " << r_lprotong.E() << " " <<  r_lprotong.M()<< " " << "4" << endl;
1351   // Vertex line - V - 1 - 0 - [1,2]
1352   DEMPOut << "V" << " " << "-1" << " " << "0" << " " << "[1,2]" << endl;
1353   // Output particles, particle line - P - Particle ID - Parent Vertex ID - PDG id - px - py - pz - energy - particle mass - status (1, undecayed physical particle)
1354   // Scattered electron
1355   DEMPOut << "P" << " " << "3" << " " << "-1" << " " << "11" << " " << r_lscatelecg.X() << " "  << r_lscatelecg.Y() << " "  << r_lscatelecg.Z() << " " << r_lscatelecg.E() << " " << r_lscatelecg.M() << " " << "1" << endl;
1356   // Produced meson
1357   DEMPOut << "P" << " " << "4" << " " << "-1" << " " << PDGtype(produced_X) << " " << r_l_Ejectile_g.X() << " "  << r_l_Ejectile_g.Y() << " "  << r_l_Ejectile_g.Z() << " " << r_l_Ejectile_g.E() << " " << r_l_Ejectile_g.M() << " " << "1" << endl;
1358   // Recoil hadron
1359   DEMPOut << "P" << " " << "5" << " " << "-1" << " " << PDGtype(recoil_hadron) << " " << l_Recoil_g.X() << " "  << l_Recoil_g.Y() << " "  << l_Recoil_g.Z() << " " << l_Recoil_g.E() << " " <<  l_Recoil_g.M() << " " << "1" << endl;
1360
1361 }
1362
1363 /*--------------------------------------------------*/
1364 bool DEMP_Reaction::SolnCheck(){
1365
1366   //  // Double Checking for solution viability
1367   //  if (TMath::Abs(f_Scat_hadron_Mass-r_l_scat_hadron_solved->M())>1){
1368   //    //cerr << "Mass Missmatch" << endl;
1369   //    //cerr << TMath::Abs(proton_mass_mev-Proton->M()) << endl;
1370   //    return false;
1371   //  }
1372   //  if (TMath::Abs(W_in()-W_out())>1){
1373   //    //cerr << "W Missmatch" << endl;
1374   //    //cerr << TMath::Abs(W_in()-W_out()) << endl;
1375   //    return false;
1376   //  }
1377   //  *Final = *r_l_scat_hadron_solved + *r_lX_solved;
1378   //
1379   //  if (TMath::Abs(Initial->Px()-Final->Px())>1){
1380   //    //cerr << "Px Missmatch" << endl;
1381   //    //cerr << TMath::Abs(Initial->Px()-Final->Px()) << endl;
1382   //    return false;
1383   //  }
1384   //
1385   //  if (TMath::Abs(Initial->Py()-Final->Py())>1){
1386   //    //cerr << "Py Missmatch" << endl;
1387   //    //cerr << TMath::Abs(Initial->Py()-Final->Py()) << endl;
1388   //    return false;
1389   //  }
1390   //
1391   //  if (TMath::Abs(Initial->Pz()-Final->Pz())>1){
1392   //    //cerr << "Pz Missmatch" << endl;
1393   //    //cerr << TMath::Abs(Initial->Pz()-Final->Pz()) << endl;
1394   //    return false;
1395   //  }
1396   //
1397   //  if (TMath::Abs(Initial->E()-Final->E())>1){
1398   //    return false;
1399   //  }
1400   return true;
1401 }
1402
1403 /*--------------------------------------------------*/
1404 double DEMP_Reaction::W_in(){
1405   return 0;
1406 }
1407
1408 /*--------------------------------------------------*/
1409 double DEMP_Reaction::W_out(){
1410   return 0;
1411 }
1412
1413 /*--------------------------------------------------*/
1414
1415 int DEMP_Reaction::Solve(){
1416
1417   VertBeamElec->SetPxPyPzE(r_lelectron.Px(), r_lelectron.Py(), r_lelectron.Pz(), r_lelectron.E());
1418   VertScatElec->SetPxPyPzE(r_lscatelec.Px(), r_lscatelec.Py(), r_lscatelec.Pz(), r_lscatelec.E());
1419   Target->SetPxPyPzE(r_lproton.Px(), r_lproton.Py(), r_lproton.Pz(), r_lproton.E());
1420   *Photon = *VertBeamElec - *VertScatElec;
1421   *Interaction = *Photon;
1422
1423   *Initial = *Interaction+*Target;
1424
1425   theta =  f_Ejectile_Theta_Col;
1426   phi   =  f_Ejectile_Phi_Col;
1427
1428   return this->Solve(theta, phi);
1429 }
1430
1431
1432 int DEMP_Reaction::Solve(double theta, double phi){
1433
1434   W_in_val = W_in();
1435
1436   if (W_in_val<0){
1437     return 0;
1438   }
1439
1440   UnitVect->SetTheta(theta);
1441   UnitVect->SetPhi(phi);
1442   UnitVect->SetMag(1);
1443
1444   double* pars = new double[9];
1445
1446   pars[0] = UnitVect->X();
1447   pars[1] = UnitVect->Y();
1448   pars[2] = UnitVect->Z();
1449   pars[3] = Initial->Px();
1450   pars[4] = Initial->Py();
1451   pars[5] = Initial->Pz();
1452   pars[6] = Initial->E();
1453   pars[7] = f_Ejectile_Mass;
1454   pars[8] = f_Recoil_Mass;
1455
1456   F->SetParameters(pars);
1457
1458   ///*--------------------------------------------------*/
1459   // Looking for the 1st Solution:
1460   //    If a solution found, then this will be the fist solution. Then we proceed to look for the 2nd solution.
1461   //    If no soluion found, then exit solve function
1462
1463   P = F->GetX(0, 0, pars[6], 0.0001, 10000);
1464
1465   if (TMath::Abs(F->Eval(P)) < 1){
1466     fSolveEvents_1Sol++;
1467   } else {
1468     fSolveEvents_0Sol++;
1469     return 0;
1470   }
1471
1472   TLorentzVector * r_l_Ejectile_solved_1_temp = new TLorentzVector();
1473   TLorentzVector * r_l_Ejectile_solved_2_temp = new TLorentzVector();
1474
1475   Float_t r_l_Ejectile_E = sqrt( pow(P*pars[0],2) + pow(P*pars[1],2) + pow(P*pars[2],2) + pow(f_Ejectile_Mass,2) );
1476   r_l_Ejectile_solved_1_temp->SetPxPyPzE(P*pars[0], P*pars[1],  P*pars[2], r_l_Ejectile_E);
1477
1478   ///*--------------------------------------------------*/
1479   // Looking for the 2nd Solution
1480
1481   P2 = F->GetX(0, P+100, pars[6], 0.0001, 10000);
1482   Float_t r_l_Ejectile_E_2 = sqrt( pow(P2 * pars[0],2) + pow(P2 * pars[1],2) + pow(P2 * pars[2],2) + pow(f_Ejectile_Mass,2) );
1483   r_l_Ejectile_solved_2_temp->SetPxPyPzE(P2 * pars[0], P2 * pars[1],  P2 * pars[2], r_l_Ejectile_E_2);
1484
1485   ///*--------------------------------------------------*/
1486   // If a valid 2nd solution is found, then we are certian that there are two solutions.
1487   //   - We then increament the counter for 2nd solution scenario
1488   //   - We then decreament the counter for the 1st solution scenario
1489
1490   if (TMath::Abs(F->Eval(P2)) < 1){
1491     fSolveEvents_2Sol++;
1492     fSolveEvents_1Sol--;
1493     if ( Int_t(CoinToss->Uniform(0,100)) < 50) {
1494       r_l_Ejectile_solved.SetPxPyPzE(r_l_Ejectile_solved_1_temp->X(), r_l_Ejectile_solved_1_temp->Y(), r_l_Ejectile_solved_1_temp->Z(), r_l_Ejectile_solved_1_temp->E());
1495     } else {
1496       r_l_Ejectile_solved.SetPxPyPzE(r_l_Ejectile_solved_2_temp->X(), r_l_Ejectile_solved_2_temp->Y(), r_l_Ejectile_solved_2_temp->Z(), r_l_Ejectile_solved_2_temp->E());
1497     }
1498   }
1499   else {
1500     r_l_Ejectile_solved.SetPxPyPzE(r_l_Ejectile_solved_1_temp->X(), r_l_Ejectile_solved_1_temp->Y(), r_l_Ejectile_solved_1_temp->Z(), r_l_Ejectile_solved_1_temp->E());
1501   }
1502
1503   ///*--------------------------------------------------*/
1504   /// Solve for the recoil information with the "solved" Ejectile informaiton
1505   TLorentzVector * r_l_hadron_temp= new TLorentzVector();
1506   *r_l_hadron_temp = *Initial- r_l_Ejectile_solved;
1507   r_l_Recoil_solved.SetPxPyPzE(r_l_hadron_temp->Px(), r_l_hadron_temp->Py(), r_l_hadron_temp->Pz(), r_l_hadron_temp->E());
1508
1509   delete r_l_Ejectile_solved_1_temp;
1510   delete r_l_Ejectile_solved_2_temp;
1511
1512   delete r_l_hadron_temp;
1513   delete[] pars;
1514
1515   return 1;
1516
1517 }