Merge pull request #29577 from vargasa/MaterialAccounting

DD4hep: SimTracker/TrackerMaterialAnalysis Migration MaterialAccounting

SimDataFormats/ValidationFormats/interface/MaterialAccountingDetector.h

@@ -15,6 +15,7 @@
 class MaterialAccountingDetector {
   friend class MaterialAccountingTrack;
   friend class TrackingMaterialAnalyser;
+  friend class DD4hep_TrackingMaterialAnalyser;
 
 public:
   MaterialAccountingDetector(void)

SimTracker/TrackerMaterialAnalysis/plugins/dd4hep/BoundingBox.cc

@@ -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); }

SimTracker/TrackerMaterialAnalysis/plugins/dd4hep/BoundingBox.h

@@ -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

SimTracker/TrackerMaterialAnalysis/plugins/dd4hep/DD4hep_Listgroups.cc

@@ -31,6 +31,8 @@
 #include <TText.h>
 #include <TColor.h>
 
+#include "DD4hep_MaterialAccountingGroup.h"
+
 class DD4hep_ListGroups : public edm::one::EDAnalyzer<> {
 public:
   DD4hep_ListGroups(const edm::ParameterSet &iConfig);
@@ -47,8 +49,10 @@ class DD4hep_ListGroups : public edm::one::EDAnalyzer<> {
   std::set<std::string_view> m_group_names;
   std::vector<unsigned int> m_color;
   std::vector<int> m_gradient;
+  std::vector<DD4hep_MaterialAccountingGroup *> m_groups;
   void fillColor();
   void fillGradient();
+  void produceAndSaveSummaryPlot(cms::DDCompactView cpv);
   std::vector<std::pair<std::shared_ptr<TLine>, std::shared_ptr<TText>>> overlayEtaReferences();
 };
 
@@ -59,6 +63,12 @@ DD4hep_ListGroups::DD4hep_ListGroups(const edm::ParameterSet &iConfig)
 
 DD4hep_ListGroups::~DD4hep_ListGroups() {}
 
+void DD4hep_ListGroups::produceAndSaveSummaryPlot(cms::DDCompactView cpv) {
+  for (auto n : m_group_names) {
+    m_groups.push_back(new DD4hep_MaterialAccountingGroup(n.data(), cpv));
+  }
+}
+
 void DD4hep_ListGroups::fillColor(void) {
   // With the introduction of the support for PhaseI and PhaseII detectors it
   // became quite difficult to maintain a list of colors that is in sync with
@@ -239,36 +249,32 @@ void DD4hep_ListGroups::analyze(const edm::Event &evt, const edm::EventSetup &se
 
   for (const auto &t : fv.specpars()) {
     m_group_names.insert(t->strValue("TrackingMaterialGroup"));
-    std::cout << t->strValue("TrackingMaterialGroup") << std::endl;
   }
 
   for (const auto &i : m_group_names) {
     cms::DDFilter filter1("TrackingMaterialGroup", {i.data(), i.size()});
     cms::DDFilteredView fv1(*cpv, filter1);
     bool firstChild = fv1.firstChild();
-    //fv1.printFilter();
-
-    std::cout << "***" << i << std::endl;
 
     for (const auto j : fv1.specpars()) {
       for (const auto k : j->paths) {
-        std::cout << k << std::endl;
         if (firstChild) {
-          std::cout << "Find children matching selection for a parent " << fv1.name() << "\n";
           std::vector<std::vector<cms::Node *>> children = fv1.children(k);
           for (auto const &path : children) {
             for (auto const &node : path) {
-              std::cout << node->GetName() << ", ";
+              edm::LogVerbatim("TrackingMaterialGroup") << node->GetName() << "/";
             }
             cms::Translation trans = fv1.translation(path);
-            std::cout << "(" << trans.x() << ", " << trans.y() << ", " << trans.z() << ")\n";
+            edm::LogVerbatim("TrackingMaterialGroup")
+                << "(" << trans.x() << ", " << trans.y() << ", " << trans.z() << ")\n";
           }
-          std::cout << "\n";
         }
       }
     }
-    std::cout << "*************" << std::endl;
   }
+
+  if (m_saveSummaryPlot)
+    produceAndSaveSummaryPlot(*cpv);
 }
 
 void DD4hep_ListGroups::endJob() {}

SimTracker/TrackerMaterialAnalysis/plugins/dd4hep/DD4hep_MaterialAccountingGroup.cc

@@ -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
+            << 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.;
+}