Create Beamline benchmark (#154)

Adds a first round of electron beamline benchmarks which demonstrate the path taken by the electron beam through the outgoing beamline magnets.

Plots show the position and shape of the beamspot, correlation between position and momentum. The position, rotation and radius of the detector planes is also output.

Co-authored-by: Dmitry Kalinkin <dmitry.kalinkin@gmail.com>

.gitlab-ci.yml

@@ -126,6 +126,7 @@ get_data:
 include: 
   - local: 'benchmarks/backgrounds/config.yml'
   - local: 'benchmarks/backwards_ecal/config.yml'
+  - local: 'benchmarks/beamline/config.yml'
   - local: 'benchmarks/calo_pid/config.yml'
   - local: 'benchmarks/ecal_gaps/config.yml'
   - local: 'benchmarks/tracking_detectors/config.yml'
@@ -159,6 +160,7 @@ deploy_results:
     - "collect_results:backwards_ecal"
     - "collect_results:barrel_ecal"
     - "collect_results:barrel_hcal"
+    - "collect_results:beamline"
     - "collect_results:calo_pid"
     - "collect_results:ecal_gaps"
     - "collect_results:lfhcal"

Snakefile

@@ -33,6 +33,7 @@ def find_epic_libraries():
 include: "benchmarks/backgrounds/Snakefile"
 include: "benchmarks/backwards_ecal/Snakefile"
 include: "benchmarks/barrel_ecal/Snakefile"
+include: "benchmarks/beamline/Snakefile"
 include: "benchmarks/calo_pid/Snakefile"
 include: "benchmarks/ecal_gaps/Snakefile"
 include: "benchmarks/material_scan/Snakefile"

benchmarks/beamline/Snakefile

+SIMOUTDIR="sim_output/beamline/"
+ANALYSISDIR=SIMOUTDIR+"analysis/"
+
+rule beamline_steering_sim:
+    input:
+        macro=workflow.source_path("beamlineGPS.mac"),
+    output:
+        SIMOUTDIR+"beamlineTest{CAMPAIGN}.edm4hep.root",
+    shell:
+        """
+            exec npsim \
+            --runType run \
+            --enableG4GPS \
+            --macroFile {input.macro} \
+            --compactFile $DETECTOR_PATH/epic_ip6_extended.xml \
+            --outputFile {output} \
+        """
+
+rule beamline_steering_analysis:
+    params:
+        xml=os.getenv("DETECTOR_PATH")+"/epic_ip6_extended.xml",
+    input:
+        script=workflow.source_path("beamlineAnalysis.C"),
+        header=workflow.source_path("shared_functions.h"),
+        data=SIMOUTDIR+"beamlineTest{CAMPAIGN}.edm4hep.root",
+    output:
+        rootfile=ANALYSISDIR+"beamlineTestAnalysis{CAMPAIGN}.root",
+        beamspot_canvas=ANALYSISDIR+"beamspot_{CAMPAIGN}.png",
+        x_px_canvas=ANALYSISDIR+"x_px_{CAMPAIGN}.png",
+        y_py_canvas=ANALYSISDIR+"y_py_{CAMPAIGN}.png",
+        fitted_position_means_stdevs_canvas=ANALYSISDIR+"fitted_position_means_stdevs_{CAMPAIGN}.png",
+        fitted_momentum_means_stdevs_canvas=ANALYSISDIR+"fitted_momentum_means_stdevs_{CAMPAIGN}.png",
+        pipe_parameter_canvas=ANALYSISDIR+"pipe_parameter_{CAMPAIGN}.png",
+    shell:
+        """
+            root -l -b -q '{input.script}("{input.data}", "{output.rootfile}", "{params.xml}",
+            "{output.beamspot_canvas}", "{output.x_px_canvas}", "{output.y_py_canvas}",
+            "{output.fitted_position_means_stdevs_canvas}", "{output.fitted_momentum_means_stdevs_canvas}",
+            "{output.pipe_parameter_canvas}")'
+        """
+
+rule beamline:
+    input:
+        ANALYSISDIR+"beamlineTestAnalysis{CAMPAIGN}.root",
+        ANALYSISDIR+"beamspot_{CAMPAIGN}.png",
+        ANALYSISDIR+"x_px_{CAMPAIGN}.png",
+        ANALYSISDIR+"y_py_{CAMPAIGN}.png",
+        ANALYSISDIR+"fitted_position_means_stdevs_{CAMPAIGN}.png",
+        ANALYSISDIR+"fitted_momentum_means_stdevs_{CAMPAIGN}.png",
+        ANALYSISDIR+"pipe_parameter_{CAMPAIGN}.png",        
+    output:
+        directory("results/beamline/{CAMPAIGN}/")
+    shell:
+        """
+            mkdir {output}        
+            cp {input} {output}
+        """
+
+rule beamline_local:
+    input:
+        "results/beamline/local/"

benchmarks/beamline/beamlineGPS.mac

+/gps/ang/type beam2d
+/gps/ang/sigma_x 0.0002017 rad # 201.7 urad
+/gps/ang/sigma_y 0.0001873 rad # 187.3 urad
+/gps/pos/type Beam
+/gps/pos/sigma_x 0.119 mm
+/gps/pos/sigma_y 0.0107 mm
+#/gps/energy 17.846 GeV
+/gps/energy 18 GeV
+/gps/particle e-
+
+/run/beamOn 10000
\ No newline at end of file

benchmarks/beamline/config.yml

+sim:beamline:
+  extends: .det_benchmark
+  stage: simulate
+  script:
+    - |
+      snakemake $SNAKEMAKE_FLAGS --cores 2 \
+        sim_output/beamline/beamlineTestlocal.edm4hep.root
+
+bench:beamline:
+  extends: .det_benchmark
+  stage: benchmarks
+  needs:
+    - ["sim:beamline"]
+  script:
+    - snakemake $SNAKEMAKE_FLAGS --cores 2 beamline_local
+
+collect_results:beamline:
+  extends: .det_benchmark
+  stage: collect
+  needs:
+    - "bench:beamline"
+  script:
+    - ls -lrht
+    - mv results{,_save}/ # move results directory out of the way to preserve it
+    - snakemake $SNAKEMAKE_FLAGS --cores 1 --delete-all-output beamline_local
+    - mv results{_save,}/

benchmarks/beamline/shared_functions.h

+#include "DD4hep/Detector.h"
+#include "DDRec/CellIDPositionConverter.h"
+#include "edm4hep/MCParticleCollection.h"
+#include "edm4hep/SimTrackerHitCollection.h"
+#include "edm4hep/SimCalorimeterHitCollection.h"
+#include "DD4hep/VolumeManager.h"
+#include "TFile.h"
+
+using RVecHits = ROOT::VecOps::RVec<edm4hep::SimTrackerHitData>;
+
+//-----------------------------------------------------------------------------------------
+// Grab Component functor
+//-----------------------------------------------------------------------------------------
+struct getSubID{
+  getSubID(std::string cname, dd4hep::Detector& det, std::string rname = "BackwardsBeamlineHits") : componentName(cname), detector(det), readoutName(rname){}
+  
+  ROOT::VecOps::RVec<int> operator()(const RVecHits& evt) {
+    auto decoder = detector.readout(readoutName).idSpec().decoder();
+    auto indexID = decoder->index(componentName);
+    ROOT::VecOps::RVec<int> result;
+    for(const auto& h: evt) {
+      result.push_back(decoder->get(h.cellID,indexID));      
+    }
+    return result;    
+  };
+  
+  void SetComponent(std::string cname){
+    componentName = cname;
+  }
+  void SetReadout(std::string rname){
+    readoutName = rname;
+  }
+  
+  private: 
+  std::string componentName;
+  dd4hep::Detector& detector;
+  std::string readoutName;
+};
+
+//-----------------------------------------------------------------------------------------
+// Transform global x,y,z position and momentum into local coordinates
+//-----------------------------------------------------------------------------------------
+struct globalToLocal{
+  globalToLocal(dd4hep::Detector& det) : detector(det){
+     volumeManager = dd4hep::VolumeManager::getVolumeManager(det);
+  }
+  
+  ROOT::VecOps::RVec<ROOT::VecOps::RVec<double>> operator()(const RVecHits& evt) {
+
+    ROOT::VecOps::RVec<ROOT::VecOps::RVec<double>> result;
+    ROOT::VecOps::RVec<double> xvec;
+    ROOT::VecOps::RVec<double> yvec;
+    ROOT::VecOps::RVec<double> zvec;
+    ROOT::VecOps::RVec<double> pxvec;
+    ROOT::VecOps::RVec<double> pyvec;
+    ROOT::VecOps::RVec<double> pzvec;
+
+    for(const auto& h: evt) {
+      auto cellID = h.cellID;
+      // dd4hep::DetElement detelement = volumeManager.lookupDetElement(cellID);
+      // auto detelement = volumeManager.lookupVolumePlacement(cellID);
+      auto detelement = volumeManager.lookupDetElement(cellID);
+      const TGeoMatrix& transform = detelement.nominal().worldTransformation();
+      // transform.Print();
+      //auto transform = volumeManager.worldTransformation(cellID);
+      //Convert position to local coordinates
+      auto pos = h.position;
+      Double_t globalPos[3] = {pos.x/10, pos.y/10, pos.z/10};
+      Double_t localPos[3];
+      transform.MasterToLocal(globalPos, localPos);
+      // std::cout << "Global: " << globalPos[0] << " " << globalPos[1] << " " << globalPos[2] << std::endl;
+      // std::cout << "Local: " << localPos[0] << " " << localPos[1] << " " << localPos[2] << std::endl;
+      //Transform global momentum to local coordinates
+      auto mom = h.momentum;
+      Double_t globalMom[3] = {mom.x, mom.y, mom.z};
+      Double_t localMom[3];
+      transform.MasterToLocalVect(globalMom, localMom);
+      // std::cout << "Global: " << globalMom[0] << " " << globalMom[1] << " " << globalMom[2] << std::endl;
+      // std::cout << "Local: " << localMom[0] << " " << localMom[1] << " " << localMom[2] << std::endl;
+
+      xvec.push_back(localPos[0]);
+      yvec.push_back(localPos[1]);
+      zvec.push_back(localPos[2]);
+      pxvec.push_back(localMom[0]);
+      pyvec.push_back(localMom[1]);
+      pzvec.push_back(localMom[2]);
+      
+    }
+
+    result.push_back(xvec);
+    result.push_back(yvec);
+    result.push_back(zvec);
+    result.push_back(pxvec);
+    result.push_back(pyvec);
+    result.push_back(pzvec);    
+
+    return result;    
+  };
+  
+  private: 
+  dd4hep::Detector& detector;
+  dd4hep::VolumeManager volumeManager;
+};
+
+struct volParams{
+  double radius;
+  double xPos;
+  double yPos;
+  double zPos;
+  double rotation;
+};
+
+// Functor to get the volume element parameters from the CellID
+struct getVolumeParametersFromCellID {
+  getVolumeParametersFromCellID(dd4hep::Detector& det) : detector(det) {
+    volumeManager = dd4hep::VolumeManager::getVolumeManager(det);
+  }
+
+  ROOT::VecOps::RVec<volParams> operator()(const RVecHits& evt) {
+    ROOT::VecOps::RVec<volParams> result;
+    // Look up the detector element using the CellID
+    for(const auto& h : evt) {
+      auto  cellID = h.cellID;
+      auto  detelement = volumeManager.lookupDetElement(cellID);
+      const TGeoMatrix& transform      = detelement.nominal().worldTransformation();
+      const Double_t*   translation    = transform.GetTranslation();
+      const Double_t*   rotationMatrix = transform.GetRotationMatrix(); // Compute rotation angle around the Y-axis
+      double rotationAngleY = std::atan2(rotationMatrix[2], rotationMatrix[8]); // atan2(r13, r33)
+      auto volume = detelement.volume();
+      dd4hep::ConeSegment cone = volume.solid();
+      volParams params{
+        cone.rMax1(),
+        translation[0],
+        translation[1],
+        translation[2],
+        rotationAngleY
+      };
+      result.push_back(params);
+    }
+    return result;
+  }
+
+private:
+  dd4hep::Detector& detector;
+  dd4hep::VolumeManager volumeManager;
+};
+
+TH1F* CreateFittedHistogram(const std::string& histName, 
+  const std::string& histTitle, 
+  const std::map<TString, double>& valueMap, 
+  const std::map<TString, double>& errorMap, 
+  const std::string& xAxisLabel) {
+  // Number of bins corresponds to the number of entries in the value map
+  int nPoints = valueMap.size();
+  TH1F* hist = new TH1F(histName.c_str(), (";" + xAxisLabel + ";" + histTitle).c_str(), nPoints, 0.5, nPoints + 0.5);
+
+  // Fill the histogram with values and errors, and set custom bin labels
+  int binIndex = 1; // Start from bin 1
+  for (const auto& [name, value] : valueMap) {
+    hist->SetBinContent(binIndex, value); // Set the bin content
+    if(errorMap.size()==valueMap.size()){
+      hist->SetBinError(binIndex, errorMap.at(name)); // Set the bin error
+    }
+    else{
+      hist->SetBinError(binIndex, 0); // Set the bin error to 0 if not provided
+    }
+    hist->GetXaxis()->SetBinLabel(binIndex, name); // Set the bin label
+    binIndex++;
+  }
+
+  // Customize the histogram
+  hist->SetMarkerStyle(20);
+  hist->SetMarkerSize(1.0);
+  hist->SetLineColor(kBlue);
+  hist->SetMarkerColor(kRed);
+
+  return hist;
+}
\ No newline at end of file