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