From 61e0f7954d12248476acf3db24d8f6b208a023d5 Mon Sep 17 00:00:00 2001 From: Wouter Deconinck <wouter.deconinck@umanitoba.ca> Date: Fri, 13 Aug 2021 09:06:30 +0000 Subject: [PATCH] Updated MRich with Murad's current version --- athena.xml | 2 +- compact/definitions.xml | 3 +- compact/mrich.xml | 129 +++++++++++++++++++++-------- compact/optical_materials.xml | 16 ++++ src/MRich_geo.cpp | 147 +++++++++++++++++++++++----------- 5 files changed, 217 insertions(+), 80 deletions(-) diff --git a/athena.xml b/athena.xml index a0b3aa93..9edd4748 100644 --- a/athena.xml +++ b/athena.xml @@ -173,7 +173,7 @@ ### PID detectors </documentation> <!--include ref="compact/dirc.xml"/--> - <!--include ref="compact/mrich.xml"/--> + <include ref="compact/mrich.xml"/> <!--include ref="compact/forward_trd.xml"/--> <include ref="compact/drich.xml"/> diff --git a/compact/definitions.xml b/compact/definitions.xml index 5b8de90a..865628c5 100644 --- a/compact/definitions.xml +++ b/compact/definitions.xml @@ -330,7 +330,8 @@ Examples: <constant name="ForwardTRD_length" value="25.0*cm"/> <constant name="ForwardTOF_length" value="ExtraTrackingEndcapP_length-ForwardTRD_length"/> - <constant name="BackwardCherenkov_length" value="20.0*cm"/> + <constant name="BackwardCherenkov_rmax" value="100.0*cm"/> + <constant name="BackwardCherenkov_length" value="40.0*cm"/> <constant name="ExtraTrackingEndcapN_length" value="5.0*cm"/> <constant name="BackwardTOF_length" value="0.0*cm"/> diff --git a/compact/mrich.xml b/compact/mrich.xml index e7ff797b..7a8e34b6 100644 --- a/compact/mrich.xml +++ b/compact/mrich.xml @@ -1,10 +1,11 @@ <lccdd> <comment> MRICH </comment> <define> - <constant name="MRICH_rmin" value="15*cm"/> - <constant name="MRICH_rmax" value="BarrelTracking_rmax"/> + <constant name="MRICH_rmin" value="10*cm"/> + <constant name="MRICH_rmax" value="BackwardCherenkov_rmax"/> <constant name="MRICH_length" value="BackwardCherenkov_length"/> - <constant name="MRICH_zmin" value="-BarrelTracking_length/2.0 - BackwardTracking_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"/> @@ -15,7 +16,9 @@ <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="MRICHPhotoDetCopper_thickness" value="0.1*mm"/> + <constant name="MRICHPhotoDetKapton_thickness" value="0.2*mm"/> + <constant name="MRICHRearExtraSpace_thickness" value="9.7*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"/> @@ -39,6 +42,8 @@ + MRICHFresnelLens_thickness + MRICHLensPhotoDet_length + MRICHGlassPhotoDet_thickness + + MRICHPhotoDetCopper_thickness + + MRICHPhotoDetKapton_thickness + MRICHRearExtraSpace_thickness "/> </define> @@ -54,35 +59,87 @@ <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> - --> + reflect="true" + projective="true" + vis="InvisibleWithDaughters" + material="Air"> + <dimensions rmin="MRICH_rmin" rmax="MRICH_rmax" length="abs(MRICH_length)" zmin="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"/> + <!--layer thickness="MRICHPhotoDetCopper_thickness" material="Copper"/--> + <!--layer thickness="MRICHPhotoDetKapton_thickness" material="Kapton"/--> + </photodet> + </module> + <comment> + Modules are only listed here for one quadrant + </comment> + <positions scale="1.03"> + <position x="-41.3250000000" y="178.6750000000"/> + <position x=" 96.0250000001" y="178.6750000000"/> + </positions> + <positions scale="1.036"> + <position x="-41.3250000000" y="316.0250000001"/> + <position x=" 96.0250000001" y="316.0250000001"/> + <position x="233.3750000002" y="178.6750000000"/> + <position x="233.3750000002" y="316.0250000001"/> + </positions> + <positions scale="1.052"> + <position x="-41.3250000000" y="453.3750000002"/> + <position x=" 96.0250000001" y="453.3750000002"/> + <position x="233.3750000002" y="453.3750000002"/> + <position x="370.7250000003" y="453.3750000002"/> + <position x="370.7250000003" y="316.0250000001"/> + <position x="370.7250000003" y="178.6750000000"/> + </positions> + <positions scale="1.078"> + <position x="-41.3250000000" y="590.7250000003"/> + <position x=" 96.0250000001" y="590.7250000003"/> + <position x="233.3750000002" y="590.7250000003"/> + <position x="370.7250000003" y="590.7250000003"/> + <position x="508.0750000004" y="590.7250000003"/> + <position x="508.0750000004" y="178.6750000000"/> + <position x="508.0750000004" y="316.0250000001"/> + <position x="508.0750000004" y="453.3750000002"/> + </positions> + <positions scale="1.09"> + <position x="-41.3250000000" y="728.0750000004"/> + <position x=" 96.0250000001" y="728.0750000004"/> + <position x="233.3750000002" y="728.0750000004"/> + <!--position x="370.7250000003" y="728.0750000004"/--> + <!--position x="508.0750000004" y="728.0750000004"/--> + <!--position x="645.4250000005" y="590.7250000003"/--> + <position x="645.4250000005" y="453.3750000002"/> + <position x="645.4250000005" y="316.0250000001"/> + <position x="645.4250000005" y="178.6750000000"/> + </positions> + <positions scale="1.095"> + <!--position x="-41.3250000000" y="865.4250000005"/--> + <!--position x=" 96.0250000001" y="865.4250000005"/--> + <!--position x="233.3750000002" y="865.4250000005"/--> + <position x="782.7750000006" y="178.6750000000"/> + <!--position x="782.7750000006" y="316.0250000001"/--> + <!--position x="782.7750000006" y="453.3750000002"/--> + </positions> </detector> </detectors> @@ -92,5 +149,11 @@ <id>system:8,module:14,sensor:8,x:32:-16,y:-16</id> </readout> </readouts> + + <!--Globals> + <Parameter Name="mrichInfo" Value="mrichmod/mrich_1_geoparams-0-0-4294967295-1527211159.xml"/> + </Globals--> + + </lccdd> diff --git a/compact/optical_materials.xml b/compact/optical_materials.xml index c510b815..242add12 100644 --- a/compact/optical_materials.xml +++ b/compact/optical_materials.xml @@ -47,6 +47,21 @@ 1240*eV/600 1.49 1240*eV/400 1.49 "/> + <matrix name="ABSLENGTH__Acrylic" coldim="2" values=" + 1240*eV/1100 2631*mm + 1240*eV/1000 2631*mm + 1240*eV/900 2631*mm + 1240*eV/800 2631*mm + 1240*eV/700 2500*mm + 1240*eV/600 2272*mm + 1240*eV/500 2000*mm + 1240*eV/400 1315*mm + 1240*eV/300 1613*mm + 1240*eV/250 740*mm + 1240*eV/225 125*mm + 1240*eV/210 10*mm + 1240*eV/200 0*mm + "/> <!-- BEGIN dRICh material properties - dumped from fun4all implementation - see https://github.com/cisbani/dRICh/blob/main/share/source/g4dRIChOptics.hh @@ -377,6 +392,7 @@ <composite n="2" ref="O"/> <composite n="8" ref="H"/> <property name="RINDEX" ref="RINDEX__Acrylic"/> + <property name="ABSLENGTH" ref="ABSLENGTH__Acrylic"/> </material> <!-- BEGIN dRICh material definitions --> <material name="C2F6_DRICH"> diff --git a/src/MRich_geo.cpp b/src/MRich_geo.cpp index 33d1059c..914e8fda 100644 --- a/src/MRich_geo.cpp +++ b/src/MRich_geo.cpp @@ -27,13 +27,25 @@ static Ref_t createDetector(Detector& description, xml::Handle_t e, SensitiveDet 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(); + // read module positions + std::vector<std::pair<double,double>> positions; + for (xml_coll_t x_positions_i(x_det, _Unicode(positions)); x_positions_i; ++x_positions_i) { + xml_comp_t x_positions = x_positions_i; + for (xml_coll_t x_position_i(x_positions, _U(position)); x_position_i; ++x_position_i) { + xml_comp_t x_position = x_position_i; + positions.push_back( + std::make_pair(x_positions.scale() * x_position.x() * mm, + x_positions.scale() * x_position.y() * mm)); + } + } + bool projective = getAttrOrDefault(x_det, _Unicode(projective), false); + bool reflect = x_det.reflect(true); PlacedVolume pv; @@ -49,6 +61,7 @@ static Ref_t createDetector(Detector& description, xml::Handle_t e, SensitiveDet auto rmax = dims.rmax(); auto length = dims.length(); auto zmin = dims.zmin(); + auto zpos = zmin + length / 2; // expect only one module (for now) xml_comp_t x_mod = x_det.child(_U(module)); @@ -110,14 +123,14 @@ static Ref_t createDetector(Detector& description, xml::Handle_t e, SensitiveDet // - 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 groove_pitch = getAttrOrDefault(x_lens, _Unicode(pitch), 0.2 * mm);// 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 eff_diameter = getAttrOrDefault(x_lens, _Unicode(effective_diameter), 152.4 * 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 center_thickness = getAttrOrDefault(x_lens, _U(thickness), 0.068 * 2.54 * cm);//2.0 * mm); double n_acrylic = 1.49; - double lens_curvature = 1.0 / (lens_f*(n_acrylic - 1.0)); + double lens_curvature = 1.0 / (lens_f*(n_acrylic - 1.0)); //confirmed double full_ring_rmax = std::min(eff_diameter / 2.0, lens_width/2.0); double N_grooves = std::ceil((full_ring_rmax) / groove_pitch); @@ -143,15 +156,12 @@ static Ref_t createDetector(Detector& description, xml::Handle_t e, SensitiveDet 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)); + lens_vol.placeVolume(lens_groove_vol); //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 @@ -173,6 +183,7 @@ static Ref_t createDetector(Detector& description, xml::Handle_t e, SensitiveDet 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); @@ -221,12 +232,32 @@ static Ref_t createDetector(Detector& description, xml::Handle_t e, SensitiveDet 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); + // sensitive pv.addPhysVolID("sensor", n_sensor); + window_vol.setSensitiveDetector(sens); + sensitives[mod_name].push_back(pv); + ++n_sensor; + + // photon detector electronics layers + double layer_zpos = window_zpos + photodet_thickness/2.0; + int i_layer = 1; + for (xml_coll_t li(x_photodet, _Unicode(layer)); li; ++li) { + xml_comp_t x_layer = li; + Material layer_mat = description.material(x_layer.materialStr()); + double layer_thickness = x_layer.thickness(); + Box layer_box(photodet_width/2.0,photodet_width/2.0,layer_thickness/2.0); + Volume layer_vol(mod_name + "_layer_" + std::to_string(i_layer), layer_box, layer_mat); + layer_zpos += layer_thickness / 2.0; + pv = m_volume.placeVolume(layer_vol,Position(0,0,layer_zpos)); + DetElement layer_de(mod_de, std::string("mod_layer_de_") + std::to_string(i_layer), 1); + layer_de.setPlacement(pv); + layer_zpos += layer_thickness / 2.0; + i_layer++; + } //for (size_t ic = 0; ic < sensVols.size(); ++ic) { // PlacedVolume sens_pv = sensVols[ic]; @@ -242,14 +273,10 @@ static Ref_t createDetector(Detector& description, xml::Handle_t e, SensitiveDet //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); @@ -263,41 +290,71 @@ static Ref_t createDetector(Detector& description, xml::Handle_t e, SensitiveDet 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. + + // place modules + int i_mod = 1; // starts at 1 + for (auto& p: positions) { + + // get positions in one quadrant + double x = p.first; + double y = p.second; + double z = -zpos; + + // and place in all quadrants (intentional shadowing) + for (auto& p: decltype(positions){{x,y}, {y,-x}, {-x,-y}, {-y,x}}) { + + // get positions (intentional shadowing) + double x = p.first; + double y = p.second; + + // get angles + double rotAngX = atan(y/z); + double rotAngY = -1.*atan(x/z); + + /* + 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. + */ + + Transform3D tr; + if(projective) { + tr = Translation3D(x, y, 0) + * RotationX(rotAngX) + * RotationY(rotAngY); + } else { + tr = Translation3D(x, y, 0) * RotationX(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++; } - // 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); + Volume motherVol = description.pickMotherVolume(sdet); + if (reflect) { + pv = motherVol.placeVolume(envVol, Transform3D(RotationZYX(0, M_PI, 0), Position(0, 0, -zpos))); + } else { + pv = motherVol.placeVolume(envVol, Transform3D(RotationZYX(0, 0, 0), Position(0, 0, +zpos))); + } + pv.addPhysVolID("system", x_det.id()); + sdet.setPlacement(pv); return sdet; } -- GitLab