diff --git a/CMakeLists.txt b/CMakeLists.txt
index 716afdc684a18e00d5b2afab6532ab4a93122488..f1e31b27a393edaad08d6662c8b60f89de466794 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -54,6 +54,12 @@ else()
endif()
find_package(fmt REQUIRED)
+if(DEFINED IRT_AUXFILE)
+ add_compile_definitions(IRT_AUXFILE)
+ find_package(IRT REQUIRED)
+ include_directories( ${CMAKE_INSTALL_PREFIX}/include ${CMAKE_INSTALL_PREFIX}/include/IRT)
+endif()
+
#-----------------------------------------------------------------------------------
set(a_lib_name ${PROJECT_NAME})
@@ -67,6 +73,9 @@ dd4hep_add_plugin(${a_lib_name}
target_link_libraries(${a_lib_name}
PUBLIC ${DD4hep_required_libraries} fmt::fmt
)
+if(DEFINED IRT_AUXFILE)
+ target_link_libraries(${a_lib_name} PUBLIC IRT)
+endif()
#-----------------------------------------------------------------------------------
# Parse jinja templates: once by default, and once for all yml files
diff --git a/compact/drich.xml b/compact/drich.xml
index c98d41301f7d33363cff28faf98d432a7b341b09..2cb26423d7747558e3a2c2c2a9969dcd6101b55f 100644
--- a/compact/drich.xml
+++ b/compact/drich.xml
@@ -32,10 +32,13 @@
0 = off
- `DRICH_debug_mirror`: 1 = draw full mirror shape for single sector; 0 = off
- `DRICH_debug_sensors`: 1 = draw full sensor sphere for a single sector; 0 = off
+- `DRICH_create_irt_file`: 1 = create an auxiliary IRT config ROOT file; 0 = off
+ - the file name is set by detector attribute `irt_filename` below
</comment>
<constant name="DRICH_debug_optics" value="0"/>
<constant name="DRICH_debug_mirror" value="0"/>
<constant name="DRICH_debug_sensors" value="0"/>
+<constant name="DRICH_create_irt_file" value="0"/>
</define>
@@ -55,6 +58,7 @@
material="Aluminum"
vis_vessel="DRICH_vessel_vis"
vis_gas="DRICH_gas_vis"
+ irt_filename="irt-drich.root"
>
diff --git a/scripts/create_IRT_auxfile.py b/scripts/create_IRT_auxfile.py
new file mode 100755
index 0000000000000000000000000000000000000000..1f0b79a3dd03b466dd0d1839678b245b4d967656
--- /dev/null
+++ b/scripts/create_IRT_auxfile.py
@@ -0,0 +1,54 @@
+# produce auxiliary ROOT file for the IRT algorithm
+
+import shutil, os, sys, argparse
+import xml.etree.ElementTree as et
+
+parser = argparse.ArgumentParser()
+parser.add_argument(
+ '-o', '--output-name', dest='outFile', default='irt-drich.root',
+ help='output auxiliary ROOT file name', type=str
+ )
+argv = parser.parse_args()
+
+import DDG4
+
+
+def run():
+
+ # compact files
+ if not 'DETECTOR_PATH' in os.environ:
+ print('ERROR: env var DETECTOR_PATH not set',file=sys.stderr)
+ exit(1)
+ mainFile = os.environ['DETECTOR_PATH'] + '/ecce.xml'
+ richFile = os.environ['DETECTOR_PATH'] + '/compact/drich.xml'
+
+ # backup original richFile, then parse
+ shutil.copy(richFile,richFile+'.bak')
+ richTree = et.parse(richFile)
+
+ # enable `DRICH_create_irt_file` mode
+ for constant in richTree.iter(tag='constant'):
+ if(constant.attrib['name']=='DRICH_create_irt_file'):
+ constant.set('value','1')
+
+ # set auxiliary file name
+ for detector in richTree.iter(tag='detector'):
+ detector.set('irt_filename',argv.outFile)
+
+ # overwrite original richFile
+ richTree.write(richFile)
+
+ # produce IRT config file
+ try:
+ kernel = DDG4.Kernel()
+ kernel.loadGeometry(f'file:{mainFile}')
+ kernel.terminate()
+ print(f'\n -> produced {argv.outFile}\n')
+ except:
+ pass
+
+ # revert to the original richFile
+ os.replace(richFile+'.bak',richFile)
+
+if __name__ == "__main__":
+ run()
diff --git a/src/DRICH_geo.cpp b/src/DRICH_geo.cpp
index 1a2bd4f2a616b71511ba4b2141c58e763fc9d672..0dfaf46fabe4520deadcf98a1d6999a9c744bb52 100644
--- a/src/DRICH_geo.cpp
+++ b/src/DRICH_geo.cpp
@@ -17,6 +17,16 @@
#include <XML/Helper.h>
+#ifdef IRT_AUXFILE
+#include <CherenkovDetectorCollection.h>
+#include <CherenkovPhotonDetector.h>
+#include <CherenkovRadiator.h>
+#include <OpticalBoundary.h>
+#include <ParametricSurface.h>
+
+#include <TFile.h>
+#endif
+
using namespace dd4hep;
using namespace dd4hep::rec;
@@ -101,6 +111,11 @@ static Ref_t createDetector(Detector& desc, xml::Handle_t handle, SensitiveDetec
long debugOpticsMode = desc.constantAsLong("DRICH_debug_optics");
bool debugMirror = desc.constantAsLong("DRICH_debug_mirror") == 1;
bool debugSensors = desc.constantAsLong("DRICH_debug_sensors") == 1;
+#ifdef IRT_AUXFILE
+ // - IRT auxiliary file
+ auto irtAuxFileName = detElem.attr<std::string>(_Unicode(irt_filename));
+ bool createIrtFile = desc.constantAsLong("DRICH_create_irt_file") == 1;
+#endif
// if debugging optics, override some settings
bool debugOptics = debugOpticsMode > 0;
@@ -132,6 +147,58 @@ static Ref_t createDetector(Detector& desc, xml::Handle_t handle, SensitiveDetec
printout(DEBUG, "DRICH_geo", "sectorMask, sectorOffset, moduleMask, moduleOffset = 0x%x %d 0x%x %d",
sectorBits.mask(), sectorBits.offset(), moduleBits.mask(), moduleBits.offset());
+#ifdef IRT_AUXFILE
+ // IRT geometry auxiliary file ===========================================================
+ /* - optionally generate an auxiliary ROOT file, storing geometry objects for IRT
+ * - use compact file variable `DRICH_create_irt_file` to control this
+ */
+ TFile* irtAuxFile = nullptr;
+ CherenkovDetectorCollection* irtGeometry;
+ CherenkovDetector* irtDetector;
+ if (createIrtFile) {
+ irtAuxFile = new TFile(irtAuxFileName.c_str(), "RECREATE");
+ printout(ALWAYS, "IRTLOG", "Producing auxiliary ROOT file for IRT: %s", irtAuxFileName.c_str());
+ irtGeometry = new CherenkovDetectorCollection();
+ irtDetector = irtGeometry->AddNewDetector(detName.c_str());
+ }
+
+ // container volume (envelope?)
+ /* FIXME: have no connection to GEANT G4LogicalVolume pointers; however all is needed
+ * is to make them unique so that std::map work internally; resort to using integers,
+ * who cares; material pointer can seemingly be '0', and effective refractive index
+ * for all radiators will be assigned at the end by hand; FIXME: should assign it on
+ * per-photon basis, at birth, like standalone GEANT code does;
+ */
+ FlatSurface* irtBoundary;
+ TVector3 normX(1, 0, 0); // normal vectors
+ TVector3 normY(0, -1, 0);
+ if (createIrtFile) {
+ irtBoundary = new FlatSurface((1 / mm) * TVector3(0, 0, vesselZmin), normX, normY);
+ for (int isec = 0; isec < nSectors; isec++) {
+ auto rad = irtGeometry->SetContainerVolume(
+ irtDetector, // Cherenkov detector
+ "GasVolume", // name
+ isec, // path
+ (G4LogicalVolume*)(0x0), // G4LogicalVolume (inaccessible? use an integer instead)
+ nullptr, // G4RadiatorMaterial (inaccessible?)
+ irtBoundary // surface
+ );
+ rad->SetAlternativeMaterialName(gasvolMat.ptr()->GetName());
+ }
+ }
+
+ // photon detector // FIXME: args (G4Solid,G4Material) inaccessible?
+ CherenkovPhotonDetector* irtPhotonDetector;
+ if (createIrtFile) {
+ irtDetector->SetReadoutCellMask(cellMask); // readout mask
+ irtPhotonDetector = new CherenkovPhotonDetector(nullptr, nullptr);
+ irtGeometry->AddPhotonDetector(irtDetector, // Cherenkov detector
+ nullptr, // G4LogicalVolume (inaccessible?)
+ irtPhotonDetector // photon detector
+ );
+ }
+#endif
+
// BUILD VESSEL ====================================================================
/* - `vessel`: aluminum enclosure, the mother volume of the dRICh
* - `gasvol`: gas volume, which fills `vessel`; all other volumes defined below
@@ -265,6 +332,47 @@ static Ref_t createDetector(Detector& desc, xml::Handle_t handle, SensitiveDetec
// filterSkin.isValid();
};
+#ifdef IRT_AUXFILE
+ // IRT aerogel + filter
+ /* AddFlatRadiator will create a pair of flat refractive surfaces internally;
+ * FIXME: should make a small gas gap at the upstream end of the gas volume;
+ */
+ FlatSurface* aerogelFlatSurface;
+ FlatSurface* filterFlatSurface;
+ if (createIrtFile) {
+ double irtAerogelZpos = vesselPos.z() + aerogelPV.position().z(); // center of aerogel, w.r.t. IP
+ double irtFilterZpos = vesselPos.z() + filterPV.position().z(); // center of filter, w.r.t. IP
+ aerogelFlatSurface = new FlatSurface((1 / mm) * TVector3(0, 0, irtAerogelZpos), normX, normY);
+ filterFlatSurface = new FlatSurface((1 / mm) * TVector3(0, 0, irtFilterZpos), normX, normY);
+ for (int isec = 0; isec < nSectors; isec++) {
+ auto aerogelFlatRadiator = irtGeometry->AddFlatRadiator(
+ irtDetector, // Cherenkov detector
+ "Aerogel", // name
+ isec, // path
+ (G4LogicalVolume*)(0x1), // G4LogicalVolume (inaccessible? use an integer instead)
+ nullptr, // G4RadiatorMaterial
+ aerogelFlatSurface, // surface
+ aerogelThickness / mm // surface thickness
+ );
+ auto filterFlatRadiator = irtGeometry->AddFlatRadiator(
+ irtDetector, // Cherenkov detector
+ "Filter", // name
+ isec, // path
+ (G4LogicalVolume*)(0x2), // G4LogicalVolume (inaccessible? use an integer instead)
+ nullptr, // G4RadiatorMaterial
+ filterFlatSurface, // surface
+ filterThickness / mm // surface thickness
+ );
+ aerogelFlatRadiator->SetAlternativeMaterialName(aerogelMat.ptr()->GetName());
+ filterFlatRadiator->SetAlternativeMaterialName(filterMat.ptr()->GetName());
+ }
+ printout(ALWAYS, "IRTLOG", "irtAerogelZpos = %f cm", irtAerogelZpos);
+ printout(ALWAYS, "IRTLOG", "irtFilterZpos = %f cm", irtFilterZpos);
+ printout(ALWAYS, "IRTLOG", "aerogel thickness = %f cm", aerogelThickness);
+ printout(ALWAYS, "IRTLOG", "filter thickness = %f cm", filterThickness);
+ }
+#endif
+
// SECTOR LOOP //////////////////////////////////////////////////////////////////////
// initialize sensor centroids (used for mirror parameterization below); this is
@@ -397,6 +505,45 @@ static Ref_t createDetector(Detector& desc, xml::Handle_t handle, SensitiveDetec
sensorSkin.isValid();
};
+#ifdef IRT_AUXFILE
+ // IRT sensors
+ FlatSurface* sensorFlatSurface;
+ if (createIrtFile) {
+ // get sensor position, w.r.t. IP
+ // - sensorGlobalPos X and Y are equivalent to sensorPV.position()
+ double sensorLocalPos[3] = {0.0, 0.0, 0.0};
+ double sensorBufferPos[3];
+ double sensorGlobalPos[3];
+ sensorPV.ptr()->LocalToMaster(sensorLocalPos, sensorBufferPos);
+ vesselPV.ptr()->LocalToMaster(sensorBufferPos, sensorGlobalPos);
+
+ // get sensor flat surface normX and normY
+ // - ignore vessel transformation, since it is a pure translation
+ double sensorLocalNormX[3] = {1.0, 0.0, 0.0};
+ double sensorLocalNormY[3] = {0.0, 1.0, 0.0};
+ double sensorGlobalNormX[3], sensorGlobalNormY[3];
+ sensorPV.ptr()->LocalToMasterVect(sensorLocalNormX, sensorGlobalNormX);
+ sensorPV.ptr()->LocalToMasterVect(sensorLocalNormY, sensorGlobalNormY);
+
+ // create the IRT sensor geometry
+ sensorFlatSurface = new FlatSurface((1 / mm) * TVector3(sensorGlobalPos), TVector3(sensorGlobalNormX),
+ TVector3(sensorGlobalNormY));
+ irtDetector->CreatePhotonDetectorInstance(isec, // sector
+ irtPhotonDetector, // CherenkovPhotonDetector
+ imodsec, // copy number
+ sensorFlatSurface // surface
+ );
+ /* // (sensor printout is verbose, uncomment to enable)
+ if(imod==0) {
+ printout(ALWAYS, "IRTLOG", "");
+ printout(ALWAYS, "IRTLOG", " SECTOR %d SENSORS:", isec);
+ }
+ printout(ALWAYS, "IRTLOG", " sensor (imodsec,x,y,z) = 0x%08x %5.2f %5.2f %5.2f cm",
+ imodsec, sensorGlobalPos[0], sensorGlobalPos[1], sensorGlobalPos[2]);
+ */
+ }
+#endif
+
// increment sensor module number
imod++;
@@ -512,8 +659,59 @@ static Ref_t createDetector(Detector& desc, xml::Handle_t handle, SensitiveDetec
SkinSurface mirrorSkin(desc, mirrorDE, Form("mirror_optical_surface%d", isec), mirrorSurf, mirrorVol);
mirrorSkin.isValid();
+#ifdef IRT_AUXFILE
+ // get mirror center coordinates, w.r.t. IP
+ /* - we have sector 0 coordinates `(zM,xM,rM)`, but here we try to access the numbers more generally,
+ * so we get the mirror centers after sectorRotation
+ * - FIXME: boolean solids make this a bit tricky, both here and from `GeoSvc`, is there an easier way?
+ */
+ SphericalSurface* mirrorSphericalSurface;
+ OpticalBoundary* mirrorOpticalBoundary;
+ if (createIrtFile) {
+ auto mirrorFinalPlacement = mirrorSectorPlacement * mirrorPlacement;
+ auto mirrorFinalCenter = vesselPos + mirrorFinalPlacement.Translation().Vect(); // w.r.t. IP
+ mirrorSphericalSurface = new SphericalSurface(
+ (1 / mm) * TVector3(mirrorFinalCenter.x(), mirrorFinalCenter.y(), mirrorFinalCenter.z()), mirrorRadius / mm);
+ mirrorOpticalBoundary = new OpticalBoundary(irtDetector->GetContainerVolume(), // CherenkovRadiator radiator
+ mirrorSphericalSurface, // surface
+ false // bool refractive
+ );
+ irtDetector->AddOpticalBoundary(isec, mirrorOpticalBoundary);
+ printout(ALWAYS, "IRTLOG", "");
+ printout(ALWAYS, "IRTLOG", " SECTOR %d MIRROR:", isec);
+ printout(ALWAYS, "IRTLOG", " mirror x = %f cm", mirrorFinalCenter.x());
+ printout(ALWAYS, "IRTLOG", " mirror y = %f cm", mirrorFinalCenter.y());
+ printout(ALWAYS, "IRTLOG", " mirror z = %f cm", mirrorFinalCenter.z());
+ printout(ALWAYS, "IRTLOG", " mirror R = %f cm", mirrorRadius);
+ }
+
+ // IRT: complete the radiator volume description; this is the rear side of the container gas volume
+ if (createIrtFile)
+ irtDetector->GetRadiator("GasVolume")->m_Borders[isec].second = mirrorSphericalSurface;
+#endif
+
}; // END SECTOR LOOP //////////////////////////
+#ifdef IRT_AUXFILE
+ // write IRT auxiliary file
+ if (createIrtFile) {
+ // set refractive indices
+ // FIXME: are these (weighted) averages? can we automate this?
+ std::map<std::string, double> rIndices;
+ rIndices.insert({"GasVolume", 1.0008});
+ rIndices.insert({"Aerogel", 1.0190});
+ rIndices.insert({"Filter", 1.5017});
+ for (auto const& [rName, rIndex] : rIndices) {
+ auto rad = irtDetector->GetRadiator(rName.c_str());
+ if (rad)
+ rad->SetReferenceRefractiveIndex(rIndex);
+ }
+ // write
+ irtGeometry->Write();
+ irtAuxFile->Close();
+ }
+#endif
+
return det;
}