From ecae0203ea6272436f3294fd3f7649b91e6d7a02 Mon Sep 17 00:00:00 2001
From: Wouter Deconinck <wdconinc@gmail.com>
Date: Fri, 12 Nov 2021 01:34:27 +0000
Subject: [PATCH] 2d material scan figure: eta vs phi, colz on integral of X/X0

---
 benchmarks/others/config.yml                  | 67 +++++++++++++----
 .../{materialScan.cxx => materialScanEta.cxx} | 27 +++++--
 benchmarks/others/materialScanEtaPhi.cxx      | 74 +++++++++++++++++++
 3 files changed, 147 insertions(+), 21 deletions(-)
 rename benchmarks/others/{materialScan.cxx => materialScanEta.cxx} (85%)
 create mode 100644 benchmarks/others/materialScanEtaPhi.cxx

diff --git a/benchmarks/others/config.yml b/benchmarks/others/config.yml
index 718b572f..5bef64db 100644
--- a/benchmarks/others/config.yml
+++ b/benchmarks/others/config.yml
@@ -1,34 +1,75 @@
-bench:materialscan:
+variables:
+  ETAMIN: "-4.5"
+  ETAMAX: "+4.5"
+  ETASTEP: "0.01"
+  PHIMIN: "0.0"
+  PHIMAX: "6.28318530718"
+  PHISTEP: "0.01"
+  TRACKING_RHOMAX: "103."
+  TRACKING_ZNMAX: "191."
+  TRACKING_ZPMAX: "350."
+  ECAL_RHOMAX: "230."
+  ECAL_ZNMAX: "355."
+  ECAL_ZPMAX: "380."
+
+bench:materialScanEta:
+  stage: benchmarks
+  extends: .det_benchmark
+  script:
+    - echo ".x benchmarks/others/materialScanEta.cxx($ETAMIN, $ETAMAX, $ETASTEP, $PHIMIN)" | materialScan ${DETECTOR_PATH}/${JUGGLER_DETECTOR}.xml -interactive
+    - mkdir results/images/
+    - mv materialScanEta.png results/images/materialScanEta.png
+    - mv materialScanEta.pdf results/images/materialScanEta.pdf
+
+bench:materialScanEta:tracking:
+  stage: benchmarks
+  extends: .det_benchmark
+  script:
+    - echo ".x benchmarks/others/materialScanEta.cxx($ETAMIN, $ETAMAX, $ETASTEP, $PHIMIN, $TRACKING_RHOMAX, $TRACKING_ZNMAX, $TRACKING_ZPMAX)" | materialScan ${DETECTOR_PATH}/${JUGGLER_DETECTOR}.xml -interactive
+    - mkdir results/images/
+    - mv materialScanEta.png results/images/materialScanEtaTracking.png
+    - mv materialScanEta.pdf results/images/materialScanEtaTracking.pdf
+
+bench:materialScanEta:ecal:
+  stage: benchmarks
+  extends: .det_benchmark
+  script:
+    - echo ".x benchmarks/others/materialScanEta.cxx($ETAMIN, $ETAMAX, $ETASTEP, $PHIMIN, $ECAL_RHOMAX, $ECAL_ZNMAX, $ECAL_ZPMAX)" | materialScan ${DETECTOR_PATH}/${JUGGLER_DETECTOR}.xml -interactive
+    - mkdir results/images/
+    - mv materialScanEta.png results/images/materialScanEtaEcal.png
+    - mv materialScanEta.pdf results/images/materialScanEtaEcal.pdf
+
+bench:materialScanEtaPhi:
   stage: benchmarks
   extends: .det_benchmark
   script:
-    - echo ".x benchmarks/others/materialScan.cxx(-4.5, +4.5, 0.01, 750., 0.0)" | materialScan ${DETECTOR_PATH}/${JUGGLER_DETECTOR}.xml -interactive
+    - echo ".x benchmarks/others/materialScanEtaPhi.cxx($ETAMIN, $ETAMAX, $ETASTEP, $PHIMIN, $PHIMAX, $PHISTEP)" | materialScan ${DETECTOR_PATH}/${JUGGLER_DETECTOR}.xml -interactive
     - mkdir results/images/
-    - mv materialScan.png results/images/materialScan.png
-    - mv materialScan.pdf results/images/materialScan.pdf
+    - mv materialScanEtaPhi.png results/images/materialScanEtaPhi.png
+    - mv materialScanEtaPhi.pdf results/images/materialScanEtaPhi.pdf
 
-bench:materialscan:tracking:
+bench:materialScanEtaPhi:tracking:
   stage: benchmarks
   extends: .det_benchmark
   script:
-    - echo ".x benchmarks/others/materialScan.cxx(-4.5, +4.5, 0.01, 750., 0.0, 103., 191., 350.)" | materialScan ${DETECTOR_PATH}/${JUGGLER_DETECTOR}.xml -interactive
+    - echo ".x benchmarks/others/materialScanEtaPhi.cxx($ETAMIN, $ETAMAX, $ETASTEP, $PHIMIN, $PHIMAX, $PHISTEP, $TRACKING_RHOMAX, $TRACKING_ZNMAX, $TRACKING_ZPMAX)" | materialScan ${DETECTOR_PATH}/${JUGGLER_DETECTOR}.xml -interactive
     - mkdir results/images/
-    - mv materialScan.png results/images/materialScanTracking.png
-    - mv materialScan.pdf results/images/materialScanTracking.pdf
+    - mv materialScanEtaPhi.png results/images/materialScanEtaPhiTracking.png
+    - mv materialScanEtaPhi.pdf results/images/materialScanEtaPhiTracking.pdf
 
-bench:materialscan:ecal:
+bench:materialScanEtaPhi:ecal:
   stage: benchmarks
   extends: .det_benchmark
   script:
-    - echo ".x benchmarks/others/materialScan.cxx(-4.5, +4.5, 0.01, 750., 0.0, 230., 355., 380.)" | materialScan ${DETECTOR_PATH}/${JUGGLER_DETECTOR}.xml -interactive
+    - echo ".x benchmarks/others/materialScanEtaPhi.cxx($ETAMIN, $ETAMAX, $ETASTEP, $PHIMIN, $PHIMAX, $PHISTEP, $ECAL_RHOMAX, $ECAL_ZNMAX, $ECAL_ZPMAX)" | materialScan ${DETECTOR_PATH}/${JUGGLER_DETECTOR}.xml -interactive
     - mkdir results/images/
-    - mv materialScan.png results/images/materialScanEcal.png
-    - mv materialScan.pdf results/images/materialScanEcal.pdf
+    - mv materialScanEtaPhi.png results/images/materialScanEtaPhiEcal.png
+    - mv materialScanEtaPhi.pdf results/images/materialScanEtaPhiEcal.pdf
 
 collect_results:materialscan:
   extends: .det_benchmark
   stage: collect
   needs:
-    - ["bench:materialscan", "bench:materialscan:tracking"]
+    - ["bench:materialScanEta", "bench:materialScanEta:tracking", "bench:materialScanEta:ecal", "bench:materialScanEtaPhi", "bench:materialScanEtaPhi:tracking", "bench:materialScanEtaPhi:ecal"]
   script:
     - ls -lrht
diff --git a/benchmarks/others/materialScan.cxx b/benchmarks/others/materialScanEta.cxx
similarity index 85%
rename from benchmarks/others/materialScan.cxx
rename to benchmarks/others/materialScanEta.cxx
index dbc7cfc8..525bd51c 100644
--- a/benchmarks/others/materialScan.cxx
+++ b/benchmarks/others/materialScanEta.cxx
@@ -40,29 +40,40 @@ Color_t color(const Material& m) {
   }
 }
 
-void materialScan(
+void materialScanEta(
     double etamin = -2.0, // minimum eta
     double etamax = +1.0, // maximum eta
     double etastep = 0.2, // steps in eta
-    double rmax = 100.0, // maximum radius to scan to
     double phi = 0.0, // phi angle for material scan
     double rhomax = 10000.0, // maximum distance from z axis
     double znmax = 10000.0, // maximum negative endcap z plane (positive number)
     double zpmax = 10000.0 // maximum positive endcap z plane (positive number)
 ) {
   // check inputs
-  if (etamin > etamax || rmax <= 0.0) {
+  if (etamin > etamax) {
     std::cout << "Error: ordered eta range required" << std::endl;
     return -1;
   }
+  if (rhomax <= 0.0) {
+    std::cout << "Error: positive rhomax required" << std::endl;
+    return -1;
+  }
+  if (znmax <= 0.0) {
+    std::cout << "Error: positive znmax required" << std::endl;
+    return -1;
+  }
+  if (zpmax <= 0.0) {
+    std::cout << "Error: positive zpmax required" << std::endl;
+    return -1;
+  }
 
   // get material scans
   size_t total{0};
   std::vector<dd4hep::rec::MaterialVec> scan;
   double x0{0}, y0{0}, z0{0};
   for (double eta = etamin; eta <= etamax + 0.5*etastep; eta += etastep) {
-    double theta = 2.0 * (atan(1) - atan(exp(-eta)));
-    double r = min((theta > 0? zpmax: -znmax) / sin(theta), min(rmax, rhomax / cos(theta)));
+    double theta = 2.0 * (atan(1) - atan(exp(-eta))); // |theta| < 90 deg, cos(theta) > 0, sin(theta) can be 0
+    double r = min((theta < 0? -znmax: zpmax) / sin(theta), rhomax / cos(theta)); // theta == 0 results in min(+inf, rhomax)
     double x = r * cos(theta) * cos(phi);
     double y = r * cos(theta) * sin(phi);
     double z = r * sin(theta); 
@@ -105,7 +116,7 @@ void materialScan(
   std::cout << histograms.size() << " histograms created" << std::endl;
 
   // plot histograms as stack
-  THStack hs("hs",Form("Material Scan (r < %.0f cm, rho < %.0f cm, -%.0f cm < z < %.0f cm)", rmax, rhomax, znmax, zpmax));
+  THStack hs("hs",Form("Material Scan (rho < %.0f cm, -%.0f cm < z < %.0f cm)", rhomax, znmax, zpmax));
   for (auto& h: histograms) {
     hs.Add(&h);
   }
@@ -116,6 +127,6 @@ void materialScan(
   hs.GetXaxis()->SetTitle("eta");
   hs.GetYaxis()->SetTitle("Fraction X0");
   hs.SetMinimum(2.5e-3);
-  cs.SaveAs("materialScan.png");
-  cs.SaveAs("materialScan.pdf");
+  cs.SaveAs("materialScanEta.png");
+  cs.SaveAs("materialScanEta.pdf");
 }
diff --git a/benchmarks/others/materialScanEtaPhi.cxx b/benchmarks/others/materialScanEtaPhi.cxx
new file mode 100644
index 00000000..dc670558
--- /dev/null
+++ b/benchmarks/others/materialScanEtaPhi.cxx
@@ -0,0 +1,74 @@
+#import <iostream>
+#import <vector>
+
+#include <DDRec/Material.h>
+#include <DDRec/MaterialScan.h>
+
+#include <TH1F.h>
+#include <THStack.h>
+
+void materialScanEtaPhi(
+    double etamin = -2.0, // minimum eta
+    double etamax = +1.0, // maximum eta
+    double etastep = 0.2, // steps in eta
+    double phimin = 0.0, // minimum phi
+    double phimax = 2.0 * M_PI, // maximum phi
+    double phistep = M_PI / 2, // steps in phi
+    double rhomax = 10000.0, // maximum distance from z axis
+    double znmax = 10000.0, // maximum negative endcap z plane (positive number)
+    double zpmax = 10000.0 // maximum positive endcap z plane (positive number)
+) {
+  // check inputs
+  if (etamin > etamax) {
+    std::cout << "Error: ordered eta range required" << std::endl;
+    return -1;
+  }
+  if (phimin > phimax) {
+    std::cout << "Error: ordered phi range required" << std::endl;
+    return -1;
+  }
+  if (rhomax <= 0.0) {
+    std::cout << "Error: positive rhomax required" << std::endl;
+    return -1;
+  }
+  if (znmax <= 0.0) {
+    std::cout << "Error: positive znmax required" << std::endl;
+    return -1;
+  }
+  if (zpmax <= 0.0) {
+    std::cout << "Error: positive zpmax required" << std::endl;
+    return -1;
+  }
+
+  // get material scans
+  size_t total{0};
+  double x0{0}, y0{0}, z0{0};
+  TH2F h2("h2","Material Scan Eta vs Phi",
+    (etamax-etamin)/etastep+1,etamin,etamax,
+    (phimax-phimin)/phistep+1,phimin,phimax
+  );
+  for (double eta = etamin; eta <= etamax + 0.5*etastep; eta += etastep) {
+    for (double phi = phimin; phi <= phimax + 0.5*phistep; phi += phistep) {
+      double theta = 2.0 * (atan(1) - atan(exp(-eta))); // |theta| < 90 deg, cos(theta) > 0, sin(theta) can be 0 
+      double r = min((theta < 0? -znmax: zpmax) / sin(theta), rhomax / cos(theta)); // theta == 0 results in min(+inf, rhomax)
+      double x = r * cos(theta) * cos(phi);
+      double y = r * cos(theta) * sin(phi);
+      double z = r * sin(theta); 
+      auto scan = gMaterialScan->scan(x0,y0,z0,x,y,z);
+      double X0{0};
+      for (auto& mat_len: scan)
+        X0 +=  mat_len.second / mat_len.first.radLength();
+      h2.Fill(eta,phi,X0);
+    }
+  }
+
+  // plot histograms as stack
+  TCanvas cs("cs","Material Scan Eta vs Phi",1920,1080);
+  auto pad = cs.cd();
+  cs.SetLogz();
+  h2.Draw("colz");
+  h2.GetXaxis()->SetTitle("eta");
+  h2.GetYaxis()->SetTitle("phi");
+  cs.SaveAs("materialScanEtaPhi.png");
+  cs.SaveAs("materialScanEtaPhi.pdf");
+}
-- 
GitLab