From 847aa69a2e5d440e5fc8be8fee1dc6ee0ac356b3 Mon Sep 17 00:00:00 2001 From: Stephen Kay Date: Fri, 25 Aug 2023 05:15:42 -0600 Subject: [PATCH 1/2] New blank Pi0 cross section files. Modified default output to be HEPMC3. Changed README accordingly. Implemented new counters in CheckLaws fn. Now reports momentum conservation failures by component too --- Config_EIC.json | 2 +- README.md | 4 +- src/eic_evgen/Pi0_sig.cc | 0 src/eic_evgen/Pi0_sig.h | 0 src/eic_evgen/eic.cc | 4 +- src/eic_evgen/eic_pim.cc | 116 ++++++++++++++---- src/eic_evgen/eic_pim.h | 9 +- .../process_routine/DEMP_Reaction.cc | 20 +-- 8 files changed, 116 insertions(+), 39 deletions(-) create mode 100644 src/eic_evgen/Pi0_sig.cc create mode 100644 src/eic_evgen/Pi0_sig.h diff --git a/Config_EIC.json b/Config_EIC.json index b3f2125..3777c99 100644 --- a/Config_EIC.json +++ b/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 diff --git a/README.md b/README.md index fe18937..cf82515 100644 --- a/README.md +++ b/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 diff --git a/src/eic_evgen/Pi0_sig.cc b/src/eic_evgen/Pi0_sig.cc new file mode 100644 index 0000000..e69de29 diff --git a/src/eic_evgen/Pi0_sig.h b/src/eic_evgen/Pi0_sig.h new file mode 100644 index 0000000..e69de29 diff --git a/src/eic_evgen/eic.cc b/src/eic_evgen/eic.cc index d800d32..0d42170 100644 --- a/src/eic_evgen/eic.cc +++ b/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{ diff --git a/src/eic_evgen/eic_pim.cc b/src/eic_evgen/eic_pim.cc index ed51bf5..8a11e8d 100644 --- a/src/eic_evgen/eic_pim.cc +++ b/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,29 +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) { @@ -1016,19 +1053,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; } diff --git a/src/eic_evgen/eic_pim.h b/src/eic_evgen/eic_pim.h index 87e3f03..99f8fc1 100644 --- a/src/eic_evgen/eic_pim.h +++ b/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>>> SigPar; diff --git a/src/eic_evgen/process_routine/DEMP_Reaction.cc b/src/eic_evgen/process_routine/DEMP_Reaction.cc index 899fc84..595eab8 100644 --- a/src/eic_evgen/process_routine/DEMP_Reaction.cc +++ b/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; } From e3b8448e993ed08ea48f611e2882e9d5137875e8 Mon Sep 17 00:00:00 2001 From: Stephen JD Kay Date: Fri, 25 Aug 2023 05:41:38 -0600 Subject: [PATCH 2/2] Removed outdated comment. --- src/eic_evgen/eic_pim.cc | 1 - 1 file changed, 1 deletion(-) diff --git a/src/eic_evgen/eic_pim.cc b/src/eic_evgen/eic_pim.cc index 8a11e8d..14c527f 100644 --- a/src/eic_evgen/eic_pim.cc +++ b/src/eic_evgen/eic_pim.cc @@ -1038,7 +1038,6 @@ int pim::CheckLaws(TLorentzVector P_E0, TLorentzVector P_t, TLorentzVector P_e, 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());