Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

dRICH large sensor for focal point region studies #351

Open
wants to merge 24 commits into
base: main
Choose a base branch
from
Open

dRICH large sensor for focal point region studies #351

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
24 commits
Select commit Hold shift + click to select a range
646fb74
Added oversized sensor for dRICH, to catch all reflected photons (for…
cpecar Jan 13, 2023
b29f29c
2 new debug optics modes: one to create only large sensor for one sec…
cpecar Jan 13, 2023
5d388f3
[pre-commit.ci] auto fixes from pre-commit.com hooks
pre-commit-ci[bot] Jan 13, 2023
8585073
fix: `airgapZpos` and `filterZpos` scoping
c-dilks Jan 13, 2023
d7e99e9
style
c-dilks Jan 13, 2023
3dd9fd7
fix: add Connor to copyright headers
c-dilks Jan 13, 2023
1cefaac
Merge remote-tracking branch 'origin/main' into dRICHlargeSensor
c-dilks Jan 25, 2023
6e233ac
Changed reading in of focal point positions to ifstream. Added variab…
cpecar Jan 27, 2023
bc4a944
changed variable type to string
cpecar Jan 27, 2023
e0df142
reset to default debug values
cpecar Jan 27, 2023
d7f8567
[pre-commit.ci] auto fixes from pre-commit.com hooks
pre-commit-ci[bot] Jan 27, 2023
68f9484
Update compact/pid/drich.xml
cpecar Jan 27, 2023
96e0efe
Update src/DRICH_geo.cpp
cpecar Jan 27, 2023
f17ea63
Update src/DRICH_geo.cpp
cpecar Jan 27, 2023
8982518
Changed from C style cast to static_cast
cpecar Jan 27, 2023
64c8e56
combining case statements
cpecar Jan 27, 2023
806ce52
debugOptics changed to debugOpticsMode (no int to bool conversion)
cpecar Jan 27, 2023
52d6306
ifstream to skip txt comments starting with '#'
cpecar Jan 31, 2023
99e0d2e
[pre-commit.ci] auto fixes from pre-commit.com hooks
pre-commit-ci[bot] Jan 31, 2023
ab00025
Merge branch 'main' into dRICHlargeSensor
cpecar May 10, 2023
6d3f5f7
fix inconsistency after merging commits from main
cpecar May 10, 2023
1493ac6
Merge branch 'main' into dRICHlargeSensor
c-dilks Jun 6, 2023
1838608
Update src/DRICH_geo.cpp
cpecar Jun 15, 2023
16b99be
Update src/DRICH_geo.cpp
cpecar Jun 15, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 4 additions & 1 deletion compact/pid/drich.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
<!-- SPDX-License-Identifier: LGPL-3.0-or-later -->
<!-- Copyright (C) 2023 Christopher Dilks, Chandradoy Chatterjee -->
<!-- Copyright (C) 2023 Christopher Dilks, Chandradoy Chatterjee, Connor Pecar -->

<lccdd>

Expand Down Expand Up @@ -37,6 +37,8 @@
- `DRICH_debug_optics`: 1 = all components become vacuum, except for mirrors; test opticalphotons from IP
2 = all components become vacuum, except for mirrors and `gasvol`, test charged particles from IP
3 = all components become vacuum, except for mirrors and sensors, test reflected trajectories' hits
4 = all components become vacuum, except for mirrors and sensors, oversized sensor created to catch all photons
5 = all components created normally, calculated focal points from mode 4 drawn if environment variable set
0 = off
- `DRICH_debug_sector`: 1 = only include one sector (will be set automatically when `DRICH_debug_optics>0`)
- `DRICH_debug_mirror`: 1 = draw full mirror shape for single sector; 0 = off
Expand All @@ -46,6 +48,7 @@
<constant name="DRICH_debug_sector" value="0"/>
<constant name="DRICH_debug_mirror" value="0"/>
<constant name="DRICH_debug_sensors" value="0"/>
<constant name="DRICH_FP_file" value="0" type="string"/>
</define>


Expand Down
232 changes: 159 additions & 73 deletions src/DRICH_geo.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
// SPDX-License-Identifier: LGPL-3.0-or-later
// Copyright (C) 2022 Christopher Dilks, Junhuai Xu
// Copyright (C) 2022, 2023 Christopher Dilks, Junhuai Xu, Connor Pecar

//==========================================================================
// dRICH: Dual Ring Imaging Cherenkov Detector
Expand All @@ -17,6 +17,7 @@
#include "DD4hep/Printout.h"
#include "DDRec/DetectorData.h"
#include "DDRec/Surface.h"
#include "DD4hep/CartesianGridXY.h"

#include <XML/Helper.h>

Expand Down Expand Up @@ -110,10 +111,10 @@ static Ref_t createDetector(Detector& desc, xml::Handle_t handle, SensitiveDetec
bool debugSector = desc.constantAsLong("DRICH_debug_sector") == 1;
bool debugMirror = desc.constantAsLong("DRICH_debug_mirror") == 1;
bool debugSensors = desc.constantAsLong("DRICH_debug_sensors") == 1;
std::string FPfile = desc.constantAsString("DRICH_FP_file");

// if debugging optics, override some settings
bool debugOptics = debugOpticsMode > 0;
if (debugOptics) {
if (debugOpticsMode != 0) {
printout(WARNING, "DRICH_geo", "DEBUGGING DRICH OPTICS");
switch (debugOpticsMode) {
case 1:
Expand All @@ -123,6 +124,8 @@ static Ref_t createDetector(Detector& desc, xml::Handle_t handle, SensitiveDetec
vesselMat = aerogelMat = filterMat = sensorMat = desc.material("VacuumOptical");
break;
case 3:
case 4:
case 5:
vesselMat = aerogelMat = filterMat = gasvolMat = desc.material("VacuumOptical");
break;
default:
Expand All @@ -131,8 +134,19 @@ static Ref_t createDetector(Detector& desc, xml::Handle_t handle, SensitiveDetec
}
}

// sensor size rescaling for large sensor (catching all photons)
double sensorRescale = 0;
if (debugOpticsMode == 4) sensorRescale = 250;
if (sensorRescale > 0) {
auto seg = static_cast<CartesianGridXY>( desc.readout(readoutName).segmentation() );
seg.setGridSizeX(sensorRescale * seg.gridSizeX());
seg.setGridSizeY(sensorRescale * seg.gridSizeY());
seg.setOffsetX(sensorRescale * seg.offsetX());
seg.setOffsetY(sensorRescale * seg.offsetY());
}

// if debugging anything, draw only one sector and adjust visibility
if (debugOptics || debugMirror || debugSensors)
if (debugOpticsMode!=0 || debugMirror || debugSensors)
debugSector = true;
if (debugSector)
gasvolVis = vesselVis = desc.invisible();
Expand Down Expand Up @@ -206,6 +220,8 @@ static Ref_t createDetector(Detector& desc, xml::Handle_t handle, SensitiveDetec
break; // `!debugOptics`
case 1:
case 3:
case 4:
case 5:
vesselSolid = vesselBox;
gasvolSolid = gasvolBox;
break;
Expand Down Expand Up @@ -278,7 +294,7 @@ static Ref_t createDetector(Detector& desc, xml::Handle_t handle, SensitiveDetec
// aerogelSkin.isValid();

// airgap and filter placement and surface properties
if (!debugOptics) {
if (debugOpticsMode == 0) {

auto airgapPlacement =
Translation3D(radiatorPos.x(), radiatorPos.y(), radiatorPos.z()) * // re-center to originFront
Expand Down Expand Up @@ -314,6 +330,7 @@ static Ref_t createDetector(Detector& desc, xml::Handle_t handle, SensitiveDetec
desc.add(Constant("DRICH_filter_material", filterMat.ptr()->GetName(), "string"));
desc.add(Constant("DRICH_gasvol_material", gasvolMat.ptr()->GetName(), "string"));


// SECTOR LOOP //////////////////////////////////////////////////////////////////////

// initialize sensor centroids (used for mirror parameterization below); this is
Expand All @@ -324,7 +341,6 @@ static Ref_t createDetector(Detector& desc, xml::Handle_t handle, SensitiveDetec
// int sensorCount = 0;

for (int isec = 0; isec < nSectors; isec++) {

// debugging filters, limiting the number of sectors
if (debugSector && isec != 0)
continue;
Expand Down Expand Up @@ -423,7 +439,7 @@ static Ref_t createDetector(Detector& desc, xml::Handle_t handle, SensitiveDetec
auto sensorSphPos = Position(sensorSphCenterX, 0., sensorSphCenterZ) + originFront;

// sensitivity
if (!debugOptics || debugOpticsMode == 3)
if (debugOpticsMode == 0 || debugOpticsMode == 3 || debugOpticsMode == 4)
sensorVol.setSensitiveDetector(sens);

// reconstruction constants
Expand Down Expand Up @@ -459,84 +475,154 @@ static Ref_t createDetector(Detector& desc, xml::Handle_t handle, SensitiveDetec

// 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);

// shift global coordinates so we can apply spherical patch cuts
double zCheck = z + sensorSphCenterZ;
double xCheck = x + sensorSphCenterX;
double yCheck = y;
double rCheck = std::hypot(xCheck, yCheck);
double phiCheck = std::atan2(yCheck, xCheck);

// patch cut
bool patchCut = std::fabs(phiCheck) < sensorSphPatchPhiw && zCheck > sensorSphPatchZmin &&
rCheck > sensorSphPatchRmin && rCheck < sensorSphPatchRmax;
if (debugSensors)
patchCut = std::fabs(phiCheck) < sensorSphPatchPhiw;
if (patchCut) {

// append sensor position to centroid calculation
// if (isec == 0) {
// sensorCentroidX += xCheck;
// sensorCentroidZ += zCheck;
// sensorCount++;
// }

// placement (note: transformations are in reverse order)
// - transformations operate on global coordinates; the corresponding
// generator coordinates are provided in the comments
auto sensorPlacement =
if(debugOpticsMode != 4) {
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);

// shift global coordinates so we can apply spherical patch cuts
double zCheck = z + sensorSphCenterZ;
double xCheck = x + sensorSphCenterX;
double yCheck = y;
double rCheck = std::hypot(xCheck, yCheck);
double phiCheck = std::atan2(yCheck, xCheck);

// patch cut
bool patchCut = std::fabs(phiCheck) < sensorSphPatchPhiw && zCheck > sensorSphPatchZmin &&
rCheck > sensorSphPatchRmin && rCheck < sensorSphPatchRmax;
if (debugSensors)
patchCut = std::fabs(phiCheck) < sensorSphPatchPhiw;
if (patchCut) {

// append sensor position to centroid calculation
// if (isec == 0) {
// sensorCentroidX += xCheck;
// sensorCentroidZ += zCheck;
// sensorCount++;
// }

// placement (note: transformations are in reverse order)
// - transformations operate on global coordinates; the corresponding
// generator coordinates are provided in the comments
auto sensorPlacement =
sectorRotation * // rotate about beam axis to sector
Translation3D(sensorSphPos.x(), sensorSphPos.y(), sensorSphPos.z()) * // move sphere to reference position
RotationX(phiGen) * // rotate about `zGen`
RotationZ(thetaGen) * // rotate about `yGen`
Translation3D(-sensorThickness / 2.0, 0.,
0.) * // pull back so sensor active surface is at spherical surface
0.) * // pull back so sensor active surface is at spherical surface
Translation3D(sensorSphRadius, 0., 0.) * // push radially to spherical surface
RotationY(M_PI / 2) * // rotate sensor to be compatible with generator coords
RotationZ(-M_PI / 2); // correction for readout segmentation mapping
auto sensorPV = gasvolVol.placeVolume(sensorVol, sensorPlacement);
auto sensorPV = gasvolVol.placeVolume(sensorVol, sensorPlacement);

// generate LUT for module number -> sensor position, for readout mapping tests
// if(isec==0) printf("%d %f %f\n",imod,sensorPV.position().x(),sensorPV.position().y());
// generate LUT for module number -> sensor position, for readout mapping tests
// if(isec==0) printf("%d %f %f\n",imod,sensorPV.position().x(),sensorPV.position().y());

// properties
sensorPV.addPhysVolID("sector", isec)
// properties
sensorPV.addPhysVolID("sector", isec)
.addPhysVolID("module", imod); // NOTE: must be consistent with `sensorIDfields`
auto imodsec = encodeSensorID(sensorPV.volIDs());
std::string modsecName = secName + "_" + std::to_string(imod);
DetElement sensorDE(det, "sensor_de_" + modsecName, imodsec);
sensorDE.setPlacement(sensorPV);
if (!debugOptics || debugOpticsMode == 3) {
SkinSurface sensorSkin(desc, sensorDE, "sensor_optical_surface_" + modsecName, sensorSurf, sensorVol);
sensorSkin.isValid();
auto imodsec = encodeSensorID(sensorPV.volIDs());
std::string modsecName = secName + "_" + std::to_string(imod);
DetElement sensorDE(det, "sensor_de_" + modsecName, imodsec);
sensorDE.setPlacement(sensorPV);
if (debugOpticsMode == 0 || debugOpticsMode == 3 || debugOpticsMode == 5) {
SkinSurface sensorSkin(desc, sensorDE, "sensor_optical_surface_" + modsecName, sensorSurf, sensorVol);
sensorSkin.isValid();
}

// increment sensor module number
imod++;

} // end patch cuts
} // end phiGen loop
} // end thetaGen loop
} // end if(debugOpticsMode != 4)

// large sensor placement (for focal point testing):
if (debugOpticsMode == 4){
Comment on lines +557 to +560
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

else ?

Box sensorSolidScaled(sensorRescale * sensorSide / 2., sensorRescale * sensorSide / 2., sensorThickness / 2.);
Volume sensorVolScaled(detName + "_sensor_" + secName, sensorSolidScaled, sensorMat);
sensorVolScaled.setVisAttributes(sensorVis);
sensorVolScaled.setSensitiveDetector(sens);

double phiGen = 0;
auto sensorPlacement =
RotationZ(sectorRotation) *
Translation3D(sensorSphPos.x()+100, sensorSphPos.y(), sensorSphPos.z()-125) *
RotationX(phiGen) *
Translation3D(sensorSphRadius, 0., 0.)*
RotationZ(0);
auto sensorPV = gasvolVol.placeVolume(sensorVolScaled, sensorPlacement);
sensorPV.addPhysVolID("sector", isec).addPhysVolID("module", imod);

DetElement sensorDE(det, Form("sensor_de%d_%d", isec, imod), 0);
sensorDE.setPlacement(sensorPV);

SkinSurface sensorSkin(desc, sensorDE, Form("sensor_optical_surface%d", isec), sensorSurf, sensorVolScaled);
sensorSkin.isValid();
} // end large sensor placement

// for debugOpticsMode 5: drawing calculated focal points
// of thrown beams of parallel photons
if( FPfile!="0" && debugOpticsMode == 5){
Cone FPsolid(1.,0.25,0.5,1,1);
Volume FPvol(detName + "_FPpos_" + secName, FPsolid, aerogelMat);

std::ifstream fptxt(FPfile);
std::string fpstr;
double fpx, fpy, fpz, dirx, diry, dirz;
int fpnum=0;

while(std::getline(fptxt,fpstr)){
if(fpstr[0]=='#') continue;
std::istringstream stringStream(fpstr);
stringStream >> fpx >> fpy >> fpz >> dirx >> diry >> dirz;
if( std::abs(fpx) < 1000 && std::abs(fpy) < 1000 && std::abs(fpz) < 1000){
double zrot = 0;
if(dirx < 0){
if(diry < 0) zrot = -1*std::atan2(dirx,diry) + M_PI;
else zrot = -1*std::atan2(dirx,diry);
}

// increment sensor module number
imod++;

} // end patch cuts
} // end phiGen loop
} // end thetaGen loop
else{
if(diry < 0) zrot = -1*std::atan2(dirx,diry) - M_PI;
else zrot = std::atan2(dirx,diry) + M_PI;
}
double xrot = std::atan2(sqrt(diry*diry+dirx*dirx),abs(dirz));
xrot = 0;
zrot = 0;
auto FPPlacement =
RotationZ(0) *
Translation3D(fpx, fpy, fpz-vesselPos.z())*
RotationX(xrot)*
RotationZ(zrot)
;
auto FPPV = gasvolVol.placeVolume(FPvol, FPPlacement);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Variable names which are all-caps are considered reserved for preprocessor definitions, absolute constants, and environment variable names.


FPPV.addPhysVolID("sector", isec).addPhysVolID("module", imod);
DetElement FPDE(det, Form("FP_de%d_%d", isec, fpnum), 0);
FPDE.setPlacement(FPPV);
fpnum++;
}
}
}

// calculate centroid sensor position
// if (isec == 0) {
Expand Down