From b50710f051ec44bd041219eb6d5a574e9c96b7ac Mon Sep 17 00:00:00 2001 From: Stephen Kay Date: Wed, 6 Sep 2023 07:07:47 -0600 Subject: [PATCH 01/21] Updated error flags in eic.cc --- src/eic_evgen/eic.cc | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/eic_evgen/eic.cc b/src/eic_evgen/eic.cc index 0d42170..5d8deb2 100644 --- a/src/eic_evgen/eic.cc +++ b/src/eic_evgen/eic.cc @@ -390,26 +390,26 @@ vector>>> ReadCrossSectionPar(TString particle, TSt string sigL_ParamFile, sigT_ParamFile; if (particle == "Pi+" && hadron == "Neutron"){ - cout << "Add Pi+/Neutron case here" << endl; + //cout << "Add Pi+/Neutron case here" << endl; } else if (particle == "Pi-" && hadron == "Proton"){ - cout << "Add Pi-/Proton case here" << endl; + //cout << "Add Pi-/Proton case here" << endl; } else if (particle == "K+" && hadron == "Lambda"){ - cout << "Add K+/Lambda case here" << endl; + //cout << "Add K+/Lambda case here" << endl; sigL_ParamFile = "../src/eic_evgen/CrossSection_Params/KPlusLambda_Param_sigL"; sigT_ParamFile = "../src/eic_evgen/CrossSection_Params/KPlusLambda_Param_sigT"; // Shouldn't really have a relative path, should look at setting a DEMPGen variable and doing this in a better way later } else if (particle == "K+" && hadron == "Sigma0"){ - cout << "Add K+/Sigma case here" << endl; + //cout << "Add K+/Sigma case here" << endl; sigL_ParamFile = "../src/eic_evgen/CrossSection_Params/KPlusSigma_Param_sigL"; sigT_ParamFile = "../src/eic_evgen/CrossSection_Params/KPlusSigma_Param_sigT"; } else if (particle == "Pi0"){ - cout << "Add Pi0 case here" << endl; + //cout << "Add Pi0 case here" << endl; } else{ - cout << "Throw some error" << endl; + cerr << " !!!!! " << endl << "Warning! Combination of specified ejectile and recoil particles not found!" << " !!!!! " << endl; } //.................................................................................................... From 9308ba3ef618af776adb43426f436777737d3201 Mon Sep 17 00:00:00 2001 From: Stephen Kay Date: Wed, 6 Sep 2023 07:14:03 -0600 Subject: [PATCH 02/21] Updated error messages in eic.cc --- src/eic_evgen/eic.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/eic_evgen/eic.cc b/src/eic_evgen/eic.cc index 5d8deb2..095e50f 100644 --- a/src/eic_evgen/eic.cc +++ b/src/eic_evgen/eic.cc @@ -409,7 +409,7 @@ vector>>> ReadCrossSectionPar(TString particle, TSt //cout << "Add Pi0 case here" << endl; } else{ - cerr << " !!!!! " << endl << "Warning! Combination of specified ejectile and recoil particles not found!" << " !!!!! " << endl; + cerr << " !!!!!! " << endl << "Warning! Combination of specified ejectile and recoil particles not found!" << " !!!!!! " << endl; } //.................................................................................................... From 13801b15200510b9f98f0834f9cb4562116e53b3 Mon Sep 17 00:00:00 2001 From: Stephen Kay Date: Wed, 6 Sep 2023 07:19:15 -0600 Subject: [PATCH 03/21] Else case now exits code as intended --- src/eic_evgen/eic.cc | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/eic_evgen/eic.cc b/src/eic_evgen/eic.cc index 095e50f..e8cc7e8 100644 --- a/src/eic_evgen/eic.cc +++ b/src/eic_evgen/eic.cc @@ -409,7 +409,8 @@ vector>>> ReadCrossSectionPar(TString particle, TSt //cout << "Add Pi0 case here" << endl; } else{ - cerr << " !!!!!! " << endl << "Warning! Combination of specified ejectile and recoil particles not found!" << " !!!!!! " << endl; + cerr << " !!!!! " << endl << "Warning!" << endl << "Combination of specified ejectile and recoil particles not found!" << "Cross section parameters cannot be read, check inputs!" << endl << "Warning!" << endl << " !!!!! " << endl; + exit(-1); } //.................................................................................................... From 3eede6696d7c7c503b3f7be85929e85b2dbb28c7 Mon Sep 17 00:00:00 2001 From: Stephen Kay Date: Wed, 6 Sep 2023 07:36:32 -0600 Subject: [PATCH 04/21] Added a new variable that sets the limit for the t cut based upon ejectile. This is fixed in eic.cc --- src/eic_evgen/eic.cc | 4 ++++ src/eic_evgen/eic_pim.cc | 4 ++-- src/eic_evgen/eic_pim.h | 1 + src/eic_evgen/process_routine/DEMP_Reaction.cc | 15 ++++----------- 4 files changed, 11 insertions(+), 13 deletions(-) diff --git a/src/eic_evgen/eic.cc b/src/eic_evgen/eic.cc index e8cc7e8..886273c 100644 --- a/src/eic_evgen/eic.cc +++ b/src/eic_evgen/eic.cc @@ -179,18 +179,22 @@ void eic(Json::Value obj) { if (particle == "pi+" || particle == "Pion+" || particle == "Pi+"){ fQsq_Min = 3.0; fQsq_Max = 35.0; fW_Min = 2.0; fW_Max = 10.2; + fT_Max = 1.3; } else if (particle == "pi0" || particle == "Pion0" || particle == "Pi0"){ fQsq_Min = 5.0; fQsq_Max = 1000.0; fW_Min = 2.0; fW_Max = 10.0; + fT_Max = 0.5; } else if (particle == "K+"){ fQsq_Min = 1.0; fQsq_Max = 35.0; fW_Min = 2.0; fW_Max = 10.0; + fT_Max = 2.0; } else{ fQsq_Min = 5.0; fQsq_Max = 35.0; fW_Min = 2.0; fW_Max = 10.0; + fT_Max = 1.3; } // SJDK - 01/06/21 diff --git a/src/eic_evgen/eic_pim.cc b/src/eic_evgen/eic_pim.cc index 14c527f..24cb9f7 100644 --- a/src/eic_evgen/eic_pim.cc +++ b/src/eic_evgen/eic_pim.cc @@ -40,8 +40,8 @@ int fWLessShell, fWLess1P9, fSDiff; 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; +// SJDK 03/04/23 - Added in Qsq Min/Max, W Min/Max and t max (06/09/23) +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, fT_Max; double fOmega_Theta_I, fOmega_Theta_F, fOmega_Theta_Col, fOmega_Phi_Col; diff --git a/src/eic_evgen/eic_pim.h b/src/eic_evgen/eic_pim.h index 99f8fc1..d0b77f8 100644 --- a/src/eic_evgen/eic_pim.h +++ b/src/eic_evgen/eic_pim.h @@ -143,6 +143,7 @@ extern double fQsq_Min; extern double fQsq_Max; extern double fW_Min; extern double fW_Max; +extern double fT_Max; extern double fMandSConserve; extern double fTop_Pion_Mom; diff --git a/src/eic_evgen/process_routine/DEMP_Reaction.cc b/src/eic_evgen/process_routine/DEMP_Reaction.cc index 595eab8..6ff27af 100644 --- a/src/eic_evgen/process_routine/DEMP_Reaction.cc +++ b/src/eic_evgen/process_routine/DEMP_Reaction.cc @@ -461,16 +461,9 @@ void DEMP_Reaction::Processing_Event() { } */ - // 31/01/23 SJDK - New limit on t, remove only events outside the parameterisation range, limits depend upon particle typ - if (rEjectile == "Pi+" && fT_GeV > 1.3 ) { - t_ev++; - return; - } - else if (rEjectile == "K+" && fT_GeV > 2.0) { - t_ev++; - return; - } - else if (rEjectile == "Pi0+" && fT_GeV > 0.5){ // 03/02/23 - SJDK - Not sure what range is used for pi0, assume < 0.5 for now, would be u in this case anyway? + // 31/01/23 SJDK - New limit on t, remove only events outside the parameterisation range + // 06/09/23 SJDK - fT_Max set in eic.cc depending upon ejectile type + if (fT_GeV > fT_Max ) { t_ev++; return; } @@ -775,7 +768,7 @@ void DEMP_Reaction::Detail_Output() { 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 -t > " << fT_Max << "GeV " << setw(30) << t_ev << endl; DEMPDetails << "Number of events with w less than threshold " << setw(20) << fWSqNeg << endl; DEMPDetails << "Number of events with Sigma negative " << setw(20) << fNSigmaNeg << endl; DEMPDetails << "Seed used for the Random Number Generator " << setw(20) << fSeed << endl; From cf52f84792d98733380bb59ae1729d4dcbc7343a Mon Sep 17 00:00:00 2001 From: Stephen Kay Date: Wed, 6 Sep 2023 07:40:20 -0600 Subject: [PATCH 05/21] Removed determination of w'. Not used. Not affect on output --- src/eic_evgen/process_routine/DEMP_Reaction.cc | 8 +------- 1 file changed, 1 insertion(+), 7 deletions(-) diff --git a/src/eic_evgen/process_routine/DEMP_Reaction.cc b/src/eic_evgen/process_routine/DEMP_Reaction.cc index 6ff27af..bbe9d10 100644 --- a/src/eic_evgen/process_routine/DEMP_Reaction.cc +++ b/src/eic_evgen/process_routine/DEMP_Reaction.cc @@ -351,13 +351,6 @@ void DEMP_Reaction::Processing_Event() { r_lw = r_lproton + r_lphoton; fW = r_lw.Mag(); - // ---------------------------------------------------------------------------------------------- - // Calculate w prime w' = (proton + photon - pion)^2 - // ---------------------------------------------------------------------------------------------- - - lwp = r_lprotong + r_lphotong - r_l_Ejectile_g; - fW_Prime_GeV = lwp.Mag(); - fsini = r_lelectron + r_lproton; fsfin = r_lscatelec + r_l_Ejectile + l_Recoil; @@ -747,6 +740,7 @@ Double_t DEMP_Reaction::Get_Total_Cross_Section() { /*--------------------------------------------------*/ /// Output generator detail +// 06/06/23 SJDK - Formatting of these is all messed up annoyingly, would be nice to get the final numbers to align. They don't currently void DEMP_Reaction::Detail_Output() { From ae1138b9960a660d2f2556cc4e9bc42ad26f98df Mon Sep 17 00:00:00 2001 From: Stephen Kay Date: Wed, 6 Sep 2023 07:43:52 -0600 Subject: [PATCH 06/21] Removed mandelstam s calculation and conservation check. Not used. --- src/eic_evgen/process_routine/DEMP_Reaction.cc | 8 -------- 1 file changed, 8 deletions(-) diff --git a/src/eic_evgen/process_routine/DEMP_Reaction.cc b/src/eic_evgen/process_routine/DEMP_Reaction.cc index bbe9d10..aa31b92 100644 --- a/src/eic_evgen/process_routine/DEMP_Reaction.cc +++ b/src/eic_evgen/process_routine/DEMP_Reaction.cc @@ -351,14 +351,6 @@ void DEMP_Reaction::Processing_Event() { r_lw = r_lproton + r_lphoton; fW = r_lw.Mag(); - fsini = r_lelectron + r_lproton; - fsfin = r_lscatelec + r_l_Ejectile + l_Recoil; - - fsinig = fsini * fm; - fsfing = fsfin * fm; - // SJDK 15/06/21 - Mandlestam S conservation check - doesn't actually seem to be utilised? - fMandSConserve = std::abs( fsinig.Mag() - fsfing.Mag() ); - // SJDK 15/06/21 - Added integer counters for conservation law check and for NaN check 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 fNaN++; From d0b568eda3d448f3ad96753709dc6d0efc09345b Mon Sep 17 00:00:00 2001 From: Stephen Kay Date: Wed, 6 Sep 2023 07:46:41 -0600 Subject: [PATCH 07/21] Removed kSConserve true/false setting. Not used, did not alter output at all --- src/eic_evgen/process_routine/DEMP_Reaction.cc | 4 ---- 1 file changed, 4 deletions(-) diff --git a/src/eic_evgen/process_routine/DEMP_Reaction.cc b/src/eic_evgen/process_routine/DEMP_Reaction.cc index aa31b92..df2ae35 100644 --- a/src/eic_evgen/process_routine/DEMP_Reaction.cc +++ b/src/eic_evgen/process_routine/DEMP_Reaction.cc @@ -356,10 +356,6 @@ void DEMP_Reaction::Processing_Event() { fNaN++; return; } - kSConserve = false; - if( std::abs( fsinig.Mag() - fsfing.Mag() ) < fDiff ) { - kSConserve = true; - } ///*--------------------------------------------------*/ //-> 10/05/23 - Love added a slimmed down, simpler to read version of the CheckLaws fn From 471ddccadf5b146dc5a862ba48f352eb14ed04c6 Mon Sep 17 00:00:00 2001 From: Stephen Kay Date: Wed, 6 Sep 2023 07:57:08 -0600 Subject: [PATCH 08/21] Commented out fermi momentum parts of DEMP_Reaction.cc for now. These should not be enabled for e/p collisions. This does alter then number of events. See more w/o fermi momentum enabled --- .../process_routine/DEMP_Reaction.cc | 34 +++++++++---------- 1 file changed, 17 insertions(+), 17 deletions(-) diff --git a/src/eic_evgen/process_routine/DEMP_Reaction.cc b/src/eic_evgen/process_routine/DEMP_Reaction.cc index df2ae35..75d35d4 100644 --- a/src/eic_evgen/process_routine/DEMP_Reaction.cc +++ b/src/eic_evgen/process_routine/DEMP_Reaction.cc @@ -108,7 +108,7 @@ void DEMP_Reaction::Init() { // cout << rNEvents << " " << fNEvents << endl; - rFermiMomentum = pd->fermiMomentum(); + //rFermiMomentum = pd->fermiMomentum(); // ---------------------------------------------------- // Proton in collider (lab) frame @@ -128,7 +128,7 @@ void DEMP_Reaction::Init() { // ---------------------------------------------------- // Electron in collider (lab) frame - cout << "Fermi momentum: " << rFermiMomentum << endl; + //cout << "Fermi momentum: " << rFermiMomentum << endl; r_lelectron = GetElectronVector_lab(); r_lelectrong = r_lelectron * fm; @@ -211,9 +211,9 @@ void DEMP_Reaction::Processing_Event() { // Considering Fermi momentum for the proton // ---------------------------------------------------- // SJDK - 31/01/23 - This doesn't seem to do anything? - if( kCalcFermi ) { - Consider_Proton_Fermi_Momentum(); - } + // if( kCalcFermi ) { + // Consider_Proton_Fermi_Momentum(); + // } // ---------------------------------------------------- // Boost vector from collider (lab) frame to protons rest frame (Fix target) @@ -629,24 +629,24 @@ TLorentzVector DEMP_Reaction::GetProtonVector_lab() { // Proton in collider (lab) frame // ---------------------------------------------------- -void DEMP_Reaction::Consider_Proton_Fermi_Momentum() { +// void DEMP_Reaction::Consider_Proton_Fermi_Momentum() { - fProton_Mom_Col = fProton_Mom_Col + rFermiMomentum; - fProton_Theta_Col = acos( fRandom->Uniform( cos(0.0) , cos(fPi) ) ); - fProton_Phi_Col = fRandom->Uniform( 0 , 360 ); +// fProton_Mom_Col = fProton_Mom_Col + rFermiMomentum; +// fProton_Theta_Col = acos( fRandom->Uniform( cos(0.0) , cos(fPi) ) ); +// fProton_Phi_Col = fRandom->Uniform( 0 , 360 ); - double px, py, pz, e; +// double px, py, pz, e; - px = fProton_Mom_Col * sin(fProton_Theta_Col) * cos(fProton_Phi_Col); - py = fProton_Mom_Col * sin(fProton_Theta_Col) * sin(fProton_Phi_Col); - pz = fProton_Mom_Col * cos(fProton_Theta_Col); - e = sqrt( pow( fProton_Mom_Col , 2 ) + pow( fProton_Mass , 2 ) ); +// px = fProton_Mom_Col * sin(fProton_Theta_Col) * cos(fProton_Phi_Col); +// py = fProton_Mom_Col * sin(fProton_Theta_Col) * sin(fProton_Phi_Col); +// pz = fProton_Mom_Col * cos(fProton_Theta_Col); +// e = sqrt( pow( fProton_Mom_Col , 2 ) + pow( fProton_Mass , 2 ) ); - r_lproton.SetPxPyPzE(px,py,pz,e); +// r_lproton.SetPxPyPzE(px,py,pz,e); - r_lprotong = r_lproton*fm; +// r_lprotong = r_lproton*fm; -} +// } // ---------------------------------------------------- // Electron in collider (lab) frame From 6d3da0c911795e279b825e0012135cb05d12ed89 Mon Sep 17 00:00:00 2001 From: Stephen Kay Date: Wed, 6 Sep 2023 08:50:00 -0600 Subject: [PATCH 09/21] Reordered details output, clarified comments, added a check condition which prints out whether total tried = total recorded + total cut --- .../process_routine/DEMP_Reaction.cc | 90 +++++++++++-------- 1 file changed, 51 insertions(+), 39 deletions(-) diff --git a/src/eic_evgen/process_routine/DEMP_Reaction.cc b/src/eic_evgen/process_routine/DEMP_Reaction.cc index 75d35d4..636963c 100644 --- a/src/eic_evgen/process_routine/DEMP_Reaction.cc +++ b/src/eic_evgen/process_routine/DEMP_Reaction.cc @@ -108,7 +108,7 @@ void DEMP_Reaction::Init() { // cout << rNEvents << " " << fNEvents << endl; - //rFermiMomentum = pd->fermiMomentum(); + rFermiMomentum = pd->fermiMomentum(); // ---------------------------------------------------- // Proton in collider (lab) frame @@ -128,7 +128,7 @@ void DEMP_Reaction::Init() { // ---------------------------------------------------- // Electron in collider (lab) frame - //cout << "Fermi momentum: " << rFermiMomentum << endl; + cout << "Fermi momentum: " << rFermiMomentum << endl; r_lelectron = GetElectronVector_lab(); r_lelectrong = r_lelectron * fm; @@ -211,9 +211,9 @@ void DEMP_Reaction::Processing_Event() { // Considering Fermi momentum for the proton // ---------------------------------------------------- // SJDK - 31/01/23 - This doesn't seem to do anything? - // if( kCalcFermi ) { - // Consider_Proton_Fermi_Momentum(); - // } + if( kCalcFermi ) { + Consider_Proton_Fermi_Momentum(); + } // ---------------------------------------------------- // Boost vector from collider (lab) frame to protons rest frame (Fix target) @@ -629,24 +629,24 @@ TLorentzVector DEMP_Reaction::GetProtonVector_lab() { // Proton in collider (lab) frame // ---------------------------------------------------- -// void DEMP_Reaction::Consider_Proton_Fermi_Momentum() { +void DEMP_Reaction::Consider_Proton_Fermi_Momentum() { -// fProton_Mom_Col = fProton_Mom_Col + rFermiMomentum; -// fProton_Theta_Col = acos( fRandom->Uniform( cos(0.0) , cos(fPi) ) ); -// fProton_Phi_Col = fRandom->Uniform( 0 , 360 ); + fProton_Mom_Col = fProton_Mom_Col + rFermiMomentum; + fProton_Theta_Col = acos( fRandom->Uniform( cos(0.0) , cos(fPi) ) ); + fProton_Phi_Col = fRandom->Uniform( 0 , 360 ); -// double px, py, pz, e; + double px, py, pz, e; -// px = fProton_Mom_Col * sin(fProton_Theta_Col) * cos(fProton_Phi_Col); -// py = fProton_Mom_Col * sin(fProton_Theta_Col) * sin(fProton_Phi_Col); -// pz = fProton_Mom_Col * cos(fProton_Theta_Col); -// e = sqrt( pow( fProton_Mom_Col , 2 ) + pow( fProton_Mass , 2 ) ); + px = fProton_Mom_Col * sin(fProton_Theta_Col) * cos(fProton_Phi_Col); + py = fProton_Mom_Col * sin(fProton_Theta_Col) * sin(fProton_Phi_Col); + pz = fProton_Mom_Col * cos(fProton_Theta_Col); + e = sqrt( pow( fProton_Mom_Col , 2 ) + pow( fProton_Mass , 2 ) ); -// r_lproton.SetPxPyPzE(px,py,pz,e); + r_lproton.SetPxPyPzE(px,py,pz,e); -// r_lprotong = r_lproton*fm; + r_lprotong = r_lproton*fm; -// } +} // ---------------------------------------------------- // Electron in collider (lab) frame @@ -728,33 +728,45 @@ Double_t DEMP_Reaction::Get_Total_Cross_Section() { /*--------------------------------------------------*/ /// Output generator detail -// 06/06/23 SJDK - Formatting of these is all messed up annoyingly, would be nice to get the final numbers to align. They don't currently - +// 06/09/23 SJDK - Formatting of these is all messed up annoyingly, would be nice to get the final numbers to align. They don't currently. +// Cuts are now ordered as they are applied in the generator void DEMP_Reaction::Detail_Output() { + DEMPDetails << "Seed used for the Random Number Generator " << setw(20) << fSeed << endl; + DEMPDetails << endl; DEMPDetails << "Total events tried " << setw(20) << fNGenerated << endl; + DEMPDetails << "Total events cut " << setw(20) << (qsq_ev + w_ev + w_neg_ev + fNaN + fConserve + t_ev + fNSigmaNeg) << 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 " << fQsq_Min << " < qsq < " << fQsq_Max << " " << setw(20) << qsq_ev << endl; - DEMPDetails << "Number of events with Meson (X) energy NaN " << setw(20) << fNaN << 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 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 > " << fT_Max << "GeV " << setw(30) << t_ev << endl; - DEMPDetails << "Number of events with w less than threshold " << setw(20) << fWSqNeg << endl; - DEMPDetails << "Number of events with Sigma negative " << setw(20) << fNSigmaNeg << endl; - DEMPDetails << "Seed used for the Random Number Generator " << setw(20) << fSeed << endl; + if (fNGenerated != (qsq_ev + w_ev + w_neg_ev + fNaN + fConserve + t_ev + fNSigmaNeg + fNRecorded)){ + DEMPDetails << "Total events cut + recorded = events tried? " << setw(20) << "NO! ERROR!" << endl; + } + else{ + DEMPDetails << "Total events cut + recorded = events tried? " << setw(22) << "Yes! :)" << endl; + } + + DEMPDetails << endl << "Cut details -" << endl; + DEMPDetails << "Events cut due to qsq < " << fQsq_Min << " or qsq > "<< fQsq_Max << " " << setw(20) << qsq_ev << endl; + DEMPDetails << "Events cut due to negative Wsq value " << setw(20) << w_neg_ev << endl; + DEMPDetails << "Events cut due to W < " << fW_Min << " or W > " << fW_Max << " " << setw(20) << w_ev << endl; + DEMPDetails << "Events cut due to ejectile (X) energy NaN " << setw(20) << fNaN << endl; + DEMPDetails << "Events cut due to conservation law check failure " << setw(20) << fConserve << endl; + DEMPDetails << "Events cut due to -t > " << fT_Max << "GeV " << setw(30) << t_ev << endl; + DEMPDetails << "Events cut due to -ve cross section value " << setw(20) << fNSigmaNeg << endl; + + DEMPDetails << endl << "Conservation law checks details -" << endl; + DEMPDetails << "Total events PASSING conservation law check with tolerance " << fDiff << setw(17) << conserve << endl; + DEMPDetails << "Events cut due to energy conservation check ONLY " << setw(20) << ene << endl; + DEMPDetails << "Events cut due to momentum conservation check ONLY " << setw(20) << mom << endl; + DEMPDetails << "Events cut due to energy AND momentum conservation checks " << setw(20) << ene_mom << endl; + DEMPDetails << "Events cut due to px conservation law check " << setw(20) << mom_px << endl; + DEMPDetails << "Events cut due to py conservation law check " << setw(20) << mom_py << endl; + DEMPDetails << "Events cut due to pz conservation law check " << setw(20) << mom_pz << endl; + DEMPDetails << "Events cut due to px and py conservation law checks " << setw(20) << mom_pxpy << endl; + DEMPDetails << "Events cut due to px and pz conservation law checks " << setw(20) << mom_pxpz << endl; + DEMPDetails << "Events cut due to py and pz conservation law checks " << setw(20) << mom_pypz << endl; + DEMPDetails << "Events cut due to px, py and pz conservation law checks " << setw(20) << mom_pxpypz << endl; + } ////*-------------------------------------------------- From 8f42989c1e7c867881acb09864784d3b912361a2 Mon Sep 17 00:00:00 2001 From: Stephen Kay Date: Wed, 6 Sep 2023 09:38:03 -0600 Subject: [PATCH 10/21] Added in pieces of solve function, commented for now. Modified config .json file. Now has a string to specify calculation method which is checked and a boolean assigned in eic.cc. Default is Analytical if not specified or entered incorrectly. Need to get other variables for solve fn in and test actual routine --- Config_EIC.json | 1 + src/eic_evgen/eic.cc | 17 +- src/eic_evgen/eic_pim.cc | 1 + src/eic_evgen/eic_pim.h | 1 + .../process_routine/DEMP_Reaction.cc | 160 ++++++++++++------ 5 files changed, 124 insertions(+), 56 deletions(-) diff --git a/Config_EIC.json b/Config_EIC.json index 3777c99..bd33599 100644 --- a/Config_EIC.json +++ b/Config_EIC.json @@ -30,6 +30,7 @@ "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 (ePIC) "det_location": "ip6", // choices: ip6 for STAR, ip8 for PHENIX + "calc_method": "Analytical", // choices: Analytical or Solve, affects how the ejectile/recoil hadron 4-vectors are calculated. Default is analytical "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 "e_Theta_Low": 60.0, // The minimum scattered electron angle (theta) that will be generated in degrees diff --git a/src/eic_evgen/eic.cc b/src/eic_evgen/eic.cc index 886273c..4a5cd4d 100644 --- a/src/eic_evgen/eic.cc +++ b/src/eic_evgen/eic.cc @@ -316,9 +316,22 @@ void eic(Json::Value obj) { fEjectileX_Theta_F = 60.0 * fDEG2RAD; cout << "Max ejectile X theta not specified in .json file, defaulting to 60 degrees." << endl; } - - SigPar = ReadCrossSectionPar(particle, hadron); + // 06/09/23 - SJDK - Added string to check method chosen + TString CalcMethod = obj["calc_method"].asString(); + if(CalcMethod == "Analytical"){ + UseSolve = false; + } + else if (CalcMethod == "Solve"){ + UseSolve = true; + } + else{ + cout << "! WARNING !" << endl << "! WARNING !- Calculation method not specified or not recognised, defaulting to Analytical - ! WARNING!" << endl << "! WARNING !" << endl; + UseSolve = false; + } + + SigPar = ReadCrossSectionPar(particle, hadron); + if(particle != "pi0"){ // Default case now Reaction* r1 = new Reaction(particle, hadron); r1->process_reaction(); diff --git a/src/eic_evgen/eic_pim.cc b/src/eic_evgen/eic_pim.cc index 24cb9f7..724b146 100644 --- a/src/eic_evgen/eic_pim.cc +++ b/src/eic_evgen/eic_pim.cc @@ -24,6 +24,7 @@ TTree *t1; int gKinematics_type; bool gPi0_decay; +bool UseSolve; string gDet_location; string gOutputType; // SJDK 12/01/22 - Added output type as a variable you can specify in the .json file string gBeamPart; // SJDK 12/01/22 - Added output type as a variable you can specify in the .json file diff --git a/src/eic_evgen/eic_pim.h b/src/eic_evgen/eic_pim.h index d0b77f8..3198c79 100644 --- a/src/eic_evgen/eic_pim.h +++ b/src/eic_evgen/eic_pim.h @@ -50,6 +50,7 @@ extern TString gfile_name; extern TString gParticle; extern TString gHadron; extern bool gPi0_decay; +extern bool UseSolve; extern std::string gDet_location; extern std::string gOutputType; extern std::string gBeamPart; diff --git a/src/eic_evgen/process_routine/DEMP_Reaction.cc b/src/eic_evgen/process_routine/DEMP_Reaction.cc index 636963c..8879234 100644 --- a/src/eic_evgen/process_routine/DEMP_Reaction.cc +++ b/src/eic_evgen/process_routine/DEMP_Reaction.cc @@ -93,13 +93,12 @@ void DEMP_Reaction::Init() { dRootTree->Branch("EventWeight", &fEventWeight, "fEventWeight/D"); - /*--------------------------------------------------*/ qsq_ev = 0, t_ev = 0, w_neg_ev = 0, w_ev = 0; rNEvents = fNEvents; rNEvent_itt = 0; - + // 02/06/21 SJDK // Set these values once the beam energies are read in fPSF = ( fEBeam * ( fScatElec_E_Hi - fScatElec_E_Lo ) *( sin( fScatElec_Theta_F ) - sin( fScatElec_Theta_I ) ) * 2 * fPI *( sin( fEjectileX_Theta_F ) - sin( fEjectileX_Theta_I ) ) * 2 * fPI ); @@ -108,7 +107,7 @@ void DEMP_Reaction::Init() { // cout << rNEvents << " " << fNEvents << endl; - rFermiMomentum = pd->fermiMomentum(); + //rFermiMomentum = pd->fermiMomentum(); // ---------------------------------------------------- // Proton in collider (lab) frame @@ -128,7 +127,8 @@ void DEMP_Reaction::Init() { // ---------------------------------------------------- // Electron in collider (lab) frame - cout << "Fermi momentum: " << rFermiMomentum << endl; + // 06/09/23 - SJDK - Commenting out for now, should be disabled for e/p collisions + //cout << "Fermi momentum: " << rFermiMomentum << endl; r_lelectron = GetElectronVector_lab(); r_lelectrong = r_lelectron * fm; @@ -177,7 +177,12 @@ void DEMP_Reaction::Init() { cout << "Produced particle in exclusive production: " << rEjectile << "; with mass: " << f_Ejectile_Mass << " MeV "<< endl; cout << fEBeam << " GeV electrons on " << fPBeam << " GeV ions" << endl; - + if(UseSolve == true){ + cout << rEjectile << " and " << rEjectile_scat_hadron << " 4-vectors calculated using Solve function" << endl; + } + else if(UseSolve == false){ + cout << rEjectile << " and " << rEjectile_scat_hadron << " 4-vectors calculated using analytical solution" << endl; + } // 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) // See slide 11 in https://indico.cern.ch/event/1072579/contributions/4796856/attachments/2456676/4210776/CAP-EIC-June-7-2022-Seryi-r2.pdf // If available in the future, this could be replaced by some fixed function @@ -202,7 +207,48 @@ void DEMP_Reaction::Init() { else{ 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; } - + + // /*--------------------------------------------------*/ + // // SJDK 03/04/22 - New set of initialisation stuff for the solve function from Ishan and Bill + + // CoinToss = new TRandom3(); + + // // F = new TF1("F", + // // "[6]-sqrt([7]**2+x**2)-sqrt([8]**2+([3]-[0]*x)**2+([4]-[1]*x)**2+([5]-[2]*x)**2)", + // // 0, 12000); + + // F = new TF1("F", + // "[6]-sqrt([7]**2+x**2)-sqrt([8]**2+([3]-[0]*x)**2+([4]-[1]*x)**2+([5]-[2]*x)**2)", + // 0, r_lproton.E()); + + // char AngleGenName[100] = "AngleGen"; + // double dummy[2] = {0,1}; + // // Changed the theta range here to match the one actually provided in the input .json file, these are already converted to radians in eic.cc (see ~ line 293) + // // double ThetaRange[2] = {fX_Theta_I, fX_Theta_F}; + + // f_Ejectile_Theta_I = fEjectileX_Theta_I; + // f_Ejectile_Theta_F = fEjectileX_Theta_F; + + // double ThetaRange[2] = {f_Ejectile_Theta_I, f_Ejectile_Theta_F}; + // double PhiRange[2] = {0, 360*TMath::DegToRad()}; + + // AngleGen = new CustomRand(AngleGenName, dummy, + // ThetaRange, PhiRange); + + // UnitVect = new TVector3(0,0,1); + + // ///*--------------------------------------------------*/ + // // Produced hadron and recoilded hadron from the solve function + + // VertBeamElec = new TLorentzVector(); + // VertScatElec = new TLorentzVector(); + + // Initial = new TLorentzVector(); + // Target = new TLorentzVector(); + // Photon = new TLorentzVector(); + // Interaction = new TLorentzVector(); + // Final = new TLorentzVector(); + } void DEMP_Reaction::Processing_Event() { @@ -210,10 +256,10 @@ void DEMP_Reaction::Processing_Event() { // ---------------------------------------------------- // Considering Fermi momentum for the proton // ---------------------------------------------------- - // SJDK - 31/01/23 - This doesn't seem to do anything? - if( kCalcFermi ) { - Consider_Proton_Fermi_Momentum(); - } + // SJDK - 06/06/23 - Commenting out for now, this increases the number of recorded events - Should have no fermi momentum for e/p collisions + // if( kCalcFermi ) { + // Consider_Proton_Fermi_Momentum(); + // } // ---------------------------------------------------- // Boost vector from collider (lab) frame to protons rest frame (Fix target) @@ -283,62 +329,68 @@ void DEMP_Reaction::Processing_Event() { return; } - // --------------------------------------------------------- - // Pion momentum in collider frame, analytic solution starts - // --------------------------------------------------------- + // SJDK - 06/09/23 - Check UseSolved boolean, process through relevant loop + if (UseSolve == false){ + // --------------------------------------------------------- + // Pion momentum in collider frame, analytic solution starts + // --------------------------------------------------------- - double fupx = sin( f_Ejectile_Theta_Col ) * cos( f_Ejectile_Phi_Col ); - double fupy = sin( f_Ejectile_Theta_Col ) * sin( f_Ejectile_Phi_Col ); - double fupz = cos( f_Ejectile_Theta_Col ); + double fupx = sin( f_Ejectile_Theta_Col ) * cos( f_Ejectile_Phi_Col ); + double fupy = sin( f_Ejectile_Theta_Col ) * sin( f_Ejectile_Phi_Col ); + double fupz = cos( f_Ejectile_Theta_Col ); - double fuqx = sin( r_lphoton.Theta() ) * cos( r_lphoton.Phi() ); - double fuqy = sin( r_lphoton.Theta() ) * sin( r_lphoton.Phi() ); - double fuqz = cos( r_lphoton.Theta() ); + double fuqx = sin( r_lphoton.Theta() ) * cos( r_lphoton.Phi() ); + double fuqy = sin( r_lphoton.Theta() ) * sin( r_lphoton.Phi() ); + double fuqz = cos( r_lphoton.Theta() ); - double fa = -(r_lphoton.Vect()).Mag() * ( fupx * fuqx + fupy * fuqy + fupz * fuqz ); - double fb = pow ( (r_lphoton.Vect()).Mag() , 2 ); - double fc = r_lphoton.E() + fProton_Mass; + double fa = -(r_lphoton.Vect()).Mag() * ( fupx * fuqx + fupy * fuqy + fupz * fuqz ); + double fb = pow ( (r_lphoton.Vect()).Mag() , 2 ); + double fc = r_lphoton.E() + fProton_Mass; - fa = ( fa - std::abs( (r_lproton.Vect()).Mag() ) * ( ( ( r_lproton.X() / (r_lproton.Vect()).Mag() ) * fupx ) + - ( ( r_lproton.Y() / (r_lproton.Vect()).Mag() ) * fupy ) + - ( ( r_lproton.Z() / (r_lproton.Vect()).Mag() ) * fupz ) ) ); + fa = ( fa - std::abs( (r_lproton.Vect()).Mag() ) * ( ( ( r_lproton.X() / (r_lproton.Vect()).Mag() ) * fupx ) + + ( ( r_lproton.Y() / (r_lproton.Vect()).Mag() ) * fupy ) + + ( ( r_lproton.Z() / (r_lproton.Vect()).Mag() ) * fupz ) ) ); - double factor = ( pow( (r_lproton.Vect()).Mag() , 2 ) + 2.0 * (r_lphoton.Vect()).Mag() * (r_lproton.Vect()).Mag() * - ( ( ( r_lproton.X() / (r_lproton.Vect()).Mag() ) * fuqx ) + - ( ( r_lproton.Y() / (r_lproton.Vect()).Mag() ) * fuqy ) + - ( ( r_lproton.Z() / (r_lproton.Vect()).Mag() ) * fuqz ) ) ); + double factor = ( pow( (r_lproton.Vect()).Mag() , 2 ) + 2.0 * (r_lphoton.Vect()).Mag() * (r_lproton.Vect()).Mag() * + ( ( ( r_lproton.X() / (r_lproton.Vect()).Mag() ) * fuqx ) + + ( ( r_lproton.Y() / (r_lproton.Vect()).Mag() ) * fuqy ) + + ( ( r_lproton.Z() / (r_lproton.Vect()).Mag() ) * fuqz ) ) ); - fb = fb + factor; - fc = r_lphoton.E() + r_lproton.E(); + fb = fb + factor; + fc = r_lphoton.E() + r_lproton.E(); - double ft = fc * fc - fb + f_Ejectile_Mass * f_Ejectile_Mass - f_Recoil_Mass * f_Recoil_Mass; + double ft = fc * fc - fb + f_Ejectile_Mass * f_Ejectile_Mass - f_Recoil_Mass * f_Recoil_Mass; - double fQA = 4.0 * ( fa * fa - fc * fc ); - double fQB = 4.0 * fc * ft; + double fQA = 4.0 * ( fa * fa - fc * fc ); + double fQB = 4.0 * fc * ft; - double fQC = -4.0 * fa * fa * f_Ejectile_Mass * f_Ejectile_Mass - ft * ft; + double fQC = -4.0 * fa * fa * f_Ejectile_Mass * f_Ejectile_Mass - ft * ft; - fradical = fQB * fQB - 4.0 * fQA * fQC; + fradical = fQB * fQB - 4.0 * fQA * fQC; - fepi1 = ( -fQB - sqrt( fradical ) ) / ( 2.0 * fQA ); - fepi2 = ( -fQB + sqrt( fradical ) ) / ( 2.0 * fQA ); - - ///--------------------------------------------------------- - /// Particle X momentum in collider frame, analytic solution - /// And obtain recoiled proton in collider (lab) frame - ///--------------------------------------------------------- - - r_l_Ejectile.SetPxPyPzE( (sqrt( pow( fepi1 , 2) - pow(f_Ejectile_Mass , 2) ) ) * sin(f_Ejectile_Theta_Col) * cos(f_Ejectile_Phi_Col), - ( sqrt( pow( fepi1 , 2) - pow(f_Ejectile_Mass , 2) ) ) * sin(f_Ejectile_Theta_Col) * sin(f_Ejectile_Phi_Col), - ( sqrt( pow( fepi1 , 2) - pow(f_Ejectile_Mass , 2) ) ) * cos(f_Ejectile_Theta_Col), - fepi1 ); + fepi1 = ( -fQB - sqrt( fradical ) ) / ( 2.0 * fQA ); + fepi2 = ( -fQB + sqrt( fradical ) ) / ( 2.0 * fQA ); + + ///--------------------------------------------------------- + /// Particle X momentum in collider frame, analytic solution + /// And obtain recoiled proton in collider (lab) frame + ///--------------------------------------------------------- + + r_l_Ejectile.SetPxPyPzE( (sqrt( pow( fepi1 , 2) - pow(f_Ejectile_Mass , 2) ) ) * sin(f_Ejectile_Theta_Col) * cos(f_Ejectile_Phi_Col), + ( sqrt( pow( fepi1 , 2) - pow(f_Ejectile_Mass , 2) ) ) * sin(f_Ejectile_Theta_Col) * sin(f_Ejectile_Phi_Col), + ( sqrt( pow( fepi1 , 2) - pow(f_Ejectile_Mass , 2) ) ) * cos(f_Ejectile_Theta_Col), + fepi1 ); - l_Recoil.SetPxPyPzE( ( r_lproton + r_lelectron - r_lscatelec - r_l_Ejectile).X(), - ( r_lproton + r_lelectron - r_lscatelec - r_l_Ejectile ).Y(), - ( r_lproton + r_lelectron - r_lscatelec - r_l_Ejectile ).Z(), - sqrt( pow( ( ( ( r_lproton + r_lelectron - r_lscatelec - r_l_Ejectile ).Vect() ).Mag()),2) + - pow( f_Recoil_Mass , 2) ) ); - + l_Recoil.SetPxPyPzE( ( r_lproton + r_lelectron - r_lscatelec - r_l_Ejectile).X(), + ( r_lproton + r_lelectron - r_lscatelec - r_l_Ejectile ).Y(), + ( r_lproton + r_lelectron - r_lscatelec - r_l_Ejectile ).Z(), + sqrt( pow( ( ( ( r_lproton + r_lelectron - r_lscatelec - r_l_Ejectile ).Vect() ).Mag()),2) + + pow( f_Recoil_Mass , 2) ) ); + } + else if (UseSolve == true){ + r_l_Ejectile.SetPxPyPzE(0,0,0,0); + l_Recoil_g.SetPxPyPzE(0,0,0,0); + } ///-------------------------------------------------- r_l_Ejectile_g = r_l_Ejectile * fm; From a0be5cb55e9e679d042ee1d1d3e06edca1b6170f Mon Sep 17 00:00:00 2001 From: Stephen Kay Date: Fri, 8 Sep 2023 07:34:42 -0600 Subject: [PATCH 11/21] Solve function merged back in after removal of particle class usage. Now selected via a flag in the .json config flag. Output .txt file also differs depending upon calculation choice --- src/eic_evgen/eic.cc | 14 +- src/eic_evgen/eic_pim.cc | 10 +- src/eic_evgen/eic_pim.h | 4 + .../process_routine/DEMP_Reaction.cc | 299 ++++++++++++++---- src/eic_evgen/reaction_routine.cc | 9 +- src/eic_evgen/reaction_routine.h | 43 +-- 6 files changed, 277 insertions(+), 102 deletions(-) diff --git a/src/eic_evgen/eic.cc b/src/eic_evgen/eic.cc index 4a5cd4d..5c2b57f 100644 --- a/src/eic_evgen/eic.cc +++ b/src/eic_evgen/eic.cc @@ -233,17 +233,17 @@ void eic(Json::Value obj) { cout << "Using LUND output format" << endl; } else if (gOutputType == "HEPMC3"){ - cout << "Using HEPMC3 output format for EPIC" << endl; + cout << "Using HEPMC3 output format for ePIC" << endl; } else{ cout << "Output type not recognised!" << endl; - cout << "Setting output type to HEPMC3 by default!" << endl; + cout << "Setting output type to HEPMC3 (ePIC) by default!" << endl; gOutputType = "HEPMC3"; } } else{ cout << "Output type not specified in .json file!" << endl; - cout << "Setting output type to HEPMC3 by default!" << endl; + cout << "Setting output type to HEPMC3 (ePIC) by default!" << endl; gOutputType = "HEPMC3"; } ///*--------------------------------------------------*/ @@ -407,23 +407,21 @@ vector>>> ReadCrossSectionPar(TString particle, TSt string sigL_ParamFile, sigT_ParamFile; if (particle == "Pi+" && hadron == "Neutron"){ - //cout << "Add Pi+/Neutron case here" << endl; + // When pion model parameterised in some way, add Pi+/Neutron case here - } else if (particle == "Pi-" && hadron == "Proton"){ - //cout << "Add Pi-/Proton case here" << endl; + // When pion model parameterised in some way, add Pi-/Proton case here - } else if (particle == "K+" && hadron == "Lambda"){ - //cout << "Add K+/Lambda case here" << endl; sigL_ParamFile = "../src/eic_evgen/CrossSection_Params/KPlusLambda_Param_sigL"; sigT_ParamFile = "../src/eic_evgen/CrossSection_Params/KPlusLambda_Param_sigT"; // Shouldn't really have a relative path, should look at setting a DEMPGen variable and doing this in a better way later } else if (particle == "K+" && hadron == "Sigma0"){ - //cout << "Add K+/Sigma case here" << endl; sigL_ParamFile = "../src/eic_evgen/CrossSection_Params/KPlusSigma_Param_sigL"; sigT_ParamFile = "../src/eic_evgen/CrossSection_Params/KPlusSigma_Param_sigT"; } else if (particle == "Pi0"){ - //cout << "Add Pi0 case here" << endl; + // When pi0 model paramterised, add it here } else{ cerr << " !!!!! " << endl << "Warning!" << endl << "Combination of specified ejectile and recoil particles not found!" << "Cross section parameters cannot be read, check inputs!" << endl << "Warning!" << endl << " !!!!! " << endl; diff --git a/src/eic_evgen/eic_pim.cc b/src/eic_evgen/eic_pim.cc index 724b146..671cda5 100644 --- a/src/eic_evgen/eic_pim.cc +++ b/src/eic_evgen/eic_pim.cc @@ -39,14 +39,14 @@ int fWLessShell, fWLess1P9, fSDiff; //long int fNEvents, fNRecorded, fNGenerated, fWSqNeg, fNSigmaNeg, fNWeightUnphys, fNWeightReject, fLundRecorded, fNFile; -unsigned long long int fNEvents, fNRecorded, fNGenerated, fWSqNeg, fNSigmaNeg, fNaN, fConserve, fNWeightUnphys, fNWeightReject, fLundRecorded, fNFile; +unsigned long long int fNEvents, fNRecorded, fNGenerated, fWSqNeg, fNSigmaNeg, fNaN, fConserve, fNWeightUnphys, fNWeightReject, fLundRecorded, fNFile, fSolveEvents_0Sol, fSolveEvents_1Sol, fSolveEvents_2Sol; // SJDK 03/04/23 - Added in Qsq Min/Max, W Min/Max and t max (06/09/23) 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, fT_Max; double fOmega_Theta_I, fOmega_Theta_F, fOmega_Theta_Col, fOmega_Phi_Col; -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 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; @@ -240,7 +240,6 @@ void pim::Initilize() { // 18/01/23 - The luminosity below is some default assumtpion, more up to date values are set in DEMP prod and depend upon beam energy combinations if they are specified // See slide 11 in https://indico.cern.ch/event/1072579/contributions/4796856/attachments/2456676/4210776/CAP-EIC-June-7-2022-Seryi-r2.pdf for more info // fLumi = 0.374e33; // Jlab design - //fLumi = 1e34; // https://eic.jlab.org/wiki/index.php/EIC_luminosity - OUTDATED fLumi = 1e33; // 18/01/23, this seems a better default based upon more up to date info, see link above fuBcm2 = 1.0e-30; fPI = 3.1415926; @@ -279,14 +278,12 @@ void pim::Initilize() { fRecoilProton_Mass_GeV = fRecoilProton_Mass/1000.0; fPion_Mass = 139.57039 ; fPion_Mass_GeV = fPion_Mass/1000.0; - fKaon_Mass = 493.677; fKaon_Mass_GeV = fKaon_Mass/1000.0; fLambda_Mass = 1115.683; fLambda_Mass_GeV = fLambda_Mass/1000.0; fSigma_Mass = 1192.642; fSigma_Mass_GeV = fSigma_Mass/1000.0; - fOmega_Mass = 782.66; fOmega_Mass_GeV = fOmega_Mass/1000.0; @@ -333,6 +330,9 @@ void pim::Initilize() { fConserve = 0; fNWeightUnphys = 0; fNWeightReject = 0; + fSolveEvents_0Sol = 0; + fSolveEvents_1Sol = 0; + fSolveEvents_2Sol = 0; fSDiff = 0; fScatElecEnergyLess = 0; fScatElecThetaLess = 0; diff --git a/src/eic_evgen/eic_pim.h b/src/eic_evgen/eic_pim.h index 3198c79..f022c9d 100644 --- a/src/eic_evgen/eic_pim.h +++ b/src/eic_evgen/eic_pim.h @@ -116,6 +116,10 @@ extern unsigned long long int fConserve; extern unsigned long long int fNWeightUnphys; extern unsigned long long int fNWeightReject; +extern unsigned long long int fSolveEvents_0Sol; +extern unsigned long long int fSolveEvents_1Sol; +extern unsigned long long int fSolveEvents_2Sol; + extern unsigned long long int fNFile; extern double fK; diff --git a/src/eic_evgen/process_routine/DEMP_Reaction.cc b/src/eic_evgen/process_routine/DEMP_Reaction.cc index 8879234..088b8d7 100644 --- a/src/eic_evgen/process_routine/DEMP_Reaction.cc +++ b/src/eic_evgen/process_routine/DEMP_Reaction.cc @@ -54,7 +54,7 @@ void DEMP_Reaction::process_reaction() { rNEvent_itt = i; fNGenerated ++; - Progress_Report(); // This is happens at each 10% of the total event is processed + Progress_Report(); // This happens at each 10% of the total event is processed Processing_Event(); } @@ -106,12 +106,14 @@ void DEMP_Reaction::Init() { fElectron_Kin_Col = fElectron_Kin_Col_GeV * 1000.0; // cout << rNEvents << " " << fNEvents << endl; - + + // 08/09/23 - SJDK - Fermi momentum commented out for now, this is not fully implemented yet + // In future, this will be enabled/disabled automatically depending upon the specified hadron beam //rFermiMomentum = pd->fermiMomentum(); // ---------------------------------------------------- // Proton in collider (lab) frame - + // 08/09/23 - SJDK - The naming needs to be adjusted to be the incoming hadron beam, not the proton. Make this generic. r_lproton = GetProtonVector_lab(); r_lprotong = GetProtonVector_lab() * fm; @@ -134,7 +136,7 @@ void DEMP_Reaction::Init() { r_lelectrong = r_lelectron * fm; ///*--------------------------------------------------*/ - /// Getting the particle mass from the data base + /// Getting the ejectile (produced meson) particle mass from the data base produced_X = ParticleEnum(rEjectile); f_Ejectile_Mass = ParticleMass(produced_X)*1000; //MeV @@ -208,46 +210,42 @@ 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; } - // /*--------------------------------------------------*/ - // // SJDK 03/04/22 - New set of initialisation stuff for the solve function from Ishan and Bill - - // CoinToss = new TRandom3(); - - // // F = new TF1("F", - // // "[6]-sqrt([7]**2+x**2)-sqrt([8]**2+([3]-[0]*x)**2+([4]-[1]*x)**2+([5]-[2]*x)**2)", - // // 0, 12000); + if(UseSolve == true){ + /*--------------------------------------------------*/ + // SJDK 03/04/22 - New set of initialisation stuff for the solve function from Ishan and Bill + + CoinToss = new TRandom3(); - // F = new TF1("F", - // "[6]-sqrt([7]**2+x**2)-sqrt([8]**2+([3]-[0]*x)**2+([4]-[1]*x)**2+([5]-[2]*x)**2)", - // 0, r_lproton.E()); + F = new TF1("F", + "[6]-sqrt([7]**2+x**2)-sqrt([8]**2+([3]-[0]*x)**2+([4]-[1]*x)**2+([5]-[2]*x)**2)", + 0, r_lproton.E()); - // char AngleGenName[100] = "AngleGen"; - // double dummy[2] = {0,1}; - // // Changed the theta range here to match the one actually provided in the input .json file, these are already converted to radians in eic.cc (see ~ line 293) - // // double ThetaRange[2] = {fX_Theta_I, fX_Theta_F}; + char AngleGenName[100] = "AngleGen"; + double dummy[2] = {0,1}; - // f_Ejectile_Theta_I = fEjectileX_Theta_I; - // f_Ejectile_Theta_F = fEjectileX_Theta_F; + f_Ejectile_Theta_I = fEjectileX_Theta_I; + f_Ejectile_Theta_F = fEjectileX_Theta_F; - // double ThetaRange[2] = {f_Ejectile_Theta_I, f_Ejectile_Theta_F}; - // double PhiRange[2] = {0, 360*TMath::DegToRad()}; + double ThetaRange[2] = {f_Ejectile_Theta_I, f_Ejectile_Theta_F}; + double PhiRange[2] = {0, 360*TMath::DegToRad()}; - // AngleGen = new CustomRand(AngleGenName, dummy, - // ThetaRange, PhiRange); + AngleGen = new CustomRand(AngleGenName, dummy, + ThetaRange, PhiRange); - // UnitVect = new TVector3(0,0,1); + UnitVect = new TVector3(0,0,1); - // ///*--------------------------------------------------*/ - // // Produced hadron and recoilded hadron from the solve function + // ///*--------------------------------------------------*/ + // // Produced hadron and recoilded hadron from the solve function - // VertBeamElec = new TLorentzVector(); - // VertScatElec = new TLorentzVector(); + VertBeamElec = new TLorentzVector(); + VertScatElec = new TLorentzVector(); - // Initial = new TLorentzVector(); - // Target = new TLorentzVector(); - // Photon = new TLorentzVector(); - // Interaction = new TLorentzVector(); - // Final = new TLorentzVector(); + Initial = new TLorentzVector(); + Target = new TLorentzVector(); + Photon = new TLorentzVector(); + Interaction = new TLorentzVector(); + Final = new TLorentzVector(); + } } @@ -276,10 +274,9 @@ void DEMP_Reaction::Processing_Event() { fScatElec_Energy_Col = fRandom->Uniform( fScatElec_E_Lo * fElectron_Energy_Col , fScatElec_E_Hi * fElectron_Energy_Col ); // ---------------------------------------------------- - // Produced Particle X in Collider frame + // Produced ejectile in Collider frame // ---------------------------------------------------- - /// The generic produced particle in the exclusive reaction is labelled as X f_Ejectile_Theta_Col = acos( fRandom->Uniform( cos(f_Ejectile_Theta_I), cos(f_Ejectile_Theta_F ) ) ); f_Ejectile_Phi_Col = fRandom->Uniform( 0 , 2.0 * fPi ); @@ -388,9 +385,13 @@ void DEMP_Reaction::Processing_Event() { pow( f_Recoil_Mass , 2) ) ); } else if (UseSolve == true){ - r_l_Ejectile.SetPxPyPzE(0,0,0,0); - l_Recoil_g.SetPxPyPzE(0,0,0,0); + if(!Solve()){ + return; + } + r_l_Ejectile = r_l_Ejectile_solved; + l_Recoil = r_l_Recoil_solved; } + ///-------------------------------------------------- r_l_Ejectile_g = r_l_Ejectile * fm; @@ -409,16 +410,18 @@ void DEMP_Reaction::Processing_Event() { return; } -///*--------------------------------------------------*/ -//-> 10/05/23 - Love added a slimmed down, simpler to read version of the CheckLaws fn -// -// To check the conservation of the energy and momentum, there two methods avalaible: -// Method 1: Give the four-vectors of the initial and final states partciles, -// tolerance factor will be defaulted 1e-6 MeV -// Method 2: Give the four-vectors of the initial and final states partciles, -// and the prefered tolerance factor. -// - + ///*--------------------------------------------------*/ + //-> 10/05/23 - Love added a slimmed down, simpler to read version of the CheckLaws fn + // + // To check the conservation of the energy and momentum, there two methods avalaible: + // Method 1: Give the four-vectors of the initial and final states partciles, + // tolerance factor will be defaulted 1e-6 MeV + // CheckLaws(e_beam, h_beam, scatt_e, ejectile, recoil) <- input 4 vectors + // Method 2: Give the four-vectors of the initial and final states partciles, + // and the prefered tolerance factor. + // CheckLaws(e_beam, h_beam, scatt_e, ejectile, recoil, tolerance) <- input 4 vectors and tolerance value in GeV + // Both functions return 1 if conservations laws are satisified + if( pd->CheckLaws(r_lelectron, r_lproton, r_lscatelec, r_l_Ejectile, l_Recoil) !=1 ){ fConserve++; return; @@ -787,19 +790,35 @@ void DEMP_Reaction::Detail_Output() { DEMPDetails << "Seed used for the Random Number Generator " << setw(20) << fSeed << endl; DEMPDetails << endl; DEMPDetails << "Total events tried " << setw(20) << fNGenerated << endl; - DEMPDetails << "Total events cut " << setw(20) << (qsq_ev + w_ev + w_neg_ev + fNaN + fConserve + t_ev + fNSigmaNeg) << endl; - DEMPDetails << "Total events recorded " << setw(20) << fNRecorded << endl; - if (fNGenerated != (qsq_ev + w_ev + w_neg_ev + fNaN + fConserve + t_ev + fNSigmaNeg + fNRecorded)){ - DEMPDetails << "Total events cut + recorded = events tried? " << setw(20) << "NO! ERROR!" << endl; + if(UseSolve == true){ + DEMPDetails << "Total events cut " << setw(20) << (qsq_ev + w_ev + w_neg_ev + fNaN + fConserve + t_ev + fNSigmaNeg + fSolveEvents_0Sol) << endl; + DEMPDetails << "Total events recorded " << setw(20) << fNRecorded << endl; + if (fNGenerated != (qsq_ev + w_ev + w_neg_ev + fNaN + fConserve + t_ev + fNSigmaNeg + fNRecorded + fSolveEvents_0Sol)){ + DEMPDetails << "Total events cut + recorded = events tried? " << setw(20) << "NO! ERROR!" << endl; + } + else{ + DEMPDetails << "Total events cut + recorded = events tried? " << setw(22) << "Yes! :)" << endl; + } + } else{ - DEMPDetails << "Total events cut + recorded = events tried? " << setw(22) << "Yes! :)" << endl; + DEMPDetails << "Total events cut " << setw(20) << (qsq_ev + w_ev + w_neg_ev + fNaN + fConserve + t_ev + fNSigmaNeg) << endl; + DEMPDetails << "Total events recorded " << setw(20) << fNRecorded << endl; + if (fNGenerated != (qsq_ev + w_ev + w_neg_ev + fNaN + fConserve + t_ev + fNSigmaNeg + fNRecorded)){ + DEMPDetails << "Total events cut + recorded = events tried? " << setw(20) << "NO! ERROR!" << endl; + } + else{ + DEMPDetails << "Total events cut + recorded = events tried? " << setw(22) << "Yes! :)" << endl; + } } DEMPDetails << endl << "Cut details -" << endl; DEMPDetails << "Events cut due to qsq < " << fQsq_Min << " or qsq > "<< fQsq_Max << " " << setw(20) << qsq_ev << endl; DEMPDetails << "Events cut due to negative Wsq value " << setw(20) << w_neg_ev << endl; DEMPDetails << "Events cut due to W < " << fW_Min << " or W > " << fW_Max << " " << setw(20) << w_ev << endl; + if(UseSolve == true){ + DEMPDetails << "Events cut due to solve function finding 0 solutions " << setw(20) << fSolveEvents_0Sol << endl; + } DEMPDetails << "Events cut due to ejectile (X) energy NaN " << setw(20) << fNaN << endl; DEMPDetails << "Events cut due to conservation law check failure " << setw(20) << fConserve << endl; DEMPDetails << "Events cut due to -t > " << fT_Max << "GeV " << setw(30) << t_ev << endl; @@ -818,6 +837,12 @@ void DEMP_Reaction::Detail_Output() { DEMPDetails << "Events cut due to py and pz conservation law checks " << setw(20) << mom_pypz << endl; DEMPDetails << "Events cut due to px, py and pz conservation law checks " << setw(20) << mom_pxpypz << endl; + if(UseSolve == true){ + DEMPDetails << endl << "Solve function, addtional info -" << endl; + DEMPDetails << "Number of events with 0 Solution " << setw(20) << fSolveEvents_0Sol << endl; + DEMPDetails << "Number of events with 1 Solution " << setw(20) << fSolveEvents_1Sol << endl; + DEMPDetails << "Number of events with 2 Solution " << setw(20) << fSolveEvents_2Sol << endl; + } } @@ -1041,7 +1066,7 @@ void DEMP_Reaction::DEMPReact_HEPMC3_Out_Init() { void DEMP_Reaction::DEMPReact_HEPMC3_Output() { - // HEPMC3 output for Athena/EPIC simulations + // HEPMC3 output for Athena/ePIC simulations // First line - E - Event# - #Vertices - #Particles DEMPOut << std::scientific << std::setprecision(15) << "E" << " " << print_itt << " " << "1" << " " << 5 << endl; @@ -1064,3 +1089,167 @@ void DEMP_Reaction::DEMPReact_HEPMC3_Output() { 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; } + +/*--------------------------------------------------*/ + +bool DEMP_Reaction::SolnCheck() +{ + // + // // Double Checking for solution viability + // if (TMath::Abs(f_Scat_hadron_Mass-r_l_scat_hadron_solved->M())>1){ + // //cerr << "Mass Missmatch" << endl; + // //cerr << TMath::Abs(proton_mass_mev-Proton->M()) << endl; + // return false; + // } + // if (TMath::Abs(W_in()-W_out())>1){ + // //cerr << "W Missmatch" << endl; + // //cerr << TMath::Abs(W_in()-W_out()) << endl; + // return false; + // } + // *Final = *r_l_scat_hadron_solved + *r_lX_solved; + // + // if (TMath::Abs(Initial->Px()-Final->Px())>1){ + // //cerr << "Px Missmatch" << endl; + // //cerr << TMath::Abs(Initial->Px()-Final->Px()) << endl; + // return false; + // } + // + // if (TMath::Abs(Initial->Py()-Final->Py())>1){ + // //cerr << "Py Missmatch" << endl; + // //cerr << TMath::Abs(Initial->Py()-Final->Py()) << endl; + // return false; + // } + // + // if (TMath::Abs(Initial->Pz()-Final->Pz())>1){ + // //cerr << "Pz Missmatch" << endl; + // //cerr << TMath::Abs(Initial->Pz()-Final->Pz()) << endl; + // return false; + // } + // + // if (TMath::Abs(Initial->E()-Final->E())>1){ + // return false; + // } + return true; +} + +/*--------------------------------------------------*/ +double DEMP_Reaction::W_in() +{ + return 0; +} + +/*--------------------------------------------------*/ +double DEMP_Reaction::W_out() +{ + return 0; +} + +/*--------------------------------------------------*/ + +int DEMP_Reaction::Solve() +{ + + + VertBeamElec->SetPxPyPzE(r_lelectron.Px(), r_lelectron.Py(), r_lelectron.Pz(), r_lelectron.E()); + VertScatElec->SetPxPyPzE(r_lscatelec.Px(), r_lscatelec.Py(), r_lscatelec.Pz(), r_lscatelec.E()); + Target->SetPxPyPzE(r_lproton.Px(), r_lproton.Py(), r_lproton.Pz(), r_lproton.E()); + *Photon = *VertBeamElec - *VertScatElec; + *Interaction = *Photon; + + *Initial = *Interaction+*Target; + + theta = f_Ejectile_Theta_Col; + phi = f_Ejectile_Phi_Col; + + return this->Solve(theta, phi); +} + + +int DEMP_Reaction::Solve(double theta, double phi) +{ + + W_in_val = W_in(); + + if (W_in_val<0){ + return 0; + } + + UnitVect->SetTheta(theta); + UnitVect->SetPhi(phi); + UnitVect->SetMag(1); + + double* pars = new double[9]; + + pars[0] = UnitVect->X(); + pars[1] = UnitVect->Y(); + pars[2] = UnitVect->Z(); + pars[3] = Initial->Px(); + pars[4] = Initial->Py(); + pars[5] = Initial->Pz(); + pars[6] = Initial->E(); + pars[7] = f_Ejectile_Mass; + pars[8] = f_Recoil_Mass; + + F->SetParameters(pars); + + + ///*--------------------------------------------------*/ + // Looking for the 1st Solution: + // If a solution found, then this will be the fist solution. Then we proceed to look for the 2nd solution. + // If no soluion found, then exit solve function + + P = F->GetX(0, 0, pars[6], 0.0001, 10000); + + if (TMath::Abs(F->Eval(P)) < 1){ + fSolveEvents_1Sol++; + } else { + fSolveEvents_0Sol++; + return 0; + } + + TLorentzVector * r_l_Ejectile_solved_1_temp = new TLorentzVector(); + TLorentzVector * r_l_Ejectile_solved_2_temp = new TLorentzVector(); + + 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) ); + r_l_Ejectile_solved_1_temp->SetPxPyPzE(P*pars[0], P*pars[1], P*pars[2], r_l_Ejectile_E); + + ///*--------------------------------------------------*/ + // Looking for the 2nd Solution + + P2 = F->GetX(0, P+100, pars[6], 0.0001, 10000); + 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) ); + r_l_Ejectile_solved_2_temp->SetPxPyPzE(P2 * pars[0], P2 * pars[1], P2 * pars[2], r_l_Ejectile_E_2); + + ///*--------------------------------------------------*/ + // If a valid 2nd solution is found, then we are certian that there are two solutions. + // - We then increament the counter for 2nd solution scenario + // - We then decreament the counter for the 1st solution scenario + + if (TMath::Abs(F->Eval(P2)) < 1){ + fSolveEvents_2Sol++; + fSolveEvents_1Sol--; + if ( Int_t(CoinToss->Uniform(0,100)) < 50) { + 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()); + } else { + 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()); + } + } + else { + 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()); + } + + ///*--------------------------------------------------*/ + /// Solve for the recoil information with the "solved" Ejectile informaiton + TLorentzVector * r_l_hadron_temp= new TLorentzVector(); + *r_l_hadron_temp = *Initial- r_l_Ejectile_solved; + 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()); + + delete r_l_Ejectile_solved_1_temp; + delete r_l_Ejectile_solved_2_temp; + + delete r_l_hadron_temp; + delete[] pars; + + return 1; + +} diff --git a/src/eic_evgen/reaction_routine.cc b/src/eic_evgen/reaction_routine.cc index 90ece05..6112507 100644 --- a/src/eic_evgen/reaction_routine.cc +++ b/src/eic_evgen/reaction_routine.cc @@ -78,17 +78,10 @@ Reaction::~Reaction() { /// void Reaction::process_reaction() { -// if (rEjectile == "Pi0") { -// // Pi0_Production* r1 = new Pi0_Production("Eta"); -// Pi0_Production* rr1 = new Pi0_Production(rEjectile); -// rr1->process_reaction(); -// delete rr1; -// } -// else{ + // SJDK - 19/12/22 - New generic DEMP reaction class, the intention is that this should be able to handle any case DEMP_Reaction* rr1 = new DEMP_Reaction(rEjectile, rRecoil); rr1->process_reaction(); delete rr1; -// } } diff --git a/src/eic_evgen/reaction_routine.h b/src/eic_evgen/reaction_routine.h index c8e388e..42ebd59 100644 --- a/src/eic_evgen/reaction_routine.h +++ b/src/eic_evgen/reaction_routine.h @@ -74,7 +74,7 @@ class DEMP_Reaction { Double_t Get_Total_Cross_Section(); -// Double_t GetPi0_CrossSection(); + // Double_t GetPi0_CrossSection(); /*--------------------------------------------------*/ // Parameters @@ -87,7 +87,7 @@ class DEMP_Reaction { std::string sTFile; /// Generator output files. For documentation and monitoring purposes std::string sLFile; /// Lund input file into the EIC simulation - std::string sDFile; /// Root dianostic plot in root file format + std::string sDFile; /// Root diagnostic plot in root file format std::ofstream DEMPOut; std::ofstream DEMPDetails; @@ -137,11 +137,8 @@ class DEMP_Reaction { TLorentzVector r_l_Ejectile; TLorentzVector r_l_Ejectile_g; -// Particle* r_l_Ejectile_solved; -// Particle* l_Recoil_solved; - TLorentzVector r_l_Ejectile_solved; - TLorentzVector l_Recoil_solved; + TLorentzVector r_l_Recoil_solved; double f_Ejectile_Mass; double f_Ejectile_Mass_GeV; @@ -220,28 +217,19 @@ class DEMP_Reaction { ///*--------------------------------------------------*/ // Rory Check algorithm + + TLorentzVector* Interaction_Solve; + TLorentzVector* Target_Solve; + + TLorentzVector* VertBeamElec; + TLorentzVector* VertScatElec; - // Particle* Interaction; - // Particle* Target; - // - // Particle* Initial; - // Particle* Final; - // - // Particle* VertBeamElec; - // Particle* VertScatElec; - // Particle* Photon; - - - TLorentzVector Interaction_Solve; - TLorentzVector Target_Solve; - - TLorentzVector Initial; - TLorentzVector Final; + TLorentzVector* Initial; + TLorentzVector* Target; + TLorentzVector* Photon; + TLorentzVector* Interaction; + TLorentzVector* Final; - TLorentzVector VertBeamElec; - TLorentzVector VertScatElec; - TLorentzVector Photon; - bool SolnCheck(); double W_in_Solve(); double W_out_Solve(); @@ -256,6 +244,9 @@ class DEMP_Reaction { int Solve(); int Solve(double theta, double phi); + double W_in(); + double W_out(); + ///*--------------------------------------------------*/ // Needed for the Solve function From b4694693b36eeb02f82dbddb16abbd334acc6e07 Mon Sep 17 00:00:00 2001 From: Stephen JD Kay Date: Fri, 8 Sep 2023 07:45:00 -0600 Subject: [PATCH 12/21] Pedantic comment typo fix --- src/eic_evgen/eic.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/eic_evgen/eic.cc b/src/eic_evgen/eic.cc index 5c2b57f..c5ea328 100644 --- a/src/eic_evgen/eic.cc +++ b/src/eic_evgen/eic.cc @@ -421,7 +421,7 @@ vector>>> ReadCrossSectionPar(TString particle, TSt sigT_ParamFile = "../src/eic_evgen/CrossSection_Params/KPlusSigma_Param_sigT"; } else if (particle == "Pi0"){ - // When pi0 model paramterised, add it here + // When pi0 model parameterised, add it here } else{ cerr << " !!!!! " << endl << "Warning!" << endl << "Combination of specified ejectile and recoil particles not found!" << "Cross section parameters cannot be read, check inputs!" << endl << "Warning!" << endl << " !!!!! " << endl; From 9ddc3e9297eb2069c230b858acf1fe49568292dd Mon Sep 17 00:00:00 2001 From: Stephen Kay Date: Wed, 13 Sep 2023 04:37:19 -0600 Subject: [PATCH 13/21] Added info to README regarding calculation methods. --- README.md | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/README.md b/README.md index cf82515..1781bd4 100644 --- a/README.md +++ b/README.md @@ -38,6 +38,18 @@ Event generator for Deep Exclusive Meson Production The file Config.json contains all the configuration options. Use this as a template for other configuration files, which may be given as an argument to the event generator. +## Ejectile calculation methods + +- The EIC module of DEMPgen has two different calculation methods that may be used to calculate the ejectile properties. +- These are referred to as the "Analytical" and "Solve" methods. Either can be used in a simulation run, just change the config file to switch between them. +- In your config .json file, select between the two versions using the + -"calc_method" argument, this can be set to Analytical or Solve, e.g. + - "calc_method": "Analytical", + - OR + - "calc_method": "Solve", + +- If you are using the shell scripts provided to run a lot of simulations, modify Config_EIC.json BEFORE running the shell scripts. + ## Output The event data is output to the configured file location relative to the build directory. The TTree in this file contains all kinematic data for all particles in the laboratory rest frame. Variables with the prefix "Vert" represent values read at the interaction vertex. Values with the prefix "Lab" represent values read after all correcting effects (multiple scattering, ionization, etc.) have been applied. From 45dc3287ce3e861ca2a02f4696535de8d558e47e Mon Sep 17 00:00:00 2001 From: Stephen JD Kay Date: Wed, 13 Sep 2023 04:39:55 -0600 Subject: [PATCH 14/21] Update README.md --- README.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/README.md b/README.md index 1781bd4..e3ae010 100644 --- a/README.md +++ b/README.md @@ -38,12 +38,12 @@ Event generator for Deep Exclusive Meson Production The file Config.json contains all the configuration options. Use this as a template for other configuration files, which may be given as an argument to the event generator. -## Ejectile calculation methods +### Ejectile calculation methods - The EIC module of DEMPgen has two different calculation methods that may be used to calculate the ejectile properties. - These are referred to as the "Analytical" and "Solve" methods. Either can be used in a simulation run, just change the config file to switch between them. - In your config .json file, select between the two versions using the - -"calc_method" argument, this can be set to Analytical or Solve, e.g. + - "calc_method" argument, this can be set to Analytical or Solve, e.g. - "calc_method": "Analytical", - OR - "calc_method": "Solve", From 6a0c0dfd4207635abd4578fc218c4b28a3b628fe Mon Sep 17 00:00:00 2001 From: Stephen JD Kay Date: Wed, 13 Sep 2023 04:42:28 -0600 Subject: [PATCH 15/21] Update README.md Clarifications of calc method in .json files. Formatting fixes. --- README.md | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/README.md b/README.md index e3ae010..5366c08 100644 --- a/README.md +++ b/README.md @@ -42,11 +42,11 @@ The file Config.json contains all the configuration options. Use this as a templ - The EIC module of DEMPgen has two different calculation methods that may be used to calculate the ejectile properties. - These are referred to as the "Analytical" and "Solve" methods. Either can be used in a simulation run, just change the config file to switch between them. -- In your config .json file, select between the two versions using the - - "calc_method" argument, this can be set to Analytical or Solve, e.g. - - "calc_method": "Analytical", - - OR - - "calc_method": "Solve", +- In your config .json file, select between the two versions using the "calc_method" argument, this can be set to Analytical or Solve, e.g. + - "calc_method": "Analytical" + - OR + - "calc_method": "Solve" + - If the method is not entered exactly as above, it will not be recognised and the generator will default to the analytical method. - If you are using the shell scripts provided to run a lot of simulations, modify Config_EIC.json BEFORE running the shell scripts. From 0f68b90d0204aa9d6362dea375f14511ce975d81 Mon Sep 17 00:00:00 2001 From: Stephen JD Kay Date: Wed, 13 Sep 2023 04:46:35 -0600 Subject: [PATCH 16/21] Update README.md Added comment to script running section to remind users to change calc method before running. --- README.md | 63 ++++++++++++++++++++++++++++--------------------------- 1 file changed, 32 insertions(+), 31 deletions(-) diff --git a/README.md b/README.md index 5366c08..c2728d2 100644 --- a/README.md +++ b/README.md @@ -1,5 +1,5 @@ # DEMPgen -Event generator for Deep Exclusive Meson Production +Event generator for Deep Exclusive Meson Production. ## Building @@ -25,7 +25,7 @@ Event generator for Deep Exclusive Meson Production ### Building on the iFarm -- Building on the JLab iFarm requires you to set up some software versions beforehand, to build successfully, I did the following +- Building on the JLab iFarm requires you to set up some software versions beforehand, to build successfully, I did the following - - Comment out any CUE or other initialisation in your .login/.cshrc scripts - Login to an ifarm node - module load cmake/3.19.4 @@ -108,28 +108,29 @@ This project uses [JsonCpp](https://github.com/open-source-parsers/jsoncpp "Json ## Processing Scripts and json examples -- There are several shell scripts in the main directory that can be used to run the generator and produce some output, information on these scripts is provided below -- To remove clutter from the main directory, example .json scripts were moved to a new folder +- There are several shell scripts in the main directory that can be used to run the generator and produce some output, information on these scripts is provided below. +- To remove clutter from the main directory, example .json scripts were moved to a new folder - - json_examples +- Remember to set the ejectile calculation method in the Config_EIC.json file *before* running these scripts if you do not want to use the default analytical method. ### Process_EIC.csh !!! NOTICE !!! -This script copies Config_EIC.json and formats a new file based upon this, DO NOT MODIFY Config_EIC.json if you want to use this script! +This script copies Config_EIC.json and formats a new file based upon this, DO NOT MODIFY Config_EIC.json (other than the calculation method) if you want to use this script! !!! NOTICE !!! - To facilitate the submission of batch jobs, I created a csh script to automatically construct .json config files and run them. This script can also be utilised to run the generator manually, without the need to go and edit a json file. - The script requires 8 arguments (which is a lot, I know), but in the K+ case, it expects 9. They are as follows - - - Arg 1 - FileNum -> For batch running, we typically run X files of Y events, this argument is just X, if you're running manually as a test, just input 1 or whatever you fancy - - 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 - 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 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 + - Arg 1 - FileNum -> For batch running, we typically run X files of Y events, this argument is just X, if you're running manually as a test, just input 1 or whatever you fancy. + - 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 - 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 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. - So as an example if you executed the following - @@ -142,33 +143,33 @@ This script copies Config_EIC.json and formats a new file based upon this, DO NO - This script creates and submits batch jobs. It is designed for use with the torque queueing system on Lark at the University of Regina. The jobs the script creates and submits all execute the Process_EIC.csh script described above. This script requries a very similar set of arguments - - - Arg 1 - NumFiles -> The batch script is designed to run X jobs of Y events, this number is just X, the number of files you want to run - - 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 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 + - Arg 1 - NumFiles -> The batch script is designed to run X jobs of Y events, this number is just X, the number of files you want to run. + - 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 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 invalif. + - 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. -- The script automatically generates a random seed itself using the /dev/urandom function +- The script automatically generates a random seed itself using the /dev/urandom function. ### Process_EIC_iFarm.csh -- This version is for use on the JLab iFarm/Farm -- This script uses the same arguments as Process_EIC.csh +- This version is for use on the JLab iFarm/Farm. +- This script uses the same arguments as Process_EIC.csh. - If you are processing interactively (i.e. on the iFarm), you will need to make sure that you have executed - - module use /group/halla/modulefiles - module load root -- You should also check the path set at the top looks OK - - By default, /eic/users/${USER}/DEMPGen is assumed +- You should also check the path set at the top looks OK. + - By default, /eic/users/${USER}/DEMPGen is assumed. ### JLab_Batch_Submission.sh -- This version should be used to submit jobs to the Farm batch queueing system (swif2) -- It uses the same arguments as Batch_Submission_EIC.sh -- You should also check the paths set throughout the script look ok (for example, in the COMMAND: ... line) - - By default, /eic/users/${USER}/DEMPGen is assumed +- This version should be used to submit jobs to the Farm batch queueing system (swif2). +- It uses the same arguments as Batch_Submission_EIC.sh. +- You should also check the paths set throughout the script look ok (for example, in the COMMAND: ... line). + - By default, /eic/users/${USER}/DEMPGen is assumed. ### json_examples From 899d3e220661c0647fed446f1bee9eb60011c5fb Mon Sep 17 00:00:00 2001 From: Stephen JD Kay Date: Wed, 13 Sep 2023 04:48:13 -0600 Subject: [PATCH 17/21] Update README.md Minor fix of a code line, removed . so it can be directly copy/pasted. --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index c2728d2..2beb7c4 100644 --- a/README.md +++ b/README.md @@ -9,7 +9,7 @@ Event generator for Deep Exclusive Meson Production. - After downloading the source create a build directory and cd to it. Take note of the location of the source directory (where CMakeLists.txt should be stored) and run the commands: - - mdkir build. + - mdkir build - cd build - cmake .. - make -j8 From c0ca9ff316093d7964fe3eda24a5d7ac37691c61 Mon Sep 17 00:00:00 2001 From: Stephen JD Kay Date: Wed, 13 Sep 2023 04:51:42 -0600 Subject: [PATCH 18/21] Update README.md Clarified config files section. This referred to a file that no longer exists. Clarified that there are specific EIC/SoLID templates for use. --- README.md | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/README.md b/README.md index 2beb7c4..6b94f9b 100644 --- a/README.md +++ b/README.md @@ -36,7 +36,10 @@ Event generator for Deep Exclusive Meson Production. ## Configuration -The file Config.json contains all the configuration options. Use this as a template for other configuration files, which may be given as an argument to the event generator. +The file Config_EIC.json and Config_SoLID.json contain all the configuration options for EIC and SoLID event generation respectively. Use these as a template for other configuration files, which may be given as an argument to the event generator. +- Note that Config_EIC.json is used in shell scripts (see below) that run large numbers of jobs on batch queueing systems. *Do not* modify this file unless you know what you're doing with it or you really need to edit something. + - Copying the template to a new file and editing that is strongly recommended in all cases. + ### Ejectile calculation methods From e587cadde73be38e193416165e60f9fa1b16def4 Mon Sep 17 00:00:00 2001 From: Stephen Kay Date: Wed, 13 Sep 2023 08:03:07 -0600 Subject: [PATCH 19/21] Renamed parts of eic.cc, eic.h and reaction routines to match more consistent conventions. Modified config template and shell scripts to match. --- .gitignore | 4 +- Batch_Submission_EIC.sh | 26 +++---- Config_EIC.json | 4 +- JLab_Batch_Submission.sh | 22 +++--- Process_EIC.csh | 28 +++---- Process_EIC_iFarm.csh | 28 +++---- README.md | 8 +- src/eic_evgen/eic.cc | 118 +++++++++++++++--------------- src/eic_evgen/eic.h | 2 +- src/eic_evgen/reaction_routine.cc | 21 +++--- src/eic_evgen/reaction_routine.h | 8 +- 11 files changed, 135 insertions(+), 134 deletions(-) diff --git a/.gitignore b/.gitignore index 2c26578..3cdaf75 100644 --- a/.gitignore +++ b/.gitignore @@ -9,4 +9,6 @@ build/ *_job.txt Config_EIC_* sjdkay_job.txt -*#* \ No newline at end of file +*#* +.cc* +.h* \ No newline at end of file diff --git a/Batch_Submission_EIC.sh b/Batch_Submission_EIC.sh index 4f54b3b..bec4b33 100755 --- a/Batch_Submission_EIC.sh +++ b/Batch_Submission_EIC.sh @@ -14,7 +14,7 @@ fi if [[ "$#" -ne 7 && "$#" -ne 8 ]]; then echo "" echo "!!! ERROR !!! - Expected 7 or 8 arguments - !!! ERROR !!!" - echo "Expect - NumFiles NumEvents EBeamE HBeamE OutputType InteractionPoint Particle Hadron(optional)" + echo "Expect - NumFiles NumEvents EBeamE HBeamE OutputType InteractionPoint Ejectile RecoilHadron(optional)" echo "See the Config_EIC.json file or the README for options and try again, exiting" echo "!!! ERROR !!! - Expected 7 or 8 arguments - !!! ERROR !!!" echo "" @@ -28,35 +28,35 @@ EBeamE=$3 HBeamE=$4 OutputType=$5 InteractionPoint=$6 -Particle=$7 +Ejectile=$7 # If K+ specified, check the 8th argument, expect this to exist for K+, if it does NOT (case 1), set a default -if [[ $Particle == "K+" && -z "$8" ]]; then +if [[ $Ejectile == "K+" && -z "$8" ]]; then echo "!!! WARNING !!! - For K+ production expect a hadron specified, defaulting to Lambda - !!! WARNING !!!" - Hadron="Lambda" -elif [[ $Particle == "K+" && ! -z "$8" ]]; then # If 8th argument is not a blank string (i.e. it exists), set the Hadron to this - Hadron=$8 -else # Any other case (non K+), set Hadron to be a blank string. We don't actually care for Pi+, Pi0 production etc. - Hadron="" + RecoilHadron="Lambda" +elif [[ $Ejectile == "K+" && ! -z "$8" ]]; then # If 8th argument is not a blank string (i.e. it exists), set the RecoilHadron to this + RecoilHadron=$8 +else # Any other case (non K+), set RecoilHadron to be a blank string. We don't actually care for Pi+, Pi0 production etc. + RecoilHadron="" fi i=1 while [[ $i -le $NumFiles ]]; do # This is the name of the job submission script the shell script creates - batch="${USER}_EICDempGen_${EBeamE}on${HBeamE}_${Particle}${Hadron}_${InteractionPoint}_${NumEvents}_${i}_Job.txt" # The name of the job submission script it'll create each time + batch="${USER}_EICDempGen_${EBeamE}on${HBeamE}_${Ejectile}${RecoilHadron}_${InteractionPoint}_${NumEvents}_${i}_Job.txt" # The name of the job submission script it'll create each time echo "Running ${batch} for file ${i}" cp /dev/null ${batch} RandomSeed=$(od -An -N3 -i /dev/urandom) echo "#!/bin/csh" >> ${batch} # Tells your job which shell to run in - echo "#PBS -N DEMPGen_${EBeamE}on${HBeamE}_${Particle}${Hadron}_${InteractionPoint}_${NumEvents}_${i}" >> ${batch} # Name your job + echo "#PBS -N DEMPGen_${EBeamE}on${HBeamE}_${Ejectile}${RecoilHadron}_${InteractionPoint}_${NumEvents}_${i}" >> ${batch} # Name your job echo "#PBS -m abe" >> ${batch} # Email you on job start, end or error #echo "#PBS -M ${USER}@jlab.org" >>${batch} # Your email address, change it to be what you like echo "#PBS -r n" >> ${batch} # Don't re-run if it crashes - echo "#PBS -o /home/${USER}/trq_output/${EBeamE}on${HBeamE}_${Particle}${Hadron}_${InteractionPoint}_${NumEvents}_${i}.out" >> ${batch} # Output directory and file name, set to what you like - echo "#PBS -e /home/${USER}/trq_output/${EBeamE}on${HBeamE}_${Particle}${Hadron}_${InteractionPoint}_${NumEvents}_${i}.err" >> ${batch} # Error output directory and file name + echo "#PBS -o /home/${USER}/trq_output/${EBeamE}on${HBeamE}_${Ejectile}${RecoilHadron}_${InteractionPoint}_${NumEvents}_${i}.out" >> ${batch} # Output directory and file name, set to what you like + echo "#PBS -e /home/${USER}/trq_output/${EBeamE}on${HBeamE}_${Ejectile}${RecoilHadron}_${InteractionPoint}_${NumEvents}_${i}.err" >> ${batch} # Error output directory and file name echo "date" >> ${batch} echo "cd /home/apps/DEMPgen/" >> ${batch} # Tell your job to go to the directory with the script you want to run - echo "./Process_EIC.csh ${i} ${NumEvents} ${EBeamE} ${HBeamE} ${RandomSeed} ${OutputType} ${InteractionPoint} ${Particle} ${Hadron}" >> ${batch} # Run your script, change this to what you like + echo "./Process_EIC.csh ${i} ${NumEvents} ${EBeamE} ${HBeamE} ${RandomSeed} ${OutputType} ${InteractionPoint} ${Ejectile} ${RecoilHadron}" >> ${batch} # Run your script, change this to what you like echo "date">>${batch} echo "exit">>${batch} # End of your job script echo "Submitting batch" diff --git a/Config_EIC.json b/Config_EIC.json index bd33599..6e7f2ea 100644 --- a/Config_EIC.json +++ b/Config_EIC.json @@ -23,8 +23,8 @@ //************************************** /// This section if for EIC simulation only "Kinematics_type" : 1, // Kinematics type (1->FF, 2->TSSA) - "particle": "Pion+", // Choices: omega, pi+, pi0, K+ - "hadron": "neutron", // Choices: Neutron, proton, Lambda or Sigma0 + "ejectile": "Pion+", // Choices: omega, pi+, pi0, K+ + "recoil_hadron": "neutron", // Choices: Neutron, proton, Lambda or Sigma0 "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 diff --git a/JLab_Batch_Submission.sh b/JLab_Batch_Submission.sh index fa1dad0..2c8158c 100755 --- a/JLab_Batch_Submission.sh +++ b/JLab_Batch_Submission.sh @@ -11,7 +11,7 @@ echo "Running as ${USER}" # Checks who you're running this as if [[ "$#" -ne 7 && "$#" -ne 8 ]]; then echo "" echo "!!! ERROR !!! - Expected 7 or 8 arguments - !!! ERROR !!!" - echo "Expect - NumFiles NumEvents EBeamE HBeamE OutputType InteractionPoint Particle Hadron(optional)" + echo "Expect - NumFiles NumEvents EBeamE HBeamE OutputType InteractionPoint Ejectile RecoilHadron(optional)" echo "See the Config_EIC.json file or the README for options and try again, exiting" echo "!!! ERROR !!! - Expected 7 or 8 arguments - !!! ERROR !!!" echo "" @@ -25,16 +25,16 @@ EBeamE=$3 HBeamE=$4 OutputType=$5 InteractionPoint=$6 -Particle=$7 +Ejectile=$7 # If K+ specified, check the 8th argument, expect this to exist for K+, if it does NOT (case 1), set a default -if [[ $Particle == "K+" && -z "$8" ]]; then +if [[ $Ejectile == "K+" && -z "$8" ]]; then echo "!!! WARNING !!! - For K+ production expect a hadron specified, defaulting to Lambda - !!! WARNING !!!" - Hadron="Lambda" -elif [[ $Particle == "K+" && ! -z "$8" ]]; then # If 8th argument is not a blank string (i.e. it exists), set the Hadron to this - Hadron=$8 -else # Any other case (non K+), set Hadron to be a blank string. We don't actually care for Pi+, Pi0 production etc. - Hadron="" + RecoilHadron="Lambda" +elif [[ $Ejectile == "K+" && ! -z "$8" ]]; then # If 8th argument is not a blank string (i.e. it exists), set the RecoilHadron to this + RecoilHadron=$8 +else # Any other case (non K+), set RecoilHadron to be a blank string. We don't actually care for Pi+, Pi0 production etc. + RecoilHadron="" fi Workflow="EIC_DEMPGen_${USER}" # Change this as desired @@ -47,16 +47,16 @@ while true; do ( while [[ $i -le $NumFiles ]]; do # This is the name of the job submission script the shell script creates - batch="${USER}_EICDempGen_${EBeamE}on${HBeamE}_${Particle}${Hadron}_${InteractionPoint}_${NumEvents}_${i}_Job.txt" # The name of the job submission script it'll create each time + batch="${USER}_EICDempGen_${EBeamE}on${HBeamE}_${Ejectile}${RecoilHadron}_${InteractionPoint}_${NumEvents}_${i}_Job.txt" # The name of the job submission script it'll create each time echo "Running ${batch} for file ${i}" cp /dev/null ${batch} RandomSeed=$(od -An -N3 -i /dev/urandom) echo "PROJECT: c-kaonlt" >> ${batch} # Is eic a valid project? echo "TRACK: analysis" >> ${batch} - echo "JOBNAME: DEMPGen_${EBeamE}on${HBeamE}_${Particle}${Hadron}_${InteractionPoint}_${NumEvents}_${i}" >> ${batch} + echo "JOBNAME: DEMPGen_${EBeamE}on${HBeamE}_${Ejectile}${RecoilHadron}_${InteractionPoint}_${NumEvents}_${i}" >> ${batch} echo "MEMORY: 2000 MB" >> ${batch} # Request 2GB RAM - probably too much echo "CPU: 1" >> ${batch} # Request 1 CPU core per job - echo "COMMAND:/group/eic/users/${USER}/DEMPGen/Process_EIC_iFarm.csh ${i} ${NumEvents} ${EBeamE} ${HBeamE} ${RandomSeed} ${OutputType} ${InteractionPoint} ${Particle} ${Hadron}" >> ${batch} + echo "COMMAND:/group/eic/users/${USER}/DEMPGen/Process_EIC_iFarm.csh ${i} ${NumEvents} ${EBeamE} ${HBeamE} ${RandomSeed} ${OutputType} ${InteractionPoint} ${Ejectile} ${RecoilHadron}" >> ${batch} echo "MAIL: ${USER}@jlab.org" >> ${batch} echo "Submitting batch" eval "swif2 add-jsub ${Workflow} -script ${batch} 2>/dev/null" # Swif2 job submission, uses old jsub scripts diff --git a/Process_EIC.csh b/Process_EIC.csh index 647028f..b0f03e2 100755 --- a/Process_EIC.csh +++ b/Process_EIC.csh @@ -4,7 +4,7 @@ echo"" echo "This file is intended to be run as part of a batch job submission, however, you can also run it on its own." -echo "Expected input is - FileNumber NumberOfEvents ElectronBeamEnergy HadronBeamEnergy RandomSeed OutputType IP Particle Hadron(Optional, for K+)" +echo "Expected input is - FileNumber NumberOfEvents ElectronBeamEnergy HadronBeamEnergy RandomSeed OutputType IP Ejectile RecoilHadron(Optional, for K+)" echo "Please see the README for more info." echo "" @@ -22,35 +22,35 @@ set HBeamE=$4 set RandomSeed=$5 set OutputType=$6 set InteractionPoint=$7 -set Particle=$8 +set Ejectile=$8 -if ($Particle == "K+" && $#argv == 8 ) then +if ($Ejectile == "K+" && $#argv == 8 ) then echo "! WARNING !" echo "! WARNING ! - For K+ production expect a hadron specified, defaulting to Lambda - ! WARNING !" echo "! WARNING !" - set Hadron="Lambda" -else if ($Particle == "K+" && $#argv == 9 ) then - set Hadron=$9 + set RecoilHadron="Lambda" +else if ($Ejectile == "K+" && $#argv == 9 ) then + set RecoilHadron=$9 else - set Hadron="" + set RecoilHadron="" endif -echo "Running target polarisation up, FF setting for file $FileNum with $NumEvents events per file for $EBeamE GeV e- on $HBeamE GeV p using random seed $RandomSeed, using $OutputType format output for $Particle $Hadron events." +echo "Running target polarisation up, FF setting for file $FileNum with $NumEvents events per file for $EBeamE GeV e- on $HBeamE GeV p using random seed $RandomSeed, using $OutputType format output for $Ejectile $RecoilHadron events." # Set the config file name based upon inputs -set ConfigFilename = 'Config_EIC_'$EBeamE'on'$HBeamE'_'$InteractionPoint'_'$Particle$Hadron'_'$NumEvents'_'$FileNum'.json' +set ConfigFilename = 'Config_EIC_'$EBeamE'on'$HBeamE'_'$InteractionPoint'_'$Ejectile$RecoilHadron'_'$NumEvents'_'$FileNum'.json' # Copy the default config file to our constructed filename cp Config_EIC.json $ConfigFilename # Use sed commands to change our config file based upon inputs -sed -i 's/"file_name" \:.*/"file_name" \: "DEMPGen_'$EBeamE'on'$HBeamE'_'$InteractionPoint'_'$Particle$Hadron'_'$NumEvents'_'$FileNum'",/' $ConfigFilename +sed -i 's/"file_name" \:.*/"file_name" \: "DEMPGen_'$EBeamE'on'$HBeamE'_'$InteractionPoint'_'$Ejectile$RecoilHadron'_'$NumEvents'_'$FileNum'",/' $ConfigFilename sed -i 's/"n_events" \:.*/"n_events" \: '$NumEvents',/' $ConfigFilename sed -i 's/"generator_seed"\:.*/"generator_seed" \: '$RandomSeed',/' $ConfigFilename sed -i 's/"ebeam"\:.*/"ebeam" \: '$EBeamE',/' $ConfigFilename sed -i 's/"hbeam"\:.*/"hbeam" \: '$HBeamE',/' $ConfigFilename -sed -i 's/"particle"\:.*/"particle" \: "'$Particle'",/' $ConfigFilename -sed -i 's/"hadron"\:.*/"hadron" \: "'$Hadron'",/' $ConfigFilename +sed -i 's/"ejectile"\:.*/"ejectile" \: "'$Ejectile'",/' $ConfigFilename +sed -i 's/"recoil_hadron"\:.*/"recoil_hadron" \: "'$RecoilHadron'",/' $ConfigFilename sed -i 's/"det_location"\:.*/"det_location" \: "'$InteractionPoint'",/' $ConfigFilename sed -i 's/"OutputType"\:.*/"OutputType"\: "'$OutputType'",/' $ConfigFilename @@ -60,8 +60,8 @@ cd data/ sleep 5 # Filename as it's created is a bit odd, so rename it -set OriginalOutput = 'eic_input_DEMPGen_'$EBeamE'on'$HBeamE'_'$InteractionPoint'_'$Particle$Hadron'_'$NumEvents'_'$FileNum'.dat' -set RenamedOutput = 'eic_DEMPGen_'$EBeamE'on'$HBeamE'_'$InteractionPoint'_'$Particle$Hadron'_'$NumEvents'_'$FileNum'.dat' +set OriginalOutput = 'eic_input_DEMPGen_'$EBeamE'on'$HBeamE'_'$InteractionPoint'_'$Ejectile$RecoilHadron'_'$NumEvents'_'$FileNum'.dat' +set RenamedOutput = 'eic_DEMPGen_'$EBeamE'on'$HBeamE'_'$InteractionPoint'_'$Ejectile$RecoilHadron'_'$NumEvents'_'$FileNum'.dat' mv "OutputFiles/"$OriginalOutput "OutputFiles/"$RenamedOutput rm -rf ../$ConfigFilename diff --git a/Process_EIC_iFarm.csh b/Process_EIC_iFarm.csh index 35cb6c3..bcea831 100755 --- a/Process_EIC_iFarm.csh +++ b/Process_EIC_iFarm.csh @@ -22,7 +22,7 @@ cd "$DEMPGENPath" echo"" echo "This file is intended to be run as part of a batch job submission, however, you can also run it on its own." -echo "Expected input is - FileNumber NumberOfEvents ElectronBeamEnergy HadronBeamEnergy RandomSeed OutputType IP Particle Hadron(Optional, for K+)" +echo "Expected input is - FileNumber NumberOfEvents ElectronBeamEnergy HadronBeamEnergy RandomSeed OutputType IP Ejectile Recoil_Hadron(Optional, for K+)" echo "Please see the README for more info." echo "" @@ -40,23 +40,23 @@ set HBeamE=$4 set RandomSeed=$5 set OutputType=$6 set InteractionPoint=$7 -set Particle=$8 +set Ejectile=$8 -if ($Particle == "K+" && $#argv == 8 ) then +if ($Ejectile == "K+" && $#argv == 8 ) then echo "! WARNING !" - echo "! WARNING ! - For K+ production expect a hadron specified, defaulting to Lambda - ! WARNING !" + echo "! WARNING ! - For K+ production expect a recoil hadron specified, defaulting to Lambda - ! WARNING !" echo "! WARNING !" - set Hadron="Lambda" -else if ($Particle == "K+" && $#argv == 9 ) then - set Hadron=$9 + set RecoilHadron="Lambda" +else if ($Ejectile == "K+" && $#argv == 9 ) then + set RecoilHadron=$9 else - set Hadron="" + set RecoilHadron="" endif -echo "Running target polarisation up, FF setting for file $FileNum with $NumEvents events per file for $EBeamE GeV e- on $HBeamE GeV p using random seed $RandomSeed, using $OutputType format output for $Particle $Hadron events." +echo "Running target polarisation up, FF setting for file $FileNum with $NumEvents events per file for $EBeamE GeV e- on $HBeamE GeV p using random seed $RandomSeed, using $OutputType format output for $Ejectile $RecoilHadron events." # Set the config file name based upon inputs -set ConfigFilename = $DEMPGENPath'/Config_EIC_'$EBeamE'on'$HBeamE'_'$InteractionPoint'_'$Particle$Hadron'_'$NumEvents'_'$FileNum'.json' +set ConfigFilename = $DEMPGENPath'/Config_EIC_'$EBeamE'on'$HBeamE'_'$InteractionPoint'_'$Ejectile$RecoilHadron'_'$NumEvents'_'$FileNum'.json' # Copy the default config file to our constructed filename cp "$DEMPGENPath/Config_EIC.json" $ConfigFilename @@ -67,8 +67,8 @@ sed -i 's/"n_events" \:.*/"n_events" \: '$NumEvents',/' $ConfigFilename sed -i 's/"generator_seed"\:.*/"generator_seed" \: '$RandomSeed',/' $ConfigFilename sed -i 's/"ebeam"\:.*/"ebeam" \: '$EBeamE',/' $ConfigFilename sed -i 's/"hbeam"\:.*/"hbeam" \: '$HBeamE',/' $ConfigFilename -sed -i 's/"particle"\:.*/"particle" \: "'$Particle'",/' $ConfigFilename -sed -i 's/"hadron"\:.*/"hadron" \: "'$Hadron'",/' $ConfigFilename +sed -i 's/"ejectile"\:.*/"ejectile" \: "'$Ejectile'",/' $ConfigFilename +sed -i 's/"recoil_hadron"\:.*/"recoil_hadron" \: "'$RecoilHadron'",/' $ConfigFilename sed -i 's/"det_location"\:.*/"det_location" \: "'$InteractionPoint'",/' $ConfigFilename sed -i 's/"OutputType"\:.*/"OutputType"\: "'$OutputType'",/' $ConfigFilename @@ -78,8 +78,8 @@ eval $DEMPGENPath/build/DEMPgen $ConfigFilename sleep 5 # Filename as it's created is a bit odd, so rename it -set OriginalOutput = $DEMPGENPath'/data/OutputFiles/eic_input_DEMPGen_'$EBeamE'on'$HBeamE'_'$InteractionPoint'_'$Particle$Hadron'_'$NumEvents'_'$FileNum'.dat' -set RenamedOutput = $DEMPGENPath'/data/OutputFiles/eic_DEMPGen_'$EBeamE'on'$HBeamE'_'$InteractionPoint'_'$Particle$Hadron'_'$NumEvents'_'$FileNum'.dat' +set OriginalOutput = $DEMPGENPath'/data/OutputFiles/eic_input_DEMPGen_'$EBeamE'on'$HBeamE'_'$InteractionPoint'_'$Ejectile$RecoilHadron'_'$NumEvents'_'$FileNum'.dat' +set RenamedOutput = $DEMPGENPath'/data/OutputFiles/eic_DEMPGen_'$EBeamE'on'$HBeamE'_'$InteractionPoint'_'$Ejectile$RecoilHadron'_'$NumEvents'_'$FileNum'.dat' mv $OriginalOutput $RenamedOutput diff --git a/README.md b/README.md index 6b94f9b..041c2c0 100644 --- a/README.md +++ b/README.md @@ -132,8 +132,8 @@ This script copies Config_EIC.json and formats a new file based upon this, DO NO - 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 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. + - Arg 8 - Ejectile -> The produced ejectile (meson) in the reaction, choose from omega, pi+, pi0 or K+. + - Arg 9 - RecoilHadron -> OPTIONAL - This only matters if you select K+ as the ejectile, in this case, choose from Lambda or Sigma0 here. If your choice is invalid (or you don't specify arg9), the default is Lambda. - So as an example if you executed the following - @@ -152,8 +152,8 @@ The jobs the script creates and submits all execute the Process_EIC.csh script d - 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 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 invalif. - - 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. + - Arg 7 - Ejectile -> The produced ejectile (meson) in the reaction, choose from omega, pi+, pi0 or K+. + - Arg 8 - RecoilHadron -> OPTIONAL - This only matters if you select K+ as the ejectile, in this case, choose from Lambda or Sigma0 here. If your choice is invalid (or you don't specify arg9), the default is Lambda. - The script automatically generates a random seed itself using the /dev/urandom function. diff --git a/src/eic_evgen/eic.cc b/src/eic_evgen/eic.cc index c5ea328..dbc6189 100644 --- a/src/eic_evgen/eic.cc +++ b/src/eic_evgen/eic.cc @@ -49,7 +49,7 @@ void eic() { /*--------------------------------------------------*/ // 18/01/23 - SJDK- This function is never used since eic() is only called with a json object as the argument. Commented out for now, delete later? -/* void eic(int event_number, int target_direction, int kinematics_type, TString file_name, int fEIC_seed, TString particle, TString hadron, TString det_location, TString OutputType, double EBeam, double HBeam) { +/* void eic(int event_number, int target_direction, int kinematics_type, TString file_name, int fEIC_seed, TString Ejectile, TString RecoilHadron, TString det_location, TString OutputType, double EBeam, double HBeam) { TString targetname; TString charge; @@ -65,45 +65,45 @@ void eic() { fNEvents = event_number; fSeed = fEIC_seed; - cout << EBeam << " elec " << HBeam << " hadrons" << endl; + cout << EBeam << " elec " << HBeam << " RecoilHadrons" << endl; fEBeam = EBeam; fPBeam = HBeam; pim* myPim = new pim(fSeed); myPim->Initilize(); - // 09/02/22 - SJDK - Special case for the kaon, if hadron not specified, default to Lambda - if (particle == "K+"){ - if (hadron != "Lambda" && hadron != "Sigma0"){ - hadron = "Lambda"; + // 09/02/22 - SJDK - Special case for the kaon, if RecoilHadron not specified, default to Lambda + if (Ejectile == "K+"){ + if (RecoilHadron != "Lambda" && RecoilHadron != "Sigma0"){ + RecoilHadron = "Lambda"; } else{ - hadron = ExtractParticle(hadron); + RecoilHadron = ExtractParticle(RecoilHadron); } - Reaction* r1 = new Reaction(particle, hadron); + Reaction* r1 = new Reaction(Ejectile, RecoilHadron); r1->process_reaction(); delete r1; } - else if (particle == "pi+" || particle == "Pion+" || particle == "Pi+"){ - hadron = "Neutron"; - particle = ExtractParticle(particle); - charge = ExtractCharge(particle); - Reaction* r1 = new Reaction(particle, hadron); + else if (Ejectile == "pi+" || Ejectile == "Pion+" || Ejectile == "Pi+"){ + RecoilHadron = "Neutron"; + Ejectile = ExtractParticle(particle); + charge = ExtractCharge(Ejectile); + Reaction* r1 = new Reaction(Ejectile, RecoilHadron); r1->process_reaction(); delete r1; } - else if (particle == "pi0" || particle == "Pion0" || particle == "Pi0"){ - hadron = "Proton"; - particle = ExtractParticle(particle); - charge = ExtractCharge(particle); - //Reaction* r1 = new Reaction(particle); - Reaction* r1 = new Reaction(particle, hadron); + else if (Ejectile == "pi0" || Ejectile == "Pion0" || Ejectile == "Pi0"){ + RecoilHadron = "Proton"; + Ejectile = ExtractParticle(Ejectile); + charge = ExtractCharge(Ejectile); + //Reaction* r1 = new Reaction(Ejectile); + Reaction* r1 = new Reaction(Ejectile, RecoilHadron); r1->process_reaction(); delete r1; } else{ - particle = ExtractParticle(particle); - charge = ExtractCharge(particle); - Reaction* r1 = new Reaction(particle); + Ejectile = ExtractParticle(Ejectile); + charge = ExtractCharge(Ejectile); + Reaction* r1 = new Reaction(Ejectile); r1->process_reaction(); delete r1; } @@ -139,54 +139,54 @@ void eic(Json::Value obj) { // TDatime dsTime; // cout << "Start Time: " << dsTime.GetHour() << ":" << dsTime.GetMinute() << endl; // 21/12/22 - SJDK - Should do a check if these are defined or not, should crash if not defined or set defaults, see other quantities below - TString particle = obj["particle"].asString(); - TString hadron = obj["hadron"].asString(); // 09/02/22 - SJDK - Added in hadron type argument for K+ - // SJDK - 08/02/22 - This is terrible, need to change this, particle should just be K+ - // Add a new flag which, hadron - where this is specified too, then add conditionals elsewhere based on this + TString Ejectile = obj["ejectile"].asString(); + TString RecoilHadron = obj["recoil_hadron"].asString(); // 09/02/22 - SJDK - Added in RecoilHadron type argument for K+ + // SJDK - 08/02/22 - This is terrible, need to change this, Ejectile should just be K+ + // Add a new flag which, RecoilHadron - where this is specified too, then add conditionals elsewhere based on this // New conditional, special case for Kaon - particle = ExtractParticle(particle); - charge = ExtractCharge(particle); - if(hadron == "Sigma" || hadron == "sigma"){ // SJDK - 31/01/23 - If hadron specified as Sigma, interpret this as Sigma0. Also correct for lower case - hadron = "Sigma0"; + Ejectile = ExtractParticle(Ejectile); + charge = ExtractCharge(Ejectile); + if(RecoilHadron == "Sigma" || RecoilHadron == "sigma"){ // SJDK - 31/01/23 - If RecoilHadron specified as Sigma, interpret this as Sigma0. Also correct for lower case + RecoilHadron = "Sigma0"; } - if (hadron == "lambda"){ // SJDK - 31/01/23 - Make Lambda selection case insensitive - hadron = "Lambda"; + if (RecoilHadron == "lambda"){ // SJDK - 31/01/23 - Make Lambda selection case insensitive + RecoilHadron = "Lambda"; } - if (particle == "K+"){ - if (hadron != "Lambda" && hadron != "Sigma0"){ - hadron = "Lambda"; + if (Ejectile == "K+"){ + if (RecoilHadron != "Lambda" && RecoilHadron != "Sigma0"){ + RecoilHadron = "Lambda"; cout << "! WARNING !" << endl; - cout << "! WARNING !- K+ production specified but hadron not recognised, deaulting to Lambda - ! WARNING!" << endl; + cout << "! WARNING !- K+ production specified but RecoilHadron not recognised, deaulting to Lambda - ! WARNING!" << endl; cout << "! WARNING !" << endl; } else{ - hadron = ExtractParticle(hadron); + RecoilHadron = ExtractParticle(RecoilHadron); } } - // SJDK - 19/12/22 - Specify hadron to neutron/proton for pi+/pi0 production, for pi0 production, may want to adjust? - else if (particle == "pi+" || particle == "Pion+" || particle == "Pi+"){ - hadron = "Neutron"; + // SJDK - 19/12/22 - Specify RecoilHadron to neutron/proton for pi+/pi0 production, for pi0 production, may want to adjust? + else if (Ejectile == "pi+" || Ejectile == "Pion+" || Ejectile == "Pi+"){ + RecoilHadron = "Neutron"; } - else if (particle == "pi0" || particle == "Pion0" || particle == "Pi0"){ - hadron = "Proton"; + else if (Ejectile == "pi0" || Ejectile == "Pion0" || Ejectile == "Pi0"){ + RecoilHadron = "Proton"; } - else { // SJDK -09/02/22 - Note that in future this could be changed to get different hadrons in other reactions if desired - hadron = ""; + else { // SJDK -09/02/22 - Note that in future this could be changed to get different RecoilHadrons in other reactions if desired + RecoilHadron = ""; } // SJDK 03/04/23 - Change to how Qsq range is set/chosen, could add as an override variable later too - // Set min/max Qsq values depending upon particle type - if (particle == "pi+" || particle == "Pion+" || particle == "Pi+"){ + // Set min/max Qsq values depending upon Ejectile type + if (Ejectile == "pi+" || Ejectile == "Pion+" || Ejectile == "Pi+"){ fQsq_Min = 3.0; fQsq_Max = 35.0; fW_Min = 2.0; fW_Max = 10.2; fT_Max = 1.3; } - else if (particle == "pi0" || particle == "Pion0" || particle == "Pi0"){ + else if (Ejectile == "pi0" || Ejectile == "Pion0" || Ejectile == "Pi0"){ fQsq_Min = 5.0; fQsq_Max = 1000.0; fW_Min = 2.0; fW_Max = 10.0; fT_Max = 0.5; } - else if (particle == "K+"){ + else if (Ejectile == "K+"){ fQsq_Min = 1.0; fQsq_Max = 35.0; fW_Min = 2.0; fW_Max = 10.0; fT_Max = 2.0; @@ -330,15 +330,15 @@ void eic(Json::Value obj) { UseSolve = false; } - SigPar = ReadCrossSectionPar(particle, hadron); + SigPar = ReadCrossSectionPar(Ejectile, RecoilHadron); - if(particle != "pi0"){ // Default case now - Reaction* r1 = new Reaction(particle, hadron); + if(Ejectile != "pi0"){ // Default case now + Reaction* r1 = new Reaction(Ejectile, RecoilHadron); r1->process_reaction(); delete r1; } else{ // Treat pi0 slightly differently for now - Reaction* r1 = new Reaction(particle); + Reaction* r1 = new Reaction(Ejectile); r1->process_reaction(); delete r1; } @@ -402,29 +402,29 @@ TString ExtractCharge(TString particle) { return charge; } -vector>>> ReadCrossSectionPar(TString particle, TString hadron){ +vector>>> ReadCrossSectionPar(TString EjectileX, TString RecoilHad){ string sigL_ParamFile, sigT_ParamFile; - if (particle == "Pi+" && hadron == "Neutron"){ + if (EjectileX == "Pi+" && RecoilHad == "Neutron"){ // When pion model parameterised in some way, add Pi+/Neutron case here - } - else if (particle == "Pi-" && hadron == "Proton"){ + else if (EjectileX == "Pi-" && RecoilHad == "Proton"){ // When pion model parameterised in some way, add Pi-/Proton case here - } - else if (particle == "K+" && hadron == "Lambda"){ + else if (EjectileX == "K+" && RecoilHad == "Lambda"){ sigL_ParamFile = "../src/eic_evgen/CrossSection_Params/KPlusLambda_Param_sigL"; sigT_ParamFile = "../src/eic_evgen/CrossSection_Params/KPlusLambda_Param_sigT"; // Shouldn't really have a relative path, should look at setting a DEMPGen variable and doing this in a better way later } - else if (particle == "K+" && hadron == "Sigma0"){ + else if (EjectileX == "K+" && RecoilHad == "Sigma0"){ sigL_ParamFile = "../src/eic_evgen/CrossSection_Params/KPlusSigma_Param_sigL"; sigT_ParamFile = "../src/eic_evgen/CrossSection_Params/KPlusSigma_Param_sigT"; } - else if (particle == "Pi0"){ + else if (EjectileX == "Pi0"){ // When pi0 model parameterised, add it here } else{ - cerr << " !!!!! " << endl << "Warning!" << endl << "Combination of specified ejectile and recoil particles not found!" << "Cross section parameters cannot be read, check inputs!" << endl << "Warning!" << endl << " !!!!! " << endl; + cerr << " !!!!! " << endl << "Warning!" << endl << "Combination of specified ejectile and recoil hadron not found!" << "Cross section parameters cannot be read, check inputs!" << endl << "Warning!" << endl << " !!!!! " << endl; exit(-1); } diff --git a/src/eic_evgen/eic.h b/src/eic_evgen/eic.h index 81b036d..c4847bf 100644 --- a/src/eic_evgen/eic.h +++ b/src/eic_evgen/eic.h @@ -58,7 +58,7 @@ void SetEICSeed(int); TString ExtractParticle(TString); TString ExtractCharge(TString); -vector>>> ReadCrossSectionPar(TString particle, TString hadron); +vector>>> ReadCrossSectionPar(TString EjectileX, TString RecoilHad); #endif diff --git a/src/eic_evgen/reaction_routine.cc b/src/eic_evgen/reaction_routine.cc index 6112507..0bec8e6 100644 --- a/src/eic_evgen/reaction_routine.cc +++ b/src/eic_evgen/reaction_routine.cc @@ -19,12 +19,11 @@ using namespace std; /*--------------------------------------------------*/ /// Reaction +Reaction::Reaction(TString ejectile_str) { -Reaction::Reaction(TString particle_str) { - - rEjectile = particle_str; - cout << "Produced particle is: " << GetParticle() << endl; - cout << "Generated process: e + p -> e' + p' + " << GetParticle() << endl; + rEjectile = ejectile_str; + cout << "Produced ejectile is: " << GetEjectile() << endl; + cout << "Generated process: e + p -> e' + p' + " << GetEjectile() << endl; tTime.Start(); cout << "/*--------------------------------------------------*/" << endl; @@ -37,13 +36,13 @@ Reaction::Reaction(TString particle_str) { } // SJDK 09/02/22 - New reaction where the particle and hadron are specified -Reaction::Reaction(TString particle_str, TString hadron_str) { +Reaction::Reaction(TString ejectile_str, TString recoil_hadron_str) { - rEjectile = particle_str; - rRecoil = hadron_str; - cout << "Produced particle is: " << GetParticle() << endl; - cout << "Produced hadron is: " << GetHadron() << endl; - cout << "Generated process: e + p -> e'+ " << GetHadron() << " + " << GetParticle() << endl; + rEjectile = ejectile_str; + rRecoil = recoil_hadron_str; + cout << "Produced ejectile is: " << GetEjectile() << endl; + cout << "Produced recoil hadron is: " << GetRecoilHadron() << endl; + cout << "Generated process: e + p -> e'+ " << GetRecoilHadron() << " + " << GetEjectile() << endl; tTime.Start(); cout << "/*--------------------------------------------------*/" << endl; diff --git a/src/eic_evgen/reaction_routine.h b/src/eic_evgen/reaction_routine.h index 42ebd59..8b3c7b0 100644 --- a/src/eic_evgen/reaction_routine.h +++ b/src/eic_evgen/reaction_routine.h @@ -29,8 +29,8 @@ class Reaction{ ~Reaction(); void process_reaction(); - TString GetParticle() {return rEjectile;}; - TString GetHadron() {return rRecoil;}; + TString GetEjectile() {return rEjectile;}; + TString GetRecoilHadron() {return rRecoil;}; protected: TStopwatch tTime; @@ -48,8 +48,8 @@ class DEMP_Reaction { ~DEMP_Reaction(); void process_reaction(); - TString GetParticle() {return rEjectile;}; - TString GetHadron() {return rRecoil;}; + TString GetEjectile() {return rEjectile;}; + TString GetRecoilHadron() {return rRecoil;}; protected: From c66a6261a3500bcf5b7d1d541c169604b6b53d59 Mon Sep 17 00:00:00 2001 From: Stephen Kay Date: Wed, 13 Sep 2023 08:16:12 -0600 Subject: [PATCH 20/21] New generic fHBeam double rather than fPBeam. In place in DEMP_reaction.cc, no changes to output. --- src/eic_evgen/eic.cc | 576 +++++++++--------- src/eic_evgen/eic_pim.cc | 3 +- src/eic_evgen/eic_pim.h | 2 + .../process_routine/DEMP_Reaction.cc | 16 +- 4 files changed, 301 insertions(+), 296 deletions(-) diff --git a/src/eic_evgen/eic.cc b/src/eic_evgen/eic.cc index dbc6189..f487d23 100644 --- a/src/eic_evgen/eic.cc +++ b/src/eic_evgen/eic.cc @@ -33,17 +33,17 @@ using namespace std; void eic() { - Int_t target_direction, kinematics_type; - Double_t EBeam, HBeam; + Int_t target_direction, kinematics_type; + Double_t EBeam, HBeam; - cout << "Target Direction (1->Up, 2->Down): "; cin >> target_direction; cout << endl; - cout << "Kinematics type (1->FF, 2->TSSA): "; cin >> kinematics_type; cout << endl; - cout << "Enter the number of events: "; cin >> fNEvents; cout << endl; - cout << "Enter the file number: "; cin >> fNFile; cout << endl; - cout << "Enter the electron beam energy: "; cin >> EBeam; cout << endl; - cout << "Enter the hadron beam energy: "; cin >> HBeam; cout << endl; + cout << "Target Direction (1->Up, 2->Down): "; cin >> target_direction; cout << endl; + cout << "Kinematics type (1->FF, 2->TSSA): "; cin >> kinematics_type; cout << endl; + cout << "Enter the number of events: "; cin >> fNEvents; cout << endl; + cout << "Enter the file number: "; cin >> fNFile; cout << endl; + cout << "Enter the electron beam energy: "; cin >> EBeam; cout << endl; + cout << "Enter the hadron beam energy: "; cin >> HBeam; cout << endl; -// eic(target_direction, kinematics_type, fNEvents); + // eic(target_direction, kinematics_type, fNEvents); } @@ -51,63 +51,63 @@ void eic() { // 18/01/23 - SJDK- This function is never used since eic() is only called with a json object as the argument. Commented out for now, delete later? /* void eic(int event_number, int target_direction, int kinematics_type, TString file_name, int fEIC_seed, TString Ejectile, TString RecoilHadron, TString det_location, TString OutputType, double EBeam, double HBeam) { - TString targetname; - TString charge; + TString targetname; + TString charge; - if( target_direction == 1 ) targetname = "up"; - if( target_direction == 2 ) targetname = "down"; + if( target_direction == 1 ) targetname = "up"; + if( target_direction == 2 ) targetname = "down"; - gKinematics_type = kinematics_type; - gfile_name = file_name; - - fNFile = 1; - - fNEvents = event_number; - - fSeed = fEIC_seed; - cout << EBeam << " elec " << HBeam << " RecoilHadrons" << endl; - fEBeam = EBeam; - fPBeam = HBeam; - - pim* myPim = new pim(fSeed); - myPim->Initilize(); - // 09/02/22 - SJDK - Special case for the kaon, if RecoilHadron not specified, default to Lambda - if (Ejectile == "K+"){ - if (RecoilHadron != "Lambda" && RecoilHadron != "Sigma0"){ - RecoilHadron = "Lambda"; - } - else{ - RecoilHadron = ExtractParticle(RecoilHadron); - } - Reaction* r1 = new Reaction(Ejectile, RecoilHadron); - r1->process_reaction(); - delete r1; - } - else if (Ejectile == "pi+" || Ejectile == "Pion+" || Ejectile == "Pi+"){ - RecoilHadron = "Neutron"; - Ejectile = ExtractParticle(particle); - charge = ExtractCharge(Ejectile); - Reaction* r1 = new Reaction(Ejectile, RecoilHadron); - r1->process_reaction(); - delete r1; - } - else if (Ejectile == "pi0" || Ejectile == "Pion0" || Ejectile == "Pi0"){ - RecoilHadron = "Proton"; - Ejectile = ExtractParticle(Ejectile); - charge = ExtractCharge(Ejectile); - //Reaction* r1 = new Reaction(Ejectile); - Reaction* r1 = new Reaction(Ejectile, RecoilHadron); - r1->process_reaction(); - delete r1; - } - else{ - Ejectile = ExtractParticle(Ejectile); - charge = ExtractCharge(Ejectile); - Reaction* r1 = new Reaction(Ejectile); - r1->process_reaction(); - delete r1; - } -} + gKinematics_type = kinematics_type; + gfile_name = file_name; + + fNFile = 1; + + fNEvents = event_number; + + fSeed = fEIC_seed; + cout << EBeam << " elec " << HBeam << " RecoilHadrons" << endl; + fEBeam = EBeam; + fPBeam = HBeam; + + pim* myPim = new pim(fSeed); + myPim->Initilize(); + // 09/02/22 - SJDK - Special case for the kaon, if RecoilHadron not specified, default to Lambda + if (Ejectile == "K+"){ + if (RecoilHadron != "Lambda" && RecoilHadron != "Sigma0"){ + RecoilHadron = "Lambda"; + } + else{ + RecoilHadron = ExtractParticle(RecoilHadron); + } + Reaction* r1 = new Reaction(Ejectile, RecoilHadron); + r1->process_reaction(); + delete r1; + } + else if (Ejectile == "pi+" || Ejectile == "Pion+" || Ejectile == "Pi+"){ + RecoilHadron = "Neutron"; + Ejectile = ExtractParticle(particle); + charge = ExtractCharge(Ejectile); + Reaction* r1 = new Reaction(Ejectile, RecoilHadron); + r1->process_reaction(); + delete r1; + } + else if (Ejectile == "pi0" || Ejectile == "Pion0" || Ejectile == "Pi0"){ + RecoilHadron = "Proton"; + Ejectile = ExtractParticle(Ejectile); + charge = ExtractCharge(Ejectile); + //Reaction* r1 = new Reaction(Ejectile); + Reaction* r1 = new Reaction(Ejectile, RecoilHadron); + r1->process_reaction(); + delete r1; + } + else{ + Ejectile = ExtractParticle(Ejectile); + charge = ExtractCharge(Ejectile); + Reaction* r1 = new Reaction(Ejectile); + r1->process_reaction(); + delete r1; + } + } */ /*--------------------------------------------------*/ @@ -115,240 +115,242 @@ void eic() { // SJDK 21/12/22 - Note that this is the one that actually gets used, reads in the .json file void eic(Json::Value obj) { - TString targetname; - TString charge; + TString targetname; + TString charge; - int target_direction = obj["Targ_dir"].asInt(); - gKinematics_type = obj["Kinematics_type"].asInt(); + int target_direction = obj["Targ_dir"].asInt(); + gKinematics_type = obj["Kinematics_type"].asInt(); - if( target_direction == 1 ) targetname = "up"; - if( target_direction == 2 ) targetname = "down"; + if( target_direction == 1 ) targetname = "up"; + if( target_direction == 2 ) targetname = "down"; - gfile_name = obj["file_name"].asString(); + gfile_name = obj["file_name"].asString(); - gPi0_decay = obj["pi0_decay"].asBool(); + gPi0_decay = obj["pi0_decay"].asBool(); - fNFile = 1; - fNEvents = obj["n_events"].asUInt64(); + fNFile = 1; + fNEvents = obj["n_events"].asUInt64(); - fSeed = obj["generator_seed"].asInt(); + fSeed = obj["generator_seed"].asInt(); - pim* myPim = new pim(fSeed); - myPim->Initilize(); + pim* myPim = new pim(fSeed); + myPim->Initilize(); -// TDatime dsTime; -// cout << "Start Time: " << dsTime.GetHour() << ":" << dsTime.GetMinute() << endl; - // 21/12/22 - SJDK - Should do a check if these are defined or not, should crash if not defined or set defaults, see other quantities below - TString Ejectile = obj["ejectile"].asString(); - TString RecoilHadron = obj["recoil_hadron"].asString(); // 09/02/22 - SJDK - Added in RecoilHadron type argument for K+ - // SJDK - 08/02/22 - This is terrible, need to change this, Ejectile should just be K+ - // Add a new flag which, RecoilHadron - where this is specified too, then add conditionals elsewhere based on this - // New conditional, special case for Kaon - Ejectile = ExtractParticle(Ejectile); - charge = ExtractCharge(Ejectile); - if(RecoilHadron == "Sigma" || RecoilHadron == "sigma"){ // SJDK - 31/01/23 - If RecoilHadron specified as Sigma, interpret this as Sigma0. Also correct for lower case - RecoilHadron = "Sigma0"; - } - if (RecoilHadron == "lambda"){ // SJDK - 31/01/23 - Make Lambda selection case insensitive - RecoilHadron = "Lambda"; - } - if (Ejectile == "K+"){ - if (RecoilHadron != "Lambda" && RecoilHadron != "Sigma0"){ - RecoilHadron = "Lambda"; - cout << "! WARNING !" << endl; - cout << "! WARNING !- K+ production specified but RecoilHadron not recognised, deaulting to Lambda - ! WARNING!" << endl; - cout << "! WARNING !" << endl; - } - else{ - RecoilHadron = ExtractParticle(RecoilHadron); - } - } - // SJDK - 19/12/22 - Specify RecoilHadron to neutron/proton for pi+/pi0 production, for pi0 production, may want to adjust? - else if (Ejectile == "pi+" || Ejectile == "Pion+" || Ejectile == "Pi+"){ - RecoilHadron = "Neutron"; - } - else if (Ejectile == "pi0" || Ejectile == "Pion0" || Ejectile == "Pi0"){ - RecoilHadron = "Proton"; - } - else { // SJDK -09/02/22 - Note that in future this could be changed to get different RecoilHadrons in other reactions if desired - RecoilHadron = ""; - } + // TDatime dsTime; + // cout << "Start Time: " << dsTime.GetHour() << ":" << dsTime.GetMinute() << endl; + // 21/12/22 - SJDK - Should do a check if these are defined or not, should crash if not defined or set defaults, see other quantities below + TString Ejectile = obj["ejectile"].asString(); + TString RecoilHadron = obj["recoil_hadron"].asString(); // 09/02/22 - SJDK - Added in RecoilHadron type argument for K+ + // SJDK - 08/02/22 - This is terrible, need to change this, Ejectile should just be K+ + // Add a new flag which, RecoilHadron - where this is specified too, then add conditionals elsewhere based on this + // New conditional, special case for Kaon + Ejectile = ExtractParticle(Ejectile); + charge = ExtractCharge(Ejectile); + if(RecoilHadron == "Sigma" || RecoilHadron == "sigma"){ // SJDK - 31/01/23 - If RecoilHadron specified as Sigma, interpret this as Sigma0. Also correct for lower case + RecoilHadron = "Sigma0"; + } + if (RecoilHadron == "lambda"){ // SJDK - 31/01/23 - Make Lambda selection case insensitive + RecoilHadron = "Lambda"; + } + if (Ejectile == "K+"){ + if (RecoilHadron != "Lambda" && RecoilHadron != "Sigma0"){ + RecoilHadron = "Lambda"; + cout << "! WARNING !" << endl; + cout << "! WARNING !- K+ production specified but RecoilHadron not recognised, deaulting to Lambda - ! WARNING!" << endl; + cout << "! WARNING !" << endl; + } + else{ + RecoilHadron = ExtractParticle(RecoilHadron); + } + } + // SJDK - 19/12/22 - Specify RecoilHadron to neutron/proton for pi+/pi0 production, for pi0 production, may want to adjust? + else if (Ejectile == "pi+" || Ejectile == "Pion+" || Ejectile == "Pi+"){ + RecoilHadron = "Neutron"; + } + else if (Ejectile == "pi0" || Ejectile == "Pion0" || Ejectile == "Pi0"){ + RecoilHadron = "Proton"; + } + else { // SJDK -09/02/22 - Note that in future this could be changed to get different RecoilHadrons in other reactions if desired + RecoilHadron = ""; + } - // SJDK 03/04/23 - Change to how Qsq range is set/chosen, could add as an override variable later too - // Set min/max Qsq values depending upon Ejectile type - if (Ejectile == "pi+" || Ejectile == "Pion+" || Ejectile == "Pi+"){ - fQsq_Min = 3.0; fQsq_Max = 35.0; - fW_Min = 2.0; fW_Max = 10.2; - fT_Max = 1.3; - } - else if (Ejectile == "pi0" || Ejectile == "Pion0" || Ejectile == "Pi0"){ - fQsq_Min = 5.0; fQsq_Max = 1000.0; - fW_Min = 2.0; fW_Max = 10.0; - fT_Max = 0.5; - } - else if (Ejectile == "K+"){ - fQsq_Min = 1.0; fQsq_Max = 35.0; - fW_Min = 2.0; fW_Max = 10.0; - fT_Max = 2.0; - } - else{ - fQsq_Min = 5.0; fQsq_Max = 35.0; - fW_Min = 2.0; fW_Max = 10.0; - fT_Max = 1.3; - } + // SJDK 03/04/23 - Change to how Qsq range is set/chosen, could add as an override variable later too + // Set min/max Qsq values depending upon Ejectile type + if (Ejectile == "pi+" || Ejectile == "Pion+" || Ejectile == "Pi+"){ + fQsq_Min = 3.0; fQsq_Max = 35.0; + fW_Min = 2.0; fW_Max = 10.2; + fT_Max = 1.3; + } + else if (Ejectile == "pi0" || Ejectile == "Pion0" || Ejectile == "Pi0"){ + fQsq_Min = 5.0; fQsq_Max = 1000.0; + fW_Min = 2.0; fW_Max = 10.0; + fT_Max = 0.5; + } + else if (Ejectile == "K+"){ + fQsq_Min = 1.0; fQsq_Max = 35.0; + fW_Min = 2.0; fW_Max = 10.0; + fT_Max = 2.0; + } + else{ + fQsq_Min = 5.0; fQsq_Max = 35.0; + fW_Min = 2.0; fW_Max = 10.0; + fT_Max = 1.3; + } - // SJDK - 01/06/21 - // Set beam energies from .json read in - if (obj.isMember("ebeam")){ - fEBeam = obj["ebeam"].asDouble(); - } - else{ - fEBeam = 5; - cout << "Electron beam energy not specified in .json file, defaulting to 5 GeV." << endl; - } - if (obj.isMember("hbeam")){ - fPBeam = obj["hbeam"].asDouble(); - } - else{ - fPBeam = 100; - cout << "Ion beam energy not specified in .json file, defaulting to 100 GeV." << endl; - } + // SJDK - 01/06/21 + // Set beam energies from .json read in + if (obj.isMember("ebeam")){ + fEBeam = obj["ebeam"].asDouble(); + } + else{ + fEBeam = 5; + cout << "Electron beam energy not specified in .json file, defaulting to 5 GeV." << endl; + } + if (obj.isMember("hbeam")){ + fPBeam = obj["hbeam"].asDouble(); + fHBeam = obj["hbeam"].asDouble(); + } + else{ + fPBeam = 100; + fHBeam = 100; + cout << "Ion beam energy not specified in .json file, defaulting to 100 GeV." << endl; + } - if (obj.isMember("hbeam_part")){ - gBeamPart = obj["hbeam_part"].asString(); - } + if (obj.isMember("hbeam_part")){ + gBeamPart = obj["hbeam_part"].asString(); + } + else { + gBeamPart = "Proton"; + } + + // SJDK - 12/01/22 + // Set output type as a .json read in + // Should be Pythia6, LUND or HEPMC3 + if (obj.isMember("OutputType")){ + gOutputType = obj["OutputType"].asString(); + if (gOutputType == "Pythia6"){ + cout << "Using Pythia6 output format for Fun4All" << endl; + } + else if (gOutputType == "LUND"){ + cout << "Using LUND output format" << endl; + } + else if (gOutputType == "HEPMC3"){ + cout << "Using HEPMC3 output format for ePIC" << endl; + } + else{ + cout << "Output type not recognised!" << endl; + cout << "Setting output type to HEPMC3 (ePIC) by default!" << endl; + gOutputType = "HEPMC3"; + } + } + else{ + cout << "Output type not specified in .json file!" << endl; + cout << "Setting output type to HEPMC3 (ePIC) by default!" << endl; + gOutputType = "HEPMC3"; + } + ///*--------------------------------------------------*/ + /// The detector selection is determined here + /// The incidence proton phi angle is + if (obj.isMember("det_location")){ + gDet_location = obj["det_location"].asString(); + if (gDet_location == "ip8") { + fProton_incidence_phi = 0.0; + } + else if (gDet_location == "ip6") { + fProton_incidence_phi = fPi; + } else { - gBeamPart = "Proton"; + fProton_incidence_phi = 0.0; + cout << "The interaction point requested is not recognized!" << endl; + cout << "Therefore ip6 is used by default." << endl; } + } + else{ // 21/12/22 - This could probably be combined with the else statement above in some way + fProton_incidence_phi = 0.0; + cout << "The interaction points was not specified in the .json file!" << endl; + cout << "Therefore ip6 is used by default" << endl; + } - // SJDK - 12/01/22 - // Set output type as a .json read in - // Should be Pythia6, LUND or HEPMC3 - if (obj.isMember("OutputType")){ - gOutputType = obj["OutputType"].asString(); - if (gOutputType == "Pythia6"){ - cout << "Using Pythia6 output format for Fun4All" << endl; - } - else if (gOutputType == "LUND"){ - cout << "Using LUND output format" << endl; - } - else if (gOutputType == "HEPMC3"){ - cout << "Using HEPMC3 output format for ePIC" << endl; - } - else{ - cout << "Output type not recognised!" << endl; - cout << "Setting output type to HEPMC3 (ePIC) by default!" << endl; - gOutputType = "HEPMC3"; - } - } - else{ - cout << "Output type not specified in .json file!" << endl; - cout << "Setting output type to HEPMC3 (ePIC) by default!" << endl; - gOutputType = "HEPMC3"; - } - ///*--------------------------------------------------*/ - /// The detector selection is determined here - /// The incidence proton phi angle is - if (obj.isMember("det_location")){ - gDet_location = obj["det_location"].asString(); - if (gDet_location == "ip8") { - fProton_incidence_phi = 0.0; - } - else if (gDet_location == "ip6") { - fProton_incidence_phi = fPi; - } - else { - fProton_incidence_phi = 0.0; - cout << "The interaction point requested is not recognized!" << endl; - cout << "Therefore ip6 is used by default." << endl; - } - } - else{ // 21/12/22 - This could probably be combined with the else statement above in some way - fProton_incidence_phi = 0.0; - cout << "The interaction points was not specified in the .json file!" << endl; - cout << "Therefore ip6 is used by default" << endl; - } - - if (obj.isMember("Ee_Low")){ - fScatElec_E_Lo = obj["Ee_Low"].asDouble(); - } - else{ - fScatElec_E_Lo = 0.5; - cout << "Minumum scattered electron energy not specified in .json file, defaulting to 0.5*EBeam." << endl; - } + if (obj.isMember("Ee_Low")){ + fScatElec_E_Lo = obj["Ee_Low"].asDouble(); + } + else{ + fScatElec_E_Lo = 0.5; + cout << "Minumum scattered electron energy not specified in .json file, defaulting to 0.5*EBeam." << endl; + } - if (obj.isMember("Ee_High")){ - fScatElec_E_Hi = obj["Ee_High"].asDouble(); - } - else{ - fScatElec_E_Hi = 2.5; - cout << "Max scattered electron energy not specified in .json file, defaulting to 2.5*EBeam." << endl; - } + if (obj.isMember("Ee_High")){ + fScatElec_E_Hi = obj["Ee_High"].asDouble(); + } + else{ + fScatElec_E_Hi = 2.5; + cout << "Max scattered electron energy not specified in .json file, defaulting to 2.5*EBeam." << endl; + } - if (obj.isMember("e_Theta_Low")){ - fScatElec_Theta_I = obj["e_Theta_Low"].asDouble() * fDEG2RAD; - } - else{ - fScatElec_Theta_I = 60.0 * fDEG2RAD; - cout << "Min scattered electron theta not specified in .json file, defaulting to 60 degrees." << endl; - } + if (obj.isMember("e_Theta_Low")){ + fScatElec_Theta_I = obj["e_Theta_Low"].asDouble() * fDEG2RAD; + } + else{ + fScatElec_Theta_I = 60.0 * fDEG2RAD; + cout << "Min scattered electron theta not specified in .json file, defaulting to 60 degrees." << endl; + } - if (obj.isMember("e_Theta_High")){ - fScatElec_Theta_F = obj["e_Theta_High"].asDouble() * fDEG2RAD; - } - else{ - fScatElec_Theta_F = 175.0 * fDEG2RAD; - cout << "Max scattered electron theta not specified in .json file, defaulting to 175 degrees." << endl; - } + if (obj.isMember("e_Theta_High")){ + fScatElec_Theta_F = obj["e_Theta_High"].asDouble() * fDEG2RAD; + } + else{ + fScatElec_Theta_F = 175.0 * fDEG2RAD; + cout << "Max scattered electron theta not specified in .json file, defaulting to 175 degrees." << endl; + } - if (obj.isMember("EjectileX_Theta_Low")){ - fEjectileX_Theta_I = obj["EjectileX_Theta_Low"].asDouble() * fDEG2RAD; - } - else{ - fEjectileX_Theta_I = 0.0 * fDEG2RAD; - cout << "Min ejectile X theta not specified in .json file, defaulting to 0 degrees." << endl; - } + if (obj.isMember("EjectileX_Theta_Low")){ + fEjectileX_Theta_I = obj["EjectileX_Theta_Low"].asDouble() * fDEG2RAD; + } + else{ + fEjectileX_Theta_I = 0.0 * fDEG2RAD; + cout << "Min ejectile X theta not specified in .json file, defaulting to 0 degrees." << endl; + } - if (obj.isMember("EjectileX_Theta_High")){ - fEjectileX_Theta_F = obj["EjectileX_Theta_High"].asDouble() * fDEG2RAD; - } - else{ - fEjectileX_Theta_F = 60.0 * fDEG2RAD; - cout << "Max ejectile X theta not specified in .json file, defaulting to 60 degrees." << endl; - } + if (obj.isMember("EjectileX_Theta_High")){ + fEjectileX_Theta_F = obj["EjectileX_Theta_High"].asDouble() * fDEG2RAD; + } + else{ + fEjectileX_Theta_F = 60.0 * fDEG2RAD; + cout << "Max ejectile X theta not specified in .json file, defaulting to 60 degrees." << endl; + } - // 06/09/23 - SJDK - Added string to check method chosen - TString CalcMethod = obj["calc_method"].asString(); - if(CalcMethod == "Analytical"){ - UseSolve = false; - } - else if (CalcMethod == "Solve"){ - UseSolve = true; - } - else{ - cout << "! WARNING !" << endl << "! WARNING !- Calculation method not specified or not recognised, defaulting to Analytical - ! WARNING!" << endl << "! WARNING !" << endl; - UseSolve = false; - } + // 06/09/23 - SJDK - Added string to check method chosen + TString CalcMethod = obj["calc_method"].asString(); + if(CalcMethod == "Analytical"){ + UseSolve = false; + } + else if (CalcMethod == "Solve"){ + UseSolve = true; + } + else{ + cout << "! WARNING !" << endl << "! WARNING !- Calculation method not specified or not recognised, defaulting to Analytical - ! WARNING!" << endl << "! WARNING !" << endl; + UseSolve = false; + } - SigPar = ReadCrossSectionPar(Ejectile, RecoilHadron); + SigPar = ReadCrossSectionPar(Ejectile, RecoilHadron); - if(Ejectile != "pi0"){ // Default case now - Reaction* r1 = new Reaction(Ejectile, RecoilHadron); - r1->process_reaction(); - delete r1; - } - else{ // Treat pi0 slightly differently for now - Reaction* r1 = new Reaction(Ejectile); - r1->process_reaction(); - delete r1; - } + if(Ejectile != "pi0"){ // Default case now + Reaction* r1 = new Reaction(Ejectile, RecoilHadron); + r1->process_reaction(); + delete r1; + } + else{ // Treat pi0 slightly differently for now + Reaction* r1 = new Reaction(Ejectile); + r1->process_reaction(); + delete r1; + } } /*--------------------------------------------------*/ /*--------------------------------------------------*/ void SetEICSeed (int seed) { - fSeed = seed; + fSeed = seed; } ///*--------------------------------------------------*/ @@ -361,27 +363,27 @@ void SetEICSeed (int seed) { TString ExtractParticle(TString particle) { - /// Make the input particle case insansitive - particle.ToLower(); - if (particle.Contains("on")) { - particle.ReplaceAll("on", ""); - }; + /// Make the input particle case insansitive + particle.ToLower(); + if (particle.Contains("on")) { + particle.ReplaceAll("on", ""); + }; - if (particle.Contains("plus")) { - particle.ReplaceAll("plus", "+"); - } + if (particle.Contains("plus")) { + particle.ReplaceAll("plus", "+"); + } - if (particle.Contains("minus")) { - particle.ReplaceAll("minus", "-"); - } + if (particle.Contains("minus")) { + particle.ReplaceAll("minus", "-"); + } - if (particle.Contains("zero")) { - particle.ReplaceAll("zero", "0"); - } + if (particle.Contains("zero")) { + particle.ReplaceAll("zero", "0"); + } - particle[0] = toupper(particle[0]); - cout << "Particle: " << particle << endl; - return particle; + particle[0] = toupper(particle[0]); + cout << "Particle: " << particle << endl; + return particle; } diff --git a/src/eic_evgen/eic_pim.cc b/src/eic_evgen/eic_pim.cc index 671cda5..93346ff 100644 --- a/src/eic_evgen/eic_pim.cc +++ b/src/eic_evgen/eic_pim.cc @@ -42,7 +42,8 @@ int fWLessShell, fWLess1P9, fSDiff; unsigned long long int fNEvents, fNRecorded, fNGenerated, fWSqNeg, fNSigmaNeg, fNaN, fConserve, fNWeightUnphys, fNWeightReject, fLundRecorded, fNFile, fSolveEvents_0Sol, fSolveEvents_1Sol, fSolveEvents_2Sol; // SJDK 03/04/23 - Added in Qsq Min/Max, W Min/Max and t max (06/09/23) -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, fT_Max; +// 13/09/23 - SJDK - New generic HBeam value (rather than proton beam) +double fK, fm, fElectron_Kin_Col_GeV, fElectron_Kin_Col, fRand, fLumi, fuBcm2, fPI, fDEG2RAD, fRAD2DEG, fEBeam, fPBeam, fHBeam, 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, fT_Max; double fOmega_Theta_I, fOmega_Theta_F, fOmega_Theta_Col, fOmega_Phi_Col; diff --git a/src/eic_evgen/eic_pim.h b/src/eic_evgen/eic_pim.h index f022c9d..06d43e7 100644 --- a/src/eic_evgen/eic_pim.h +++ b/src/eic_evgen/eic_pim.h @@ -134,6 +134,8 @@ extern double fDEG2RAD; extern double fRAD2DEG; extern double fEBeam; extern double fPBeam; +// 13/09/23 - SJDK - New generic HBeam value (rather than proton beam) +extern double fHBeam; extern double fScatElec_Theta_I; extern double fScatElec_Theta_F; extern double fPion_Theta_I; // SJDK 19/12/22 - These should be removed in future, specific to pion reaction cases. Should be generic MesonX diff --git a/src/eic_evgen/process_routine/DEMP_Reaction.cc b/src/eic_evgen/process_routine/DEMP_Reaction.cc index 088b8d7..041bf58 100644 --- a/src/eic_evgen/process_routine/DEMP_Reaction.cc +++ b/src/eic_evgen/process_routine/DEMP_Reaction.cc @@ -178,7 +178,7 @@ void DEMP_Reaction::Init() { f_Ejectile_Theta_F = fEjectileX_Theta_F; cout << "Produced particle in exclusive production: " << rEjectile << "; with mass: " << f_Ejectile_Mass << " MeV "<< endl; - cout << fEBeam << " GeV electrons on " << fPBeam << " GeV ions" << endl; + cout << fEBeam << " GeV electrons on " << fHBeam << " GeV ions" << endl; if(UseSolve == true){ cout << rEjectile << " and " << rEjectile_scat_hadron << " 4-vectors calculated using Solve function" << endl; } @@ -188,22 +188,22 @@ void DEMP_Reaction::Init() { // 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) // See slide 11 in https://indico.cern.ch/event/1072579/contributions/4796856/attachments/2456676/4210776/CAP-EIC-June-7-2022-Seryi-r2.pdf // If available in the future, this could be replaced by some fixed function - if ((fEBeam == 5.0 ) && (fPBeam == 41.0) ){ + if ((fEBeam == 5.0 ) && (fHBeam == 41.0) ){ fLumi = 0.44e33; } - else if ((fEBeam == 5.0 ) && (fPBeam == 100.0) ){ + else if ((fEBeam == 5.0 ) && (fHBeam == 100.0) ){ fLumi = 3.68e33; } - else if ((fEBeam == 10.0 ) && (fPBeam == 100.0) ){ + else if ((fEBeam == 10.0 ) && (fHBeam == 100.0) ){ fLumi = 4.48e33; } - else if ((fEBeam == 18.0 ) && (fPBeam == 275.0) ){ + else if ((fEBeam == 18.0 ) && (fHBeam == 275.0) ){ fLumi = 1.54e33; } - else if ((fEBeam == 3.5 ) && (fPBeam == 20) ){ // EICC optimal beam energy combination + else if ((fEBeam == 3.5 ) && (fHBeam == 20) ){ // EICC optimal beam energy combination fLumi = 2e33; } - else if ((fEBeam == 2.8 ) && (fPBeam == 13) ){ // EICC lowest beam energy combination + else if ((fEBeam == 2.8 ) && (fHBeam == 13) ){ // EICC lowest beam energy combination fLumi = 0.7e33; } else{ @@ -666,7 +666,7 @@ TLorentzVector DEMP_Reaction::GetProtonVector_lab() { // fProton_Phi_Col = fPi; fProton_Phi_Col = fProton_incidence_phi; - fProton_Mom_Col = fPBeam * 1e3; + fProton_Mom_Col = fHBeam * 1e3; fVertex_X = 0.; fVertex_Y = 0.; fVertex_Z = 0.; From cef672554c6d3d9a99cdf20bb4e95c9bbe533322 Mon Sep 17 00:00:00 2001 From: Stephen JD Kay Date: Thu, 14 Sep 2023 09:23:33 -0600 Subject: [PATCH 21/21] Update DEMP_Reaction.cc Added % calculations for each cut. These are shown in a new column in the details (.txt) output. --- .../process_routine/DEMP_Reaction.cc | 46 +++++++++---------- 1 file changed, 23 insertions(+), 23 deletions(-) diff --git a/src/eic_evgen/process_routine/DEMP_Reaction.cc b/src/eic_evgen/process_routine/DEMP_Reaction.cc index 041bf58..9446a10 100644 --- a/src/eic_evgen/process_routine/DEMP_Reaction.cc +++ b/src/eic_evgen/process_routine/DEMP_Reaction.cc @@ -791,8 +791,8 @@ void DEMP_Reaction::Detail_Output() { DEMPDetails << endl; DEMPDetails << "Total events tried " << setw(20) << fNGenerated << endl; if(UseSolve == true){ - DEMPDetails << "Total events cut " << setw(20) << (qsq_ev + w_ev + w_neg_ev + fNaN + fConserve + t_ev + fNSigmaNeg + fSolveEvents_0Sol) << endl; - DEMPDetails << "Total events recorded " << setw(20) << fNRecorded << endl; + DEMPDetails << "Total events cut " << setw(20) << (qsq_ev + w_ev + w_neg_ev + fNaN + fConserve + t_ev + fNSigmaNeg + fSolveEvents_0Sol) << setw(20) << ((double) (qsq_ev + w_ev + w_neg_ev + fNaN + fConserve + t_ev + fNSigmaNeg + fSolveEvents_0Sol)/(double)fNGenerated)*100 << " %" << endl; + DEMPDetails << "Total events recorded " << setw(20) << fNRecorded << setw(20) << ((double)fNRecorded/(double)fNGenerated)*100 << " %" << endl; if (fNGenerated != (qsq_ev + w_ev + w_neg_ev + fNaN + fConserve + t_ev + fNSigmaNeg + fNRecorded + fSolveEvents_0Sol)){ DEMPDetails << "Total events cut + recorded = events tried? " << setw(20) << "NO! ERROR!" << endl; } @@ -802,8 +802,8 @@ void DEMP_Reaction::Detail_Output() { } else{ - DEMPDetails << "Total events cut " << setw(20) << (qsq_ev + w_ev + w_neg_ev + fNaN + fConserve + t_ev + fNSigmaNeg) << endl; - DEMPDetails << "Total events recorded " << setw(20) << fNRecorded << endl; + DEMPDetails << "Total events cut " << setw(20) << (qsq_ev + w_ev + w_neg_ev + fNaN + fConserve + t_ev + fNSigmaNeg) << setw(20) << ((double) (qsq_ev + w_ev + w_neg_ev + fNaN + fConserve + t_ev + fNSigmaNeg)/(double)fNGenerated)*100 << " %" << endl; + DEMPDetails << "Total events recorded " << setw(20) << fNRecorded << setw(20) << ((double)fNRecorded/(double)fNGenerated)*100 << " %" << endl; if (fNGenerated != (qsq_ev + w_ev + w_neg_ev + fNaN + fConserve + t_ev + fNSigmaNeg + fNRecorded)){ DEMPDetails << "Total events cut + recorded = events tried? " << setw(20) << "NO! ERROR!" << endl; } @@ -813,29 +813,29 @@ void DEMP_Reaction::Detail_Output() { } DEMPDetails << endl << "Cut details -" << endl; - DEMPDetails << "Events cut due to qsq < " << fQsq_Min << " or qsq > "<< fQsq_Max << " " << setw(20) << qsq_ev << endl; - DEMPDetails << "Events cut due to negative Wsq value " << setw(20) << w_neg_ev << endl; - DEMPDetails << "Events cut due to W < " << fW_Min << " or W > " << fW_Max << " " << setw(20) << w_ev << endl; + DEMPDetails << "Events cut due to qsq < " << fQsq_Min << " or qsq > "<< fQsq_Max << " " << setw(20) << qsq_ev << setw(20) << ((double)qsq_ev/(double)fNGenerated)*100 << " %" << endl; + DEMPDetails << "Events cut due to negative Wsq value " << setw(20) << w_neg_ev << setw(20) << ((double)w_neg_ev/(double)fNGenerated)*100 << " %" << endl; + DEMPDetails << "Events cut due to W < " << fW_Min << " or W > " << fW_Max << " " << setw(20) << w_ev << setw(20) << ((double)w_ev/(double)fNGenerated)*100 << " %" << endl; if(UseSolve == true){ - DEMPDetails << "Events cut due to solve function finding 0 solutions " << setw(20) << fSolveEvents_0Sol << endl; + DEMPDetails << "Events cut due to solve function finding 0 solutions " << setw(20) << fSolveEvents_0Sol << setw(20) << ((double)fSolveEvents_0Sol/(double)fNGenerated)*100 << " %" << endl; } - DEMPDetails << "Events cut due to ejectile (X) energy NaN " << setw(20) << fNaN << endl; - DEMPDetails << "Events cut due to conservation law check failure " << setw(20) << fConserve << endl; - DEMPDetails << "Events cut due to -t > " << fT_Max << "GeV " << setw(30) << t_ev << endl; - DEMPDetails << "Events cut due to -ve cross section value " << setw(20) << fNSigmaNeg << endl; + DEMPDetails << "Events cut due to ejectile (X) energy NaN " << setw(20) << fNaN << setw(20) << ((double)fNaN/(double)fNGenerated)*100 << " %" << endl; + DEMPDetails << "Events cut due to conservation law check failure " << setw(20) << fConserve << setw(20) << ((double)fConserve/(double)fNGenerated)*100 << " %" << endl; + DEMPDetails << "Events cut due to -t > " << fT_Max << "GeV " << setw(30) << t_ev << setw(20) << ((double)t_ev/(double)fNGenerated)*100 << " %" << endl; + DEMPDetails << "Events cut due to -ve cross section value " << setw(20) << fNSigmaNeg << setw(20) << ((double)fNSigmaNeg/(double)fNGenerated)*100 << " %" << endl; DEMPDetails << endl << "Conservation law checks details -" << endl; - DEMPDetails << "Total events PASSING conservation law check with tolerance " << fDiff << setw(17) << conserve << endl; - DEMPDetails << "Events cut due to energy conservation check ONLY " << setw(20) << ene << endl; - DEMPDetails << "Events cut due to momentum conservation check ONLY " << setw(20) << mom << endl; - DEMPDetails << "Events cut due to energy AND momentum conservation checks " << setw(20) << ene_mom << endl; - DEMPDetails << "Events cut due to px conservation law check " << setw(20) << mom_px << endl; - DEMPDetails << "Events cut due to py conservation law check " << setw(20) << mom_py << endl; - DEMPDetails << "Events cut due to pz conservation law check " << setw(20) << mom_pz << endl; - DEMPDetails << "Events cut due to px and py conservation law checks " << setw(20) << mom_pxpy << endl; - DEMPDetails << "Events cut due to px and pz conservation law checks " << setw(20) << mom_pxpz << endl; - DEMPDetails << "Events cut due to py and pz conservation law checks " << setw(20) << mom_pypz << endl; - DEMPDetails << "Events cut due to px, py and pz conservation law checks " << setw(20) << mom_pxpypz << endl; + DEMPDetails << "Total events PASSING conservation law check with tolerance " << fDiff << setw(17) << conserve << endl; + DEMPDetails << "Events cut due to energy conservation check ONLY " << setw(20) << ene << setw(20) << ((double)ene/(double)fNGenerated)*100 << " %" << endl; + DEMPDetails << "Events cut due to momentum conservation check ONLY " << setw(20) << mom << setw(20) << ((double)mom/(double)fNGenerated)*100 << " %" << endl; + DEMPDetails << "Events cut due to energy AND momentum conservation checks " << setw(20) << ene_mom << setw(20) << ((double)ene_mom/(double)fNGenerated)*100 << " %" << endl; + DEMPDetails << "Events cut due to px conservation law check " << setw(20) << mom_px << setw(20) << ((double)mom_px/(double)fNGenerated)*100 << " %" << endl; + DEMPDetails << "Events cut due to py conservation law check " << setw(20) << mom_py << setw(20) << ((double)mom_py/(double)fNGenerated)*100 << " %" << endl; + DEMPDetails << "Events cut due to pz conservation law check " << setw(20) << mom_pz << setw(20) << ((double)mom_pz/(double)fNGenerated)*100 << " %" << endl; + DEMPDetails << "Events cut due to px and py conservation law checks " << setw(20) << mom_pxpy << setw(20) << ((double)mom_pxpy/(double)fNGenerated)*100 << " %" << endl; + DEMPDetails << "Events cut due to px and pz conservation law checks " << setw(20) << mom_pxpz << setw(20) << ((double)mom_pxpz/(double)fNGenerated)*100 << " %" << endl; + DEMPDetails << "Events cut due to py and pz conservation law checks " << setw(20) << mom_pypz << setw(20) << ((double)mom_pypz/(double)fNGenerated)*100 << " %" << endl; + DEMPDetails << "Events cut due to px, py and pz conservation law checks " << setw(20) << mom_pxpypz << setw(20) << ((double)mom_pxpypz/(double)fNGenerated)*100 << " %" << endl; if(UseSolve == true){ DEMPDetails << endl << "Solve function, addtional info -" << endl;