Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-09-15 08:17:22 

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

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!