IRT auxfile creation

modified:   CMakeLists.txt

include IRT headers

add IRT aux file output; use Detector constant accessors; formatting

new file:   scripts/createIRTconfig.py

add CI job for IRT aux file creation

remove WITH_IRT

revert comments

add build option and macro `IRT_AUXFILE`

clang-format

modified:   src/DRICH_geo.cpp

modified:   src/DRICH_geo.cpp

modified:   src/DRICH_geo.cpp

switch to raw pointers, since ROOT is smart enough

create IRT objects iff `createIrtFile==1`

set refractive indices for IRT

clang-format

convince clang-format to not mess up complicated `Transform3D`s

modify comments

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

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"
   >
 
 

scripts/create_IRT_auxfile.py

+# 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()

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;
 }