plugins/dd4hep/DD4hep_TrackingMaterialAnalyser.cc

0001 #include <iostream>  // FIXME: switch to MessagLogger & friends
0002 #include <fstream>
0003 #include <sstream>
0004 #include <vector>
0005 #include <string>
0006 #include <stdexcept>
0007 #include <cstring>
0008 #include <cstdlib>
0009 #include <fmt/printf.h>
0010
0011 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0012 #include "FWCore/Framework/interface/Event.h"
0013 #include "FWCore/Framework/interface/EventSetup.h"
0014 #include "FWCore/Framework/interface/ESTransientHandle.h"
0015 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0016 #include "FWCore/Utilities/interface/EDMException.h"
0017
0018 #include "DetectorDescription/DDCMS/interface/DDCompactView.h"
0019 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0020
0021 #include "SimDataFormats/ValidationFormats/interface/MaterialAccountingStep.h"
0022 #include "SimDataFormats/ValidationFormats/interface/MaterialAccountingTrack.h"
0023 #include "DD4hep_TrackingMaterialAnalyser.h"
0024 #include "DD4hep_TrackingMaterialPlotter.h"
0025
0026 //-------------------------------------------------------------------------
0027 DD4hep_TrackingMaterialAnalyser::DD4hep_TrackingMaterialAnalyser(const edm::ParameterSet& iPSet) {
0028   m_materialToken =
0029       consumes<std::vector<MaterialAccountingTrack> >(iPSet.getParameter<edm::InputTag>("MaterialAccounting"));
0030   m_groupNames = iPSet.getParameter<std::vector<std::string> >("Groups");
0031   const std::string& splitmode = iPSet.getParameter<std::string>("SplitMode");
0032   if (strcasecmp(splitmode.c_str(), "NearestLayer") == 0) {
0033     m_splitMode = NEAREST_LAYER;
0034   } else if (strcasecmp(splitmode.c_str(), "InnerLayer") == 0) {
0035     m_splitMode = INNER_LAYER;
0036   } else if (strcasecmp(splitmode.c_str(), "OuterLayer") == 0) {
0037     m_splitMode = OUTER_LAYER;
0038   } else {
0039     m_splitMode = UNDEFINED;
0040     throw edm::Exception(edm::errors::LogicError)
0041         << "Invalid SplitMode \"" << splitmode
0042         << "\". Acceptable values are \"NearestLayer\", \"InnerLayer\", \"OuterLayer\".";
0043   }
0044   m_dddToken = esConsumes(iPSet.getParameter<edm::ESInputTag>("DDDetector"));
0045   m_skipAfterLastDetector = iPSet.getParameter<bool>("SkipAfterLastDetector");
0046   m_skipBeforeFirstDetector = iPSet.getParameter<bool>("SkipBeforeFirstDetector");
0047   m_saveSummaryPlot = iPSet.getParameter<bool>("SaveSummaryPlot");
0048   m_saveDetailedPlots = iPSet.getParameter<bool>("SaveDetailedPlots");
0049   m_saveParameters = iPSet.getParameter<bool>("SaveParameters");
0050   m_saveXml = iPSet.getParameter<bool>("SaveXML");
0051   m_isHGCal = iPSet.getParameter<bool>("isHGCal");
0052   m_isHFNose = iPSet.getParameter<bool>("isHFNose");
0053   if (m_saveSummaryPlot) {
0054     if (m_isHGCal) {
0055       m_plotter = std::make_unique<DD4hep_TrackingMaterialPlotter>(550., 300., 10);
0056     } else if (m_isHFNose) {
0057       m_plotter = std::make_unique<DD4hep_TrackingMaterialPlotter>(1200., 350., 10);
0058     } else {
0059       m_plotter = std::make_unique<DD4hep_TrackingMaterialPlotter>(300., 120., 10);
0060     }  // 10x10 points per cm2
0061   } else {
0062     m_plotter = nullptr;
0063   }
0064 }
0065
0066 //-------------------------------------------------------------------------
0067 DD4hep_TrackingMaterialAnalyser::~DD4hep_TrackingMaterialAnalyser(void) {}
0068
0069 //-------------------------------------------------------------------------
0070 void DD4hep_TrackingMaterialAnalyser::saveParameters(const char* name) {
0071   std::ofstream parameters(name);
0072   edm::LogVerbatim("TrackerMaterialAnalysis") << "TrackingMaterialAnalyser" << std::endl;
0073   for (unsigned int i = 0; i < m_groups.size(); ++i) {
0074     DD4hep_MaterialAccountingGroup& layer = *(m_groups[i]);
0075     edm::LogVerbatim("TrackerMaterialAnalysis") << "TrackingMaterialAnalyser" << layer.name() << std::endl;
0076     edm::LogVerbatim("TrackerMaterialAnalysis")
0077         << "TrackingMaterialAnalyser" << fmt::sprintf("\tnumber of hits:               %9d", layer.tracks())
0078         << std::endl;
0079     edm::LogVerbatim("TrackerMaterialAnalysis")
0080         << "TrackingMaterialAnalyser"
0081         << fmt::sprintf("\tnormalized segment length:    %9.1f ± %9.1f cm", layer.averageLength(), layer.sigmaLength())
0082         << std::endl;
0083     edm::LogVerbatim("TrackerMaterialAnalysis") << "TrackingMaterialAnalyser"
0084                                                 << fmt::sprintf("\tnormalized radiation lengths: %9.3f ± %9.3f",
0085                                                                 layer.averageRadiationLengths(),
0086                                                                 layer.sigmaRadiationLengths())
0087                                                 << std::endl;
0088     edm::LogVerbatim("TrackerMaterialAnalysis")
0089         << "TrackingMaterialAnalyser"
0090         << fmt::sprintf("\tnormalized energy loss:       %6.5fe-03 ± %6.5fe-03 GeV",
0091                         layer.averageEnergyLoss(),
0092                         layer.sigmaEnergyLoss())
0093         << std::endl;
0094     parameters << fmt::sprintf("%-20s\t%7d\t%5.1f ± %5.1f cm\t%6.4f ± %6.4f \t%6.4fe-03 ± %6.4fe-03 GeV",
0095                                layer.name(),
0096                                layer.tracks(),
0097                                layer.averageLength(),
0098                                layer.sigmaLength(),
0099                                layer.averageRadiationLengths(),
0100                                layer.sigmaRadiationLengths(),
0101                                layer.averageEnergyLoss(),
0102                                layer.sigmaEnergyLoss())
0103                << std::endl;
0104   }
0105   edm::LogVerbatim("TrackerMaterialAnalysis") << "TrackingMaterialAnalyser" << std::endl;
0106
0107   parameters.close();
0108 }
0109
0110 //-------------------------------------------------------------------------
0111 void DD4hep_TrackingMaterialAnalyser::saveXml(const char* name) {
0112   std::ofstream xml(name);
0113   xml << "<?xml version=\"1.0\" encoding=\"utf-8\"?>" << std::endl;
0114   xml << "<Groups>" << std::endl;
0115   for (unsigned int i = 0; i < m_groups.size(); ++i) {
0116     DD4hep_MaterialAccountingGroup& layer = *(m_groups[i]);
0117     xml << "  <Group name=\"" << layer.name() << "\">\n"
0118         << "    <Parameter name=\"TrackerRadLength\" value=\"" << layer.averageRadiationLengths() << "\"/>\n"
0119         << "    <Parameter name=\"TrackerXi\" value=\"" << layer.averageEnergyLoss() << "\"/>\n"
0120         << "  </Group>\n"
0121         << std::endl;
0122   }
0123   xml << "</Groups>" << std::endl;
0124 }
0125
0126 //-------------------------------------------------------------------------
0127 void DD4hep_TrackingMaterialAnalyser::saveLayerPlots() {
0128   for (unsigned int i = 0; i < m_groups.size(); ++i) {
0129     DD4hep_MaterialAccountingGroup& layer = *(m_groups[i]);
0130     layer.savePlots();
0131   }
0132 }
0133
0134 //-------------------------------------------------------------------------
0135 void DD4hep_TrackingMaterialAnalyser::endJob(void) {
0136   if (m_saveParameters)
0137     saveParameters("parameters");
0138
0139   if (m_saveXml)
0140     saveXml("parameters.xml");
0141
0142   if (m_saveDetailedPlots)
0143     saveLayerPlots();
0144
0145   if (m_saveSummaryPlot and m_plotter) {
0146     m_plotter->normalize();
0147     m_plotter->draw();
0148   }
0149 }
0150
0151 //-------------------------------------------------------------------------
0152
0153 //-------------------------------------------------------------------------
0154 void DD4hep_TrackingMaterialAnalyser::analyze(const edm::Event& event, const edm::EventSetup& setup) {
0155   using namespace edm;
0156   auto hDDD = setup.getTransientHandle(m_dddToken);
0157
0158   m_groups.reserve(m_groupNames.size());
0159   // Initialize m_groups iff it has size equal to zero, so that we are
0160   // sure it will never be repopulated with the same entries over and
0161   // over again in the eventloop, at each call of the analyze method.
0162   if (m_groups.empty()) {
0163     for (unsigned int i = 0; i < m_groupNames.size(); ++i)
0164       m_groups.emplace_back(new DD4hep_MaterialAccountingGroup(m_groupNames[i], *hDDD));
0165
0166     edm::LogVerbatim("TrackerMaterialAnalysis")
0167         << "TrackingMaterialAnalyser: List of the tracker groups: " << std::endl;
0168     for (unsigned int i = 0; i < m_groups.size(); ++i)
0169       edm::LogVerbatim("TrackerMaterialAnalysis")
0170           << i << " TrackingMaterialAnalyser:\t" << m_groups[i]->info() << std::endl;
0171   }
0172   Handle<std::vector<MaterialAccountingTrack> > h_tracks;
0173   event.getByToken(m_materialToken, h_tracks);
0174
0175   for (std::vector<MaterialAccountingTrack>::const_iterator t = h_tracks->begin(), end = h_tracks->end(); t != end;
0176        ++t) {
0177     MaterialAccountingTrack track(*t);
0178     split(track);
0179   }
0180 }
0181
0182 //-------------------------------------------------------------------
0183 // split a track in segments, each associated to a sensitive detector
0184 // in a DetLayer; then, associate each step to one segment, splitting
0185 // the steps across the segment boundaries
0186 //
0187 // Nota Bene: this implementation assumes that the steps stored along
0188 // each track are consecutive and adjacent, and that no step can span
0189 // across 3 layers, since all steps should split at layer boundaries
0190
0191 void DD4hep_TrackingMaterialAnalyser::split(MaterialAccountingTrack& track) {
0192   using namespace edm;
0193   // group sensitive detectors by their DetLayer
0194   std::vector<int> group(track.detectors().size());
0195   for (unsigned int i = 0; i < track.detectors().size(); ++i)
0196     group[i] = findLayer(track.detectors()[i]);
0197
0198   for (unsigned int i = 0; i < group.size(); ++i)
0199     if (group[i] > 0)
0200       edm::LogVerbatim("TrackerMaterialAnalysis") << "TrackingMaterialAnalyser:\n"
0201                                                   << "For detector i: " << i << " index: " << group[i]
0202                                                   << " R-ranges: " << m_groups[group[i] - 1]->getBoundingR().first
0203                                                   << ", " << m_groups[group[i] - 1]->getBoundingR().second << group[i]
0204                                                   << " Z-ranges: " << m_groups[group[i] - 1]->getBoundingZ().first
0205                                                   << ", " << m_groups[group[i] - 1]->getBoundingZ().second << std::endl;
0206
0207   unsigned int detectors = track.detectors().size();
0208   if (detectors == 0) {
0209     // the track doesn't cross any active detector:
0210     // keep al material as unassigned
0211     if (m_plotter)
0212       for (unsigned int i = 1; i < track.steps().size(); ++i)
0213         m_plotter->plotSegmentUnassigned(track.steps()[i]);
0214   } else {
0215     const double TOLERANCE = 0.0001;  // 1 um tolerance
0216     std::vector<double> limits(detectors + 2);
0217
0218     // define the trivial limits
0219     if (m_skipBeforeFirstDetector)
0220       limits[0] = track.detectors()[0].m_curvilinearIn - TOLERANCE;
0221     else
0222       limits[0] = -TOLERANCE;
0223     if (m_skipAfterLastDetector)
0224       limits[detectors] = track.detectors()[detectors - 1].m_curvilinearOut + TOLERANCE;
0225     else
0226       limits[detectors] = track.summary().length() + TOLERANCE;
0227     limits[detectors + 1] = INFINITY;  // this is probably no more needed, but doesn't harm...
0228
0229     // pick the algorithm to define the non-trivial limits
0230     switch (m_splitMode) {
0231       // assign each segment to the the nearest layer
0232       // e.g. the material between pixel barrel 3 and TIB 1 will be split among the two
0233       case NEAREST_LAYER:
0234         for (unsigned int i = 1; i < detectors; ++i) {
0235           limits[i] = (track.detectors()[i - 1].m_curvilinearOut + track.detectors()[i].m_curvilinearIn) / 2.;
0236         }
0237         break;
0238
0239       // assign each segment to the the inner layer
0240       // e.g. all material between pixel barrel 3 and TIB 1 will go into the pixel barrel
0241       case INNER_LAYER:
0242         for (unsigned int i = 1; i < detectors; ++i)
0243           limits[i] = track.detectors()[i].m_curvilinearIn - TOLERANCE;
0244         break;
0245
0246       // assign each segment to the the outer layer
0247       // e.g. all material between pixel barrel 3 and TIB 1 will go into the TIB
0248       case OUTER_LAYER:
0249         for (unsigned int i = 1; i < detectors; ++i)
0250           limits[i] = track.detectors()[i - 1].m_curvilinearOut + TOLERANCE;
0251         break;
0252
0253       case UNDEFINED:
0254         [[fallthrough]];
0255
0256       default:
0257         // throw something
0258         throw edm::Exception(edm::errors::LogicError) << "Invalid SplitMode";
0259     }
0260
0261     double begin = 0.;   // beginning of step, along the track
0262     double end = 0.;     // end of step, along the track
0263     unsigned int i = 1;  // step conter
0264
0265     // skip the material before the first layer
0266     while (end < limits[0]) {
0267       const MaterialAccountingStep& step = track.steps()[i++];
0268       end = begin + step.length();
0269
0270       // do not account material before the first layer
0271       if (m_plotter)
0272         m_plotter->plotSegmentUnassigned(step);
0273
0274       begin = end;
0275     }
0276
0277     unsigned int index = 0;  // which detector
0278     while (i < track.steps().size()) {
0279       const MaterialAccountingStep& step = track.steps()[i++];
0280
0281       end = begin + step.length();
0282
0283       if (begin > limits[detectors]) {
0284         // segment after last layer and skipping requested in configuation
0285         if (m_plotter)
0286           m_plotter->plotSegmentUnassigned(step);
0287         begin = end;
0288         continue;
0289       }
0290
0291       // from here onwards we should be in the accountable region, either completely in a single layer:
0292       //   limits[index] <= begin < end <= limits[index+1]
0293       // or possibly split between 2 layers
0294       //   limits[index] < begin < limits[index+1] < end <  limits[index+2]
0295       if (begin < limits[index] or end > limits[index + 2]) {
0296         // sanity check
0297         edm::LogError("TrackerMaterialAnalysis")
0298             << "TrackingMaterialAnalyser::split(): ERROR: internal logic error, expected " << limits[index] << " < "
0299             << begin << " < " << limits[index + 1] << std::endl;
0300         break;
0301       }
0302
0303       if (limits[index] <= begin and end <= limits[index + 1]) {
0304         // step completely inside current detector range
0305         track.detectors()[index].account(step, begin, end);
0306         if (m_plotter)
0307           m_plotter->plotSegmentInLayer(step, group[index]);
0308       } else {
0309         // step shared beteewn two detectors, transition at limits[index+1]
0310         double fraction = (limits[index + 1] - begin) / (end - begin);
0311         assert(fraction < 1.);
0312         std::pair<MaterialAccountingStep, MaterialAccountingStep> parts = step.split(fraction);
0313
0314         if (m_plotter) {
0315           if (index > 0)
0316             m_plotter->plotSegmentInLayer(parts.first, group[index]);
0317           else
0318             // track outside acceptance, keep as unassocated
0319             m_plotter->plotSegmentUnassigned(parts.first);
0320
0321           if (index + 1 < detectors)
0322             m_plotter->plotSegmentInLayer(parts.second, group[index + 1]);
0323           else
0324             // track outside acceptance, keep as unassocated
0325             m_plotter->plotSegmentUnassigned(parts.second);
0326         }
0327
0328         track.detectors()[index].account(parts.first, begin, limits[index + 1]);
0329         ++index;  // next layer
0330         if (index < detectors)
0331           track.detectors()[index].account(parts.second, limits[index + 1], end);
0332       }
0333       begin = end;
0334     }
0335   }
0336
0337   // add the material from each detector to its layer (if there is one and only one)
0338   for (unsigned int i = 0; i < track.detectors().size(); ++i)
0339     if (group[i] != 0)
0340       m_groups[group[i] - 1]->addDetector(track.detectors()[i]);
0341
0342   // end of track: commit internal buffers and reset the m_groups internal state for a new track
0343   for (unsigned int i = 0; i < m_groups.size(); ++i)
0344     m_groups[i]->endOfTrack();
0345 }
0346
0347 //-------------------------------------------------------------------------
0348 // find the layer index (0: none, 1-3: PixelBarrel,
0349 //                       4-7: TID, 8-13: TOB, 14-15,28-29: PixelEndcap,
0350 //                       16-18,30-32: TID, 19-27,33-41: TEC)
0351 int DD4hep_TrackingMaterialAnalyser::findLayer(const MaterialAccountingDetector& detector) {
0352   int index = 0;
0353   size_t inside = 0;
0354   for (size_t i = 0; i < m_groups.size(); ++i)
0355     if (m_groups[i]->isInside(detector)) {
0356       ++inside;
0357       index = i + 1;
0358     }
0359   if (inside == 0) {
0360     index = 0;
0361     edm::LogError("TrackerMaterialAnalysis")
0362         << "TrackingMaterialAnalyser::findLayer(...): ERROR: detector does not belong to any DetLayer" << std::endl;
0363     edm::LogError("TrackerMaterialAnalysis")
0364         << "TrackingMaterialAnalyser::findLayer(...): detector position: " << std::fixed
0365         << " (r: " << std::setprecision(1) << std::setw(5) << detector.position().perp()
0366         << ", z: " << std::setprecision(1) << std::setw(6) << detector.position().z()
0367         << ", phi: " << std::setprecision(3) << std::setw(6) << detector.position().phi() << ")" << std::endl;
0368   }
0369   if (inside > 1) {
0370     index = 0;
0371     edm::LogError("TrackerMaterialAnalysis") << "TrackingMaterialAnalyser::findLayer(...): ERROR: detector belongs to "
0372                                              << inside << " DetLayers" << std::endl;
0373     edm::LogError("TrackerMaterialAnalysis")
0374         << "TrackingMaterialAnalyser::findLayer(...): detector position: " << std::fixed
0375         << " (r: " << std::setprecision(1) << std::setw(5) << detector.position().perp()
0376         << ", z: " << std::setprecision(1) << std::setw(6) << detector.position().z()
0377         << ", phi: " << std::setprecision(3) << std::setw(6) << detector.position().phi() << ")" << std::endl;
0378   }
0379
0380   return index;
0381 }
0382
0383 //-------------------------------------------------------------------------
0384 // define as a plugin
0385 #include "FWCore/PluginManager/interface/ModuleDef.h"
0386 #include "FWCore/Framework/interface/MakerMacros.h"
0387 DEFINE_FWK_MODULE(DD4hep_TrackingMaterialAnalyser);