ATLAS Offline Software
Loading...
Searching...
No Matches
RHadronMasses.py
Go to the documentation of this file.
1# Copyright (C) 2002-2023 CERN for the benefit of the ATLAS collaboration
2
3# This file contains a number of helper functions for defining R-hadron mass spectra
4# A large data table at the top is then used in several of the helper functions
5
6
7"""
8The offset options. Dictionary of PDG IDs with offsets.
9Anti-particles by definition have the same mass as standard particles.
10Columns are:
11 PDG ID
12 PDG ID of SUSY constituent
13 Has an anti-particle (i.e. need to define a particle with an opposite-signed PDG ID)
14 Name (for PDG table)
15 Charge (for PDG table)
16
17 Mass Offset in the Pythia6 simulation configuration
18 8 mass options documented in the INT note of the stopped particle analysis:
19 https://cds.cern.ch/record/2665511
20 First 3 are the meson fits; 4th is the first specrum shifted to make the gluino R-glueball heaver
21 than the lightest R-mesons, but with the same relative splitting otherwise.
22 Second set of four follow the same pattern, but for the baryon fits.
23
24The list of possible R-hadrons comes from the Pythia8 code, in src/RHadrons.cc (at the top)
25"""
26
27first_mass_set = 4
28offset_options = {
29# Fundamental SUSY particles
30 1000005 : [ 0 , False , '~b ' , -1./3. , 0.000 , 0.000 , 0.000 , 0.000 , 0.000 , 0.000 , 0.000 , 0.000 , 0.000 ] ,
31 1000006 : [ 0 , False , '~t ' , 2./3. , 0.000 , 0.000 , 0.000 , 0.000 , 0.000 , 0.000 , 0.000 , 0.000 , 0.000 ] ,
32 1000021 : [ 0 , False , '~g ' , 0 , 0.000 , 0.000 , 0.000 , 0.000 , 0.000 , 0.000 , 0.000 , 0.000 , 0.000 ] ,
33 1000022 : [ 0 , False , '~chi10 ' , 0 , 0.000 , 0.000 , 0.000 , 0.000 , 0.000 , 0.000 , 0.000 , 0.000 , 0.000 ] ,
34 1000039 : [ 0 , False , '~Gr ' , 0 , 0.000 , 0.000 , 0.000 , 0.000 , 0.000 , 0.000 , 0.000 , 0.000 , 0.000 ] ,
35# Gluino R-glueball
36 1000993 : [ 1000021 , False , '~g_glueball ' , 0 , 0.700 , 0.700 , 0.700 , 0.700 , 0.700 , 0.700 , 0.700 , 0.700 , 0.700 ] ,
37# Gluino q-qbar R-mesons
38 1009113 : [ 1000021 , False , '~g_rho ' , 0 , 0.650 , 0.610 , 0.403 , 0.717 , 0.000 , 0.717 , 0.403 , 0.717 , 0.000 ] ,
39 1009223 : [ 1000021 , False , '~g_omega ' , 0 , 0.650 , 0.610 , 0.403 , 0.717 , 0.000 , 0.717 , 0.403 , 0.717 , 0.000 ] ,
40 1009333 : [ 1000021 , False , '~g_phi ' , 0 , 1.800 , 0.924 , 0.828 , 1.074 , 0.000 , 1.054 , 0.828 , 1.074 , 0.000 ] ,
41 1009443 : [ 1000021 , False , '~g_J/Psi ' , 0 , 3.400 , 3.259 , 3.099 , 3.419 , 0.000 , 3.399 , 3.099 , 3.419 , 0.000 ] ,
42 1009553 : [ 1000021 , False , '~g_Upsilon ' , 0 , 9.460 , 9.460 , 9.460 ,11.000 , 0.000 , 9.460 , 9.460 ,11.000 ,10.000 ] ,
43# Light-flavor Gluino R-mesons
44 1009213 : [ 1000021 , True , '~g_rho+ ' , 1 , 0.650 , 0.610 , 0.403 , 0.717 , 0.000 , 0.717 , 0.403 , 0.717 , 0.000 ] ,
45# Strange Gluino R-mesons
46 1009313 : [ 1000021 , True , '~g_K*0 ' , 0 , 0.825 , 0.768 , 0.620 , 0.891 , 0.000 , 0.886 , 0.620 , 0.891 , 0.000 ] ,
47 1009323 : [ 1000021 , True , '~g_K*+ ' , 1 , 0.825 , 0.768 , 0.620 , 0.891 , 0.000 , 0.886 , 0.620 , 0.891 , 0.000 ] ,
48# Charm Gluino R-mesons
49 1009413 : [ 1000021 , True , '~g_D*+ ' , 1 , 2.000 , 1.940 , 1.843 , 2.067 , 0.000 , 2.062 , 1.843 , 2.067 , 0.000 ] ,
50 1009423 : [ 1000021 , True , '~g_D*0 ' , 0 , 2.000 , 1.940 , 1.843 , 2.067 , 0.000 , 2.062 , 1.843 , 2.067 , 0.000 ] ,
51 1009433 : [ 1000021 , True , '~g_D*s+ ' , 1 , 2.200 , 2.094 , 2.034 , 2.248 , 0.000 , 2.228 , 2.034 , 2.248 , 0.000 ] ,
52# Bottom Gluino R-mesons
53 1009513 : [ 1000021 , True , '~g_B*0 ' , 0 , 5.000 , 5.043 , 5.039 , 5.859 , 0.000 , 5.094 , 5.039 , 5.859 , 0.000 ] ,
54 1009523 : [ 1000021 , True , '~g_B*+ ' , 1 , 5.000 , 5.043 , 5.039 , 5.859 , 0.000 , 5.094 , 5.039 , 5.859 , 0.000 ] ,
55 1009533 : [ 1000021 , True , '~g_B*s0 ' , 0 , 5.200 , 5.195 , 5.195 , 6.039 , 0.000 , 5.259 , 5.195 , 6.039 , 0.000 ] ,
56 1009543 : [ 1000021 , True , '~g_B*c+ ' , 1 , 7.000 , 6.360 , 6.280 , 7.210 , 0.000 , 6.430 , 6.280 , 7.210 , 0.000 ] ,
57# Light-flavor singlet Gluino R-baryons
58 1093214 : [ 1000021 , True , '~g_Lambda0 ' , 0 , 1.150 , 0.562 , 0.280 , 0.715 , 0.000 , 0.715 , 0.280 , 0.715 , 0.000 ] ,
59# Charm singlet Gluino R-baryons
60 1094214 : [ 1000021 , True , '~g_Sigmac*+ ' , 1 , 2.300 , 1.726 , 1.489 , 1.885 , 0.000 , 1.885 , 1.489 , 1.885 , 0.000 ] ,
61 1094314 : [ 1000021 , True , '~g_Xic*0 ' , 0 , 2.300 , 1.878 , 1.688 , 2.050 , 0.000 , 2.050 , 1.688 , 2.050 , 0.000 ] ,
62 1094324 : [ 1000021 , True , '~g_Xic*+ ' , 1 , 2.300 , 1.878 , 1.688 , 2.050 , 0.000 , 2.050 , 1.688 , 2.050 , 0.000 ] ,
63# Bottom singlet Gluino R-baryons
64 1095214 : [ 1000021 , True , '~g_Sigmab*0 ' , 0 , 5.600 , 4.826 , 4.796 , 5.660 , 0.000 , 4.915 , 4.796 , 5.660 , 0.000 ] ,
65 1095314 : [ 1000021 , True , '~g_Xib*- ' , -1 , 5.750 , 4.978 , 4.970 , 5.840 , 0.000 , 5.080 , 4.970 , 5.840 , 0.000 ] ,
66 1095324 : [ 1000021 , True , '~g_Xib*0 ' , 0 , 5.750 , 4.978 , 4.970 , 5.840 , 0.000 , 5.080 , 4.970 , 5.840 , 0.000 ] ,
67# Light flavor Gluino R-baryons
68 1091114 : [ 1000021 , True , '~g_Delta- ' , -1 , 0.975 , 0.812 , 0.530 , 0.965 , 0.000 , 0.965 , 0.530 , 0.965 , 0.000 ] ,
69 1092114 : [ 1000021 , True , '~g_Delta0 ' , 0 , 0.975 , 0.812 , 0.530 , 0.965 , 0.000 , 0.965 , 0.530 , 0.965 , 0.000 ] ,
70 1092214 : [ 1000021 , True , '~g_Delta+ ' , 1 , 0.975 , 0.812 , 0.530 , 0.965 , 0.000 , 0.965 , 0.530 , 0.965 , 0.000 ] ,
71 1092224 : [ 1000021 , True , '~g_Delta++ ' , 2 , 0.975 , 0.812 , 0.530 , 0.965 , 0.000 , 0.965 , 0.530 , 0.965 , 0.000 ] ,
72# Strange Gluino R-baryons
73 1093114 : [ 1000021 , True , '~g_Sigma*- ' , -1 , 1.150 , 1.094 , 0.858 , 1.260 , 0.000 , 1.260 , 0.859 , 1.260 , 0.000 ] ,
74 1093224 : [ 1000021 , True , '~g_Sigma*+ ' , 1 , 1.150 , 1.094 , 0.858 , 1.260 , 0.000 , 1.260 , 0.859 , 1.260 , 0.000 ] ,
75 1093314 : [ 1000021 , True , '~g_Xi*- ' , -1 , 1.300 , 1.246 , 1.058 , 1.440 , 0.000 , 1.425 , 1.058 , 1.440 , 0.000 ] ,
76 1093324 : [ 1000021 , True , '~g_Xi*0 ' , 0 , 1.300 , 1.246 , 1.058 , 1.440 , 0.000 , 1.425 , 1.058 , 1.440 , 0.000 ] ,
77 1093334 : [ 1000021 , True , '~g_Omega- ' , -1 , 1.600 , 1.398 , 1.257 , 1.620 , 0.000 , 1.590 , 1.257 , 1.620 , 0.000 ] ,
78# Charm Gluino R-baryons
79 1094114 : [ 1000021 , True , '~g_Sigmac*0 ' , 0 , 2.300 , 2.258 , 2.068 , 2.430 , 0.000 , 2.430 , 2.068 , 2.430 , 0.000 ] ,
80 1094224 : [ 1000021 , True , '~g_Sigmac*++' , 2 , 2.300 , 2.258 , 2.068 , 2.430 , 0.000 , 2.430 , 2.068 , 2.430 , 0.000 ] ,
81 1094334 : [ 1000021 , True , '~g_Omegac*0 ' , 0 , 2.300 , 2.562 , 2.466 , 2.790 , 0.000 , 2.760 , 2.466 , 2.790 , 0.000 ] ,
82# Bottom Gluino R-baryons
83 1095114 : [ 1000021 , True , '~g_Sigmab*- ' , -1 , 5.600 , 5.358 , 5.350 , 6.220 , 0.000 , 5.460 , 5.350 , 6.220 , 0.000 ] ,
84 1095224 : [ 1000021 , True , '~g_Sigmab*+ ' , 1 , 5.600 , 5.358 , 5.350 , 6.220 , 0.000 , 5.460 , 5.350 , 6.220 , 0.000 ] ,
85 1095334 : [ 1000021 , True , '~g_Omegab*- ' , -1 , 5.900 , 5.662 , 5.662 , 6.580 , 0.000 , 5.790 , 5.662 , 6.580 , 0.000 ] ,
86
87# Sbottom R-mesons
88 1000512 : [ 1000005 , True , '~B0 ' , 0 , 0.325 , 0.314 , 0.220 , 0.365 , 0.000 , 0.365 , 0.220 , 0.365 , 0.000 ] ,
89 1000522 : [ 1000005 , True , '~B- ' , -1 , 0.325 , 0.314 , 0.220 , 0.365 , 0.000 , 0.365 , 0.220 , 0.365 , 0.000 ] ,
90 1000532 : [ 1000005 , True , '~Bs0 ' , 0 , 0.500 , 0.466 , 0.419 , 0.540 , 0.000 , 0.530 , 0.419 , 0.540 , 0.000 ] ,
91 1000542 : [ 1000005 , True , '~Bc- ' , -1 , 1.500 , 1.630 , 1.550 , 1.710 , 0.000 , 1.700 , 1.550 , 1.710 , 0.000 ] ,
92 1000552 : [ 1000005 , True , '~etab0 ' , 0 , 4.800 , 4.730 , 4.730 , 5.500 , 0.000 , 4.730 , 4.730 , 5.500 , 0.000 ] ,
93# Sbottom R-baryons
94 1005113 : [ 1000005 , True , '~Sigmab- ' , -1 , 0.650 , 0.672 , 0.530 , 0.763 , 0.000 , 0.763 , 0.530 , 0.763 , 0.000 ] ,
95 1005211 : [ 1000005 , True , '~Sigmab0 ' , 0 , 0.650 , 0.496 , 0.171 , 0.632 , 0.000 , 0.632 , 0.171 , 0.632 , 0.000 ] ,
96 1005213 : [ 1000005 , True , '~Sigmab*0 ' , 0 , 0.650 , 0.672 , 0.530 , 0.763 , 0.000 , 0.763 , 0.530 , 0.763 , 0.000 ] ,
97 1005223 : [ 1000005 , True , '~Sigmab+ ' , 1 , 0.650 , 0.672 , 0.530 , 0.763 , 0.000 , 0.763 , 0.530 , 0.763 , 0.000 ] ,
98 1005311 : [ 1000005 , True , '~Xib- ' , -1 , 0.825 , 0.691 , 0.498 , 0.833 , 0.000 , 0.828 , 0.498 , 0.833 , 0.000 ] ,
99 1005313 : [ 1000005 , True , '~Xib*- ' , -1 , 0.825 , 0.810 , 0.686 , 0.922 , 0.000 , 0.917 , 0.686 , 0.922 , 0.000 ] ,
100 1005321 : [ 1000005 , True , '~Xib0 ' , 0 , 0.825 , 0.691 , 0.498 , 0.833 , 0.000 , 0.828 , 0.498 , 0.833 , 0.000 ] ,
101 1005323 : [ 1000005 , True , '~Xib*0 ' , 0 , 0.825 , 0.810 , 0.686 , 0.922 , 0.000 , 0.917 , 0.686 , 0.922 , 0.000 ] ,
102 1005333 : [ 1000005 , True , '~Omegab- ' , -1 , 1.000 , 0.952 , 0.863 , 1.095 , 0.000 , 1.075 , 0.863 , 1.095 , 0.000 ] ,
103# Stop R-mesons
104 1000612 : [ 1000006 , True , '~T+ ' , 1 , 0.325 , 0.314 , 0.220 , 0.365 , 0.000 , 0.365 , 0.220 , 0.365 , 0.000 ] ,
105 1000622 : [ 1000006 , True , '~T0 ' , 0 , 0.325 , 0.314 , 0.220 , 0.365 , 0.000 , 0.365 , 0.220 , 0.365 , 0.000 ] ,
106 1000632 : [ 1000006 , True , '~Ts+ ' , 1 , 0.500 , 0.466 , 0.419 , 0.540 , 0.000 , 0.530 , 0.419 , 0.540 , 0.000 ] ,
107 1000642 : [ 1000006 , True , '~Tc0 ' , 0 , 1.500 , 1.630 , 1.550 , 1.710 , 0.000 , 1.700 , 1.550 , 1.710 , 0.000 ] ,
108 1000652 : [ 1000006 , True , '~etat+ ' , 1 , 4.800 , 4.730 , 4.730 , 5.500 , 0.000 , 4.730 , 4.730 , 5.500 , 0.000 ] ,
109# Stop R-baryons
110 1006113 : [ 1000006 , True , '~Sigmat0 ' , 0 , 0.650 , 0.672 , 0.530 , 0.763 , 0.000 , 0.763 , 0.530 , 0.763 , 0.000 ] ,
111 1006211 : [ 1000006 , True , '~Sigmat+ ' , 1 , 0.650 , 0.496 , 0.171 , 0.632 , 0.000 , 0.632 , 0.171 , 0.632 , 0.000 ] ,
112 1006213 : [ 1000006 , True , '~Sigmat*+ ' , 1 , 0.650 , 0.672 , 0.530 , 0.763 , 0.000 , 0.763 , 0.530 , 0.763 , 0.000 ] ,
113 1006223 : [ 1000006 , True , '~Sigmat++ ' , 2 , 0.650 , 0.672 , 0.530 , 0.763 , 0.000 , 0.763 , 0.530 , 0.763 , 0.000 ] ,
114 1006311 : [ 1000006 , True , '~Xit0 ' , 0 , 0.825 , 0.691 , 0.498 , 0.833 , 0.000 , 0.828 , 0.498 , 0.833 , 0.000 ] ,
115 1006313 : [ 1000006 , True , '~Xit*0 ' , 0 , 0.825 , 0.810 , 0.686 , 0.922 , 0.000 , 0.917 , 0.686 , 0.922 , 0.000 ] ,
116 1006321 : [ 1000006 , True , '~Xit+ ' , 1 , 0.825 , 0.691 , 0.498 , 0.833 , 0.000 , 0.828 , 0.498 , 0.833 , 0.000 ] ,
117 1006323 : [ 1000006 , True , '~Xit*+ ' , 1 , 0.825 , 0.810 , 0.686 , 0.922 , 0.000 , 0.917 , 0.686 , 0.922 , 0.000 ] ,
118 1006333 : [ 1000006 , True , '~Omegat0 ' , 0 , 1.000 , 0.952 , 0.863 , 1.095 , 0.000 , 1.075 , 0.863 , 1.095 , 0.000 ] ,
119 }
120
121# Now programmatically calculate the missing spectra
122# These are designed to just flop the rho above or below the gluinoball
123gb_offset = offset_options[1009113][first_mass_set+5]-offset_options[1009113][first_mass_set+1]
124for pid in offset_options:
125 # Skip fundamental SUSY particles and R-glueball
126 if offset_options[pid][0] == 0 or pid == 1000993: continue
127 # Setup #4 to be mass set #1 but with 1009113 matching mass set 5
128 offset_options[pid][first_mass_set+4] = offset_options[pid][first_mass_set+1]+gb_offset
129 # Setup #8 to be mass set #5 but with 1009113 matching mass set 1
130 offset_options[pid][first_mass_set+8] = offset_options[pid][first_mass_set+5]-gb_offset
131
132
133
134def charge( c ):
135 """ Function to return a PDG table formatted charge given an integer charge
136 Input is the charge either as a string or number
137 """
138 n = int(c)
139 if n==0: return ' 0'
140 if n==1: return ' +'
141 if n==2: return '++'
142 if n==-1: return ' -'
143 if n==-2: return '--'
144 raise RuntimeError('Unexpected charge: '+str(n))
145
146
147def get_quarks( y ):
148 """ Function to return a list of quarks in a hadron
149 """
150 x = abs(y)
151 # For stop/sbottom mesons, just the last quark!
152 if '000' in str(x): return str(x)[5:6]
153 # For mesons, just two quarks
154 if '00' in str(x): return str(x)[4:6]
155 # For baryons, three quarks
156 return str(x)[3:6]
157
158
159def is_baryon( x ):
160 # 1000993, gluinoball, is also not a baryon
161 b_n = 0
162 if '009' in str(x): b_n=0 # gluino meson
163 elif '09' in str(x): b_n=1 # gluino baryon
164 elif '0006' in str(x): b_n=0 # stop meson
165 elif '0005' in str(x): b_n=0 # sbottom meson
166 elif '006' in str(x): b_n=1 # stop baryon
167 elif '005' in str(x): b_n=1 # sbottom baryon
168 else: # Otherwise, what on earth was this??
169 raise RuntimeError('is_baryon ERROR Unknown PDG ID: '+str(x))
170 if int(x)<0: return -b_n
171 return b_n
172
173
174def anti_name( x ):
175 """ Function to turn a particle name into an anti-particle name (mostly charge issues)
176 """
177 # These look a little funny to get both + and ++
178 if '*+' in x: return x.replace('*','bar*').replace('+','-')
179 if '*-' in x: return x.replace('*','bar*').replace('-','+')
180 if '++' in x: return x.replace('++','bar--')
181 if '--' in x: return x.replace('--','bar++')
182 if '+' in x: return x.replace('+','bar-')
183 if '-' in x: return x.replace('-','bar+')
184 if '*0' in x: return x.replace('*0','bar*0')
185 return x.replace('0','bar0')
186
187
188def get_gluino_Rhadron_masses(input_file, mass_spectrum=1):
189 """ Function to return a dictionary of PDG IDs and masses based on an input param/SLHA/LHE file
190 First parameter: input file (string or file handle)
191 Second parameter: mass spectrum (enumeration value)
192 """
193 # Expect a string file name or file handle
194 if isinstance(input_file, str):
195 in_file = open(input_file,'r')
196 else:
197 in_file = input_file
198
199 # Expect SLHA file format. Read for mass block, then look for relevant masses, then exit
200 masses = {}
201 mass_block = False
202 for l in in_file.readlines():
203 # Are we entering the mass block?
204 if 'BLOCK' in l.upper().split('#')[0].split() and 'MASS' in l.upper().split('#')[0].split():
205 mass_block = True
206 continue
207 # Otherwise, if we've not yet entered the mass block, keep reading
208 elif not mass_block: continue
209 # If we're past the mass block, then stop reading
210 if 'BLOCK' in l.upper().split('#')[0].split(): break
211 # Skip empty lines, comments, etc
212 if len(l.split('#')[0].split())<2: continue
213 pdg_id = int(l.split('#')[0].split()[0])
214 # Set the input masses
215 if pdg_id in offset_options:
216 mass = float(l.split('#')[0].split()[1])
217 # If we have decoupled the thing, don't include it!
218 if mass > 7e3: continue
219 # Otherwise, it goes in the list
220 masses[pdg_id] = mass
221 # Not an ID we care about otherwise! Skip!
222 # Done reading file; close if it's our responsibility
223 if isinstance(input_file, str):
224 in_file.close()
225
226 # Set the remainder of the masses
227 had_rhadron=False
228 for pid in offset_options:
229 # Skip fundamental particles - they should be read from the input file!
230 if offset_options[pid][0] == 0: continue
231 # Check if the constituent is in there (e.g. skip stop R-hadrons for gluino production)
232 if offset_options[pid][0] not in masses: continue
233 # Check if the mass spectrum is available
234 if mass_spectrum<0 or first_mass_set+mass_spectrum>len(offset_options[pid]):
235 raise RuntimeError("Unknown mass set requested: "+str(mass_spectrum)+" > number of options ("+str(len(offset_options[pid])-first_mass_set+1)+") for PID "+str(pid))
236 # Add 'normal' R-hadron
237 masses[pid] = masses[ offset_options[pid][0] ] + offset_options[pid][first_mass_set+mass_spectrum]
238 # If needed, add anti-R-hadron
239 if offset_options[pid][1]:
240 masses[-pid] = masses[ offset_options[pid][0] ] + offset_options[pid][first_mass_set+mass_spectrum]
241 had_rhadron = True
242
243 # Make sure we generated some R-hadrons
244 if not had_rhadron:
245 raise RuntimeError('No R-hadrons generated!')
246
247 # Return the dictionary
248 return masses
249
250
251def update_PDG_table(input_file, pdg_table, mass_spectrum=1):
252 """ Function to update a PDG table with R-hadron masses
253 First input parameter: input file (string or file handle)
254 Second input parameter: output PDG table (string or file handle)
255 Third input parameter: mass spectrum (enumeration value)
256 Gets R-hadron masses based on get_gluino_Rhadron_masses()
257 """
258 # Check that we had the right output file type
259 # Get the masses that we need
260 masses = get_gluino_Rhadron_masses(input_file,mass_spectrum)
261 # Get the output file ready
262 # Open for appending (assume that's what was done if given a file handle)
263 lines = None
264 if isinstance(pdg_table, str):
265 lines = open(pdg_table).readlines()
266 else:
267 lines = pdg_table.readlines()
268 # Add all our R-hadrons to the table!
269 pdgcodes = []
270 for pid in masses:
271 pdgcodes += [pid]
272 # For the PDG table, we only write positive-signed PDG ID particles
273 if pid<0: continue
274 # Note that we follow the Pythia6 convention of *including* fundamental SUSY particles
275 # The format is VERY specific; needs mass and width (we always set the width to 0)
276 # Mass is in MeV here, rather than GeV as in the dictionary
277 lines.append('M %i %11.7E +0.0E+00 -0.0E+00 %s %s'%(pid,masses[pid]*1000.,offset_options[pid][2],charge(offset_options[pid][3])) + '\n')
278 lines.append('W %i %11.7E +0.0E+00 -0.0E+00 %s %s'%(pid,0.E+00,offset_options[pid][2],charge(offset_options[pid][3])) + '\n')
279
280 update = open('PDGTABLE.MeV', 'w')
281 update.write(''.join(lines))
282 update.close()
283
284 from ExtraParticles.PDGHelpers import updateExtraParticleAcceptList
285 updateExtraParticleAcceptList('G4particle_acceptlist_ExtraParticles.txt', pdgcodes)
286 # Nothing to return
287
288
289def update_particle_table(input_file, particle_table='particles.txt', mass_spectrum=1):
290 """ Function to update a particle table with R-hadron masses
291 First input parameter: input file (string or file handle)
292 Second input parameter: output particle table (string or file handle)
293 Third input parameter: mass spectrum (enumeration value)
294 Gets R-hadron masses based on get_gluino_Rhadron_masses()
295 """
296 # Get the masses that we need
297 masses = get_gluino_Rhadron_masses(input_file,mass_spectrum)
298 # Get the output file ready
299 # Open for appending (assume that's what was done if given a file handle)
300 if isinstance (particle_table, str):
301 out_file = open(particle_table,'a')
302 else:
303 out_file = particle_table
304 # Add all our R-hadrons to the table!
305 # Note that we MUST write the primary first, followed by the compound particles
306 primaries = []
307 extras = []
308 for pid in masses:
309 if offset_options[abs(pid)][0]==0: extras += [pid]
310 elif not offset_options[abs(pid)][0] in primaries: primaries += [ offset_options[abs(pid)][0] ]
311
312 # Rounds per primary
313 for p in primaries:
314 # Note that we follow the old convention of *including* fundamental SUSY particles
315 # The format is VERY specific; needs mass and width (we always set the width to 0)
316 # Mass is in MeV here, rather than GeV as in the dictionary
317 if p>0: out_file.write(' %i %04.3f # %s\n'%(p,masses[p],offset_options[abs(p)][2]))
318 # For the anti-particle, also need the anti-name
319 else: out_file.write(' %i %04.3f # %s\n'%(p,masses[p],anti_name(offset_options[abs(p)][2])))
320 # Now include the secondaries
321 for pid in masses:
322 if offset_options[abs(pid)][0]!=p: continue
323 # Note that we follow the old convention of *including* fundamental SUSY particles
324 # The format is VERY specific; needs mass and width (we always set the width to 0)
325 # Mass is in MeV here, rather than GeV as in the dictionary
326 if pid>0: out_file.write(' %i %04.3f # %s\n'%(pid,masses[pid],offset_options[abs(pid)][2]))
327 # For the anti-particle, also need the anti-name
328 else: out_file.write(' %i %04.3f # %s\n'%(pid,masses[pid],anti_name(offset_options[abs(pid)][2])))
329 # Done with secondaries for this primary
330
331 for p in extras:
332 if p in primaries: continue
333 if p>0: out_file.write(' %i %04.3f # %s\n'%(p,masses[p],offset_options[abs(p)][2]))
334 # For the anti-particle, also need the anti-name
335 else: out_file.write(' %i %04.3f # %s\n'%(p,masses[p],anti_name(offset_options[abs(p)][2])))
336
337 # Done writing all the lines! Clean up if necessary
338 if isinstance(particle_table, str):
339 out_file.close()
340
341 # Nothing to return
342
343
344def get_Pythia8_commands(input_file, mass_spectrum=1):
345 """ Function to return a list of Pythia8 commands to set up an R-hadron mass spectrum.
346 First input parameter: input file (string or file handle)
347 Second input parameter: mass spectrum (enumeration value)
348 """
349 # Get the masses for this configuration
350 masses = get_gluino_Rhadron_masses(input_file,mass_spectrum)
351 # Tell Pythia8 we are going to use our own masses
352 commands = ['RHadrons:setMasses = off']
353
354 # Add commands to set all the masses
355 for pid in masses:
356 # Only set masses for particles (not anti-particles)
357 if pid<0: continue
358 # Actual command takes the form PDGID:m0 = somemass
359 commands += [ str(pid)+':m0 = '+str(masses[pid]) ]
360
361 # All done!
362 return commands
363
364
365def get_interaction_list(input_file, interaction_file='ProcessList.txt', mass_spectrum=1):
366 """ Function to write all possible interactiosn that we need
367 First input parameter: input file (string or file handle)
368 Second input parameter: output PDG table (string or file handle)
369 Third input parameter: mass spectrum (enumeration value)
370 Gets R-hadron masses based on get_gluino_Rhadron_masses()
371 """
372 # Get the masses that we need. Note that we don't really need masses, just PDG IDs
373 masses = get_gluino_Rhadron_masses(input_file,mass_spectrum)
374 # Get the output file ready
375 # Open for appending (assume that's what was done if given a file handle)
376 if isinstance (interaction_file, str):
377 out_file = open(interaction_file,'a')
378 else:
379 out_file = interaction_file
380
381 # Helpful lists to move us along
382 sm_particles = {
383 # Name : Charge , Baryon # , Strangeness
384 'pi0' : [ 0 , 0 , 0 ],
385 'pi+' : [ 1 , 0 , 0 ],
386 'pi-' : [ -1 , 0 , 0 ],
387 'neutron' : [ 0 , 1 , 0 ],
388 'proton' : [ 1 , 1 , 0 ],
389 'kaon0' : [ 0 , 0 , 1 ],
390 'anti_kaon0' : [ 0 , 0 , -1 ],
391 'kaon+' : [ 1 , 0 , 1 ],
392 'kaon-' : [ -1 , 0 , -1 ]
393 }
394 targets = [ 'proton' , 'neutron' ]
395
396 incoming_rhadrons = {}
397 outgoing_rhadrons = {}
398 for pid in masses:
399 # Only for bound states
400 if offset_options[abs(pid)][0]==0: continue
401 # All of them are on the list of incoming RHadrons
402 # Deal with strangeness
403 # Approximation! Bottom number -> -Charm number -> Strangeness
404 # Approximation needed because outgoing SM charms are not treated in G4 at the moment
405 s_number = 0
406 my_q = get_quarks(pid)
407 if '3' in my_q or '4' in my_q or '5' in my_q:
408 if len(my_q)>2:
409 # Gluino R-baryons
410 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')
411 elif len(my_q)>1 and '9' in str(pid):
412 # Gluino R-mesons
413 if my_q in ['33','44','55','35']: s_number=0 # 33, 44, 55, 35 - one is anti-quark, so they cancel
414 # By convention both 43 and 53 have charge +1, which means c-sbar or c-bbar
415 elif my_q in ['43','53']: s_number = 2 if pid>0 else -2
416 # Only one of bottom / charm / strange. Deal with neutral convention first
417 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
418 elif offset_options[abs(pid)][3]==0 and '4' in my_q: s_number=1 if pid<0 else -1
419 # Now charged convention
420 elif '3' in my_q or '5' in my_q: s_number=offset_options[abs(pid)][3]
421 elif '4' in my_q: s_number=-offset_options[abs(pid)][3]
422 elif len(my_q)>1:
423 # Squark R-baryons
424 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')
425 else:
426 # Squark R-mesons
427 s_number = my_q.count('3') - my_q.count('4') + my_q.count('5')
428 s_number = s_number if pid>0 else -s_number
429 else: s_number=0
430 # Build the dictionary
431 pid_name = offset_options[pid][2].strip() if pid>0 else anti_name(offset_options[abs(pid)][2]).strip()
432 charge = offset_options[abs(pid)][3] if pid>0 else -offset_options[abs(pid)][3]
433 incoming_rhadrons[pid_name] = [ charge , is_baryon(pid) , s_number ]
434 # Smaller list of outgoing rhadrons.
435 # No charm or bottom
436 if '4' in my_q or '5' in my_q: continue
437 outgoing_rhadrons[pid_name] = [ charge , is_baryon(pid) , s_number ]
438
439 # Add all our R-hadrons to the table
440 for proj in incoming_rhadrons:
441 # Loop over targets
442 for t in targets:
443 # Loop over possible outgoing R-hadrons
444 for orhad in outgoing_rhadrons:
445 # Possible 2>2 reactions
446 for osm1 in sm_particles:
447 # Check for charge conservation
448 total_charge = incoming_rhadrons[proj][0]+sm_particles[t][0]-outgoing_rhadrons[orhad][0]-sm_particles[osm1][0]
449 # Check for baryon number conservation
450 total_bnumber = incoming_rhadrons[proj][1]+sm_particles[t][1]-outgoing_rhadrons[orhad][1]-sm_particles[osm1][1]
451 # Check for strangeness conservation
452 total_snumber = incoming_rhadrons[proj][2]+sm_particles[t][2]-outgoing_rhadrons[orhad][2]-sm_particles[osm1][2]
453 # Check if it's an allowed reaction
454 if total_charge==0 and total_bnumber==0 and total_snumber==0:
455 out_file.write( ' # '.join([str(proj),str(t),str(orhad),str(osm1)])+'\n' )
456 # Now loop over possible 2>3 reactions
457 for osm2 in sm_particles:
458 # Check for charge conservation
459 total_charge = incoming_rhadrons[proj][0]+sm_particles[t][0]-outgoing_rhadrons[orhad][0]-sm_particles[osm1][0]-sm_particles[osm2][0]
460 # Check for baryon number conservation
461 total_bnumber = incoming_rhadrons[proj][1]+sm_particles[t][1]-outgoing_rhadrons[orhad][1]-sm_particles[osm1][1]-sm_particles[osm2][1]
462 # Check for strangeness conservation
463 total_snumber = incoming_rhadrons[proj][2]+sm_particles[t][2]-outgoing_rhadrons[orhad][2]-sm_particles[osm1][2]-sm_particles[osm2][2]
464 # Check if it's an allowed reaction
465 if total_charge==0 and total_bnumber==0 and total_snumber==0:
466 out_file.write( ' # '.join([str(proj),str(t),str(orhad),str(osm1),str(osm2)])+'\n' )
467 # Wrote out the reaction
468 # Loop over 2>3
469 # Loop over 2>2
470 # Loop over outgoing RHadrons
471 # Loop over targets
472 # Loop over projectiles
473
474 # Done writing all the lines! Clean up if necessary
475 if isinstance(interaction_file, str):
476 out_file.close()
477
478 # Nothing to return
479
480
481def print_masses(spectrum=-1):
482 """ Print the mass spectra.
483 Input parameter: spectrum number. If -1, print all spectra.
484 """
485 for i in sorted(offset_options.keys()):
486 s= str(offset_options[i][2])+' '+str(i)
487 if spectrum<0:
488 from past.builtins import range # Temporary workaround for python3 compatibility use range in CA-based config
489 for j in range(first_mass_set,len(offset_options[i])): s+=' '+str(offset_options[i][j])
490 else:
491 if first_mass_set+spectrum>len(offset_options[i]):
492 raise RuntimeError('Spectrum #'+str(spectrum)+' not known for PID '+str(i))
493 else:
494 s+=' '+str(offset_options[i][spectrum+first_mass_set])
495 print (s)
496 # Done!
std::string replace(std::string s, const std::string &s2, const std::string &s3)
Definition hcg.cxx:310
std::vector< std::string > split(const std::string &s, const std::string &t=":")
Definition hcg.cxx:177
print_masses(spectrum=-1)
update_particle_table(input_file, particle_table='particles.txt', mass_spectrum=1)
get_gluino_Rhadron_masses(input_file, mass_spectrum=1)
get_Pythia8_commands(input_file, mass_spectrum=1)
get_interaction_list(input_file, interaction_file='ProcessList.txt', mass_spectrum=1)
update_PDG_table(input_file, pdg_table, mass_spectrum=1)