New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
DD4hep: SimTracker/TrackerMaterialAnalysis Migration MaterialAccounting #29577
Merged
Merged
Changes from all commits
Commits
Show all changes
6 commits
Select commit
Hold shift + click to select a range
d037ad2
Cleanup
vargasa 43cbb3b
Migrating MaterialAccountingGroup,BoundingBox
vargasa c7b7b3b
Use shared_ptrs instead
vargasa 8bc27cd
scram build code-format
vargasa b498780
Mostly intact
vargasa c3dd904
Add migrated class as friend
vargasa File filter
Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
21 changes: 21 additions & 0 deletions
21
SimTracker/TrackerMaterialAnalysis/plugins/dd4hep/BoundingBox.cc
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,21 @@ | ||
#include "BoundingBox.h" | ||
|
||
void BoundingBox::grow(double r, double z) { | ||
if (r < r_min) | ||
r_min = r; | ||
if (r > r_max) | ||
r_max = r; | ||
if (z < z_min) | ||
z_min = z; | ||
if (z > z_max) | ||
z_max = z; | ||
} | ||
|
||
void BoundingBox::grow(double skin) { | ||
r_min -= skin; // yes, we allow r_min to go negative | ||
r_max += skin; | ||
z_min -= skin; | ||
z_max += skin; | ||
} | ||
|
||
bool BoundingBox::inside(double r, double z) const { return (r >= r_min and r <= r_max and z >= z_min and z <= z_max); } |
30 changes: 30 additions & 0 deletions
30
SimTracker/TrackerMaterialAnalysis/plugins/dd4hep/BoundingBox.h
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,30 @@ | ||
#ifndef BoundingBox_h | ||
#define BoundingBox_h | ||
|
||
#include <utility> | ||
|
||
class BoundingBox { | ||
private: | ||
double r_min; | ||
double r_max; | ||
double z_min; | ||
double z_max; | ||
|
||
public: | ||
BoundingBox() : r_min(0.), r_max(0.), z_min(0.), z_max(0.) {} | ||
|
||
BoundingBox(double min_r, double max_r, double min_z, double max_z) | ||
: r_min(min_r), r_max(max_r), z_min(min_z), z_max(max_z) {} | ||
|
||
void grow(double r, double z); | ||
|
||
void grow(double skin); | ||
|
||
bool inside(double r, double z) const; | ||
|
||
std::pair<double, double> range_r() const { return std::make_pair(r_min, r_max); } | ||
|
||
std::pair<double, double> range_z() const { return std::make_pair(z_min, z_max); } | ||
}; | ||
|
||
#endif |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
236 changes: 236 additions & 0 deletions
236
SimTracker/TrackerMaterialAnalysis/plugins/dd4hep/DD4hep_MaterialAccountingGroup.cc
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,236 @@ | ||
#include <sstream> | ||
#include <iomanip> | ||
#include <string> | ||
#include <stdexcept> | ||
|
||
#include <TFile.h> | ||
#include <TH1F.h> | ||
#include <TProfile.h> | ||
#include <TCanvas.h> | ||
#include <TFrame.h> | ||
|
||
#include "DataFormats/GeometryVector/interface/GlobalPoint.h" | ||
#include "DataFormats/Math/interface/Vector3D.h" | ||
#include "DetectorDescription/DDCMS/interface/Filter.h" | ||
#include "DetectorDescription/DDCMS/interface/DDFilteredView.h" | ||
#include "SimDataFormats/ValidationFormats/interface/MaterialAccountingStep.h" | ||
#include "SimDataFormats/ValidationFormats/interface/MaterialAccountingDetector.h" | ||
#include "DD4hep_MaterialAccountingGroup.h" | ||
|
||
DD4hep_MaterialAccountingGroup::DD4hep_MaterialAccountingGroup(const std::string& name, const cms::DDCompactView& geometry) | ||
: m_name(name), | ||
m_elements(), | ||
m_boundingbox(), | ||
m_accounting(), | ||
m_errors(), | ||
m_tracks(0), | ||
m_counted(false), | ||
m_file(nullptr) { | ||
cms::DDFilter filter("TrackingMaterialGroup", {name.data(), name.size()}); | ||
cms::DDFilteredView fv(geometry, filter); | ||
bool firstChild = fv.firstChild(); | ||
|
||
edm::LogVerbatim("TrackingMaterialAnalysis") << "Elements within: " << name; | ||
|
||
for (const auto j : fv.specpars()) { | ||
for (const auto k : j->paths) { | ||
if (firstChild) { | ||
std::vector<std::vector<cms::Node*>> children = fv.children(k); | ||
for (auto const& path : children) { | ||
cms::Translation trans = fv.translation(path); | ||
GlobalPoint gp = GlobalPoint(trans.x(), trans.y(), trans.z()); | ||
m_elements.push_back(gp); | ||
std::cout << "Adding element at (r,z) " << gp.perp() << "," << gp.z() << std::endl; | ||
} | ||
} | ||
} | ||
} | ||
|
||
for (unsigned int i = 0; i < m_elements.size(); ++i) { | ||
m_boundingbox.grow(m_elements[i].perp(), m_elements[i].z()); | ||
} | ||
|
||
m_boundingbox.grow(s_tolerance); | ||
|
||
std::cout << "Final BBox r_range: " << m_boundingbox.range_r().first << ", " << m_boundingbox.range_r().second | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. also here |
||
<< std::endl | ||
<< "Final BBox z_range: " << m_boundingbox.range_z().first << ", " << m_boundingbox.range_z().second | ||
<< std::endl; | ||
|
||
// initialize the histograms | ||
m_dedx_spectrum = | ||
std::make_shared<TH1F>(TH1F((m_name + "_dedx_spectrum").c_str(), "Energy loss spectrum", 1000, 0., 1.)); | ||
m_radlen_spectrum = | ||
std::make_shared<TH1F>(TH1F((m_name + "_radlen_spectrum").c_str(), "Radiation lengths spectrum", 1000, 0., 1.)); | ||
m_dedx_vs_eta = | ||
std::make_shared<TProfile>(TProfile((m_name + "_dedx_vs_eta").c_str(), "Energy loss vs. eta", 600, -3., 3.)); | ||
m_dedx_vs_z = | ||
std::make_shared<TProfile>(TProfile((m_name + "_dedx_vs_z").c_str(), "Energy loss vs. Z", 6000, -300., 300.)); | ||
m_dedx_vs_r = | ||
std::make_shared<TProfile>(TProfile((m_name + "_dedx_vs_r").c_str(), "Energy loss vs. R", 1200, 0., 120.)); | ||
m_radlen_vs_eta = std::make_shared<TProfile>( | ||
TProfile((m_name + "_radlen_vs_eta").c_str(), "Radiation lengths vs. eta", 600, -3., 3.)); | ||
m_radlen_vs_z = std::make_shared<TProfile>( | ||
TProfile((m_name + "_radlen_vs_z").c_str(), "Radiation lengths vs. Z", 6000, -300., 300.)); | ||
m_radlen_vs_r = std::make_shared<TProfile>( | ||
TProfile((m_name + "_radlen_vs_r").c_str(), "Radiation lengths vs. R", 1200, 0., 120.)); | ||
|
||
m_dedx_spectrum->SetDirectory(nullptr); | ||
m_radlen_spectrum->SetDirectory(nullptr); | ||
m_dedx_vs_eta->SetDirectory(nullptr); | ||
m_dedx_vs_z->SetDirectory(nullptr); | ||
m_dedx_vs_r->SetDirectory(nullptr); | ||
m_radlen_vs_eta->SetDirectory(nullptr); | ||
m_radlen_vs_z->SetDirectory(nullptr); | ||
m_radlen_vs_r->SetDirectory(nullptr); | ||
} | ||
|
||
bool DD4hep_MaterialAccountingGroup::isInside(const MaterialAccountingDetector& detector) const { | ||
const GlobalPoint& position = detector.position(); | ||
|
||
edm::LogVerbatim("MaterialAccountingGroup") | ||
<< "Testing position: (x, y, z, r) = " << position.x() << ", " << position.y() << ", " << position.z() << ", " | ||
<< position.perp() << std::endl; | ||
|
||
if (not m_boundingbox.inside(position.perp(), position.z())) { | ||
edm::LogVerbatim("MaterialAccountingGroup") | ||
<< "r outside of: (" << m_boundingbox.range_r().first << ", " << m_boundingbox.range_r().second | ||
<< "), Z ouside of: (" << m_boundingbox.range_z().first << ", " << m_boundingbox.range_z().second << ")" | ||
<< std::endl; | ||
return false; | ||
} else { | ||
edm::LogVerbatim("MaterialAccountingGroup") | ||
<< "r within: (" << m_boundingbox.range_r().first << ", " << m_boundingbox.range_r().second << "), Z within: (" | ||
<< m_boundingbox.range_z().first << ", " << m_boundingbox.range_z().second << ")" << std::endl; | ||
|
||
for (unsigned int i = 0; i < m_elements.size(); ++i) { | ||
edm::LogVerbatim("MaterialAccountingGroup") | ||
<< "Closest testing agains(x, y, z, r): (" << m_elements[i].x() << ", " << m_elements[i].y() << ", " | ||
<< m_elements[i].z() << ", " << m_elements[i].perp() << ") --> " << (position - m_elements[i]).mag() | ||
<< " vs tolerance: " << s_tolerance << std::endl; | ||
if ((position - m_elements[i]).mag2() < (s_tolerance * s_tolerance)) | ||
return true; | ||
} | ||
return false; | ||
} | ||
} | ||
|
||
bool DD4hep_MaterialAccountingGroup::addDetector(const MaterialAccountingDetector& detector) { | ||
if (!isInside(detector)) | ||
return false; | ||
|
||
m_buffer += detector.material(); | ||
m_counted = true; | ||
|
||
return true; | ||
} | ||
|
||
void DD4hep_MaterialAccountingGroup::endOfTrack(void) { | ||
if (m_counted) { | ||
m_accounting += m_buffer; | ||
m_errors += m_buffer * m_buffer; | ||
++m_tracks; | ||
|
||
GlobalPoint average((m_buffer.in().x() + m_buffer.out().x()) / 2., | ||
(m_buffer.in().y() + m_buffer.out().y()) / 2., | ||
(m_buffer.in().z() + m_buffer.out().z()) / 2.); | ||
m_dedx_spectrum->Fill(m_buffer.energyLoss()); | ||
m_radlen_spectrum->Fill(m_buffer.radiationLengths()); | ||
m_dedx_vs_eta->Fill(average.eta(), m_buffer.energyLoss(), 1.); | ||
m_dedx_vs_z->Fill(average.z(), m_buffer.energyLoss(), 1.); | ||
m_dedx_vs_r->Fill(average.perp(), m_buffer.energyLoss(), 1.); | ||
m_radlen_vs_eta->Fill(average.eta(), m_buffer.radiationLengths(), 1.); | ||
m_radlen_vs_z->Fill(average.z(), m_buffer.radiationLengths(), 1.); | ||
m_radlen_vs_r->Fill(average.perp(), m_buffer.radiationLengths(), 1.); | ||
} | ||
m_counted = false; | ||
m_buffer = MaterialAccountingStep(); | ||
} | ||
|
||
void DD4hep_MaterialAccountingGroup::savePlots(void) { | ||
m_file = new TFile((m_name + ".root").c_str(), "RECREATE"); | ||
savePlot(m_dedx_spectrum, m_name + "_dedx_spectrum"); | ||
savePlot(m_radlen_spectrum, m_name + "_radlen_spectrum"); | ||
savePlot(m_dedx_vs_eta, averageEnergyLoss(), m_name + "_dedx_vs_eta"); | ||
savePlot(m_dedx_vs_z, averageEnergyLoss(), m_name + "_dedx_vs_z"); | ||
savePlot(m_dedx_vs_r, averageEnergyLoss(), m_name + "_dedx_vs_r"); | ||
savePlot(m_radlen_vs_eta, averageRadiationLengths(), m_name + "_radlen_vs_eta"); | ||
savePlot(m_radlen_vs_z, averageRadiationLengths(), m_name + "_radlen_vs_z"); | ||
savePlot(m_radlen_vs_r, averageRadiationLengths(), m_name + "_radlen_vs_r"); | ||
m_file->Write(); | ||
m_file->Close(); | ||
delete m_file; | ||
} | ||
|
||
void DD4hep_MaterialAccountingGroup::savePlot(std::shared_ptr<TH1F> plot, const std::string& name) { | ||
TCanvas canvas(name.c_str(), plot->GetTitle(), 1280, 1024); | ||
plot->SetFillColor(15); // grey | ||
plot->SetLineColor(1); // black | ||
plot->Draw("c e"); | ||
canvas.GetFrame()->SetFillColor(kWhite); | ||
canvas.Draw(); | ||
canvas.SaveAs((name + ".png").c_str(), ""); | ||
plot->SetDirectory(m_file); | ||
} | ||
|
||
void DD4hep_MaterialAccountingGroup::savePlot(std::shared_ptr<TProfile> plot, float average, const std::string& name) { | ||
std::unique_ptr<TH1F> line = std::make_unique<TH1F>( | ||
TH1F((name + "_par").c_str(), "Parametrization", 1, plot->GetXaxis()->GetXmin(), plot->GetXaxis()->GetXmax())); | ||
|
||
line->SetBinContent(1, average); | ||
|
||
TCanvas canvas(name.c_str(), plot->GetTitle(), 1280, 1024); | ||
plot->SetFillColor(15); // grey | ||
plot->SetLineColor(1); // black | ||
plot->SetLineWidth(2); | ||
plot->Draw("c e6"); | ||
line->SetLineColor(2); // red | ||
line->SetLineWidth(2); | ||
line->Draw("same"); | ||
canvas.GetFrame()->SetFillColor(kWhite); | ||
canvas.Draw(); | ||
canvas.SaveAs((name + ".png").c_str(), ""); | ||
plot->SetDirectory(m_file); | ||
line->SetDirectory(m_file); | ||
line->Write(); | ||
} | ||
|
||
std::string DD4hep_MaterialAccountingGroup::info(void) const { | ||
std::stringstream out; | ||
out << std::setw(48) << std::left << m_name << std::right << std::fixed; | ||
|
||
out << "BBox: " << std::setprecision(1) << std::setw(6) << m_boundingbox.range_z().first << " < Z < " | ||
<< std::setprecision(1) << std::setw(6) << m_boundingbox.range_z().second; | ||
out << ", " << std::setprecision(1) << std::setw(5) << m_boundingbox.range_r().first << " < R < " | ||
<< std::setprecision(1) << std::setw(5) << m_boundingbox.range_r().second; | ||
out << " Elements: " << std::setw(6) << m_elements.size(); | ||
return out.str(); | ||
} | ||
|
||
MaterialAccountingStep DD4hep_MaterialAccountingGroup::average(void) const { | ||
return m_tracks ? m_accounting / m_tracks : MaterialAccountingStep(); | ||
} | ||
|
||
double DD4hep_MaterialAccountingGroup::averageLength(void) const { return m_tracks ? m_accounting.length() / m_tracks : 0.; } | ||
|
||
double DD4hep_MaterialAccountingGroup::averageEnergyLoss(void) const { | ||
return m_tracks ? m_accounting.energyLoss() / m_tracks : 0.; | ||
} | ||
|
||
double DD4hep_MaterialAccountingGroup::sigmaLength(void) const { | ||
return m_tracks ? std::sqrt(m_errors.length() / m_tracks - averageLength() * averageLength()) : 0.; | ||
} | ||
|
||
double DD4hep_MaterialAccountingGroup::sigmaRadiationLengths(void) const { | ||
return m_tracks | ||
? std::sqrt(m_errors.radiationLengths() / m_tracks - averageRadiationLengths() * averageRadiationLengths()) | ||
: 0.; | ||
} | ||
|
||
double DD4hep_MaterialAccountingGroup::sigmaEnergyLoss(void) const { | ||
return m_tracks ? std::sqrt(m_errors.energyLoss() / m_tracks - averageEnergyLoss() * averageEnergyLoss()) : 0.; | ||
} | ||
|
||
double DD4hep_MaterialAccountingGroup::averageRadiationLengths(void) const { | ||
return m_tracks ? m_accounting.radiationLengths() / m_tracks : 0.; | ||
} |
Oops, something went wrong.
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Shouldn't all these
cout
statements be changed toLogVerbatim
? Is this program only for testing?There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
@cvuosalo This specific part has not changed from the existing code. Due to time constraints I'm focusing on having it working and we can address these details in future iterations