MaterialAccountingGroup.cc

Go to the documentation of this file.

#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/Core/interface/DDFilteredView.h"
#include "DetectorDescription/Core/interface/DDCompactView.h"
#include "SimDataFormats/ValidationFormats/interface/MaterialAccountingStep.h"
#include "SimDataFormats/ValidationFormats/interface/MaterialAccountingDetector.h"
#include "MaterialAccountingGroup.h"

double const MaterialAccountingGroup::s_tolerance =
    0.01;  // 100um should be small enough that no elements from different layers/groups are so close

MaterialAccountingGroup::MaterialAccountingGroup(const std::string& name, const DDCompactView& geometry)
    : m_name(name),
      m_elements(),
      m_boundingbox(),
      m_accounting(),
      m_errors(),
      m_tracks(0),
      m_counted(false),
      m_file(nullptr) {
  // retrieve the elements from DDD
  DDValue namevalue;
  if (TString(name.c_str()).Contains("Tracker")) {
    namevalue = DDValue("TrackingMaterialGroup", name);
  } else if (TString(name.c_str()).Contains("HGCal")) {
    namevalue = DDValue("Volume", name);
  } else if (TString(name.c_str()).Contains("HFNose")) {
    namevalue = DDValue("Volume", name);
  } else {
    LogTrace("MaterialAccountingGroup") << name << std::endl;
    edm::LogError("MaterialAccountingGroup") << "Only Tracker , HGCal and HFNose are supported" << std::endl;
  }

  DDSpecificsMatchesValueFilter filter{namevalue};
  DDFilteredView fv(geometry, filter);
  LogTrace("MaterialAccountingGroup") << "Elements within: " << name << std::endl;
  while (fv.next()) {
    // DD3Vector and DDTranslation are the same type as math::XYZVector
    math::XYZVector position = fv.translation() / 10.;  // mm -> cm
    LogTrace("MaterialAccountingGroup") << "Adding element at(r,z): ("
                                        << GlobalPoint(position.x(), position.y(), position.z()).perp() << ", "
                                        << GlobalPoint(position.x(), position.y(), position.z()).z() << ") cm"
                                        << std::endl;
    LogTrace("MaterialAccountingGroup") << "Name of added element: " << fv.logicalPart().toString() << std::endl;
    m_elements.push_back(GlobalPoint(position.x(), position.y(), position.z()));
  }

  // grow the bounding box
  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);
  LogTrace("MaterialAccountingGroup") << "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 = new TH1F((m_name + "_dedx_spectrum").c_str(), "Energy loss spectrum", 1000, 0, 1);
  m_radlen_spectrum = new TH1F((m_name + "_radlen_spectrum").c_str(), "Radiation lengths spectrum", 1000, 0, 1);
  m_dedx_vs_eta = new TProfile((m_name + "_dedx_vs_eta").c_str(), "Energy loss vs. eta", 600, -3, 3);
  m_dedx_vs_z = new TProfile((m_name + "_dedx_vs_z").c_str(), "Energy loss vs. Z", 6000, -300, 300);
  m_dedx_vs_r = new TProfile((m_name + "_dedx_vs_r").c_str(), "Energy loss vs. R", 1200, 0, 120);
  m_radlen_vs_eta = new TProfile((m_name + "_radlen_vs_eta").c_str(), "Radiation lengths vs. eta", 600, -3, 3);
  m_radlen_vs_z = new TProfile((m_name + "_radlen_vs_z").c_str(), "Radiation lengths vs. Z", 6000, -300, 300);
  m_radlen_vs_r = new 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);
}

MaterialAccountingGroup::~MaterialAccountingGroup(void) {
  delete m_dedx_spectrum;
  delete m_dedx_vs_eta;
  delete m_dedx_vs_z;
  delete m_dedx_vs_r;
  delete m_radlen_spectrum;
  delete m_radlen_vs_eta;
  delete m_radlen_vs_z;
  delete m_radlen_vs_r;
}

// TODO the inner check could be sped up in many ways
// (sorting the m_elements, partitioning the bounding box, ...)
// but is it worth?
// especially with the segmentation of the layers ?
bool MaterialAccountingGroup::inside(const MaterialAccountingDetector& detector) const {
  const GlobalPoint& position = detector.position();
  // first check to see if the point is inside the bounding box
  LogTrace("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())) {
    LogTrace("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 {
    // now check if the point is actually close enough to any element
    LogTrace("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) {
      LogTrace("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 MaterialAccountingGroup::addDetector(const MaterialAccountingDetector& detector) {
  if (not inside(detector))
    return false;

  // multiple hits in the same layer (overlaps, etc.) from a single track still count as one for averaging,
  // since the energy deposits from the track have been already split between the different detectors
  m_buffer += detector.material();
  m_counted = true;

  return true;
}

void MaterialAccountingGroup::endOfTrack(void) {
  // add a detector
  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 MaterialAccountingGroup::savePlot(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(), "");

  // store te plot into m_file
  plot->SetDirectory(m_file);
}

void MaterialAccountingGroup::savePlot(TProfile* plot, float average, const std::string& name) {
  // Nota Bene:
  // these "line" plots are not deleted explicitly since
  //   - deleting them before saving them to a TFile will not save them
  //   - deleting them after the TFile they're stored into results in a SEGV
  // ROOT is probably "taking care" (read: messing things up) somehow...
  TH1F* line =
      new 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(), "");

  // store te plots into m_file
  plot->SetDirectory(m_file);
  line->SetDirectory(m_file);
}

std::string 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();
}

void 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;
}