#include <TrackingMaterialAnalyser.h>

Inheritance diagram for TrackingMaterialAnalyser:

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

virtual	~TrackingMaterialAnalyser ()

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

	EDAnalyzer ()

ModuleDescription const &	moduleDescription () const

std::string	workerType () const

virtual	~EDAnalyzer ()

Public Member Functions inherited from edm::EDConsumerBase
std::vector< ConsumesInfo >	consumesInfo () const

	EDConsumerBase ()

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

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

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

std::vector < ProductHolderIndexAndSkipBit > const &	itemsToGetFromEvent () const

void	labelsForToken (EDGetToken iToken, Labels &oLabels) const

void	modulesDependentUpon (const std::string &iProcessName, std::vector< const char * > &oModuleLabels) const

void	modulesWhoseProductsAreConsumed (std::vector< ModuleDescription const * > &modules, ProductRegistry const &preg, std::map< std::string, ModuleDescription const * > const &labelsToDesc, std::string const &processName) const

bool	registeredToConsume (ProductHolderIndex, bool, BranchType) const

bool	registeredToConsumeMany (TypeID const &, BranchType) 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 ()

void	endJob ()

int	findLayer (const MaterialAccountingDetector &detector)

void	saveLayerPlots ()

void	saveParameters (const char *name)

void	saveXml (const char *name)

void	split (MaterialAccountingTrack &track)

Private Attributes
std::vector< std::string >	m_groupNames

std::vector < MaterialAccountingGroup * >	m_groups

edm::EDGetTokenT< std::vector < MaterialAccountingTrack > >	m_materialToken

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

Public Types inherited from edm::EDConsumerBase
typedef ProductLabels	Labels

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::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 16 of file TrackingMaterialAnalyser.h.

Member Enumeration Documentation

enum TrackingMaterialAnalyser::SplitMode

private

Enumerator
NEAREST_LAYER
INNER_LAYER
OUTER_LAYER
UNDEFINED

Definition at line 23 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 Exception, edm::ParameterSet::getParameter(), INNER_LAYER, edm::errors::LogicError, m_groupNames, m_materialToken, 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_materialToken           = consumes<std::vector<MaterialAccountingTrack> >(
       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 59 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 137 of file TrackingMaterialAnalyser.cc.

References end, edm::EventSetup::get(), i, LogDebug, m_groupNames, m_groups, m_materialToken, split(), and tree::t.

 {
   edm::ESTransientHandle<DDCompactView> hDDD;
   setup.get<IdealGeometryRecord>().get( hDDD );
 
   m_groups.reserve( m_groupNames.size() );
   // Initialize m_groups iff it has size equal to zero, so that we are
   // sure it will never be repopulated with the same entries over and
   // over again in the eventloop, at each call of the analyze method.
   if (m_groups.size() == 0) {
     for (unsigned int i = 0; i < m_groupNames.size(); ++i)
       m_groups.push_back( new MaterialAccountingGroup( m_groupNames[i], * hDDD) ); 
 
     LogDebug("TrackingMaterialAnalyser")
         << "TrackingMaterialAnalyser: List of the tracker groups: " << std::endl;
     for (unsigned int i = 0; i < m_groups.size(); ++i)
       LogDebug("TrackingMaterialAnalyser")
           << "TrackingMaterialAnalyser:\t" << m_groups[i]->info() << std::endl;
   }
   edm::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) {
     MaterialAccountingTrack track(*t);
     split( track );
   }
 }

void TrackingMaterialAnalyser::beginJob ( void )

inlineprivatevirtual

Reimplemented from edm::EDAnalyzer.

Definition at line 31 of file TrackingMaterialAnalyser.h.

31 {}

void TrackingMaterialAnalyser::endJob ( void )

privatevirtual

Reimplemented from edm::EDAnalyzer.

Definition at line 117 of file TrackingMaterialAnalyser.cc.

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

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

int TrackingMaterialAnalyser::findLayer ( const MaterialAccountingDetector & detector )

private

Definition at line 330 of file TrackingMaterialAnalyser.cc.

References ecal_dqm_sourceclient-live_cfg::cerr, i, cmsHarvester::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;
   }
 
   // if (index > 0)
   //   std::cout << m_groups[index-1]->info() << " " << index << std::endl;
   return index;
 }

void TrackingMaterialAnalyser::saveLayerPlots ( )

private

Definition at line 108 of file TrackingMaterialAnalyser.cc.

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

Referenced by endJob().

 {
   for (unsigned int i = 0; i < m_groups.size(); ++i) {
     MaterialAccountingGroup & layer = *(m_groups[i]);
     layer.savePlots();
   }
 }

void TrackingMaterialAnalyser::saveParameters ( const char * name )

private

Definition at line 66 of file TrackingMaterialAnalyser.cc.

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

Referenced by endJob().

 {
   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:       %6.5fe-03 ± %6.5fe-03 GeV") % 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 91 of file TrackingMaterialAnalyser.cc.

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

Referenced by endJob().

 {
   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 174 of file TrackingMaterialAnalyser.cc.

References begin, ecal_dqm_sourceclient-live_cfg::cerr, MuonGeometrySanityCheck_cfi::detectors(), MaterialAccountingTrack::detectors(), end, Exception, findLayer(), HLT_25ns14e33_v1_cff::fraction, watchdog::group, i, cmsHarvester::index, INNER_LAYER, MaterialAccountingStep::length(), edm::errors::LogicError, m_groups, m_plotter, m_skipAfterLastDetector, m_skipBeforeFirstDetector, m_splitMode, NEAREST_LAYER, or, OUTER_LAYER, contentValuesFiles::parts, TrackingMaterialPlotter::plotSegmentInLayer(), TrackingMaterialPlotter::plotSegmentUnassigned(), MaterialAccountingStep::split(), relval_parameters_module::step, MaterialAccountingTrack::steps(), MaterialAccountingTrack::summary(), and UNDEFINED.

Referenced by analyze().

 {
   // group sensitive detectors by their DetLayer
   std::vector<int> group( track.detectors().size() );
   for (unsigned int i = 0; i < track.detectors().size(); ++i)
     group[i] = findLayer( track.detectors()[i] );
 
   unsigned int detectors = track.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.steps().size(); ++i)
         m_plotter->plotSegmentUnassigned( track.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.detectors()[0].m_curvilinearIn - TOLERANCE;
     else
       limits[0] = - TOLERANCE;
     if (m_skipAfterLastDetector)
       limits[detectors] = track.detectors()[detectors-1].m_curvilinearOut + TOLERANCE;
     else
       limits[detectors] = track.summary().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.detectors()[i-1].m_curvilinearOut + track.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.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.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.;          // beginning 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.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;
 
     unsigned int index = 0;     // which detector
     while (i < track.steps().size()) {
       const MaterialAccountingStep & step = track.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.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.detectors()[index].account( parts.first, begin, limits[index+1] );
         ++index;          // next layer
         //std::cout << '!' << std::endl;
         // std::cout << "next layer (" << index << "): "
         //           << " old det: " << group[index-1] << " new det: " << group[index]
         //           << " " << limits[index] << ".." << limits[index+1] << std::endl;
         if (index < detectors)
           track.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.detectors().size(); ++i)
     if (group[i] != 0)
       m_groups[group[i]-1]->addDetector( track.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();
 }