From 05867ecae93617154bb1ff4e9eb7b3d5cb42d282 Mon Sep 17 00:00:00 2001
From: Zach Marshall <ZLMarshall@lbl.gov>
Date: Sun, 5 May 2019 19:17:09 +0200
Subject: [PATCH] RHadrons updates for stops and sbottoms in master
This merge request fixes two sets of problems.
First, there were problems in the logic that assigned hadronic
interactions to R-hadrons. The problem was essentially that one quark in
the light quark system was ignored when calculating the R-hadron
properties. This is (surprisingly) a quite good approximation as long as
you aren't making extremely heavy R-hadrons, and as long as you don't
have stop or sbottom R-mesons, where there is only one light quark in
the R-hadron. That's fixed, the list of R-hadron interactions has
changed somewhat, and the results change a little bit.
Second, we needed a few minor adjustments to get stop and sbottom
R-hadron simulation working. That should all be done. Validation plots
have been made and everything.
This MR gets master back in sync with 21.0 following !23144
---
.../RHadrons/python/RHadronMasses.py | 106 ++++++++++--------
.../RHadrons/src/G4ProcessHelper.cxx | 54 ++++-----
.../RHadrons/src/G4ProcessHelper.hh | 5 +-
.../preInclude.RHadronsPythia8.py | 6 +-
4 files changed, 91 insertions(+), 80 deletions(-)
diff --git a/Simulation/G4Extensions/RHadrons/python/RHadronMasses.py b/Simulation/G4Extensions/RHadrons/python/RHadronMasses.py
index 3d5370f5bf4..add95cfd479 100644
--- a/Simulation/G4Extensions/RHadrons/python/RHadronMasses.py
+++ b/Simulation/G4Extensions/RHadrons/python/RHadronMasses.py
@@ -26,6 +26,8 @@ The list of possible R-hadrons comes from the Pythia8 code, in src/RHadrons.cc (
first_mass_set = 4
offset_options = {
# Fundamental SUSY particles
+ 1000005 : [ 0 , False , '~b ' , -1./3. , 0.000 , 0.000 , 0.000 , 0.000 , 0.000 , 0.000 , 0.000 , 0.000 , 0.000 ] ,
+ 1000006 : [ 0 , False , '~t ' , 2./3. , 0.000 , 0.000 , 0.000 , 0.000 , 0.000 , 0.000 , 0.000 , 0.000 , 0.000 ] ,
1000021 : [ 0 , False , '~g ' , 0 , 0.000 , 0.000 , 0.000 , 0.000 , 0.000 , 0.000 , 0.000 , 0.000 , 0.000 ] ,
1000022 : [ 0 , False , '~chi10 ' , 0 , 0.000 , 0.000 , 0.000 , 0.000 , 0.000 , 0.000 , 0.000 , 0.000 , 0.000 ] ,
1000039 : [ 0 , False , '~Gr ' , 0 , 0.000 , 0.000 , 0.000 , 0.000 , 0.000 , 0.000 , 0.000 , 0.000 , 0.000 ] ,
@@ -82,47 +84,49 @@ offset_options = {
1095334 : [ 1000021 , True , '~g_Omegab*- ' , -1 , 5.900 , 5.662 , 5.662 , 6.580 , 0.000 , 5.790 , 5.662 , 6.580 , 0.000 ] ,
# Sbottom R-mesons
- 1000512 : [ 1000005 , True , '~B0 ' , 0 , 0.325 , 0.314 , 0.220 , 0.365 , 0.000 , 0.365 , 0.220 , 0.365 , 0.000 ] ,
- 1000522 : [ 1000005 , True , '~B- ' , -1 , 0.325 , 0.314 , 0.220 , 0.365 , 0.000 , 0.365 , 0.220 , 0.365 , 0.000 ] ,
- 1000532 : [ 1000005 , True , '~Bs0 ' , 0 , 0.500 , 0.466 , 0.419 , 0.540 , 0.000 , 0.530 , 0.419 , 0.540 , 0.000 ] ,
- 1000542 : [ 1000005 , True , '~Bc- ' , -1 , 1.500 , 1.630 , 1.550 , 1.710 , 0.000 , 1.700 , 1.550 , 1.710 , 0.000 ] ,
- 1000552 : [ 1000005 , True , '~etab0 ' , 0 , 4.800 , 4.730 , 4.730 , 5.500 , 0.000 , 4.730 , 4.730 , 5.500 , 0.000 ] ,
+ 1000512 : [ 1000005 , True , '~B0 ' , 0 , 0.325 , 0.314 , 0.220 , 0.365 , 0.000 , 0.365 , 0.220 , 0.365 , 0.000 ] ,
+ 1000522 : [ 1000005 , True , '~B- ' , -1 , 0.325 , 0.314 , 0.220 , 0.365 , 0.000 , 0.365 , 0.220 , 0.365 , 0.000 ] ,
+ 1000532 : [ 1000005 , True , '~Bs0 ' , 0 , 0.500 , 0.466 , 0.419 , 0.540 , 0.000 , 0.530 , 0.419 , 0.540 , 0.000 ] ,
+ 1000542 : [ 1000005 , True , '~Bc- ' , -1 , 1.500 , 1.630 , 1.550 , 1.710 , 0.000 , 1.700 , 1.550 , 1.710 , 0.000 ] ,
+ 1000552 : [ 1000005 , True , '~etab0 ' , 0 , 4.800 , 4.730 , 4.730 , 5.500 , 0.000 , 4.730 , 4.730 , 5.500 , 0.000 ] ,
# Sbottom R-baryons
- 1005113 : [ 1000005 , True , '~Sigmab- ' , -1 , 0.650 , 0.672 , 0.530 , 0.763 , 0.000 , 0.763 , 0.530 , 0.763 , 0.000 ] ,
- 1005211 : [ 1000005 , True , '~Sigmab0 ' , 0 , 0.650 , 0.496 , 0.171 , 0.632 , 0.000 , 0.632 , 0.171 , 0.632 , 0.000 ] ,
- 1005213 : [ 1000005 , True , '~Sigmab*0 ' , 0 , 0.650 , 0.672 , 0.530 , 0.763 , 0.000 , 0.763 , 0.530 , 0.763 , 0.000 ] ,
- 1005223 : [ 1000005 , True , '~Sigmab+ ' , 1 , 0.650 , 0.672 , 0.530 , 0.763 , 0.000 , 0.763 , 0.530 , 0.763 , 0.000 ] ,
- 1005311 : [ 1000005 , True , '~Xib- ' , -1 , 0.825 , 0.691 , 0.498 , 0.833 , 0.000 , 0.828 , 0.498 , 0.833 , 0.000 ] ,
- 1005313 : [ 1000005 , True , '~Xib*- ' , -1 , 0.825 , 0.810 , 0.686 , 0.922 , 0.000 , 0.917 , 0.686 , 0.922 , 0.000 ] ,
- 1005321 : [ 1000005 , True , '~Xib0 ' , 0 , 0.825 , 0.691 , 0.498 , 0.833 , 0.000 , 0.828 , 0.498 , 0.833 , 0.000 ] ,
- 1005323 : [ 1000005 , True , '~Xib*0 ' , 0 , 0.825 , 0.810 , 0.686 , 0.922 , 0.000 , 0.917 , 0.686 , 0.922 , 0.000 ] ,
- 1005333 : [ 1000005 , True , '~Omegab- ' , -1 , 1.000 , 0.952 , 0.863 , 1.095 , 0.000 , 1.075 , 0.863 , 1.095 , 0.000 ] ,
+ 1005113 : [ 1000005 , True , '~Sigmab- ' , -1 , 0.650 , 0.672 , 0.530 , 0.763 , 0.000 , 0.763 , 0.530 , 0.763 , 0.000 ] ,
+ 1005211 : [ 1000005 , True , '~Sigmab0 ' , 0 , 0.650 , 0.496 , 0.171 , 0.632 , 0.000 , 0.632 , 0.171 , 0.632 , 0.000 ] ,
+ 1005213 : [ 1000005 , True , '~Sigmab*0 ' , 0 , 0.650 , 0.672 , 0.530 , 0.763 , 0.000 , 0.763 , 0.530 , 0.763 , 0.000 ] ,
+ 1005223 : [ 1000005 , True , '~Sigmab+ ' , 1 , 0.650 , 0.672 , 0.530 , 0.763 , 0.000 , 0.763 , 0.530 , 0.763 , 0.000 ] ,
+ 1005311 : [ 1000005 , True , '~Xib- ' , -1 , 0.825 , 0.691 , 0.498 , 0.833 , 0.000 , 0.828 , 0.498 , 0.833 , 0.000 ] ,
+ 1005313 : [ 1000005 , True , '~Xib*- ' , -1 , 0.825 , 0.810 , 0.686 , 0.922 , 0.000 , 0.917 , 0.686 , 0.922 , 0.000 ] ,
+ 1005321 : [ 1000005 , True , '~Xib0 ' , 0 , 0.825 , 0.691 , 0.498 , 0.833 , 0.000 , 0.828 , 0.498 , 0.833 , 0.000 ] ,
+ 1005323 : [ 1000005 , True , '~Xib*0 ' , 0 , 0.825 , 0.810 , 0.686 , 0.922 , 0.000 , 0.917 , 0.686 , 0.922 , 0.000 ] ,
+ 1005333 : [ 1000005 , True , '~Omegab- ' , -1 , 1.000 , 0.952 , 0.863 , 1.095 , 0.000 , 1.075 , 0.863 , 1.095 , 0.000 ] ,
# Stop R-mesons
- 1000612 : [ 1000006 , True , '~T+ ' , 1 , 0.325 , 0.314 , 0.220 , 0.365 , 0.000 , 0.365 , 0.220 , 0.365 , 0.000 ] ,
- 1000622 : [ 1000006 , True , '~T0 ' , 0 , 0.325 , 0.314 , 0.220 , 0.365 , 0.000 , 0.365 , 0.220 , 0.365 , 0.000 ] ,
- 1000632 : [ 1000006 , True , '~Ts+ ' , 1 , 0.500 , 0.466 , 0.419 , 0.540 , 0.000 , 0.530 , 0.419 , 0.540 , 0.000 ] ,
- 1000642 : [ 1000006 , True , '~Tc0 ' , 0 , 1.500 , 1.630 , 1.550 , 1.710 , 0.000 , 1.700 , 1.550 , 1.710 , 0.000 ] ,
- 1000652 : [ 1000006 , True , '~etat+ ' , 1 , 4.800 , 4.730 , 4.730 , 5.500 , 0.000 , 4.730 , 4.730 , 5.500 , 0.000 ] ,
+ 1000612 : [ 1000006 , True , '~T+ ' , 1 , 0.325 , 0.314 , 0.220 , 0.365 , 0.000 , 0.365 , 0.220 , 0.365 , 0.000 ] ,
+ 1000622 : [ 1000006 , True , '~T0 ' , 0 , 0.325 , 0.314 , 0.220 , 0.365 , 0.000 , 0.365 , 0.220 , 0.365 , 0.000 ] ,
+ 1000632 : [ 1000006 , True , '~Ts+ ' , 1 , 0.500 , 0.466 , 0.419 , 0.540 , 0.000 , 0.530 , 0.419 , 0.540 , 0.000 ] ,
+ 1000642 : [ 1000006 , True , '~Tc0 ' , 0 , 1.500 , 1.630 , 1.550 , 1.710 , 0.000 , 1.700 , 1.550 , 1.710 , 0.000 ] ,
+ 1000652 : [ 1000006 , True , '~etat+ ' , 1 , 4.800 , 4.730 , 4.730 , 5.500 , 0.000 , 4.730 , 4.730 , 5.500 , 0.000 ] ,
# Stop R-baryons
- 1006113 : [ 1000006 , True , '~Sigmat0 ' , 0 , 0.650 , 0.672 , 0.530 , 0.763 , 0.000 , 0.763 , 0.530 , 0.763 , 0.000 ] ,
- 1006211 : [ 1000006 , True , '~Sigmat+ ' , 1 , 0.650 , 0.496 , 0.171 , 0.632 , 0.000 , 0.632 , 0.171 , 0.632 , 0.000 ] ,
- 1006213 : [ 1000006 , True , '~Sigmat*+ ' , 1 , 0.650 , 0.672 , 0.530 , 0.763 , 0.000 , 0.763 , 0.530 , 0.763 , 0.000 ] ,
- 1006223 : [ 1000006 , True , '~Sigmat++ ' , 2 , 0.650 , 0.672 , 0.530 , 0.763 , 0.000 , 0.763 , 0.530 , 0.763 , 0.000 ] ,
- 1006311 : [ 1000006 , True , '~Xit0 ' , 0 , 0.825 , 0.691 , 0.498 , 0.833 , 0.000 , 0.828 , 0.498 , 0.833 , 0.000 ] ,
- 1006313 : [ 1000006 , True , '~Xit*0 ' , 0 , 0.825 , 0.810 , 0.686 , 0.922 , 0.000 , 0.917 , 0.686 , 0.922 , 0.000 ] ,
- 1006321 : [ 1000006 , True , '~Xit+ ' , 1 , 0.825 , 0.691 , 0.498 , 0.833 , 0.000 , 0.828 , 0.498 , 0.833 , 0.000 ] ,
- 1006323 : [ 1000006 , True , '~Xit*+ ' , 1 , 0.825 , 0.810 , 0.686 , 0.922 , 0.000 , 0.917 , 0.686 , 0.922 , 0.000 ] ,
- 1006333 : [ 1000006 , True , '~Omegat0 ' , 0 , 1.000 , 0.952 , 0.863 , 1.095 , 0.000 , 1.075 , 0.863 , 1.095 , 0.000 ] ,
+ 1006113 : [ 1000006 , True , '~Sigmat0 ' , 0 , 0.650 , 0.672 , 0.530 , 0.763 , 0.000 , 0.763 , 0.530 , 0.763 , 0.000 ] ,
+ 1006211 : [ 1000006 , True , '~Sigmat+ ' , 1 , 0.650 , 0.496 , 0.171 , 0.632 , 0.000 , 0.632 , 0.171 , 0.632 , 0.000 ] ,
+ 1006213 : [ 1000006 , True , '~Sigmat*+ ' , 1 , 0.650 , 0.672 , 0.530 , 0.763 , 0.000 , 0.763 , 0.530 , 0.763 , 0.000 ] ,
+ 1006223 : [ 1000006 , True , '~Sigmat++ ' , 2 , 0.650 , 0.672 , 0.530 , 0.763 , 0.000 , 0.763 , 0.530 , 0.763 , 0.000 ] ,
+ 1006311 : [ 1000006 , True , '~Xit0 ' , 0 , 0.825 , 0.691 , 0.498 , 0.833 , 0.000 , 0.828 , 0.498 , 0.833 , 0.000 ] ,
+ 1006313 : [ 1000006 , True , '~Xit*0 ' , 0 , 0.825 , 0.810 , 0.686 , 0.922 , 0.000 , 0.917 , 0.686 , 0.922 , 0.000 ] ,
+ 1006321 : [ 1000006 , True , '~Xit+ ' , 1 , 0.825 , 0.691 , 0.498 , 0.833 , 0.000 , 0.828 , 0.498 , 0.833 , 0.000 ] ,
+ 1006323 : [ 1000006 , True , '~Xit*+ ' , 1 , 0.825 , 0.810 , 0.686 , 0.922 , 0.000 , 0.917 , 0.686 , 0.922 , 0.000 ] ,
+ 1006333 : [ 1000006 , True , '~Omegat0 ' , 0 , 1.000 , 0.952 , 0.863 , 1.095 , 0.000 , 1.075 , 0.863 , 1.095 , 0.000 ] ,
}
# Now programmatically calculate the missing spectra
+# These are designed to just flop the rho above or below the gluinoball
+gb_offset = offset_options[1009113][first_mass_set+5]-offset_options[1009113][first_mass_set+1]
for pid in offset_options:
# Skip fundamental SUSY particles and R-glueball
if offset_options[pid][0] == 0 or pid == 1000993: continue
# Setup #4 to be mass set #1 but with 1009113 matching mass set 5
- offset_options[pid][first_mass_set+4] = offset_options[pid][first_mass_set+1]+0.088
- # Setup #8 to be mass set #5 but with 1009113 matches mass set 2
- offset_options[pid][first_mass_set+8] = offset_options[pid][first_mass_set+5]-0.088
+ offset_options[pid][first_mass_set+4] = offset_options[pid][first_mass_set+1]+gb_offset
+ # Setup #8 to be mass set #5 but with 1009113 matching mass set 1
+ offset_options[pid][first_mass_set+8] = offset_options[pid][first_mass_set+5]-gb_offset
@@ -139,13 +143,16 @@ def charge( c ):
raise RuntimeError('Unexpected charge: '+str(n))
-def get_quarks( x ):
+def get_quarks( y ):
""" Function to return a list of quarks in a hadron
"""
+ x = abs(y)
+ # For stop/sbottom mesons, just the last quark!
+ if '000' in str(x): return str(x)[5:6]
# For mesons, just two quarks
- if '00' in str(x): return str(x)[4:5]
+ if '00' in str(x): return str(x)[4:6]
# For baryons, three quarks
- return str(x)[3:5]
+ return str(x)[3:6]
def is_baryon( x ):
@@ -254,8 +261,8 @@ def update_PDG_table(input_file, pdg_table, mass_spectrum=1):
# Note that we follow the Pythia6 convention of *including* fundamental SUSY particles
# The format is VERY specific; needs mass and width (we always set the width to 0)
# Mass is in MeV here, rather than GeV as in the dictionary
- out_file.write('\nM %i %11.5E +0.0E+00 -0.0E+00 %s %s'%(pid,masses[pid]*1000.,offset_options[pid][2],charge(offset_options[pid][3])))
- out_file.write('\nW %i %11.5E +0.0E+00 -0.0E+00 %s %s'%(pid,0.E+00,offset_options[pid][2],charge(offset_options[pid][3])))
+ out_file.write('\nM %i %11.7E +0.0E+00 -0.0E+00 %s %s'%(pid,masses[pid]*1000.,offset_options[pid][2],charge(offset_options[pid][3])))
+ out_file.write('\nW %i %11.7E +0.0E+00 -0.0E+00 %s %s'%(pid,0.E+00,offset_options[pid][2],charge(offset_options[pid][3])))
# Done writing all the lines! Clean up if necessary
if type(pdg_table) is str: out_file.close()
@@ -377,24 +384,31 @@ def get_interaction_list(input_file, interaction_file='ProcessList.txt', mass_sp
# All of them are on the list of incoming RHadrons
# Deal with strangeness
# Approximation! Bottom number -> -Charm number -> Strangeness
- # Approximation needed because I don't know how outgoing SM charms are treated in G4 at the moment
+ # Approximation needed because outgoing SM charms are not treated in G4 at the moment
s_number = 0
- if '3' in get_quarks(pid) or '4' in get_quarks(pid) or '5' in get_quarks(pid):
- if len(get_quarks(pid))>2:
- s_number = -(get_quarks(pid).count('3')+get_quarks(pid).count('4')+get_quarks(pid).count('5')) if pid>0 else get_quarks(pid).count('3')+get_quarks(pid).count('4')+get_quarks(pid).count('5')
+ my_q = get_quarks(pid)
+ if '3' in my_q or '4' in my_q or '5' in my_q:
+ if len(my_q)>2:
+ # Gluino R-baryons
+ s_number = -(my_q.count('3')+my_q.count('4')+my_q.count('5')) if pid>0 else my_q.count('3')+my_q.count('4')+my_q.count('5')
+ elif len(my_q)>1:
+ # Squark R-baryons or Gluino R-mesons
+ if my_q.count('3') + my_q.count('4') + my_q.count('5')>1: s_number=0
+ elif offset_options[abs(pid)][3]==0 and ('3' in my_q or '5' in my_q): s_number=1 if pid>0 else -1
+ elif offset_options[abs(pid)][3]==0 and '4' in my_q: s_number=1 if pid<0 else -1
+ elif '3' in my_q or '5' in my_q: s_number=offset_options[abs(pid)][3]
+ elif '4' in my_q: s_number=-offset_options[abs(pid)][3]
else:
- if get_quarks(pid).count('3') + get_quarks(pid).count('4') + get_quarks(pid).count('5')>1: s_number=0
- elif offset_options[abs(pid)][3]==0 and ('3' in get_quarks(pid) or '5' in get_quarks(pid)): s_number=1 if pid>0 else -1
- elif offset_options[abs(pid)][3]==0 and '4' in get_quarks(pid): s_number=1 if pid<0 else -1
- elif '3' in get_quarks(pid) or '5' in get_quarks(pid): s_number=offset_options[abs(pid)][3]
- elif '4' in get_quarks(pid): s_number=-offset_options[abs(pid)][3]
+ # Squark R-mesons
+ s_number = my_q.count('3') + my_q.count('4') + my_q.count('5')
+ s_number = s_number if pid>0 else -s_number
else: s_number=0
# Build the dictionary
pid_name = offset_options[pid][2].strip() if pid>0 else anti_name(offset_options[abs(pid)][2]).strip()
incoming_rhadrons[pid_name] = [ offset_options[abs(pid)][3] , is_baryon(pid) , s_number ]
# Smaller list of outgoing rhadrons.
# No charm or bottom
- if '4' in get_quarks(pid) or '5' in get_quarks(pid): continue
+ if '4' in my_q or '5' in my_q: continue
outgoing_rhadrons[pid_name] = [ offset_options[abs(pid)][3] , is_baryon(pid) , s_number ]
# Add all our R-hadrons to the table
diff --git a/Simulation/G4Extensions/RHadrons/src/G4ProcessHelper.cxx b/Simulation/G4Extensions/RHadrons/src/G4ProcessHelper.cxx
index 7fccf9e08da..4b3e490764a 100644
--- a/Simulation/G4Extensions/RHadrons/src/G4ProcessHelper.cxx
+++ b/Simulation/G4Extensions/RHadrons/src/G4ProcessHelper.cxx
@@ -179,15 +179,11 @@ G4bool G4ProcessHelper::ApplicabilityTester(const G4ParticleDefinition& aPart){
G4double G4ProcessHelper::GetInclusiveCrossSection(const G4DynamicParticle *aParticle,
const G4Element *anElement){
-
//We really do need a dedicated class to handle the cross sections. They might not always be constant
-
//Disassemble the PDG-code
-
G4int thePDGCode = aParticle->GetDefinition()->GetPDGEncoding();
double boost = (aParticle->GetKineticEnergy()+aParticle->GetMass())/aParticle->GetMass();
- // G4cout<<"thePDGCode: "<<thePDGCode<<G4endl;
G4double theXsec = 0;
G4String name = aParticle->GetDefinition()->GetParticleName();
@@ -256,6 +252,11 @@ G4double G4ProcessHelper::GetInclusiveCrossSection(const G4DynamicParticle *aPar
}
ReactionProduct G4ProcessHelper::GetFinalState(const G4Track& aTrack, G4ParticleDefinition*& aTarget){
+ return GetFinalStateInternal(aTrack,aTarget,false);
+}
+
+// Version where we know if we baryonize already
+ReactionProduct G4ProcessHelper::GetFinalStateInternal(const G4Track& aTrack,G4ParticleDefinition*& aTarget, const bool baryonize_failed) {
const G4DynamicParticle* aDynamicParticle = aTrack.GetDynamicParticle();
@@ -304,16 +305,17 @@ ReactionProduct G4ProcessHelper::GetFinalState(const G4Track& aTrack, G4Particle
bool baryonise=false;
- if(reggemodel
- &&CLHEP::RandFlat::shoot()>0.9
- &&(
+ if(!baryonize_failed
+ && reggemodel
+ && CLHEP::RandFlat::shoot()>0.9
+ && (
(CustomPDGParser::s_isMesonino(theIncidentPDG)&&theIncidentPDG>0)
||
CustomPDGParser::s_isRMeson(theIncidentPDG)
)
- )
+ ){
baryonise=true;
-
+ }
//Making a pointer directly to the ReactionProductList we are looking at. Makes life easier :-)
ReactionProductList* aReactionProductList = &((*theReactionMap)[theIncidentPDG]);
@@ -362,14 +364,17 @@ ReactionProduct G4ProcessHelper::GetFinalState(const G4Track& aTrack, G4Particle
} else {
G4cerr << "ReactionProduct has unsupported number of secondaries: "<<secondaries<<G4endl;
}
- } /*else {
- G4cout<<"There was an impossible process"<<G4endl;
- }*/
+ }
}
- // G4cout<<"The size of the ReactionProductList is: "<<theReactionProductList.size()<<G4endl;
- if (theReactionProductList.size()==0) G4Exception("G4ProcessHelper", "NoProcessPossible", FatalException,
- "GetFinalState: No process could be selected from the given list.");
+ if (theReactionProductList.size()==0 && baryonize_failed){
+ G4Exception("G4ProcessHelper", "NoProcessPossible", FatalException,
+ "GetFinalState: No process could be selected from the given list.");
+ } else if (theReactionProductList.size()==0 && !baryonize_failed) {
+ // Baryonization had not yet failed -- try again
+ G4cout << "G4ProcessHelper::GetFinalStateInternal WARNING Could not select an appropriate process in first pass" << G4endl;
+ return GetFinalStateInternal(aTrack,aTarget,true);
+ }
// For the Regge model no phase space considerations. We pick a process at random
if(reggemodel)
@@ -403,19 +408,15 @@ ReactionProduct G4ProcessHelper::GetFinalState(const G4Track& aTrack, G4Particle
TwotoThreeFlag.push_back(true);
}
Probabilities.push_back(CumulatedProbability);
- // G4cout<<"Pushing back cumulated probability: "<<CumulatedProbability<<G4endl;
}
//Renormalising probabilities
- // G4cout<<"Probs: ";
for (std::vector<G4double>::iterator it = Probabilities.begin();
it != Probabilities.end();
it++)
{
*it /= CumulatedProbability;
- // G4cout<<*it<<" ";
}
- // G4cout<<G4endl;
// Choosing ReactionProduct
@@ -428,14 +429,10 @@ ReactionProduct G4ProcessHelper::GetFinalState(const G4Track& aTrack, G4Particle
while(!selected && tries < 100){
i=0;
G4double dice = CLHEP::RandFlat::shoot();
- // G4cout<<"What's the dice?"<<dice<<G4endl;
while(dice>Probabilities[i] && i<theReactionProductList.size()){
- // G4cout<<"i: "<<i<<G4endl;
i++;
}
- // G4cout<<"Chosen i: "<<i<<G4endl;
-
if(!TwotoThreeFlag[i]) {
// 2 -> 2 processes are chosen immediately
selected = true;
@@ -459,9 +456,6 @@ ReactionProduct G4ProcessHelper::GetFinalState(const G4Track& aTrack, G4Particle
}
if(tries>=100) G4cerr<<"Could not select process!!!!"<<G4endl;
- // G4cout<<"So far so good"<<G4endl;
- // G4cout<<"Sec's: "<<theReactionProductList[i].size()<<G4endl;
-
//Updating checkfraction:
if (theReactionProductList[i].size()==2) {
n_22++;
@@ -476,7 +470,7 @@ ReactionProduct G4ProcessHelper::GetFinalState(const G4Track& aTrack, G4Particle
return theReactionProductList[i];
}
-G4double G4ProcessHelper::ReactionProductMass(const ReactionProduct& aReaction,const G4DynamicParticle* aDynamicParticle) const {
+G4double G4ProcessHelper::ReactionProductMass(const ReactionProduct& aReaction,const G4DynamicParticle* aDynamicParticle) const{
// Incident energy:
G4double E_incident = aDynamicParticle->GetTotalEnergy();
//G4cout<<"Total energy: "<<E_incident<<" Kinetic: "<<aDynamicParticle->GetKineticEnergy()<<G4endl;
@@ -492,7 +486,7 @@ G4double G4ProcessHelper::ReactionProductMass(const ReactionProduct& aReaction,c
//G4cout<<"Mass contrib: "<<(particleTable->FindParticle(*r_it)->GetPDGMass())/CLHEP::MeV<<" MeV"<<G4endl;
M_after += particleTable->FindParticle(*r_it)->GetPDGMass();
}
- //G4cout<<"Intending to return this ReactionProductMass: "<<(sqrts - M_after)/CLHEP::MeV<<" MeV"<<G4endl;
+ //G4cout<<"Intending to return this ReactionProductMass: " << sqrts << " - " << M_after << " MeV"<<G4endl;
return sqrts - M_after;
}
@@ -507,7 +501,7 @@ G4bool G4ProcessHelper::ReactionGivesBaryon(const ReactionProduct& aReaction) co
return false;
}
-G4double G4ProcessHelper::PhaseSpace(const ReactionProduct& aReaction,const G4DynamicParticle* aDynamicParticle) const{
+G4double G4ProcessHelper::PhaseSpace(const ReactionProduct& aReaction,const G4DynamicParticle* aDynamicParticle) const {
G4double qValue = ReactionProductMass(aReaction,aDynamicParticle);
// Eq 4 of https://arxiv.org/pdf/hep-ex/0404001.pdf
G4double phi = sqrt(1+qValue/(2*0.139*CLHEP::GeV))*pow(qValue/(1.1*CLHEP::GeV),3./2.);
@@ -516,7 +510,7 @@ G4double G4ProcessHelper::PhaseSpace(const ReactionProduct& aReaction,const G4Dy
void G4ProcessHelper::ReadAndParse(const G4String& str,
std::vector<G4String>& tokens,
- const G4String& delimiters) const
+ const G4String& delimiters)
{
// Skip delimiters at beginning.
G4String::size_type lastPos = str.find_first_not_of(delimiters, 0);
diff --git a/Simulation/G4Extensions/RHadrons/src/G4ProcessHelper.hh b/Simulation/G4Extensions/RHadrons/src/G4ProcessHelper.hh
index f833bcf6cdf..ee30bd1aff7 100644
--- a/Simulation/G4Extensions/RHadrons/src/G4ProcessHelper.hh
+++ b/Simulation/G4Extensions/RHadrons/src/G4ProcessHelper.hh
@@ -54,6 +54,9 @@ private:
static G4ProcessHelper* pinstance;
+ // Version where we know if we baryonize already
+ ReactionProduct GetFinalStateInternal(const G4Track& aTrack,G4ParticleDefinition*& aTarget, const bool baryonize_failed);
+
G4double Regge(const double boost) const;
G4double Pom(const double boost) const;
@@ -65,7 +68,7 @@ private:
G4bool ReactionGivesBaryon(const ReactionProduct& aReaction) const;
void ReadAndParse(const G4String& str,
std::vector<G4String>& tokens,
- const G4String& delimiters = " ") const;
+ const G4String& delimiters = " ");
//Map of applicable particles
std::map<const G4ParticleDefinition*,G4bool> known_particles;
diff --git a/Simulation/SimulationJobOptions/share/specialConfig/preInclude.RHadronsPythia8.py b/Simulation/SimulationJobOptions/share/specialConfig/preInclude.RHadronsPythia8.py
index 54de5f8ec60..37bcd0a172d 100644
--- a/Simulation/SimulationJobOptions/share/specialConfig/preInclude.RHadronsPythia8.py
+++ b/Simulation/SimulationJobOptions/share/specialConfig/preInclude.RHadronsPythia8.py
@@ -66,9 +66,9 @@ try:
# MC Channel Number. Try the standard two spots, and fall back to the run number for evgen
if 'mc_channel_number' in f.infos and len(f.infos['mc_channel_number'])>0:
runNumber = f.infos['mc_channel_number'][0]
- if runNumber<=0 and 'mc_channel_number' in f.infos['tag_info']:
+ elif 'mc_channel_number' in f.infos['tag_info']:
runNumber = f.infos['tag_info']['mc_channel_number']
- if runNumber<=0:
+ else:
runNumber = f.infos['run_number'][0]
# This is also used for digitization, so protect in case we're there
if "StreamHITS" in f.infos["stream_names"]:
@@ -181,7 +181,7 @@ if lifetime>0.:
else:
addLineToPhysicsConfiguration("DoDecays","1")
addLineToPhysicsConfiguration("HadronLifeTime", str(lifetime))
- # If we reading particle records, and the lifetime is short, stop them as well
+ # If we are reading particle records, and the lifetime is short, stop them as well
if lifetime<1. and hasattr(runArgs,'inputEVNT_TRFile'):
addLineToPhysicsConfiguration("DoDecays","1")
addLineToPhysicsConfiguration("HadronLifeTime", 0.000001)
--
GitLab