From 69a4591bdaf273f41ac4bb3cb4b1e9fdae83be7c Mon Sep 17 00:00:00 2001 From: Whitney Armstrong <warmstrong@anl.gov> Date: Thu, 15 Jul 2021 18:49:39 +0000 Subject: [PATCH] Resolve "Detailed MRICH" --- athena.xml | 93 +++++++++- bin/debug_mrich | 8 + compact/ce_mrich.xml | 25 --- compact/central_tracker.xml | 2 +- compact/definitions.xml | 6 +- compact/display.xml | 2 + compact/mrich.xml | 96 ++++++++++ macro/mrich_vis.mac | 50 +++++ macro/vis.mac | 29 +++ src/GeometryHelpers.h | 13 +- src/MRich.cpp | 351 ++++++++++++++++++++++++++++++++++++ src/ce_MRICH.cpp | 101 ----------- 12 files changed, 635 insertions(+), 141 deletions(-) create mode 100755 bin/debug_mrich delete mode 100644 compact/ce_mrich.xml create mode 100644 compact/mrich.xml create mode 100644 macro/mrich_vis.mac create mode 100644 macro/vis.mac create mode 100644 src/MRich.cpp delete mode 100644 src/ce_MRICH.cpp diff --git a/athena.xml b/athena.xml index d54c46d6..f8a75669 100644 --- a/athena.xml +++ b/athena.xml @@ -4,11 +4,14 @@ <!-- Some information about detector --> <info name="Athena Detector" title="Athena Detector" - author="Jihee Kim, Whitney Armstrong" + author="Athena Collaboration" url="https://eicweb.phy.anl.gov/EIC/detectors/athena.git" status="development" version="v1 2021-03-16"> - <comment>Athena Detector</comment> + <comment>Athena Detector + - https://eicweb.phy.anl.gov/EIC/detectors/athena.git + - https://eicweb.phy.anl.gov/EIC/detectors/ip6.git + </comment> </info> <define> @@ -39,11 +42,31 @@ 4.0*eV 1.5 5.1*eV 1.5 "/> + <matrix name="ABSLENGTH__Pyrex" coldim="2" values=" + 1.0*eV 10.0*cm + 4.0*eV 10.0*cm + 5.1*eV 10.0*cm + "/> <matrix name= "REFLECTIVITY_mirror" coldim="2" values=" 1.0*eV 0.9 4.0*eV 0.9 5.1*eV 0.9 - "/> + "/> + <matrix name="RINDEX__Aerogel" coldim="2" values=" + 1.0*eV 1.030 + 4.0*eV 1.030 + 5.1*eV 1.030 + "/> + <matrix name="ABSLENGTH__Aerogel" coldim="2" values=" + 1.0*eV 4.0*cm + 4.0*eV 4.0*cm + 5.1*eV 4.0*cm + "/> + <matrix name="RINDEX__Acrylic" coldim="2" values=" + 1240*eV/1100 1.49 + 1240*eV/600 1.49 + 1240*eV/400 1.49 + "/> </properties> <includes> @@ -57,37 +80,87 @@ <fraction n="0.754" ref="N"/> <fraction n="0.234" ref="O"/> <fraction n="0.012" ref="Ar"/> - <property name="RINDEX" ref="RINDEX__Vacuum"/> + <property name="RINDEX" ref="RINDEX__Air"/> + <property name="ABSLENGTH" coldim="2" values="1*eV 200*m 5*eV 200*m"/> </material> <material name="N2cherenkov"> <D type="density" value="0.00125" unit="g/cm3"/> <composite n="1" ref="N"/> <property name="RINDEX" ref="RINDEX__N2"/> </material> - <material name="PyrexGlass"> + <material name="PyrexGlassOptical"> <D type="density" value="2.23" unit="g/cm3"/> <fraction n="0.806" ref="SiliconOxide"/> <fraction n="0.130" ref="BoronOxide"/> <fraction n="0.040" ref="SodiumOxide"/> <fraction n="0.023" ref="AluminumOxide"/> <property name="RINDEX" ref="RINDEX__Pyrex"/> + <property name="ABSLENGTH" ref="ABSLENGTH__Pyrex"/> + </material> + <material name="AerogelOptical"> + <D value="0.2" unit="g / cm3"/> + <fraction n="0.625" ref="SiliconOxide"/> + <fraction n="0.374" ref="SiliconOxide"/> + <fraction n="0.1" ref="C"/> + <property name="RINDEX" ref="RINDEX__Aerogel"/> + <property name="ABSLENGTH" ref="ABSLENGTH__Aerogel"/> + </material> + <material name="AcrylicOptical"> + <D type="density" value="1.18" unit="g/cm3"/> + <composite n="5" ref="C"/> + <composite n="2" ref="O"/> + <composite n="8" ref="H"/> + <property name="RINDEX" ref="RINDEX__Acrylic"/> </material> </materials> <surfaces> <comment> For the values of "finish", model and type, see TGeoOpticalSurface.h ! </comment> + <opticalsurface finish="polished" model="glisur" name="MirrorOpticalSurface" type="dielectric_metal" value="0"> <property name="REFLECTIVITY" ref="REFLECTIVITY_mirror"/> <property name="RINDEX" coldim="2" values="1.034*eV 1.5 4.136*eV 1.5"/> <!--<property name="EFFICIENCY" ref="EFFICIENCY0x8b77240"/>--> </opticalsurface> + + <!-- <opticalsurface name="mirror2" finish="polished" model="glisur" type="dielectric_dielectric"> <property name="REFLECTIVITY" coldim="2" values="1.034*eV 0.8 4.136*eV 0.9"/> <property name="EFFICIENCY" coldim="2" values="2.034*eV 0.8 4.136*eV 1.0"/> <property name="RINDEX" coldim="2" values="1.034*eV 1.5 4.136*eV 1.5"/> </opticalsurface> - <!-- Define the dimensions of the world volume --> + --> + + <opticalsurface finish="polished" model="unified" name="MRICH_MirrorOpticalSurface" type="dielectric_metal" value="0"> + </opticalsurface> + <opticalsurface finish="polished" model="unified" name="MRICH_LensOpticalSurface" type="dielectric_dielectric" value="0"> + <property name="REFLECTIVITY" coldim="2" values="1240*eV/1100 0.08 1240*eV/400 0.08"/> + <!-- + <property name="RINDEX" coldim="2" values="2.034*eV 1.56 4.136*eV 1.56"/> + <property name="SPECULARLOBECONSTANT" coldim="2" values="2.034*eV 0.3 4.136*eV 0.3 "/> + <property name="SPECULARSPIKECONSTANT" coldim="2" values="2.034*eV 0.2 4.136*eV 0.2 "/> + <property name="BACKSCATTERCONSTANT" coldim="2" values="2.034*eV 0.2 4.136*eV 0.2 "/> + --> + </opticalsurface> + + <opticalsurface finish="polished" model="unified" name="MRICH_PhotoSensorOpticalSurface" type="dielectric_dielectric" value="0"> + <!-- + <property name="RINDEX" coldim="2" values="2.034*eV 1.35 4.136*eV 1.40"/> + <property name="SPECULARLOBECONSTANT" coldim="2" values="2.034*eV 0.3 4.136*eV 0.3 "/> + <property name="SPECULARSPIKECONSTANT" coldim="2" values="2.034*eV 0.2 4.136*eV 0.2 "/> + <property name="BACKSCATTERCONSTANT" coldim="2" values="2.034*eV 0.2 4.136*eV 0.2 "/> + --> + </opticalsurface> + + <opticalsurface finish="polished" model="unified" name="MRICH_AerogelOpticalSurface" type="dielectric_dielectric" value="0"> + <!-- + <property name="RINDEX" coldim="2" values="2.034*eV 1.010 4.136*eV 1.010"/> + <property name="SPECULARLOBECONSTANT" coldim="2" values="2.034*eV 0.3 4.136*eV 0.3 "/> + <property name="SPECULARSPIKECONSTANT" coldim="2" values="2.034*eV 0.2 4.136*eV 0.2 "/> + <property name="BACKSCATTERCONSTANT" coldim="2" values="2.034*eV 0.2 4.136*eV 0.2 "/> + --> + </opticalsurface> </surfaces> <limits> @@ -189,7 +262,9 @@ <include ref="compact/tof_barrel.xml"/> - <!--include ref="compact/rwell_tracker_barrel.xml"/--> + <!-- + <include ref="compact/rwell_tracker_barrel.xml"/> + --> <include ref="compact/cb_DIRC.xml"/> <!-- When changing magnet, also select dimensions in definitions.xml. --> @@ -206,7 +281,7 @@ <include ref="compact/hcal.xml"/> <!--include ref="compact/ce_GEM.xml"/--> <!--include ref="compact/gem_tracker_endcap.xml"/--> - <include ref="compact/ce_mrich.xml"/> + <include ref="compact/mrich.xml"/> <include ref="compact/tof_endcap.xml"/> <include ref="compact/forward_trd.xml"/> <include ref="compact/gaseous_rich.xml"/> @@ -216,12 +291,12 @@ <include ref="ip6/far_forward_romanpots.xml"/> <include ref="ip6/far_forward_detectors.xml"/> + <!-- <include ref="compact/mm_tracker_barrel.xml"/> <include ref="compact/cb_VTX_Barrel.xml"/> --> - <readouts> </readouts> diff --git a/bin/debug_mrich b/bin/debug_mrich new file mode 100755 index 00000000..30bff695 --- /dev/null +++ b/bin/debug_mrich @@ -0,0 +1,8 @@ +#!/bin/bash + +npsim --enableG4GPS \ + --compact athena.xml \ + --runType vis \ + --macro macro/mrich_vis.mac \ + -O derp.root + diff --git a/compact/ce_mrich.xml b/compact/ce_mrich.xml deleted file mode 100644 index 7b8860aa..00000000 --- a/compact/ce_mrich.xml +++ /dev/null @@ -1,25 +0,0 @@ -<?xml version="1.0" encoding="UTF-8"?> -<lccdd> - <define> - <constant name="ce_MRICHRMin" value="15*cm"/> - <constant name="ce_MRICHRMax" value="BarrelTracking_rmax"/> - <constant name="ce_MRICHLength" value="BackwardCherenkov_length"/> - <constant name="ce_MRICHZMin" value="-BarrelTracking_length/2.0 - BackwardTracking_length"/> - - </define> - - <detectors> - <detector id="ce_MRICH_ID" name="ce_MRICH" type="refdet_ce_MRICH" readout="MRICHHits" vis="BlueVis" material="Air"> - <dimensions rmin="ce_MRICHRMin" rmax="ce_MRICHRMax" length="ce_MRICHLength" z="ce_MRICHZMin"/> - <modules material="Quartz" thickness="10*cm" width="10*cm" gap="1*cm" /> - </detector> - </detectors> - - <readouts> - <readout name="MRICHHits"> - <segmentation type="CartesianGridXY" grid_size_x="3*mm" grid_size_y="3*mm" /> - <id>system:8,sector:8,module:14,x:32:-16,y:-16</id> - </readout> - </readouts> - -</lccdd> diff --git a/compact/central_tracker.xml b/compact/central_tracker.xml index 2ac1351d..5785b111 100644 --- a/compact/central_tracker.xml +++ b/compact/central_tracker.xml @@ -57,7 +57,7 @@ <constant name="TrackerEndcapInnerLayer_dz" value="TrackerEndcapInner_length/TrackerEndcapInner_nLayers"/> <constant name="TrackerEndcapOuter_zmin" value="TrackerBarrelOuter_length/2.0"/> - <constant name="TrackerEndcapOuter_zmax" value="TrackerEndcapOuter_zmin + 72"/> + <constant name="TrackerEndcapOuter_zmax" value="TrackerEndcapOuter_zmin + 70"/> <constant name="TrackerEndcapOuter_length" value="TrackerEndcapOuter_zmax- TrackerEndcapOuter_zmin"/> <constant name="TrackerEndcapOuterLayer_dz" value="TrackerEndcapOuter_length/TrackerEndcapOuter_nLayers"/> diff --git a/compact/definitions.xml b/compact/definitions.xml index 5d4bf1fa..eb58560d 100644 --- a/compact/definitions.xml +++ b/compact/definitions.xml @@ -5,8 +5,8 @@ <constant name="world_z" value="100*m"/> <constant name="Pi" value="3.14159265359"/> - <constant name="mil" value="0.0254*mm"/> + <constant name="inch" value="2.54*cm"/> <comment> ------------ @@ -232,7 +232,7 @@ Unused IDs: 132-139 </comment> - <constant name="ce_MRICH_ID" value="130"/> + <constant name="MRICH_ID" value="130"/> <constant name="ce_GEM_ID" value="131"/> <comment> @@ -340,7 +340,7 @@ <constant name="ForwardTRD_length" value="10.0*cm"/> <constant name="ForwardTOF_length" value="3.0*cm"/> - <constant name="BackwardCherenkov_length" value="15.0*cm"/> + <constant name="BackwardCherenkov_length" value="20.0*cm"/> <constant name="BackwardTOF_length" value="3.0*cm"/> <comment> Total length of PID detectors above</comment> diff --git a/compact/display.xml b/compact/display.xml index 27e6d71b..bd358cc3 100644 --- a/compact/display.xml +++ b/compact/display.xml @@ -6,9 +6,11 @@ <vis name="AnlOrange" alpha="1" r="255/256" g="121/256" b="0"/> <vis name="AnlRed" alpha="1" r="205/256" g="32/256" b="44/256"/> <vis name="AnlGold" alpha="1" r="248/256" g="178/256" b="0/256"/> + <vis name="AnlGold_1" alpha="1" r="248/256" g="178/256" b="0/256" visible="true" drawingStyle="wireframe"/> <vis name="AnlBlue" alpha="1" r="0/256" g="96/256" b="156/256"/> <vis name="AnlTeal" alpha="1" r="0/256" g="161/256" b="156/256"/> <vis name="AnlGray" alpha="1" r="102/256" g="102/256" b="102/256"/> + <vis name="AnlGray_1" alpha="1" r="102/256" g="102/256" b="102/256" visible="true" drawingStyle="wireframe"/> <vis name="AnlLight_Gray" alpha="1" r="209/256" g="209/256" b="209/256"/> <vis name="AnlOff_White" alpha="1" r="242/256" g="242/256" b="242/256"/> <vis name="AnlDelta_Red" alpha="1" r="161/256" g="43/256" b="47/256"/> diff --git a/compact/mrich.xml b/compact/mrich.xml new file mode 100644 index 00000000..e7ff797b --- /dev/null +++ b/compact/mrich.xml @@ -0,0 +1,96 @@ +<lccdd> + <comment> MRICH </comment> + <define> + <constant name="MRICH_rmin" value="15*cm"/> + <constant name="MRICH_rmax" value="BarrelTracking_rmax"/> + <constant name="MRICH_length" value="BackwardCherenkov_length"/> + <constant name="MRICH_zmin" value="-BarrelTracking_length/2.0 - BackwardTracking_length"/> + <constant name="MRICHAerogel_thickness" value="30.0*mm"/> + <constant name="MRICHAerogel_width" value="126.5*mm"/> + <constant name="MRICHFoam_thickness" value="2*mm"/> + <constant name="MRICHFresnelLens_thickness" value="0.06*inch"/> + <constant name="MRICHAerogelLensGap_thickness" value="2*mm"/> + <constant name="MRICHLensMirrorGap_thickness" value="2*mm"/> + <constant name="MRICHPhotoDet_thickness" value="1.5*mm"/> + <constant name="MRICHPhotoDet_length" value="48.5*mm"/> + <constant name="MRICHGlassWindow_width" value="103.5*mm"/> + <constant name="MRICHGlassPhotoDet_thickness" value="2.0*mm"/> + <constant name="MRICHRearExtraSpace_thickness" value="10.0*mm"/> + <constant name="MRICHLensPhotoDet_length" value="136.4*mm"/> + <constant name="MRICHMirror_thickness" value="2.0*mm"/> + <constant name="MRICHMirror_length" value="MRICHLensPhotoDet_length - MRICHLensMirrorGap_thickness"/> + <constant name="MRICHMirror_width1" value="MRICHAerogel_width "/> + <constant name="MRICHMirror_width2" value="MRICHGlassWindow_width"/> + + <constant name="MRICHFresnelLensEffectiveDiameter" value="6.0*inch"/> + <constant name="MRICHFresnelLensGroove_pitch" value="inch/125"/> + + <constant name="MRICHCarbonFrame_thickness" value="1.0*mm"/> + <constant name="MRICHCarbonFrame_width" value="MRICHAerogel_width+2.0*MRICHFoam_thickness + 2.0*MRICHCarbonFrame_thickness"/> + + <constant name="MRICHModules_nx" value="floor((MRICH_rmax-MRICH_rmin)/MRICHCarbonFrame_width)"/> + <constant name="MRICHModules_ny" value="floor((MRICH_rmax-MRICH_rmin)/MRICHCarbonFrame_width)"/> + + <constant name="MRICHCarbonFrame_length" + value="MRICHAerogel_thickness + + 2.0*MRICHCarbonFrame_thickness + + 2.0*MRICHFoam_thickness + + MRICHAerogelLensGap_thickness + + MRICHFresnelLens_thickness + + MRICHLensPhotoDet_length + + MRICHGlassPhotoDet_thickness + + MRICHRearExtraSpace_thickness "/> + </define> + + <limits> + </limits> + + <regions> + </regions> + + <display> + </display> + + <detectors> + <detector id="MRICH_ID" name="MRICH" type="athena_MRICH" + readout="MRICHHits" + projective="false" + vis="InvisibleWithDaughters" material="Air"> + <dimensions rmin="MRICH_rmin" rmax="MRICH_rmax" length="MRICHCarbonFrame_length" zmin="-1.0*abs(MRICH_zmin)"/> + <module name="MRICH_module1" vis="InvisibleWithDaughters" + width="MRICHCarbonFrame_width" + height="MRICHCarbonFrame_width" + length="MRICHCarbonFrame_length"> + <frame vis="AnlGray" thickness="MRICHCarbonFrame_thickness" material="CarbonFiber"/> + <aerogel vis="AnlTeal" length="MRICHAerogel_thickness" width="MRICHAerogel_width" material="AerogelOptical"> + <frame vis="AnlGold_1" thickness="MRICHFoam_thickness" material="PolystyreneFoam" /> + </aerogel> + <lens vis="AnlViolet" thickness="MRICHFresnelLens_thickness" + pitch="MRICHFresnelLensGroove_pitch" focal_length="6.0*inch" + effective_diameter="MRICHFresnelLensEffectiveDiameter" + width="MRICHAerogel_width" + material="AcrylicOptical"/> + <mirror vis="AnlGray" x1="MRICHMirror_width1" x2="MRICHMirror_width2" length="MRICHMirror_length" + surface="MRICH_MirrorOpticalSurface" thickness="MRICHMirror_thickness" material="AluminumOxide"/> + <photodet width="MRICHGlassWindow_width" thickness="MRICHGlassPhotoDet_thickness" material="PyrexGlassOptical"> + <sensor nx="2" ny="2" thickness="MRICHPhotoDet_thickness" width="MRICHPhotoDet_length" material="SiliconOxide"/> + </photodet> + </module> + <!-- + <layer id="1"> + <array nx="4" ny="4" module="MRICH_module1"> + <position x="0" y="0" z="0"/> + </array> + </layer> + --> + </detector> + </detectors> + + <readouts> + <readout name="MRICHHits"> + <segmentation type="CartesianGridXY" grid_size_x="3*mm" grid_size_y="3*mm" /> + <id>system:8,module:14,sensor:8,x:32:-16,y:-16</id> + </readout> + </readouts> + +</lccdd> diff --git a/macro/mrich_vis.mac b/macro/mrich_vis.mac new file mode 100644 index 00000000..63ee5d45 --- /dev/null +++ b/macro/mrich_vis.mac @@ -0,0 +1,50 @@ +/run/initialize + +#/vis/open OGL 800x800-0+0 +/vis/open OGLIQt 800x800-0+0 +# +/vis/drawVolume +/vis/viewer/set/viewpointThetaPhi 20 30 +/vis/viewer/zoom 1.1 +#/vis/viewer/set/style wireframe +#/vis/scene/add/axes 0 0 0 1 m +/vis/scene/add/trajectories rich smooth +/vis/modeling/trajectories/create/drawByCharge +#/vis/modeling/trajectories/drawByCharge-0/default/setDrawStepPts true +#/vis/modeling/trajectories/drawByCharge-0/default/setStepPtsSize 2 +/vis/scene/add/hits +/vis/scene/endOfEventAction accumulate 2 + +#/vis/filtering/trajectories/create/particleFilter +#/vis/filtering/trajectories/particleFilter-0/add gamma +#/vis/filtering/trajectories/particleFilter-0/invert true +#/vis/filtering/trajectories/particleFilter-0/verbose true +#/vis/filtering/trajectories/particleFilter-0/active true + + +/vis/ogl/set/displayListLimit 50000 + +/vis/viewer/flush + + +/gps/verbose 2 +/gps/particle pi- +/gps/number 1 + +/gps/ene/type Gauss +/gps/ene/mono 8.0 GeV +/gps/ene/sigma 1.0 GeV + +/gps/pos/type Volume +/gps/pos/shape Cylinder +/gps/pos/centre 0.0 0.0 0.0 cm +/gps/pos/radius 0.001 cm +/gps/pos/halfz 1 cm +/gps/position 0 0 0 cm + +#/gps/direction 0 0.1 1.0 +/gps/ang/type iso +/gps/ang/mintheta 10 degree +/gps/ang/maxtheta 20 degree + +/run/beamOn 1 diff --git a/macro/vis.mac b/macro/vis.mac new file mode 100644 index 00000000..3823c07c --- /dev/null +++ b/macro/vis.mac @@ -0,0 +1,29 @@ +/vis/open OGL 800x800-0+0 + +/vis/drawVolume +/vis/viewer/set/viewpointThetaPhi 30 30 +#/vis/viewer/zoom 30. +#/vis/viewer/set/style wireframe +#/vis/scene/add/axes 0 0 0 1 m +/vis/scene/add/trajectories rich smooth +/vis/modeling/trajectories/create/drawByCharge +#/vis/modeling/trajectories/drawByCharge-0/default/setDrawStepPts true +#/vis/modeling/trajectories/drawByCharge-0/default/setStepPtsSize 2 +/vis/scene/add/hits +/vis/scene/endOfEventAction accumulate 2 + +#/vis/filtering/trajectories/create/particleFilter +#/vis/filtering/trajectories/particleFilter-0/add gamma +#/vis/filtering/trajectories/particleFilter-0/invert true +#/vis/filtering/trajectories/particleFilter-0/verbose true +#/vis/filtering/trajectories/particleFilter-0/active true + + +/vis/ogl/set/displayListLimit 500000 + +/vis/viewer/flush + +/run/beamOn 1 + +#/control/execute macro/gun.mac + diff --git a/src/GeometryHelpers.h b/src/GeometryHelpers.h index 5ec55025..5194e100 100644 --- a/src/GeometryHelpers.h +++ b/src/GeometryHelpers.h @@ -5,9 +5,18 @@ // some utility functions that can be shared namespace athena::geo { -typedef ROOT::Math::XYPoint Point; +using Point = ROOT::Math::XYPoint; -// fill rectangles in a ring +/** Fill rectangles in a ring (disk). + * + * @param ref 2D reference point. + * @param sx x side length + * @param sy y side length + * @param rmin inner radius of disk + * @param rmax outer radius of disk to fill + * @param phmin phi min + * @param phmax phi max + */ std::vector<Point> fillRectangles(Point ref, double sx, double sy, double rmin, double rmax, double phmin = -M_PI, double phmax = M_PI); // fill squares in a ring diff --git a/src/MRich.cpp b/src/MRich.cpp new file mode 100644 index 00000000..33d1059c --- /dev/null +++ b/src/MRich.cpp @@ -0,0 +1,351 @@ +// +// Author : Whit Armstrong (warmstrong@anl.gov) +// +#include <XML/Helper.h> +#include "TMath.h" +#include "TString.h" +#include "DDRec/Surface.h" +#include "DDRec/DetectorData.h" +#include "DD4hep/OpticalSurfaces.h" +#include "DD4hep/DetFactoryHelper.h" +#include "DD4hep/Printout.h" +#include "GeometryHelpers.h" +#include "Math/Vector3D.h" +#include "Math/AxisAngle.h" +#include "Math/VectorUtil.h" + +using namespace std; +using namespace dd4hep; +using namespace dd4hep::rec; + +using Placements = vector<PlacedVolume>; + + + +static Ref_t createDetector(Detector& description, xml::Handle_t e, SensitiveDetector sens){ + xml_det_t x_det = e; + Material air = description.material("AirOptical"); + Material vacuum = description.vacuum(); + string det_name = x_det.nameStr(); + //xml::Component pos = x_det.position(); + DetElement sdet(det_name, x_det.id()); + Assembly assembly(det_name); + sens.setType("photoncounter"); + OpticalSurfaceManager surfMgr = description.surfaceManager(); + + bool projective = getAttrOrDefault(x_det, _Unicode(projective), false); + + PlacedVolume pv; + + map<string, Volume> modules; + map<string, Placements> sensitives; + map<string, Volume> module_assemblies; + std::map<std::string,DetElement> module_assembly_delements; + + int n_sensor = 1; + + xml::Component dims = x_det.dimensions(); + auto rmin = dims.rmin(); + auto rmax = dims.rmax(); + auto length = dims.length(); + auto zmin = dims.zmin(); + + // expect only one module (for now) + xml_comp_t x_mod = x_det.child(_U(module)); + string mod_name = x_mod.nameStr(); + double mod_width = getAttrOrDefault(x_mod, _U(width), 130.0 * mm); + double mod_height = getAttrOrDefault(x_mod, _U(height), 130.0 * mm); + double mod_length = getAttrOrDefault(x_mod, _U(length), 130.0 * mm); + + // various components + xml_comp_t x_frame = x_mod.child(_Unicode(frame)); + xml_comp_t x_aerogel = x_mod.child(_Unicode(aerogel)); + xml_comp_t x_lens = x_mod.child(_Unicode(lens)); + xml_comp_t x_mirror = x_mod.child(_Unicode(mirror)); + xml_comp_t x_photodet = x_mod.child(_Unicode(photodet)); + + // module + Box m_solid(mod_width / 2.0, mod_height / 2.0, mod_length / 2.0); + Volume m_volume(mod_name, m_solid, air); + m_volume.setVisAttributes(description.visAttributes(x_mod.visStr())); + DetElement mod_de( mod_name + std::string("_mod_") + std::to_string(1), 1); + + // todo module frame + double frame_thickness = getAttrOrDefault(x_frame, _U(thickness), 2.0 * mm); + + // aerogel box + xml_comp_t x_aerogel_frame = x_aerogel.child(_Unicode(frame)); + double aerogel_width = getAttrOrDefault(x_aerogel, _U(width), 130.0 * mm); + double aerogel_length = getAttrOrDefault(x_aerogel, _U(length), 130.0 * mm); + double foam_thickness = getAttrOrDefault(x_aerogel_frame, _U(thickness), 2.0 * mm); + Material foam_mat = description.material(x_aerogel_frame.materialStr()); + Material aerogel_mat = description.material(x_aerogel.materialStr()); + auto aerogel_vis = getAttrOrDefault<std::string>(x_aerogel, _U(vis), std::string("InvisibleWithDaughters")); + auto foam_vis = getAttrOrDefault<std::string>(x_aerogel_frame, _U(vis), std::string("RedVis")); + + // aerogel foam frame + Box foam_box(aerogel_width / 2.0 + foam_thickness, aerogel_width / 2.0 + foam_thickness, (aerogel_length + foam_thickness) / 2.0); + Box aerogel_sub_box(aerogel_width / 2.0, aerogel_width / 2.0, (aerogel_length + foam_thickness) / 2.0); + SubtractionSolid foam_frame_solid(foam_box, aerogel_sub_box, Position(0, 0, foam_thickness)); + Volume foam_vol(mod_name+"_aerogel_frame", foam_frame_solid, foam_mat); + foam_vol.setVisAttributes(description.visAttributes(foam_vis)); + double foam_frame_zpos = -mod_length / 2.0 + frame_thickness + (aerogel_length + foam_thickness) / 2.0; + m_volume.placeVolume(foam_vol,Position(0,0,foam_frame_zpos)); + + // aerogel + Box aerogel_box(aerogel_width / 2.0, aerogel_width / 2.0, (aerogel_length) / 2.0); + Volume aerogel_vol(mod_name+"_aerogel", aerogel_box, aerogel_mat); + aerogel_vol.setVisAttributes(description.visAttributes(aerogel_vis)); + double aerogel_zpos = foam_frame_zpos + foam_thickness / 2.0; + pv = m_volume.placeVolume(aerogel_vol,Position(0,0,aerogel_zpos)); + DetElement aerogel_de(mod_de, mod_name + std::string("_aerogel_de") + std::to_string(1), 1); + aerogel_de.setPlacement(pv); + + auto aerogel_surf = surfMgr.opticalSurface(dd4hep::getAttrOrDefault(x_aerogel, _Unicode(surface), "MRICH_AerogelOpticalSurface")); + SkinSurface skin0(description, aerogel_de, Form("MRICH_aerogel_skin_surface_%d", 1), aerogel_surf, aerogel_vol); + skin0.isValid(); + + // Fresnel Lens + // - The lens has a constant groove pitch (delta r) as opposed to fixing the groove height. + // - The lens area outside of the effective diamtere is flat. + // - The grooves are not curved, rather they are polycone shaped, ie a flat approximating the curvature. + auto lens_vis = getAttrOrDefault<std::string>(x_lens, _U(vis), std::string("AnlBlue")); + double groove_pitch = getAttrOrDefault(x_lens, _Unicode(pitch), 0.5 * mm); + double lens_f = getAttrOrDefault(x_lens, _Unicode(focal_length), 6.0*2.54*cm); + double eff_diameter = getAttrOrDefault(x_lens, _Unicode(effective_diameter), 130.0 * mm); + double lens_width = getAttrOrDefault(x_lens, _Unicode(width), 6.7*2.54*cm); + double center_thickness = getAttrOrDefault(x_lens, _U(thickness), 2.0 * mm); + + double n_acrylic = 1.49; + double lens_curvature = 1.0 / (lens_f*(n_acrylic - 1.0)); + double full_ring_rmax = std::min(eff_diameter / 2.0, lens_width/2.0); + + double N_grooves = std::ceil((full_ring_rmax) / groove_pitch); + double groove_last_rmin = (N_grooves - 1) * groove_pitch; + double groove_last_rmax = N_grooves * groove_pitch; + + auto groove_sagitta = [&](double r) { return lens_curvature * std::pow(r, 2) / (1.0 + 1.0); }; + double lens_thickness = groove_sagitta(groove_last_rmax) - groove_sagitta(groove_last_rmin) + center_thickness; + + Material lens_mat = description.material(x_lens.materialStr()); + Box lens_box(lens_width / 2.0, lens_width / 2.0, (center_thickness) / 2.0); + SubtractionSolid flat_lens(lens_box, Tube(0.0, full_ring_rmax, 2 * center_thickness)); + + Assembly lens_vol(mod_name + "_lens"); + Volume flatpart_lens_vol( "flatpart_lens", flat_lens, lens_mat); + lens_vol.placeVolume(flatpart_lens_vol);//,Position(0,0,lens_zpos)); + + Solid fresnel_lens_solid; + + int i_groove = 0; + double groove_rmax = groove_pitch; + double groove_rmin = 0; + + while ( groove_rmax <= full_ring_rmax ) { + double dZ = groove_sagitta(groove_rmax) - groove_sagitta(groove_rmin); + //std::cout << " dZ = " << dZ << ", lens_thickness = " << lens_thickness << "\n"; + //std::cout << "groove_rmin = " << groove_rmin << "\n"; + //std::cout << "groove_rmax = " << groove_rmax << "\n"; + Polycone groove_solid(0, 2.0 * M_PI, + {groove_rmin, groove_rmin, groove_rmin}, + {groove_rmax, groove_rmax, groove_rmin}, + {-lens_thickness/2.0, lens_thickness/2.0-dZ, lens_thickness/2.0}); + Volume lens_groove_vol("lens_groove_" + std::to_string(i_groove), groove_solid, lens_mat); + lens_vol.placeVolume(lens_groove_vol); //,Position(0,0,lens_zpos)); + //Volume groove_vol(groove_solid, lens_mat, par->name.c_str(), 0, 0, 0); + //new G4PVPlacement(0, par->pos, Groove_log[i], par->name.c_str(), motherLV, false, 0, OverlapCheck()); + //phi1 = phi1 + halfpi; //g4 pre-defined: halfpi=pi/2 + //Tube sub_cylinder(r0, r1, 3*eff_diameter); + //IntersectionSolid groove_solid(lens_box,lens_sphere, Position(0,0,-eff_diameter/2.0 + lens_thickness/2.0+(t-lens_t)/2.0 )); + //IntersectionSolid lens_ring(groove_solid, sub_cylinder); + //if (i_groove == 0) { + // fresnel_lens_solid = groove_solid; + //} else { + // fresnel_lens_solid = UnionSolid(fresnel_lens_solid, groove_solid); + //} + //r0 = r1; + //if(i_groove > 3) { + // SubtractionSolid flat_lens(lens_box,Tube(0.0, r0, 3*eff_diameter)); + // fresnel_lens_solid = UnionSolid(fresnel_lens_solid, flat_lens); + // break; // temporary + //} + i_groove++; + groove_rmin = (i_groove )*groove_pitch; + groove_rmax = (i_groove+1)*groove_pitch; + } + //fresnel_lens_solid = UnionSolid(fresnel_lens_solid, flat_lens); + //Volume lens_vol(mod_name + "_lens", fresnel_lens_solid, lens_mat); + + lens_vol.setVisAttributes(description.visAttributes(lens_vis)); + double lens_zpos = aerogel_zpos +aerogel_length/ 2.0 + foam_thickness + lens_thickness/2.0; + pv = m_volume.placeVolume(lens_vol,Position(0,0,lens_zpos)); + DetElement lens_de(mod_de, mod_name + std::string("_lens_de") + std::to_string(1), 1); + lens_de.setPlacement(pv); + + auto surf = surfMgr.opticalSurface(dd4hep::getAttrOrDefault(x_lens, _Unicode(surface), "MRICH_LensOpticalSurface")); + SkinSurface skin(description, lens_de, Form("MRichFresnelLens_skin_surface_%d", 1), surf, lens_vol); + skin.isValid(); + + // mirror + auto mirror_vis = getAttrOrDefault<std::string>(x_mirror, _U(vis), std::string("AnlGray")); + double mirror_x1 = getAttrOrDefault(x_mirror, _U(x1), 100.0 * mm); + double mirror_x2 = getAttrOrDefault(x_mirror, _U(x2), 80.0 * mm); + double mirror_length = getAttrOrDefault(x_mirror, _U(length), 130.0 * mm); + double mirror_thickness = getAttrOrDefault(x_mirror, _U(thickness), 2.0 * mm); + double outer_x1 = (mirror_x1+mirror_thickness)/2.0; + double outer_x2 = (mirror_x2+mirror_thickness)/2.0; + Trd2 outer_mirror_trd(outer_x1, outer_x2, outer_x1, outer_x2, mirror_length/2.0); + Trd2 inner_mirror_trd(mirror_x1 / 2.0, mirror_x2 / 2.0, mirror_x1 / 2.0,mirror_x2 / 2.0, mirror_length/2.0+0.1*mm); + SubtractionSolid mirror_solid(outer_mirror_trd, inner_mirror_trd); + Material mirror_mat = description.material(x_mirror.materialStr()); + Volume mirror_vol(mod_name+"_mirror", mirror_solid, mirror_mat); + double mirror_zpos = lens_zpos + lens_thickness/2.0 + foam_thickness + mirror_length/2.0; + pv = m_volume.placeVolume(mirror_vol,Position(0,0,mirror_zpos)); + DetElement mirror_de(mod_de, mod_name + std::string("_mirror_de") + std::to_string(1), 1); + mirror_de.setPlacement(pv); + + auto mirror_surf = surfMgr.opticalSurface(dd4hep::getAttrOrDefault(x_mirror, _Unicode(surface), "MRICH_MirrorOpticalSurface")); + SkinSurface skin1(description, mirror_de, Form("MRICH_mirror_skin_surface_%d", 1), mirror_surf, mirror_vol); + skin1.isValid(); + + // photon detector + xml_comp_t x_photodet_sensor = x_photodet.child(_Unicode(sensor)); + auto photodet_vis = getAttrOrDefault<std::string>(x_photodet, _U(vis), std::string("AnlRed")); + double photodet_width = getAttrOrDefault(x_photodet, _U(width), 130.0 * mm); + double photodet_thickness = getAttrOrDefault(x_photodet, _U(thickness), 2.0 * mm); + double sensor_thickness = getAttrOrDefault(x_photodet_sensor, _U(thickness), 2.0 * mm); + Material photodet_mat = description.material(x_photodet.materialStr()); + Material sensor_mat = description.material(x_photodet_sensor.materialStr()); + int sensor_nx = getAttrOrDefault(x_photodet_sensor, _Unicode(nx), 2); + int sensor_ny = getAttrOrDefault(x_photodet_sensor, _Unicode(ny), 2); + + Box window_box(photodet_width/2.0,photodet_width/2.0,photodet_thickness/2.0); + Volume window_vol(mod_name+"_window", window_box, photodet_mat); + window_vol.setSensitiveDetector(sens); + double window_zpos = mirror_zpos + mirror_length/2.0+photodet_thickness/2.0; + pv = m_volume.placeVolume(window_vol,Position(0,0,window_zpos)); + DetElement comp_de(mod_de, std::string("mod_sensor_de_") + std::to_string(1) , 1); + comp_de.setPlacement(pv); + pv.addPhysVolID("sensor", n_sensor); + + //for (size_t ic = 0; ic < sensVols.size(); ++ic) { + // PlacedVolume sens_pv = sensVols[ic]; + // DetElement comp_de(mod_de, std::string("de_") + sens_pv.volume().name(), ic + 1); + // comp_de.setPlacement(sens_pv); + // // Acts::ActsExtension* sensorExtension = new Acts::ActsExtension(); + // //// sensorExtension->addType("sensor", "detector"); + // // comp_de.addExtension<Acts::ActsExtension>(sensorExtension); + // //// comp_de.setAttributes(description, sens_pv.volume(), x_layer.regionStr(), + // //// x_layer.limitsStr(), + // //// xml_det_t(xmleles[m_nam]).visStr()); + //} + //DetElement window_de(sdet, mod_name + std::string("_window_de") + std::to_string(1), 1); + //window_de.setPlacement(pv); + + window_vol.setSensitiveDetector(sens); + sensitives[mod_name].push_back(pv); + ++n_sensor; + modules[mod_name] = m_volume; + module_assembly_delements[mod_name] = mod_de; + // end module + + int i_mod = 1; + // detector envelope + Tube envShape(rmin, rmax, length / 2., 0., 2 * M_PI); + Volume envVol("MRICH_Envelope", envShape, air); + envVol.setVisAttributes(description.visAttributes(x_det.visStr())); + + // place modules in the sectors (disk) + auto points = athena::geo::fillSquares({0., 0.}, mod_width, rmin, rmax); + + // mod_name = ... + Placements& sensVols = sensitives[mod_name]; + auto mod_v = modules[mod_name]; + // determine module direction, always facing z = 0 + double roty = dims.zmin() < 0. ? -M_PI : 0 ; + int imod = 1; + for (auto& p : points) { + ROOT::Math::XYZVector x_location(p.x(), p.y(), zmin+std::signbit(zmin)*mod_length/2.0); + ROOT::Math::XYZVector z_dir(0, 0, 1); + ROOT::Math::XYZVector x_dir(1, 0, 0); + ROOT::Math::XYZVector rot_axis = x_location.Cross(z_dir); + double rot_angle = ROOT::Math::VectorUtil::Angle(z_dir,x_location); + ROOT::Math::AxisAngle proj_rot(rot_axis,-1.0*rot_angle); + ROOT::Math::AxisAngle grid_fix_rot(x_location,0.0*rot_angle); + auto new_x_dir = grid_fix_rot*x_dir; + // operations are inversely ordered + Transform3D tr = Translation3D(p.x(), p.y(), 0.) // move to position + * RotationY(roty); // facing z = 0. + + if(projective) { + tr = Translation3D(p.x(), p.y(), 0.) // move to position + * grid_fix_rot // keep the modules oriented vertially + * proj_rot // projective rotation + * RotationY(roty); // facing z = 0. + } + // mod placement + pv = envVol.placeVolume(mod_v, tr); + pv.addPhysVolID("module", i_mod); + + auto mod_det_element = module_assembly_delements[mod_name].clone(mod_name + "__" + std::to_string(i_mod)); + mod_det_element.setPlacement(pv); + sdet.add(mod_det_element); + i_mod++; + } + + // place envelope + Volume motherVol = description.pickMotherVolume(sdet); + PlacedVolume envPV = motherVol.placeVolume(envVol, Position(0, 0, zmin)); + envPV.addPhysVolID("system", x_det.id()); + sdet.setPlacement(envPV); + return sdet; +} + + +//void addModules(Volume &mother, xml::DetElement &detElem, Detector &description, SensitiveDetector &sens) +//{ +// xml::Component dims = detElem.dimensions(); +// xml::Component mods = detElem.child(_Unicode(modules)); +// +// auto rmin = dims.rmin(); +// auto rmax = dims.rmax(); +// +// auto mThick = mods.attr<double>(_Unicode(thickness)); +// auto mWidth = mods.attr<double>(_Unicode(width)); +// auto mGap = mods.attr<double>(_Unicode(gap)); +// +// auto modMat = description.material(mods.materialStr()); +// auto gasMat = description.material("AirOptical"); +// +// // single module +// Box mShape(mWidth/2., mWidth/2., mThick/2. - 0.1*mm); +// Volume mVol("ce_MRICH_mod_Solid", mShape, modMat); +// +// // a thin gas layer to detect optical photons +// Box modShape(mWidth/2., mWidth/2., mThick/2.); +// Volume modVol("ce_MRICH_mod_Solid_v", modShape, gasMat); +// // thin gas layer is on top (+z) of the material +// modVol.placeVolume(mVol, Position(0., 0., -0.1*mm)); +// +// modVol.setVisAttributes(description.visAttributes(mods.visStr())); +// sens.setType("photoncounter"); +// modVol.setSensitiveDetector(sens); +// +// // place modules in the sectors (disk) +// auto points = ref::utils::fillSquares({0., 0.}, mWidth + mGap, rmin - mGap, rmax + mGap); +// +// // determine module direction, always facing z = 0 +// double roty = dims.z() > 0. ? M_PI/2. : -M_PI/2.; +// int imod = 1; +// for (auto &p : points) { +// // operations are inversely ordered +// Transform3D tr = Translation3D(p.x(), p.y(), 0.) // move to position +// * RotationY(roty); // facing z = 0. +// auto modPV = mother.placeVolume(modVol, tr); +// modPV.addPhysVolID("sector", 1).addPhysVolID("module", imod ++); +// } +//} + +// clang-format off +DECLARE_DETELEMENT(athena_MRICH, createDetector) + diff --git a/src/ce_MRICH.cpp b/src/ce_MRICH.cpp deleted file mode 100644 index 5920ddc3..00000000 --- a/src/ce_MRICH.cpp +++ /dev/null @@ -1,101 +0,0 @@ - -#include <XML/Helper.h> -#include "TMath.h" -#include "TString.h" -#include "DDRec/Surface.h" -#include "DDRec/DetectorData.h" -#include "DD4hep/OpticalSurfaces.h" -#include "DD4hep/DetFactoryHelper.h" -#include "DD4hep/Printout.h" -#include "GeometryHelpers.h" - -using namespace std; -using namespace dd4hep; -using namespace dd4hep::rec; - - -void addModules(Volume &mother, xml::DetElement &detElem, Detector &desc, SensitiveDetector &sens); - - -// 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(); - // xml::Component rads = detElem.child(_Unicode(radiator)); - - auto rmin = dims.rmin(); - auto rmax = dims.rmax(); - auto length = dims.length(); - auto z0 = dims.z(); - - auto gasMat = desc.material("AirOptical"); - - // detector envelope - Tube envShape(rmin, rmax, length/2., 0., 2*M_PI); - Volume envVol("ce_MRICH_GVol", envShape, gasMat); - envVol.setVisAttributes(desc.visAttributes(detElem.visStr())); - - // modules - addModules(envVol, detElem, desc, sens); - - // place envelope - Volume motherVol = desc.pickMotherVolume(det); - PlacedVolume envPV = motherVol.placeVolume(envVol, Position(0, 0, z0)); - envPV.addPhysVolID("system", detID); - det.setPlacement(envPV); - return det; -} - - -void addModules(Volume &mother, xml::DetElement &detElem, Detector &desc, SensitiveDetector &sens) -{ - xml::Component dims = detElem.dimensions(); - xml::Component mods = detElem.child(_Unicode(modules)); - - auto rmin = dims.rmin(); - auto rmax = dims.rmax(); - - auto mThick = mods.attr<double>(_Unicode(thickness)); - auto mWidth = mods.attr<double>(_Unicode(width)); - auto mGap = mods.attr<double>(_Unicode(gap)); - - auto modMat = desc.material(mods.materialStr()); - auto gasMat = desc.material("AirOptical"); - - // single module - Box mShape(mWidth/2., mWidth/2., mThick/2. - 0.1*mm); - Volume mVol("ce_MRICH_mod_Solid", mShape, modMat); - - // a thin gas layer to detect optical photons - Box modShape(mWidth/2., mWidth/2., mThick/2.); - Volume modVol("ce_MRICH_mod_Solid_v", modShape, gasMat); - // thin gas layer is on top (+z) of the material - modVol.placeVolume(mVol, Position(0., 0., -0.1*mm)); - - modVol.setVisAttributes(desc.visAttributes(mods.visStr())); - sens.setType("photoncounter"); - modVol.setSensitiveDetector(sens); - - // place modules in the sectors (disk) - auto points = athena::geo::fillSquares({0., 0.}, mWidth + mGap, rmin - mGap, rmax + mGap); - - // determine module direction, always facing z = 0 - double roty = dims.z() > 0. ? M_PI/2. : -M_PI/2.; - int imod = 1; - for (auto &p : points) { - // operations are inversely ordered - Transform3D tr = Translation3D(p.x(), p.y(), 0.) // move to position - * RotationY(roty); // facing z = 0. - auto modPV = mother.placeVolume(modVol, tr); - modPV.addPhysVolID("sector", 1).addPhysVolID("module", imod ++); - } -} - -// clang-format off -DECLARE_DETELEMENT(refdet_ce_MRICH, createDetector) - -- GitLab