From 0fb63eacf8f584f8a7c6fab35589900b680f998d Mon Sep 17 00:00:00 2001
From: Christopher Dilks <dilks@jlab.org>
Date: Fri, 6 Aug 2021 19:49:08 +0000
Subject: [PATCH] Resolve "Detailed Dual RICH implementation"
---
athena.xml | 2 +-
compact/definitions.xml | 4 +-
compact/display.xml | 9 +
compact/drich.xml | 180 ++++++++++++++++
compact/materials.xml | 21 +-
compact/optical_materials.xml | 383 ++++++++++++++++++++++++++++++++
src/DRich.cpp | 395 ++++++++++++++++++++++++++++++++++
7 files changed, 989 insertions(+), 5 deletions(-)
create mode 100644 compact/drich.xml
create mode 100644 src/DRich.cpp
diff --git a/athena.xml b/athena.xml
index 89f8c620..8592b4af 100644
--- a/athena.xml
+++ b/athena.xml
@@ -164,7 +164,7 @@
<!--include ref="compact/dirc.xml"/-->
<!--include ref="compact/mrich.xml"/-->
<include ref="compact/forward_trd.xml"/>
- <include ref="compact/gaseous_rich.xml"/>
+ <include ref="compact/drich.xml"/>
<documentation level="10">
## Central calorimetry
diff --git a/compact/definitions.xml b/compact/definitions.xml
index f44c6e99..dc7901d3 100644
--- a/compact/definitions.xml
+++ b/compact/definitions.xml
@@ -324,7 +324,7 @@ Examples:
### PID Detector Region Parameters
</documentation>
- <constant name="ForwardRICH_length" value="1.1*m"/>
+ <constant name="ForwardRICH_length" value="180.0*cm"/>
<constant name="ForwardTRD_length" value="10.0*cm"/>
<constant name="ForwardTOF_length" value="3.0*cm"/>
@@ -353,7 +353,7 @@ Examples:
Generic tracking space allocations
</documentation>
- <constant name="ForwardTracking_length" value="0.0*cm"/>
+ <constant name="ForwardTracking_length" value="25.0*cm"/>
<documentation>
`BackwardTracking_length` and `ForwardTracking_length` compensate for the asymmetry of the setup
</documentation>
diff --git a/compact/display.xml b/compact/display.xml
index d20038be..2aeb5571 100644
--- a/compact/display.xml
+++ b/compact/display.xml
@@ -77,6 +77,15 @@
<vis name="ci_GEMVis" r= "0.8" g="0.4" b="0.3" alpha="0.8" showDaughters="true" visible="true"/>
<vis name="ci_HCALVis" r= "0.6" g="0" b="0.6" alpha="1.0" showDaughters="true" visible="true"/>
+ <vis name="DRICH_vessel_vis" alpha="1.0" r="1.0" g="1.0" b="1.0" showDaughters="true" visible="true" />
+ <vis name="DRICH_gas_vis" alpha="1.0" r="1.0" g="0.0" b="0.0" showDaughters="true" visible="true" />
+ <vis name="DRICH_aerogel_vis" alpha="1.0" r="0.0" g="1.0" b="1.0" showDaughters="true" visible="true" />
+ <vis name="DRICH_filter_vis" alpha="1.0" r="1.0" g="1.0" b="0.0" showDaughters="true" visible="true" />
+ <vis name="DRICH_mirror_vis" alpha="1.0" r="0.5" g="0.5" b="0.5" showDaughters="true" visible="true" />
+ <vis name="DRICH_sensor_vis" alpha="1.0" r="0.0" g="1.0" b="0.0" showDaughters="true" visible="true" />
+
+ </display>
+
<comment>
Deprecated colors.
vis name="GreenVis" alpha="1.0" r= "0.0" g="1.0" b="0.0" showDaughters="true" visible="true"/
diff --git a/compact/drich.xml b/compact/drich.xml
new file mode 100644
index 00000000..6a700900
--- /dev/null
+++ b/compact/drich.xml
@@ -0,0 +1,180 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<lccdd>
+
+ <define>
+ <!-- TODO [low priority]: some of these, viz. radii, could be parameterized
+ with respect to other variables; for now they are hard coded in case
+ other detectors' parameters are changing -->
+ <!-- parameters for re-scaling fun4all design to ATHENA -->
+ <constant name="DRICH_scale" value="0.963"/> <!-- overall scale factor from fun4all to ATHENA -->
+ <constant name="DRICH_f4a_length" value="161.0*cm"/> <!-- z-length of fun4all design -->
+ <!-- vessel (=snout+tank) geometry -->
+ <constant name="DRICH_zmin" value="BarrelTracking_length/2.0 + ForwardTracking_length "/> <!-- vessel front -->
+ <constant name="DRICH_Length" value="ForwardRICH_length"/> <!-- overall vessel length (including snout) -->
+ <constant name="DRICH_rmin0" value="ForwardPID_rmin1"/> <!-- bore radius at dRICh vessel frontplane -->
+ <constant name="DRICH_rmin1" value="19.0*cm"/> <!-- bore radius at dRICh vessel backplane -->
+ <constant name="DRICH_wall_thickness" value="0.5*cm"/> <!-- thickness of radial walls -->
+ <constant name="DRICH_window_thickness" value="0.1*cm"/> <!-- thickness of entrance and exit walls -->
+ <!-- snout geometry: cone with front radius rmax0 and back radius of rmax1 -->
+ <constant name="DRICH_SnoutLength" value="50.0*cm"/>
+ <constant name="DRICH_rmax0" value="110.0*cm"/>
+ <constant name="DRICH_rmax1" value="125.0*cm"/>
+ <!-- tank geometry: cylinder, holding the majority of detector components -->
+ <constant name="DRICH_rmax2" value="200*cm"/> <!-- cylinder radius; 20 cm gap between dRICh and HCalBarrel -->
+ <!-- additional parameters -->
+ <constant name="DRICH_aerogel_thickness" value="4.0*cm"/> <!-- aerogel thickness -->
+ </define>
+
+ <detectors>
+
+ <detector
+ id="ForwardRICH_ID"
+ name="DRICH"
+ type="athena_DRICH"
+ readout="DRICHHits"
+ gas="C2F6_DRICH"
+ material="Aluminum"
+ vis_vessel="DRICH_vessel_vis"
+ vis_gas="DRICH_gas_vis"
+ >
+
+ <!-- envelope dimensions (see above)
+ - `wall_thickness`: thickness of radial walls
+ - `window_thickness`: thickness of entrance and exit disks
+ -->
+ <dimensions
+ z0="DRICH_zmin"
+ length="DRICH_Length"
+ snout_length="DRICH_SnoutLength"
+ rmin0="DRICH_rmin0"
+ rmin1="DRICH_rmin1"
+ rmax0="DRICH_rmax0"
+ rmax1="DRICH_rmax1"
+ rmax2="DRICH_rmax2"
+ nsectors="6"
+ wall_thickness="DRICH_wall_thickness"
+ window_thickness="DRICH_window_thickness"
+ />
+
+ <!-- radiator defined in a wedge of azimuthal space
+ - `phiw` is phi width of wedge
+ - `thickness` defined separately for aerogel and filter
+ - `frontplane` is the front of the aerogel, w.r.t. front plane of vessel envelope
+ - `pitch` controls the angle of the radiator (0=vertical)
+ - filter is applied to backplane of aerogel
+ -->
+ <radiator
+ rmin="DRICH_rmin0 + DRICH_wall_thickness + 2.0*cm"
+ rmax="DRICH_rmax0 - DRICH_wall_thickness - 2.0*cm"
+ phiw="56*degree"
+ frontplane="DRICH_window_thickness + 0.5*DRICH_aerogel_thickness"
+ pitch="0*degree"
+ >
+ <aerogel
+ material="Aerogel_DRICH"
+ vis="DRICH_aerogel_vis"
+ thickness="DRICH_aerogel_thickness"
+ />
+ <filter
+ material="Acrylic_DRICH"
+ vis="DRICH_filter_vis"
+ thickness="0.3*mm"
+ />
+ </radiator>
+
+
+ <!-- spherical mirror is part of a sphere
+ - `rmin` and `rmax` provide polar angle boundaries
+ - `phiw` is the azimuthal width
+ - `radius` is the radius of the sphere
+ - `centerx` is the transverse position of the center
+ of the sphere, for the sector on the +x axis
+ - the back of the mirror will pass through `backplane`
+ - set `debug` to 1 so draw reference sphere instead, view with y-clipping
+ -->
+ <mirror
+ material="Acrylic_DRICH"
+ surface="MirrorSurface_DRICH"
+ vis="DRICH_mirror_vis"
+ backplane="DRICH_Length-2.0*cm"
+ thickness="0.2*cm"
+ radius="290*DRICH_scale*cm"
+ centerx="145*DRICH_scale*cm"
+ rmin="DRICH_rmin1 + DRICH_wall_thickness + 0.0*cm"
+ rmax="DRICH_rmax2 - DRICH_wall_thickness - 2.0*cm"
+ phiw="54*degree"
+ debug="0"
+ />
+
+ <sensors>
+ <!-- geometry for a single square sensor
+ - based on Hamamatsu H13700 MAPMT
+ (https://www.hamamatsu.com/us/en/product/type/H13700/index.html)
+ - N.B. not ideal for a magnetic field, SiPM matrix would be better
+ - effective area: 48.5x48.5mm
+ - enclosure size: 52x52mm
+ - 16x16 channel matrix (see readout segmentation below)
+ - pixel size: 3x3mm
+ - `side` is the side length of the square
+ - `thickness` is the depth of the sensor
+ - `gap` provides room between the squares, to help
+ prevent them from overlapping
+ - the value of `side` will determine how many sensors there are,
+ since the sensor placement algorithm will try to place as many
+ as it can in the specified patch below
+ -->
+ <module
+ material="Silicon"
+ surface="SensorSurface_DRICH"
+ vis="DRICH_sensor_vis"
+ side="48*mm"
+ thickness="35*mm"
+ gap="0.5*(52-48)*mm + 2*mm"
+ />
+ <!-- sensors will be tiled on this sphere
+ - `center{x,y,z} is defined for sector on +x axis, defined w.r.t. snout frontplane,
+ and `radius` is the sphere radius; the first term of each of these comes from
+ the fun4all design
+ - these attributes were determined from a spherical fit to the
+ sensor placement in the fun4all port
+ - set `debug` to 1 so draw reference sphere instead, view with y-clipping
+ -->
+ <sphere
+ radius="159.76*DRICH_scale*cm"
+ centerx="144.91*DRICH_scale*cm"
+ centery="0*DRICH_scale*cm"
+ centerz="-197.25*DRICH_scale*cm + DRICH_Length - 0.5*DRICH_scale*DRICH_f4a_length"
+ debug="0"
+ />
+ <!-- sensors will be limited to a patch of the sphere
+ - `thetamin` and `thetamax` define pseudorapidity coverage
+ - `widthfactor` controls the azimuthal coverage, where lower=wider
+ - `taper` defines half the angle between the azimuthal boundaries
+ - the size of the sensor controls how many sensors are placed
+ -->
+ <sphericalpatch
+ thetamin="-10*degree"
+ thetamax="22*degree"
+ widthfactor="1.8"
+ taper="56*degree"
+ />
+ </sensors>
+
+ </detector>
+ </detectors>
+
+ <readouts>
+ <readout name="DRICHHits">
+ <segmentation type="CartesianGridXY" grid_size_x="3*mm" grid_size_y="3*mm"/>
+ <!-- cellID: 64bits
+ - bits 1-8: dRICh ID
+ - bits 9-16: sector number
+ - bits 17-32: photosensor number
+ - bits 33-48: x pixel (signed)
+ - bits 49-64: y pixel (signed)
+ -->
+ <id>system:8,sector:8,module:16,x:32:-16,y:-16</id>
+ </readout>
+ </readouts>
+
+</lccdd>
diff --git a/compact/materials.xml b/compact/materials.xml
index e305c5fa..a1ab3a2a 100644
--- a/compact/materials.xml
+++ b/compact/materials.xml
@@ -24,7 +24,7 @@
<composite n="15" ref="C"/>
<composite n="7" ref="O"/>
</material>
- <material name="Quartz">
+ <material name="Quartz"> <!-- different density than `SiliconDioxide` -->
<D type="density" value="2.2" unit="g/cm3"/>
<composite n="1" ref="Si"/>
<composite n="2" ref="O"/>
@@ -45,11 +45,16 @@
<fraction n="0.998" ref="Fe"/>
<fraction n=".002" ref="C"/>
</material>
- <material name="SiliconOxide">
+ <material name="SiliconOxide"> <!-- different density than `SiliconDioxide` -->
<D type="density" value="2.65" unit="g/cm3"/>
<composite n="1" ref="Si"/>
<composite n="2" ref="O"/>
</material>
+ <material name="SiliconDioxide"> <!-- density from `G4_SILICON_DIOXIDE` (NIST DB) -->
+ <D type="density" value="2.32" unit="g/cm3"/>
+ <composite n="1" ref="Si"/>
+ <composite n="2" ref="O"/>
+ </material>
<material name="BoronOxide">
<D type="density" value="2.46" unit="g/cm3"/>
<composite n="2" ref="B"/>
@@ -258,6 +263,18 @@
<D type="density" value="2.21" unit="g/cm3"/>
<composite n="1" ref="C"/>
</material>
+ <material name="PolyvinylAcetate">
+ <D type="density" value="1.19" unit="g/cm3"/>
+ <composite n="4" ref="C"/>
+ <composite n="6" ref="H"/>
+ <composite n="2" ref="O"/>
+ </material>
+ <material name="Plexiglass">
+ <D type="density" value="1.19" unit="g/cm3"/>
+ <composite n="5" ref="C"/>
+ <composite n="8" ref="H"/>
+ <composite n="2" ref="O"/>
+ </material>
<material name="StainlessSteel">
<D type="density" value="8.3" unit="g / cm3"/>
<fraction n="0.74" ref="Fe"/>
diff --git a/compact/optical_materials.xml b/compact/optical_materials.xml
index 75daac07..c510b815 100644
--- a/compact/optical_materials.xml
+++ b/compact/optical_materials.xml
@@ -47,6 +47,298 @@
1240*eV/600 1.49
1240*eV/400 1.49
"/>
+ <!-- BEGIN dRICh material properties
+ - dumped from fun4all implementation
+ - see https://github.com/cisbani/dRICh/blob/main/share/source/g4dRIChOptics.hh
+ -->
+ <!-- dRICh gas -->
+ <matrix name="RINDEX__C2F6_DRICH" coldim="2" values="
+ 1.7712*eV 1.00075
+ 1.92389*eV 1.00075
+ 2.10539*eV 1.00075
+ 2.3247*eV 1.00075
+ 2.59502*eV 1.00076
+ 2.93647*eV 1.00076
+ 3.38139*eV 1.00077
+ 3.98521*eV 1.00078
+ 4.85156*eV 1.0008
+ 6.19921*eV 1.00084
+ "/>
+ <matrix name="ABSLENGTH__C2F6_DRICH" coldim="2" values="
+ 1.7712*eV 1000*mm
+ 1.92389*eV 1000*mm
+ 2.10539*eV 1000*mm
+ 2.3247*eV 1000*mm
+ 2.59502*eV 1000*mm
+ 2.93647*eV 1000*mm
+ 3.38139*eV 1000*mm
+ 3.98521*eV 1000*mm
+ 4.85156*eV 1000*mm
+ 6.19921*eV 1000*mm
+ "/>
+ <!-- dRICh aerogel -->
+ <matrix name="RINDEX__Aerogel_DRICH" coldim="2" values="
+ 1.87855*eV 1.01638
+ 1.96673*eV 1.01642
+ 2.0549*eV 1.01647
+ 2.14308*eV 1.01652
+ 2.23126*eV 1.01657
+ 2.31943*eV 1.01662
+ 2.40761*eV 1.01667
+ 2.49579*eV 1.01673
+ 2.58396*eV 1.01679
+ 2.67214*eV 1.01685
+ 2.76032*eV 1.01691
+ 2.84849*eV 1.01698
+ 2.93667*eV 1.01705
+ 3.02485*eV 1.01712
+ 3.11302*eV 1.01719
+ 3.2012*eV 1.01727
+ 3.28938*eV 1.01734
+ 3.37755*eV 1.01742
+ 3.46573*eV 1.01751
+ 3.55391*eV 1.01759
+ 3.64208*eV 1.01768
+ 3.73026*eV 1.01777
+ 3.81844*eV 1.01787
+ 3.90661*eV 1.01796
+ 3.99479*eV 1.01806
+ 4.08297*eV 1.01816
+ 4.17114*eV 1.01827
+ 4.25932*eV 1.01838
+ 4.3475*eV 1.01849
+ 4.43567*eV 1.0186
+ 4.52385*eV 1.01872
+ 4.61203*eV 1.01884
+ 4.7002*eV 1.01897
+ 4.78838*eV 1.01909
+ 4.87656*eV 1.01922
+ 4.96473*eV 1.01936
+ 5.05291*eV 1.0195
+ 5.14109*eV 1.01964
+ 5.22927*eV 1.01979
+ 5.31744*eV 1.01994
+ 5.40562*eV 1.02009
+ 5.4938*eV 1.02025
+ 5.58197*eV 1.02041
+ 5.67015*eV 1.02057
+ 5.75833*eV 1.02074
+ 5.8465*eV 1.02092
+ 5.93468*eV 1.0211
+ 6.02286*eV 1.02128
+ 6.11103*eV 1.02147
+ 6.19921*eV 1.02167
+ "/>
+ <matrix name="ABSLENGTH__Aerogel_DRICH" coldim="2" values="
+ 1.87855*eV 154*mm
+ 1.96673*eV 156.17*mm
+ 2.0549*eV 158.154*mm
+ 2.14308*eV 159.974*mm
+ 2.23126*eV 161.651*mm
+ 2.31943*eV 163.2*mm
+ 2.40761*eV 164.636*mm
+ 2.49579*eV 165.97*mm
+ 2.58396*eV 167.213*mm
+ 2.67214*eV 168.374*mm
+ 2.76032*eV 169.461*mm
+ 2.84849*eV 170.481*mm
+ 2.93667*eV 171.439*mm
+ 3.02485*eV 172.342*mm
+ 3.11302*eV 173.193*mm
+ 3.2012*eV 173.998*mm
+ 3.28938*eV 174.759*mm
+ 3.37755*eV 175.481*mm
+ 3.46573*eV 176.165*mm
+ 3.55391*eV 176.817*mm
+ 3.64208*eV 162.708*mm
+ 3.73026*eV 140.953*mm
+ 3.81844*eV 122.516*mm
+ 3.90661*eV 106.833*mm
+ 3.99479*eV 93.4428*mm
+ 4.08297*eV 81.9694*mm
+ 4.17114*eV 72.1072*mm
+ 4.25932*eV 63.6011*mm
+ 4.3475*eV 56.2434*mm
+ 4.43567*eV 49.8599*mm
+ 4.52385*eV 44.3054*mm
+ 4.61203*eV 39.4601*mm
+ 4.7002*eV 35.2211*mm
+ 4.78838*eV 31.5049*mm
+ 4.87656*eV 28.2374*mm
+ 4.96473*eV 25.359*mm
+ 5.05291*eV 22.8166*mm
+ 5.14109*eV 20.5674*mm
+ 5.22927*eV 18.5724*mm
+ 5.31744*eV 16.7992*mm
+ 5.40562*eV 15.2205*mm
+ 5.4938*eV 13.8125*mm
+ 5.58197*eV 12.5541*mm
+ 5.67015*eV 11.4277*mm
+ 5.75833*eV 10.4166*mm
+ 5.8465*eV 9.50928*mm
+ 5.93468*eV 8.69176*mm
+ 6.02286*eV 7.95608*mm
+ 6.11103*eV 7.29168*mm
+ 6.19921*eV 6.69064*mm
+ "/>
+ <matrix name="RAYLEIGH__Aerogel_DRICH" coldim="2" values="
+ 1.87855*eV 309.218*mm
+ 1.96673*eV 257.383*mm
+ 2.0549*eV 215.968*mm
+ 2.14308*eV 182.559*mm
+ 2.23126*eV 155.367*mm
+ 2.31943*eV 133.053*mm
+ 2.40761*eV 114.606*mm
+ 2.49579*eV 99.2482*mm
+ 2.58396*eV 86.3795*mm
+ 2.67214*eV 75.5291*mm
+ 2.76032*eV 66.3313*mm
+ 2.84849*eV 58.4918*mm
+ 2.93667*eV 51.777*mm
+ 3.02485*eV 45.998*mm
+ 3.11302*eV 41.0045*mm
+ 3.2012*eV 36.6696*mm
+ 3.28938*eV 32.8931*mm
+ 3.37755*eV 29.5904*mm
+ 3.46573*eV 26.6917*mm
+ 3.55391*eV 24.1402*mm
+ 3.64208*eV 21.8856*mm
+ 3.73026*eV 19.8884*mm
+ 3.81844*eV 18.1144*mm
+ 3.90661*eV 16.533*mm
+ 3.99479*eV 15.1206*mm
+ 4.08297*eV 13.856*mm
+ 4.17114*eV 12.7208*mm
+ 4.25932*eV 11.7005*mm
+ 4.3475*eV 10.7791*mm
+ 4.43567*eV 9.94752*mm
+ 4.52385*eV 9.1938*mm
+ 4.61203*eV 8.51136*mm
+ 4.7002*eV 7.88964*mm
+ 4.78838*eV 7.32468*mm
+ 4.87656*eV 6.80988*mm
+ 4.96473*eV 6.33864*mm
+ 5.05291*eV 5.907*mm
+ 5.14109*eV 5.51232*mm
+ 5.22927*eV 5.14932*mm
+ 5.31744*eV 4.81668*mm
+ 5.40562*eV 4.51044*mm
+ 5.4938*eV 4.22796*mm
+ 5.58197*eV 3.9666*mm
+ 5.67015*eV 3.72504*mm
+ 5.75833*eV 3.50196*mm
+ 5.8465*eV 3.29604*mm
+ 5.93468*eV 3.10464*mm
+ 6.02286*eV 2.92644*mm
+ 6.11103*eV 2.76144*mm
+ 6.19921*eV 2.607*mm
+ "/>
+ <!-- dRICh acrylic -->
+ <matrix name="RINDEX__Acrylic_DRICH" coldim="2" values="
+ 4.13281*eV 1.5017
+ 4.22099*eV 1.5017
+ 4.30916*eV 1.5017
+ 4.39734*eV 1.5017
+ 4.48552*eV 1.5017
+ 4.57369*eV 1.5017
+ 4.66187*eV 1.5017
+ 4.75005*eV 1.5017
+ 4.83822*eV 1.5017
+ 4.9264*eV 1.5017
+ 5.01458*eV 1.5017
+ 5.10275*eV 1.5017
+ 5.19093*eV 1.5017
+ 5.27911*eV 1.5017
+ 5.36728*eV 1.5017
+ 5.45546*eV 1.5017
+ 5.54364*eV 1.5017
+ 5.63181*eV 1.5017
+ 5.71999*eV 1.5017
+ 5.80817*eV 1.5017
+ 5.89634*eV 1.5017
+ 5.98452*eV 1.5017
+ 6.0727*eV 1.5017
+ 6.16087*eV 1.5017
+ 6.24905*eV 1.5017
+ 6.33723*eV 1.5017
+ 6.4254*eV 1.5017
+ 6.51358*eV 1.5017
+ 6.60176*eV 1.5017
+ 6.68993*eV 1.5017
+ 6.77811*eV 1.5017
+ 6.86629*eV 1.5017
+ 6.95446*eV 1.5017
+ 7.04264*eV 1.5017
+ 7.13082*eV 1.5017
+ 7.21899*eV 1.5017
+ 7.30717*eV 1.5017
+ 7.39535*eV 1.5017
+ 7.48353*eV 1.5017
+ 7.5717*eV 1.5017
+ 7.65988*eV 1.5017
+ 7.74806*eV 1.5017
+ 7.83623*eV 1.5017
+ 7.92441*eV 1.5017
+ 8.01259*eV 1.5017
+ 8.10076*eV 1.5017
+ 8.18894*eV 1.5017
+ 8.27712*eV 1.5017
+ 8.36529*eV 1.5017
+ 8.45347*eV 1.5017
+ "/>
+ <matrix name="ABSLENGTH__Acrylic_DRICH" coldim="2" values="
+ 4.13281*eV 82.0704*mm
+ 4.22099*eV 36.9138*mm
+ 4.30916*eV 13.3325*mm
+ 4.39734*eV 5.03627*mm
+ 4.48552*eV 2.3393*mm
+ 4.57369*eV 1.36177*mm
+ 4.66187*eV 0.933192*mm
+ 4.75005*eV 0.708268*mm
+ 4.83822*eV 0.573082*mm
+ 4.9264*eV 0.483641*mm
+ 5.01458*eV 0.420282*mm
+ 5.10275*eV 0.373102*mm
+ 5.19093*eV 0.33662*mm
+ 5.27911*eV 0.307572*mm
+ 5.36728*eV 0.283902*mm
+ 5.45546*eV 0.264235*mm
+ 5.54364*eV 0.247641*mm
+ 5.63181*eV 0.233453*mm
+ 5.71999*eV 0.221177*mm
+ 5.80817*eV 0.210456*mm
+ 5.89634*eV 0.201012*mm
+ 5.98452*eV 0.192627*mm
+ 6.0727*eV 0.185134*mm
+ 6.16087*eV 0.178399*mm
+ 6.24905*eV 0.172309*mm
+ 6.33723*eV 0.166779*mm
+ 6.4254*eV 0.166779*mm
+ 6.51358*eV 0.166779*mm
+ 6.60176*eV 0.166779*mm
+ 6.68993*eV 0.166779*mm
+ 6.77811*eV 0.166779*mm
+ 6.86629*eV 0.166779*mm
+ 6.95446*eV 0.166779*mm
+ 7.04264*eV 0.166779*mm
+ 7.13082*eV 0.166779*mm
+ 7.21899*eV 0.166779*mm
+ 7.30717*eV 0.166779*mm
+ 7.39535*eV 0.166779*mm
+ 7.48353*eV 0.166779*mm
+ 7.5717*eV 0.166779*mm
+ 7.65988*eV 0.166779*mm
+ 7.74806*eV 0.166779*mm
+ 7.83623*eV 0.166779*mm
+ 7.92441*eV 0.166779*mm
+ 8.01259*eV 0.166779*mm
+ 8.10076*eV 0.166779*mm
+ 8.18894*eV 0.166779*mm
+ 8.27712*eV 0.166779*mm
+ 8.36529*eV 0.166779*mm
+ 8.45347*eV 0.166779*mm
+ "/>
+ <!-- END dRICh optical properties -->
</properties>
<materials>
<material name="AirOptical">
@@ -86,6 +378,32 @@
<composite n="8" ref="H"/>
<property name="RINDEX" ref="RINDEX__Acrylic"/>
</material>
+ <!-- BEGIN dRICh material definitions -->
+ <material name="C2F6_DRICH">
+ <D type="density" value="0.005734" unit="g/cm3"/>
+ <composite n="2" ref="C"/>
+ <composite n="6" ref="F"/>
+ <property name="RINDEX" ref="RINDEX__C2F6_DRICH"/>
+ <property name="ABSLENGTH" ref="ABSLENGTH__C2F6_DRICH"/>
+ </material>
+ <material name="Aerogel_DRICH">
+ <D type="density" value="0.1" unit="g/cm3"/>
+ <comment> n_air = [dens(Si02)-dens(aerogel)] / [dens(Si02)-dens(Air) ] </comment>
+ <fraction n=" (2.32-0.1) / (2.32-0.0012)" ref="Air"/>
+ <fraction n="1 - (2.32-0.1) / (2.32-0.0012)" ref="SiliconDioxide"/>
+ <property name="RINDEX" ref="RINDEX__Aerogel_DRICH"/>
+ <property name="ABSLENGTH" ref="ABSLENGTH__Aerogel_DRICH"/>
+ <property name="RAYLEIGH" ref="RAYLEIGH__Aerogel_DRICH"/>
+ </material>
+ <material name="Acrylic_DRICH">
+ <D type="density" value="1.19" unit="g/cm3"/>
+ <comment> TO BE IMPROVED </comment>
+ <fraction n="0.99" ref="Plexiglass"/>
+ <fraction n="0.01" ref="PolyvinylAcetate"/>
+ <property name="RINDEX" ref="RINDEX__Acrylic_DRICH"/>
+ <property name="ABSLENGTH" ref="ABSLENGTH__Acrylic_DRICH"/>
+ </material>
+ <!-- END dRICh material definitions -->
</materials>
<surfaces>
<opticalsurface finish="polished" model="glisur" name="MirrorOpticalSurface" type="dielectric_metal" value="0">
@@ -131,5 +449,70 @@
-->
</opticalsurface>
+ <!-- BEGIN dRICh surface definitions -->
+ <opticalsurface name="MirrorSurface_DRICH" model="unified" finish="polishedfrontpainted" type="dielectric_dielectric">
+ <property name="REFLECTIVITY" coldim="2" values="
+ 2.04358*eV 0.867812
+ 2.0664*eV 0.865156
+ 2.09046*eV 0.863906
+ 2.14023*eV 0.86375
+ 2.16601*eV 0.864062
+ 2.20587*eV 0.864531
+ 2.23327*eV 0.864375
+ 2.26137*eV 0.865625
+ 2.31972*eV 0.865313
+ 2.35005*eV 0.865
+ 2.38116*eV 0.864844
+ 2.41313*eV 0.863828
+ 2.44598*eV 0.863516
+ 2.47968*eV 0.863125
+ 2.53081*eV 0.862187
+ 2.58354*eV 0.861719
+ 2.6194*eV 0.861328
+ 2.69589*eV 0.861016
+ 2.73515*eV 0.861094
+ 2.79685*eV 0.861602
+ 2.86139*eV 0.862305
+ 2.95271*eV 0.86375
+ 3.04884*eV 0.865586
+ 3.12665*eV 0.867383
+ 3.2393*eV 0.870059
+ 3.39218*eV 0.874199
+ 3.52508*eV 0.878105
+ 3.66893*eV 0.88252
+ 3.82396*eV 0.887617
+ 3.99949*eV 0.893721
+ 4.13281*eV 0.898184
+ 4.27679*eV 0.902744
+ 4.48244*eV 0.907837
+ 4.65057*eV 0.9102
+ 4.89476*eV 0.909316
+ 5.02774*eV 0.906174
+ 5.16816*eV 0.900422
+ 5.31437*eV 0.891521
+ 5.63821*eV 0.859954
+ 5.90401*eV 0.820831
+ 6.19921*eV 0.762502
+ "/>
+ </opticalsurface>
+ <opticalsurface name="SensorSurface_DRICH" model="glisur" finish="polished" type="dielectric_metal">
+ <property name="EFFICIENCY" coldim="2" values="
+ 1*eV 1
+ 4*eV 1
+ 7*eV 1
+ "/>
+ <property name="IMAGINARYRINDEX" coldim="2" values="
+ 1*eV 1.69
+ 4*eV 1.69
+ 7*eV 1.69
+ "/>
+ <property name="REALRINDEX" coldim="2" values="
+ 1*eV 1.92
+ 4*eV 1.92
+ 7*eV 1.92
+ "/>
+ </opticalsurface>
+ <!-- END dRICh surface definitions -->
+
</surfaces>
</lccdd>
diff --git a/src/DRich.cpp b/src/DRich.cpp
new file mode 100644
index 00000000..5d9ebd83
--- /dev/null
+++ b/src/DRich.cpp
@@ -0,0 +1,395 @@
+//==========================================================================
+// dRICh: Dual Ring Imaging Cherenkov Detector
+//--------------------------------------------------------------------------
+//
+// Author: Christopher Dilks (Duke University)
+//
+// - Design Adapted from Standalone Fun4all and GEMC implementations
+// [ Evaristo Cisbani, Cristiano Fanelli, Alessio Del Dotto, et al. ]
+//
+//==========================================================================
+
+#include <XML/Helper.h>
+#include "TMath.h"
+#include "TString.h"
+#include "GeometryHelpers.h"
+#include "Math/Point2D.h"
+#include "DDRec/Surface.h"
+#include "DDRec/DetectorData.h"
+#include "DD4hep/OpticalSurfaces.h"
+#include "DD4hep/DetFactoryHelper.h"
+#include "DD4hep/Printout.h"
+
+using namespace std;
+using namespace dd4hep;
+using namespace dd4hep::rec;
+
+// create the detector
+static Ref_t createDetector(Detector& desc, xml::Handle_t handle, SensitiveDetector sens) {
+ xml::DetElement detElem = handle;
+ std::string detName = detElem.nameStr();
+ int detID = detElem.id();
+
+ DetElement det(detName, detID);
+ xml::Component dims = detElem.dimensions();
+ OpticalSurfaceManager surfMgr = desc.surfaceManager();
+
+ // attributes -----------------------------------------------------------
+ // TODO [low priority]: some attributes have default values, some do not;
+ // make this more consistent
+ // - vessel
+ double vesselZ0 = dims.attr<double>(_Unicode(z0));
+ double vesselLength = dims.attr<double>(_Unicode(length));
+ double vesselRmin0 = dims.attr<double>(_Unicode(rmin0));
+ double vesselRmin1 = dims.attr<double>(_Unicode(rmin1));
+ double vesselRmax0 = dims.attr<double>(_Unicode(rmax0));
+ double vesselRmax1 = dims.attr<double>(_Unicode(rmax1));
+ double vesselRmax2 = dims.attr<double>(_Unicode(rmax2));
+ double snoutLength = dims.attr<double>(_Unicode(snout_length));
+ int nSectors = getAttrOrDefault<int>(dims, _Unicode(nsectors), 6);
+ double wallThickness = getAttrOrDefault<double>(dims, _Unicode(wall_thickness), 0.5*cm);
+ double windowThickness = getAttrOrDefault<double>(dims, _Unicode(window_thickness), 0.1*cm);
+ auto vesselMat = desc.material(getAttrOrDefault(detElem, _Unicode(material), "Aluminum"));
+ auto gasvolMat = desc.material(getAttrOrDefault(detElem, _Unicode(gas), "AirOptical"));
+ auto vesselVis = desc.visAttributes(detElem.attr<std::string>(_Unicode(vis_vessel)));
+ auto gasvolVis = desc.visAttributes(detElem.attr<std::string>(_Unicode(vis_gas)));
+ // - radiator (applies to aerogel and filter)
+ auto radiatorElem = detElem.child(_Unicode(radiator));
+ double radiatorRmin = getAttrOrDefault<double>(radiatorElem, _Unicode(rmin), 10.*cm);
+ double radiatorRmax = getAttrOrDefault<double>(radiatorElem, _Unicode(rmax), 80.*cm);
+ double radiatorPhiw = getAttrOrDefault<double>(radiatorElem, _Unicode(phiw), 2*M_PI/nSectors);
+ double radiatorPitch = getAttrOrDefault<double>(radiatorElem, _Unicode(pitch), 0.*degree);
+ double radiatorFrontplane = getAttrOrDefault<double>(radiatorElem, _Unicode(frontplane), 2.5*cm);
+ // - aerogel
+ auto aerogelElem = radiatorElem.child(_Unicode(aerogel));
+ auto aerogelMat = desc.material(aerogelElem.attr<std::string>(_Unicode(material)));
+ auto aerogelVis = desc.visAttributes(aerogelElem.attr<std::string>(_Unicode(vis)));
+ double aerogelThickness = getAttrOrDefault<double>(aerogelElem, _Unicode(thickness), 2.5*cm);
+ // - filter
+ auto filterElem = radiatorElem.child(_Unicode(filter));
+ auto filterMat = desc.material(filterElem.attr<std::string>(_Unicode(material)));
+ auto filterVis = desc.visAttributes(filterElem.attr<std::string>(_Unicode(vis)));
+ double filterThickness = getAttrOrDefault<double>(filterElem, _Unicode(thickness), 2.5*cm);
+ // - mirror
+ auto mirrorElem = detElem.child(_Unicode(mirror));
+ auto mirrorMat = desc.material(mirrorElem.attr<std::string>(_Unicode(material)));
+ auto mirrorVis = desc.visAttributes(mirrorElem.attr<std::string>(_Unicode(vis)));
+ auto mirrorSurf = surfMgr.opticalSurface(getAttrOrDefault(mirrorElem, _Unicode(surface), "MirrorOpticalSurface"));
+ double mirrorBackplane = getAttrOrDefault<double>(mirrorElem, _Unicode(backplane), 240.*cm);
+ double mirrorThickness = getAttrOrDefault<double>(mirrorElem, _Unicode(thickness), 2.*mm);
+ double mirrorRadius = getAttrOrDefault<double>(mirrorElem, _Unicode(radius), 190*cm);
+ double mirrorCenterX = getAttrOrDefault<double>(mirrorElem, _Unicode(centerx), 95*cm);
+ double mirrorRmin = getAttrOrDefault<double>(mirrorElem, _Unicode(rmin), 10.*cm);
+ double mirrorRmax = getAttrOrDefault<double>(mirrorElem, _Unicode(rmax), 150.*cm);
+ double mirrorPhiw = getAttrOrDefault<double>(mirrorElem, _Unicode(phiw), 2*M_PI/nSectors);
+ int mirrorDebug = getAttrOrDefault<int>(mirrorElem, _Unicode(debug), 0);
+ // - sensor module
+ auto sensorElem = detElem.child(_Unicode(sensors)).child(_Unicode(module));
+ auto sensorMat = desc.material(sensorElem.attr<std::string>(_Unicode(material)));
+ auto sensorVis = desc.visAttributes(sensorElem.attr<std::string>(_Unicode(vis)));
+ auto sensorSurf = surfMgr.opticalSurface(getAttrOrDefault(sensorElem, _Unicode(surface), "MirrorOpticalSurface"));
+ double sensorSide = sensorElem.attr<double>(_Unicode(side));
+ double sensorGap = sensorElem.attr<double>(_Unicode(gap));
+ double sensorThickness = sensorElem.attr<double>(_Unicode(thickness));
+ // - sensor sphere
+ auto sensorSphElem = detElem.child(_Unicode(sensors)).child(_Unicode(sphere));
+ double sensorSphRadius = sensorSphElem.attr<double>(_Unicode(radius));
+ double sensorSphCenterX = sensorSphElem.attr<double>(_Unicode(centerx));
+ double sensorSphCenterY = sensorSphElem.attr<double>(_Unicode(centery));
+ double sensorSphCenterZ = sensorSphElem.attr<double>(_Unicode(centerz));
+ int sensorSphDebug = getAttrOrDefault<int>(sensorSphElem, _Unicode(debug), 0);
+ // - sensor sphere patch cuts
+ auto sensorSphPatchElem = detElem.child(_Unicode(sensors)).child(_Unicode(sphericalpatch));
+ double sensorSphPatchThetaMin = sensorSphPatchElem.attr<double>(_Unicode(thetamin));
+ double sensorSphPatchThetaMax = sensorSphPatchElem.attr<double>(_Unicode(thetamax));
+ double sensorSphPatchWidthFactor = sensorSphPatchElem.attr<double>(_Unicode(widthfactor));
+ double sensorSphPatchTaper = sensorSphPatchElem.attr<double>(_Unicode(taper));
+
+
+ // BUILD VESSEL ====================================================================
+ /* - `vessel`: aluminum enclosure, the mother volume of the dRICh
+ * - `gasvol`: gas volume, which fills `vessel`; all other volumes defined below
+ * are children of `gasvol`
+ * - the dRICh vessel geometry has two regions: the snout refers to the conic region
+ * in the front, housing the aerogel, while the tank refers to the cylindrical
+ * region, housing the rest of the detector components
+ */
+
+ // snout solids
+ double boreDelta = vesselRmin1 - vesselRmin0;
+ Cone vesselSnout(
+ snoutLength/2.0,
+ vesselRmin0,
+ vesselRmax0,
+ vesselRmin0 + boreDelta * snoutLength / vesselLength,
+ vesselRmax1
+ );
+ Cone gasvolSnout(
+ /* note: `gasvolSnout` extends a bit into the tank, so it touches `gasvolTank`
+ * - the extension distance is equal to the tank `windowThickness`, so the
+ * length of `gasvolSnout` == length of `vesselSnout`
+ * - the extension backplane radius is calculated using similar triangles
+ */
+ snoutLength/2.0,
+ vesselRmin0 + wallThickness,
+ vesselRmax0 - wallThickness,
+ vesselRmin0 + boreDelta * (snoutLength-windowThickness) / vesselLength + wallThickness,
+ vesselRmax1 - wallThickness + windowThickness * (vesselRmax1 - vesselRmax0) / snoutLength
+ );
+
+ // tank solids
+ Cone vesselTank(
+ (vesselLength - snoutLength)/2.0,
+ vesselSnout.rMin2(),
+ vesselRmax2,
+ vesselRmin1,
+ vesselRmax2
+ );
+ Cone gasvolTank(
+ (vesselLength - snoutLength)/2.0 - windowThickness,
+ gasvolSnout.rMin2(),
+ vesselRmax2 - wallThickness,
+ vesselRmin1 + wallThickness,
+ vesselRmax2 - wallThickness
+ );
+
+ // snout + tank solids
+ UnionSolid vesselSolid(
+ vesselTank,
+ vesselSnout,
+ Position(0., 0., -vesselLength/2.)
+ );
+ UnionSolid gasvolSolid(
+ gasvolTank,
+ gasvolSnout,
+ Position(0., 0., -vesselLength/2. + windowThickness)
+ );
+
+
+ // volumes
+ Volume vesselVol(detName, vesselSolid, vesselMat);
+ Volume gasvolVol(detName + "_gas", gasvolSolid, gasvolMat);
+ vesselVol.setVisAttributes(vesselVis);
+ gasvolVol.setVisAttributes(gasvolVis);
+
+ // reference position
+ auto snoutFront = Position(0., 0., -(vesselLength + snoutLength)/2.);
+
+
+ // sensitive detector type
+ sens.setType("photoncounter");
+
+
+ // SECTOR LOOP //////////////////////////////////
+ for(int isec=0; isec<nSectors; isec++) {
+
+ if(mirrorDebug*isec>0 || sensorSphDebug*isec>0) continue; // if debugging, draw only 1 sector
+ double sectorRotation = isec * 360/nSectors * degree; // sector rotation about z axis
+
+ // BUILD RADIATOR ====================================================================
+
+ // derived attributes
+ auto radiatorPos = Position(0., 0., radiatorFrontplane) + snoutFront;
+
+ // solid and volume: create aerogel and filter sectors
+ Tube aerogelSolid(radiatorRmin, radiatorRmax, aerogelThickness/2, -radiatorPhiw/2.0, radiatorPhiw/2.0);
+ Tube filterSolid( radiatorRmin, radiatorRmax, filterThickness/2, -radiatorPhiw/2.0, radiatorPhiw/2.0);
+ Volume aerogelVol("aerogel_v", aerogelSolid, aerogelMat);
+ Volume filterVol( "filter_v", filterSolid, filterMat);
+ aerogelVol.setVisAttributes(aerogelVis);
+ filterVol.setVisAttributes(filterVis);
+
+ // placement
+ auto aerogelPV = gasvolVol.placeVolume(aerogelVol,
+ RotationZ(sectorRotation) // rotate about beam axis to sector
+ * Translation3D(radiatorPos.x(), radiatorPos.y(), radiatorPos.z()) // re-center to snoutFront
+ * RotationY(radiatorPitch) // change polar angle to specified pitch
+ );
+ auto filterPV = gasvolVol.placeVolume(filterVol,
+ RotationZ(sectorRotation) // rotate about beam axis to sector
+ * Translation3D(radiatorPos.x(), radiatorPos.y(), radiatorPos.z()) // re-center to snoutFront
+ * RotationY(radiatorPitch) // change polar angle
+ * Translation3D(0., 0., (aerogelThickness+filterThickness)/2.) // move to aerogel backplane
+ );
+
+ // properties
+ // TODO [critical]: define skin properties for aerogel and filter
+ DetElement aerogelDE(det, Form("aerogel_de%d", isec), isec);
+ aerogelDE.setPlacement(aerogelPV);
+ //SkinSurface aerogelSkin(desc, aerogelDE, Form("mirror_optical_surface%d", isec), aerogelSurf, aerogelVol);
+ //aerogelSkin.isValid();
+ DetElement filterDE(det, Form("filter_de%d", isec), isec);
+ filterDE.setPlacement(filterPV);
+ //SkinSurface filterSkin(desc, filterDE, Form("mirror_optical_surface%d", isec), filterSurf, filterVol);
+ //filterSkin.isValid();
+
+
+ // BUILD MIRROR ====================================================================
+
+ // derived attributes
+ auto mirrorPos = Position(mirrorCenterX, 0., mirrorBackplane) + snoutFront;
+ double mirrorThetaRot = std::asin(mirrorCenterX/mirrorRadius);
+ double mirrorTheta1 = mirrorThetaRot - std::asin((mirrorCenterX-mirrorRmin)/mirrorRadius);
+ double mirrorTheta2 = mirrorThetaRot + std::asin((mirrorRmax-mirrorCenterX)/mirrorRadius);
+
+ // if debugging, draw full sphere
+ if(mirrorDebug>0) { mirrorTheta1=0; mirrorTheta2=M_PI; mirrorPhiw=2*M_PI; };
+
+ // solid and volume: create sphere at origin, with specified angular limits
+ Sphere mirrorSolid(
+ mirrorRadius-mirrorThickness,
+ mirrorRadius,
+ mirrorTheta1,
+ mirrorTheta2,
+ -mirrorPhiw/2.0,
+ mirrorPhiw/2.0
+ );
+ Volume mirrorVol("mirror_v", mirrorSolid, mirrorMat);
+ mirrorVol.setVisAttributes(mirrorVis);
+
+ // placement (note: transformations are in reverse order)
+ auto mirrorPV = gasvolVol.placeVolume(mirrorVol,
+ RotationZ(sectorRotation) // rotate about beam axis to sector
+ * Translation3D(0,0,-mirrorRadius) // move longitudinally so it intersects snoutFront
+ * Translation3D(mirrorPos.x(), mirrorPos.y(), mirrorPos.z()) // re-center to snoutFront
+ * RotationY(-mirrorThetaRot) // rotate about origin
+ );
+
+ // properties
+ DetElement mirrorDE(det, Form("mirror_de%d", isec), isec);
+ mirrorDE.setPlacement(mirrorPV);
+ SkinSurface mirrorSkin(desc, mirrorDE, Form("mirror_optical_surface%d", isec), mirrorSurf, mirrorVol);
+ mirrorSkin.isValid();
+
+
+ // BUILD SENSORS ====================================================================
+
+ // if debugging sphere properties, restrict number of sensors drawn
+ if(sensorSphDebug>0) { sensorSide = 2*M_PI*sensorSphRadius / 64; };
+
+ // solid and volume: single sensor module
+ Box sensorSolid(sensorSide/2., sensorSide/2., sensorThickness/2.);
+ Volume sensorVol("sensor_v", sensorSolid, sensorMat);
+ sensorVol.setVisAttributes(sensorVis);
+
+ // sensitivity
+ sensorVol.setSensitiveDetector(sens);
+
+ // SENSOR MODULE LOOP ------------------------
+ /* ALGORITHM: generate sphere of positions
+ * - NOTE: there are two coordinate systems here:
+ * - "global" the main ATHENA coordinate system
+ * - "generator" (vars end in `Gen`) is a local coordinate system for
+ * generating points on a sphere; it is related to the global system by
+ * a rotation; we do this so the "patch" (subset of generated
+ * positions) of sensors we choose to build is near the equator, where
+ * point distribution is more uniform
+ * - PROCEDURE: loop over `thetaGen`, with subloop over `phiGen`, each divided evenly
+ * - the number of points to generate depends how many sensors (+`sensorGap`)
+ * can fit within each ring of constant `thetaGen` or `phiGen`
+ * - we divide the relevant circumference by the sensor
+ * size(+`sensorGap`), and this number is allowed to be a fraction,
+ * because likely we don't care about generating a full sphere and
+ * don't mind a "seam" at the overlap point
+ * - if we pick a patch of the sphere near the equator, and not near
+ * the poles or seam, the sensor distribution will appear uniform
+ */
+
+ // initialize module number for this sector
+ int imod=1;
+
+ // thetaGen loop: iterate less than "0.5 circumference / sensor size" times
+ double nTheta = M_PI*sensorSphRadius / (sensorSide+sensorGap);
+ for(int t=0; t<(int)(nTheta+0.5); t++) {
+ double thetaGen = t/((double)nTheta) * M_PI;
+
+ // phiGen loop: iterate less than "circumference at this latitude / sensor size" times
+ double nPhi = 2*M_PI * sensorSphRadius * std::sin(thetaGen) / (sensorSide+sensorGap);
+ for(int p=0; p<(int)(nPhi+0.5); p++) {
+ double phiGen = p/((double)nPhi) * 2*M_PI - M_PI; // shift to [-pi,pi]
+
+ // determine global phi and theta
+ // - convert {radius,thetaGen,phiGen} -> {xGen,yGen,zGen}
+ double xGen = sensorSphRadius * std::sin(thetaGen) * std::cos(phiGen);
+ double yGen = sensorSphRadius * std::sin(thetaGen) * std::sin(phiGen);
+ double zGen = sensorSphRadius * std::cos(thetaGen);
+ // - convert {xGen,yGen,zGen} -> global {x,y,z} via rotation
+ double x = zGen;
+ double y = xGen;
+ double z = yGen;
+ // - convert global {x,y,z} -> global {phi,theta}
+ double phi = std::atan2(y,x);
+ double theta = std::acos(z/sensorSphRadius);
+
+ // cut spherical patch
+ // TODO [low priority]: instead of cutting a patch with complicated parameters,
+ // can we use something simpler such as Union or
+ // Intersection of solids?
+ // - applied on global coordinates
+ // - theta cuts are signed, since we will offset the patch along +x;
+ // from the offset position, toward barrel is positive, toward beam is negative
+ double thetaSigned = (x<0?-1:1) * theta;
+ // - position of yz planes, associated to theta cuts
+ double xmin = sensorSphRadius * std::sin(sensorSphPatchThetaMin/2);
+ double xmax = sensorSphRadius * std::sin(sensorSphPatchThetaMax/2);
+ // - we want a phi wedge, but offset from the origin to allow more width, so
+ // define phiCheck to account for the offset; the amount of the offset,
+ // and hence the width, is controlled by `sensorSphPatchWidthFactor`
+ double phiCheck = std::atan2(y,(x+sensorSphCenterX)/sensorSphPatchWidthFactor);
+ // - apply cuts (only if not debugging)
+ // - NOTE: use `x<xmax` for straight cuts, or `theta<sensorSphPatchThetaMax` for
+ // rounded cuts (which allows for more sensors)
+ if( ( std::fabs(phiCheck)<sensorSphPatchTaper && x>xmin && theta<sensorSphPatchThetaMax && z>0 )
+ || sensorSphDebug>0
+ ) {
+
+ // placement (note: transformations are in reverse order)
+ // - transformations operate on global coordinates; the corresponding
+ // generator coordinates are provided in the comments
+ auto sensorPV = gasvolVol.placeVolume(sensorVol,
+ RotationZ(sectorRotation) // rotate about beam axis to sector
+ * Translation3D(sensorSphCenterX, sensorSphCenterY, sensorSphCenterZ) // move sphere to specified center
+ * Translation3D(snoutFront.x(), snoutFront.y(), snoutFront.z()) // move sphere to reference position
+ * RotationX(phiGen) // rotate about `zGen`
+ * RotationZ(thetaGen) // rotate about `yGen`
+ * Translation3D(sensorSphRadius, 0., 0.) // push radially to spherical surface
+ * RotationY(M_PI/2) // rotate sensor to be compatible with generator coords
+ );
+
+ // properties
+ sensorPV.addPhysVolID("sector", isec).addPhysVolID("module", imod);
+ DetElement sensorDE(det, Form("sensor_de%d_%d", isec, imod), 10000*isec+imod);
+ sensorDE.setPlacement(sensorPV);
+ SkinSurface sensorSkin(desc, sensorDE, Form("sensor_optical_surface%d", isec), sensorSurf, sensorVol);
+ sensorSkin.isValid();
+
+ // increment sensor module number
+ imod++;
+
+ }; // end patch cuts
+ }; // end phiGen loop
+ }; // end thetaGen loop
+ // END SENSOR MODULE LOOP ------------------------
+
+
+ }; // END SECTOR LOOP //////////////////////////
+
+
+ // place gas volume
+ PlacedVolume gasvolPV = vesselVol.placeVolume(gasvolVol,Position(0, 0, 0));
+ DetElement gasvolDE(det, "gasvol_de", 0);
+ gasvolDE.setPlacement(gasvolPV);
+
+ // place mother volume (vessel)
+ Volume motherVol = desc.pickMotherVolume(det);
+ PlacedVolume vesselPV = motherVol.placeVolume(vesselVol,
+ Position(0, 0, vesselZ0) - snoutFront
+ );
+ vesselPV.addPhysVolID("system", detID);
+ det.setPlacement(vesselPV);
+
+ return det;
+};
+
+// clang-format off
+DECLARE_DETELEMENT(athena_DRICH, createDetector)
--
GitLab