File indexing completed on 2024-04-06 12:31:03
0001 #include <sstream>
0002 #include <iomanip>
0003 #include <string>
0004 #include <stdexcept>
0005
0006 #include <TFile.h>
0007 #include <TH1F.h>
0008 #include <TProfile.h>
0009 #include <TCanvas.h>
0010 #include <TFrame.h>
0011
0012 #include <DD4hep/DD4hepUnits.h>
0013
0014 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0015
0016 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0017 #include "DataFormats/Math/interface/Vector3D.h"
0018 #include "DetectorDescription/DDCMS/interface/DDFilteredView.h"
0019 #include "SimDataFormats/ValidationFormats/interface/MaterialAccountingStep.h"
0020 #include "SimDataFormats/ValidationFormats/interface/MaterialAccountingDetector.h"
0021 #include "DD4hep_MaterialAccountingGroup.h"
0022
0023 DD4hep_MaterialAccountingGroup::DD4hep_MaterialAccountingGroup(const std::string& name,
0024 const cms::DDCompactView& geometry)
0025 : m_name(name),
0026 m_elements(),
0027 m_boundingbox(),
0028 m_accounting(),
0029 m_errors(),
0030 m_tracks(0),
0031 m_counted(false),
0032 m_file(nullptr) {
0033 cms::DDFilter filter("TrackingMaterialGroup", {name.data(), name.size()});
0034 cms::DDFilteredView fv(geometry, filter);
0035 bool firstChild = fv.firstChild();
0036
0037 edm::LogVerbatim("TrackingMaterialAnalysis") << "Elements within: " << name;
0038
0039 for (const auto& j : fv.specpars()) {
0040 for (const auto& k : j.second->paths) {
0041 if (firstChild) {
0042 std::vector<std::vector<cms::Node*>> children = fv.children(k);
0043 for (auto const& path : children) {
0044 cms::Translation trans = fv.translation(path) / dd4hep::cm;
0045 GlobalPoint gp = GlobalPoint(trans.x(), trans.y(), trans.z());
0046 m_elements.emplace_back(gp);
0047 edm::LogVerbatim("TrackerMaterialAnalysis")
0048 << "MaterialAccountingGroup:\t"
0049 << "Adding element at (r,z) " << gp.perp() << "," << gp.z() << std::endl;
0050 }
0051 }
0052 }
0053 }
0054
0055 for (unsigned int i = 0; i < m_elements.size(); ++i) {
0056 m_boundingbox.grow(m_elements[i].perp(), m_elements[i].z());
0057 }
0058
0059 m_boundingbox.grow(s_tolerance);
0060
0061 edm::LogVerbatim("TrackerMaterialAnalysis")
0062 << "MaterialAccountingGroup:\t"
0063 << "Final BBox r_range: " << m_boundingbox.range_r().first << ", " << m_boundingbox.range_r().second << std::endl
0064 << "Final BBox z_range: " << m_boundingbox.range_z().first << ", " << m_boundingbox.range_z().second << std::endl;
0065
0066 m_dedx_spectrum = std::make_shared<TH1F>((m_name + "_dedx_spectrum").c_str(), "Energy loss spectrum", 1000, 0., 1.);
0067 m_radlen_spectrum =
0068 std::make_shared<TH1F>((m_name + "_radlen_spectrum").c_str(), "Radiation lengths spectrum", 1000, 0., 1.);
0069 m_dedx_vs_eta = std::make_shared<TProfile>((m_name + "_dedx_vs_eta").c_str(), "Energy loss vs. eta", 600, -3., 3.);
0070 m_dedx_vs_z = std::make_shared<TProfile>((m_name + "_dedx_vs_z").c_str(), "Energy loss vs. Z", 6000, -300., 300.);
0071 m_dedx_vs_r = std::make_shared<TProfile>((m_name + "_dedx_vs_r").c_str(), "Energy loss vs. R", 1200, 0., 120.);
0072 m_radlen_vs_eta =
0073 std::make_shared<TProfile>((m_name + "_radlen_vs_eta").c_str(), "Radiation lengths vs. eta", 600, -3., 3.);
0074 m_radlen_vs_z =
0075 std::make_shared<TProfile>((m_name + "_radlen_vs_z").c_str(), "Radiation lengths vs. Z", 6000, -300., 300.);
0076 m_radlen_vs_r =
0077 std::make_shared<TProfile>((m_name + "_radlen_vs_r").c_str(), "Radiation lengths vs. R", 1200, 0., 120.);
0078
0079 m_dedx_spectrum->SetDirectory(nullptr);
0080 m_radlen_spectrum->SetDirectory(nullptr);
0081 m_dedx_vs_eta->SetDirectory(nullptr);
0082 m_dedx_vs_z->SetDirectory(nullptr);
0083 m_dedx_vs_r->SetDirectory(nullptr);
0084 m_radlen_vs_eta->SetDirectory(nullptr);
0085 m_radlen_vs_z->SetDirectory(nullptr);
0086 m_radlen_vs_r->SetDirectory(nullptr);
0087 }
0088
0089 bool DD4hep_MaterialAccountingGroup::isInside(const MaterialAccountingDetector& detector) const {
0090 const GlobalPoint& position = detector.position();
0091
0092 edm::LogVerbatim("MaterialAccountingGroup")
0093 << "Testing position: (x, y, z, r) = " << position.x() << ", " << position.y() << ", " << position.z() << ", "
0094 << position.perp() << std::endl;
0095
0096 if (not m_boundingbox.inside(position.perp(), position.z())) {
0097 edm::LogVerbatim("MaterialAccountingGroup")
0098 << "r outside of: (" << m_boundingbox.range_r().first << ", " << m_boundingbox.range_r().second
0099 << "), Z ouside of: (" << m_boundingbox.range_z().first << ", " << m_boundingbox.range_z().second << ")"
0100 << std::endl;
0101 return false;
0102 } else {
0103 edm::LogVerbatim("MaterialAccountingGroup")
0104 << "r within: (" << m_boundingbox.range_r().first << ", " << m_boundingbox.range_r().second << "), Z within: ("
0105 << m_boundingbox.range_z().first << ", " << m_boundingbox.range_z().second << ")" << std::endl;
0106
0107 for (unsigned int i = 0; i < m_elements.size(); ++i) {
0108 edm::LogVerbatim("MaterialAccountingGroup")
0109 << "Closest testing agains(x, y, z, r): (" << m_elements[i].x() << ", " << m_elements[i].y() << ", "
0110 << m_elements[i].z() << ", " << m_elements[i].perp() << ") --> " << (position - m_elements[i]).mag()
0111 << " vs tolerance: " << s_tolerance << std::endl;
0112 if ((position - m_elements[i]).mag2() < (s_tolerance * s_tolerance))
0113 return true;
0114 }
0115 return false;
0116 }
0117 }
0118
0119 bool DD4hep_MaterialAccountingGroup::addDetector(const MaterialAccountingDetector& detector) {
0120 if (!isInside(detector))
0121 return false;
0122
0123 m_buffer += detector.material();
0124 m_counted = true;
0125
0126 return true;
0127 }
0128
0129 void DD4hep_MaterialAccountingGroup::endOfTrack(void) {
0130 if (m_counted) {
0131 m_accounting += m_buffer;
0132 m_errors += m_buffer * m_buffer;
0133 ++m_tracks;
0134
0135 GlobalPoint average((m_buffer.in().x() + m_buffer.out().x()) / 2.,
0136 (m_buffer.in().y() + m_buffer.out().y()) / 2.,
0137 (m_buffer.in().z() + m_buffer.out().z()) / 2.);
0138 m_dedx_spectrum->Fill(m_buffer.energyLoss());
0139 m_radlen_spectrum->Fill(m_buffer.radiationLengths());
0140 m_dedx_vs_eta->Fill(average.eta(), m_buffer.energyLoss(), 1.);
0141 m_dedx_vs_z->Fill(average.z(), m_buffer.energyLoss(), 1.);
0142 m_dedx_vs_r->Fill(average.perp(), m_buffer.energyLoss(), 1.);
0143 m_radlen_vs_eta->Fill(average.eta(), m_buffer.radiationLengths(), 1.);
0144 m_radlen_vs_z->Fill(average.z(), m_buffer.radiationLengths(), 1.);
0145 m_radlen_vs_r->Fill(average.perp(), m_buffer.radiationLengths(), 1.);
0146 }
0147 m_counted = false;
0148 m_buffer = MaterialAccountingStep();
0149 }
0150
0151 void DD4hep_MaterialAccountingGroup::savePlots(void) {
0152 m_file = std::make_unique<TFile>((m_name + ".root").c_str(), "RECREATE");
0153 savePlot(m_dedx_spectrum, m_name + "_dedx_spectrum");
0154 savePlot(m_radlen_spectrum, m_name + "_radlen_spectrum");
0155 savePlot(m_dedx_vs_eta, averageEnergyLoss(), m_name + "_dedx_vs_eta");
0156 savePlot(m_dedx_vs_z, averageEnergyLoss(), m_name + "_dedx_vs_z");
0157 savePlot(m_dedx_vs_r, averageEnergyLoss(), m_name + "_dedx_vs_r");
0158 savePlot(m_radlen_vs_eta, averageRadiationLengths(), m_name + "_radlen_vs_eta");
0159 savePlot(m_radlen_vs_z, averageRadiationLengths(), m_name + "_radlen_vs_z");
0160 savePlot(m_radlen_vs_r, averageRadiationLengths(), m_name + "_radlen_vs_r");
0161 m_file->Write();
0162 m_file->Close();
0163 }
0164
0165 void DD4hep_MaterialAccountingGroup::savePlot(std::shared_ptr<TH1F>& plot, const std::string& name) {
0166 TCanvas canvas(name.c_str(), plot->GetTitle(), 1280, 1024);
0167 plot->SetFillColor(15);
0168 plot->SetLineColor(1);
0169 plot->Draw("c e");
0170 canvas.GetFrame()->SetFillColor(kWhite);
0171 canvas.Draw();
0172 canvas.SaveAs((name + ".png").c_str(), "");
0173 plot->SetDirectory(m_file.get());
0174 }
0175
0176 void DD4hep_MaterialAccountingGroup::savePlot(std::shared_ptr<TProfile>& plot, float average, const std::string& name) {
0177 std::unique_ptr<TH1F> line = std::make_unique<TH1F>(
0178 (name + "_par").c_str(), "Parametrization", 1, plot->GetXaxis()->GetXmin(), plot->GetXaxis()->GetXmax());
0179
0180 line->SetBinContent(1, average);
0181
0182 TCanvas canvas(name.c_str(), plot->GetTitle(), 1280, 1024);
0183 plot->SetFillColor(15);
0184 plot->SetLineColor(1);
0185 plot->SetLineWidth(2);
0186 plot->Draw("c e6");
0187 line->SetLineColor(2);
0188 line->SetLineWidth(2);
0189 line->Draw("same");
0190 canvas.GetFrame()->SetFillColor(kWhite);
0191 canvas.Draw();
0192 canvas.SaveAs((name + ".png").c_str(), "");
0193 plot->SetDirectory(m_file.get());
0194 line->SetDirectory(m_file.get());
0195 line->Write();
0196 }
0197
0198 std::string DD4hep_MaterialAccountingGroup::info(void) const {
0199 std::stringstream out;
0200 out << std::setw(48) << std::left << m_name << std::right << std::fixed;
0201
0202 out << "BBox: " << std::setprecision(1) << std::setw(6) << m_boundingbox.range_z().first << " < Z < "
0203 << std::setprecision(1) << std::setw(6) << m_boundingbox.range_z().second;
0204 out << ", " << std::setprecision(1) << std::setw(5) << m_boundingbox.range_r().first << " < R < "
0205 << std::setprecision(1) << std::setw(5) << m_boundingbox.range_r().second;
0206 out << " Elements: " << std::setw(6) << m_elements.size();
0207 return out.str();
0208 }
0209
0210 MaterialAccountingStep DD4hep_MaterialAccountingGroup::average(void) const {
0211 return m_tracks ? m_accounting / m_tracks : MaterialAccountingStep();
0212 }
0213
0214 double DD4hep_MaterialAccountingGroup::averageLength(void) const {
0215 return m_tracks ? m_accounting.length() / m_tracks : 0.;
0216 }
0217
0218 double DD4hep_MaterialAccountingGroup::averageEnergyLoss(void) const {
0219 return m_tracks ? m_accounting.energyLoss() / m_tracks : 0.;
0220 }
0221
0222 double DD4hep_MaterialAccountingGroup::sigmaLength(void) const {
0223 return m_tracks ? std::sqrt(m_errors.length() / m_tracks - averageLength() * averageLength()) : 0.;
0224 }
0225
0226 double DD4hep_MaterialAccountingGroup::sigmaRadiationLengths(void) const {
0227 return m_tracks
0228 ? std::sqrt(m_errors.radiationLengths() / m_tracks - averageRadiationLengths() * averageRadiationLengths())
0229 : 0.;
0230 }
0231
0232 double DD4hep_MaterialAccountingGroup::sigmaEnergyLoss(void) const {
0233 return m_tracks ? std::sqrt(m_errors.energyLoss() / m_tracks - averageEnergyLoss() * averageEnergyLoss()) : 0.;
0234 }
0235
0236 double DD4hep_MaterialAccountingGroup::averageRadiationLengths(void) const {
0237 return m_tracks ? m_accounting.radiationLengths() / m_tracks : 0.;
0238 }
