#include <TrackingMaterialAnalyser.h>

Inheritance diagram for TrackingMaterialAnalyser:

Public Member Functions
	TrackingMaterialAnalyser (const edm::ParameterSet &)

virtual	~TrackingMaterialAnalyser ()

Public Member Functions inherited from edm::EDAnalyzer
	EDAnalyzer ()

std::string	workerType () const

virtual	~EDAnalyzer ()

Public Member Functions inherited from edm::EDConsumerBase
	EDConsumerBase ()

ProductHolderIndex	indexFrom (EDGetToken, BranchType, TypeID const &) const

void	itemsMayGet (BranchType, std::vector< ProductHolderIndex > &) const

void	itemsToGet (BranchType, std::vector< ProductHolderIndex > &) const

void	labelsForToken (EDGetToken iToken, Labels &oLabels) const

void	updateLookup (BranchType iBranchType, ProductHolderIndexHelper const &)

virtual	~EDConsumerBase ()

Private Types
enum	SplitMode { NEAREST_LAYER, INNER_LAYER, OUTER_LAYER, UNDEFINED }

Private Member Functions
void	analyze (const edm::Event &, const edm::EventSetup &)

void	beginJob (const edm::EventSetup &)

void	endJob ()

int	findLayer (const MaterialAccountingDetector &detector)

void	save ()

void	saveLayerPlots (const char *name)

void	saveParameters (const char *name)

void	saveXml (const char *name)

void	split (MaterialAccountingTrack &track)

Private Attributes
edm::ESWatcher < IdealGeometryRecord >	m_geometryWatcher

std::vector< std::string >	m_groupNames

std::vector < MaterialAccountingGroup * >	m_groups

edm::InputTag	m_material

TrackingMaterialPlotter *	m_plotter

bool	m_saveDetailedPlots

bool	m_saveParameters

bool	m_saveSummaryPlot

bool	m_saveXml

bool	m_skipAfterLastDetector

bool	m_skipBeforeFirstDetector

SplitMode	m_splitMode

Additional Inherited Members
Public Types inherited from edm::EDAnalyzer
typedef EDAnalyzer	ModuleType

typedef WorkerT< EDAnalyzer >	WorkerType

Static Public Member Functions inherited from edm::EDAnalyzer
static const std::string &	baseType ()

static void	fillDescriptions (ConfigurationDescriptions &descriptions)

static void	prevalidate (ConfigurationDescriptions &)

Protected Member Functions inherited from edm::EDAnalyzer
void	callWhenNewProductsRegistered (std::function< void(BranchDescription const &)> const &func)

CurrentProcessingContext const *	currentContext () const

Protected Member Functions inherited from edm::EDConsumerBase
template<typename ProductType , BranchType B = InEvent>
EDGetTokenT< ProductType >	consumes (edm::InputTag const &tag)

EDGetToken	consumes (const TypeToGet &id, edm::InputTag const &tag)

template<BranchType B>
EDGetToken	consumes (TypeToGet const &id, edm::InputTag const &tag)

ConsumesCollector	consumesCollector ()
	Use a ConsumesCollector to gather consumes information from helper functions. More...

template<typename ProductType , BranchType B = InEvent>
void	consumesMany ()

void	consumesMany (const TypeToGet &id)

template<BranchType B>
void	consumesMany (const TypeToGet &id)

template<typename ProductType , BranchType B = InEvent>
EDGetTokenT< ProductType >	mayConsume (edm::InputTag const &tag)

EDGetToken	mayConsume (const TypeToGet &id, edm::InputTag const &tag)

template<BranchType B>
EDGetToken	mayConsume (const TypeToGet &id, edm::InputTag const &tag)

Detailed Description

Definition at line 18 of file TrackingMaterialAnalyser.h.

Member Enumeration Documentation

enum TrackingMaterialAnalyser::SplitMode

private

Enumerator
NEAREST_LAYER
INNER_LAYER
OUTER_LAYER
UNDEFINED

Definition at line 25 of file TrackingMaterialAnalyser.h.

                  {
     NEAREST_LAYER,
     INNER_LAYER,
     OUTER_LAYER,
     UNDEFINED
   };

Constructor & Destructor Documentation

TrackingMaterialAnalyser::TrackingMaterialAnalyser ( const edm::ParameterSet & iPSet )

Definition at line 30 of file TrackingMaterialAnalyser.cc.

References edm::hlt::Exception, edm::ParameterSet::getParameter(), INNER_LAYER, edm::errors::LogicError, m_groupNames, m_material, m_plotter, m_saveDetailedPlots, m_saveParameters, m_saveSummaryPlot, m_saveXml, m_skipAfterLastDetector, m_skipBeforeFirstDetector, m_splitMode, NEAREST_LAYER, NULL, OUTER_LAYER, AlCaHLTBitMon_QueryRunRegistry::string, and UNDEFINED.

 {
   m_material                = iPSet.getParameter<edm::InputTag>("MaterialAccounting");
   m_groupNames              = iPSet.getParameter<std::vector<std::string> >("Groups");
   const std::string & splitmode = iPSet.getParameter<std::string>("SplitMode");
   if (strcasecmp(splitmode.c_str(), "NearestLayer") == 0) {
     m_splitMode = NEAREST_LAYER;
   } else if (strcasecmp(splitmode.c_str(), "InnerLayer") == 0) {
     m_splitMode = INNER_LAYER;
   } else if (strcasecmp(splitmode.c_str(), "OuterLayer") == 0) {
     m_splitMode = OUTER_LAYER;
   } else {
     m_splitMode = UNDEFINED;
     throw edm::Exception(edm::errors::LogicError) << "Invalid SplitMode \"" << splitmode << "\". Acceptable values are \"NearestLayer\", \"InnerLayer\", \"OuterLayer\".";
   }
   m_skipAfterLastDetector   = iPSet.getParameter<bool>("SkipAfterLastDetector");
   m_skipBeforeFirstDetector = iPSet.getParameter<bool>("SkipBeforeFirstDetector");
   m_saveSummaryPlot         = iPSet.getParameter<bool>("SaveSummaryPlot");
   m_saveDetailedPlots       = iPSet.getParameter<bool>("SaveDetailedPlots");
   m_saveParameters          = iPSet.getParameter<bool>("SaveParameters");
   m_saveXml                 = iPSet.getParameter<bool>("SaveXML");
   if (m_saveSummaryPlot)
     m_plotter               = new TrackingMaterialPlotter( 300., 120., 10 );      // 10x10 points per cm2
   else
     m_plotter               = NULL;
 }

TrackingMaterialAnalyser::~TrackingMaterialAnalyser ( void )

virtual

Definition at line 58 of file TrackingMaterialAnalyser.cc.

References m_plotter.

 {
   if (m_plotter)
     delete m_plotter;
 }

Member Function Documentation

void TrackingMaterialAnalyser::analyze	(	const edm::Event &	event,
		const edm::EventSetup &	setup
	)

privatevirtual

Implements edm::EDAnalyzer.

Definition at line 144 of file TrackingMaterialAnalyser.cc.

References edm::ESWatcher< T >::check(), gather_cfg::cout, end, edm::EventSetup::get(), i, info, m_geometryWatcher, m_groupNames, m_groups, m_material, m_plotter, nullptr, TrackingMaterialPlotter::reset(), save(), split(), and lumiQTWidget::t.

 {
   if (m_geometryWatcher.check(setup)) {
     edm::ESTransientHandle<DDCompactView> hDDD;
     setup.get<IdealGeometryRecord>().get( hDDD );
 
     if (m_groups.empty()) {
       // initialise the layers for the first time
       m_groups.resize( m_groupNames.size(), nullptr );
     } else {
       // the geometry has changed: save the current parameters/xml/plots and re-initialise the layers
       save();
       for (unsigned int i = 0; i < m_groups.size(); ++i)
         delete m_groups[i];
       if (m_plotter)
         m_plotter->reset();
     }
 
     for (unsigned int i = 0; i < m_groups.size(); ++i)
       m_groups[i] = new MaterialAccountingGroup( m_groupNames[i], * hDDD);
 
     // INFO
     std::cout << "TrackingMaterialAnalyser: List of the tracker groups: " << std::endl;
     for (unsigned int i = 0; i < m_groups.size(); ++i)
       std::cout << '\t' << m_groups[i]->info() << std::endl;
     std::cout << std::endl;
   }
 
   edm::Handle< std::vector<MaterialAccountingTrack> > h_tracks;
   event.getByLabel(m_material, h_tracks);
 
   for (std::vector<MaterialAccountingTrack>::const_iterator t = h_tracks->begin(), end = h_tracks->end(); t != end; ++t) {
     MaterialAccountingTrack track(*t);
     split( track );
   }
 }

void TrackingMaterialAnalyser::beginJob ( const edm::EventSetup & )

private

void TrackingMaterialAnalyser::endJob ( void )

privatevirtual

Reimplemented from edm::EDAnalyzer.

Definition at line 138 of file TrackingMaterialAnalyser.cc.

References save().

 {
   save();
 }

int TrackingMaterialAnalyser::findLayer ( const MaterialAccountingDetector & detector )

private

Definition at line 364 of file TrackingMaterialAnalyser.cc.

References dtNoiseDBValidation_cfg::cerr, i, getHLTprescales::index, m_groups, PV3DBase< T, PVType, FrameType >::perp(), PV3DBase< T, PVType, FrameType >::phi(), MaterialAccountingDetector::position(), and PV3DBase< T, PVType, FrameType >::z().

Referenced by split().

 {
   int    index  = 0;
   size_t inside = 0;
   for (size_t i = 0; i < m_groups.size(); ++i)
     if (m_groups[i]->inside(detector)) {
       ++inside;
       index = i+1;
     }
   if (inside == 0) {
     index = 0;
     std::cerr << "TrackingMaterialAnalyser::findLayer(...): ERROR: detector does not belong to any DetLayer" << std::endl;
     std::cerr << "TrackingMaterialAnalyser::findLayer(...): detector position: " << std::fixed
               << " (r: " << std::setprecision(1) << std::setw(5) << detector.position().perp()
               << ", z: " << std::setprecision(1) << std::setw(6) << detector.position().z()
               << ", phi: " << std::setprecision(3) << std::setw(6) << detector.position().phi() << ")" 
               << std::endl;
   }
   if (inside > 1) {
     index = 0;
     std::cerr << "TrackingMaterialAnalyser::findLayer(...): ERROR: detector belongs to " << inside << "DetLayers" << std::endl;
     std::cerr << "TrackingMaterialAnalyser::findLayer(...): detector position: " << std::fixed
               << " (r: " << std::setprecision(1) << std::setw(5) << detector.position().perp()
               << ", z: " << std::setprecision(1) << std::setw(6) << detector.position().z()
               << ", phi: " << std::setprecision(3) << std::setw(6) << detector.position().phi() << ")" 
               << std::endl;
   }
 
   return index;
 }

void TrackingMaterialAnalyser::save ( )

private

Definition at line 120 of file TrackingMaterialAnalyser.cc.

References m_plotter, m_saveDetailedPlots, m_saveParameters, m_saveSummaryPlot, m_saveXml, saveLayerPlots(), saveParameters(), and saveXml().

Referenced by analyze(), and endJob().

 {
   if (m_saveParameters)
     saveParameters("parameters");
 
   if (m_saveXml)
     saveXml("parameters.xml");
 
   if (m_saveDetailedPlots)
     saveLayerPlots("layers");
 
   if (m_saveSummaryPlot and m_plotter) {
     m_plotter->normalize();
     m_plotter->draw();
   }
 }

void TrackingMaterialAnalyser::saveLayerPlots ( const char * name )

private

Definition at line 107 of file TrackingMaterialAnalyser.cc.

References i, m_groups, and MaterialAccountingGroup::savePlots().

Referenced by save().

 {
   if (boost::filesystem::is_directory(name))
     boost::filesystem::remove_all(name);
   boost::filesystem::create_directory(name);
   for (unsigned int i = 0; i < m_groups.size(); ++i) {
     MaterialAccountingGroup & layer = *(m_groups[i]);
     layer.savePlots(name);
   }
 }

void TrackingMaterialAnalyser::saveParameters ( const char * name )

private

Definition at line 65 of file TrackingMaterialAnalyser.cc.

References MaterialAccountingGroup::averageEnergyLoss(), MaterialAccountingGroup::averageLength(), MaterialAccountingGroup::averageRadiationLengths(), gather_cfg::cout, diffTreeTool::format(), i, m_groups, MaterialAccountingGroup::name(), Parameters::parameters, MaterialAccountingGroup::sigmaEnergyLoss(), MaterialAccountingGroup::sigmaLength(), MaterialAccountingGroup::sigmaRadiationLengths(), and MaterialAccountingGroup::tracks().

Referenced by save().

 {
   std::ofstream parameters(name);
   std::cout << std::endl;
   for (unsigned int i = 0; i < m_groups.size(); ++i) {
     MaterialAccountingGroup & layer = *(m_groups[i]);
     std::cout << layer.name() << std::endl;
     std::cout << boost::format("\tnumber of hits:               %9d") % layer.tracks() << std::endl;
     std::cout << boost::format("\tnormalized segment length:    %9.1f ± %9.1f cm")  % layer.averageLength()           % layer.sigmaLength()           << std::endl;
     std::cout << boost::format("\tnormalized radiation lengths: %9.3f ± %9.3f")     % layer.averageRadiationLengths() % layer.sigmaRadiationLengths() << std::endl;
     std::cout << boost::format("\tnormalized energy loss:       %9.3f ± %9.3f MeV") % layer.averageEnergyLoss()       % layer.sigmaEnergyLoss()       << std::endl;
     parameters << boost::format("%-20s\t%7d\t%5.1f ± %5.1f cm\t%6.4f ± %6.4f \t%6.4fe-03 ± %6.4fe-03 GeV")
                                 % layer.name()
                                 % layer.tracks()
                                 % layer.averageLength()               % layer.sigmaLength()
                                 % layer.averageRadiationLengths()     % layer.sigmaRadiationLengths()
                                 % layer.averageEnergyLoss()           % layer.sigmaEnergyLoss()
                << std::endl;
   }
   std::cout << std::endl;
 
   parameters.close();
 }

void TrackingMaterialAnalyser::saveXml ( const char * name )

private

Definition at line 90 of file TrackingMaterialAnalyser.cc.

References MaterialAccountingGroup::averageEnergyLoss(), MaterialAccountingGroup::averageRadiationLengths(), i, m_groups, and MaterialAccountingGroup::name().

Referenced by save().

 {
   std::ofstream xml(name);
   xml << "<?xml version=\"1.0\" encoding=\"utf-8\"?>" << std::endl;
   xml << "<Groups>" << std::endl;
   for (unsigned int i = 0; i < m_groups.size(); ++i) {
     MaterialAccountingGroup & layer = *(m_groups[i]);
     xml << "  <Group name=\"" << layer.name() << "\">\n"
         << "    <Parameter name=\"TrackerRadLength\" value=\"" << layer.averageRadiationLengths() << "\"/>\n"
         << "    <Parameter name=\"TrackerXi\" value=\"" << layer.averageEnergyLoss() << "\"/>\n"
         << "  </Group>\n" 
         << std::endl;
   }
   xml << "</Groups>" << std::endl;
 }

void TrackingMaterialAnalyser::split ( MaterialAccountingTrack & track )

private

Definition at line 188 of file TrackingMaterialAnalyser.cc.

References begin, dtNoiseDBValidation_cfg::cerr, MuonGeometrySanityCheck_cfi::detectors(), end, findLayer(), i, getHLTprescales::index, INNER_LAYER, MaterialAccountingStep::length(), edm::errors::LogicError, MaterialAccountingTrack::m_detectors, m_groups, m_plotter, m_skipAfterLastDetector, m_skipBeforeFirstDetector, m_splitMode, MaterialAccountingTrack::m_steps, MaterialAccountingTrack::m_total, NEAREST_LAYER, or, OUTER_LAYER, CfgNavigationSchool_cfi::parts, TrackingMaterialPlotter::plotSegmentInLayer(), TrackingMaterialPlotter::plotSegmentUnassigned(), MaterialAccountingStep::split(), relval_parameters_module::step, and UNDEFINED.

Referenced by analyze().

 {
   // group sensitive detectors by their DetLayer
   std::vector<int> group( track.m_detectors.size() );
   for (unsigned int i = 0; i < track.m_detectors.size(); ++i)
     group[i] = findLayer( track.m_detectors[i] );
 
   unsigned int detectors = track.m_detectors.size();
   if (detectors == 0) {
     // the track doesn't cross any active detector:
     // keep al material as unassigned
     if (m_plotter)
       for (unsigned int i = 1; i < track.m_steps.size(); ++i)
         m_plotter->plotSegmentUnassigned( track.m_steps[i] );
   } else {
     const double TOLERANCE = 0.0001;    // 1 um tolerance
     std::vector<double> limits(detectors + 2);
 
     // define the trivial limits
     if (m_skipBeforeFirstDetector)
       limits[0] = track.m_detectors[0].m_curvilinearIn - TOLERANCE;
     else
       limits[0] = - TOLERANCE;
     if (m_skipAfterLastDetector)
       limits[detectors] = track.m_detectors[detectors-1].m_curvilinearOut + TOLERANCE;
     else
       limits[detectors] = track.m_total.length() + TOLERANCE;
     limits[detectors+1] = INFINITY;     // this is probably no more needed, but doesn't harm...
 
     // pick the algorithm to define the non-trivial limits
     switch (m_splitMode) {
       // assign each segment to the the nearest layer
       // e.g. the material between pixel barrel 3 and TIB 1 will be split among the two
       case NEAREST_LAYER:
         for (unsigned int i = 1; i < detectors; ++i)
           limits[i] = (track.m_detectors[i-1].m_curvilinearOut + track.m_detectors[i].m_curvilinearIn) / 2.;
         break;
 
       // assign each segment to the the inner layer
       // e.g. all material between pixel barrel 3 and TIB 1 will go into the pixel barrel
       case INNER_LAYER:
         for (unsigned int i = 1; i < detectors; ++i)
           limits[i] = track.m_detectors[i].m_curvilinearIn - TOLERANCE;
         break;
 
       // assign each segment to the the outer layer
       // e.g. all material between pixel barrel 3 and TIB 1 will go into the TIB
       case OUTER_LAYER:
         for (unsigned int i = 1; i < detectors; ++i)
           limits[i] = track.m_detectors[i-1].m_curvilinearOut + TOLERANCE;
         break;
 
       case UNDEFINED:
       default:
         // throw something
         throw edm::Exception(edm::errors::LogicError) << "Invalid SplitMode";
     }
 
     //for (unsigned int i = 0; i < detectors; ++i)
     //  std::cout << "MaterialAccountingTrack::split(): detector region boundaries: [" << limits[i] << ", " << limits[i+1] << "] along track" << std::endl;
 
     double begin = 0.;          // begginning of step, along the track
     double end   = 0.;          // end of step, along the track
     unsigned int i = 1;         // step conter
 
     // skip the material before the first layer
     //std::cout << "before first layer, skipping" << std::endl;
     while (end < limits[0]) {
       const MaterialAccountingStep & step = track.m_steps[i++];
       end = begin + step.length();
 
       // do not account material before the first layer
       if (m_plotter)
         m_plotter->plotSegmentUnassigned( step );
 
       begin = end;
       //std::cout << '.';
     }
     //std::cout << std::endl;
 
     // optionally split a step across the first layer boundary
     //std::cout << "first layer (0): " << limits[0] << ".." << limits[1] << std::endl;
     if (begin < limits[0] and end > limits[0]) {
       const MaterialAccountingStep & step = track.m_steps[i++];
       end = begin + step.length();
 
       double fraction = (limits[0] - begin) / (end - begin);
       std::pair<MaterialAccountingStep, MaterialAccountingStep> parts = step.split(fraction);
 
       //std::cout << '!' << std::endl;
       track.m_detectors[0].account( parts.second, limits[1], end );
 
       if (m_plotter) {
         // step partially before first layer, keep first part as unassocated
         m_plotter->plotSegmentUnassigned( parts.first );
 
         // associate second part to first layer
         m_plotter->plotSegmentInLayer( parts.second,  group[0] );
       }
       begin = end;
     }
 
     unsigned int index = 0;     // which detector
     while (i < track.m_steps.size()) {
       const MaterialAccountingStep & step = track.m_steps[i++];
 
       end = begin + step.length();
 
       if (begin > limits[detectors]) {
         // segment after last layer and skipping requested in configuation
         if (m_plotter)
           m_plotter->plotSegmentUnassigned( step );
         begin = end;
         continue;
       }
 
       // from here onwards we should be in the accountable region, either completely in a single layer:
       //   limits[index] <= begin < end <= limits[index+1]
       // or possibly split between 2 layers
       //   limits[index] < begin < limits[index+1] < end <  limits[index+2]
       if (begin < limits[index] or end > limits[index+2]) {
         // sanity check
         std::cerr << "MaterialAccountingTrack::split(): ERROR: internal logic error, expected " << limits[index] << " < " << begin << " < " << limits[index+1] << std::endl;
         break;
       }
 
       //std::cout << '.';
       if (limits[index] <= begin and end <= limits[index+1]) {
         // step completely inside current detector range
         track.m_detectors[index].account( step, begin, end );
         if (m_plotter)
           m_plotter->plotSegmentInLayer( step, group[index] );
       } else {
         // step shared beteewn two detectors, transition at limits[index+1]
         double fraction = (limits[index+1] - begin) / (end - begin);
         std::pair<MaterialAccountingStep, MaterialAccountingStep> parts = step.split(fraction);
 
         if (m_plotter) {
           if (index > 0)
             m_plotter->plotSegmentInLayer( parts.first, group[index] );
           else
             // track outside acceptance, keep as unassocated
             m_plotter->plotSegmentUnassigned( parts.first );
 
           if (index+1 < detectors)
             m_plotter->plotSegmentInLayer( parts.second,  group[index+1] );
           else
             // track outside acceptance, keep as unassocated
             m_plotter->plotSegmentUnassigned( parts.second );
         }
 
         track.m_detectors[index].account( parts.first, begin, limits[index+1] );
         ++index;          // next layer
         //std::cout << '!' << std::endl;
         //std::cout << "next layer (" << index << "): " << limits[index] << ".." << limits[index+1] << std::endl;
         if (index < detectors)
           track.m_detectors[index].account( parts.second, limits[index+1], end );
       }
       begin = end;
     }
 
   }
   //std::cout << std::endl;
 
   // add the material from each detector to its layer (if there is one and only one)
   for (unsigned int i = 0; i < track.m_detectors.size(); ++i)
     if (group[i] != 0)
       m_groups[group[i]-1]->addDetector( track.m_detectors[i] );
 
   // end of track: commit internal buffers and reset the m_groups internal state for a new track
   for (unsigned int i = 0; i < m_groups.size(); ++i)
     m_groups[i]->endOfTrack();
 }