Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 23:31:27

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);  // grey
0168   plot->SetLineColor(1);   // black
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);  // grey
0184   plot->SetLineColor(1);   // black
0185   plot->SetLineWidth(2);
0186   plot->Draw("c e6");
0187   line->SetLineColor(2);  // red
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 }