TrackerMaterialAnalysis/plugins/MaterialAccountingGroup.cc

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;  // 100um should be small enough that no elements from different layers/groups are so close
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   // retrieve the elements from DDD
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     // DD3Vector and DDTranslation are the same type as math::XYZVector
0050     math::XYZVector position = fv.translation() / 10.;  // mm -> cm
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   // grow the bounding box
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   // initialize the histograms
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 // TODO the inner check could be sped up in many ways
0100 // (sorting the m_elements, partitioning the bounding box, ...)
0101 // but is it worth?
0102 // especially with the segmentation of the layers ?
0103 bool MaterialAccountingGroup::inside(const MaterialAccountingDetector& detector) const {
0104   const GlobalPoint& position = detector.position();
0105   // first check to see if the point is inside the bounding box
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     // now check if the point is actually close enough to any element
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   // multiple hits in the same layer (overlaps, etc.) from a single track still count as one for averaging,
0137   // since the energy deposits from the track have been already split between the different detectors
0138   m_buffer += detector.material();
0139   m_counted = true;
0140
0141   return true;
0142 }
0143
0144 void MaterialAccountingGroup::endOfTrack(void) {
0145   // add a detector
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);  // grey
0170   plot->SetLineColor(1);   // black
0171   plot->Draw("c e");
0172   canvas.GetFrame()->SetFillColor(kWhite);
0173   canvas.Draw();
0174   canvas.SaveAs((name + ".png").c_str(), "");
0175
0176   // store te plot into m_file
0177   plot->SetDirectory(m_file);
0178 }
0179
0180 void MaterialAccountingGroup::savePlot(TProfile* plot, float average, const std::string& name) {
0181   // Nota Bene:
0182   // these "line" plots are not deleted explicitly since
0183   //   - deleting them before saving them to a TFile will not save them
0184   //   - deleting them after the TFile they're stored into results in a SEGV
0185   // ROOT is probably "taking care" (read: messing things up) somehow...
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);  // grey
0192   plot->SetLineColor(1);   // black
0193   plot->SetLineWidth(2);
0194   plot->Draw("c e6");
0195   line->SetLineColor(2);  // red
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   // store te plots into m_file
0203   plot->SetDirectory(m_file);
0204   line->SetDirectory(m_file);
0205 }
0206
0207 /// get some infos
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 }