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