generate_lbne10kT.pl
Go to the documentation of this file.
1 #!/usr/bin/perl
2 
3 # contact tylerdalion@gmail.com for any GDML/generate questions
4 # I would love to help!
5 
6 # This is essentially the same script as the one that generates 34kt,
7 # just with a different set of parameters/dimensions. Eventually these
8 # differing parameters will be handled by passing different xml inputs
9 # to a single script, thereby setting the parameters.
10 
11 # For now, the differences are:
12 #--------------------- 34kt ------------------ 10kt -----------------
13 # MaxDrift 385 228
14 # nAPALong 18 10
15 
16 # The only structural differnce in this script is that 10kt places
17 # cryostats side by side, whereas 34kt is end to end.
18 
19 
20 # Each subroutine generates a fragment GDML file, and the last subroutine
21 # creates an XML file that make_gdml.pl will use to appropriately arrange
22 # the fragment GDML files to create the final desired DUNE GDML file,
23 # to be named by make_gdml output command
24 
25 # If you are playing with different geometries, you can use the
26 # suffix command to help organize your work.
27 
28 use Math::Trig;
29 use Getopt::Long;
30 use Math::BigFloat;
31 Math::BigFloat->precision(-16);
32 
33 GetOptions( "help|h" => \$help,
34  "suffix|s:s" => \$suffix,
35  "output|o:s" => \$output,
36  "wires|w:s" => \$wires,
37  "helpcube|c" => \$helpcube);
38 
39 if ( defined $help )
40 {
41  # If the user requested help, print the usage notes and exit.
42  usage();
43  exit;
44 }
45 
46 if ( ! defined $suffix )
47 {
48  # The user didn't supply a suffix, so append nothing to the file
49  # names.
50  $suffix = "";
51 }
52 else
53 {
54  # Otherwise, stick a "-" before the suffix, so that a suffix of
55  # "test" applied to filename.gdml becomes "filename-test.gdml".
56  $suffix = "-" . $suffix;
57 }
58 
59 
60 #++++++++++++++++++++++++ Begin defining variables +++++++++++++++++++++++++
61 
62 # Define detector geometry variables - later to be put in a parameters
63 # XML file to be parsed as an input?
64 
65 # set wires on to be the default, unless given an input by the user
66 $wires_on = 1; # 1=on, 0=off
67 if (defined $wires)
68 {
69 $wires_on = $wires
70 }
71 
72 $tpc_on=1;
73 $inch = 2.54;
74 
75 
76 ###### ALL PARAMETERS FROM DocDb-3383 #######
77 
78 
79 ##################################################################
80 ##################### wire plane parameters ######################
81 
82 $UWirePitch = 0.49;
83 $VWirePitch = 0.5;
84 $XWirePitch = 0.45;
85 
86 #$UAngle = 45.7;
87 #$VAngle = 44.3;
88 $UAngle = 36;
89 $VAngle = 36;
90 
91 $SinUAngle = sin( deg2rad($UAngle) );
92 $CosUAngle = cos( deg2rad($UAngle) );
93 $TanUAngle = tan( deg2rad($UAngle) );
94 
95 $SinVAngle = sin( deg2rad($VAngle) );
96 $CosVAngle = cos( deg2rad($VAngle) );
97 $TanVAngle = tan( deg2rad($VAngle) );
98 
99 $UWire_yint = $UWirePitch/$SinUAngle;
100 $UWire_zint = $UWirePitch/$CosUAngle;
101 
102 $VWire_yint = $VWirePitch/$SinVAngle;
103 $VWire_zint = $VWirePitch/$CosVAngle;
104 
105 $TPCWireThickness = 0.015;
106 
107 $TPCWirePlaneThickness = $TPCWireThickness;
108 #height and length defined lower
109 
110 
111 
112 
113 
114 ###########################################################################
115 ############## modular APA dimension and spacing parameters ###############
116 
117 $nCryos = 2;
118 $nAPAWide = 3;
119 $nAPAHigh = 2;
120 $nAPALong = 10;
121 
122 # 4 APA testing geometry
123 #$nCryos = 1;
124 #$nAPAWide = 1;
125 #$nAPAHigh = 2;
126 #$nAPALong = 2;
127 
128 $CPAThickness = 5.1;
129 $APAFrame_x = 2*$inch; # this does not include the wire spacing
130 $APAWirePlaneSpacing = 0.476; # spacing between all of the wire planes (u, v, and x)
131 $MaxDrift = 228;
132  #MaxDrift is the distance form the edge of the CPA to the edge of the first wire plane
133  #TODO: the implementation and value of MaxDrift will have to be re-evaluated if adding the grid plane.
134 
135 $APALongGap = 1.5; # separation between APAs along the incident beam axis
136 $APAVerticalGap = 2.5; # separation between APAs along the vertical axis
137 
138 #Cryostat space with LAr outside of entire fiducial volume
139 $SpaceCPAtoCryoWall = 85;
140 $SpaceAPAToFloor = 50;
141 $SpaceAPAToTopLAr = 50;
142 $UpstreamLArPadding = 250;
143 $DownstreamLArPadding = 50;
144 
145 
146 $APAWidth = 2*$MaxDrift
147  + 4*$APAWirePlaneSpacing
148  + 2*$TPCWirePlaneThickness
149  + $APAFrame_x;
150  # This is the distance edge of CPA to edge of CPA.
151  # MaxDrift covers the spacing involving the grid plane.
152  # Unlike the other dimensions, this does not describe
153  # a physical dimension of an APA, but rather is used
154  # to position the "module" one would think of as having
155  # the APA and the corresponding drift distances.
156 
157 $APAHeight = 700; # doesn't include front-end boards or hanger posts
158 $APALength = 252; # doesn't include the dead space on the sides of APA frame for wrapping wires
159 
160 
161 
162 
163 
164 
165 ####################################################################
166 ################# APA Frame and Paddle Dimensions ##################
167 
168 $APAFrameZSide_x = $APAFrame_x;
169 $APAFrameZSide_y = 4*$inch;
170 $APAFrameZSide_z = $APALength;
171 
172 
173 $LightPaddleWidth = 0.476;
174 $LightPaddleHeight = 4*$inch;
175 $LightPaddleLength = 225-0.001;
176 $nLightPaddlesPerAPA = 10; # 10, or 20 for double coverage (for now)
177 $PaddleYInterval = (2*$APAHeight+$APAVerticalGap-$LightPaddleHeight-2*$APAFrameZSide_y)
178  /(2*$nLightPaddlesPerAPA-1);
179 $FrameToPaddleSpace = ($PaddleYInterval-$APAVerticalGap)/2;
180 
181 $SiPM_z = 0;
182 
183 # $PaddleYInterval is defined so that the center-to-center distance in the
184 # y direction between paddles is uniform between vertically stacked APAs.
185 # $FrameToPaddleSpace is from the BOTTOM of the APA frame (4" in y direction)
186 # to the CENTER of a paddle, including the 4" part of the frame. This variable's
187 # primary purpose is to position the lowest paddle in each APA.
188 
189 
190 $APAFrameZSide_x = $APAFrame_x;
191 $APAFrameZSide_y = 4*$inch;
192 $APAFrameZSide_z = $APALength;
193 
194 $APAFrameYSide_x = $APAFrame_x;
195 $APAFrameYSide_y = $APAHeight-2*$APAFrameZSide_y;
196 $APAFrameYSide_z = 4*$inch;
197 
198 # Two outer Y supports will sandwich the inner Y supports and light paddles
199 $APAFrameYOuterSupport_x = ($APAFrame_x-$LightPaddleWidth)/2;
200 $APAFrameYOuterSupport_y = $APAHeight-2*$APAFrameZSide_y;
201 $APAFrameYOuterSupport_z = 4*$inch;
202 
203 $APAFrameYInnerSupport_x = $LightPaddleWidth;
204 $APAFrameYInnerSupport_y = $PaddleYInterval-$LightPaddleHeight;
205 $APAFrameYInnerSupport_z = 4*$inch;
206 
207 $APAFrameZHalfSupport_x = $APAFrame_x;
208 $APAFrameZHalfSupport_y = 2*$inch;
209 $APAFrameZHalfSupport_z = $APAFrameZSide_z
210  - 2*$APAFrameYSide_z
211  - $APAFrameYSupport_z;
212 
213 
214 $EdgeFrameSteelThickness = 0.12*$inch;
215 $InnerFrameSteelThickness = 0.062*$inch;
216 
217 
218 
219 
220 
221 
222 
223 ##############################################################
224 ############## Cryo and TPC relevant dimensions #############
225 
226  # TODO: These fiducial parameters were useful in placement within the Cryostat,
227  # but now that it is only symmetric in the x direction, the only placement where
228  # these parameters are useful is in the x direction. possibly deprecate so as
229  # not to be misleading or overly complicated here.
230 $FiducialWidth = $APAWidth*$nAPAWide + $CPAThickness*($nAPAWide+1);
231 $FiducialHeight = $APAHeight*$nAPAHigh + ($nAPAHigh - 1)*$APAVerticalGap;
232 $FiducialLength = $APALength*$nAPALong + ($nAPALong - 1)*$APALongGap;
233 
234 $SteelThickness = 0.5*$inch; #half inch
235 $HeightGaseousAr = 50;
236 
237 
238 $ArgonWidth = $FiducialWidth
239  + 2*$SpaceCPAtoCryoWall;
240 
241 $ArgonHeight = $FiducialHeight
242  + $SpaceAPAToFloor + $SpaceAPAToTopLAr
243  + $HeightGaseousAr;
244  # both liquid AND gaseous argon
245 
246 $ArgonLength = $FiducialLength
247  + $UpstreamLArPadding + $DownstreamLArPadding;
248 
249 $CryostatWidth = $ArgonWidth + 2*$SteelThickness;
250 $CryostatHeight = $ArgonHeight + 2*$SteelThickness;
251 $CryostatLength = $ArgonLength + 2*$SteelThickness;
252 
253 $TPCWidth = ($APAWidth-$APAFrame_x)/2;
254  # this distance is the distance from edge of
255  # the APA frame to the edge of the CPA
256 $TPCHeight = $APAHeight + $APAVerticalGap;
257 $TPCLength = $APALength + $APALongGap;
258 
259 
260  # $TPCWirePlaneThickness now defined higher up
261 $TPCWirePlaneHeight = $APAHeight;
262  # the wire plane region is the full height of the APA since
263  # the previous number doesn't have the front-end boards, etc.
264 $TPCWirePlaneLength = $APALength;
265  # the APA Length doesn't have the spacing on between the two
266  # APAs so the Wire Plane Length is the full length
267 
268 
269 ### TODO: MAY NEED ADJUSTMENT:
270 #TPC Active Variables -- apply cuts here
271 
272 $TPCActiveWidth = $TPCWidth-(3*$APAWirePlaneSpacing);
273 $TPCActiveHeight = $TPCWirePlaneHeight;
274 $TPCActiveLength = $TPCWirePlaneLength;
275 
276 $posTPCActive_X = $TPCWidth/2-$TPCActiveWidth/2;
277 $posTPCActive_Y = 0;
278 $posTPCActive_Z = 0;
279 
280 
281 
282 
283 
284 
285 
286 
287 ##################################################################
288 ############## DetEnc and World relevant parameters #############
289 
290 $ArToAr = 300;
291  # x distance between the LAr in each side by side cryo
292 $ConcretePadding = 50;
293 $FoamPadding = 80;
294 $TotalPadding = $ConcretePadding+$FoamPadding;
295 $DetEncWidth = 2*$CryostatWidth+2*$TotalPadding+2*$FoamPadding + $ArToAr;
296 $DetEncHeight = $CryostatHeight+$ConcretePadding;
297  # no foam on bottom or top, no concrete on top
298 $DetEncLength = $CryostatLength+2*$TotalPadding;
299 
300 
301  # We want the world origin to be at the very front of the fiducial volume.
302  # move it to the front of the enclosure, then back it up through the concrete/foam,
303  # then through the Cryostat shell, then through the upstream dead LAr (including the
304  # dead LAr on the edge of the TPC, but this is covered in $UpstreamLArPadding).
305  # This is to be added to the z position of every volume in volWorld
306 
307 $OriginZSet = $DetEncLength/2
308  - $TotalPadding
309  - $SteelThickness
310  - $UpstreamLArPadding;
311 
312  # We want the world origin to be vertically centered between the stacked APAs.
313  # the cryostat sits on top of concrete padding, move the detector enclosure back
314  # This is to be added to the y position of every volume in volWorld
315 
316 $OriginYSet = $HeightGaseousAr/2
317  - $ConcretePadding/2;
318 
319  # similar X set variable may be set if needed later
320 
321 
322 # TODO: Needs work from here on... wait until design stabilizes
323 
324 $RockThickness = 3000;
325 
326 $WorldWidth = 3*$RockThickness;
327 $WorldHeight = 3*$RockThickness;
328 $WorldLength = 3*$RockThickness;
329 
330 
331 
332 
333 
334 
335 
336 
337 #############################################################
338 ############## Service Building and Hillside ###############
339 
340 $HighBayTopRockThickness = 320;
341 $LowBayTopRockThickness = 250;
342 
343 $ServiceBuildingExtraWidth = 890;
344  # for now 890 is the difference between 4244 and 3354, which are
345  # the spec'd service building width and detector enclosure width 11/1/2012
346 
347 #LowBayDefinitions
348 $LowBayInsideWidth = $DetEncWidth + $ServiceBuildingExtraWidth;
349 $LowBayInsideHeight = 350;
350 $LowBayInsideLength = 2666;
351 
352 $LowBayCeilingThickness = 2*$ConcretePadding;
353 $LowBaySideWallThickness = $ConcretePadding;
354 $LowBayBackWallThickness = $ConcretePadding;
355 
356 #HighBayDefinitions
357 $HighBayInsideWidth = $DetEncWidth + $ServiceBuildingExtraWidth;
358 $HighBayInsideHeight = 1048;
359 $HighBayInsideLength = 1719;
360 
361 $HighBayCeilingThickness = 2*$ConcretePadding;
362 $HighBayFrontWallThickness = $ConcretePadding;
363 $HighBaySideWallThickness = $ConcretePadding;
364 $HighBayBackWallThickness = 2*$ConcretePadding;
365 
366 $HighBayOverlap = 500; #the opening for access to the DetEncl and position of the HighBay
367 
368 $GabionThickness = 400;
369 
370 $SlopeRockFillHeight = ($HighBayInsideHeight + $HighBayCeilingThickness +$HighBayTopRockThickness - $LowBayInsideHeight - $LowBayCeilingThickness - $LowBayTopRockThickness);
371 $SlopeRockFillLength = $SlopeRockFillHeight*2;
372 
373 $HillSideAngle = 20; #in degrees
374 $HillSideLength = $RockThickness + $DetEncLength;
375 $HillSideHeight = $HillSideLength * tan( deg2rad($HillSideAngle));
376 $HillSideWidth = $RockThickness - ($ServiceBuildingExtraWidth/2 + $LowBaySideWallThickness + $GabionThickness);
377 
378 $HillSideBoxLength = $RockThickness - ($LowBayInsideLength +$LowBayBackWallThickness + $HighBayOverlap +$HighBayBackWallThickness - $DetEncLength );
379 $HillSideBoxWidth = $LowBayInsideWidth + 2*$LowBaySideWallThickness + 2*$GabionThickness;
380 $HillSideBoxHeight = $LowBayInsideHeight +$LowBayCeilingThickness +$LowBayTopRockThickness;
381 
382 $HillSideMiddleLength = $RockThickness + $DetEncLength - $HighBayOverlap - $HighBayBackWallThickness - $SlopeRockFillLength;
383 $HillSideMiddleHeight = $HillSideMiddleLength * tan( deg2rad($HillSideAngle));
384 $HillSideMiddleWidth = $HillSideBoxWidth;
385 
386 
387 
388 
389 
390 #################### Help get your bearings when drawing a symmetric detector. only in world.
391 
392 $PosDirCubeSide = 0;
393 if (defined $helpcube)
394 {
395 $PosDirCubeSide = $ArToAr; #seems to be a good proportion
396 }
397 
398 
399 
400 
401 
402 #+++++++++++++++++++++++++ End defining variables ++++++++++++++++++++++++++
403 
404 # run the sub routines that generate the fragments
405 
406 gen_Define(); # generates definitions at beginning of GDML
407 gen_Materials(); # generates materials to be used
408 
409 gen_TPC(); # generates wires, wire planes, and puts them in volTPC
410  # This is the bulk of the code, and has zero wires option
411 gen_Cryostat(); # places (2*nAPAWide x nAPAHigh x nAPALong) volTPC,
412  # half rotated 180 about Y
413 gen_Enclosure(); # places two cryostats and concrete volumes
414 
415 #gen_ServiceBuilding(); #puts the service building over top of the Enclosure, note that the floor of the service
416  #building is built here but carved out of the surrounding rock
417 
418 #gen_HillSide(); #puts the rock around the detector enclosure and rock fill and hillside around the service building
419 
420 gen_World(); # places the enclosure among DUSEL Rock
421 
422 
423 write_fragments(); # writes the XML input for make_gdml.pl
424  # which zips together the final GDML
425 exit;
426 
427 
428 #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
429 #+++++++++++++++++++++++++++++++++++++++++ usage +++++++++++++++++++++++++++++++++++++++++
430 #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
431 
432 sub usage()
433 {
434  print "Usage: $0 [-h|--help] [-o|--output <fragments-file>] [-s|--suffix <string>]\n";
435  print " if -o is omitted, output goes to STDOUT; <fragments-file> is input to make_gdml.pl\n";
436  print " -s <string> appends the string to the file names; useful for multiple detector versions\n";
437  print " -h prints this message, then quits\n";
438 }
439 
440 
441 
442 #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
443 #++++++++++++++++++++++++++++++++++++++ gen_Define +++++++++++++++++++++++++++++++++++++++
444 #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
445 
446 sub gen_Define()
447 {
448 
449 # Create the <define> fragment file name,
450 # add file to list of fragments,
451 # and open it
452  $DEF = "dune_10kT_Def" . $suffix . ".gdml";
453  push (@gdmlFiles, $DEF);
454  $DEF = ">" . $DEF;
455  open(DEF) or die("Could not open file $DEF for writing");
456 
457 
458 print DEF <<EOF;
459 <?xml version='1.0'?>
460 <gdml>
461 <define>
462  <position name="posTPCActive" unit="cm" x="$posTPCActive_X" y="$posTPCActive_Y" z="$posTPCActive_Z"/>
463  <position name="posCenter" unit="cm" x="0" y="0" z="0"/>
464  <rotation name="rPlus90AboutX" unit="deg" x="90" y="0" z="0"/>
465  <rotation name="rMinus90AboutY" unit="deg" x="0" y="270" z="0"/>
466  <rotation name="rMinus90AboutYMinus90AboutX" unit="deg" x="270" y="270" z="0"/>
467  <rotation name="rPlusUAngleAboutX" unit="deg" x="90-$UAngle" y="0" z="0"/>
468  <rotation name="rPlusVAngleAboutX" unit="deg" x="90+$VAngle" y="0" z="0"/>
469  <rotation name="rPlus180AboutY" unit="deg" x="0" y="180" z="0"/>
470  <rotation name="rIdentity" unit="deg" x="0" y="0" z="0"/>
471 </define>
472 </gdml>
473 EOF
474  close (DEF);
475 }
476 
477 
478 
479 
480 #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
481 #+++++++++++++++++++++++++++++++++++++ gen_Materials +++++++++++++++++++++++++++++++++++++
482 #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
483 
484 sub gen_Materials()
485 {
486 
487 # Create the <materials> fragment file name,
488 # add file to list of output GDML fragments,
489 # and open it
490  $MAT = "dune_10kT_Materials" . $suffix . ".gdml";
491  push (@gdmlFiles, $MAT);
492  $MAT = ">" . $MAT;
493  open(MAT) or die("Could not open file $MAT for writing");
494 
495 
496  print MAT <<EOF;
497 <materials>
498  <element name="videRef" formula="VACUUM" Z="1"> <atom value="1"/> </element>
499  <element name="bromine" formula="Br" Z="35"> <atom value="79.904"/> </element>
500  <element name="hydrogen" formula="H" Z="1"> <atom value="1.0079"/> </element>
501  <element name="nitrogen" formula="N" Z="7"> <atom value="14.0067"/> </element>
502  <element name="oxygen" formula="O" Z="8"> <atom value="15.999"/> </element>
503  <element name="aluminum" formula="Al" Z="13"> <atom value="26.9815"/> </element>
504  <element name="silicon" formula="Si" Z="14"> <atom value="28.0855"/> </element>
505  <element name="carbon" formula="C" Z="6"> <atom value="12.0107"/> </element>
506  <element name="potassium" formula="K" Z="19"> <atom value="39.0983"/> </element>
507  <element name="chromium" formula="Cr" Z="24"> <atom value="51.9961"/> </element>
508  <element name="iron" formula="Fe" Z="26"> <atom value="55.8450"/> </element>
509  <element name="nickel" formula="Ni" Z="28"> <atom value="58.6934"/> </element>
510  <element name="calcium" formula="Ca" Z="20"> <atom value="40.078"/> </element>
511  <element name="magnesium" formula="Mg" Z="12"> <atom value="24.305"/> </element>
512  <element name="sodium" formula="Na" Z="11"> <atom value="22.99"/> </element>
513  <element name="titanium" formula="Ti" Z="22"> <atom value="47.867"/> </element>
514  <element name="argon" formula="Ar" Z="18"> <atom value="39.9480"/> </element>
515  <element name="sulphur" formula="S" Z="16"> <atom value="32.065"/> </element>
516  <element name="phosphorus" formula="P" Z="15"> <atom value="30.973"/> </element>
517 
518  <material name="Vacuum" formula="Vacuum">
519  <D value="1.e-25" unit="g/cm3"/>
520  <fraction n="1.0" ref="videRef"/>
521  </material>
522 
523  <material name="ALUMINUM_Al" formula="ALUMINUM_Al">
524  <D value="2.6990" unit="g/cm3"/>
525  <fraction n="1.0000" ref="aluminum"/>
526  </material>
527 
528  <material name="SILICON_Si" formula="SILICON_Si">
529  <D value="2.3300" unit="g/cm3"/>
530  <fraction n="1.0000" ref="silicon"/>
531  </material>
532 
533  <material name="epoxy_resin" formula="C38H40O6Br4">
534  <D value="1.1250" unit="g/cm3"/>
535  <composite n="38" ref="carbon"/>
536  <composite n="40" ref="hydrogen"/>
537  <composite n="6" ref="oxygen"/>
538  <composite n="4" ref="bromine"/>
539  </material>
540 
541  <material name="SiO2" formula="SiO2">
542  <D value="2.2" unit="g/cm3"/>
543  <composite n="1" ref="silicon"/>
544  <composite n="2" ref="oxygen"/>
545  </material>
546 
547  <material name="Al2O3" formula="Al2O3">
548  <D value="3.97" unit="g/cm3"/>
549  <composite n="2" ref="aluminum"/>
550  <composite n="3" ref="oxygen"/>
551  </material>
552 
553  <material name="Fe2O3" formula="Fe2O3">
554  <D value="5.24" unit="g/cm3"/>
555  <composite n="2" ref="iron"/>
556  <composite n="3" ref="oxygen"/>
557  </material>
558 
559  <material name="CaO" formula="CaO">
560  <D value="3.35" unit="g/cm3"/>
561  <composite n="1" ref="calcium"/>
562  <composite n="1" ref="oxygen"/>
563  </material>
564 
565  <material name="MgO" formula="MgO">
566  <D value="3.58" unit="g/cm3"/>
567  <composite n="1" ref="magnesium"/>
568  <composite n="1" ref="oxygen"/>
569  </material>
570 
571  <material name="Na2O" formula="Na2O">
572  <D value="2.27" unit="g/cm3"/>
573  <composite n="2" ref="sodium"/>
574  <composite n="1" ref="oxygen"/>
575  </material>
576 
577  <material name="TiO2" formula="TiO2">
578  <D value="4.23" unit="g/cm3"/>
579  <composite n="1" ref="titanium"/>
580  <composite n="2" ref="oxygen"/>
581  </material>
582 
583  <material name="FeO" formula="FeO">
584  <D value="5.745" unit="g/cm3"/>
585  <composite n="1" ref="iron"/>
586  <composite n="1" ref="oxygen"/>
587  </material>
588 
589  <material name="CO2" formula="CO2">
590  <D value="1.562" unit="g/cm3"/>
591  <composite n="1" ref="iron"/>
592  <composite n="2" ref="oxygen"/>
593  </material>
594 
595  <material name="P2O5" formula="P2O5">
596  <D value="1.562" unit="g/cm3"/>
597  <composite n="2" ref="phosphorus"/>
598  <composite n="5" ref="oxygen"/>
599  </material>
600 
601  <material formula=" " name="DUSEL_Rock">
602  <D value="2.82" unit="g/cm3"/>
603  <fraction n="0.5267" ref="SiO2"/>
604  <fraction n="0.1174" ref="FeO"/>
605  <fraction n="0.1025" ref="Al2O3"/>
606  <fraction n="0.0473" ref="MgO"/>
607  <fraction n="0.0422" ref="CO2"/>
608  <fraction n="0.0382" ref="CaO"/>
609  <fraction n="0.0240" ref="carbon"/>
610  <fraction n="0.0186" ref="sulphur"/>
611  <fraction n="0.0053" ref="Na2O"/>
612  <fraction n="0.00070" ref="P2O5"/>
613  <fraction n="0.0771" ref="oxygen"/>
614  </material>
615 
616  <material name="fibrous_glass">
617  <D value="2.74351" unit="g/cm3"/>
618  <fraction n="0.600" ref="SiO2"/>
619  <fraction n="0.118" ref="Al2O3"/>
620  <fraction n="0.001" ref="Fe2O3"/>
621  <fraction n="0.224" ref="CaO"/>
622  <fraction n="0.034" ref="MgO"/>
623  <fraction n="0.010" ref="Na2O"/>
624  <fraction n="0.013" ref="TiO2"/>
625  </material>
626 
627  <material name="FR4">
628  <D value="1.98281" unit="g/cm3"/>
629  <fraction n="0.47" ref="epoxy_resin"/>
630  <fraction n="0.53" ref="fibrous_glass"/>
631  </material>
632 
633  <material name="STEEL_STAINLESS_Fe7Cr2Ni" formula="STEEL_STAINLESS_Fe7Cr2Ni">
634  <D value="7.9300" unit="g/cm3"/>
635  <fraction n="0.0010" ref="carbon"/>
636  <fraction n="0.1792" ref="chromium"/>
637  <fraction n="0.7298" ref="iron"/>
638  <fraction n="0.0900" ref="nickel"/>
639  </material>
640 
641  <material name="LAr" formula="LAr">
642  <D value="1.40" unit="g/cm3"/>
643  <fraction n="1.0000" ref="argon"/>
644  </material>
645 
646  <material name="ArGas" formula="ArGas">
647  <D value="0.00166" unit="g/cm3"/>
648  <fraction n="1.0" ref="argon"/>
649  </material>
650 
651  <material formula=" " name="Air">
652  <D value="0.001205" unit="g/cm3"/>
653  <fraction n="0.781154" ref="nitrogen"/>
654  <fraction n="0.209476" ref="oxygen"/>
655  <fraction n="0.00934" ref="argon"/>
656  </material>
657 
658  <material formula=" " name="G10">
659  <D value="1.7" unit="g/cm3"/>
660  <fraction n="0.2805" ref="silicon"/>
661  <fraction n="0.3954" ref="oxygen"/>
662  <fraction n="0.2990" ref="carbon"/>
663  <fraction n="0.0251" ref="hydrogen"/>
664  </material>
665 
666  <material formula=" " name="Granite">
667  <D value="2.7" unit="g/cm3"/>
668  <fraction n="0.438" ref="oxygen"/>
669  <fraction n="0.257" ref="silicon"/>
670  <fraction n="0.222" ref="sodium"/>
671  <fraction n="0.049" ref="aluminum"/>
672  <fraction n="0.019" ref="iron"/>
673  <fraction n="0.015" ref="potassium"/>
674  </material>
675 
676  <material formula=" " name="ShotRock">
677  <D value="2.7*0.6" unit="g/cm3"/>
678  <fraction n="0.438" ref="oxygen"/>
679  <fraction n="0.257" ref="silicon"/>
680  <fraction n="0.222" ref="sodium"/>
681  <fraction n="0.049" ref="aluminum"/>
682  <fraction n="0.019" ref="iron"/>
683  <fraction n="0.015" ref="potassium"/>
684  </material>
685 
686  <material formula=" " name="Dirt">
687  <D value="1.7" unit="g/cm3"/>
688  <fraction n="0.438" ref="oxygen"/>
689  <fraction n="0.257" ref="silicon"/>
690  <fraction n="0.222" ref="sodium"/>
691  <fraction n="0.049" ref="aluminum"/>
692  <fraction n="0.019" ref="iron"/>
693  <fraction n="0.015" ref="potassium"/>
694  </material>
695 
696  <material formula=" " name="Concrete">
697  <D value="2.3" unit="g/cm3"/>
698  <fraction n="0.530" ref="oxygen"/>
699  <fraction n="0.335" ref="silicon"/>
700  <fraction n="0.060" ref="calcium"/>
701  <fraction n="0.015" ref="sodium"/>
702  <fraction n="0.020" ref="iron"/>
703  <fraction n="0.040" ref="aluminum"/>
704  </material>
705 
706  <material formula="H2O" name="Water">
707  <D value="1.0" unit="g/cm3"/>
708  <fraction n="0.1119" ref="hydrogen"/>
709  <fraction n="0.8881" ref="oxygen"/>
710  </material>
711 
712  <material formula="Ti" name="Titanium">
713  <D value="4.506" unit="g/cm3"/>
714  <fraction n="1." ref="titanium"/>
715  </material>
716 
717  <material name="TPB" formula="TPB">
718  <D value="1.40" unit="g/cm3"/>
719  <fraction n="1.0000" ref="argon"/>
720  </material>
721 
722  <material name="Glass">
723  <D value="2.74351" unit="g/cm3"/>
724  <fraction n="0.600" ref="SiO2"/>
725  <fraction n="0.118" ref="Al2O3"/>
726  <fraction n="0.001" ref="Fe2O3"/>
727  <fraction n="0.224" ref="CaO"/>
728  <fraction n="0.034" ref="MgO"/>
729  <fraction n="0.010" ref="Na2O"/>
730  <fraction n="0.013" ref="TiO2"/>
731  </material>
732 
733  <material name="Acrylic">
734  <D value="1.19" unit="g/cm3"/>
735  <fraction n="0.600" ref="carbon"/>
736  <fraction n="0.320" ref="oxygen"/>
737  <fraction n="0.080" ref="hydrogen"/>
738  </material>
739 
740 </materials>
741 EOF
742 
743 close(MAT);
744 }
745 
746 
747 
748 
749 #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
750 #++++++++++++++++++++++++++++++++++++++++ gen_TPC ++++++++++++++++++++++++++++++++++++++++
751 #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
752 
753 
754 sub gen_TPC()
755 {
756 
757 #constructs everything inside volTPC, namely
758 # (moving from left to right, or from +x to -x)
759 # -volTPCPlaneU: with wires angled from vertical slightly different than in V
760 # -volTPCPlaneV: with wires angled from vertical slightly differently than in U
761 # -volTPCPlaneX: with vertical wires
762 
763 
764 # Create the TPC fragment file name,
765 # add file to list of output GDML fragments,
766 # and open it
767  $TPC = "dune_10kT_TPC" . $suffix . ".gdml";
768  push (@gdmlFiles, $TPC);
769  $TPC = ">" . $TPC;
770  open(TPC) or die("Could not open file $TPC for writing");
771 
772 
773 # The standard XML prefix and starting the gdml
774  print TPC <<EOF;
775 <?xml version='1.0'?>
776 <gdml>
777 EOF
778 
779 
780 # All the TPC solids save the wires.
781 print TPC <<EOF;
782 <solids>
783  <box name="TPC" lunit="cm"
784  x="$TPCWidth"
785  y="$TPCHeight"
786  z="$TPCLength"/>
787  <box name="TPCPlane" lunit="cm"
788  x="$TPCWirePlaneThickness"
789  y="$TPCWirePlaneHeight"
790  z="$TPCWirePlaneLength"/>
791  <box name="TPCActive" lunit="cm"
792  x="$TPCActiveWidth"
793  y="$TPCActiveHeight"
794  z="$TPCActiveLength"/>
795 EOF
796 
797 
798 #++++++++++++++++++++++++++++ Wire Solids ++++++++++++++++++++++++++++++
799 
800 print TPC <<EOF;
801 
802  <tube name="TPCWireVert"
803  rmax="0.5*$TPCWireThickness"
804  z="$TPCWirePlaneHeight"
805  deltaphi="360"
806  aunit="deg"
807  lunit="cm"/>
808 EOF
809 
810 # Set number of wires to default to zero, when $wires_on = 0, for a low memory
811 # version. But if $wires_on = 1, calculate the number of wires on each side of each
812 # plane to be used in the for loops
813 
814 my $NumberCornerUWires = 0;
815 my $NumberSideUWires = 0;
816 my $NumberCommonUWires = 0;
817 my $NumberCornerVWires = 0;
818 my $NumberSideVWires = 0;
819 my $NumberCommonVWires = 0;
820 my $NumberVerticalWires = 0;
821 
822 if ($wires_on == 1)
823 {
824  # Number of wires in one corner
825 $NumberCornerUWires = int( $TPCWirePlaneLength/($UWirePitch/$CosUAngle) );
826 
827 $NumberCornerVWires = int( $TPCWirePlaneLength/($VWirePitch/$CosVAngle) );
828 
829  # Total number of wires touching one vertical (longer) side
830  # Note that the total number of wires per plane is this + another set of corner wires
831 $NumberSideUWires = int( $TPCWirePlaneHeight/($UWirePitch/$SinUAngle) );
832 
833 $NumberSideVWires = int( $TPCWirePlaneHeight/($VWirePitch/$SinVAngle) );
834 
835  # Number of wires per side that aren't cut off by the corner
836 $NumberCommonUWires = $NumberSideUWires - $NumberCornerUWires;
837 
838 $NumberCommonVWires = $NumberSideVWires - $NumberCornerVWires;
839 
840  # number of wires on the vertical plane
841 $NumberVerticalWires = int( ($TPCWirePlaneLength-$TPCWireThickness)/$XWirePitch );
842 }
843 
844 # These XML comments throughout make the GDML file easier to navigate
845 print TPC <<EOF;
846 
847 <!--+++++++++++++++++++ U Wire Solids ++++++++++++++++++++++-->
848 
849 EOF
850 
851 # The corner wires for the U plane
852 if ($wires_on==1)
853 {
854  for ($i = 0; $i < $NumberCornerUWires; ++$i)
855  {
856  # Subtraction to avoid corners of wires overlapping the TPCPlane sides,
857  # equal to 0.5*TCPWireThickness*($TanUAngle+1/$TanUAngle),
858  # allowing for 30deg<UAngle
859 
860  print TPC <<EOF;
861  <tube name="TPCWireU$i"
862  rmax="0.5*$TPCWireThickness"
863  z="$UWirePitch*($TanUAngle+1/$TanUAngle)*($i+1)-0.01732"
864  deltaphi="360"
865  aunit="deg"
866  lunit="cm"/>
867 EOF
868 
869  }
870  # Next, the wire used many times in the middle of the U plane.
871  # Subtraction again to avoid wire corners overlapping, equal to
872  # 0.5*TCPWireThickness*2/$TanVAngle, allowing for 30deg<VAngle
873 
874  print TPC <<EOF;
875  <tube name="TPCWireUCommon"
876  rmax="0.5*$TPCWireThickness"
877  z="$TPCWirePlaneLength/$SinUAngle-0.02598"
878  deltaphi="360"
879  aunit="deg"
880  lunit="cm"/>
881 EOF
882 
883 } else {
884 #inform the gdml user why no wires show up when -w=0 is used
885 
886 print TPC <<EOF;
887 
888 
889  <!-- The command -w=0 has been used when running generate_dune_gdml-NEW.pl -->
890 
891  <!-- This GDML version has no wires and uses much less memory -->
892 
893 EOF
894 
895 }
896 
897 print TPC <<EOF;
898 
899 
900 
901 
902 
903 
904 
905 
906 
907 <!--+++++++++++++++++++ V Wire Solids ++++++++++++++++++++++-->
908 
909 
910 EOF
911 
912 # The corner wires for the V plane
913 if ($wires_on==1)
914 {
915  for ($i = 0; $i < $NumberCornerVWires; ++$i)
916  {
917  # Same subtraction to avoid corners of wires overlapping
918  # the TPCPlane sides
919 
920  print TPC <<EOF;
921 
922  <tube name="TPCWireV$i"
923  rmax="0.5*$TPCWireThickness"
924  z="$VWirePitch*($TanVAngle+1/$TanVAngle)*($i+1)-0.01732"
925  deltaphi="360"
926  aunit="deg"
927  lunit="cm"/>
928 
929 EOF
930 
931  }
932 
933  # The wire used many times in the middle of the V plane
934  # Same subtraction as U common
935 
936  print TPC <<EOF;
937  <tube name="TPCWireVCommon"
938  rmax="0.5*$TPCWireThickness"
939  z="$TPCWirePlaneLength/$SinVAngle-0.02598"
940  deltaphi="360"
941  aunit="deg"
942  lunit="cm"/>
943 EOF
944 
945 } else {
946 #inform the gdml user why no wires show up when -w=0 is used
947 
948 print TPC <<EOF;
949 
950 
951  <!-- no wires in this GDML -->
952 
953 EOF
954 
955 }
956 
957 
958 
959 # Begin structure and create the vertical wire logical volume
960 print TPC <<EOF;
961 </solids>
962 <structure>
963  <volume name="volTPCActive">
964  <materialref ref="LAr"/>
965  <solidref ref="TPCActive"/>
966  </volume>
967 
968 
969 
970 
971 
972 
973 
974 
975 
976 <!--+++++++++++++++++ Wire Logical Volumes ++++++++++++++++++++-->
977 
978 EOF
979 
980 
981 if ($wires_on==1)
982 {
983  print TPC <<EOF;
984  <volume name="volTPCWireVert">
985  <materialref ref="STEEL_STAINLESS_Fe7Cr2Ni" />
986  <solidref ref="TPCWireVert" />
987  </volume>
988 EOF
989 
990  # Corner U wires logical volumes
991  for ($i = 0; $i < $NumberCornerUWires; ++$i)
992  {
993  print TPC <<EOF;
994  <volume name="volTPCWireU$i">
995  <materialref ref="STEEL_STAINLESS_Fe7Cr2Ni" />
996  <solidref ref="TPCWireU$i" />
997  </volume>
998 EOF
999 
1000  }
1001 
1002  # Common U wire logical volume, referenced many times
1003  print TPC <<EOF;
1004  <volume name="volTPCWireUCommon">
1005  <materialref ref="STEEL_STAINLESS_Fe7Cr2Ni" />
1006  <solidref ref="TPCWireUCommon" />
1007  </volume>
1008 EOF
1009 
1010  # Corner V wires logical volumes
1011  for ($i = 0; $i < $NumberCornerVWires; ++$i)
1012  {
1013  print TPC <<EOF;
1014  <volume name="volTPCWireV$i">
1015  <materialref ref="STEEL_STAINLESS_Fe7Cr2Ni" />
1016  <solidref ref="TPCWireV$i" />
1017  </volume>
1018 EOF
1019 
1020  }
1021 
1022  # Common V wire logical volume, referenced many times
1023  print TPC <<EOF;
1024  <volume name="volTPCWireVCommon">
1025  <materialref ref="STEEL_STAINLESS_Fe7Cr2Ni" />
1026  <solidref ref="TPCWireVCommon" />
1027  </volume>
1028 EOF
1029 
1030 } else {
1031 #inform the gdml user why no wires show up when -w=0 is used
1032 
1033 print TPC <<EOF;
1034 
1035 
1036  <!-- The command -w=0 has been used when running generate_dune_gdml-NEW.pl -->
1037 
1038  <!-- This GDML version has no wires and uses much less memory -->
1039 
1040 EOF
1041 
1042 }
1043 
1044 
1045 
1046 
1047 #+++++++++++++++++++++++++ Position physical wires ++++++++++++++++++++++++++
1048 
1049 # ++++++++++++++++++++++ U Plane +++++++++++++++++++++++
1050 
1051 # Create U plane logical volume
1052 print TPC <<EOF;
1053 
1054 
1055 
1056 
1057 
1058 
1059 
1060 <!--+++++++++++++++++++++ U Plane ++++++++++++++++++++++++-->
1061 
1062 
1063  <volume name="volTPCPlaneU">
1064  <materialref ref="LAr"/>
1065  <solidref ref="TPCPlane"/>
1066 EOF
1067 
1068 if ($wires_on==0)
1069 {
1070 print TPC <<EOF;
1071 
1072  <!-- no wires -->
1073 
1074 EOF
1075 
1076 } else {
1077 
1078 # Starting with the bottom left corner wires:
1079  # x=0 to center the wires in the plane
1080  # y positioning: (-0.5*$TPCWirePlaneHeight) starts the incremental increase
1081  # from the bottom of the plane, and trigonometry gives the increment
1082  # z positioning: Looking at the plane from the positive x direction,
1083  # (0.5*$TPCWirePlaneLength) starts the incremental increase from
1084  # the lower left corner.
1085  # rotation: same as common wire in code below
1086 
1087 for ($i = 0; $i < $NumberCornerUWires; ++$i)
1088 {
1089 my $ypos = (-0.5*$TPCWirePlaneHeight)+0.5*($i+1)*$UWire_yint;
1090 my $zpos = (0.5*$TPCWirePlaneLength)-0.5*($i+1)*$UWire_zint;
1091 
1092 my $diff=(0.5*$TPCWirePlaneLength)-0.5*($NumberCornerUWires)*$UWire_zint;
1093 my $zpos=$zpos-$diff;
1094 
1095 print TPC <<EOF;
1096  <physvol>
1097  <volumeref ref="volTPCWireU$i"/>
1098  <position name="posTPCWireU$i" unit="cm" x="0" y="$ypos " z="$zpos"/>
1099  <rotationref ref="rPlusUAngleAboutX"/>
1100  </physvol>
1101 EOF
1102 
1103 }
1104 
1105 
1106 # Moving upwards to the common wires:
1107  # x and z are zero to center the wires along a vertical axis
1108  # y positioning: The trick is positioning the lowest common wire so that the pitch
1109  # is consistent, then the increment is double the increment of
1110  # the corner wires since there is no z incriment.
1111  # rotation: wires in \\\\ direction, so +90deg to bring them to vertical and
1112  # +UAngle counterclockwise to arrive at proper orientation
1113 # Note that the counter maintains wire number (in pos. name) counting bottom to top
1114 
1115 for ($i = $NumberCornerUWires; $i < $NumberSideUWires; ++$i)
1116 {
1117 my $ypos = (-0.5*$TPCWirePlaneHeight)+0.5*($NumberCornerUWires)*$UWire_yint+($i+1-$NumberCornerUWires)*$UWire_yint;
1118 
1119 print TPC <<EOF;
1120  <physvol>
1121  <volumeref ref="volTPCWireUCommon"/>
1122  <position name="posTPCWireU$i" unit="cm" x="0" y="$ypos " z="0"/>
1123  <rotationref ref="rPlusUAngleAboutX"/>
1124  </physvol>
1125 EOF
1126 
1127 }
1128 
1129 
1130 # Finally moving to the corner wires on the top right:
1131  # x=0 to center the wires in the plane
1132  # y positioning: plug wire number into same equation
1133  # z positioning: start at z=0 and go negatively at the same z increment
1134  # rotation: same as common wire in code above
1135 # note that the counter maintains wire number shown in the position name
1136 
1137 for ($i = $NumberSideUWires; $i < $NumberSideUWires+$NumberCornerUWires-1; ++$i)
1138 {
1139  # Make a counter to recall the right logical volume reference:
1140  # We want the last U wire in this loop (the highest wire) to be the
1141  # first wire in the logical volume loop for U wires.
1142 
1143 $j = $NumberSideUWires+$NumberCornerUWires - $i - 2;
1144 
1145  # Note that since we are referencing the same logical volumes/same solids for
1146  # the top wires as well as the bottom, the pattern of "stacking" wire on top of wire
1147  # with an incremental separation is likely to cause the top corner wires to be a
1148  # a little shorter than they can be, but never any longer. There is no immediately
1149  # elegant way to fix this, but at 5mm pitch and around 45deg wire orientation, the
1150  # wires can be at most 1cm shorter than possible which is negligible until the top
1151  # 20 wires or so where 1cm is >5% of their length. This also means that there
1152  # could be one more space for a wire left over, but that is highly unlikely.
1153 
1154 my $ypos = (-0.5*$TPCWirePlaneHeight)+0.5*($NumberCornerUWires)*$UWire_yint+($NumberCommonUWires)*$UWire_yint+0.5*($i+1-$NumberSideUWires)*$UWire_yint;
1155 my $zpos = -0.5*($i+1-$NumberSideUWires)*$UWire_zint;
1156 
1157 print TPC <<EOF;
1158  <physvol>
1159  <volumeref ref="volTPCWireU$j"/>
1160  <position name="posTPCWireU$i" unit="cm" x="0" y="$ypos " z="$zpos"/>
1161  <rotationref ref="rPlusUAngleAboutX"/>
1162  </physvol>
1163 EOF
1164 
1165 }
1166 
1167 } #ends else
1168 
1169 
1170 # ++++++++++++++++++++++ V Plane +++++++++++++++++++++++
1171 
1172 # End U plane and create V plane logical volume
1173 print TPC <<EOF;
1174  </volume>
1175 
1176 
1177 
1178 
1179 
1180 
1181 
1182 <!--+++++++++++++++++++++ V Plane ++++++++++++++++++++++++-->
1183 
1184 
1185  <volume name="volTPCPlaneV">
1186  <materialref ref="LAr"/>
1187  <solidref ref="TPCPlane"/>
1188 EOF
1189 
1190 if ($wires_on==0)
1191 {
1192 print TPC <<EOF;
1193 
1194  <!-- no wires -->
1195 
1196 EOF
1197 
1198 } else {
1199 
1200 
1201 # Starting with the bottom right corner wires:
1202  # x=0 to center the wires in the plane
1203  # y positioning: (-0.5*$TPCWirePlaneHeight) starts the incremental increase
1204  # from the bottom of the plane, and trigonometry gives the increment
1205  # z positioning: Looking at the plane from the positive x direction,
1206  # (-0.5*$TPCWirePlaneLength) starts the incremental increase from
1207  # the lower right corner.
1208  # rotation: same as common wire in code below
1209 
1210 for ($i = 0; $i < $NumberCornerVWires; ++$i)
1211 {
1212 my $ypos = (-0.5*$TPCWirePlaneHeight)+0.5*($i+1)*$VWire_yint;
1213 my $zpos = (-0.5*$TPCWirePlaneLength)+0.5*($i+1)*$VWire_zint;
1214 
1215 my $diff=(-0.5*$TPCWirePlaneLength)+0.5*($NumberCornerVWires)*$VWire_zint;
1216 my $zpos=$zpos-$diff;
1217 
1218 print TPC <<EOF;
1219  <physvol>
1220  <volumeref ref="volTPCWireV$i"/>
1221  <position name="posTPCWireV$i" unit="cm" x="0" y="$ypos " z="$zpos"/>
1222  <rotationref ref="rPlusVAngleAboutX"/>
1223  </physvol>
1224 EOF
1225 
1226 }
1227 
1228 
1229 # Moving upwards to the common wires:
1230  # x and z are zero to center the wires along a vertical axis
1231  # y positioning: Plug wire number into the same corner ypos equation
1232  # rotation: wires in //// direction, so +90deg to bring them to vertical and
1233  # --VAngle counterclockwise to arrive at proper orientation
1234 # Note that the counter maintains wire number in the position name
1235 
1236 for ($i = $NumberCornerVWires; $i < $NumberSideVWires; ++$i)
1237 {
1238 my $ypos = (-0.5*$TPCWirePlaneHeight)-0.5*($NumberCornerVWires)*$VWire_yint+($i+1)*$VWire_yint;
1239 
1240 print TPC <<EOF;
1241  <physvol>
1242  <volumeref ref="volTPCWireVCommon"/>
1243  <position name="posTPCWireV$i" unit="cm" x="0" y="$ypos " z="0"/>
1244  <rotationref ref="rPlusVAngleAboutX"/>
1245  </physvol>
1246 EOF
1247 
1248 }
1249 
1250 
1251 # Finally moving to the corner wires on the top right:
1252  # x=0 to center the wires in the y
1253  # plane positioning: plug wire number into same equation
1254  # z positioning: start at z=0 and go positively at the same z increment
1255  # rotation: same as common wire in code above
1256 # note that the counter maintains wire number shown in the position name
1257 
1258 for ($i = $NumberSideVWires; $i < $NumberSideVWires+$NumberCornerVWires-1; ++$i)
1259 {
1260  # Make a counter to recall the right logical volume reference where the last
1261  # wire in this loop is the smallest, first wire in the logical volume loop, just as in U
1262 
1263 $j = $NumberSideVWires+$NumberCornerVWires - $i - 2;
1264 
1265  # Note that since we are referencing the same logical volumes/same solids for
1266  # the top wires as well as the bottom, the pattern of "stacking" wire on top of wire
1267  # with an incremental separation is likely to cause the top corner wires to be a
1268  # a little shorter than they can be, but never any longer. Just as in U
1269 
1270 my $ypos = (-0.5*$TPCWirePlaneHeight)+0.5*($NumberCornerVWires)*$VWire_yint+($NumberCommonVWires)*$VWire_yint+0.5*($i+1-$NumberSideVWires)*$VWire_yint;
1271 my $zpos = 0.5*($i+1-$NumberSideVWires)*$VWire_zint;
1272 
1273 print TPC <<EOF;
1274  <physvol>
1275  <volumeref ref="volTPCWireV$j"/>
1276  <position name="posTPCWireV$i" unit="cm" x="0" y="$ypos " z="$zpos"/>
1277  <rotationref ref="rPlusVAngleAboutX"/>
1278  </physvol>
1279 EOF
1280 }
1281 
1282 } #ends else
1283 
1284 
1285 
1286 # ++++++++++++++++++++++ X Plane +++++++++++++++++++++++
1287 
1288 # End V plane and create X plane logical volume
1289 print TPC <<EOF;
1290  </volume>
1291 
1292 
1293 
1294 
1295 
1296 
1297 <!--+++++++++++++++++++++ X Plane ++++++++++++++++++++++++-->
1298 
1299 
1300  <volume name="volTPCPlaneX">
1301  <materialref ref="LAr"/>
1302  <solidref ref="TPCPlane"/>
1303 EOF
1304 
1305 if ($wires_on==0)
1306 {
1307 print TPC <<EOF;
1308 
1309  <!-- no wires -->
1310 
1311 EOF
1312 
1313 } else {
1314 
1315 # This is the simplest plane, one loop creates all of the wires
1316  # x and y position at zero to center the wires
1317  # z position: moving from front of detector to back, in the positive z direction,
1318  # starting at (-0.5*$TPCWirePlaneLength), the right side looking from
1319  # the +x direction
1320 
1321 for ($i=0; $i<$NumberVerticalWires; ++$i)
1322 {
1323 
1324  # DocDb-6464 pg. 7 says that the center of the first wire is half of the pitch from the edge
1325 my $zpos = (-0.5*$TPCWirePlaneLength)+$TPCWireThickness/2+$XWirePitch*($i+0.5);
1326 
1327 print TPC <<EOF;
1328  <physvol>
1329  <volumeref ref="volTPCWireVert"/>
1330  <position name="posTPCWireX$i" unit="cm" x="0" y="0 " z="$zpos"/>
1331  <rotationref ref="rPlus90AboutX" />
1332  </physvol>
1333 EOF
1334 
1335 }
1336 
1337 } #ends else
1338 
1339 print TPC <<EOF;
1340  </volume>
1341 EOF
1342 
1343 #+++++++++++++++++++++ Position physical wires Above +++++++++++++++++++++
1344 
1345 #wrap up the TPC file
1346 print TPC <<EOF;
1347  <volume name="volTPC">
1348  <materialref ref="LAr" />
1349  <solidref ref="TPC" />
1350  <physvol>
1351  <volumeref ref="volTPCPlaneU" />
1352  <position name="posTPCPlaneU" unit="cm" x="(-$TPCWidth/2)+3*$APAWirePlaneSpacing" y="0" z="0" />
1353  </physvol>
1354  <physvol>
1355  <volumeref ref="volTPCPlaneV" />
1356  <position name="posTPCPlaneV" unit="cm" x="(-$TPCWidth/2)+2*$APAWirePlaneSpacing" y="0" z="0" />
1357  </physvol>
1358  <physvol>
1359  <volumeref ref="volTPCPlaneX" />
1360  <position name="posTPCPlaneX" unit="cm" x="(-$TPCWidth/2)+$APAWirePlaneSpacing" y="0" z="0" />
1361  </physvol>
1362  <physvol>
1363  <volumeref ref="volTPCActive"/>
1364  <positionref ref="posTPCActive"/>
1365  </physvol>
1366  </volume>
1367 </structure>
1368 </gdml>
1369 EOF
1370 
1371  close(GDML);
1372 
1373 } #end of sub gen_TPC
1374 
1375 
1376 
1377 
1378 
1379 #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1380 #++++++++++++++++++++++++++++++++++++++ gen_Cryostat +++++++++++++++++++++++++++++++++++++
1381 #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1382 
1383 sub gen_Cryostat()
1384 {
1385 
1386 # Create the cryostat fragment file name,
1387 # add file to list of output GDML fragments,
1388 # and open it
1389  $CRYO = "dune_10kT_Cryostat" . $suffix . ".gdml";
1390  push (@gdmlFiles, $CRYO);
1391  $CRYO = ">" . $CRYO;
1392  open(CRYO) or die("Could not open file $CRYO for writing");
1393 
1394 
1395 # The standard XML prefix and starting the gdml
1396  print CRYO <<EOF;
1397 <?xml version='1.0'?>
1398 <gdml>
1399 EOF
1400 
1401 
1402 # All the cryostat solids.
1403 print CRYO <<EOF;
1404 <solids>
1405  <box name="Cryostat" lunit="cm"
1406  x="$CryostatWidth"
1407  y="$CryostatHeight"
1408  z="$CryostatLength"/>
1409  <box name="ArgonInterior" lunit="cm"
1410  x="$ArgonWidth"
1411  y="$ArgonHeight"
1412  z="$ArgonLength"/>
1413  <box name="GaseousArgon" lunit="cm"
1414  x="$ArgonWidth"
1415  y="$HeightGaseousAr"
1416  z="$ArgonLength"/>
1417  <subtraction name="SteelShell">
1418  <first ref="Cryostat"/>
1419  <second ref="ArgonInterior"/>
1420  </subtraction>
1421 
1422  <box name="Cathode" lunit="cm"
1423  x="$CPAThickness"
1424  y="$TPCHeight"
1425  z="$TPCLength"/>
1426 
1427  <box name="LightPaddle" lunit="cm"
1428  x="$LightPaddleWidth"
1429  y="$LightPaddleHeight"
1430  z="$LightPaddleLength + $SiPM_z"/>
1431 
1432  <box name="APAFrameYSideHollow" lunit="cm"
1433  x="$APAFrameYSide_x-2*$EdgeFrameSteelThickness"
1434  y="$APAFrameYSide_y-2*$EdgeFrameSteelThickness"
1435  z="$APAFrameYSide_z" />
1436  <box name="APAFrameYSideShell" lunit="cm"
1437  x="$APAFrameYSide_x"
1438  y="$APAFrameYSide_y"
1439  z="$APAFrameYSide_z" />
1440  <subtraction name="APAFrameYSide">
1441  <first ref="APAFrameYSideShell" />
1442  <second ref="APAFrameYSideHollow"/>
1443  <positionref ref="posCenter" />
1444  <rotationref ref="rIdentity" />
1445  </subtraction>
1446 
1447  <box name="APAFrameZSideHollow" lunit="cm"
1448  x="$APAFrameZSide_x-2*$EdgeFrameSteelThickness"
1449  y="$APAFrameZSide_y-2*$EdgeFrameSteelThickness"
1450  z="$APAFrameZSide_z" />
1451  <box name="APAFrameZSideShell" lunit="cm"
1452  x="$APAFrameZSide_x"
1453  y="$APAFrameZSide_y"
1454  z="$APAFrameZSide_z" />
1455  <subtraction name="APAFrameZSide">
1456  <first ref="APAFrameZSideShell" />
1457  <second ref="APAFrameZSideHollow"/>
1458  <positionref ref="posCenter" />
1459  <rotationref ref="rIdentity" />
1460  </subtraction>
1461 
1462  <box name="APAFrameYOuterSupport" lunit="cm"
1463  x="$EdgeFrameSteelThickness"
1464  y="$APAFrameYOuterSupport_y"
1465  z="$APAFrameYOuterSupport_z" />
1466 
1467  <box name="APAFrameYInnerSupport" lunit="cm"
1468  x="$APAFrameYInnerSupport_x"
1469  y="$APAFrameYInnerSupport_y"
1470  z="$APAFrameYInnerSupport_z" />
1471 
1472  <box name="APAFrameZHalfSupport" lunit="cm"
1473  x="$APAFrameZHalfSupport_x"
1474  y="$APAFrameZHalfSupport_y"
1475  z="$APAFrameZHalfSupport_z" />
1476 </solids>
1477 EOF
1478 
1479 # Cryostat structure
1480 print CRYO <<EOF;
1481 <structure>
1482  <volume name="volSteelShell">
1483  <materialref ref="STEEL_STAINLESS_Fe7Cr2Ni" />
1484  <solidref ref="SteelShell" />
1485  </volume>
1486  <volume name="volGaseousArgon">
1487  <materialref ref="ArGas"/>
1488  <solidref ref="GaseousArgon"/>
1489  </volume>
1490 
1491  <volume name="volCathode">
1492  <materialref ref="STEEL_STAINLESS_Fe7Cr2Ni" />
1493  <solidref ref="Cathode" />
1494  </volume>
1495 
1496  <volume name="volOpDetSensitive_Bar">
1497  <materialref ref="LAr"/>
1498  <solidref ref="LightPaddle"/>
1499  </volume>
1500 
1501  <volume name="volAPAFrameYSide">
1502  <materialref ref="STEEL_STAINLESS_Fe7Cr2Ni" />
1503  <solidref ref="APAFrameYSide" />
1504  </volume>
1505 
1506  <volume name="volAPAFrameZSide">
1507  <materialref ref="STEEL_STAINLESS_Fe7Cr2Ni" />
1508  <solidref ref="APAFrameZSide" />
1509  </volume>
1510 
1511  <volume name="volAPAFrameYOuterSupport">
1512  <materialref ref="STEEL_STAINLESS_Fe7Cr2Ni" />
1513  <solidref ref="APAFrameYOuterSupport" />
1514  </volume>
1515 
1516  <volume name="volAPAFrameYInnerSupport">
1517  <materialref ref="STEEL_STAINLESS_Fe7Cr2Ni" />
1518  <solidref ref="APAFrameYInnerSupport" />
1519  </volume>
1520 
1521  <volume name="volAPAFrameZHalfSupport">
1522  <materialref ref="STEEL_STAINLESS_Fe7Cr2Ni" />
1523  <solidref ref="APAFrameZHalfSupport" />
1524  </volume>
1525 
1526  <volume name="volCryostat">
1527  <materialref ref="LAr" />
1528  <solidref ref="Cryostat" />
1529  <physvol>
1530  <volumeref ref="volGaseousArgon"/>
1531  <position name="posGaseousArgon" unit="cm" x="0" y="$ArgonHeight/2-$HeightGaseousAr/2" z="0"/>
1532  </physvol>
1533  <physvol>
1534  <volumeref ref="volSteelShell"/>
1535  <position name="posSteelShell" unit="cm" x="0" y="0" z="0"/>
1536  </physvol>
1537 EOF
1538 
1539 # nested for loops to place the non-rotated AND rotated volTPC
1540  # x loop rotation: There are six drift volumes. Looking into the
1541  # detector from incident direction, and counting from right (-x) to
1542  # left (+x), odd number volumes need to be rotated 180deg about Y in
1543  # order for the cathode to be on the right of the APA.
1544 
1545 if ($tpc_on==1) {
1546 
1547 
1548  for($k=0 ; $k<$nAPALong ; $k++)
1549  {
1550 
1551  for($i=0 ; $i<$nAPAWide+1 ; $i++)
1552  {
1553 
1554  for($j=0 ; $j<$nAPAHigh ; $j++)
1555  {
1556 
1557 
1558  $CAT_X = - $FiducialWidth/2
1559  + $CPAThickness/2
1560  + $i*($APAWidth);
1561 
1562  $APACenter_x = - $FiducialWidth/2
1563  + $CPAThickness/2
1564  + ($i+0.5)*($APAWidth);
1565 
1566  $APACenter_y = - $ArgonHeight/2 + $SpaceAPAToFloor
1567  + $APAHeight/2
1568  + $j*($APAHeight+$APAVerticalGap);
1569 
1570  $APACenter_z = - $ArgonLength/2
1571  + $UpstreamLArPadding
1572  + $APALength/2
1573  + $k*($APALength+$APALongGap);
1574 
1575 
1576  print CRYO <<EOF;
1577 
1578  <physvol>
1579  <volumeref ref="volCathode" />
1580  <position name="posCathode\-$i\-$j\-$k" unit="cm" x="$CAT_X" y="$APACenter_y" z="$APACenter_z" />
1581  </physvol>
1582 
1583 EOF
1584 
1585  # This if statement is to stop placement of this last set of TPCs/LightPaddles when placing
1586  # the last set of CPAs that do not have any drift volume beyond them in +x
1587 
1588  if ($i<$nAPAWide){
1589 
1590  # place APA volumes around this center: Frame, TPCs, paddles
1591  make_APA($APACenter_x, $APACenter_y, $APACenter_z, $i, $j, $k);
1592 
1593  for ($paddle = 0; $paddle<$nLightPaddlesPerAPA; $paddle++)
1594  {
1595 
1596  # All Light Paddle centers will have the same
1597  # X coordinate as the center of the current APA
1598  # Z coordinate as the current TPC pair
1599  # The Y coordinate must be looped over:
1600 
1601  #the multiplication by j here is a temporary dirty trick to get around some strange behavior
1602 
1603  $Paddle_Y = $APACenter_y
1604  - $APAHeight/2
1605  + $j*$FrameToPaddleSpace
1606  + (1-$j)*($LightPaddleHeight/2 + $APAFrameZSide_y)
1607  + $PaddleYInterval*$paddle;
1608 
1609  print CRYO <<EOF;
1610 
1611  <physvol>
1612  <volumeref ref="volOpDetSensitive"/>
1613  <position name="posPaddle\-$paddle\-TPC\-$i\-$j\-$k" unit="cm" x="$APACenter_x" y="$Paddle_Y" z="$APACenter_z + $SiPM_z/2"/>
1614  <rotationref ref="rIdentity"/>
1615  </physvol>
1616 
1617 EOF
1618 
1619  }#end Paddle for-loop
1620 
1621  }#end if
1622 
1623  } #high
1624  } #wide
1625  } #long
1626 
1627 }# end if tpc
1628 
1629 print CRYO <<EOF;
1630  </volume>
1631 </structure>
1632 </gdml>
1633 EOF
1634 
1635 close(CRYO);
1636 }
1637 
1638 
1639 #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1640 #+++++++++++++++++++++++++++++++++++++++ make_APA ++++++++++++++++++++++++++++++++++++++++
1641 #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1642 
1643 # Must be called only within gen_Cryostat(),
1644 # with APA center x, y, and z as three arguments
1645 
1646 sub make_APA()
1647 {
1648 
1649  $TPC_X_rot = $_[0] - ($TPCWidth+$APAFrame_x)/2;
1650  $TPC_X = $_[0] + ($TPCWidth+$APAFrame_x)/2;
1651  #for the rotation, remember the wires are on the negative side of volTPC
1652 
1653  print CRYO <<EOF;
1654  <physvol>
1655  <volumeref ref="volAPAFrameYOuterSupport"/>
1656  <position name="posAPAFrameYOuterSupportNeg\-$_[3]\-$_[4]\-$_[5]" unit="cm"
1657  x="$_[0] - ($APAFrameYOuterSupport_x + $APAFrameYInnerSupport_x/2 - $EdgeFrameSteelThickness/2)"
1658  y="$_[1]"
1659  z="$_[2]"/>
1660  <rotationref ref="rIdentity"/>
1661  </physvol>
1662  <physvol>
1663  <volumeref ref="volAPAFrameYOuterSupport"/>
1664  <position name="posAPAFrameYOuterSupportPos\-$_[3]\-$_[4]\-$_[5]" unit="cm"
1665  x="$_[0] + ($APAFrameYOuterSupport_x + $APAFrameYInnerSupport_x/2 - $EdgeFrameSteelThickness/2)"
1666  y="$_[1]"
1667  z="$_[2]"/>
1668  <rotationref ref="rIdentity"/>
1669  </physvol>
1670  <physvol>
1671  <volumeref ref="volAPAFrameYSide"/>
1672  <position name="posAPAFrameYSideNeg\-$_[3]\-$_[4]\-$_[5]" unit="cm"
1673  x="$_[0]"
1674  y="$_[1]"
1675  z="$_[2] - $APALength/2 + $APAFrameYSide_z/2"/>
1676  <rotationref ref="rIdentity"/>
1677  </physvol>
1678  <physvol>
1679  <volumeref ref="volAPAFrameYSide"/>
1680  <position name="posAPAFrameYSidePos\-$_[3]\-$_[4]\-$_[5]" unit="cm"
1681  x="$_[0]"
1682  y="$_[1]"
1683  z="$_[2] + $APALength/2 - $APAFrameYSide_z/2"/>
1684  <rotationref ref="rIdentity"/>
1685  </physvol>
1686  <physvol>
1687  <volumeref ref="volAPAFrameZSide"/>
1688  <position name="posAPAFrameZSidePos\-$_[3]\-$_[4]\-$_[5]" unit="cm"
1689  x="$_[0]"
1690  y="$_[1] + $APAHeight/2 - $APAFrameZSide_y/2"
1691  z="$_[2]"/>
1692  <rotationref ref="rIdentity"/>
1693  </physvol>
1694  <physvol>
1695  <volumeref ref="volAPAFrameZSide"/>
1696  <position name="posAPAFrameZSideNeg\-$_[3]\-$_[4]\-$_[5]" unit="cm"
1697  x="$_[0]"
1698  y="$_[1] - $APAHeight/2 + $APAFrameZSide_y/2"
1699  z="$_[2]"/>
1700  <rotationref ref="rIdentity"/>
1701  </physvol>
1702 
1703  <physvol>
1704  <volumeref ref="volTPC"/>
1705  <position name="posTPCrot\-$_[3]\-$_[4]\-$_[5]" unit="cm" x="$TPC_X_rot" y="$_[1]" z="$_[2]"/>
1706  <rotationref ref="rPlus180AboutY"/>
1707  </physvol>
1708  <physvol>
1709  <volumeref ref="volTPC"/>
1710  <position name="posTPC\-$_[3]\-$_[4]\-$_[5]" unit="cm" x="$TPC_X" y="$_[1]" z="$_[2]"/>
1711  <rotationref ref="rIdentity"/>
1712  </physvol>
1713 EOF
1714 
1715 }
1716 
1717 
1718 
1719 #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1720 #+++++++++++++++++++++++++++++++++++++ gen_Enclosure +++++++++++++++++++++++++++++++++++++
1721 #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1722 
1723 sub gen_Enclosure()
1724 {
1725 
1726 # Create the detector enclosure fragment file name,
1727 # add file to list of output GDML fragments,
1728 # and open it
1729  $ENCL = "dune_10kT_DetEnclosure" . $suffix . ".gdml";
1730  push (@gdmlFiles, $ENCL);
1731  $ENCL = ">" . $ENCL;
1732  open(ENCL) or die("Could not open file $ENCL for writing");
1733 
1734 
1735 # The standard XML prefix and starting the gdml
1736  print ENCL <<EOF;
1737 <?xml version='1.0'?>
1738 <gdml>
1739 EOF
1740 
1741 
1742 # All the detector enclosure solids.
1743 print ENCL <<EOF;
1744 <solids>
1745 
1746  <box name="FoamPadBlock" lunit="cm"
1747  x="$CryostatWidth+2*$FoamPadding"
1748  y="$CryostatHeight"
1749  z="$CryostatLength+2*$FoamPadding" />
1750 
1751  <subtraction name="FoamPadding">
1752  <first ref="FoamPadBlock"/>
1753  <second ref="Cryostat"/>
1754  <position name="posCryoInFoam" x="0" y="0" z="0"/>
1755  </subtraction>
1756 
1757  <box name="ConcreteWall" lunit="cm"
1758  x="$ArToAr - 2*$FoamPadding"
1759  y="$CryostatHeight+$ConcretePadding"
1760  z="$CryostatLength+2*$TotalPadding"/>
1761 
1762  <box name="DetEnclosure" lunit="cm"
1763  x="$DetEncWidth"
1764  y="$DetEncHeight"
1765  z="$DetEncLength"/>
1766 
1767 </solids>
1768 EOF
1769 
1770 
1771 
1772 # Detector enclosure structure
1773  print ENCL <<EOF;
1774 <structure>
1775  <volume name="volFoamPadding">
1776  <materialref ref="fibrous_glass"/>
1777  <solidref ref="FoamPadding"/>
1778  </volume>
1779 
1780  <volume name="volConcreteWall">
1781  <materialref ref="Concrete"/>
1782  <solidref ref="ConcreteWall"/>
1783  </volume>
1784 
1785  <volume name="volDetEnclosure">
1786  <materialref ref="Concrete"/>
1787  <solidref ref="DetEnclosure"/>
1788 EOF
1789 
1790  # option for one cryostat
1791 if($nCryos==1){
1792  print ENCL <<EOF;
1793  <physvol>
1794  <volumeref ref="volFoamPadding"/>
1795  <position name="posNegFoamCryo" unit="cm" x="-$CryostatWidth/2-$ArToAr/2" y="$ConcretePadding/2" z="0"/>
1796  </physvol>
1797  <physvol>
1798  <volumeref ref="volCryostat"/>
1799  <position name="posNegCryo" unit="cm" x="-$CryostatWidth/2-$ArToAr/2" y="$ConcretePadding/2" z="0" />
1800  </physvol>
1801 EOF
1802  } elsif ($nCryos==2) {
1803  print ENCL <<EOF;
1804  <physvol>
1805  <volumeref ref="volFoamPadding"/>
1806  <position name="posNegFoamCryo" unit="cm" x="-$CryostatWidth/2-$ArToAr/2" y="$ConcretePadding/2" z="0"/>
1807  </physvol>
1808  <physvol>
1809  <volumeref ref="volCryostat"/>
1810  <position name="posNegCryo" unit="cm" x="-$CryostatWidth/2-$ArToAr/2" y="$ConcretePadding/2" z="0" />
1811  </physvol>
1812  <physvol>
1813  <volumeref ref="volConcreteWall"/>
1814  <position name="posConcreteWall" unit="cm" x="0" y="0" z="0"/>
1815  </physvol>
1816  <physvol>
1817  <volumeref ref="volFoamPadding"/>
1818  <position name="posPosFoamCryo" unit="cm" x="$CryostatWidth/2+$ArToAr/2" y="$ConcretePadding/2" z="0"/>
1819  </physvol>
1820  <physvol>
1821  <volumeref ref="volCryostat"/>
1822  <position name="posPosCryo" unit="cm" x="$CryostatWidth/2+$ArToAr/2" y="$ConcretePadding/2" z="0" />
1823  </physvol>
1824 EOF
1825  }
1826 
1827 
1828  print ENCL <<EOF;
1829  </volume>
1830 EOF
1831 
1832 print ENCL <<EOF;
1833 </structure>
1834 </gdml>
1835 EOF
1836 
1837 close(ENCL);
1838 }
1839 
1840 
1841 
1842 
1843 
1844 #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1845 #+++++++++++++++++++++++++++++++++++++ gen_ServiceBuilding +++++++++++++++++++++++++++++++++++++
1846 #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1847 
1848 sub gen_ServiceBuilding()
1849 {
1850 
1851 # Create the service building above the detector enclosure
1852 # add file to list of output GDML fragments,
1853 # and open it
1854  $SRVBUILD = "dune_10kT_DetServiceBuilding" . $suffix . ".gdml";
1855  push (@gdmlFiles, $SRVBUILD);
1856  $SRVBUILD = ">" . $SRVBUILD;
1857  open(SRVBUILD) or die("Could not open file $SRVBUILD for writing");
1858 
1859 
1860 # The standard XML prefix and starting the gdml
1861  print SRVBUILD <<EOF;
1862 <?xml version='1.0'?>
1863 <gdml>
1864 EOF
1865 
1866 
1867 # All the detector enclosure solids.
1868 print SRVBUILD <<EOF;
1869 <solids>
1870 
1871  <box name="LowBaySolid" lunit="cm"
1872  x="$LowBayInsideWidth + 2*$LowBaySideWallThickness"
1873  y="$LowBayInsideHeight + $LowBayCeilingThickness"
1874  z="$LowBayInsideLength + $LowBayBackWallThickness"/>
1875 
1876  <box name="LowBayInsideSpace" lunit="cm"
1877  x="$LowBayInsideWidth"
1878  y="$LowBayInsideHeight"
1879  z="$LowBayInsideLength" />
1880 
1881  <subtraction name="LowBayWalls">
1882  <first ref="LowBaySolid"/>
1883  <second ref="LowBayInsideSpace"/>
1884  <position name="posLowBayInsideSpace" unit="cm" x="0" y="-$LowBayCeilingThickness/2-0.001" z="$LowBayBackWallThickness/2+0.01"/>
1885  </subtraction>
1886 
1887  <box name="HighBaySolid" lunit="cm"
1888  x="$HighBayInsideWidth + 2*$HighBaySideWallThickness"
1889  y="$HighBayInsideHeight + $HighBayCeilingThickness"
1890  z="$HighBayInsideLength + $HighBayBackWallThickness + $HighBayFrontWallThickness"/>
1891 
1892  <box name="HighBayInsideSpace" lunit="cm"
1893  x="$HighBayInsideWidth"
1894  y="$HighBayInsideHeight"
1895  z="$HighBayInsideLength" />
1896 
1897  <box name="HighBayRearOpening" lunit="cm" x="$DetEncWidth"
1898  y="$LowBayInsideHeight" z="$HighBayBackWallThickness+0.01"/>
1899 
1900  <subtraction name="HighBayTemp">
1901  <first ref="HighBaySolid"/>
1902  <second ref="HighBayInsideSpace"/>
1903  <position name="posHighBayInsideSpace" unit="cm" x="0" y="-$HighBayCeilingThickness/2-0.001" z="$HighBayBackWallThickness/4"/>
1904  </subtraction>
1905 
1906  <subtraction name="HighBayWalls">
1907  <first ref="HighBayTemp"/>
1908  <second ref="HighBayRearOpening"/>
1909  <position name="posHighBayRearOpening" unit="cm" x="0" y="-(($HighBayInsideHeight+$HighBayCeilingThickness)/2 - $LowBayInsideHeight/2)" z="-($HighBayBackWallThickness/4+$HighBayInsideLength/2)"/>
1910  </subtraction>
1911 
1912  <box name="FrontGabion" lunit="cm"
1913  x="$HighBayInsideWidth + 2*$HighBaySideWallThickness + 2*$GabionThickness"
1914  y="$HighBayInsideHeight + $HighBayCeilingThickness +$HighBayTopRockThickness"
1915  z="$GabionThickness"/>
1916 
1917  <box name="SideGabion" lunit="cm"
1918  x="$GabionThickness"
1919  y="$HighBayInsideHeight + $HighBayCeilingThickness +$HighBayTopRockThickness"
1920  z="$HighBayInsideLength + $HighBayFrontWallThickness +$HighBayBackWallThickness"/>
1921 
1922  <box name="RearGabion" lunit="cm"
1923  x="$GabionThickness"
1924  y="$LowBayInsideHeight + $LowBayCeilingThickness +$LowBayTopRockThickness"
1925  z="$LowBayInsideLength +$LowBayBackWallThickness"/>
1926 
1927  <box name="HighBayTopRock" lunit="cm"
1928  x="$HighBayInsideWidth+2*$HighBaySideWallThickness"
1929  y="$HighBayTopRockThickness"
1930  z="$HighBayInsideLength + $HighBayFrontWallThickness + $HighBayBackWallThickness"/>
1931 
1932  <box name="LowBayTopRock" lunit="cm"
1933  x="$LowBayInsideWidth+2*$LowBaySideWallThickness"
1934  y="$LowBayTopRockThickness"
1935  z="$LowBayInsideLength + $LowBayBackWallThickness"/>
1936 
1937  <arb8 name="SlopeRockFill" v1x="0" v1y="0" v2x="$SlopeRockFillLength" v2y="0" v3x="0.01" v3y="$SlopeRockFillHeight" v4x="0" v4y="$SlopeRockFillHeight" v5x="0" v5y="0" v6x="$SlopeRockFillLength" v6y="0" v7x="0.01" v7y="$SlopeRockFillHeight" v8x="0" v8y="$SlopeRockFillHeight" dz="$LowBayInsideWidth/2+$LowBaySideWallThickness + $GabionThickness" lunit="cm"/>
1938 
1939  <arb8 name="HillSide" v1x="0" v1y="0" v2x="$HillSideHeight" v2y="0" v3x="0.001" v3y="$HillSideLength" v4x="0" v4y="$HillSideLength" v5x="0" v5y="0" v6x="$HillSideHeight" v6y="0" v7x="0.001" v7y="$HillSideLength" v8x="0" v8y="$HillSideLength" dz="$HillSideWidth/2" lunit="cm"/>
1940 
1941  <arb8 name="HillSideMiddle" v1x="0" v1y="0" v2x="$HillSideMiddleHeight" v2y="0" v3x="0.001" v3y="$HillSideMiddleLength" v4x="0" v4y="$HillSideMiddleLength" v5x="0" v5y="0" v6x="$HillSideMiddleHeight" v6y="0" v7x="0.001" v7y="$HillSideMiddleLength" v8x="0" v8y="$HillSideMiddleLength" dz="$HillSideMiddleWidth/2" lunit="cm"/>
1942 
1943  <box name="HillSideBox" lunit="cm"
1944  x="$HillSideBoxWidth"
1945  y="$HillSideBoxHeight"
1946  z="$HillSideBoxLength"/>
1947 
1948 </solids>
1949 EOF
1950 
1951 
1952 # ServiceBuilding structure
1953  print SRVBUILD <<EOF;
1954 <structure>
1955  <volume name="volHighBay">
1956  <materialref ref="Concrete"/>
1957  <solidref ref="HighBayWalls"/>
1958  </volume>
1959  <volume name="volLowBay">
1960  <materialref ref="Concrete"/>
1961  <solidref ref="LowBayWalls"/>
1962  </volume>
1963  <volume name="volFrontGabion">
1964  <materialref ref="DUSEL_Rock"/>
1965  <solidref ref="FrontGabion"/>
1966  </volume>
1967  <volume name="volSideGabion1">
1968  <materialref ref="DUSEL_Rock"/>
1969  <solidref ref="SideGabion"/>
1970  </volume>
1971  <volume name="volSideGabion2">
1972  <materialref ref="DUSEL_Rock"/>
1973  <solidref ref="SideGabion"/>
1974  </volume>
1975 
1976  <volume name="volRearGabion1">
1977  <materialref ref="DUSEL_Rock"/>
1978  <solidref ref="RearGabion"/>
1979  </volume>
1980  <volume name="volRearGabion2">
1981  <materialref ref="DUSEL_Rock"/>
1982  <solidref ref="RearGabion"/>
1983  </volume>
1984 
1985  <volume name="volHighBayTopRock">
1986  <materialref ref="DUSEL_Rock"/>
1987  <solidref ref="HighBayTopRock"/>
1988  </volume>
1989 
1990  <volume name="volLowBayTopRock">
1991  <materialref ref="DUSEL_Rock"/>
1992  <solidref ref="LowBayTopRock"/>
1993  </volume>
1994 
1995  <volume name="volSlopeRockFill">
1996  <materialref ref="DUSEL_Rock"/>
1997  <solidref ref="SlopeRockFill"/>
1998  </volume>
1999 
2000  <volume name="volHillSideBox">
2001  <materialref ref="DUSEL_Rock"/>
2002  <solidref ref="HillSideBox"/>
2003  </volume>
2004 
2005  <volume name="volHillSideMiddle">
2006  <materialref ref="DUSEL_Rock"/>
2007  <solidref ref="HillSideMiddle"/>
2008  </volume>
2009 
2010  <volume name="volHillSide1">
2011  <materialref ref="DUSEL_Rock"/>
2012  <solidref ref="HillSide"/>
2013  </volume>
2014 
2015  <volume name="volHillSide2">
2016  <materialref ref="DUSEL_Rock"/>
2017  <solidref ref="HillSide"/>
2018  </volume>
2019 
2020 EOF
2021 
2022 print SRVBUILD <<EOF;
2023 </structure>
2024 </gdml>
2025 EOF
2026 
2027 close(SRVBUILD);
2028 }
2029 
2030 #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2031 #+++++++++++++++++++++++++++++++++++++++ gen_HillSide +++++++++++++++++++++++++++++++++++++++
2032 #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2033 
2034 sub gen_HillSide()
2035 {
2036 # Create the WORLD fragment file name,
2037 # add file to list of output GDML fragments,
2038 # and open it
2039  $HILLSIDE = "dune_10kT_HillSide" . $suffix . ".gdml";
2040  push (@gdmlFiles, $HILLSIDE);
2041  $HILLSIDE = ">" . $HILLSIDE;
2042  open(HILLSIDE) or die("Could not open file $HILLSIDE for writing");
2043 
2044 
2045 # The standard XML prefix and starting the gdml
2046  print HILLSIDE <<EOF;
2047 <?xml version='1.0'?>
2048 <gdml>
2049 EOF
2050 
2051 
2052 # All the World solids.
2053 print HILLSIDE <<EOF;
2054 <solids>
2055 
2056  <box name="Ground" lunit="cm"
2057  x="$DetEncWidth+2*$RockThickness"
2058  y="$DetEncHeight+$RockThickness"
2059  z="$DetEncLength+2*$RockThickness"/>
2060  <box name="GroundRock" lunit="cm"
2061  x="$DetEncWidth+2*$RockThickness"
2062  y="$DetEncHeight+$RockThickness"
2063  z="$DetEncLength+2*$RockThickness"/>
2064 
2065  <box name="DetectorCavern" lunit="cm"
2066  x="$DetEncWidth"
2067  y="$DetEncHeight+0.001"
2068  z="$DetEncLength"/>
2069 
2070  <box name="HighBayFloorSpace" lunit="cm"
2071  x="$HighBayInsideWidth + $HighBaySideWallThickness"
2072  y="$ConcretePadding"
2073  z="$HighBayInsideLength - $HighBayOverlap"/>
2074 
2075  <subtraction name="GroundRockTemp">
2076  <first ref="GroundRock"/>
2077  <second ref="HighBayFloorSpace"/>
2078  <position name="posCavernInGround" unit="cm" x="0" y="$RockThickness/2" z="0"/>
2079  </subtraction>
2080 
2081  <subtraction name="GroundRockWithCavern">
2082  <first ref="GroundRockTemp"/>
2083  <second ref="DetectorCavern"/>
2084  <position name="posCavernInGround" unit="cm" x="0" y="$RockThickness/2" z="0"/>
2085  </subtraction>
2086 
2087 </solids>
2088 EOF
2089 
2090 # World structure
2091 print HILLSIDE <<EOF;
2092 <structure>
2093 
2094  <volume name="volGroundRockWithCavern">
2095  <materialref ref="DUSEL_Rock"/>
2096  <solidref ref="GroundRockWithCavern"/>
2097  </volume>
2098 
2099  <volume name="volGround">
2100  <materialref ref="Air"/>
2101  <solidref ref="Ground"/>
2102  <physvol>
2103  <volumeref ref="volGroundRockWithCavern"/>
2104  <position name="posGroundRockWithCavern" unit="cm" x="0" y="-($RockThickness+$DetEncHeight)/2 + $DetEncHeight/2" z="0"/>
2105  </physvol>
2106  <physvol>
2107  <volumeref ref="volFrontGabion"/>
2108  <position name="posFrontGabion" unit="cm" x="0" y="$DetEncHeight/2 +($HighBayInsideHeight+$HighBayCeilingThickness)/2 + $HighBayTopRockThickness/2" z="($DetEncLength/2+($HighBayInsideLength+$HighBayFrontWallThickness) -$HighBayOverlap) +$GabionThickness/2"/>
2109  </physvol>
2110  <physvol>
2111  <volumeref ref="volSideGabion1"/>
2112  <position name="posSideGabion" unit="cm" x="-($DetEncWidth+$ServiceBuildingExtraWidth+2*$HighBaySideWallThickness+$GabionThickness)/2" y="$DetEncHeight/2 +($HighBayInsideHeight+$HighBayCeilingThickness)/2 + $HighBayTopRockThickness/2" z="($DetEncLength/2+($HighBayInsideLength+$HighBayBackWallThickness+$HighBayFrontWallThickness)/2 -$HighBayBackWallThickness - $HighBayOverlap)"/>
2113  </physvol>
2114  <physvol>
2115  <volumeref ref="volSideGabion2"/>
2116  <position name="posSideGabion" unit="cm" x="($DetEncWidth+$ServiceBuildingExtraWidth+2*$HighBaySideWallThickness+$GabionThickness)/2" y="$DetEncHeight/2 +($HighBayInsideHeight+$HighBayCeilingThickness)/2 + $HighBayTopRockThickness/2" z="($DetEncLength/2+($HighBayInsideLength+$HighBayBackWallThickness+$HighBayFrontWallThickness)/2 -$HighBayBackWallThickness - $HighBayOverlap)"/>
2117  </physvol>
2118  <physvol>
2119  <volumeref ref="volRearGabion1"/>
2120  <position name="posRearGabion" unit="cm" x="-($DetEncWidth+$ServiceBuildingExtraWidth+2*$LowBaySideWallThickness+$GabionThickness)/2" y="$DetEncHeight/2 +($LowBayInsideHeight+$LowBayCeilingThickness)/2 + $LowBayTopRockThickness/2" z="-(($LowBayInsideLength+$LowBayBackWallThickness)/2 + $HighBayOverlap+$HighBayBackWallThickness - $DetEncLength/2)"/>
2121  </physvol>
2122  <physvol>
2123  <volumeref ref="volRearGabion2"/>
2124  <position name="posRearGabion" unit="cm" x="($DetEncWidth+$ServiceBuildingExtraWidth+2*$LowBaySideWallThickness+$GabionThickness)/2" y="$DetEncHeight/2 +($LowBayInsideHeight+$LowBayCeilingThickness)/2 + $LowBayTopRockThickness/2" z="-(($LowBayInsideLength+$LowBayBackWallThickness)/2 + $HighBayOverlap+$HighBayBackWallThickness - $DetEncLength/2)"/>
2125  </physvol>
2126  <physvol>
2127  <volumeref ref="volHighBayTopRock"/>
2128  <position name="posHighBayTopRock" unit="cm" x="0" y="$DetEncHeight/2 + $HighBayInsideHeight + $HighBayCeilingThickness + $HighBayTopRockThickness/2" z="($DetEncLength/2+($HighBayInsideLength+$HighBayBackWallThickness+$HighBayFrontWallThickness)/2 -$HighBayBackWallThickness - $HighBayOverlap)"/>
2129  </physvol>
2130  <physvol>
2131  <volumeref ref="volLowBayTopRock"/>
2132  <position name="posLowBayTopRock" unit="cm" x="0" y="$DetEncHeight/2 + $LowBayInsideHeight + $LowBayCeilingThickness + $LowBayTopRockThickness/2" z="-(($LowBayInsideLength+$LowBayBackWallThickness)/2 + $HighBayOverlap+$HighBayBackWallThickness - $DetEncLength/2)"/>
2133  </physvol>
2134  <physvol>
2135  <volumeref ref="volSlopeRockFill"/>
2136  <position name="posSlopeRockFill" unit="cm" x="0" y="$DetEncHeight/2+($LowBayInsideHeight+$LowBayCeilingThickness) +$SlopeRockFillHeight/3" z="-($HighBayOverlap+$HighBayBackWallThickness - $DetEncLength/2)"/>
2137  <rotationref ref="rMinus90AboutY"/>
2138  </physvol>
2139  <physvol>
2140  <volumeref ref="volHillSideBox"/>
2141  <position name="posHillSideBox" unit="cm" x="0" y="-$ConcretePadding/2+$DetEncHeight/2 + $HillSideBoxHeight/2" z="$DetEncLength/2-($DetEncLength/2 +$HillSideBoxLength/2 + ($LowBayInsideLength +$LowBayBackWallThickness + $HighBayBackWallThickness + $HighBayOverlap - $DetEncLength ) )"/>
2142  </physvol>
2143 
2144  <physvol>
2145  <volumeref ref="volHillSideMiddle"/>
2146  <position name="posHillSideMiddle" unit="cm" x="0" y="$DetEncHeight/2 + $HillSideBoxHeight" z="-($DetEncLength/2 +$HillSideBoxLength + ($LowBayInsideLength +$LowBayBackWallThickness + $HighBayBackWallThickness + $HighBayOverlap - $DetEncLength ) )"/>
2147  <rotationref ref="rMinus90AboutYMinus90AboutX"/>
2148  </physvol>
2149 
2150  <physvol>
2151  <volumeref ref="volHillSide1"/>
2152  <position name="posHillSide" unit="cm" x="-($DetEncWidth+$ServiceBuildingExtraWidth+2*$LowBaySideWallThickness+2*$GabionThickness+$HillSideWidth)/2" y="$DetEncHeight/2" z="$RockThickness-$DetEncLength/2"/>
2153  <rotationref ref="rMinus90AboutYMinus90AboutX"/>
2154  </physvol>
2155  <physvol>
2156  <volumeref ref="volHillSide2"/>
2157  <position name="posHillSide" unit="cm" x="($DetEncWidth+$ServiceBuildingExtraWidth+2*$LowBaySideWallThickness+2*$GabionThickness+$HillSideWidth)/2" y="$DetEncHeight/2" z="$RockThickness-$DetEncLength/2"/>
2158  <rotationref ref="rMinus90AboutYMinus90AboutX"/>
2159  </physvol>
2160  </volume>
2161 </structure>
2162 </gdml>
2163 EOF
2164 
2165 # make_gdml.pl will take care of <setup/>
2166 
2167 close(HILLSIDE);
2168 }
2169 
2170 
2171 #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2172 #+++++++++++++++++++++++++++++++++++++++ gen_World +++++++++++++++++++++++++++++++++++++++
2173 #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2174 
2175 sub gen_World()
2176 {
2177 
2178 # Create the WORLD fragment file name,
2179 # add file to list of output GDML fragments,
2180 # and open it
2181  $WORLD = "dune_10kT_World" . $suffix . ".gdml";
2182  push (@gdmlFiles, $WORLD);
2183  $WORLD = ">" . $WORLD;
2184  open(WORLD) or die("Could not open file $WORLD for writing");
2185 
2186 
2187 # The standard XML prefix and starting the gdml
2188  print WORLD <<EOF;
2189 <?xml version='1.0'?>
2190 <gdml>
2191 EOF
2192 
2193 
2194 # All the World solids.
2195 print WORLD <<EOF;
2196 <solids>
2197  <box name="World" lunit="cm"
2198  x="$DetEncWidth+3*$RockThickness"
2199  y="$DetEncHeight+3*$RockThickness"
2200  z="$DetEncLength+3*$RockThickness"/>
2201 </solids>
2202 EOF
2203 
2204 # World structure
2205 print WORLD <<EOF;
2206 <structure>
2207  <volume name="volWorld" >
2208  <materialref ref="Air"/>
2209  <solidref ref="World"/>
2210 <!-- <physvol>
2211  <volumeref ref="volHighBay"/>
2212  <position name="posHighBay" unit="cm" x="0" y="$OriginYSet+$DetEncHeight/2+($HighBayInsideHeight+$HighBayCeilingThickness)/2" z="$OriginZSet+($DetEncLength/2+($HighBayInsideLength+$HighBayBackWallThickness+$HighBayFrontWallThickness)/2 -$HighBayBackWallThickness - $HighBayOverlap)"/>
2213  </physvol>
2214  <physvol>
2215  <volumeref ref="volLowBay"/>
2216  <position name="posLowBay" unit="cm" x="0" y="$OriginYSet+$DetEncHeight/2+($LowBayInsideHeight+$LowBayCeilingThickness)/2" z="$OriginZSet-(($LowBayInsideLength+$LowBayBackWallThickness)/2 + $HighBayOverlap+$HighBayBackWallThickness - $DetEncLength/2)"/>
2217  </physvol>
2218  <physvol>
2219  <volumeref ref="volGround"/>
2220  <position name="posGround" unit="cm" x="0" y="0" z="0"/>
2221  </physvol> -->
2222  <physvol>
2223  <volumeref ref="volDetEnclosure"/>
2224 
2225  <position name="posDetEnclosure" unit="cm" x="0" y="$OriginYSet" z="$OriginZSet"/>
2226  </physvol>
2227  </volume>
2228 </structure>
2229 </gdml>
2230 EOF
2231 
2232 # make_gdml.pl will take care of <setup/>
2233 
2234 close(WORLD);
2235 }
2236 
2237 
2238 
2239 #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2240 #++++++++++++++++++++++++++++++++++++ write_fragments ++++++++++++++++++++++++++++++++++++
2241 #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2242 
2243 sub write_fragments()
2244 {
2245  # This subroutine creates an XML file that summarizes the the subfiles output
2246  # by the other sub routines - it is the input file for make_gdml.pl which will
2247  # give the final desired GDML file. Specify its name with the output option.
2248  # (you can change the name when running make_gdml)
2249 
2250  # This code is taken straigh from the similar MicroBooNE generate script, Thank you.
2251 
2252  if ( ! defined $output )
2253  {
2254  $output = "-"; # write to STDOUT
2255  }
2256 
2257  # Set up the output file.
2258  $OUTPUT = ">" . $output;
2259  open(OUTPUT) or die("Could not open file $OUTPUT");
2260 
2261  print OUTPUT <<EOF;
2262 <?xml version='1.0'?>
2263 
2264 <!-- Input to Geometry/gdml/make_gdml.pl; define the GDML fragments
2265  that will be zipped together to create a detector description.
2266  -->
2267 
2268 <config>
2269 
2270  <constantfiles>
2271 
2272  <!-- These files contain GDML <constant></constant>
2273  blocks. They are read in separately, so they can be
2274  interpreted into the remaining GDML. See make_gdml.pl for
2275  more information.
2276  -->
2277 
2278 EOF
2279 
2280  foreach $filename (@defFiles)
2281  {
2282  print OUTPUT <<EOF;
2283  <filename> $filename </filename>
2284 EOF
2285  }
2286 
2287  print OUTPUT <<EOF;
2288 
2289  </constantfiles>
2290 
2291  <gdmlfiles>
2292 
2293  <!-- The GDML file fragments to be zipped together. -->
2294 
2295 EOF
2296 
2297  foreach $filename (@gdmlFiles)
2298  {
2299  print OUTPUT <<EOF;
2300  <filename> $filename </filename>
2301 EOF
2302  }
2303 
2304  print OUTPUT <<EOF;
2305 
2306  </gdmlfiles>
2307 
2308 </config>
2309 EOF
2310 
2311  close(OUTPUT);
2312 }
Generated by   doxygen 1.8.11