generate_dunevd10kt_3view_30deg_v1_refactored.pl
Go to the documentation of this file.
1 #!/usr/bin/perl
2 
3 #
4 #
5 # First attempt to make a GDML fragment generator for the DUNE vertical drift
6 # 10kt detector geometry with 3 views: +/- Xdeg for induction and 90 deg collection
7 # The lower chamber is not added yet.
8 # !!!NOTE!!!: the readout is on a positive Y plane (drift along horizontal X)
9 # due to current reco limitations)
10 # No photon detectors declared
11 # Simplified treatment of inter-module dead spaces
12 #
13 # Created: Thu Oct 1 16:45:27 CEST 2020
14 # Vyacheslav Galymov <vgalymov@ipnl.in2p3.fr>
15 #
16 # Modified:
17 # VG: Added defs to enable use in the refactored sim framework
18 # VG: 23.02.21 Adjust plane dimensions to fit a given number of ch per side
19 # VG: 23.02.21 Group CRUs in CRPs
20 # VG: 02.03.21 The length for the ROP is force to be the lenght
21 # given by nch_collection x pitch_collection
22 #
23 #################################################################################
24 
25 # Each subroutine generates a fragment GDML file, and the last subroutine
26 # creates an XML file that make_gdml.pl will use to appropriately arrange
27 # the fragment GDML files to create the final desired DUNE GDML file,
28 # to be named by make_gdml output command
29 
30 ##################################################################################
31 
32 
33 #use warnings;
34 use gdmlMaterials;
35 use Math::Trig;
36 use Getopt::Long;
37 use Math::BigFloat;
38 Math::BigFloat->precision(-16);
39 
40 ###
41 GetOptions( "help|h" => \$help,
42  "suffix|s:s" => \$suffix,
43  "output|o:s" => \$output,
44  "wires|w:s" => \$wires,
45  "workspace|k:s" => \$wkspc);
46 
47 my $FieldCage_switch="off";
48 my $Cathode_switch="off";
49 
50 if ( defined $help )
51 {
52  # If the user requested help, print the usage notes and exit.
53  usage();
54  exit;
55 }
56 
57 if ( ! defined $suffix )
58 {
59  # The user didn't supply a suffix, so append nothing to the file
60  # names.
61  $suffix = "";
62 }
63 else
64 {
65  # Otherwise, stick a "-" before the suffix, so that a suffix of
66  # "test" applied to filename.gdml becomes "filename-test.gdml".
67  $suffix = "-" . $suffix;
68 }
69 
70 
71 $workspace = 0;
72 if(defined $wkspc )
73 {
74  $workspace = $wkspc;
75 }
76 elsif ( $workspace != 0 )
77 {
78  print "\t\tCreating smaller workspace geometry.\n";
79 }
80 
81 # set wires on to be the default, unless given an input by the user
82 $wires_on = 1; # 1=on, 0=off
83 if (defined $wires)
84 {
85  $wires_on = $wires;
86 }
87 
88 $tpc_on = 1;
89 $basename="_";
90 
91 
92 ##################################################################
93 ############## Parameters for One Readout Panel ##################
94 
95 # parameters for 1.5 x 1.7 sub-unit Charge Readout Module / Unit
96 #$widthPCBActive = 169.0; # cm
97 #$lengthPCBActive = 150.0; # cm
98 
99 # views and channel counts
100 %nChans = ('Ind1', 298, 'Ind1Bot', 100, 'Ind2', 298, 'Col', 304);
101 $nViews = keys %nChans;
102 #print "$nViews %nChans\n";
103 
104 # first induction view
105 $wirePitchU = 0.7335; # cm
106 $wireAngleU = 150.0; # deg
107 
108 # second induction view
109 $wirePitchV = 0.7335; # cm
110 $wireAngleV = 30.0; # deg
111 
112 
113 # last collection view
114 $wirePitchZ = 0.489; # cm
115 
116 # force length to be equal to collection nch x pitch
117 $lengthPCBActive = $wirePitchZ * $nChans{'Col'};
118 $widthPCBActive = 167.7006;
119 
120 #
121 $borderCRM = 0.0; # border space aroud each CRM
122 
123 $widthCRM_active = $widthPCBActive;
124 $lengthCRM_active = $lengthPCBActive;
125 
126 $widthCRM = $widthPCBActive + 2 * $borderCRM;
127 $lengthCRM = $lengthPCBActive + 2 * $borderCRM;
128 
129 $borderCRP = 0.5; # cm
130 
131 # number of CRMs in y and z
132 $nCRM_x = 4 * 2;
133 $nCRM_z = 20 * 2;
134 
135 # create a smaller geometry
136 if( $workspace == 1 )
137 {
138  $nCRM_x = 1 * 2;
139  $nCRM_z = 1 * 2;
140 }
141 
142 # create a smaller geometry
143 if( $workspace == 2 )
144 {
145  $nCRM_x = 2 * 2;
146  $nCRM_z = 2 * 2;
147 }
148 
149 # create a smaller geometry
150 if( $workspace == 3 )
151 {
152  $nCRM_x = 3 * 2;
153  $nCRM_z = 3 * 2;
154 }
155 
156 
157 # calculate tpc area based on number of CRMs and their dimensions
158 # each CRP should have a 2x2 CRMs
159 $widthTPCActive = $nCRM_x * $widthCRM + $nCRM_x * $borderCRP; # around 1200 for full module
160 $lengthTPCActive = $nCRM_z * $lengthCRM + $nCRM_z * $borderCRP; # around 6000 for full module
161 
162 # active volume dimensions
163 $driftTPCActive = 650.0;
164 
165 # model anode strips as wires of some diameter
166 $padWidth = 0.02;
167 $ReadoutPlane = $nViews * $padWidth; # 3 readout planes (no space b/w)!
168 
169 ##################################################################
170 ############## Parameters for TPC and inner volume ###############
171 
172 # inner volume dimensions of the cryostat
173 $Argon_x = 1510;
174 $Argon_y = 1510;
175 $Argon_z = 6200;
176 
177 # width of gas argon layer on top
178 $HeightGaseousAr = 100;
179 
180 if( $workspace != 0 )
181 {
182  #active tpc + 1.0 m buffer on each side
183  $Argon_x = $driftTPCActive + $HeightGaseousAr + $ReadoutPlane + 100;
184  $Argon_y = $widthTPCActive + 200;
185  $Argon_z = $lengthTPCActive + 200;
186 }
187 
188 
189 # size of liquid argon buffer
190 $xLArBuffer = $Argon_x - $driftTPCActive - $HeightGaseousAr - $ReadoutPlane;
191 $yLArBuffer = 0.5 * ($Argon_y - $widthTPCActive);
192 $zLArBuffer = 0.5 * ($Argon_z - $lengthTPCActive);
193 
194 # cryostat
195 $SteelThickness = 0.12; # membrane
196 
197 $Cryostat_x = $Argon_x + 2*$SteelThickness;
198 $Cryostat_y = $Argon_y + 2*$SteelThickness;
199 $Cryostat_z = $Argon_z + 2*$SteelThickness;
200 
201 ##################################################################
202 ############## DetEnc and World relevant parameters #############
203 ##################################################################
204 ############## DetEnc and World relevant parameters #############
205 
206 $SteelSupport_x = 100;
207 $SteelSupport_y = 100;
208 $SteelSupport_z = 100;
209 $FoamPadding = 80;
210 $FracMassOfSteel = 0.5; #The steel support is not a solid block, but a mixture of air and steel
211 $FracMassOfAir = 1 - $FracMassOfSteel;
212 
213 
214 $SpaceSteelSupportToWall = 100;
215 $SpaceSteelSupportToCeiling = 100;
216 
217 $DetEncX = $Cryostat_x
218  + 2*($SteelSupport_x + $FoamPadding) + $SpaceSteelSupportToCeiling;
219 
220 $DetEncY = $Cryostat_y
221  + 2*($SteelSupport_y + $FoamPadding) + 2*$SpaceSteelSupportToWall;
222 
223 $DetEncZ = $Cryostat_z
224  + 2*($SteelSupport_z + $FoamPadding) + 2*$SpaceSteelSupportToWall;
225 
226 $posCryoInDetEnc_x = - $DetEncX/2 + $SteelSupport_x + $FoamPadding + $Cryostat_x/2;
227 
228 
229 $RockThickness = 4000;
230 
231  # We want the world origin to be vertically centered on active TPC
232  # This is to be added to the x and y position of every volume in volWorld
233 
234 $OriginXSet = $DetEncX/2.0
235  - $SteelSupport_x
236  - $FoamPadding
237  - $SteelThickness
238  - $xLArBuffer
239  - $driftTPCActive/2.0;
240 
241 $OriginYSet = $DetEncY/2.0
242  - $SpaceSteelSupportToWall
243  - $SteelSupport_y
244  - $FoamPadding
245  - $SteelThickness
246  - $yLArBuffer
247  - $widthTPCActive/2.0;
248 
249  # We want the world origin to be at the very front of the fiducial volume.
250  # move it to the front of the enclosure, then back it up through the concrete/foam,
251  # then through the Cryostat shell, then through the upstream dead LAr (including the
252  # dead LAr on the edge of the TPC)
253  # This is to be added to the z position of every volume in volWorld
254 
255 $OriginZSet = $DetEncZ/2.0
256  - $SpaceSteelSupportToWall
257  - $SteelSupport_z
258  - $FoamPadding
259  - $SteelThickness
260  - $zLArBuffer
261  - $borderCRM;
262 
263 ##################################################################
264 ############## Field Cage Parameters ###############
265 $FieldShaperLongTubeLength = $lengthTPCActive;
266 $FieldShaperShortTubeLength = $widthTPCActive;
267 $FieldShaperInnerRadius = 1.485;
268 $FieldShaperOuterRadius = 1.685;
269 $FieldShaperTorRad = 1.69;
270 
271 $FieldShaperLength = $FieldShaperLongTubeLength + 2*$FieldShaperOuterRadius+ 2*$FieldShaperTorRad;
272 $FieldShaperWidth = $FieldShaperShortTubeLength + 2*$FieldShaperOuterRadius+ 2*$FieldShaperTorRad;
273 
274 $FieldShaperSeparation = 5.0;
275 $NFieldShapers = ($driftTPCActive/$FieldShaperSeparation) - 1;
276 
277 $FieldCageSizeX = $FieldShaperSeparation*$NFieldShapers+2;
278 $FieldCageSizeY = $FieldShaperWidth+2;
279 $FieldCageSizeZ = $FieldShaperLength+2;
280 
281 
282 #+++++++++++++++++++++++++ End defining variables ++++++++++++++++++++++++++
283 
284 
285 #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
286 #+++++++++++++++++++++++++++++++++++++++++ usage +++++++++++++++++++++++++++++++++++++++++
287 #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
288 
289 sub usage()
290 {
291  print "Usage: $0 [-h|--help] [-o|--output <fragments-file>] [-s|--suffix <string>]\n";
292  print " if -o is omitted, output goes to STDOUT; <fragments-file> is input to make_gdml.pl\n";
293  print " -s <string> appends the string to the file names; useful for multiple detector versions\n";
294  print " -h prints this message, then quits\n";
295 }
296 
297 
298 sub gen_Extend()
299 {
300 
301 # Create the <define> fragment file name,
302 # add file to list of fragments,
303 # and open it
304  $DEF = $basename."_Ext" . $suffix . ".gdml";
305  push (@gdmlFiles, $DEF);
306  $DEF = ">" . $DEF;
307  open(DEF) or die("Could not open file $DEF for writing");
308 
309 print DEF <<EOF;
310 <?xml version='1.0'?>
311 <gdml>
312 <extension>
313  <color name="magenta" R="0.0" G="1.0" B="0.0" A="1.0" />
314  <color name="green" R="0.0" G="1.0" B="0.0" A="1.0" />
315  <color name="red" R="1.0" G="0.0" B="0.0" A="1.0" />
316  <color name="blue" R="0.0" G="0.0" B="1.0" A="1.0" />
317  <color name="yellow" R="1.0" G="1.0" B="0.0" A="1.0" />
318 </extension>
319 </gdml>
320 EOF
321  close (DEF);
322 }
323 
324 #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
325 #++++++++++++++++++++++++++++++++++++++ gen_Define +++++++++++++++++++++++++++++++++++++++
326 #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
327 
328 sub gen_Define()
329 {
330 
331 # Create the <define> fragment file name,
332 # add file to list of fragments,
333 # and open it
334  $DEF = $basename."_Def" . $suffix . ".gdml";
335  push (@gdmlFiles, $DEF);
336  $DEF = ">" . $DEF;
337  open(DEF) or die("Could not open file $DEF for writing");
338 
339 
340 print DEF <<EOF;
341 <?xml version='1.0'?>
342 <gdml>
343 <define>
344 
345 <!--
346 
347 
348 
349 -->
350 
351  <position name="posCryoInDetEnc" unit="cm" x="$posCryoInDetEnc_x" y="0" z="0"/>
352  <position name="posCenter" unit="cm" x="0" y="0" z="0"/>
353  <rotation name="rUWireAboutX" unit="deg" x="$wireAngleU" y="0" z="0"/>
354  <rotation name="rVWireAboutX" unit="deg" x="$wireAngleV" y="0" z="0"/>
355  <rotation name="rPlus90AboutX" unit="deg" x="90" y="0" z="0"/>
356  <rotation name="rPlus90AboutY" unit="deg" x="0" y="90" z="0"/>
357  <rotation name="rPlus90AboutXPlus90AboutY" unit="deg" x="90" y="90" z="0"/>
358  <rotation name="rMinus90AboutX" unit="deg" x="270" y="0" z="0"/>
359  <rotation name="rMinus90AboutY" unit="deg" x="0" y="270" z="0"/>
360  <rotation name="rMinus90AboutYMinus90AboutX" unit="deg" x="270" y="270" z="0"/>
361  <rotation name="rPlus180AboutX" unit="deg" x="180" y="0" z="0"/>
362  <rotation name="rPlus180AboutY" unit="deg" x="0" y="180" z="0"/>
363  <rotation name="rPlus180AboutXPlus180AboutY" unit="deg" x="180" y="180" z="0"/>
364  <rotation name="rIdentity" unit="deg" x="0" y="0" z="0"/>
365 </define>
366 </gdml>
367 EOF
368  close (DEF);
369 }
370 
371 #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
372 #+++++++++++++++++++++++++++++++++++++ gen_Materials +++++++++++++++++++++++++++++++++++++
373 #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
374 
375 sub gen_Materials()
376 {
377 
378 # Create the <materials> fragment file name,
379 # add file to list of output GDML fragments,
380 # and open it
381  $MAT = $basename."_Materials" . $suffix . ".gdml";
382  push (@gdmlFiles, $MAT);
383  $MAT = ">" . $MAT;
384 
385  open(MAT) or die("Could not open file $MAT for writing");
386 
387  # Add any materials special to this geometry by defining a mulitline string
388  # and passing it to the gdmlMaterials::gen_Materials() function.
389 my $asmix = <<EOF;
390  <!-- preliminary values -->
391  <material name="AirSteelMixture" formula="AirSteelMixture">
392  <D value=" 0.001205*(1-$FracMassOfSteel) + 7.9300*$FracMassOfSteel " unit="g/cm3"/>
393  <fraction n="$FracMassOfSteel" ref="STEEL_STAINLESS_Fe7Cr2Ni"/>
394  <fraction n="$FracMassOfAir" ref="Air"/>
395  </material>
396  <material name="vm2000" formula="vm2000">
397  <D value="1.2" unit="g/cm3"/>
398  <composite n="2" ref="carbon"/>
399  <composite n="4" ref="hydrogen"/>
400  </material>
401 EOF
402 
403  # add the general materials used anywere
404  print MAT gdmlMaterials::gen_Materials( $asmix );
405 
406  close(MAT);
407 }
408 
409 
410 #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
411 #++++++++++++++++++++++++++++++++++++++++ gen_TPC ++++++++++++++++++++++++++++++++++++++++
412 #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
413 # line clip on the rectangle boundary
414 sub lineClip {
415  my $x0 = $_[0];
416  my $y0 = $_[1];
417  my $nx = $_[2];
418  my $ny = $_[3];
419  my $rcl = $_[4];
420  my $rcw = $_[5];
421 
422  my $tol = 1.0E-4;
423  my @endpts = ();
424  if( abs( nx ) < tol ){
425  push( @endpts, ($x0, 0) );
426  push( @endpts, ($x0, $rcw) );
427  return @endpts;
428  }
429  if( abs( ny ) < tol ){
430  push( @endpts, (0, $y0) );
431  push( @endpts, ($rcl, $y0) );
432  return @endpts;
433  }
434 
435  # left border at x = 0
436  my $y = $y0 - $x0 * $ny/$nx;
437  if( $y >= 0 && $y <= $rcw ){
438  push( @endpts, (0, $y) );
439  }
440 
441  # right border at x = l
442  $y = $y0 + ($rcl-$x0) * $ny/$nx;
443  if( $y >= 0 && $y <= $rcw ){
444  push( @endpts, ($rcl, $y) );
445  if( scalar(@endpts) == 4 ){
446  return @endpts;
447  }
448  }
449 
450  # bottom border at y = 0
451  my $x = $x0 - $y0 * $nx/$ny;
452  if( $x >= 0 && $x <= $rcl ){
453  push( @endpts, ($x, 0) );
454  if( scalar(@endpts) == 4 ){
455  return @endpts;
456  }
457  }
458 
459  # top border at y = w
460  $x = $x0 + ($rcw-$y0)* $nx/$ny;
461  if( $x >= 0 && $x <= $rcl ){
462  push( @endpts, ($x, $rcw) );
463  }
464 
465  return @endpts;
466 }
467 
468 sub gen_Wires
469 {
470  my $length = $_[0]; #
471  my $width = $_[1]; #
472  my $nch = $_[2]; #
473  my $nchb = $_[3]; # nch per bottom side
474  my $pitch = $_[4]; #
475  my $theta = $_[5]; # deg
476  my $dia = $_[6]; #
477 
478  $theta = $theta * pi()/180.0;
479  my @dirw = (cos($theta), sin($theta));
480  my @dirp = (cos($theta - pi()/2), sin($theta - pi()/2));
481 
482  # calculate
483  my $alpha = $theta;
484  if( $alpha > pi()/2 ){
485  $alpha = pi() - $alpha;
486  }
487  my $dX = $pitch / sin( $alpha );
488  my $dY = $pitch / sin( pi()/2 - $alpha );
489  if( $length <= 0 ){
490  $length = $dX * $nchb;
491  }
492  if( $width <= 0 ){
493  $width = $dY * ($nch - $nchb);
494  }
495 
496  my @orig = (0, 0);
497  if( $dirp[0] < 0 ){
498  $orig[0] = $length;
499  }
500  if( $dirp[1] < 0 ){
501  $orig[1] = $width;
502  }
503 
504  #print "origin : @orig\n";
505  #print "pitch dir : @dirp\n";
506  #print "wire dir : @dirw\n";
507  #print "$length x $width cm2\n";
508 
509  # gen wires
510  my @winfo = ();
511  my $offset = $pitch/2;
512  foreach my $ch (0..$nch-1){
513  #print "Processing $ch\n";
514 
515  # calculate reference point for this strip
516  my @wcn = (0, 0);
517  $wcn[0] = $orig[0] + $offset * $dirp[0];
518  $wcn[1] = $orig[1] + $offset * $dirp[1];
519 
520  # line clip on the rectangle boundary
521  @endpts = lineClip( $wcn[0], $wcn[1], $dirw[0], $dirw[1], $length, $width );
522 
523  if( scalar(@endpts) != 4 ){
524  print "Could not find end points for wire $ch : @endpts\n";
525  $offset = $offset + $pitch;
526  next;
527  }
528 
529  # re-center on the mid-point
530  $endpts[0] -= $length/2;
531  $endpts[2] -= $length/2;
532  $endpts[1] -= $width/2;
533  $endpts[3] -= $width/2;
534 
535  # calculate the strip center in the rectangle of CRU
536  $wcn[0] = ($endpts[0] + $endpts[2])/2;
537  $wcn[1] = ($endpts[1] + $endpts[3])/2;
538 
539  # calculate the length
540  my $dx = $endpts[0] - $endpts[2];
541  my $dy = $endpts[1] - $endpts[3];
542  my $wlen = sqrt($dx**2 + $dy**2);
543 
544  # put all info together
545  my @wire = ($ch, $wcn[0], $wcn[1], $wlen);
546  push( @wire, @endpts );
547  push( @winfo, \@wire);
548  $offset = $offset + $pitch;
549  #last;
550  }
551  return @winfo;
552 }
553 
554 #
555 sub gen_TPC()
556 {
557  # CRM active volume
558  my $TPCActive_x = $driftTPCActive;
559  my $TPCActive_y = $widthCRM_active;
560  my $TPCActive_z = $lengthCRM_active;
561 
562  # CRM total volume
563  my $TPC_x = $TPCActive_x + $ReadoutPlane;
564  my $TPC_y = $widthCRM;
565  my $TPC_z = $lengthCRM;
566 
567  print " TPC dimensions : $TPC_x x $TPC_y x $TPC_z\n";
568 
569  $TPC = $basename."_TPC" . $suffix . ".gdml";
570  push (@gdmlFiles, $TPC);
571  $TPC = ">" . $TPC;
572  open(TPC) or die("Could not open file $TPC for writing");
573 
574  # The standard XML prefix and starting the gdml
575 print TPC <<EOF;
576  <?xml version='1.0'?>
577  <gdml>
578 EOF
579 
580  # compute wires for 1st induction
581  my @winfoU = ();
582  my @winfoV = ();
583  if( $wires_on == 1 ){
584  @winfoU = gen_Wires( $TPCActive_z, 0, # force length
585  $nChans{'Ind1'}, $nChans{'Ind1Bot'},
586  $wirePitchU, $wireAngleU, $padWidth );
587  @winfoV = gen_Wires( $TPCActive_z, 0, # force length
588  $nChans{'Ind2'}, $nChans{'Ind1Bot'},
589  $wirePitchV, $wireAngleV, $padWidth );
590 
591  }
592 
593  # All the TPC solids save the wires.
594 print TPC <<EOF;
595  <solids>
596 EOF
597 
598 print TPC <<EOF;
599  <box name="CRM"
600  x="$TPC_x"
601  y="$TPC_y"
602  z="$TPC_z"
603  lunit="cm"/>
604  <box name="CRMUPlane"
605  x="$padWidth"
606  y="$TPCActive_y"
607  z="$TPCActive_z"
608  lunit="cm"/>
609  <box name="CRMVPlane"
610  x="$padWidth"
611  y="$TPCActive_y"
612  z="$TPCActive_z"
613  lunit="cm"/>
614  <box name="CRMZPlane"
615  x="$padWidth"
616  y="$TPCActive_y"
617  z="$TPCActive_z"
618  lunit="cm"/>
619  <box name="CRMActive"
620  x="$TPCActive_x"
621  y="$TPCActive_y"
622  z="$TPCActive_z"
623  lunit="cm"/>
624 EOF
625 
626 #++++++++++++++++++++++++++++ Wire Solids ++++++++++++++++++++++++++++++
627 if($wires_on==1){
628 
629  foreach my $wire (@winfoU) {
630  my $wid = $wire->[0];
631  my $wln = $wire->[3];
632 print TPC <<EOF;
633  <tube name="CRMWireU$wid"
634  rmax="0.5*$padWidth"
635  z="$wln"
636  deltaphi="360"
637  aunit="deg" lunit="cm"/>
638 EOF
639  }
640 
641  foreach my $wire (@winfoV) {
642  my $wid = $wire->[0];
643  my $wln = $wire->[3];
644 print TPC <<EOF;
645  <tube name="CRMWireV$wid"
646  rmax="0.5*$padWidth"
647  z="$wln"
648  deltaphi="360"
649  aunit="deg" lunit="cm"/>
650 EOF
651  }
652 
653 
654 print TPC <<EOF;
655  <tube name="CRMWireZ"
656  rmax="0.5*$padWidth"
657  z="$TPCActive_y"
658  deltaphi="360"
659  aunit="deg" lunit="cm"/>
660 EOF
661 }
662 print TPC <<EOF;
663 </solids>
664 EOF
665 
666 
667 # Begin structure and create wire logical volumes
668 print TPC <<EOF;
669 <structure>
670  <volume name="volTPCActive">
671  <materialref ref="LAr"/>
672  <solidref ref="CRMActive"/>
673  <auxiliary auxtype="SensDet" auxvalue="SimEnergyDeposit"/>
674  <auxiliary auxtype="StepLimit" auxunit="cm" auxvalue="0.5208*cm"/>
675  <auxiliary auxtype="Efield" auxunit="V/cm" auxvalue="500*V/cm"/>
676  <colorref ref="blue"/>
677  </volume>
678 EOF
679 
680 if($wires_on==1)
681 {
682  foreach my $wire (@winfoU)
683  {
684  my $wid = $wire->[0];
685 print TPC <<EOF;
686  <volume name="volTPCWireU$wid">
687  <materialref ref="Copper_Beryllium_alloy25"/>
688  <solidref ref="CRMWireU$wid"/>
689  </volume>
690 EOF
691  }
692 
693  foreach my $wire (@winfoV)
694  {
695  my $wid = $wire->[0];
696 print TPC <<EOF;
697  <volume name="volTPCWireV$wid">
698  <materialref ref="Copper_Beryllium_alloy25"/>
699  <solidref ref="CRMWireV$wid"/>
700  </volume>
701 EOF
702  }
703 
704 print TPC <<EOF;
705  <volume name="volTPCWireZ">
706  <materialref ref="Copper_Beryllium_alloy25"/>
707  <solidref ref="CRMWireZ"/>
708  </volume>
709 EOF
710 }
711  # 1st induction plane
712 print TPC <<EOF;
713  <volume name="volTPCPlaneU">
714  <materialref ref="LAr"/>
715  <solidref ref="CRMUPlane"/>
716 EOF
717 if ($wires_on==1) # add wires to U plane
718 {
719  # the coordinates were computed with a corner at (0,0)
720  # so we need to move to plane coordinates
721  my $offsetZ = 0; #-0.5 * $TPCActive_z;
722  my $offsetY = 0; #-0.5 * $TPCActive_y;
723 
724  foreach my $wire (@winfoU) {
725  my $wid = $wire->[0];
726  my $zpos = $wire->[1] + $offsetZ;
727  my $ypos = $wire->[2] + $offsetY;
728 print TPC <<EOF;
729  <physvol>
730  <volumeref ref="volTPCWireU$wid"/>
731  <position name="posWireU$wid" unit="cm" x="0" y="$ypos" z="$zpos"/>
732  <rotationref ref="rUWireAboutX"/>
733  </physvol>
734 EOF
735  }
736 }
737 print TPC <<EOF;
738  </volume>
739 EOF
740 
741 # 2nd induction plane
742 print TPC <<EOF;
743  <volume name="volTPCPlaneV">
744  <materialref ref="LAr"/>
745  <solidref ref="CRMVPlane"/>
746 EOF
747 
748 if ($wires_on==1) # add wires to V plane (plane with wires reading y position)
749  {
750  # the coordinates were computed with a corner at (0,0)
751  # so we need to move to plane coordinates
752  my $offsetZ = 0; #-0.5 * $TPCActive_z;
753  my $offsetY = 0; #-0.5 * $TPCActive_y;
754 
755  foreach my $wire (@winfoV) {
756  my $wid = $wire->[0];
757  my $zpos = $wire->[1] + $offsetZ;
758  my $ypos = $wire->[2] + $offsetY;
759 print TPC <<EOF;
760  <physvol>
761  <volumeref ref="volTPCWireV$wid"/>
762  <position name="posWireV$wid" unit="cm" x="0" y="$ypos" z="$zpos"/>
763  <rotationref ref="rVWireAboutX"/>
764  </physvol>
765 EOF
766  }
767 }
768 print TPC <<EOF;
769  </volume>
770 EOF
771 
772 # collection plane
773 print TPC <<EOF;
774  <volume name="volTPCPlaneZ">
775  <materialref ref="LAr"/>
776  <solidref ref="CRMZPlane"/>
777 EOF
778 if ($wires_on==1) # add wires to Z plane (plane with wires reading z position)
779  {
780  for($i=0;$i<$nChans{'Col'};++$i)
781  {
782  #my $zpos = -0.5 * $TPCActive_z + ($i+0.5)*$wirePitchZ + 0.5*$padWidth;
783  my $zpos = ($i + 0.5 - $nChans{'Col'}/2)*$wirePitchZ;
784  if( (0.5 * $TPCActive_z - abs($zpos)) < 0 ){
785  die "Cannot place wire $i in view Z, as plane is too small\n";
786  }
787 print TPC <<EOF;
788  <physvol>
789  <volumeref ref="volTPCWireZ"/>
790  <position name="posWireZ$i" unit="cm" x="0" y="0" z="$zpos"/>
791  <rotationref ref="rPlus90AboutX"/>
792  </physvol>
793 EOF
794  }
795 }
796 print TPC <<EOF;
797  </volume>
798 EOF
799 
800 
801 $posUplane[0] = 0.5*$TPC_x - 2.5*$padWidth;
802 $posUplane[1] = 0;
803 $posUplane[2] = 0;
804 
805 $posVplane[0] = 0.5*$TPC_x - 1.5*$padWidth;
806 $posVplane[1] = 0;
807 $posVplane[2] = 0;
808 
809 $posZplane[0] = 0.5*$TPC_x - 0.5*$padWidth;
810 $posZplane[1] = 0;
811 $posZplane[2] = 0;
812 
813 $posTPCActive[0] = -$ReadoutPlane/2;
814 $posTPCActive[1] = 0;
815 $posTPCActive[2] = 0;
816 
817 
818 #wrap up the TPC file
819 print TPC <<EOF;
820  <volume name="volTPC">
821  <materialref ref="LAr"/>
822  <solidref ref="CRM"/>
823  <physvol>
824  <volumeref ref="volTPCPlaneU"/>
825  <position name="posPlaneU" unit="cm"
826  x="$posUplane[0]" y="$posUplane[1]" z="$posUplane[2]"/>
827  <rotationref ref="rIdentity"/>
828  </physvol>
829  <physvol>
830  <volumeref ref="volTPCPlaneV"/>
831  <position name="posPlaneY" unit="cm"
832  x="$posVplane[0]" y="$posVplane[1]" z="$posVplane[2]"/>
833  <rotationref ref="rIdentity"/>
834  </physvol>
835  <physvol>
836  <volumeref ref="volTPCPlaneZ"/>
837  <position name="posPlaneZ" unit="cm"
838  x="$posZplane[0]" y="$posZplane[1]" z="$posZplane[2]"/>
839  <rotationref ref="rIdentity"/>
840  </physvol>
841  <physvol>
842  <volumeref ref="volTPCActive"/>
843  <position name="posActive" unit="cm"
844  x="$posTPCActive[0]" y="$posTPCAtive[1]" z="$posTPCActive[2]"/>
845  <rotationref ref="rIdentity"/>
846  </physvol>
847  </volume>
848 EOF
849 
850 print TPC <<EOF;
851  </structure>
852  </gdml>
853 EOF
854 
855  close(TPC);
856 }
857 
858 #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
859 #++++++++++++++++++++++++++++++++++++++ gen_FieldCage ++++++++++++++++++++++++++++++++++++
860 #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
861 
862 sub gen_FieldCage {
863 
864  $FieldCage = $basename."_FieldCage" . $suffix . ".gdml";
865  push (@gdmlFiles, $FieldCage);
866  $FieldCage = ">" . $FieldCage;
867  open(FieldCage) or die("Could not open file $FieldCage for writing");
868 
869 # The standard XML prefix and starting the gdml
870 print FieldCage <<EOF;
871  <?xml version='1.0'?>
872  <gdml>
873 EOF
874 # The printing solids used in the Field Cage
875 #print "lengthTPCActive : $lengthTPCActive \n";
876 #print "widthTPCActive : $widthTPCActive \n";
877 
878 
879 print FieldCage <<EOF;
880 <solids>
881  <torus name="FieldShaperCorner" rmin="$FieldShaperInnerRadius" rmax="$FieldShaperOuterRadius" rtor="$FieldShaperTorRad" deltaphi="90" startphi="0" aunit="deg" lunit="cm"/>
882  <tube name="FieldShaperLongtube" rmin="$FieldShaperInnerRadius" rmax="$FieldShaperOuterRadius" z="$FieldShaperLongTubeLength" deltaphi="360" startphi="0" aunit="deg" lunit="cm"/>
883  <tube name="FieldShaperShorttube" rmin="$FieldShaperInnerRadius" rmax="$FieldShaperOuterRadius" z="$FieldShaperShortTubeLength" deltaphi="360" startphi="0" aunit="deg" lunit="cm"/>
884 
885  <union name="FSunion1">
886  <first ref="FieldShaperLongtube"/>
887  <second ref="FieldShaperCorner"/>
888  <position name="esquinapos1" unit="cm" x="@{[-$FieldShaperTorRad]}" y="0" z="@{[0.5*$FieldShaperLongTubeLength]}"/>
889  <rotation name="rot1" unit="deg" x="90" y="0" z="0" />
890  </union>
891 
892  <union name="FSunion2">
893  <first ref="FSunion1"/>
894  <second ref="FieldShaperShorttube"/>
895  <position name="esquinapos2" unit="cm" x="@{[-0.5*$FieldShaperShortTubeLength-$FieldShaperTorRad]}" y="0" z="@{[+0.5*$FieldShaperLongTubeLength+$FieldShaperTorRad]}"/>
896  <rotation name="rot2" unit="deg" x="0" y="90" z="0" />
897  </union>
898 
899  <union name="FSunion3">
900  <first ref="FSunion2"/>
901  <second ref="FieldShaperCorner"/>
902  <position name="esquinapos3" unit="cm" x="@{[-$FieldShaperShortTubeLength-$FieldShaperTorRad]}" y="0" z="@{[0.5*$FieldShaperLongTubeLength]}"/>
903  <rotation name="rot3" unit="deg" x="90" y="270" z="0" />
904  </union>
905 
906  <union name="FSunion4">
907  <first ref="FSunion3"/>
908  <second ref="FieldShaperLongtube"/>
909  <position name="esquinapos4" unit="cm" x="@{[-$FieldShaperShortTubeLength-2*$FieldShaperTorRad]}" y="0" z="0"/>
910  </union>
911 
912  <union name="FSunion5">
913  <first ref="FSunion4"/>
914  <second ref="FieldShaperCorner"/>
915  <position name="esquinapos5" unit="cm" x="@{[-$FieldShaperShortTubeLength-$FieldShaperTorRad]}" y="0" z="@{[-0.5*$FieldShaperLongTubeLength]}"/>
916  <rotation name="rot5" unit="deg" x="90" y="180" z="0" />
917  </union>
918 
919  <union name="FSunion6">
920  <first ref="FSunion5"/>
921  <second ref="FieldShaperShorttube"/>
922  <position name="esquinapos6" unit="cm" x="@{[-0.5*$FieldShaperShortTubeLength-$FieldShaperTorRad]}" y="0" z="@{[-0.5*$FieldShaperLongTubeLength-$FieldShaperTorRad]}"/>
923  <rotation name="rot6" unit="deg" x="0" y="90" z="0" />
924  </union>
925 
926  <union name="FieldShaperSolid">
927  <first ref="FSunion6"/>
928  <second ref="FieldShaperCorner"/>
929  <position name="esquinapos7" unit="cm" x="@{[-$FieldShaperTorRad]}" y="0" z="@{[-0.5*$FieldShaperLongTubeLength]}"/>
930  <rotation name="rot7" unit="deg" x="90" y="90" z="0" />
931  </union>
932 </solids>
933 
934 EOF
935 
936 print FieldCage <<EOF;
937 
938 <structure>
939 <volume name="volFieldShaper">
940  <materialref ref="Al2O3"/>
941  <solidref ref="FieldShaperSolid"/>
942 </volume>
943 </structure>
944 
945 EOF
946 
947 print FieldCage <<EOF;
948 
949 </gdml>
950 EOF
951 close(FieldCage);
952 }
953 
954 
955 
956 #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
957 #++++++++++++++++++++++++++++++++++++++ gen_Cryostat +++++++++++++++++++++++++++++++++++++
958 #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
959 
960 sub gen_Cryostat()
961 {
962 
963 # Create the cryostat fragment file name,
964 # add file to list of output GDML fragments,
965 # and open it
966  $CRYO = $basename."_Cryostat" . $suffix . ".gdml";
967  push (@gdmlFiles, $CRYO);
968  $CRYO = ">" . $CRYO;
969  open(CRYO) or die("Could not open file $CRYO for writing");
970 
971 
972 # The standard XML prefix and starting the gdml
973  print CRYO <<EOF;
974 <?xml version='1.0'?>
975 <gdml>
976 EOF
977 
978 # All the cryostat solids.
979 print CRYO <<EOF;
980 <solids>
981  <box name="Cryostat" lunit="cm"
982  x="$Cryostat_x"
983  y="$Cryostat_y"
984  z="$Cryostat_z"/>
985 
986  <box name="ArgonInterior" lunit="cm"
987  x="$Argon_x"
988  y="$Argon_y"
989  z="$Argon_z"/>
990 
991  <box name="GaseousArgon" lunit="cm"
992  x="$HeightGaseousAr"
993  y="$Argon_y"
994  z="$Argon_z"/>
995 
996  <subtraction name="SteelShell">
997  <first ref="Cryostat"/>
998  <second ref="ArgonInterior"/>
999  </subtraction>
1000 
1001 </solids>
1002 EOF
1003 
1004 # Cryostat structure
1005 print CRYO <<EOF;
1006 <structure>
1007  <volume name="volSteelShell">
1008  <materialref ref="STEEL_STAINLESS_Fe7Cr2Ni" />
1009  <solidref ref="SteelShell" />
1010  </volume>
1011  <volume name="volGaseousArgon">
1012  <materialref ref="ArGas"/>
1013  <solidref ref="GaseousArgon"/>
1014 EOF
1015 
1016  print CRYO <<EOF;
1017  </volume>
1018 
1019  <volume name="volCryostat">
1020  <materialref ref="LAr" />
1021  <solidref ref="Cryostat" />
1022  <physvol>
1023  <volumeref ref="volGaseousArgon"/>
1024  <position name="posGaseousArgon" unit="cm" x="@{[$Argon_x/2-$HeightGaseousAr/2]}" y="0" z="0"/>
1025  </physvol>
1026  <physvol>
1027  <volumeref ref="volSteelShell"/>
1028  <position name="posSteelShell" unit="cm" x="0" y="0" z="0"/>
1029  </physvol>
1030 EOF
1031 
1032 if ($tpc_on==1) # place TPC inside croysotat offsetting each pair of CRMs by borderCRP
1033 {
1034  $posX = $Argon_x/2 - $HeightGaseousAr - 0.5*($driftTPCActive + $ReadoutPlane);
1035  $idx = 0;
1036  my $posZ = -0.5*$Argon_z + $zLArBuffer + 0.5*$lengthCRM;
1037  for(my $ii=0;$ii<$nCRM_z;$ii++)
1038  {
1039  if( $ii % 2 == 0 ){
1040  $posZ += $borderCRP;
1041  if( $ii>0 ){
1042  $posZ += $borderCRP;
1043  }
1044  }
1045  my $posY = -0.5*$Argon_y + $yLArBuffer + 0.5*$widthCRM;
1046  for(my $jj=0;$jj<$nCRM_x;$jj++)
1047  {
1048  if( $jj % 2 == 0 ){
1049  $posY += $borderCRP;
1050  if( $jj>0 ){
1051  $posY += $borderCRP;
1052  }
1053  }
1054  print CRYO <<EOF;
1055  <physvol>
1056  <volumeref ref="volTPC"/>
1057  <position name="posTPC\-$idx" unit="cm"
1058  x="$posX" y="$posY" z="$posZ"/>
1059  </physvol>
1060 EOF
1061  $idx++;
1062  $posY += $widthCRM;
1063  }
1064 
1065  $posZ += $lengthCRM;
1066  }
1067 }
1068 
1069 #The +50 in the x positions must depend on some other parameter
1070  if ( $FieldCage_switch eq "on" ) {
1071  for ( $i=0; $i<$NFieldShapers; $i=$i+1 ) { # pmts with coating
1072 $posX = $Argon_x/2 - $HeightGaseousAr - 0.5*($driftTPCActive + $ReadoutPlane);
1073  print CRYO <<EOF;
1074  <physvol>
1075  <volumeref ref="volFieldShaper"/>
1076  <position name="posFieldShaper$i" unit="cm" x="@{[-$OriginXSet+50+($i-$NFieldShapers*0.5)*$FieldShaperSeparation]}" y="@{[-0.5*$FieldShaperShortTubeLength-$FieldShaperTorRad]}" z="0" />
1077  <rotation name="rotFS$i" unit="deg" x="0" y="0" z="90" />
1078  </physvol>
1079 EOF
1080  }
1081  }
1082 
1083 
1084 $CathodePosX =-$OriginXSet+50+(-1-$NFieldShapers*0.5)*$FieldShaperSeparation;
1085 $CathodePosY = 0;
1086 $CathodePosZ = 0;
1087  if ( $Cathode_switch eq "on" )
1088  {
1089  print CRYO <<EOF;
1090  <physvol>
1091  <volumeref ref="volGroundGrid"/>
1092  <position name="posGroundGrid01" unit="cm" x="$CathodePosX" y="@{[-$CathodePosY]}" z="@{[$CathodePosZ]}"/>
1093  <rotation name="rotGG01" unit="deg" x="0" y="0" z="90" />
1094  </physvol>
1095 
1096 EOF
1097  }
1098 
1099 
1100  print CRYO <<EOF;
1101  </volume>
1102 </structure>
1103 </gdml>
1104 EOF
1105 
1106 close(CRYO);
1107 }
1108 
1109 
1110 
1111 #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1112 #+++++++++++++++++++++++++++++++++++++ gen_Enclosure +++++++++++++++++++++++++++++++++++++
1113 #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1114 
1115 sub gen_Enclosure()
1116 {
1117 
1118 # Create the detector enclosure fragment file name,
1119 # add file to list of output GDML fragments,
1120 # and open it
1121  $ENCL = $basename."_DetEnclosure" . $suffix . ".gdml";
1122  push (@gdmlFiles, $ENCL);
1123  $ENCL = ">" . $ENCL;
1124  open(ENCL) or die("Could not open file $ENCL for writing");
1125 
1126 
1127 # The standard XML prefix and starting the gdml
1128  print ENCL <<EOF;
1129 <?xml version='1.0'?>
1130 <gdml>
1131 EOF
1132 
1133 
1134 # All the detector enclosure solids.
1135 print ENCL <<EOF;
1136 <solids>
1137 
1138  <box name="FoamPadBlock" lunit="cm"
1139  x="@{[$Cryostat_x + 2*$FoamPadding]}"
1140  y="@{[$Cryostat_y + 2*$FoamPadding]}"
1141  z="@{[$Cryostat_z + 2*$FoamPadding]}" />
1142 
1143  <subtraction name="FoamPadding">
1144  <first ref="FoamPadBlock"/>
1145  <second ref="Cryostat"/>
1146  <positionref ref="posCenter"/>
1147  </subtraction>
1148 
1149  <box name="SteelSupportBlock" lunit="cm"
1150  x="@{[$Cryostat_x + 2*$FoamPadding + 2*$SteelSupport_x]}"
1151  y="@{[$Cryostat_y + 2*$FoamPadding + 2*$SteelSupport_y]}"
1152  z="@{[$Cryostat_z + 2*$FoamPadding + 2*$SteelSupport_z]}" />
1153 
1154  <subtraction name="SteelSupport">
1155  <first ref="SteelSupportBlock"/>
1156  <second ref="FoamPadBlock"/>
1157  <positionref ref="posCenter"/>
1158  </subtraction>
1159 
1160  <box name="DetEnclosure" lunit="cm"
1161  x="$DetEncX"
1162  y="$DetEncY"
1163  z="$DetEncZ"/>
1164 
1165 </solids>
1166 EOF
1167 
1168 
1169 # Detector enclosure structure
1170  print ENCL <<EOF;
1171 <structure>
1172  <volume name="volFoamPadding">
1173  <materialref ref="fibrous_glass"/>
1174  <solidref ref="FoamPadding"/>
1175  </volume>
1176 
1177  <volume name="volSteelSupport">
1178  <materialref ref="AirSteelMixture"/>
1179  <solidref ref="SteelSupport"/>
1180  </volume>
1181 
1182  <volume name="volDetEnclosure">
1183  <materialref ref="Air"/>
1184  <solidref ref="DetEnclosure"/>
1185 
1186  <physvol>
1187  <volumeref ref="volFoamPadding"/>
1188  <positionref ref="posCryoInDetEnc"/>
1189  </physvol>
1190  <physvol>
1191  <volumeref ref="volSteelSupport"/>
1192  <positionref ref="posCryoInDetEnc"/>
1193  </physvol>
1194  <physvol>
1195  <volumeref ref="volCryostat"/>
1196  <positionref ref="posCryoInDetEnc"/>
1197  </physvol>
1198 EOF
1199 
1200 
1201 print ENCL <<EOF;
1202  </volume>
1203 EOF
1204 
1205 print ENCL <<EOF;
1206 </structure>
1207 </gdml>
1208 EOF
1209 
1210 close(ENCL);
1211 }
1212 
1213 
1214 
1215 #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1216 #+++++++++++++++++++++++++++++++++++++++ gen_World +++++++++++++++++++++++++++++++++++++++
1217 #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1218 
1219 sub gen_World()
1220 {
1221 
1222 # Create the WORLD fragment file name,
1223 # add file to list of output GDML fragments,
1224 # and open it
1225  $WORLD = $basename."_World" . $suffix . ".gdml";
1226  push (@gdmlFiles, $WORLD);
1227  $WORLD = ">" . $WORLD;
1228  open(WORLD) or die("Could not open file $WORLD for writing");
1229 
1230 
1231 # The standard XML prefix and starting the gdml
1232  print WORLD <<EOF;
1233 <?xml version='1.0'?>
1234 <gdml>
1235 EOF
1236 
1237 
1238 # All the World solids.
1239 print WORLD <<EOF;
1240 <solids>
1241  <box name="World" lunit="cm"
1242  x="@{[$DetEncX+2*$RockThickness]}"
1243  y="@{[$DetEncY+2*$RockThickness]}"
1244  z="@{[$DetEncZ+2*$RockThickness]}"/>
1245 </solids>
1246 EOF
1247 
1248 # World structure
1249 print WORLD <<EOF;
1250 <structure>
1251  <volume name="volWorld" >
1252  <materialref ref="DUSEL_Rock"/>
1253  <solidref ref="World"/>
1254 
1255  <physvol>
1256  <volumeref ref="volDetEnclosure"/>
1257  <position name="posDetEnclosure" unit="cm" x="$OriginXSet" y="$OriginYSet" z="$OriginZSet"/>
1258  </physvol>
1259 
1260  </volume>
1261 </structure>
1262 </gdml>
1263 EOF
1264 
1265 # make_gdml.pl will take care of <setup/>
1266 
1267 close(WORLD);
1268 }
1269 
1270 
1271 
1272 #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1273 #++++++++++++++++++++++++++++++++++++ write_fragments ++++++++++++++++++++++++++++++++++++
1274 #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1275 
1276 sub write_fragments()
1277 {
1278  # This subroutine creates an XML file that summarizes the the subfiles output
1279  # by the other sub routines - it is the input file for make_gdml.pl which will
1280  # give the final desired GDML file. Specify its name with the output option.
1281  # (you can change the name when running make_gdml)
1282 
1283  # This code is taken straigh from the similar MicroBooNE generate script, Thank you.
1284 
1285  if ( ! defined $output )
1286  {
1287  $output = "-"; # write to STDOUT
1288  }
1289 
1290  # Set up the output file.
1291  $OUTPUT = ">" . $output;
1292  open(OUTPUT) or die("Could not open file $OUTPUT");
1293 
1294  print OUTPUT <<EOF;
1295 <?xml version='1.0'?>
1296 
1297 <!-- Input to Geometry/gdml/make_gdml.pl; define the GDML fragments
1298  that will be zipped together to create a detector description.
1299  -->
1300 
1301 <config>
1302 
1303  <constantfiles>
1304 
1305  <!-- These files contain GDML <constant></constant>
1306  blocks. They are read in separately, so they can be
1307  interpreted into the remaining GDML. See make_gdml.pl for
1308  more information.
1309  -->
1310 
1311 EOF
1312 
1313  foreach $filename (@defFiles)
1314  {
1315  print OUTPUT <<EOF;
1316  <filename> $filename </filename>
1317 EOF
1318  }
1319 
1320  print OUTPUT <<EOF;
1321 
1322  </constantfiles>
1323 
1324  <gdmlfiles>
1325 
1326  <!-- The GDML file fragments to be zipped together. -->
1327 
1328 EOF
1329 
1330  foreach $filename (@gdmlFiles)
1331  {
1332  print OUTPUT <<EOF;
1333  <filename> $filename </filename>
1334 EOF
1335  }
1336 
1337  print OUTPUT <<EOF;
1338 
1339  </gdmlfiles>
1340 
1341 </config>
1342 EOF
1343 
1344  close(OUTPUT);
1345 }
1346 
1347 
1348 print "Some of the principal parameters for this TPC geometry (unit cm unless noted otherwise)\n";
1349 print " CRM active area : $widthCRM_active x $lengthCRM_active\n";
1350 print " CRM total area : $widthCRM x $lengthCRM\n";
1351 print " Wire pitch in U, V, Z : $wirePitchU, $wirePitchV, $wirePitchZ\n";
1352 print " TPC active volume : $driftTPCActive x $widthTPCActive x $lengthTPCActive\n";
1353 print " Argon volume : ($Argon_x, $Argon_y, $Argon_z) \n";
1354 print " Argon buffer : ($xLArBuffer, $yLArBuffer, $zLArBuffer) \n";
1355 print " Detector enclosure : $DetEncX x $DetEncY x $DetEncZ\n";
1356 print " TPC Origin : ($OriginXSet, $OriginYSet, $OriginZSet) \n";
1357 print " Field Cage : $FieldCage_switch \n";
1358 print " Cathode : $Cathode_switch \n";
1359 print " Workspace : $workspace \n";
1360 print " Wires : $wires_on \n";
1361 
1362 # run the sub routines that generate the fragments
1363 if ( $FieldCage_switch eq "on" ) { gen_FieldCage(); }
1364 #if ( $Cathode_switch eq "on" ) { gen_Cathode(); } #Cathode for now has the same geometry as the Ground Grid
1365 
1366 gen_Extend(); # generates the GDML color extension for the refactored geometry
1367 gen_Define(); # generates definitions at beginning of GDML
1368 gen_Materials(); # generates materials to be used
1369 gen_TPC(); # generate TPC for a given unit CRM
1370 gen_Cryostat(); #
1371 gen_Enclosure(); #
1372 gen_World(); # places the enclosure among DUSEL Rock
1373 write_fragments(); # writes the XML input for make_gdml.pl
1374  # which zips together the final GDML
1375 print "--- done\n\n\n";
1376 exit;