ATLAS Offline Software
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 """
8 The offset options. Dictionary of PDG IDs with offsets.
9 Anti-particles by definition have the same mass as standard particles.
10 Columns 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 
24 The list of possible R-hadrons comes from the Pythia8 code, in src/RHadrons.cc (at the top)
25 """
26 
27 first_mass_set = 4
28 offset_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
123 gb_offset = offset_options[1009113][first_mass_set+5]-offset_options[1009113][first_mass_set+1]
124 for 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 
134 def 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 
147 def 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 
159 def 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 
174 def 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 
188 def 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 
251 def 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 
289 def 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 
344 def 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 
365 def 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 
481 def 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!
DerivationFramework::TriggerMatchingUtils::sorted
std::vector< typename R::value_type > sorted(const R &r, PROJ proj={})
Helper function to create a sorted vector from an unsorted range.
replace
std::string replace(std::string s, const std::string &s2, const std::string &s3)
Definition: hcg.cxx:307
PDGHelpers.updateExtraParticleAcceptList
def updateExtraParticleAcceptList(listName='G4particle_acceptlist_ExtraParticles.txt', pdgcodes=[])
Definition: PDGHelpers.py:57
RHadronMasses.update_particle_table
def update_particle_table(input_file, particle_table='particles.txt', mass_spectrum=1)
Definition: RHadronMasses.py:289
RHadronMasses.get_gluino_Rhadron_masses
def get_gluino_Rhadron_masses(input_file, mass_spectrum=1)
Definition: RHadronMasses.py:188
RHadronMasses.get_interaction_list
def get_interaction_list(input_file, interaction_file='ProcessList.txt', mass_spectrum=1)
Definition: RHadronMasses.py:365
RHadronMasses.anti_name
def anti_name(x)
Definition: RHadronMasses.py:174
RHadronMasses.is_baryon
def is_baryon(x)
Definition: RHadronMasses.py:159
plotBeamSpotVxVal.range
range
Definition: plotBeamSpotVxVal.py:194
RHadronMasses.get_quarks
def get_quarks(y)
Definition: RHadronMasses.py:147
RHadronMasses.charge
def charge(c)
Definition: RHadronMasses.py:134
TCS::join
std::string join(const std::vector< std::string > &v, const char c=',')
Definition: Trigger/TrigT1/L1Topo/L1TopoCommon/Root/StringUtils.cxx:10
Trk::open
@ open
Definition: BinningType.h:40
RHadronMasses.get_Pythia8_commands
def get_Pythia8_commands(input_file, mass_spectrum=1)
Definition: RHadronMasses.py:344
python.CaloAddPedShiftConfig.int
int
Definition: CaloAddPedShiftConfig.py:45
RHadronMasses.print_masses
def print_masses(spectrum=-1)
Definition: RHadronMasses.py:481
RHadronMasses.update_PDG_table
def update_PDG_table(input_file, pdg_table, mass_spectrum=1)
Definition: RHadronMasses.py:251
str
Definition: BTagTrackIpAccessor.cxx:11
Trk::split
@ split
Definition: LayerMaterialProperties.h:38
python.LArMinBiasAlgConfig.float
float
Definition: LArMinBiasAlgConfig.py:65
Generated on Tue Sep 30 2025 21:16:53 for ATLAS Offline Software by   doxygen 1.8.18