CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
MaterialAccountingGroup.cc
Go to the documentation of this file.
1 #include <sstream>
2 #include <iomanip>
3 #include <string>
4 #include <stdexcept>
5 
6 #include <TFile.h>
7 #include <TH1F.h>
8 #include <TProfile.h>
9 #include <TCanvas.h>
10 #include <TFrame.h>
11 
19 
20 double const MaterialAccountingGroup::s_tolerance = 0.01; // 100um should be small enough that no elements from different layers/groups are so close
21 
23  m_name( name ),
24  m_elements(),
25  m_boundingbox(),
26  m_accounting(),
27  m_errors(),
28  m_tracks( 0 ),
29  m_counted( false ),
30  m_file( 0 )
31 {
32  // retrieve the elements from DDD
33  DDFilteredView fv( geometry );
35  filter.setCriteria(DDValue("TrackingMaterialGroup", name), DDSpecificsFilter::equals);
36  fv.addFilter(filter);
37  while (fv.next()) {
38  // DD3Vector and DDTranslation are the same type as math::XYZVector
39  math::XYZVector position = fv.translation() / 10.; // mm -> cm
40  m_elements.push_back( GlobalPoint(position.x(), position.y(), position.z()) );
41  }
42 
43  // grow the bounding box
44  for (unsigned int i = 0; i < m_elements.size(); ++i) {
46  }
48 
49  // initialize the histograms
50  m_dedx_spectrum = new TH1F((m_name + "_dedx_spectrum").c_str(), "Energy loss spectrum", 1000, 0, 1);
51  m_radlen_spectrum = new TH1F((m_name + "_radlen_spectrum").c_str(), "Radiation lengths spectrum", 1000, 0, 1);
52  m_dedx_vs_eta = new TProfile((m_name + "_dedx_vs_eta").c_str(), "Energy loss vs. eta", 600, -3, 3);
53  m_dedx_vs_z = new TProfile((m_name + "_dedx_vs_z").c_str(), "Energy loss vs. Z", 6000, -300, 300);
54  m_dedx_vs_r = new TProfile((m_name + "_dedx_vs_r").c_str(), "Energy loss vs. R", 1200, 0, 120);
55  m_radlen_vs_eta = new TProfile((m_name + "_radlen_vs_eta").c_str(), "Radiation lengths vs. eta", 600, -3, 3);
56  m_radlen_vs_z = new TProfile((m_name + "_radlen_vs_z").c_str(), "Radiation lengths vs. Z", 6000, -300, 300);
57  m_radlen_vs_r = new TProfile((m_name + "_radlen_vs_r").c_str(), "Radiation lengths vs. R", 1200, 0, 120);
58  m_dedx_spectrum->SetDirectory( 0 );
59  m_radlen_spectrum->SetDirectory( 0 );
60  m_dedx_vs_eta->SetDirectory( 0 );
61  m_dedx_vs_z->SetDirectory( 0 );
62  m_dedx_vs_r->SetDirectory( 0 );
63  m_radlen_vs_eta->SetDirectory( 0 );
64  m_radlen_vs_z->SetDirectory( 0 );
65  m_radlen_vs_r->SetDirectory( 0 );
66 }
67 
69 {
70  delete m_dedx_spectrum;
71  delete m_dedx_vs_eta;
72  delete m_dedx_vs_z;
73  delete m_dedx_vs_r;
74  delete m_radlen_spectrum;
75  delete m_radlen_vs_eta;
76  delete m_radlen_vs_z;
77  delete m_radlen_vs_r;
78 }
79 
80 // TODO the inner check could be sped up in many ways
81 // (sorting the m_elements, partitioning the bounding box, ...)
82 // but is it worth?
83 // especially with the segmentation of the layers ?
85 {
86  const GlobalPoint & position = detector.position();
87  // first check to see if the point is inside the bounding box
88  if (not m_boundingbox.inside(position.perp(), position.z())) {
89  return false;
90  } else {
91  // now check if the point is actually close enough to any element
92  for (unsigned int i = 0; i < m_elements.size(); ++i)
93  if ((position - m_elements[i]).mag2() < (s_tolerance * s_tolerance))
94  return true;
95  return false;
96  }
97 }
98 
100 {
101  if (not inside(detector))
102  return false;
103 
104  // multiple hits in the same layer (overlaps, etc.) from a single track still count as one for averaging,
105  // since the energy deposits from the track have been already split between the different detectors
106  m_buffer += detector.material();
107  m_counted = true;
108 
109  return true;
110 }
111 
113  // add a detector
114  if (m_counted) {
117  ++m_tracks;
118 
119  GlobalPoint average( (m_buffer.in().x() + m_buffer.out().x()) / 2.,
120  (m_buffer.in().y() + m_buffer.out().y()) / 2.,
121  (m_buffer.in().z() + m_buffer.out().z()) / 2. );
122  m_dedx_spectrum->Fill( m_buffer.energyLoss() );
123  m_radlen_spectrum->Fill( m_buffer.radiationLengths() );
124  m_dedx_vs_eta->Fill( average.eta(), m_buffer.energyLoss(), 1. );
125  m_dedx_vs_z->Fill( average.z(), m_buffer.energyLoss(), 1. );
126  m_dedx_vs_r->Fill( average.perp(), m_buffer.energyLoss(), 1. );
127  m_radlen_vs_eta->Fill( average.eta(), m_buffer.radiationLengths(), 1. );
128  m_radlen_vs_z->Fill( average.z(), m_buffer.radiationLengths(), 1. );
129  m_radlen_vs_r->Fill( average.perp(), m_buffer.radiationLengths(), 1. );
130  }
131  m_counted = false;
133 }
134 
136 {
137  TCanvas canvas(name.c_str(), plot->GetTitle(), 1280, 1024);
138  plot->SetFillColor(15); // grey
139  plot->SetLineColor(1); // black
140  plot->Draw("c e");
141  canvas.GetFrame()->SetFillColor(kWhite);
142  canvas.Draw();
143  canvas.SaveAs((name + ".png").c_str(), "");
144 
145  // store te plot into m_file
146  plot->SetDirectory(m_file);
147 }
148 
150 {
151  // Nota Bene:
152  // these "line" plots are not deleted explicitly since
153  // - deleting them before saving them to a TFile will not save them
154  // - deleting them after the TFile they're stored into results in a SEGV
155  // ROOT is probably "taking care" (read: messing things up) somehow...
156  TH1F * line = new TH1F((name + "_par").c_str(), "Parametrization", 1, plot->GetXaxis()->GetXmin(), plot->GetXaxis()->GetXmax());
157  line->SetBinContent(1, average);
158 
159  TCanvas canvas(name.c_str(), plot->GetTitle(), 1280, 1024);
160  plot->SetFillColor(15); // grey
161  plot->SetLineColor(1); // black
162  plot->SetLineWidth(2);
163  plot->Draw("c e6");
164  line->SetLineColor(2); // red
165  line->SetLineWidth(2);
166  line->Draw("same");
167  canvas.GetFrame()->SetFillColor(kWhite);
168  canvas.Draw();
169  canvas.SaveAs((name + ".png").c_str(), "");
170 
171  // store te plots into m_file
172  plot->SetDirectory(m_file);
173  line->SetDirectory(m_file);
174 }
175 
178 {
179  std::stringstream out;
180  out << std::setw(48) << std::left << m_name << std::right << std::fixed;;
181  out << "BBox: " << std::setprecision(1) << std::setw(6) << m_boundingbox.range_z().first << " < Z < " << std::setprecision(1) << std::setw(6) << m_boundingbox.range_z().second;
182  out << ", " << std::setprecision(1) << std::setw(5) << m_boundingbox.range_r().first << " < R < " << std::setprecision(1) << std::setw(5) << m_boundingbox.range_r().second;
183  out << " Elements: " << std::setw(6) << m_elements.size();
184  return out.str();
185 }
186 
187 
189 {
190  m_file = new TFile((m_name + ".root").c_str(), "RECREATE");
191  savePlot(m_dedx_spectrum, m_name + "_dedx_spectrum");
192  savePlot(m_radlen_spectrum, m_name + "_radlen_spectrum");
193  savePlot(m_dedx_vs_eta, averageEnergyLoss(), m_name + "_dedx_vs_eta");
194  savePlot(m_dedx_vs_z, averageEnergyLoss(), m_name + "_dedx_vs_z");
195  savePlot(m_dedx_vs_r, averageEnergyLoss(), m_name + "_dedx_vs_r");
196  savePlot(m_radlen_vs_eta, averageRadiationLengths(), m_name + "_radlen_vs_eta");
199  m_file->Write();
200  m_file->Close();
201 
202  delete m_file;
203 }
int i
Definition: DBlmapReader.cc:9
void addFilter(const DDFilter &, log_op op=AND)
MaterialAccountingGroup(const std::string &name, const DDCompactView &geometry)
explicit constructors
T perp() const
Definition: PV3DBase.h:72
double averageRadiationLengths(void) const
return the average normalized number of radiation lengths
void endOfTrack(void)
commit the buffer and reset the &quot;already hit by this track&quot; flag
MaterialAccountingStep m_errors
MaterialAccountingStep m_buffer
std::vector< GlobalPoint > m_elements
Global3DPoint GlobalPoint
Definition: GlobalPoint.h:10
std::pair< double, double > range_r() const
std::string info(void) const
get some infos
type of data representation of DDCompactView
Definition: DDCompactView.h:77
float float float z
def canvas
Definition: svgfig.py:481
double averageEnergyLoss(void) const
return the average normalized energy loss density factor for Bethe-Bloch
bool next()
set current node to the next node in the filtered tree
void savePlot(TH1F *plot, const std::string &name)
T z() const
Definition: PV3DBase.h:64
MaterialAccountingStep average(void) const
return the average normalized material accounting informations
std::pair< double, double > range_z() const
const MaterialAccountingStep & material() const
tuple out
Definition: dbtoconf.py:99
XYZVectorD XYZVector
spatial vector with cartesian internal representation
Definition: Vector3D.h:30
void savePlots(void)
save the plots
const GlobalPoint & position() const
int average
Definition: PDRates.py:137
T eta() const
Definition: PV3DBase.h:76
ESHandle< TrackerGeometry > geometry
T perp() const
Magnitude of transverse component.
bool inside(double r, double z) const
static int position[264][3]
Definition: ReadPGInfo.cc:509
bool addDetector(const MaterialAccountingDetector &detector)
buffer material from a detector, if the detector is inside the DetLayer bounds
~MaterialAccountingGroup(void)
destructor
void setCriteria(const DDValue &nameVal, comp_op, log_op l=AND, bool asString=true, bool merged=true)
Definition: DDFilter.cc:285
const DDTranslation & translation() const
The absolute translation of the current node.
MaterialAccountingStep m_accounting
volatile std::atomic< bool > shutdown_flag false
bool inside(const MaterialAccountingDetector &detector) const
check if detector is inside any part of this layer
The DDGenericFilter is a runtime-parametrized Filter looking on DDSpecifcs.
Definition: DDFilter.h:37