Merge pull request #16 from sjdkay/develop

Merging in these fairly minor changes now, I'm going to be pusing and testing several major updates to my own repo that I don't want to include in the main repo at the moment.

Config_EIC.json

@@ -28,7 +28,7 @@
     "ebeam": 5,                         // Electron beam energy in GeV
     "hbeam": 100,                        // Hadron beam energy in GeV
     "hbeam_part":"Proton",           // Hadron beam particle, proton, deut or helium3, work in progress, will need to be added to shell script later on as an argument
-    "OutputType": "HEPMC3",		   // choices: LUND (Docker), Pythia6 (ECCE Fun4All) or HEPMC3 (Athena)   
+    "OutputType": "HEPMC3",		   // choices: LUND (Docker), Pythia6 (ECCE Fun4All) or HEPMC3 (ePIC)   
     "det_location": "ip6",		   // choices: ip6 for STAR, ip8 for PHENIX
     "Ee_Low": 0.5, // The minimum scattered electron energy that will be generated as a fraction of the electron beam energy
     "Ee_High": 2.5, // The maximum scattered electron energy that will be generated as a fraction of the electron beam energy

README.md

@@ -114,7 +114,7 @@ This script copies Config_EIC.json and formats a new file based upon this, DO NO
   - Arg 3 - EBeamE -> The electron beam energy, set this to whatever you want, typically, we use 5, 10 or 18 (the nominal max for the EIC)  
   - Arg 4 - HBeamE -> The hadron beam energy, again, set this to whatevr you want. Typically we use 41, 100 or 275 (41 and 275 being the nominal min/max)  
   - Arg 5 - RandomSeed -> The random seed, self explanatory. Set this however you like, the batch submission job randomly generates a random seed to feed in here  
-  - Arg 6 - OutputType -> The format of the output file, select from LUND, Pythia6 (for Fun4All) or HEPMC3 (for ATHENA), the default is Pythia6 if your choice is invalid  
+  - Arg 6 - OutputType -> The format of the output file, select from LUND, Pythia6 (for ECCE/Fun4All) or HEPMC3 (for ePIC), the default is HEPMC3 if your choice is invalid  
   - Arg 7 - InteractionPoint -> The interaction point, choose from ip6 or ip8. The default is ip6 if your choice is invalid  
   - Arg 8 - Particle -> The produced particle (meson) in the reaction, choose from omega, pi+, pi0 or K+  
   - Arg 9 - Hadron -> OPTIONAL - This only matters if you select K+ as the particle, in this case, choose from Lambda or Sigma0 here. If your choice is invalid (or you don't specify arg9), the default is Lambda  
@@ -134,7 +134,7 @@ The jobs the script creates and submits all execute the Process_EIC.csh script d
   - Arg 2 - NumEvents -> The number of events thrown for this file, set this to whatever you want to run. For reference, with the Pi+/K+ generator, 1B files takes ~1 hour  
   - Arg 3 - EBeamE -> The electron beam energy, set this to whatever you want, typically, we use 5, 10 or 18 (the nominal max for the EIC)  
   - Arg 4 - HBeamE -> The hadron beam energy, again, set this to whatevr you want. Typically we use 41, 100 or 275 (41 and 275 being the nominal min/max)  
-  - Arg 5 - OutputType -> The format of the output file, select from LUND, Pythia6 (for Fun4All) or HEPMC3 (for ATHENA), the default is Pythia6 if your choice is invalid  
+  - Arg 5 - OutputType -> The format of the output file, select from LUND, Pythia6 (for ECCE/Fun4All) or HEPMC3 (for ePIC), the default is HEPMC3 if your choice is invalid  
   - Arg 6 - InteractionPoint -> The interaction point, choose from ip6 or ip8. The default is ip6 if your choice is invalid  
   - Arg 7 - Particle -> The produced particle (meson) in the reaction, choose from omega, pi+, pi0 or K+  
   - Arg 8 - Hadron -> OPTIONAL - This only matters if you select K+ as the particle, in this case, choose from Lambda or Sigma0 here. If your choice is invalid (or you don't specify arg9), the default is Lambda  

src/eic_evgen/eic.cc

@@ -233,8 +233,8 @@ void eic(Json::Value obj) {
 	  }
 	  else{
 	    cout << "Output type not recognised!" << endl;
-	    cout << "Setting output type to Pythia6 by default!" << endl;
-	    gOutputType = "Pythia6";
+	    cout << "Setting output type to HEPMC3 by default!" << endl;
+	    gOutputType = "HEPMC3";
 	  }
 	}
 	else{

src/eic_evgen/eic_pim.cc

@@ -36,16 +36,16 @@ int fSeed;
 bool allset, print, kCalcFermi, kCalcBremss, kCalcIon, kCalcBremssEle, kCalcIonEle, kSConserve, kFSI, kMSele, kMS;
 int fWLessShell, fWLess1P9, fSDiff;
 
-//long int fNEvents, fNRecorded, fNGenerated, fWSqNeg, fNMomConserve, fNSigmaNeg, fNWeightUnphys, fNWeightReject, fLundRecorded, fNFile; 
+//long int fNEvents, fNRecorded, fNGenerated, fWSqNeg, fNSigmaNeg, fNWeightUnphys, fNWeightReject, fLundRecorded, fNFile; 
 
-unsigned long long int fNEvents, fNRecorded, fNGenerated, fWSqNeg, fNMomConserve, fNSigmaNeg, fNaN, fConserve, fNWeightUnphys, fNWeightReject, fLundRecorded, fNFile;
+unsigned long long int fNEvents, fNRecorded, fNGenerated, fWSqNeg, fNSigmaNeg, fNaN, fConserve, fNWeightUnphys, fNWeightReject, fLundRecorded, fNFile;
 
 // SJDK 03/04/23 - Added in Qsq Min/Max and W Min/Max
 double fK, fm, fElectron_Kin_Col_GeV, fElectron_Kin_Col, fRand, fLumi, fuBcm2, fPI, fDEG2RAD, fRAD2DEG, fEBeam, fPBeam, fScatElec_Theta_I, fScatElec_Theta_F, fPion_Theta_I, fPion_Theta_F, fEjectileX_Theta_I, fEjectileX_Theta_F, fScatElec_E_Hi, fScatElec_E_Lo, fPSF, fQsq_Min, fQsq_Max, fW_Min, fW_Max; 
 
 double fOmega_Theta_I, fOmega_Theta_F, fOmega_Theta_Col, fOmega_Phi_Col;
 
-double fDiff_E, conserve, ene, mom, ene_mom;     // 18/06/21 AU -> New variables to count envents passing/not passing conservation laws
+double fDiff_E, conserve, ene, mom, ene_mom, mom_px, mom_py, mom_pz, mom_pxpy, mom_pxpz, mom_pypz, mom_pxpypz;     // 18/06/21 AU -> New variables to count envents passing/not passing conservation laws
 
 double fMandSConserve, fTop_Pion_Mom, fBot_Pion_Mom, fPion_Mom_Same, fEnergyConserve, fXMomConserve, fYMomConserve, fZMomConserve, fXMomConserve_RF, fYMomConserve_RF, fZMomConserve_RF, fEnergyConserve_RF; 
 
@@ -321,14 +321,12 @@ void pim::Initilize() {
     fPion_Alpha                                 = 0;
     fPion_Beta                                  = 0;
     fNRecorded                                  = 0;
-    fLundRecorded                               = 0;
     fNGenerated                                 = 0;
     fRatio                                      = 0;
     fWLessShell                                 = 0;
     fWLess1P9                                   = 0;
     fWSqNeg                                     = 0;
     fNSigmaNeg                                  = 0;
-    fNMomConserve                               = 0;
     // SJDK 15/06/21 - Integer counters to check number returning NaN and failing conservation laws added
     fNaN                                        = 0;
     fConserve                                   = 0;
@@ -938,6 +936,13 @@ void pim::Initilize() {
     ene                                         = 0; // The number that FAIL due to failing energy conservation check ONLY
     mom                                         = 0; // The number that FAIL due to failing the momentum conservation check ONLY
     ene_mom                                     = 0; // The number that FAIL BOTH energy and momentum check
+    mom_px                                      = 0; // Fail due to px only
+    mom_py                                      = 0; // Fail due to py only
+    mom_pz                                      = 0; // Fail due to pz only
+    mom_pxpy                                    = 0; // Fail due to px and py
+    mom_pxpz                                    = 0; // Fail due to px and pz
+    mom_pypz                                    = 0; // Fail due to py and pz
+    mom_pxpypz                                  = 0; // Fail due to px, py and pz
 
 }
 
@@ -978,30 +983,61 @@ int pim::CheckLaws(TLorentzVector P_E0, TLorentzVector P_t, TLorentzVector P_e,
   double px_check =(P_t.Px() + P_E0.Px()) - (P_e.Px()+P_pim.Px()+P_pro.Px());
   double py_check =(P_t.Py() + P_E0.Py()) - (P_e.Py()+P_pim.Py()+P_pro.Py());
   double pz_check =(P_t.Pz() + P_E0.Pz()) - (P_e.Pz()+P_pim.Pz()+P_pro.Pz());
-                                                                                                 
+
   Int_t err = -1;
   if( fabs( energy_check ) < fDiff  && fabs( px_check ) < fDiff  &&  fabs( py_check ) < fDiff && fabs( pz_check ) < fDiff){  // if both momentum components and energy pass the conservation check (simultaneously)
     conserve++;
     err = 1;
   }
-
+  
   else{
-    if((fabs( px_check ) < fDiff  &&  fabs( py_check ) < fDiff && fabs( pz_check ) < fDiff) == false ){ // If momentum check fails, check if energy also failed, add counters accordingly
-      if ( (fabs( energy_check ) < fDiff) == false ){
+    if( (fabs( px_check ) > fDiff) || (fabs( py_check ) > fDiff) || (fabs( pz_check ) > fDiff) ){ // If momentum check fails, check if energy also failed, add counters accordingly  	
+
+      // Check components
+      if( fabs( px_check ) > fDiff){
+	// px failed, check py
+	if( fabs( py_check ) > fDiff){
+	  // py failed, check pz
+	  if( fabs( pz_check ) > fDiff){
+	    // pz failed, all 3 failed
+	    mom_pxpypz++; // All failed
+	  }
+	  else mom_pxpy++; // px and py failed
+	}
+	else if( fabs( py_check ) < fDiff){
+	  // py passed, check pz
+	  if( fabs( pz_check ) > fDiff){
+	    mom_pxpz++; // px and pz failed
+	  }
+	  else mom_px++; // px failed, py and pz passd
+	}
+      }
+
+      else if ( abs( px_check ) < fDiff){
+	// px passed, check py
+	if( fabs( py_check ) > fDiff){
+	  // py failed, check pz
+	  if( fabs( pz_check ) > fDiff){
+	    mom_pypz++; // py and pz failed
+	  }
+	  else mom_py++; // Just py failed
+	}
+	else mom_pz++; // Only option left (since we know one of them failed) is that pz failed
+      }
+
+      if( fabs( energy_check ) > fDiff){
 	ene_mom++; // Both failed
       }
-      else{
+      else if( abs( energy_check ) < fDiff){
 	mom++; // Only momentum failed
       }
-    }
-    else{ // If check did not pass, but it wasn't momentum, must have been energy
-      ene++; // Energy failed, add to counter
-    }
+    }    
+    else ene++; // If check did not pass, but it wasn't momentum, must have been energy, add to counter
   }
-
-  return err;
+  
+    return err;
 }
-// SJDK - 01/06/23 - This should actually be even more flexible, add an fDiff_mom too - Set momentum and energy differences separately
+
 int pim::CheckLaws(TLorentzVector P_E0, TLorentzVector P_t, TLorentzVector P_e, TLorentzVector P_pim, TLorentzVector P_pro, double fDiff_E) {
 
   double energy_check = (P_t.E() + P_E0.E()) - (P_e.E()+P_pim.E()+P_pro.E());
@@ -1016,19 +1052,50 @@ int pim::CheckLaws(TLorentzVector P_E0, TLorentzVector P_t, TLorentzVector P_e,
   }
 
   else{
-    if((fabs( px_check ) < fDiff_E  &&  fabs( py_check ) < fDiff_E && fabs( pz_check ) < fDiff_E) == false ){ // If momentum check fails, check if energy also failed, add counters accordingly
-      if ( (fabs( energy_check ) < fDiff_E) == false ){
+    if( (fabs( px_check ) > fDiff_E) || (fabs( py_check ) > fDiff_E) || (fabs( pz_check ) > fDiff_E) ){ // If momentum check fails, check if energy also failed, add counters accordingly  	
+
+      // Check components
+      if( fabs( px_check ) > fDiff_E){
+	// px failed, check py
+	if( fabs( py_check ) > fDiff_E){
+	  // py failed, check pz
+	  if( fabs( pz_check ) > fDiff_E){
+	    // pz failed, all 3 failed
+	    mom_pxpypz++; // All failed
+	  }
+	  else mom_pxpy++; // px and py failed
+	}
+	else if( fabs( py_check ) < fDiff_E){
+	  // py passed, check pz
+	  if( fabs( pz_check ) > fDiff_E){
+	    mom_pxpz++; // px and pz failed
+	  }
+	  else mom_px++; // px failed, py and pz passd
+	}
+      }
+
+      else if ( abs( px_check ) < fDiff_E){
+	// px passed, check py
+	if( fabs( py_check ) > fDiff_E){
+	  // py failed, check pz
+	  if( fabs( pz_check ) > fDiff_E){
+	    mom_pypz++; // py and pz failed
+	  }
+	  else mom_py++; // Just py failed
+	}
+	else mom_pz++; // Only option left (since we know one of them failed) is that pz failed
+      }
+
+      if( fabs( energy_check ) > fDiff_E){
 	ene_mom++; // Both failed
       }
-      else{
+      else if( abs( energy_check ) < fDiff_E){
 	mom++; // Only momentum failed
       }
-    }
-    else{ // If check did not pass, but it wasn't momentum, must have been energy
-      ene++; // Energy failed, add to counter
-    }
+    }    
+    else ene++; // If check did not pass, but it wasn't momentum, must have been energy, add to counter
   }
-
+  
   return err;
 }
 

src/eic_evgen/eic_pim.h

@@ -115,7 +115,6 @@ extern unsigned long long int fConserve;
 extern unsigned long long int fNWeightUnphys;
 extern unsigned long long int fNWeightReject;
 
-extern unsigned long long int fLundRecorded;
 extern unsigned long long int fNFile;
 
 extern double fK;
@@ -841,6 +840,14 @@ extern double ene; // The number that FAIL due to failing energy conservation ch
 extern double mom; // The number that FAIL due to failing the momentum conservation check ONLY
 extern double ene_mom; // The number that FAIL BOTH energy and momentum check
 
+extern double mom_px; // Number that fail due to px only
+extern double mom_py; // Number that fail due to py only
+extern double mom_pz; // Number that fail due to pz only
+extern double mom_pxpy; // Number that fail due to px and py
+extern double mom_pxpz; // Number that fail due to px and pz
+extern double mom_pypz; // Number that fail due to py and pz
+extern double mom_pxpypz; // Number that fail due to px, py and pz
+
 // 27/01/22 - Love Preet - Adding in vector of cross section parameters
 extern vector<vector<vector<vector<double>>>> SigPar;
 

src/eic_evgen/process_routine/DEMP_Reaction.cc

@@ -203,7 +203,6 @@ void DEMP_Reaction::Init() {
     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;
   }
 
-
 }
 
 void DEMP_Reaction::Processing_Event() {
@@ -593,11 +592,10 @@ void DEMP_Reaction::Processing_Event() {
   fEventWeight = abs(fSigma_Col * fPSF * fuBcm2 * fLumi / fNEvents);   // in Hz
 
   fNRecorded++;
-  fLundRecorded++;
   fRatio = fNRecorded / fNGenerated;
 
   if (gOutputType == "Pythia6"){
-      DEMPReact_Pythia6_Output();
+    DEMPReact_Pythia6_Output();
   }
   else if (gOutputType == "LUND"){
     Lund_Output();
@@ -762,20 +760,24 @@ void DEMP_Reaction::Detail_Output() {
   DEMPDetails << "Total events tried                                           " << setw(20) << fNGenerated   << endl;
   DEMPDetails << "Total events recorded                                        " << setw(20) << fNRecorded    << endl;
   DEMPDetails << "Number of events with wsq negative                           " << setw(20) << w_neg_ev      << endl;
-  DEMPDetails << "Number of events with " << fW_Min << " < w < " << fW_Max       << "                             " << setw(20) << w_ev << endl;
+  DEMPDetails << "Number of events with " << fW_Min << " < w < " << fW_Max       << "                           " << setw(20) << w_ev << endl;
   DEMPDetails << "Number of events with " << fQsq_Min << " < qsq < " << fQsq_Max << "                           " << setw(20) << qsq_ev << endl;
   DEMPDetails << "Number of events with Meson (X) energy NaN                   " << setw(20) << fNaN          << endl;
-  DEMPDetails << "Number of events failing conservation law check              " << setw(20) << fConserve     << endl;
   DEMPDetails << "Total events passing conservation law check with tolerance " << fDiff << setw(17) << conserve   << endl;
+  DEMPDetails << "Total events failing conservation law checks                 " << setw(20) << fConserve     << endl;
   DEMPDetails << "Total events failing energy conservation check ONLY          " << setw(20) << ene   << endl; 
   DEMPDetails << "Total events failing momentum conservation check ONLY        " << setw(20) << mom   << endl;
-  DEMPDetails << "Total events failing momentum and energy conservation checks " << setw(20) << ene_mom   << endl;
+  DEMPDetails << "Total events failing energy AND momentum conservation checks " << setw(20) << ene_mom   << endl;
+  DEMPDetails << "Total events failing px conservation law check               " << setw(20) << mom_px   << endl;
+  DEMPDetails << "Total events failing py conservation law check               " << setw(20) << mom_py   << endl;
+  DEMPDetails << "Total events failing pz conservation law check               " << setw(20) << mom_pz   << endl;
+  DEMPDetails << "Total events failing px and py conservation law checks       " << setw(20) << mom_pxpy   << endl;
+  DEMPDetails << "Total events failing px and pz conservation law checks       " << setw(20) << mom_pxpz   << endl;
+  DEMPDetails << "Total events failing py and pz conservation law checks       " << setw(20) << mom_pypz   << endl;
+  DEMPDetails << "Total events failing px, py and pz conservation law checks   " << setw(20) << mom_pxpypz   << endl;
   DEMPDetails << "Number of events with -t > 2 (K+) or -t > 1.3 (Pi+) GeV      " << setw(20) << t_ev          << endl;
   DEMPDetails << "Number of events with w less than threshold                  " << setw(20) << fWSqNeg       << endl;
-  DEMPDetails << "Number of events with mom not conserve                       " << setw(20) << fNMomConserve << endl;
   DEMPDetails << "Number of events with Sigma negative                         " << setw(20) << fNSigmaNeg    << endl;
-  DEMPDetails << "Number of lund events                                        " << setw(20) << fLundRecorded << endl;
-
   DEMPDetails << "Seed used for the Random Number Generator                    " << setw(20) << fSeed         << endl;
 
 }