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