Improve reco material code

SimTracker/TrackerMaterialAnalysis/plugins/ListGroups.cc

@@ -173,16 +173,19 @@ void ListGroups::fillGradient(void)
   m_gradient.push_back(kRed); // Overflow highest bin
 }
 
-void ListGroups::fillColor(void)
-{
+void 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
+  // the real number of grouping used in the different scenarios. We therefore
+  // define some reasonable set and loop over it in case the number of grouping
+  // is larger than the number of colors.
+
   m_color.push_back(kBlack);          // unassigned
 
   m_color.push_back(kAzure);          // PixelBarrelLayer0_Z0
   m_color.push_back(kAzure - 1);      // PixelBarrelLayer0_Z20
   m_color.push_back(kAzure + 1) ;     // Layer1_Z0
   m_color.push_back(kAzure + 2) ;     // Layer1_Z20
-//  m_color.push_back(kAzure + 3) ;     // Layer2_Z0
-//  m_color.push_back(kAzure + 10);     // Layer2_Z15
 
   m_color.push_back(kGreen);          // EndCapDisk1_R0
   m_color.push_back(kGreen + 2);      // EndcapDisk1_R11
@@ -357,7 +360,6 @@ void ListGroups::produceAndSaveSummaryPlot(const edm::EventSetup &setup) {
       new TProfile2D( "OverallDifferencesEnergyLoss", "OverallDifferencesEnergyLoss",
                       600., -300., 300, 120., 0., 120.));
 
-//  assert(m_color.size() - 1 == m_groups.size());
   for (auto g : m_groups) {
     m_plots.push_back(
         new TH2F( g->name().c_str(), g->name().c_str(),

SimTracker/TrackerMaterialAnalysis/plugins/ListIds.cc

@@ -1,6 +1,7 @@
 #include <string>
 #include <vector>
 #include <iostream>
+#include <algorithm>
 
 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
 #include "FWCore/Framework/interface/Event.h"
@@ -62,7 +63,7 @@ std::string dddGetString(const std::string & s, const DDFilteredView & view) {
     return std::string();
 }
 
-static inline 
+static inline
 std::ostream & operator<<(std::ostream & out, const math::XYZVector & v) {
   return out << "(" << v.rho() << ", " << v.z() << ", " << v.phi() << ")";
 }
@@ -72,14 +73,22 @@ class ListIds : public edm::one::EDAnalyzer<>
 public:
   ListIds(const edm::ParameterSet &);
   virtual ~ListIds();
-  
+
 private:
   void analyze(const edm::Event &, const edm::EventSetup &) override;
   void beginJob() override {}
   void endJob() override;
+  // List of material names used to select specific detectors.
+  // Names are matched literally, w/o any usage of regexp.
+  // Names should also be specified with the correct namespace,
+  // otherwise the matching will fail.
+  bool printMaterial_;
+  std::vector<std::string> materials_;
 };
 
-ListIds::ListIds(const edm::ParameterSet &) {
+ListIds::ListIds(const edm::ParameterSet & pset)
+ : printMaterial_(pset.getUntrackedParameter<bool>("printMaterial")),
+   materials_(pset.getUntrackedParameter<std::vector<std::string> >("materials")) {
 }
 
 ListIds::~ListIds() {
@@ -95,34 +104,44 @@ ListIds::analyze(const edm::Event& evt, const edm::EventSetup& setup){
   CmsTrackerStringToEnum theCmsTrackerStringToEnum;
   DDSpecificsFilter filter;
   filter.setCriteria(DDValue(attribute, "any", 0), DDCompOp::not_equals);
-  DDFilteredView fv(*hDdd); 
+  DDFilteredView fv(*hDdd);
   fv.addFilter(filter);
   if (theCmsTrackerStringToEnum.type(dddGetString(attribute, fv)) != GeometricDet::Tracker) {
     fv.firstChild();
     if (theCmsTrackerStringToEnum.type(dddGetString(attribute, fv)) != GeometricDet::Tracker)
       throw cms::Exception("Configuration") << "The first child of the DDFilteredView is not what is expected \n"
                                             << dddGetString(attribute, fv);
   }
-   
+
   std::cout << std::fixed << std::setprecision(3);
+  bool printAnyMaterial = (std::find(materials_.begin(), materials_.end(), "ANY") != materials_.end());
   do {
-    // print the full hierarchy of all Silicon items
-    if (fv.logicalPart().material().name() == "materials:Silicon") {
-
+    // print the full hierarchy of all elements whose material
+    // has been specified by the user. An empty list of
+    // materials will print no elements. The special
+    // keyword ANY (in any location of the vector)
+    // will select all elements.
+    if (printAnyMaterial ||
+        (std::find(materials_.begin(),
+                   materials_.end(),
+                   fv.logicalPart().material().name().fullname()) != materials_.end())) {
+
       // start from 2 to skip the leading /OCMS[0]/CMSE[1] part
       const DDGeoHistory & history = fv.geoHistory();
       std::cout << '/';
       for (unsigned int h = 2; h < history.size(); ++h) {
-        std::cout << '/' << history[h].logicalPart().name().name() << '[' << history[h].copyno() << ']';
+        std::cout << '/' << history[h].logicalPart().name().ns()
+                  << ":" << history[h].logicalPart().name().name() << '[' << history[h].copyno() << ']';
       }
-
+      if (printMaterial_)
+        std::cout << " Material: |" << fv.logicalPart().material().name() << "|";
       // DD3Vector and DDTranslation are the same type as math::XYZVector
       math::XYZVector position = fv.translation() / 10.;  // mm -> cm
       std::cout << "\t" << position << std::endl;
     }
   } while (fv.next());
   std::cout << std::endl;
-  
+
   std::cout << "______________________________ std::vector<GeomDet*> from TrackerGeometry::dets() ______________________________" << std::endl;
   edm::ESHandle<TrackerGeometry> hGeo;
   setup.get<TrackerDigiGeometryRecord>().get( hGeo );
@@ -136,8 +155,8 @@ ListIds::analyze(const edm::Event& evt, const edm::EventSetup& setup){
     const Surface::PositionType & p = det.position();
     math::XYZVector position(p.x(), p.y(), p.z());
 
-    std::cout << det.subDetector() << '\t' 
-              << det.geographicalId().det() << '\t' 
+    std::cout << det.subDetector() << '\t'
+              << det.geographicalId().det() << '\t'
               << det.geographicalId().subdetId() << '\t'
               << det.geographicalId().rawId() << "\t"
               << position;
@@ -150,7 +169,7 @@ ListIds::analyze(const edm::Event& evt, const edm::EventSetup& setup){
     }
     std::cout << std::endl;
   }
-    
+
 }
 
 void

SimTracker/TrackerMaterialAnalysis/plugins/MaterialAccountingGroup.cc

@@ -34,9 +34,14 @@ MaterialAccountingGroup::MaterialAccountingGroup( const std::string & name, cons
   DDSpecificsFilter filter;
   filter.setCriteria(DDValue("TrackingMaterialGroup", name), DDCompOp::equals);
   fv.addFilter(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;
     m_elements.push_back( GlobalPoint(position.x(), position.y(), position.z()) );
   }
 
@@ -45,8 +50,14 @@ MaterialAccountingGroup::MaterialAccountingGroup( const std::string & name, cons
     m_boundingbox.grow(m_elements[i].perp(), m_elements[i].z());
   }
   m_boundingbox.grow(s_tolerance);
-
-  // initialize the histograms 
+  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);
@@ -79,24 +90,47 @@ MaterialAccountingGroup::~MaterialAccountingGroup(void)
 
 // TODO the inner check could be sped up in many ways
 // (sorting the m_elements, partitioning the bounding box, ...)
-// but is it worth? 
+// 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
-    for (unsigned int i = 0; i < m_elements.size(); ++i)
+    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 ) 
+bool MaterialAccountingGroup::addDetector( const MaterialAccountingDetector& detector )
 {
   if (not inside(detector))
     return false;
@@ -117,7 +151,7 @@ void MaterialAccountingGroup::endOfTrack(void) {
     ++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().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() );

SimTracker/TrackerMaterialAnalysis/plugins/TrackingMaterialAnalyser.cc

@@ -136,7 +136,8 @@ void TrackingMaterialAnalyser::endJob(void)
 //-------------------------------------------------------------------------
 void TrackingMaterialAnalyser::analyze(const edm::Event& event, const edm::EventSetup& setup)
 {
-  edm::ESTransientHandle<DDCompactView> hDDD;
+  using namespace edm;
+  ESTransientHandle<DDCompactView> hDDD;
   setup.get<IdealGeometryRecord>().get( hDDD );
 
   m_groups.reserve( m_groupNames.size() );
@@ -147,13 +148,13 @@ void TrackingMaterialAnalyser::analyze(const edm::Event& event, const edm::Event
     for (unsigned int i = 0; i < m_groupNames.size(); ++i)
       m_groups.push_back( new MaterialAccountingGroup( m_groupNames[i], * hDDD) );
 
-    LogDebug("TrackingMaterialAnalyser")
+    LogInfo("TrackingMaterialAnalyser")
         << "TrackingMaterialAnalyser: List of the tracker groups: " << std::endl;
     for (unsigned int i = 0; i < m_groups.size(); ++i)
-      LogDebug("TrackingMaterialAnalyser")
+      LogInfo("TrackingMaterialAnalyser")
         << i << " TrackingMaterialAnalyser:\t" << m_groups[i]->info() << std::endl;
   }
-  edm::Handle< std::vector<MaterialAccountingTrack> > h_tracks;
+  Handle< std::vector<MaterialAccountingTrack> > h_tracks;
   event.getByToken(m_materialToken, h_tracks);
 
   for (std::vector<MaterialAccountingTrack>::const_iterator t = h_tracks->begin(), end = h_tracks->end(); t != end; ++t) {

SimTracker/TrackerMaterialAnalysis/plugins/TrackingMaterialProducer.cc

@@ -30,11 +30,12 @@
 #include "TrackingMaterialProducer.h"
 
 // Uncomment the following #define directive to have the full list of
-// volumes known to G4 printed to LogDebug("TrackingMaterialProducer")
+// volumes known to G4 printed to LogInfo("TrackingMaterialProducer")
 
 #define DEBUG_G4_VOLUMES
 
 using namespace CLHEP;
+using edm::LogInfo;
 
 // missing from GEANT4 < 9.0 : G4LogicalVolumeStore::GetVolume( name )
 static
@@ -43,7 +44,7 @@ const G4LogicalVolume* GetVolume(const std::string& name) {
 
 #ifdef DEBUG_G4_VOLUMES
   for (G4LogicalVolumeStore::const_iterator volume = lvs->begin(); volume != lvs->end(); ++volume)
-    LogDebug("TrackingMaterialProducer") << "TrackingMaterialProducer: G4 registered volumes "
+    LogInfo("TrackingMaterialProducer") << "TrackingMaterialProducer: G4 registered volumes "
                                          << (*volume)->GetName() << std::endl;
 #endif
 
@@ -108,12 +109,12 @@ void TrackingMaterialProducer::update(const EndOfJob* event)
 void TrackingMaterialProducer::update(const BeginOfJob* event)
 {
   // INFO
-  LogDebug("TrackingMaterialProducer") << "TrackingMaterialProducer: List of the selected volumes: " << std::endl;
+  LogInfo("TrackingMaterialProducer") << "TrackingMaterialProducer: List of the selected volumes: " << std::endl;
   for (std::vector<std::string>::const_iterator volume_name = m_selectedNames.begin();
        volume_name != m_selectedNames.end(); ++volume_name) {
     const G4LogicalVolume* volume = GetVolume(*volume_name);
     if (volume) {
-      LogDebug("TrackingMaterialProducer") << "TrackingMaterialProducer: " << *volume_name << std::endl;
+      LogInfo("TrackingMaterialProducer") << "TrackingMaterialProducer: " << *volume_name << std::endl;
       m_selectedVolumes.push_back( volume );
     } else {
       // FIXME: throw an exception ?
@@ -160,7 +161,7 @@ void TrackingMaterialProducer::update(const G4Step* step)
 {
   const G4TouchableHistory* touchable = (G4TouchableHistory*)(step->GetTrack()->GetTouchable());
   if (not isSelectedFast( touchable )) {
-    LogDebug("TrackingMaterialProducer") << "TrackingMaterialProducer:\t[...] skipping "
+    LogInfo("TrackingMaterialProducer") << "TrackingMaterialProducer:\t[...] skipping "
                                          << touchable->GetVolume()->GetName() << std::endl;
     return;
   }
@@ -224,23 +225,30 @@ void TrackingMaterialProducer::update(const G4Step* step)
     m_track.leaveDetector( sensitive, cosThetaPost );
 
   if (sensitive)
-    LogDebug("TrackingMaterialProducer") << "Track was near sensitive     volume "
+    LogInfo("TrackingMaterialProducer") << "Track was near sensitive     volume "
                                          << sensitive->GetName() << std::endl;
   else
-    LogDebug("TrackingMaterialProducer") << "Track was near non-sensitive volume "
+    LogInfo("TrackingMaterialProducer") << "Track was near non-sensitive volume "
                                          << touchable->GetVolume()->GetName() << std::endl;
-  LogDebug("TrackingMaterialProducer") << "Step length:             "
-                                       << length << " cm" << std::endl;
-  LogDebug("TrackingMaterialProducer") << "Radiation lengths:       "
+  LogInfo("TrackingMaterialProducer")  << "Step length:             "
+                                       << length << " cm\n"
+                                       << "globalPreStep(r,z): (" << globalPositionIn.perp()
+                                       << ", " << globalPositionIn.z() << ") cm\n"
+                                       << "globalPostStep(r,z): (" << globalPositionOut.perp()
+                                       << ", " << globalPositionOut.z() << ") cm\n"
+                                       << "position(r,z): ("
+                                       << position.perp()
+                                       << ", " << position.z() << ") cm\n"
+                                       << "Radiation lengths:       "
                                        << radiationLengths << " \t\t(X0: "
-                                       << X0 << " cm)" << std::endl;
-  LogDebug("TrackingMaterialProducer") << "Energy loss:             "
+                                       << X0 << " cm)\n"
+                                       << "Energy loss:             "
                                        << energyLoss << " MeV  \t(Xi: "
-                                       << Xi << " MeV/cm)" << std::endl;
-  LogDebug("TrackingMaterialProducer") << "Track was " << (enter_sensitive ? "entering " : "in none ")
-                                       << "sensitive volume" << std::endl;
-  LogDebug("TrackingMaterialProducer") << "Track was " << (leave_sensitive ? "leaving  " : "in none ")
-                                       << "sensitive volume" << std::endl;
+                                       << Xi << " MeV/cm)\n"
+                                       << "Track was " << (enter_sensitive ? "entering " : "in none ")
+                                       << "sensitive volume\n"
+                                       << "Track was " << (leave_sensitive ? "leaving  " : "in none ")
+                                       << "sensitive volume\n";
 
 }
 
@@ -255,18 +263,18 @@ void TrackingMaterialProducer::update(const EndOfTrack* event)
   radLen_vs_eta_->Fill(track->GetMomentum().eta(), m_track.summary().radiationLengths());
   m_tracks->push_back(m_track);
 
-  // LogDebug
-  LogDebug("TrackingMaterialProducer") << "TrackingMaterialProducer: this track took "
+  // LogInfo
+  LogInfo("TrackingMaterialProducer") << "TrackingMaterialProducer: this track took "
                                        << m_track.steps().size()
                                        << " steps, and passed through "
                                        << m_track.detectors().size()
                                        << " sensitive detectors" << std::endl;
-  LogDebug("TrackingMaterialProducer") << "TrackingMaterialProducer: track length:       "
+  LogInfo("TrackingMaterialProducer") << "TrackingMaterialProducer: track length:       "
                                        << m_track.summary().length()
                                        << " cm" << std::endl;
-  LogDebug("TrackingMaterialProducer") << "TrackingMaterialProducer: radiation lengths: "
+  LogInfo("TrackingMaterialProducer") << "TrackingMaterialProducer: radiation lengths: "
                                        << m_track.summary().radiationLengths() << std::endl;
-  LogDebug("TrackingMaterialProducer") << "TrackingMaterialProducer: energy loss:        "
+  LogInfo("TrackingMaterialProducer") << "TrackingMaterialProducer: energy loss:        "
                                        << m_track.summary().energyLoss()
                                        << " MeV" << std::endl;
 }

SimTracker/TrackerMaterialAnalysis/python/trackingMaterialAnalyser_ForPhaseII_cff.py

@@ -0,0 +1,8 @@
+import FWCore.ParameterSet.Config as cms
+
+# Define arbitrary tracker material groups
+from RecoTracker.GeometryESProducer.TrackerRecoGeometryESProducer_cfi import *
+from SimTracker.TrackerMaterialAnalysis.trackingMaterialGroups_ForPhaseII_cff import *
+
+# Analyze and plot the tracking material
+from SimTracker.TrackerMaterialAnalysis.trackingMaterialAnalyser_ForPhaseII_cfi import *