TrackerMaterialAnalysis/plugins/MaterialAccountingGroup.h

0001 #ifndef MaterialAccountingGroup_h
0002 #define MaterialAccountingGroup_h
0003
0004 #include <string>
0005 #include <stdexcept>
0006 #include <utility>
0007
0008 #include "Geometry/CommonDetUnit/interface/GeomDetEnumerators.h"
0009
0010 #include "SimDataFormats/ValidationFormats/interface/MaterialAccountingStep.h"
0011 #include "SimDataFormats/ValidationFormats/interface/MaterialAccountingDetector.h"
0012
0013 class TH1;
0014 class TH1F;
0015 class TProfile;
0016 class TFile;
0017 class DDCompactView;
0018
0019 class MaterialAccountingGroup {
0020 private:
0021   // quick implementation of a bounding cilinder
0022   class BoundingBox {
0023   public:
0024     BoundingBox(double min_r, double max_r, double min_z, double max_z)
0025         : r_min(min_r), r_max(max_r), z_min(min_z), z_max(max_z) {}
0026
0027     BoundingBox() : r_min(0.), r_max(0.), z_min(0.), z_max(0.) {}
0028
0029     void grow(double r, double z) {
0030       if (r < r_min)
0031         r_min = r;
0032       if (r > r_max)
0033         r_max = r;
0034       if (z < z_min)
0035         z_min = z;
0036       if (z > z_max)
0037         z_max = z;
0038     }
0039
0040     void grow(double skin) {
0041       r_min -= skin;  // yes, we allow r_min to go negative
0042       r_max += skin;
0043       z_min -= skin;
0044       z_max += skin;
0045     }
0046
0047     bool inside(double r, double z) const { return (r >= r_min and r <= r_max and z >= z_min and z <= z_max); }
0048
0049     std::pair<double, double> range_r() const { return std::make_pair(r_min, r_max); }
0050
0051     std::pair<double, double> range_z() const { return std::make_pair(z_min, z_max); }
0052
0053   private:
0054     double r_min;
0055     double r_max;
0056     double z_min;
0057     double z_max;
0058   };
0059
0060 public:
0061   /// explicit constructors
0062   MaterialAccountingGroup(const std::string& name, const DDCompactView& geometry);
0063
0064   /// stop default copy ctor
0065   MaterialAccountingGroup(const MaterialAccountingGroup& layer) = delete;
0066
0067   /// stop default assignment operator
0068   MaterialAccountingGroup& operator=(const MaterialAccountingGroup& layer) = delete;
0069
0070   /// destructor
0071   ~MaterialAccountingGroup(void);
0072
0073 public:
0074   /// buffer material from a detector, if the detector is inside the DetLayer bounds
0075   bool addDetector(const MaterialAccountingDetector& detector);
0076
0077   /// commit the buffer and reset the "already hit by this track" flag
0078   void endOfTrack(void);
0079
0080   /// check if detector is inside any part of this layer
0081   bool inside(const MaterialAccountingDetector& detector) const;
0082
0083   /// Return the bouding limit in R for the hosted Group
0084   std::pair<double, double> getBoundingR() const { return m_boundingbox.range_r(); };
0085
0086   /// Return the bouding limit in Z for the hosted Group
0087   std::pair<double, double> getBoundingZ() const { return m_boundingbox.range_z(); };
0088
0089   /// return the average normalized material accounting informations
0090   MaterialAccountingStep average(void) const { return m_tracks ? m_accounting / m_tracks : MaterialAccountingStep(); }
0091
0092   /// return the average normalized layer thickness
0093   double averageLength(void) const { return m_tracks ? m_accounting.length() / m_tracks : 0.; }
0094
0095   /// return the average normalized number of radiation lengths
0096   double averageRadiationLengths(void) const { return m_tracks ? m_accounting.radiationLengths() / m_tracks : 0.; }
0097
0098   /// return the average normalized energy loss density factor for Bethe-Bloch
0099   double averageEnergyLoss(void) const { return m_tracks ? m_accounting.energyLoss() / m_tracks : 0.; }
0100
0101   /// return the sigma of the normalized layer thickness
0102   double sigmaLength(void) const {
0103     return m_tracks ? std::sqrt(m_errors.length() / m_tracks - averageLength() * averageLength()) : 0.;
0104   }
0105
0106   /// return the sigma of the normalized number of radiation lengths
0107   double sigmaRadiationLengths(void) const {
0108     return m_tracks ? std::sqrt(m_errors.radiationLengths() / m_tracks -
0109                                 averageRadiationLengths() * averageRadiationLengths())
0110                     : 0.;
0111   }
0112
0113   /// return the sigma of the normalized energy loss density factor for Bethe-Bloch
0114   double sigmaEnergyLoss(void) const {
0115     return m_tracks ? std::sqrt(m_errors.energyLoss() / m_tracks - averageEnergyLoss() * averageEnergyLoss()) : 0.;
0116   }
0117
0118   /// return the number of tracks that hit this layer
0119   unsigned int tracks(void) const { return m_tracks; }
0120
0121   /// get the layer name
0122   const std::string& name(void) const { return m_name; }
0123
0124   /// get some infos
0125   std::string info(void) const;
0126
0127   /// save the plots
0128   void savePlots(void);
0129
0130   // get directly the internal vector of global points that
0131   // correspond, roughly, to the center of each subcomponent of this
0132   // group.
0133
0134   const std::vector<GlobalPoint>& elements(void) const { return m_elements; }
0135
0136 private:
0137   void savePlot(TH1F* plot, const std::string& name);
0138   void savePlot(TProfile* plot, float average, const std::string& name);
0139
0140   std::string m_name;
0141   std::vector<GlobalPoint> m_elements;
0142   BoundingBox m_boundingbox;
0143   MaterialAccountingStep m_accounting;
0144   MaterialAccountingStep m_errors;
0145   unsigned int m_tracks;
0146   bool m_counted;
0147   MaterialAccountingStep m_buffer;
0148
0149   // plots of material effects distribution
0150   TH1F* m_dedx_spectrum;
0151   TH1F* m_radlen_spectrum;
0152   // plots of material effects vs. η / Z / R
0153   TProfile* m_dedx_vs_eta;
0154   TProfile* m_dedx_vs_z;
0155   TProfile* m_dedx_vs_r;
0156   TProfile* m_radlen_vs_eta;
0157   TProfile* m_radlen_vs_z;
0158   TProfile* m_radlen_vs_r;
0159
0160   // file to store plots into
0161   mutable TFile* m_file;
0162
0163   // 100um should be small enough that no elements from different layers/groups are so close
0164   static double const s_tolerance;
0165 };
0166
0167 #endif  // MaterialAccountingGroup_h