Skip to content
Snippets Groups Projects
Commit ecae0203 authored by Wouter Deconinck's avatar Wouter Deconinck
Browse files

2d material scan figure: eta vs phi, colz on integral of X/X0

parent bf1ae78b
No related branches found
No related tags found
1 merge request!1012d material scan figure: eta vs phi, colz on integral of X/X0
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 stage: benchmarks
extends: .det_benchmark extends: .det_benchmark
script: 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/ - mkdir results/images/
- mv materialScan.png results/images/materialScan.png - mv materialScanEtaPhi.png results/images/materialScanEtaPhi.png
- mv materialScan.pdf results/images/materialScan.pdf - mv materialScanEtaPhi.pdf results/images/materialScanEtaPhi.pdf
bench:materialscan:tracking: bench:materialScanEtaPhi:tracking:
stage: benchmarks stage: benchmarks
extends: .det_benchmark extends: .det_benchmark
script: 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/ - mkdir results/images/
- mv materialScan.png results/images/materialScanTracking.png - mv materialScanEtaPhi.png results/images/materialScanEtaPhiTracking.png
- mv materialScan.pdf results/images/materialScanTracking.pdf - mv materialScanEtaPhi.pdf results/images/materialScanEtaPhiTracking.pdf
bench:materialscan:ecal: bench:materialScanEtaPhi:ecal:
stage: benchmarks stage: benchmarks
extends: .det_benchmark extends: .det_benchmark
script: 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/ - mkdir results/images/
- mv materialScan.png results/images/materialScanEcal.png - mv materialScanEtaPhi.png results/images/materialScanEtaPhiEcal.png
- mv materialScan.pdf results/images/materialScanEcal.pdf - mv materialScanEtaPhi.pdf results/images/materialScanEtaPhiEcal.pdf
collect_results:materialscan: collect_results:materialscan:
extends: .det_benchmark extends: .det_benchmark
stage: collect stage: collect
needs: needs:
- ["bench:materialscan", "bench:materialscan:tracking"] - ["bench:materialScanEta", "bench:materialScanEta:tracking", "bench:materialScanEta:ecal", "bench:materialScanEtaPhi", "bench:materialScanEtaPhi:tracking", "bench:materialScanEtaPhi:ecal"]
script: script:
- ls -lrht - ls -lrht
...@@ -40,29 +40,40 @@ Color_t color(const Material& m) { ...@@ -40,29 +40,40 @@ Color_t color(const Material& m) {
} }
} }
void materialScan( void materialScanEta(
double etamin = -2.0, // minimum eta double etamin = -2.0, // minimum eta
double etamax = +1.0, // maximum eta double etamax = +1.0, // maximum eta
double etastep = 0.2, // steps in 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 phi = 0.0, // phi angle for material scan
double rhomax = 10000.0, // maximum distance from z axis double rhomax = 10000.0, // maximum distance from z axis
double znmax = 10000.0, // maximum negative endcap z plane (positive number) double znmax = 10000.0, // maximum negative endcap z plane (positive number)
double zpmax = 10000.0 // maximum positive endcap z plane (positive number) double zpmax = 10000.0 // maximum positive endcap z plane (positive number)
) { ) {
// check inputs // check inputs
if (etamin > etamax || rmax <= 0.0) { if (etamin > etamax) {
std::cout << "Error: ordered eta range required" << std::endl; std::cout << "Error: ordered eta range required" << std::endl;
return -1; 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 // get material scans
size_t total{0}; size_t total{0};
std::vector<dd4hep::rec::MaterialVec> scan; std::vector<dd4hep::rec::MaterialVec> scan;
double x0{0}, y0{0}, z0{0}; double x0{0}, y0{0}, z0{0};
for (double eta = etamin; eta <= etamax + 0.5*etastep; eta += etastep) { for (double eta = etamin; eta <= etamax + 0.5*etastep; eta += etastep) {
double theta = 2.0 * (atan(1) - atan(exp(-eta))); 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? zpmax: -znmax) / sin(theta), min(rmax, rhomax / cos(theta))); 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 x = r * cos(theta) * cos(phi);
double y = r * cos(theta) * sin(phi); double y = r * cos(theta) * sin(phi);
double z = r * sin(theta); double z = r * sin(theta);
...@@ -105,7 +116,7 @@ void materialScan( ...@@ -105,7 +116,7 @@ void materialScan(
std::cout << histograms.size() << " histograms created" << std::endl; std::cout << histograms.size() << " histograms created" << std::endl;
// plot histograms as stack // 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) { for (auto& h: histograms) {
hs.Add(&h); hs.Add(&h);
} }
...@@ -116,6 +127,6 @@ void materialScan( ...@@ -116,6 +127,6 @@ void materialScan(
hs.GetXaxis()->SetTitle("eta"); hs.GetXaxis()->SetTitle("eta");
hs.GetYaxis()->SetTitle("Fraction X0"); hs.GetYaxis()->SetTitle("Fraction X0");
hs.SetMinimum(2.5e-3); hs.SetMinimum(2.5e-3);
cs.SaveAs("materialScan.png"); cs.SaveAs("materialScanEta.png");
cs.SaveAs("materialScan.pdf"); cs.SaveAs("materialScanEta.pdf");
} }
#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");
}
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment