CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_3/src/SimTracker/TrackerMaterialAnalysis/plugins/MaterialAccountingGroup.cc

Go to the documentation of this file.
00001 #include <sstream>
00002 #include <iomanip>
00003 #include <string>
00004 #include <stdexcept>
00005 
00006 #include <TFile.h>
00007 #include <TH1F.h>
00008 #include <TProfile.h>
00009 #include <TCanvas.h>
00010 #include <TFrame.h>
00011 
00012 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
00013 #include "DataFormats/Math/interface/Vector3D.h"
00014 #include "DetectorDescription/Core/interface/DDFilteredView.h"
00015 #include "DetectorDescription/Core/interface/DDCompactView.h"
00016 #include "SimDataFormats/ValidationFormats/interface/MaterialAccountingStep.h"
00017 #include "SimDataFormats/ValidationFormats/interface/MaterialAccountingDetector.h"
00018 #include "MaterialAccountingGroup.h"
00019 
00020 double MaterialAccountingGroup::s_tolerance = 0.01; // 100um should be small enough that no elements from different layers/groups are so close
00021 
00022 MaterialAccountingGroup::MaterialAccountingGroup( const std::string & name, const DDCompactView & geometry ) :
00023   m_name( name ),
00024   m_elements(),
00025   m_boundingbox(),
00026   m_accounting(),
00027   m_errors(),
00028   m_tracks( 0 ),
00029   m_counted( false ),
00030   m_file( 0 )
00031 {
00032   // retrieve the elements from DDD
00033   DDFilteredView fv( geometry );
00034   DDSpecificsFilter filter;
00035   filter.setCriteria(DDValue("TrackingMaterialGroup", name), DDSpecificsFilter::equals);
00036   fv.addFilter(filter);
00037   while (fv.next()) {
00038     // DD3Vector and DDTranslation are the same type as math::XYZVector
00039     math::XYZVector position = fv.translation() / 10.;  // mm -> cm
00040     m_elements.push_back( GlobalPoint(position.x(), position.y(), position.z()) );
00041   }
00042 
00043   // grow the bounding box
00044   for (unsigned int i = 0; i < m_elements.size(); ++i) {
00045     m_boundingbox.grow(m_elements[i].perp(), m_elements[i].z());
00046   }
00047   m_boundingbox.grow(s_tolerance);
00048 
00049   // initialize the histograms 
00050   m_dedx_spectrum   = new TH1F((m_name + "_dedx_spectrum").c_str(),     "Energy loss spectrum",       1000,    0,   1);
00051   m_radlen_spectrum = new TH1F((m_name + "_radlen_spectrum").c_str(),   "Radiation lengths spectrum", 1000,    0,   1);
00052   m_dedx_vs_eta     = new TProfile((m_name + "_dedx_vs_eta").c_str(),   "Energy loss vs. eta",         600,   -3,   3);
00053   m_dedx_vs_z       = new TProfile((m_name + "_dedx_vs_z").c_str(),     "Energy loss vs. Z",          6000, -300, 300);
00054   m_dedx_vs_r       = new TProfile((m_name + "_dedx_vs_r").c_str(),     "Energy loss vs. R",          1200,    0, 120);
00055   m_radlen_vs_eta   = new TProfile((m_name + "_radlen_vs_eta").c_str(), "Radiation lengths vs. eta",   600,   -3,   3);
00056   m_radlen_vs_z     = new TProfile((m_name + "_radlen_vs_z").c_str(),   "Radiation lengths vs. Z",    6000, -300, 300);
00057   m_radlen_vs_r     = new TProfile((m_name + "_radlen_vs_r").c_str(),   "Radiation lengths vs. R",    1200,    0, 120);
00058   m_dedx_spectrum->SetDirectory( 0 );
00059   m_radlen_spectrum->SetDirectory( 0 );
00060   m_dedx_vs_eta->SetDirectory( 0 );
00061   m_dedx_vs_z->SetDirectory( 0 );
00062   m_dedx_vs_r->SetDirectory( 0 );
00063   m_radlen_vs_eta->SetDirectory( 0 );
00064   m_radlen_vs_z->SetDirectory( 0 );
00065   m_radlen_vs_r->SetDirectory( 0 );
00066 }
00067 
00068 MaterialAccountingGroup::~MaterialAccountingGroup(void)
00069 {
00070   delete m_dedx_spectrum;
00071   delete m_dedx_vs_eta;
00072   delete m_dedx_vs_z;
00073   delete m_dedx_vs_r;
00074   delete m_radlen_spectrum;
00075   delete m_radlen_vs_eta;
00076   delete m_radlen_vs_z;
00077   delete m_radlen_vs_r;
00078 }
00079 
00080 // TODO the inner check could be sped up in many ways
00081 // (sorting the m_elements, partitioning the bounding box, ...)
00082 // but is it worth? 
00083 // especially with the segmentation of the layers ?
00084 bool MaterialAccountingGroup::inside( const MaterialAccountingDetector& detector ) const
00085 {
00086   const GlobalPoint & position = detector.position();
00087   // first check to see if the point is inside the bounding box
00088   if (not m_boundingbox.inside(position.perp(), position.z())) {
00089     return false;
00090   } else {
00091     // now check if the point is actually close enough to any element
00092     for (unsigned int i = 0; i < m_elements.size(); ++i)
00093       if ((position - m_elements[i]).mag2() < (s_tolerance * s_tolerance))
00094         return true;
00095     return false;
00096   }
00097 }
00098 
00099 bool MaterialAccountingGroup::addDetector( const MaterialAccountingDetector& detector ) 
00100 {
00101   if (not inside(detector))
00102     return false;
00103 
00104   // multiple hits in the same layer (overlaps, etc.) from a single track still count as one for averaging,
00105   // since the energy deposits from the track have been already split between the different detectors
00106   m_buffer += detector.material();
00107   m_counted = true;
00108 
00109   return true;
00110 }
00111 
00112 void MaterialAccountingGroup::endOfTrack(void) {
00113   // add a detector
00114   if (m_counted) {
00115     m_accounting += m_buffer;
00116     m_errors     += m_buffer * m_buffer;
00117     ++m_tracks;
00118 
00119     GlobalPoint average( (m_buffer.in().x() + m_buffer.out().x()) / 2.,
00120                          (m_buffer.in().y() + m_buffer.out().y()) / 2., 
00121                          (m_buffer.in().z() + m_buffer.out().z()) / 2. );
00122     m_dedx_spectrum->Fill(   m_buffer.energyLoss() );
00123     m_radlen_spectrum->Fill( m_buffer.radiationLengths() );
00124     m_dedx_vs_eta->Fill(   average.eta(),  m_buffer.energyLoss(),       1. );
00125     m_dedx_vs_z->Fill(     average.z(),    m_buffer.energyLoss(),       1. );
00126     m_dedx_vs_r->Fill(     average.perp(), m_buffer.energyLoss(),       1. );
00127     m_radlen_vs_eta->Fill( average.eta(),  m_buffer.radiationLengths(), 1. );
00128     m_radlen_vs_z->Fill(   average.z(),    m_buffer.radiationLengths(), 1. );
00129     m_radlen_vs_r->Fill(   average.perp(), m_buffer.radiationLengths(), 1. );
00130   }
00131   m_counted = false;
00132   m_buffer  = MaterialAccountingStep();
00133 }
00134 
00135 void MaterialAccountingGroup::savePlot(TH1F * plot, const std::string & name)
00136 {
00137   TCanvas canvas(name.c_str(), plot->GetTitle(), 1280, 1024);
00138   plot->SetFillColor(15);       // grey
00139   plot->SetLineColor(1);        // black
00140   plot->Draw("c e");
00141   canvas.GetFrame()->SetFillColor(kWhite);
00142   canvas.Draw();
00143   canvas.SaveAs((name + ".png").c_str(), "");
00144 
00145   // store te plot into m_file
00146   plot->SetDirectory(m_file);
00147 }
00148 
00149 void MaterialAccountingGroup::savePlot(TProfile * plot, float average, const std::string & name)
00150 {
00151   // Nota Bene:
00152   // these "line" plots are not deleted explicitly since
00153   //   - deleting them before saving them to a TFile will not save them
00154   //   - deleting them after the TFile they're stored into results in a SEGV
00155   // ROOT is probably "taking care" (read: messing things up) somehow...
00156   TH1F * line = new TH1F((name + "_par").c_str(), "Parametrization", 1, plot->GetXaxis()->GetXmin(), plot->GetXaxis()->GetXmax());
00157   line->SetBinContent(1, average);
00158 
00159   TCanvas canvas(name.c_str(), plot->GetTitle(), 1280, 1024);
00160   plot->SetFillColor(15);       // grey
00161   plot->SetLineColor(1);        // black
00162   plot->SetLineWidth(2);
00163   plot->Draw("c e6");
00164   line->SetLineColor(2);        // red
00165   line->SetLineWidth(2);
00166   line->Draw("same");
00167   canvas.GetFrame()->SetFillColor(kWhite);
00168   canvas.Draw();
00169   canvas.SaveAs((name + ".png").c_str(), "");
00170 
00171   // store te plots into m_file
00172   plot->SetDirectory(m_file);
00173   line->SetDirectory(m_file);
00174 }
00175 
00177 std::string MaterialAccountingGroup::info(void) const
00178 {
00179   std::stringstream out;
00180   out << std::setw(48) << std::left << m_name << std::right << std::fixed;;
00181   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;
00182   out << ", "     << std::setprecision(1) << std::setw(5) << m_boundingbox.range_r().first << " < R < " << std::setprecision(1) << std::setw(5) << m_boundingbox.range_r().second;
00183   out << "   Elements: " << std::setw(6) << m_elements.size();
00184   return out.str();
00185 }
00186 
00187 
00188 void MaterialAccountingGroup::savePlots(void)
00189 {
00190   m_file = new TFile((m_name + ".root").c_str(), "RECREATE");
00191   savePlot(m_dedx_spectrum,   m_name + "_dedx_spectrum");
00192   savePlot(m_radlen_spectrum, m_name + "_radlen_spectrum");
00193   savePlot(m_dedx_vs_eta,     averageEnergyLoss(),       m_name + "_dedx_vs_eta");
00194   savePlot(m_dedx_vs_z,       averageEnergyLoss(),       m_name + "_dedx_vs_z");
00195   savePlot(m_dedx_vs_r,       averageEnergyLoss(),       m_name + "_dedx_vs_r");
00196   savePlot(m_radlen_vs_eta,   averageRadiationLengths(), m_name + "_radlen_vs_eta");
00197   savePlot(m_radlen_vs_z,     averageRadiationLengths(), m_name + "_radlen_vs_z");
00198   savePlot(m_radlen_vs_r,     averageRadiationLengths(), m_name + "_radlen_vs_r");
00199   m_file->Write();
00200   m_file->Close();
00201 
00202   delete m_file;
00203 }