All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
gen_crt_frags.py
Go to the documentation of this file.
1 '''Generate the GDML for the SBND CRT geometry.
2 
3 We define a global XML DOM which the function for each volume will append to.
4 Volume-building functions are called hierarchically, so if you call module(),
5 it will construct the module and all the parts that make it up, so you end
6 up with complete GDML for one module.
7 
8 Each physical volume has a corresponding unique logical volume, as required
9 by LArG4 to keep track of energy depositions. The solids, however, can safely
10 be referenced many times, and so are stored only once (using a hash keyed on
11 the the linear dimensions).
12 
13 The output of this code is a file "crt.gdml" which contains the GDML snippets
14 to paste into the full SBND geometry.
15 
16 Created by A. Mastbaum <mastbaum@uchicago.edu>, 2016/10/27
17 Downloaded by C. Hilgenberg <Chris.Hilgenberg@colostate.edu> 2017/10/25
18 Modified by C. Hilgenberg 2017/10/26
19 '''
20 
21 import csv
22 import math
23 import xml.etree.cElementTree as ET
24 from xml.dom import minidom
25 
26 #set true to generate standalone CRT shell, false for normal production
27 testmode = False
28 printModIds = False
29 
30 #################### Parameters #####################
31 #warm vessel (cm), increased as of 7/19/2019 to
32 #reflect full extent of WV profile in detector hall
33 WVWIDTH = 1031.8 #previously 972.0
34 WVHEIGHT = 627.4 #previously 614.0, current value includes support feet, foot height is 30.0cm
35 WVLENGTH = 2268.8 #previously 2209.0
36 WVFOOTELEVATION = 10.16 #height of bottom of WV foot w.r.t. pit floor
37 ISLANDWIDTH = 118.0 #width of ~square WV support feet islands
38 
39 IVLENGTH = 1996.0 #length of interior of cold vessel, the inactive + active LAr volume
40 
41 #OVEROPENZ = 2758.5 - 30.48
42 TOPLEDGERISERTOFLOOR = 991.9 #pit floor to top of ledge riser upon which the overburden blocks sit
43 LEDGERISERHEIGTH = 46.99 #wide flange type, W18 x 71
44 TOPCRTBEAMTOFLOOR = 970.8
45 CRTBEAMSPACING = 92.71 #horizontal center to center spacing between top CRT support beams
46 
47 BOTTOMCRTROLLERHEIGHT = 3.02 #distance between bottom CRT module and pit floor
48 SIDECRTWVOFFSET = 4.13 #set by fiberglass Unistrut standoffs (given as ~4cm, but it's Unitstrut so I assume it was rounded
49 SIDECRTPOSTWIDTH = 4.13 #Unistruit vertical support posts, dimension normal to side CRT plane
50 SIDECRTPOSTSPACING = 4.13 #set by Unistruit bracket shelf
51 SIDECRTSTACKDEPTH = 3*SIDECRTPOSTWIDTH + 2*SIDECRTPOSTSPACING
52 SIDECRTSHELFTHICK = 0.56 #to be verified, could also be 0.64 depending on steel type
53 
54 overlap = 0.0065 # in order remove the overlap between two cut module
55 gap = 0.005 # added gap between two horizontal south wal cut modules
56 
57 # south wall
58 crosstwomodule = 10.7 # set an offset on each horizontal cut module to allow overlap of 10.7 cm
59 separatetwomodule = 0.6 # as two modules are present in one self, we need need to add a gap between two
60  # horizontal module to avoid the overlaps.
61 offsetontopmodule = 48.26 # 19 inches offset (as there are aluminium pipe on the way we can not avoid).
62 verticalmodule = 485.14 # 191 inches, as this module was added later it is not same as other vertical module of 400 cm.
63 
64 #dimensions of top CRT support beams, wide flange W10 x 49
65 #true area is larger than that calculated assuming perfect I shape
66 #91.55 cm^2 vs 92.90 cm^2, adjust web and flange thickness to match true mass
67 #bottom of CRT beam offset vertically + 0.635cm from bottom of ledge riser
68 CRTBEAMLENGTH = 1132.84
69 CRTBEAMHEIGHT = 25.3
70 CRTBEAMWIDTH = 25.4
71 CRTBEAMFLANGETHICK = 1.434 #true 1.422, adjusted to match true area
72 CRTBEAMWEBTHICK = 0.89407 #true 0.86, adjusted to match true area
73 CRTBEAMMASSDENS = 73 #kg/m, density 7849 kg/m^3
74 
75 #strip width, length, thickness (cm)
76 YM = 4.1 #MINOS
77 ZM = 800.0
78 XM = 1.0
79 MINOSSTRAIGHTSNOUT = 37.5
80 MINOSBENDSNOUT = 26.5
81 MINOSSNOUTLENGTH = 0.5*(MINOSSTRAIGHTSNOUT+MINOSBENDSNOUT) #guess, to be verified of length along module that snout, where fiber routing occurs, no scintillator, extends
82 
83 
84 XC = 23.0 #CERN
85 ZC = 184.0
86 YCTOP = 1.0
87 YCBOT = 1.5
88 
89 XD = 5.0 #Double Chooz
90 ZD = 322.5
91 YD = 1.0
92 
93 #Strips per layer
94 NXM = 20
95 NXC = 8
96 NXD = 32
97 
98 # cm padding between strips and module (Al thickness)
99 PADM = 0.05
100 PADC = 0.1
101 PADD = 0.05
102 PADModule = 0.1
103 PADStrip = 0.01
104 PADTagger = 0.001
105 
106 mModL = ZM+2*PADM+2*PADStrip
107 mModW = YM*NXM+2*PADM+(NXM+1)*PADStrip
108 mModH = XM+2*PADM+2*PADStrip
109 cModW = XC*NXC+2*PADC+(NXC+1)*PADStrip #same as length
110 cModH = YCTOP+YCBOT+2*PADC+3*PADStrip
111 dModH = YD*2+2*PADD+3*PADStrip
112 
113 #MINOS mounting
114 LAYERSPACE = 8.27 #MINOS edge-to-edge distance between adjacent layers (cm), prevously 10cm
115 NMODSTACK = 9 #number of lateral MINOS modules in a single layer, single stack
116 NMODSTACKSOUTHY = 10 #number of MINOS modules in a single layer of south vertical side crt,
117 SLIDERSPACE = 18.5+mModH #MINOS center-to-center distance between fixed and sliding stacks' nearest modules (cm), previously 25.0cm
118 STACKOVERLAP = 50.0 #MINOS stack horizontal overlap (cm)
119 SIDECRTROLLOFFSET = 44.29 #offset from outmost extend of WV to center fo rolling stack (E/W sides)
120 SIDECRTNORTHWALLTOFLOOR = 152.2
121 SIDECRTSOUTHWALLLATOFFSET = 1.1*MINOSSNOUTLENGTH
122 
123 #DC mounting
124 DCSPACER=32.6 #foam spacer between DC modules in rows of 5 (strip normal to drift direction) (cm)
125 LONGOFF5=(3*ISLANDWIDTH+481.8)*0.5+181.8
126 LONGOFF2= (ISLANDWIDTH+181.8)*0.5
127 
128 #CERN mounting
129 CERNMODSPACE = 0.2
130 NTOPX=6
131 NTOPZ=14
132 NSLOPELAT=14
133 NSLOPEFRONT=6
134 SLOPEINCLINATION=90.0 #degrees w.r.t. vertical, previously 60 deg
135 CERNROOFL = NTOPZ*cModW+(NTOPZ-1)*CERNMODSPACE
136 
137 #CRT shell
138 SHELLY = 1.1*cModH+TOPCRTBEAMTOFLOOR-BOTTOMCRTROLLERHEIGHT*0.9
139 #MINOS sections positions
140 MINOSSOUTHY = -0.5*SHELLY+0.5*(NMODSTACKSOUTHY*mModW+(NMODSTACKSOUTHY-1)*SIDECRTSHELFTHICK+2*PADTagger)+WVFOOTELEVATION
141 MINOSLATFIXY = MINOSSOUTHY
142 MINOSLATROLLY = MINOSLATFIXY-0.5*mModW+10
143 MINOSLATSOUTHACTIVEOVERHANG = 2*MINOSSNOUTLENGTH
144 MINOSLATSOUTHZ = -0.5*WVLENGTH + 0.5*mModL - 0.5*MINOSLATSOUTHACTIVEOVERHANG
145 MINOSLATNORTHZ = 0.5*WVLENGTH - 0.5*mModL
146 MINOSLATCENTZ = 0.5*(-0.5*MINOSLATSOUTHACTIVEOVERHANG+MINOSLATSOUTHZ+MINOSLATNORTHZ)
147 
148 #DC section positions
149 posDCInDetEncl = (0,-480.135, 0)
150 
151 #CERN sections positions
152 CERNRIMSWVOFFSET = 34.0
153 CERNTOPY = SHELLY*0.5 - 0.6*cModH
154 CERNRIMSY = CERNTOPY - 0.5*cModH - CRTBEAMHEIGHT - 2.54 - 0.5*cModW
155 CERNRIMSZ = -0.5*WVLENGTH - CERNRIMSWVOFFSET
156 CERNRIMNY = CERNRIMSY + 29.0
157 CERNRIMNZ = CERNROOFL + CERNRIMSZ + 0.6*cModH + CRTBEAMSPACING
158 CERNTOPZ = CERNRIMSZ + 0.5*CERNROOFL #roof center assuming south edge of roof aligned with south rim center in z
159 CERNRIMLATX = 0.5*WVWIDTH + 0.5*cModH + 38.0
160 CERNRIMLATY = CERNRIMSY
161 CERNRIMLATZ = CERNTOPZ
162 
163 #previously 1.01*...
164 SHELLZ = 1.01*(0.5*cModH+CERNRIMNZ-min(CERNRIMSZ+0.5*cModH,MINOSLATSOUTHZ-(mModL+MINOSLATSOUTHACTIVEOVERHANG+SIDECRTSOUTHWALLLATOFFSET)*0.5)) #2*(CERNROOFL - 0.5*IVLENGTH - cModW)
165 SHELLWVOFFSET = SHELLZ*0.5/1.01 - WVLENGTH*0.5
166 if CERNRIMSZ+0.5*cModH < MINOSLATSOUTHZ-(mModL+MINOSLATSOUTHACTIVEOVERHANG+SIDECRTSOUTHWALLLATOFFSET)*0.5:
167  SHELLWVOFFSET-= CERNRIMSWVOFFSET
168 else :
169  SHELLWVOFFSET-= MINOSLATSOUTHACTIVEOVERHANG
170 #print('SHELL-WV OFFSET: '+str(SHELLWVOFFSET))
171 
172 #cut MINOS module lengths including snout, index is row number starting from the bottom
173 minosCutModLengthNorth = (256.54, 309.9, 309.9, 508.19, 508.19, 508.19) #6 rows
174 minosCutModLengthSoutheast = (497.84, 497.84, 497.84, 497.84, 497.84, 497.84, 497.84, 497.84, 497.84) #9 rows, 196 inch
175 minosCutModLengthSouthwest = (497.84, 497.84, 497.84, 497.84, 497.84, 497.84, 325.12, 325.12, 325.12) # 6 rows, 196 inch, 3 row 128 inch
176 MINOSNORTHY = -0.5*SHELLY+0.5*(len(minosCutModLengthNorth)*mModW+(len(minosCutModLengthNorth)-1)*SIDECRTSHELFTHICK+2*PADTagger)-PADTagger+SIDECRTNORTHWALLTOFLOOR-BOTTOMCRTROLLERHEIGHT*0.9
177 
178 ########################################################
179 
180 gdml = ET.Element('gdml')
181 if testmode: materials = ET.SubElement(gdml, 'materials')
182 solids = ET.SubElement(gdml, 'solids')
183 structure = ET.SubElement(gdml, 'structure')
184 solids_store = {}
185 modToFeb = dict()
186 feb_id = 0
187 mod_id = -1
188 beam_id = 0
189 nModM = 0
190 nModC = 0
191 nModD = 0
192 
193 def get_mod_id(style='m'):
194  global mod_id
195  global nModM
196  global nModC
197  global nModD
198 
199  mod_id += 1
200  if style == 'm':
201  nModM += 1
202  if style == 'c':
203  nModC += 1
204  if style == 'd':
205  nModD += 1
206  return mod_id
207 
208 def get_mod_id_num():
209  global mod_id
210  return str(mod_id)
211 
212 def beam():
213  '''build one wide flange beam for top CRT support'''
214 
215  global beam_id
216  beam_id += 1
217 
218  xx = str(CRTBEAMWIDTH)
219  yy = str(CRTBEAMHEIGHT)
220  zz = str(CRTBEAMLENGTH)
221  xxsub = str(0.5*(CRTBEAMWIDTH-CRTBEAMWEBTHICK))
222  yysub = str(CRTBEAMHEIGHT-2*CRTBEAMFLANGETHICK)
223  zzsub = zz
224  xsubpos = 0.5*(CRTBEAMWIDTH-0.5*(CRTBEAMWIDTH-CRTBEAMWEBTHICK))
225  area = CRTBEAMWIDTH*CRTBEAMHEIGHT-(CRTBEAMWIDTH-CRTBEAMWEBTHICK)*(CRTBEAMHEIGHT-2*CRTBEAMFLANGETHICK)
226  if beam_id==1: print('modeled - true beam areas (cm^2): '+str(area-92.90304))
227 
228  sname = 'TopCRTSupportBeam'
229  vname = 'vol'+sname+'_'+str(beam_id)
230 
231  if not sname in solids_store:
232  snameext = sname+'_external'
233  snameint = sname+'_internal'
234  snamesub = sname+'_firstsubtraction'
235  sexternal= ET.SubElement(solids, 'box', name=snameext, lunit="cm", x=xx, y=yy, z=zz)
236  sinternal = ET.SubElement(solids, 'box', name=snameint, lunit="cm", x=xxsub, y=yysub, z=zzsub)
237  ssub = ET.SubElement(solids, 'subtraction', name=snamesub)
238  #s = ET.SubElement(solids, 'subtraction', name=sname)
239  ET.SubElement(ssub, 'first', ref=snameext)
240  ET.SubElement(ssub, 'second', ref=snameint)
241  ET.SubElement(ssub, 'position', name='beamsubpos1', unit='cm', x=str(xsubpos), y='0', z='0')
242  s = ET.SubElement(solids, 'subtraction', name=sname)
243  ET.SubElement(s, 'first', ref=snamesub)
244  ET.SubElement(s, 'second', ref=snameint)
245  ET.SubElement(s, 'position', name='beamsubpos2', unit='cm', x=str(-1*xsubpos), y='0', z='0')
246  solids_store[sname] = s
247 
248  else:
249  s = solids_store[sname]
250 
251  v = ET.SubElement(structure, 'volume', name=vname) #Logical volume
252  ET.SubElement(v, 'materialref', ref='STEEL_A992')
253  ET.SubElement(v, 'solidref', ref=sname)
254 
255  return (s,v)
256 
257 def beamVol():
258 
259  padding = 0.01
260  nbeam = 29
261  xx = str(CRTBEAMLENGTH+padding)
262  yy = str(CRTBEAMHEIGHT+padding)
263  zz = (nbeam-1)*CRTBEAMSPACING+CRTBEAMWIDTH+padding
264 
265  sname = 'TopCRTSupportBeamEnclosure'
266  vname = 'vol'+sname
267 
268  beams = []
269  for i in range(nbeam):
270  beams.append(beam())
271 
272  s = ET.SubElement(solids, 'box', name=sname, lunit="cm", x=xx, y=yy, z=str(zz))
273  v = ET.SubElement(structure, 'volume', name=vname)
274  ET.SubElement(v, 'materialref', ref='Air')
275  ET.SubElement(v, 'solidref', ref=sname)
276 
277  for i, (sbeam,vbeam) in enumerate(beams):
278 
279  pv = ET.SubElement(v, 'physvol')
280  ET.SubElement(pv, 'volumeref', ref=vbeam.attrib['name'])
281 
282  dz = str(0.5*(padding-zz+CRTBEAMWIDTH) + i*CRTBEAMSPACING)
283 
284  posname = 'pos' + vbeam.attrib['name']
285  ET.SubElement(pv, 'position', name=posname,unit="cm", x='0', y='0', z=dz)
286  posname = 'rot' + vbeam.attrib['name']
287  ET.SubElement(pv, 'rotation', name=posname, unit="deg", x='0', y='90', z='0')
288 
289  return (s,v)
290 
291 
292 def strip(style="m", modnum=0, stripnum=0, length=0):
293  '''Build one scintillator strip.'''
294 
295  if style=="m":
296  x=XM
297  y=YM
298  if length==0:
299  z=ZM
300  else:
301  z=length
302  name = 'MINOS'
303  if style=="c":
304  name = 'CERN'
305  x=XC
306  if stripnum < NXC:
307  y=YCTOP
308  else:
309  y=YCBOT
310  z=ZC
311  if style=="d":
312  x=XD
313  y=YD
314  z=ZD
315  name = 'DC'
316 
317  xx = str(x)
318  yy = str(y)
319  zz = str(z)
320 
321  sname = 'AuxDetSensitive_' + name
322  vname = 'volAuxDetSensitive_' + name + '_module_'
323 
324  if modnum < 10:
325  vname += '00'
326  elif modnum < 100:
327  vname += '0'
328 
329  vname += str(modnum)+'_'
330  if style=='m' and length!=0:
331  sname+='_cut'+str(int(length))+'_'
332  vname+='cut'+str(int(length))+'_'
333 
334  if style=='c':
335  if stripnum < NXC:
336  sname += '_top_'
337  vname += 'top_'
338  else:
339  sname += '_bot_'
340  vname += 'bot_'
341 
342  vname+='strip_'
343  if stripnum < 10:
344  vname += '0'
345  vname += str(stripnum)
346 
347  if not sname in solids_store:
348  s = ET.SubElement(solids, 'box', name=sname, lunit="cm", x=xx, y=yy, z=zz) #g4 solid
349  solids_store[sname] = s
350 
351  else:
352  s = solids_store[sname]
353 
354  #vname = 'volAuxDetSensitive_' + name
355  v = ET.SubElement(structure, 'volume', name=vname) #Logical volume
356  ET.SubElement(v, 'materialref', ref='Polystyrene')
357  ET.SubElement(v, 'solidref', ref=sname)
358 
359  #print("strip produced!")
360 
361  return s, v #return solid, logical volumes
362 
363 def module(style="c", reg='tt', length=0):
364  '''Build an edge-to-edge array of scintillator strips.'''
365 
366  if style=="m":
367  ny=NXM
368  x=XM
369  y=YM
370  if length==0:
371  z=ZM
372  zz = str(mModL)
373  zzsub = str(mModL-2*PADM)
374  else:
375  if length>ZM:
376  print('WARNING: LENGTH SPECIFIED FOR CUT MINOS MODULE EXCEEDES FULL LENGTH!')
377  z = length
378  zz = str(mModL-ZM+length)
379  zzsub = str(mModL-2*PADM-ZM+length)
380  xx = str(mModH)
381  yy = str(mModW)
382  xxsub = str(mModH-2*PADM)
383  yysub = str(mModW-2*PADM)
384  name = "MINOS"
385 
386  if style=="c":
387  x=XC
388  z=ZC
389  ny=NXC
390  xx = str(cModW)
391  yy = str(cModH)
392  zz = xx
393  xxsub = str(cModW - 2*PADC)
394  yysub = str(cModH-2*PADC)
395  zzsub = xxsub
396  name = "CERN"
397 
398  if style=="d":
399  x=XD
400  y=YD
401  z=ZD
402  ny=NXD
403  xx = str(x*(ny+0.5)+2*PADD+(ny+2)*PADStrip)
404  yy = str(dModH)
405  zz = str(z+2*PADD+2*PADStrip)
406  xxsub = str(x*(ny+0.5)+(ny+2)*PADStrip)
407  yysub = str(dModH-2*PADM)
408  zzsub = str(z+2*PADStrip)
409  name = "DC"
410 
411  modnum = get_mod_id(style)
412  stripnum = 0
413 
414  sname = 'AuxDet_' + name + '_module_'
415  vname = 'vol' + sname
416 
417  if modnum < 10:
418  vname += '00'
419  elif modnum < 100:
420  vname += '0'
421  vname += str(modnum)+'_'
422 
423  if style=='m' and length!=0:
424  sname += 'cut'+str(int(length))
425  vname += 'cut'+str(int(length))+'_'
426 
427  if reg=='tt': vname += 'Top'
428  if reg=='rn': vname += 'RimNorth'
429  if reg=='rs': vname += 'RimSouth'
430  if reg=='rw': vname += 'RimWest'
431  if reg=='re': vname += 'RimEast'
432  if reg=='ss': vname += 'South'
433  if reg=='nn': vname += 'North'
434  if reg=='ws': vname += 'WestSouth'
435  if reg=='wc': vname += 'WestCenter'
436  if reg=='wn': vname += 'WestNorth'
437  if reg=='es': vname += 'EastSouth'
438  if reg=='ec': vname += 'EastCenter'
439  if reg=='en': vname += 'EastNorth'
440  if reg=='bt': vname += 'Bottom'
441 
442  snamein = sname+'_inner'
443 
444  if not sname in solids_store:
445 
446  s = ET.SubElement(solids, 'box', name=sname, lunit="cm", x=xx, y=yy, z=zz)
447  sin = ET.SubElement(solids, 'box', name=snamein, lunit="cm", x=xxsub, y=yysub, z=zzsub)
448 
449  solids_store[sname] = s
450  solids_store[snamein] = sin
451 
452  else:
453  s = solids_store[sname]
454  sin = solids_store[snamein]
455 
456  strips = []
457  strips2 = []
458 
459  #generate strips, top layer for c or d type
460  for i in range(ny):
461  if style=='m' and length>0:
462  strips.append(strip(style, modnum, stripnum,length))
463  else:
464  strips.append(strip(style, modnum, stripnum,0))
465  stripnum += 1
466 
467  #generate bottom layer strips
468  if style=='d' or style=='c':
469  for i in range(ny):
470  strips2.append(strip(style, modnum, stripnum))
471  stripnum += 1
472 
473  vnamein = vname + '_inner'
474  vin = ET.SubElement(structure, 'volume', name=vnamein)
475  ET.SubElement(vin, 'materialref', ref='Air')
476  ET.SubElement(vin, 'solidref', ref=snamein)
477 
478  #place first layer of strips (only layer for m modules)
479  #top layer for c or d modules
480  for i, (es, ev) in enumerate(strips):
481  pv = ET.SubElement(vin, 'physvol')
482  ET.SubElement(pv, 'volumeref', ref=ev.attrib['name'])
483 
484  if style=='m':
485  dy = (2*i - ny + 1)* 0.5 * (y+PADStrip)
486  dx=0
487  if style=='c':
488  dx = (2*i - ny + 1)* 0.5 * (x+PADStrip)
489  dy=0.5*(YCBOT+PADStrip)
490 
491  if style=='d':
492  dy= 0.5*(y+PADStrip)
493  dx=(i - 0.5*ny + 0.25) * (x+PADStrip)
494 
495  posname = 'pos' + ev.attrib['name']
496  ET.SubElement(pv, 'position', name=posname, unit="cm", x=str(dx), y=str(dy), z='0')
497 
498  #place bottom layers
499  if style=='c':
500  for i, (es, ev) in enumerate(strips2):
501  pv = ET.SubElement(vin, 'physvol')
502  ET.SubElement(pv, 'volumeref', ref=ev.attrib['name'])
503 
504  dy= -0.5*(YCTOP+PADStrip)
505  dz=(2*i - ny + 1)* 0.5 * (x+PADStrip)
506 
507  posname = 'pos' + ev.attrib['name']
508  ET.SubElement(pv, 'position', name=posname,
509  unit="cm", x='0', y=str(dy), z=str(dz))
510  posname = 'rot' + ev.attrib['name']
511  ET.SubElement(pv, 'rotation', name=posname, unit="deg", x='0', y='90', z='0')
512 
513 
514  if style=='d':
515  for i, (es, ev) in enumerate(strips2):
516  pv = ET.SubElement(vin, 'physvol')
517  ET.SubElement(pv, 'volumeref', ref=ev.attrib['name'])
518 
519  dy= -0.5*(y+PADStrip)
520  dx=(i - 0.5*ny + 0.75) * (x+PADStrip)
521 
522  posname = 'pos' + ev.attrib['name']
523  ET.SubElement(pv, 'position', name=posname, unit="cm", x=str(dx), y=str(dy), z='0')
524  #DC strips centered (FIX ME!!)
525 
526  v = ET.SubElement(structure, 'volume', name=vname)
527  ET.SubElement(v, 'materialref', ref='ALUMINUM_Al')
528  ET.SubElement(v, 'solidref', ref=sname)
529 
530  pv = ET.SubElement(v, 'physvol')
531  ET.SubElement(pv, 'volumeref', ref=vin.attrib['name'])
532 
533  #print ('module produced!')
534 
535  return s, v
536 
537 #build one stack of MINOS modules for east or west side CRT walls
538 #pos specifies one of 3 stacks on either side: 's, 'c', 'n'
539 def minosSideTagger(side='e', pos='n'):
540  ''' Build a side tagger (1 stacks)
541  '''
542  if side!='e' and side !='w':
543  print('bad value passed to function, minosSideTagger: side='+side)
544  if pos!='n' and pos!='c' and pos != 's':
545  print('bad value passed to function, minosSideTagger: pos='+pos)
546 
547  coords = []
548  modules = []
549  nstack=NMODSTACK
550 
551  if(pos=='c'):
552  nstack-=1
553 
554  z = mModL+2*PADTagger
555  if(pos=='s'):
556  z+=MINOSLATSOUTHACTIVEOVERHANG
557 
558  xx = str(SIDECRTSTACKDEPTH)
559  yy = str(nstack*mModW+(nstack-1)*SIDECRTSHELFTHICK+2*PADTagger)
560  zz = str(z)
561 
562  if pos=='s':
563  xxsub = str(SIDECRTSTACKDEPTH)
564  yysub = str(mModW+PADTagger)
565  zzsub = str(MINOSLATSOUTHACTIVEOVERHANG+PADTagger)
566  xpsub = '0'
567  ypsub = str(0.5*((nstack-1)*(mModW+SIDECRTSHELFTHICK)+PADTagger))
568  zpsub = str(-0.5*(z-MINOSLATSOUTHACTIVEOVERHANG))
569 
570  #loop over stacks
571  for layer in range (2): #6):
572 
573  dx = ((-1)**layer)*0.5*LAYERSPACE
574  if side=='w': dx*= -1
575  dz=0
576  #loop over modules in stack
577  for i in range(nstack):
578  dy = 0.5*(2*i+1-nstack)*(mModW+SIDECRTSHELFTHICK)
579  #if pos=='c':
580  # dy-=0.5*mModW
581  if pos=='s' and i==nstack-1 : dz = 0.5*MINOSLATSOUTHACTIVEOVERHANG
582  elif pos=='s': dz = -0.5*MINOSLATSOUTHACTIVEOVERHANG
583  coords.append((dx,dy,dz))
584 
585  sname = 'tagger_SideLat_'
586  if pos=='c': sname+='Center'
587  if pos=='s':
588  sname+='South'
589  snameint = sname + '_internal'
590  snameext = sname + '_external'
591  if pos=='n': sname+='North'
592 
593  if not sname in solids_store:
594  if pos != 's': stagger = ET.SubElement(solids, 'box', name=sname, lunit="cm", x=xx, y=yy, z=zz)
595  if pos=='s':
596  sext = ET.SubElement(solids, 'box', name=snameext, lunit="cm", x=xx, y=yy, z=zz)
597  sint = ET.SubElement(solids, 'box', name=snameint, lunit="cm", x=xxsub, y=yysub, z=zzsub)
598  stagger = ET.SubElement(solids, 'subtraction', name=sname)
599  ET.SubElement(stagger, 'first', ref=snameext)
600  ET.SubElement(stagger, 'second', ref=snameint)
601  ET.SubElement(stagger, 'position', name='crtsouthtaggersubpos', unit='cm', x=xpsub,y=ypsub,z=zpsub)
602  solids_store[snameext] = sext
603  solids_store[snameint] = sint
604  solids_store[sname] = stagger
605  else:
606  if pos=='s':
607  sext = solids_store[snameext]
608  sint = solids_store[snameint]
609  stagger = solids_store[sname]
610 
611  vname = 'vol_'+ sname+'_'
612 
613  global feb_id
614  global printModIds
615  fmod = 0
616 
617  if printModIds: print('MINOS tagger, '+side+', pos '+pos+' first module: '+str(mod_id+1)+', first FEB: '+str(feb_id+1))
618  feb_id+=2
619 
620  for i in range(len(coords)):
621  if side=='w':
622  if pos=='s':
623  modules.append(module('m','ws'))
624  if i==0: vname+='WestSouth'
625  if pos=='c':
626  modules.append(module('m','wc'))
627  if i==0: vname+='WestCenter'
628  if pos=='n':
629  modules.append(module('m','wn'))
630  if i==0: vname+='WestNorth'
631  if side=='e':
632  if pos=='s':
633  modules.append(module('m','es'))
634  if i==0: vname+='EastSouth'
635  if pos=='c':
636  modules.append(module('m','ec'))
637  if i==0: vname+='EastCenter'
638  if pos=='n':
639  modules.append(module('m','en'))
640  if i==0: vname+='EastNorth'
641  fmod+=1
642  modToFeb[mod_id] = ((feb_id-1,fmod),(feb_id,fmod))
643  if fmod==3:
644  fmod=0
645  if i != len(coords)-1: feb_id+=2
646 
647  if printModIds: print(' last module: '+str(mod_id)+', last FEB: '+str(feb_id))
648  #print('dictionary generated:')
649  #for k in modToFeb.keys():
650  # print('module '+str(k)+': '+str(modToFeb[k]))
651 
652  vtagger = ET.SubElement(structure, 'volume', name=vname)
653  ET.SubElement(vtagger, 'materialref', ref='Air')
654  ET.SubElement(vtagger, 'solidref', ref=sname)
655 
656  #place left side module phy. vol.s
657  for i, (xc,yc,zc) in enumerate(coords):
658 
659  (s,v)=modules[i]
660  pv = ET.SubElement(vtagger, 'physvol')
661  ET.SubElement(pv, 'volumeref', ref=v.attrib['name'])
662 
663  posname = 'pos' + v.attrib['name']
664  ET.SubElement(pv, 'position', name=posname, unit="cm", x=str(xc), y=str(yc), z=str(zc))
665 
666  return stagger, vtagger
667 
668 def minosNorthTagger():
669  ''' Build north MINOS tagger (2 layers cut modules in X-X) on upstream face
670  '''
671 
672  coords = []
673  modules = []
674  ny = len(minosCutModLengthNorth)
675 
676  x = 2*(mModL-ZM+max(minosCutModLengthNorth))+PADM+2*PADTagger
677  y = ny*mModW+(ny-1)*SIDECRTSHELFTHICK+2*PADTagger
678  xx = str(x)
679  yy = str(y)
680  zz = str(SIDECRTSTACKDEPTH)
681 
682 
683  global feb_id
684  global printModIds
685  if printModIds: print('MINOS tagger North, first module: '+str(mod_id+1)+', FEB: '+str(feb_id+1))
686  fmod = 0
687  feb_id+=4
688 
689  #loop over rows starting from bottom going up
690  for row in range(ny):
691  zin = 0.5*LAYERSPACE
692  xleft = 0.5*x - PADTagger - 0.5*(mModL-ZM+minosCutModLengthNorth[row])
693  yrow = -0.5*y + PADTagger + (row+0.5)*mModW + row*SIDECRTSHELFTHICK
694 
695  # (row+0.5)*mModW --> center of each module
696  # the row coordinates are the module center coordinates
697  # four each of the four modules in a given row (two columns of cut modules x 2 layers)
698  # switching layers in north wall: zin -> -zin
699  # switching columns in north wall: xleft -> - xleft
700  rowcoords = ( (xleft,yrow,zin),(-xleft,yrow,zin),(xleft,yrow,-zin),(-xleft,yrow,-zin) )
701  coords.append(rowcoords)
702  rowmodules = []
703  fmod+=4
704 
705  for i in range(4):
706  rowmodules.append(module('m','nn',minosCutModLengthNorth[row]))
707  modToFeb[mod_id] = (feb_id-3+i,fmod/4)
708 
709  modules.append(rowmodules)
710  if fmod==12:
711  fmod=0
712  if row != ny-1: feb_id+=4
713 
714  if printModIds: print(' last module: '+str(mod_id)+', FEB: '+str(feb_id))
715 
716  sname = 'tagger_SideNorth'
717  vname = 'vol_'+ sname
718 
719  stagger = ET.SubElement(solids, 'box', name=sname, lunit="cm", x=xx, y=yy, z=zz)
720  vtagger = ET.SubElement(structure, 'volume', name=vname)
721  ET.SubElement(vtagger, 'materialref', ref='Air')
722  ET.SubElement(vtagger, 'solidref', ref=sname)
723 
724  #print('no. of modules in the north side CRT:', len(coords))
725 
726  #place left side module phy. vol.s
727  for row in range(len(coords)):
728  for mod, (xc,yc,zc) in enumerate(coords[row]):
729 
730  (s,v)=modules[row][mod]
731  pv = ET.SubElement(vtagger, 'physvol')
732  ET.SubElement(pv, 'volumeref', ref=v.attrib['name'])
733 
734  posname = 'pos' + v.attrib['name']
735  ET.SubElement(pv, 'position', name=posname, unit="cm", x=str(xc), y=str(yc), z=str(zc))
736 
737  if xc>0:
738  posname = 'rotplus' + v.attrib['name']
739  ET.SubElement(pv, 'rotation', name=posname, unit="deg", x='0', y='90', z='0')
740  else:
741  posname = 'rotneg' + v.attrib['name']
742  ET.SubElement(pv, 'rotation', name=posname, unit="deg", x='0', y='-90', z='0')
743 
744  return stagger, vtagger
745 
746 def minosSouthTagger():
747  ''' Build front MINOS tagger (2 layers in X-Y) on upstream face
748  '''
749  nmody = 10
750  coords = []
751  modules = []
752 
753  # x = 2*(mModL-ZM+max(minosCutModLengthSoutheast))+overlap+2*PADTagger #0.0065 added to remove overlap problem
754  x = 2*(mModL-ZM+max(minosCutModLengthSoutheast)-crosstwomodule)+overlap+2*PADTagger #0.0065 added to remove overlap problem
755  y = mModL + 2*PADTagger
756  z = SIDECRTSTACKDEPTH+mModH+PADTagger
757 
758  xx = str(x)
759  yy = str(y)
760  zz = str(z)
761  offset = 30.48 # vertical modules are shifted from east side comparing to horizontal module (12 inch)
762  #18 inch (45.72 cm) from west side
763  for i in range(2*nmody+1):
764 
765  if i < nmody: #bottom row
766 
767  dx = -0.5*x + offset + PADTagger + (i+0.5)*mModW + i*PADModule
768  dy = -0.5*y + PADTagger + 0.5*(mModL-0.5*ZM)
769  dz = -0.5*z + PADTagger + 0.5*mModH
770 
771  #print('MINOS tagger South, first module: ', i, ', dx: ', dx,', dy: ',dy,', dz: ', dz)
772 
773  if i == nmody:
774  dx = -0.5*x + offset + PADTagger + (i+0.5)*mModW + i*PADModule
775  dy = -0.5*y + PADTagger + 0.5*(mModL-ZM+verticalmodule)
776  #dy = - 0.5*(mModL-ZM+verticalmodule)
777  dz = -0.5*z + PADTagger + 0.5*mModH
778  #dy = -0.5*y + PADTagger + 0.5*verticalmodule
779  #print('MINOS tagger South, first module: ', i, ', dx: ', dx,', dy: ',dy,', dz: ', dz)
780  #else: #top row
781  if i > nmody:
782  dx = -0.5*x+offset + PADTagger + (i+0.5-nmody-1)*mModW + (i-nmody-1)*PADModule # -ve sign of x means opposite side, switch from west to east side.
783  dy = 0.5*y - PADTagger - 0.5*(mModL - 0.5*ZM)
784  dz = -0.5*z +PADTagger + mModH + 1.5*SIDECRTPOSTWIDTH
785  #print('MINOS tagger South, first module: ', i, ', dx: ', dx,', dy: ',dy,', dz: ', dz)
786  coords.append((dx,dy,dz,1)) #x,y,z,vert=true
787 
788 
789  for i in range(NMODSTACK):
790 
791  dxeast = 0.5*(mModL-ZM+minosCutModLengthSoutheast[i]) + gap - crosstwomodule # 0.005 added for space between two module neck to neck
792  dy = -0.5*y+PADTagger+(i+0.5)*mModW + i*SIDECRTSHELFTHICK
793  dz = 0.5*z - 1.5*SIDECRTPOSTWIDTH - separatetwomodule
794 
795  if i == (NMODSTACK -1):
796  dxeast = 0.5*(mModL-ZM+minosCutModLengthSoutheast[i]) + gap - crosstwomodule - offsetontopmodule - separatetwomodule
797  dy = -0.5*y+PADTagger+(i+0.5)*mModW + i*SIDECRTSHELFTHICK
798  dz = 0.5*z - 1.5*SIDECRTPOSTWIDTH - separatetwomodule
799  coords.append((-dxeast,dy,dz,0)) #x,y,z,vert=false, east side
800 
801  for i in range(NMODSTACK):
802 
803  dxwest = 0.5*(mModL-ZM+minosCutModLengthSouthwest[i]) + gap - crosstwomodule
804  dy = -0.5*y+PADTagger+(i+0.5)*mModW + i*SIDECRTSHELFTHICK
805  dz = 0.5*z - 1.5*SIDECRTPOSTWIDTH - separatetwomodule
806  coords.append((dxwest,dy,dz,0)) #x,y,z,vert=false, west side
807 
808  global feb_id
809  global printModIds
810  if printModIds: print('MINOS tagger South, first module: '+str(mod_id+1)+', FEB: '+str(feb_id+1))
811  fmod = 0
812  feb_id+=1
813 
814 # print('no. of modules in the south side CRT:', len(coords))
815 
816  for i in range(len(coords)):
817  if i<2*nmody+1:
818  #modules.append(module('m','ss',0.5*ZM))
819  if (i == nmody):
820  modules.append(module('m','ss',verticalmodule))
821  else:
822  modules.append(module('m','ss',0.5*ZM))
823  fmod+=1
824  modToFeb[mod_id] = (feb_id,fmod)
825  #print(i, ' mod no: '+str(fmod), 'moduleid: '+str(mod_id)+', FEB: '+str(feb_id))
826  if fmod==3:
827  fmod=0
828  #feb_id+=1
829  if i!=2*nmody: feb_id+=1
830  else: feb_id+=1
831  #print(i,' vertical module: '+str(mod_id)+', FEB: '+str(feb_id))
832  if i >= 2*nmody+1 and i < 2*nmody+1+9 :
833  modules.append(module('m','ss',minosCutModLengthSoutheast[i-(2*nmody+1)]))
834  fmod+=1
835  #modToFeb[mod_id] = ((feb_id-1,fmod),(feb_id,fmod))
836  modToFeb[mod_id] = (feb_id,fmod)
837  if fmod==3:
838  fmod=0
839  if i!= (2*nmody+9): feb_id+=1
840  else: feb_id+=1
841  #print(' horizontal module: '+str(mod_id)+', FEB: '+str(feb_id))
842  if i >= 2*nmody+1+9 :
843  modules.append(module('m','ss',minosCutModLengthSouthwest[i-(2*nmody+1+9)]))
844  fmod+=1
845  #modToFeb[mod_id] = ((feb_id-1,fmod),(feb_id,fmod))
846  modToFeb[mod_id] = (feb_id,fmod)
847  if fmod==3:
848  fmod=0
849  if i!= (2*nmody+18): feb_id+=1
850  #else: feb_id+=1
851  #print(' horizontal module: '+str(mod_id)+', FEB: '+str(feb_id))
852  if printModIds: print(' last module: '+str(mod_id)+', FEB: '+str(feb_id))
853 
854  sname = 'tagger_SideSouth'
855  vname = 'vol_'+ sname
856 
857  stagger = ET.SubElement(solids, 'box', name=sname, lunit="cm", x=xx, y=yy, z=zz)
858  vtagger = ET.SubElement(structure, 'volume', name=vname)
859  ET.SubElement(vtagger, 'materialref', ref='Air')
860  ET.SubElement(vtagger, 'solidref', ref=sname)
861 
862  #place left side module phy. vol.s
863  for i, (xc,yc,zc,r) in enumerate(coords):
864 
865  (s,v)=modules[i]
866  pv = ET.SubElement(vtagger, 'physvol')
867  ET.SubElement(pv, 'volumeref', ref=v.attrib['name'])
868 
869  posname = 'pos' + v.attrib['name']
870  ET.SubElement(pv, 'position', name=posname, unit="cm", x=str(xc), y=str(yc), z=str(zc))
871 
872  posname = 'rot' + v.attrib['name']
873  if r==1 : ET.SubElement(pv, 'rotation', name=posname, unit="deg", x='90', y='0', z='90')
874  if r==0 :
875  if xc>0:
876  posname = 'rotplus' + v.attrib['name']
877  ET.SubElement(pv, 'rotation', name=posname, unit="deg", x='0', y='90', z='0')
878  else:
879  posname = 'rotneg' + v.attrib['name']
880  ET.SubElement(pv, 'rotation', name=posname, unit="deg", x='0', y='-90', z='0')
881 
882  #ET.SubElement(pv, 'rotation', name=posname, unit="deg", x='0', y='90', z='0')
883 
884  return stagger, vtagger
885 
886 def DCTagger():
887  ''' Build bottom tagger
888  '''
889  modwidth = XD*(NXD+0.5)+2*PADD+(NXD+2)*PADStrip
890  xx = str(modwidth*5 + DCSPACER*4 + 2*PADTagger)
891  yy = str(2*(YD+PADD+PADTagger)+3*PADStrip)
892  zz = str(WVLENGTH)
893 
894  coords = []
895  modules = []
896  rot = 0
897  for i in range(14):
898 
899  if (i<5):
900  dx = (2*i-5+1)*0.5*(modwidth+ DCSPACER)
901  dz = -1*LONGOFF5
902  if (i==5 or i==6):
903  dx = (ZD + 2*(PADD+PADStrip))*0.5*(-1)**i
904  dz = -1*LONGOFF2
905  if (i==7 or i==8):
906  dx = (ZD + 2*(PADD+PADStrip))*0.5*(-1)**i
907  dz = LONGOFF2
908  if (i>8):
909  dx = (2*(i-9)-5+1)*0.5*(modwidth+ DCSPACER)
910  dz = LONGOFF5
911 
912  if (i>4 and i<9):
913  rot = 1
914  else :
915  rot = 0
916 
917  coords.append((dx,0,dz,rot))
918 
919  global feb_id
920  global printModIds
921  if printModIds: print('DC tagger, first module: '+str(mod_id+1)+', FEB: '+str(feb_id+1))
922 
923  for i in range(len(coords)):
924  modules.append(module('d','bt'))
925  feb_id+=1
926  modToFeb[mod_id] = (feb_id,1)
927 
928  if printModIds: print(' last module: '+str(mod_id)+', FEB: '+str(feb_id))
929 
930  sname = 'tagger_Bottom'
931  vname = 'vol_'+ sname
932 
933  stagger = ET.SubElement(solids, 'box', name=sname, lunit="cm", x=xx, y=yy, z=zz)
934  vtagger = ET.SubElement(structure, 'volume', name=vname)
935  ET.SubElement(vtagger, 'materialref', ref='Air')
936  ET.SubElement(vtagger, 'solidref', ref=sname)
937 
938  #place left side module phy. vol.s
939  for i, (xc,yc,zc,r) in enumerate(coords):
940 
941  (s,v)=modules[i]
942  pv = ET.SubElement(vtagger, 'physvol')
943  ET.SubElement(pv, 'volumeref', ref=v.attrib['name'])
944 
945  posname = 'pos' + v.attrib['name']
946  ET.SubElement(pv, 'position', name=posname, unit="cm", x=str(xc), y=str(yc), z=str(zc))
947 
948  if r==1 :
949  posname = 'rot' + v.attrib['name']
950  ET.SubElement(pv, 'rotation', name=posname, unit="deg", x='0', y='90', z='0')
951 
952  return stagger, vtagger
953 
954 
955 def cernTopTagger():
956  ''' Build top CERN tagger (1 layer of modules)
957  '''
958  modwidth = ZC + 2*PADC + (NXC+1)*PADStrip
959  xx = str(NTOPX*modwidth+2*PADTagger+(NTOPX-1)*PADModule)
960  yy = str(YCTOP+YCBOT+3*PADStrip+2*PADC+2*PADTagger)
961  zz = str(NTOPZ*modwidth + 2*PADTagger + (NTOPZ-1)*PADModule)
962 
963  coords = []
964  modules = []
965 
966  dz = 0.5*(modwidth+PADModule)*(1 - NTOPZ)
967  dx = 0.5*(modwidth+PADModule)*(1 - NTOPX)
968 
969  for i in range(NTOPX*NTOPZ):
970 
971  coords.append((dx,0,dz))
972 
973  if (i+1)%NTOPZ == 0:
974  dx+= modwidth + PADModule
975  dz = 0.5*(modwidth+PADModule)*(1 - NTOPZ)
976  else: dz+= modwidth + PADModule
977 
978  global feb_id
979  global printModIds
980  if printModIds: print('CERN tagger Top, first module: '+str(mod_id+1)+', FEB: '+str(feb_id+1))
981 
982  for i in range(len(coords)):
983  modules.append(module('c','tt'))
984  feb_id+=1
985  modToFeb[mod_id] = (feb_id,1)
986 
987  if printModIds: print(' last module: '+str(mod_id)+', FEB: '+str(feb_id))
988 
989  sname = 'tagger_Top'
990  vname = 'vol_'+ sname
991 
992  stagger = ET.SubElement(solids, 'box', name=sname, lunit="cm", x=xx, y=yy, z=zz)
993  vtagger = ET.SubElement(structure, 'volume', name=vname)
994  ET.SubElement(vtagger, 'materialref', ref='Air')
995  ET.SubElement(vtagger, 'solidref', ref=sname)
996 
997  #place left side module phy. vol.s
998  for i, (xc,yc,zc) in enumerate(coords):
999 
1000  (s,v)=modules[i]
1001  pv = ET.SubElement(vtagger, 'physvol')
1002  ET.SubElement(pv, 'volumeref', ref=v.attrib['name'])
1003 
1004  posname = 'pos' + v.attrib['name']
1005  ET.SubElement(pv, 'position', name=posname, unit="cm", x=str(xc), y=str(yc), z=str(zc))
1006 
1007  return stagger, vtagger
1008 
1009 def cernLatRimTagger(side='L'):
1010  ''' Build east(side='R') or west(side='L') CERN rim tagger
1011  '''
1012  modwidth = ZC + 2*PADC + (NXC+1)*PADStrip
1013  xx = str(modwidth+2*PADTagger)
1014  yy = str(YCTOP+YCBOT+3*PADStrip+2*PADC+2*PADTagger)
1015  zz = str(NSLOPELAT*modwidth + 2*PADTagger + (NSLOPELAT-1)*PADModule)
1016 
1017  coords = []
1018  modules = []
1019 
1020  dz = 0.5*(modwidth+PADModule)*(1 - NSLOPELAT)
1021 
1022  for i in range(NSLOPELAT):
1023 
1024  coords.append((0,0,dz))
1025 
1026  dz+= modwidth+PADModule
1027 
1028  global feb_id
1029  global printModIds
1030  if printModIds: print('CERN tagger Lat, side '+side+' first module: '+str(mod_id+1)+', FEB: '+str(feb_id+1))
1031 
1032  for i in range(len(coords)):
1033  if side == 'L':
1034  modules.append(module('c','rw'))
1035  if side == 'R':
1036  modules.append(module('c','re'))
1037  feb_id+=1
1038  modToFeb[mod_id] = (feb_id,1)
1039 
1040  if printModIds: print(' last module: '+str(mod_id)+', FEB: '+str(feb_id))
1041 
1042  sname = 'tagger_'
1043  if side == 'L':
1044  sname += 'RimWest'
1045  if side == 'R':
1046  sname += 'RimEast'
1047  vname = 'vol_'+ sname
1048 
1049  stagger = ET.SubElement(solids, 'box', name=sname, lunit="cm", x=xx, y=yy, z=zz)
1050  vtagger = ET.SubElement(structure, 'volume', name=vname)
1051  ET.SubElement(vtagger, 'materialref', ref='Air')
1052  ET.SubElement(vtagger, 'solidref', ref=sname)
1053 
1054  #place left side module phy. vol.s
1055  for i, (xc,yc,zc) in enumerate(coords):
1056 
1057  (s,v)=modules[i]
1058  pv = ET.SubElement(vtagger, 'physvol')
1059  ET.SubElement(pv, 'volumeref', ref=v.attrib['name'])
1060 
1061  posname = 'pos' + v.attrib['name']
1062  ET.SubElement(pv, 'position', name=posname, unit="cm", x=str(xc), y=str(yc), z=str(zc))
1063 
1064  return stagger, vtagger
1065 
1066 
1067 def cernLongRimTagger(side='U'):
1068  ''' Build north(side='D') or south(side='U') CERN rim tagger
1069  '''
1070  modwidth = ZC + 2*PADC + (NXC+1)*PADStrip
1071  xx = str(NSLOPEFRONT*modwidth+2*PADTagger+(NSLOPEFRONT-1)*PADModule)
1072  yy = str(YCTOP+YCBOT+3*PADStrip+2*PADC+2*PADTagger)
1073  zz = str(modwidth + 2*PADTagger)
1074 
1075  coords = []
1076  modules = []
1077 
1078  dx = 0.5*(modwidth+PADModule)*(1 - NSLOPEFRONT)
1079 
1080  for i in range(NSLOPEFRONT):
1081 
1082  coords.append((dx,0,0))
1083  dx+= modwidth+PADModule
1084 
1085  global feb_id
1086  global printModIds
1087  if printModIds: print('CERN tagger Long, side '+side+' first module: '+str(mod_id+1)+', FEB: '+str(feb_id+1))
1088 
1089  for i in range(len(coords)):
1090  if side == 'U':
1091  modules.append(module('c','rs'))
1092  if side == 'D':
1093  modules.append(module('c','rn'))
1094  feb_id+=1
1095  modToFeb[mod_id] = (feb_id,1)
1096 
1097  if printModIds: print(' last module: '+str(mod_id)+', FEB: '+str(feb_id))
1098 
1099  sname = 'tagger_'
1100  if side == 'U':
1101  sname += 'RimSouth'
1102  if side == 'D':
1103  sname += 'RimNorth'
1104  vname = 'vol_'+ sname
1105 
1106  stagger = ET.SubElement(solids, 'box', name=sname, lunit="cm", x=xx, y=yy, z=zz)
1107  vtagger = ET.SubElement(structure, 'volume', name=vname)
1108  ET.SubElement(vtagger, 'materialref', ref='Air')
1109  ET.SubElement(vtagger, 'solidref', ref=sname)
1110 
1111  #place left side module phy. vol.s
1112  for i, (xc,yc,zc) in enumerate(coords):
1113 
1114  (s,v)=modules[i]
1115  pv = ET.SubElement(vtagger, 'physvol')
1116  ET.SubElement(pv, 'volumeref', ref=v.attrib['name'])
1117 
1118  posname = 'pos' + v.attrib['name']
1119  ET.SubElement(pv, 'position', name=posname, unit="cm", x=str(xc), y=str(yc), z=str(zc))
1120 
1121  return stagger, vtagger
1122 
1123 def detectorEnclosure():
1124 
1125  #shell outer and void dimensions
1126  WVPADY = 25
1127  xxint = str(WVWIDTH + 2*SIDECRTWVOFFSET)
1128  yyint = str(WVHEIGHT+1.0+WVPADY)
1129  zzint = str(WVLENGTH+2*SIDECRTWVOFFSET)
1130 
1131  xxext = str(WVWIDTH + 2*SIDECRTROLLOFFSET + 1.1*SIDECRTSTACKDEPTH)
1132  yyext = str(SHELLY)
1133  zzext = str(SHELLZ)
1134 
1135  #generate all of the tagger volumes, CRT modules, and strips
1136  (s,vws) = minosSideTagger('w','s') #MINOS west wall, south stack
1137  (s,vwc) = minosSideTagger('w','c') #MINOS west wall, center stack
1138  (s,vwn) = minosSideTagger('w','n') #MINOS west wall, north stack
1139  (s,ves) = minosSideTagger('e','s') #MINOS east wall, south stack
1140  (s,vec) = minosSideTagger('e','c') #MINOS east wall, center stack
1141  (s,ven) = minosSideTagger('e','n') #MINOS east wall, north stack
1142  (s,vss) = minosSouthTagger()#'U',0,0,0) #MINOS south
1143  (s,vnn) = minosNorthTagger() #MINOS north
1144  (s,vbt) = DCTagger() #DC Bottom
1145  (s,vtt) = cernTopTagger() #CERN top
1146  (s,vrw) = cernLatRimTagger('L') #CERN RimWest
1147  (s,vre) = cernLatRimTagger('R') #CERN RimEast
1148  (s,vrs) = cernLongRimTagger('U') #CERN RimSouth
1149  (s,vrn) = cernLongRimTagger('D') #CERN RimNorth
1150  (s,vbeam) = beamVol()
1151 
1152  #CRT Shell containing all of the tagger volumes and a void to cointain the warm vessel
1153  sname = 'CRT_Shell'
1154  snameext = sname+'_external'
1155  snameint = sname+'_internal'
1156  sexternal = ET.SubElement(solids, 'box', name=snameext, lunit="cm", x=xxext, y=yyext, z=zzext)
1157  sinternal = ET.SubElement(solids, 'box', name=snameint, lunit="cm", x=xxint, y=yyint, z=zzint)
1158  sshell = ET.SubElement(solids, 'subtraction', name=sname)
1159  ET.SubElement(sshell, 'first', ref=snameext)
1160  ET.SubElement(sshell, 'second', ref=snameint)
1161  ET.SubElement(sshell, 'position', name='crtshellsubpos', unit='cm', x='0',y=str(-0.5*SHELLY+WVFOOTELEVATION+0.5*WVHEIGHT+0.5*WVPADY),z=str(-SHELLWVOFFSET))
1162 
1163  vname = 'vol'+sname
1164  vshell = ET.SubElement(structure, 'volume', name=vname)
1165  ET.SubElement(vshell, 'materialref', ref='Air')
1166  ET.SubElement(vshell, 'solidref', ref=sname)
1167 
1168  #crt support beam
1169  pv = ET.SubElement(vshell, 'physvol')
1170  ET.SubElement(pv, 'volumeref', ref=vbeam.attrib['name'])
1171  posname = 'pos' + vbeam.attrib['name']
1172  ET.SubElement(pv, 'position', name=posname, unit="cm", x='0', y=str(CERNTOPY-CRTBEAMHEIGHT*0.5-cModH*0.6), z=str(0))
1173 
1174  #position MINOS west, south
1175  pv = ET.SubElement(vshell, 'physvol')
1176  ET.SubElement(pv, 'volumeref', ref=vws.attrib['name'])
1177 
1178  xc = 0.5* WVWIDTH + SIDECRTWVOFFSET + 0.5*SIDECRTSTACKDEPTH
1179  yc = MINOSLATFIXY
1180  zc = MINOSLATSOUTHZ - SHELLWVOFFSET
1181 
1182  posname = 'pos' + vws.attrib['name']
1183  ET.SubElement(pv, 'position', name=posname, unit="cm", x=str(xc), y=str(yc), z=str(zc))
1184 
1185  #position MINOS west, center
1186  pv = ET.SubElement(vshell, 'physvol')
1187  ET.SubElement(pv, 'volumeref', ref=vwc.attrib['name'])
1188 
1189  xc = 0.5* WVWIDTH + SIDECRTROLLOFFSET
1190  yc = MINOSLATROLLY
1191  zc = MINOSLATCENTZ - SHELLWVOFFSET
1192 
1193  posname = 'pos' + vwc.attrib['name']
1194  ET.SubElement(pv, 'position', name=posname, unit="cm", x=str(xc), y=str(yc), z=str(zc))
1195 
1196  #position MINOS west, north
1197  pv = ET.SubElement(vshell, 'physvol')
1198  ET.SubElement(pv, 'volumeref', ref=vwn.attrib['name'])
1199 
1200  xc = 0.5* WVWIDTH + SIDECRTWVOFFSET + 0.5*SIDECRTSTACKDEPTH
1201  yc = MINOSLATFIXY
1202  zc = MINOSLATNORTHZ - SHELLWVOFFSET
1203 
1204  posname = 'pos' + vwn.attrib['name']
1205  ET.SubElement(pv, 'position', name=posname, unit="cm", x=str(xc), y=str(yc), z=str(zc))
1206 
1207  #Position MINOS east, south
1208  pv = ET.SubElement(vshell, 'physvol')
1209  ET.SubElement(pv, 'volumeref', ref=ves.attrib['name'])
1210 
1211  xc = -0.5* WVWIDTH - SIDECRTWVOFFSET - 0.5*SIDECRTSTACKDEPTH
1212  yc = MINOSLATFIXY
1213  zc = MINOSLATSOUTHZ - SHELLWVOFFSET
1214 
1215  posname = 'pos' + ves.attrib['name']
1216  ET.SubElement(pv, 'position', name=posname, unit="cm", x=str(xc), y=str(yc), z=str(zc))
1217 
1218  #Position MINOS east, center
1219  pv = ET.SubElement(vshell, 'physvol')
1220  ET.SubElement(pv, 'volumeref', ref=vec.attrib['name'])
1221 
1222  xc = -0.5* WVWIDTH - SIDECRTROLLOFFSET
1223  yc = MINOSLATROLLY
1224  zc = MINOSLATCENTZ - SHELLWVOFFSET
1225 
1226  posname = 'pos' + vec.attrib['name']
1227  ET.SubElement(pv, 'position', name=posname, unit="cm", x=str(xc), y=str(yc), z=str(zc))
1228 
1229  #Position MINOS east, north
1230  pv = ET.SubElement(vshell, 'physvol')
1231  ET.SubElement(pv, 'volumeref', ref=ven.attrib['name'])
1232 
1233  xc = -0.5* WVWIDTH - SIDECRTWVOFFSET - 0.5*SIDECRTSTACKDEPTH
1234  yc = MINOSLATFIXY
1235  zc = MINOSLATNORTHZ - SHELLWVOFFSET
1236 
1237  posname = 'pos' + ven.attrib['name']
1238  ET.SubElement(pv, 'position', name=posname, unit="cm", x=str(xc), y=str(yc), z=str(zc))
1239 
1240  #position MINOS north (downstream)
1241  pv = ET.SubElement(vshell, 'physvol')
1242  ET.SubElement(pv, 'volumeref', ref=vnn.attrib['name'])
1243 
1244  xc = 0.0
1245  yc = MINOSNORTHY
1246  zc = 0.5* WVLENGTH + SIDECRTWVOFFSET + 0.5*SIDECRTSTACKDEPTH - SHELLWVOFFSET
1247 
1248  posname = 'pos' + vnn.attrib['name']
1249  ET.SubElement(pv, 'position', name=posname, unit="cm", x=str(xc), y=str(yc), z=str(zc))
1250 
1251  #position MINOS south (upstream)
1252  pv = ET.SubElement(vshell, 'physvol')
1253  ET.SubElement(pv, 'volumeref', ref=vss.attrib['name'])
1254 
1255  xc = 0.0
1256  yc = MINOSSOUTHY
1257  zc = -1*(0.5* WVLENGTH + SIDECRTWVOFFSET + 0.5*(SIDECRTSTACKDEPTH+PADTagger+mModH)) - SHELLWVOFFSET
1258 
1259  posname = 'pos' + vss.attrib['name']
1260  ET.SubElement(pv, 'position', name=posname, unit="cm", x=str(xc), y=str(yc), z=str(zc))
1261 
1262  #position DC Bottom
1263  pv = ET.SubElement(vshell, 'physvol')
1264  ET.SubElement(pv, 'volumeref', ref=vbt.attrib['name'])
1265 
1266  xc = posDCInDetEncl[0]
1267  yc = posDCInDetEncl[1]
1268  zc = posDCInDetEncl[2] - SHELLWVOFFSET
1269 
1270  posname = 'pos' + vbt.attrib['name']
1271  ET.SubElement(pv, 'position', name=posname, unit="cm", x=str(xc), y=str(yc), z=str(zc))
1272 
1273  #position CERN Top
1274  pv = ET.SubElement(vshell, 'physvol')
1275  ET.SubElement(pv, 'volumeref', ref=vtt.attrib['name'])
1276 
1277  xc = 0.0
1278  yc = CERNTOPY
1279  zc = CERNTOPZ - SHELLWVOFFSET
1280 
1281  posname = 'pos' + vtt.attrib['name']
1282  ET.SubElement(pv, 'position', name=posname, unit="cm", x=str(xc), y=str(yc), z=str(zc))
1283 
1284  #position CERN west rim
1285  pv = ET.SubElement(vshell, 'physvol')
1286  ET.SubElement(pv, 'volumeref', ref=vrw.attrib['name'])
1287 
1288  xc = CERNRIMLATX
1289  yc = CERNRIMLATY
1290  zc = CERNRIMLATZ - SHELLWVOFFSET
1291 
1292  posname = 'pos' + vrw.attrib['name']
1293  ET.SubElement(pv, 'position', name=posname, unit="cm", x=str(xc), y=str(yc), z=str(zc))
1294 
1295  posname = 'rot' + vrw.attrib['name']
1296  ET.SubElement(pv, 'rotation', name=posname, unit="deg", x='0', y='0', z=str(SLOPEINCLINATION))
1297 
1298  #position CERN east rim
1299  pv = ET.SubElement(vshell, 'physvol')
1300  ET.SubElement(pv, 'volumeref', ref=vre.attrib['name'])
1301 
1302  xc = -1*CERNRIMLATX
1303  yc = CERNRIMLATY
1304  zc = CERNRIMLATZ - SHELLWVOFFSET
1305 
1306  posname = 'pos' + vre.attrib['name']
1307  ET.SubElement(pv, 'position', name=posname, unit="cm", x=str(xc), y=str(yc), z=str(zc))
1308 
1309  posname = 'rot' + vre.attrib['name']
1310  ET.SubElement(pv, 'rotation', name=posname, unit="deg", x='0', y='0', z=str(-1*SLOPEINCLINATION))
1311 
1312  #position CERN south rim
1313  pv = ET.SubElement(vshell, 'physvol')
1314  ET.SubElement(pv, 'volumeref', ref=vrs.attrib['name'])
1315 
1316  xc = 0.0
1317  yc = CERNRIMSY
1318  zc = CERNRIMSZ - SHELLWVOFFSET
1319 
1320  posname = 'pos' + vrs.attrib['name']
1321  ET.SubElement(pv, 'position', name=posname, unit="cm", x=str(xc), y=str(yc), z=str(zc))
1322 
1323  posname = 'rot' + vrs.attrib['name']
1324  ET.SubElement(pv, 'rotation', name=posname, unit="deg", x=str(SLOPEINCLINATION), y='0', z='0')
1325 
1326  #position CERN north rim
1327  pv = ET.SubElement(vshell, 'physvol')
1328  ET.SubElement(pv, 'volumeref', ref=vrn.attrib['name'])
1329 
1330  xc = 0.0
1331  yc = CERNRIMNY
1332  zc = CERNRIMNZ - SHELLWVOFFSET
1333 
1334  posname = 'pos' + vrn.attrib['name']
1335  ET.SubElement(pv, 'position', name=posname, unit="cm", x=str(xc), y=str(yc), z=str(zc))
1336 
1337  posname = 'rot' + vrn.attrib['name']
1338  ET.SubElement(pv, 'rotation', name=posname, unit="deg", x=str(-1*SLOPEINCLINATION), y='0', z='0')
1339 
1340  return sshell,vshell
1341 
1342 #######################################################################
1343 # now if test mode generate materials, CRT shell, world, gdml header
1344 # else just generate CRT shell for use with master geometry generator script
1345 
1346 if testmode:
1347  m = ET.SubElement(materials, 'element', name='aluminum', formula='Al', Z='13')
1348  ET.SubElement(m, 'atom', value='26.9815')
1349  m = ET.SubElement(materials, 'element', name='nitrogen', formula='N', Z='7')
1350  ET.SubElement(m, 'atom', value='14.0067')
1351  m = ET.SubElement(materials, 'element', name='oxygen', formula='O', Z='8')
1352  ET.SubElement(m, 'atom', value='15.999')
1353  m = ET.SubElement(materials, 'element', name='argon', formula='Ar', Z='18')
1354  ET.SubElement(m, 'atom', value='39.9480')
1355  m = ET.SubElement(materials, 'element', name='hydrogen', formula='H', Z='1')
1356  ET.SubElement(m, 'atom', value='1.0079')
1357  m = ET.SubElement(materials, 'element', name='carbon', formula='C', Z='6')
1358  ET.SubElement(m, 'atom', value='12.0107')
1359  m = ET.SubElement(materials, 'element', name='iron', formula='Fe', Z='26')
1360  ET.SubElement(m, 'atom', value='55.8450')
1361  m = ET.SubElement(materials, 'element', name='niobium', formula='Nb', Z='41') #add
1362  ET.SubElement(m, 'atom', value='92.90637')
1363  m = ET.SubElement(materials, 'element', name='copper', formula='Cu', Z='29') #add
1364  ET.SubElement(m, 'atom', value='63.5463')
1365  m = ET.SubElement(materials, 'element', name='manganese', formula='Mn', Z='25') #add
1366  ET.SubElement(m, 'atom', value='54.938043')
1367  m = ET.SubElement(materials, 'element', name='molybdenum', formula='Mo', Z='42') #add
1368  ET.SubElement(m, 'atom', value='95.951')
1369  m = ET.SubElement(materials, 'element', name='nickel', formula='Ni', Z='28')
1370  ET.SubElement(m, 'atom', value='58.6934')
1371  m = ET.SubElement(materials, 'element', name='phosphorus', formula='P', Z='15')
1372  ET.SubElement(m, 'atom', value='30.973')
1373  m = ET.SubElement(materials, 'element', name='silicon', formula='Si', Z='14')
1374  ET.SubElement(m, 'atom', value='28.0855')
1375  m = ET.SubElement(materials, 'element', name='sulfur', formula='S', Z='16')
1376  ET.SubElement(m, 'atom', value='32.065')
1377  m = ET.SubElement(materials, 'element', name='vanadium', formula='V', Z='23')
1378  ET.SubElement(m, 'atom', value='50.94151')
1379 
1380 
1381  m = ET.SubElement(materials, 'material', name='ALUMINUM_Al', formula='ALUMINUM_Al')
1382  ET.SubElement(m, 'D', value='2.6990', unit='g/cm3')
1383  ET.SubElement(m, 'fraction', n='1.000', ref='aluminum')
1384 
1385  m = ET.SubElement(materials, 'material', name='Air')
1386  ET.SubElement(m, 'D', value='0.001205', unit='g/cm3')
1387  ET.SubElement(m, 'fraction', n='0.781154', ref='nitrogen')
1388  ET.SubElement(m, 'fraction', n='0.209476', ref='oxygen')
1389  ET.SubElement(m, 'fraction', n='0.00934', ref='argon')
1390 
1391  m = ET.SubElement(materials, 'material', name='Polystyrene')
1392  ET.SubElement(m, 'D', value='1.19', unit='g/cm3')
1393  ET.SubElement(m, 'fraction', n='0.077418', ref='hydrogen')
1394  ET.SubElement(m, 'fraction', n='0.922582', ref='carbon')
1395 
1396  m = ET.SubElement(materials, 'material', name='STEEL_A992')
1397  ET.SubElement(m, 'D', value='7.85', unit='g/cm3')
1398  ET.SubElement(m, 'fraction', n='0.0022', ref='carbon')
1399  ET.SubElement(m, 'fraction', n='0.0004', ref='niobium')
1400  ET.SubElement(m, 'fraction', n='0.005', ref='copper')
1401  ET.SubElement(m, 'fraction', n='0.01', ref='manganese')
1402  ET.SubElement(m, 'fraction', n='0.0014', ref='molybdenum')
1403  ET.SubElement(m, 'fraction', n='0.0044', ref='nickel')
1404  ET.SubElement(m, 'fraction', n='0.00034', ref='phosphorus')
1405  ET.SubElement(m, 'fraction', n='0.0039', ref='silicon')
1406  ET.SubElement(m, 'fraction', n='0.00044', ref='sulfur')
1407  ET.SubElement(m, 'fraction', n='0.001', ref='vanadium')
1408  ET.SubElement(m, 'fraction', n='0.97092', ref='iron')
1409 
1410 #crt shell volume
1411 (s,v) = detectorEnclosure()
1412 #(s,v) = module()
1413 
1414 if testmode:
1415  ws = ET.SubElement(solids, 'box', name='World', lunit="cm", x='1500', y='1500', z='3000')
1416  w = ET.SubElement(structure, 'volume', name='volWorld')
1417  ET.SubElement(w, 'materialref', ref='Air')
1418  ET.SubElement(w, 'solidref', ref='World')
1419  pv = ET.SubElement(w, 'physvol')
1420  ET.SubElement(pv, 'volumeref', ref=v.attrib['name'])
1421  posname = 'pos' + v.attrib['name']
1422  ET.SubElement(pv, 'position', name=posname, unit="cm", x='0', y='0', z='0')
1423  setup = ET.SubElement(gdml, 'setup', name='Default', version='1.0')
1424  ET.SubElement(setup, 'world', ref='volWorld')
1425 
1426 #no. modules generated for each subsystem
1427 print('MINOS modules generated: '+str(nModM))
1428 print('CERN modules generated: '+str(nModC))
1429 print('DblCh modules generated: '+str(nModD))
1430 
1431 #write to file
1432 with open('icarus_crt.gdml', 'w') as f:
1433  f.write(minidom.parseString(ET.tostring(gdml)).toprettyxml(indent='\t'))
1434 
1435 print('Writing dictionary to file...')
1436 with open('feb_map.txt','w') as f:
1437  writer = csv.writer(f,delimiter=',',quotechar='"',quoting=csv.QUOTE_MINIMAL)
1438  for mod in modToFeb.keys():
1439  tmp = modToFeb[mod]
1440  if type(tmp[0])==int:
1441  writer.writerow([str(mod),str(tmp[0]),str(tmp[1])])
1442  else:
1443  writer.writerow([str(mod),str(tmp[0][0]),str(tmp[0][1]),str(tmp[1][0]),str(tmp[1][1])])
if
then if[["$THISISATEST"==1]]
Definition: neoSmazza.sh:95
gen_crt_frags.beam
def beam
Definition: gen_crt_frags.py:212
gen_crt_frags.beamVol
def beamVol
Definition: gen_crt_frags.py:257
print
do one_file $F done echo for F in find $TOP name CMakeLists txt print
Definition: update_sources.sh:94
gen_crt_frags.detectorEnclosure
def detectorEnclosure
Definition: gen_crt_frags.py:1123
gen_crt_frags.minosSouthTagger
def minosSouthTagger
Definition: gen_crt_frags.py:746
gen_crt_frags.get_mod_id
def get_mod_id
Definition: gen_crt_frags.py:193
util::enumerate
auto enumerate(Iterables &&...iterables)
Range-for loop helper tracking the number of iteration.
Definition: enumerate.h:69
gen_crt_frags.cernTopTagger
def cernTopTagger
Definition: gen_crt_frags.py:955
gen_crt_frags.cernLatRimTagger
def cernLatRimTagger
Definition: gen_crt_frags.py:1009
gen_crt_frags.cernLongRimTagger
def cernLongRimTagger
Definition: gen_crt_frags.py:1067
gen_crt_frags.strip
def strip
Definition: gen_crt_frags.py:292
gen_crt_frags.DCTagger
def DCTagger
Definition: gen_crt_frags.py:886
gen_crt_frags.minosNorthTagger
def minosNorthTagger
Definition: gen_crt_frags.py:668
gen_crt_frags.minosSideTagger
def minosSideTagger
Definition: gen_crt_frags.py:539
gen_crt_frags.get_mod_id_num
def get_mod_id_num
Definition: gen_crt_frags.py:208
gen_crt_frags.module
def module
Definition: gen_crt_frags.py:363
open
open(RACETRACK) or die("Could not open file $RACETRACK for writing")
Generated by   doxygen 1.8.5