Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:31:04 

0001 #include <iomanip>
0002 #include <iostream>
0003 #include <map>
0004 #include <string>
0005 #include <vector>
0006 
0007 #include <fmt/printf.h>
0008 
0009 // ROOT
0010 #include <TCanvas.h>
0011 #include <TColor.h>
0012 #include <TFrame.h>
0013 #include <TLegend.h>
0014 #include <TLegendEntry.h>
0015 #include <TLine.h>
0016 #include <TProfile2D.h>
0017 #include <TROOT.h>
0018 #include <TStyle.h>
0019 #include <TText.h>
0020 
0021 #include "DataFormats/DetId/interface/DetId.h"
0022 #include "DataFormats/Math/interface/Vector3D.h"
0023 #include "DetectorDescription/Core/interface/DDCompactView.h"
0024 #include "DetectorDescription/Core/interface/DDFilteredView.h"
0025 #include "DetectorDescription/Core/interface/DDMaterial.h"
0026 #include "FWCore/Framework/interface/ESTransientHandle.h"
0027 #include "FWCore/Framework/interface/Event.h"
0028 #include "FWCore/Framework/interface/EventSetup.h"
0029 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0030 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0031 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0032 #include "FWCore/ParameterSet/interface/types.h"
0033 #include "FWCore/Utilities/interface/InputTag.h"
0034 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0035 
0036 #include "MaterialAccountingGroup.h"
0037 #include "TrackingMaterialPlotter.h"
0038 
0039 static bool dddGetStringRaw(const DDFilteredView &view, const std::string &name, std::string &value) {
0040   std::vector<const DDsvalues_type *> result;
0041   view.specificsV(result);
0042   for (std::vector<const DDsvalues_type *>::iterator it = result.begin(); it != result.end(); ++it) {
0043     DDValue parameter(name);
0044     if (DDfetch(*it, parameter)) {
0045       if (parameter.strings().size() == 1) {
0046         value = parameter.strings().front();
0047         return true;
0048       } else {
0049         throw cms::Exception("Configuration") << " ERROR: multiple " << name << " tags encountered";
0050         return false;
0051       }
0052     }
0053   }
0054   return false;
0055 }
0056 
0057 /*
0058 static inline
0059 double dddGetDouble(const std::string & s, DDFilteredView const & view) {
0060   std::string value;
0061   if (dddGetStringRaw(view, s, value))
0062     return double(::atof(value.c_str()));
0063   else
0064     return NAN;
0065 }
0066 */
0067 
0068 static inline std::string dddGetString(const std::string &s, DDFilteredView const &view) {
0069   std::string value;
0070   if (dddGetStringRaw(view, s, value))
0071     return value;
0072   else
0073     return std::string();
0074 }
0075 
0076 class ListGroups : public edm::one::EDAnalyzer<> {
0077 public:
0078   ListGroups(const edm::ParameterSet &);
0079   ~ListGroups() override;
0080 
0081 private:
0082   void analyze(const edm::Event &, const edm::EventSetup &) override;
0083   void beginJob() override {}
0084   void endJob() override;
0085   void fillColor();
0086   void fillMaterialDifferences();
0087   void fillGradient();
0088   std::vector<std::pair<std::shared_ptr<TLine>, std::shared_ptr<TText> > > overlayEtaReferences();
0089   void produceAndSaveSummaryPlot(const edm::EventSetup &);
0090 
0091   const bool m_saveSummaryPlot;
0092   const edm::ESGetToken<DDCompactView, IdealGeometryRecord> ddToken_;
0093   std::vector<TH2F *> m_plots;
0094   std::set<std::string> m_group_names;
0095   std::vector<MaterialAccountingGroup *> m_groups;
0096   std::vector<unsigned int> m_color;
0097   std::vector<int> m_gradient;
0098 
0099   // The following maps are automatically filled by the script
0100   // dumpFullXML, when run with -c,--compare flag, and are injected in
0101   // this code via the header ListGroupsMaterialDifference.h, which is
0102   // included below. The first value of the pair in m_diff represents
0103   // the relative difference ((new - old)/old * 100, in %) of
0104   // radiation length, while the second referes to the energy loss (in
0105   // GeV/cm) changes. The values in m_values are ordered in the very
0106   // same way, ie. they contain the new values for radiation length
0107   // and energy loss, respectively.
0108   std::map<std::string, std::pair<float, float> > m_diff;
0109   std::map<std::string, std::pair<float, float> > m_values;
0110 };
0111 
0112 ListGroups::ListGroups(const edm::ParameterSet &iPSet)
0113     : m_saveSummaryPlot(iPSet.getUntrackedParameter<bool>("SaveSummaryPlot")),
0114       ddToken_(esConsumes<DDCompactView, IdealGeometryRecord>()) {
0115   m_plots.clear();
0116   m_groups.clear();
0117   TColor::InitializeColors();
0118   fillColor();
0119   fillMaterialDifferences();
0120   fillGradient();
0121 }
0122 
0123 ListGroups::~ListGroups() {
0124   for (auto plot : m_plots)
0125     delete plot;
0126 
0127   if (!m_groups.empty())
0128     for (auto g : m_groups)
0129       delete g;
0130 }
0131 
0132 #include "ListGroupsMaterialDifference.h"
0133 
0134 void ListGroups::fillGradient(void) {
0135   m_gradient.reserve(200);
0136   unsigned int steps = 100;
0137   // if no index was given, find the highest used one and start from that plus one
0138   unsigned int index = ((TObjArray *)gROOT->GetListOfColors())->GetLast() + 1;
0139 
0140   float r1, g1, b1, r2, g2, b2;
0141   static_cast<TColor *>(gROOT->GetListOfColors()->At(kBlue + 1))->GetRGB(r1, g1, b1);
0142   static_cast<TColor *>(gROOT->GetListOfColors()->At(kAzure + 10))->GetRGB(r2, g2, b2);
0143   float delta_r = (r2 - r1) / (steps - 1);
0144   float delta_g = (g2 - g1) / (steps - 1);
0145   float delta_b = (b2 - b1) / (steps - 1);
0146 
0147   m_gradient.push_back(kBlue + 4);  // Underflow lowest bin
0148   unsigned int ii = 0;
0149   for (unsigned int i = 0; i < steps; ++i, ++ii) {
0150     new TColor(static_cast<Int_t>(index + ii), r1 + delta_r * i, g1 + delta_g * i, b1 + delta_b * i);
0151     m_gradient.push_back(index + ii);
0152   }
0153 
0154   m_gradient.push_back(kWhite);  // 0 level perfectly white
0155 
0156   static_cast<TColor *>(gROOT->GetListOfColors()->At(kOrange))->GetRGB(r1, g1, b1);
0157   static_cast<TColor *>(gROOT->GetListOfColors()->At(kOrange + 7))->GetRGB(r2, g2, b2);
0158   delta_r = (r2 - r1) / (steps - 1);
0159   delta_g = (g2 - g1) / (steps - 1);
0160   delta_b = (b2 - b1) / (steps - 1);
0161   for (unsigned int i = 0; i < steps; ++i, ++ii) {
0162     new TColor(static_cast<Int_t>(index + ii), r1 + delta_r * i, g1 + delta_g * i, b1 + delta_b * i);
0163     m_gradient.push_back(index + ii);
0164   }
0165   m_gradient.push_back(kRed);  // Overflow highest bin
0166 }
0167 
0168 void ListGroups::fillColor(void) {
0169   // With the introduction of the support for PhaseI and PhaseII detectors it
0170   // became quite difficult to maintain a list of colors that is in sync with
0171   // the real number of grouping used in the different scenarios. We therefore
0172   // define some reasonable set and loop over it in case the number of grouping
0173   // is larger than the number of colors.
0174 
0175   m_color.push_back(kBlack);  // unassigned
0176 
0177   m_color.push_back(kAzure);      // PixelBarrelLayer0_Z0
0178   m_color.push_back(kAzure - 1);  // PixelBarrelLayer0_Z20
0179   m_color.push_back(kAzure + 1);  // Layer1_Z0
0180   m_color.push_back(kAzure + 2);  // Layer1_Z20
0181 
0182   m_color.push_back(kGreen);       // EndCapDisk1_R0
0183   m_color.push_back(kGreen + 2);   // EndcapDisk1_R11
0184   m_color.push_back(kGreen + 4);   // EndcapDisk1_R7
0185   m_color.push_back(kSpring + 9);  // EndcapDisk2_R0
0186   m_color.push_back(kSpring + 4);  // EndcapDisk2_R7
0187   m_color.push_back(kSpring);      // EndcapDisk2_R7
0188 
0189   m_color.push_back(kRed);          // TECDisk0_R20
0190   m_color.push_back(kRed + 2);      // TECDisk0_R40
0191   m_color.push_back(kRed - 7);      // TECDisk0_R50
0192   m_color.push_back(kRed - 5);      // TECDisk0_R60
0193   m_color.push_back(kRed - 10);     // TECDisk0_R90
0194   m_color.push_back(kRed - 1);      // TECDisk1_Inner
0195   m_color.push_back(kRed - 2);      // TECDisk1_Outer
0196   m_color.push_back(kRed - 3);      // TECDisk1_R20
0197   m_color.push_back(kPink - 2);     // TECDisk2_Inner
0198   m_color.push_back(kPink - 3);     // TECDisk2_Outer
0199   m_color.push_back(kPink - 4);     // TECDisk2_R20
0200   m_color.push_back(kPink + 9);     // TECDisk3_Inner
0201   m_color.push_back(kPink + 8);     // TECDisk3_Outer
0202   m_color.push_back(kPink + 7);     // TECDisk3
0203   m_color.push_back(kMagenta - 2);  // TECDisk4_Inner
0204   m_color.push_back(kMagenta - 3);  // TECDisk4_Outer
0205   m_color.push_back(kMagenta - 4);  // TECDisk4_R33
0206   m_color.push_back(kMagenta - 5);  // TECDisk5_Inner
0207   m_color.push_back(kMagenta - 6);  // TECDisk5_Outer
0208   m_color.push_back(kMagenta - 7);  // TECDisk5_R33
0209   m_color.push_back(kRed);          // TECDisk6
0210   m_color.push_back(kMagenta - 9);  // TECDisk7_R40
0211   m_color.push_back(kViolet);       // TECDisk8
0212 
0213   m_color.push_back(kOrange + 9);   // TIBLayer0_Z0
0214   m_color.push_back(kOrange + 7);   // TIBLayer0_Z20
0215   m_color.push_back(kOrange + 5);   // TIBLayer0_Z40
0216   m_color.push_back(kOrange - 2);   // TIBLayer1_Z0
0217   m_color.push_back(kOrange - 3);   // TIBLayer1_Z30
0218   m_color.push_back(kOrange - 6);   // TIBLayer1_Z60
0219   m_color.push_back(kOrange + 4);   // TIBLayer2_Z0
0220   m_color.push_back(kOrange - 7);   // TIBLayer2_Z40
0221   m_color.push_back(kOrange);       // TIBLayer3_Z0
0222   m_color.push_back(kOrange + 10);  // TIBLayer3_Z50
0223 
0224   m_color.push_back(kViolet + 10);  // TIDDisk1_R0
0225   m_color.push_back(kViolet + 6);   // TIDDisk1_R30
0226   m_color.push_back(kViolet + 3);   // TIDDisk1_R40
0227   m_color.push_back(kViolet - 7);   // TIDDisk2_R25
0228   m_color.push_back(kViolet - 1);   // TIDDisk2_R30
0229   m_color.push_back(kViolet + 9);   // TIDDisk2_R40
0230   m_color.push_back(kViolet - 5);   // TIDDisk3_R24
0231   m_color.push_back(kViolet - 3);   // TIDDisk3_R30
0232   m_color.push_back(kViolet);       // TIDDisk3_R40
0233 
0234   m_color.push_back(kAzure);      // TOBLayer0_Z0
0235   m_color.push_back(kAzure + 8);  // TOBLayer0_Z20
0236   m_color.push_back(kAzure + 2);  // TOBLayer0_Z70
0237   m_color.push_back(kAzure + 4);  // TOBLayer0_Z80
0238   m_color.push_back(kCyan + 1);   // TOBLayer1_Z0
0239   m_color.push_back(kCyan - 9);   // TOBLayer1_Z20
0240   m_color.push_back(kCyan + 3);   // TOBLayer1_Z80
0241   m_color.push_back(kCyan + 4);   // TOBLayer1_Z90
0242   m_color.push_back(kAzure);      // TOBLayer2_Z0
0243   m_color.push_back(kAzure + 8);  // TOBLayer2_Z25
0244   m_color.push_back(kAzure + 2);  // TOBLayer2_Z80
0245   m_color.push_back(kAzure + 5);  // TOBLayer2_Z90
0246   m_color.push_back(kCyan + 1);   // TOBLayer3_Z0
0247   m_color.push_back(kCyan - 9);   // TOBLayer3_Z25
0248   m_color.push_back(kCyan + 3);   // TOBLayer3_Z80
0249   m_color.push_back(kCyan + 4);   // TOBLayer3_Z90
0250   m_color.push_back(kAzure);      // TOBLayer4_Z0
0251   m_color.push_back(kAzure + 8);  // TOBLayer4_Z25
0252   m_color.push_back(kAzure + 2);  // TOBLayer4_Z80
0253   m_color.push_back(kAzure + 5);  // TOBLayer4_Z90
0254   m_color.push_back(kCyan + 1);   // TOBLayer5_Z0
0255   m_color.push_back(kCyan - 9);   // TOBLayer5_Z25
0256   m_color.push_back(kCyan + 3);   // TOBLayer5_Z80
0257   m_color.push_back(kCyan + 4);   // TOBLayer5_Z90
0258 }
0259 
0260 std::vector<std::pair<std::shared_ptr<TLine>, std::shared_ptr<TText> > > ListGroups::overlayEtaReferences() {
0261   std::vector<std::pair<std::shared_ptr<TLine>, std::shared_ptr<TText> > > lines;
0262 
0263   lines.reserve(40);
0264   std::pair<float, float> deltaZ(293, 298);
0265   std::pair<float, float> deltaR(115, 118);
0266   float text_size = 0.033;
0267 
0268   for (float eta = 0.; eta <= 3.8; eta += 0.2) {
0269     float theta = 2. * atan(exp(-eta));
0270     if (eta >= 1.8) {
0271       lines.push_back(std::make_pair<std::shared_ptr<TLine>, std::shared_ptr<TText> >(
0272           std::make_shared<TLine>(deltaZ.first, deltaZ.first * tan(theta), deltaZ.second, deltaZ.second * tan(theta)),
0273           std::make_shared<TText>(deltaZ.first, deltaZ.first * tan(theta), fmt::sprintf("%2.1f", eta).c_str())));
0274       lines.back().second->SetTextFont(42);
0275       lines.back().second->SetTextSize(text_size);
0276       lines.back().second->SetTextAlign(33);
0277       lines.push_back(std::make_pair<std::shared_ptr<TLine>, std::shared_ptr<TText> >(
0278           std::make_shared<TLine>(-deltaZ.first, deltaZ.first * tan(theta), -deltaZ.second, deltaZ.second * tan(theta)),
0279           std::make_shared<TText>(-deltaZ.first, deltaZ.first * tan(theta), fmt::sprintf("-%2.1f", eta).c_str())));
0280       lines.back().second->SetTextFont(42);
0281       lines.back().second->SetTextSize(text_size);
0282       lines.back().second->SetTextAlign(13);
0283     } else {
0284       lines.push_back(std::make_pair<std::shared_ptr<TLine>, std::shared_ptr<TText> >(
0285           std::make_shared<TLine>(deltaR.first / tan(theta), deltaR.first, deltaR.second / tan(theta), deltaR.second),
0286           std::make_shared<TText>(deltaR.first / tan(theta), deltaR.first, fmt::sprintf("%2.1f", eta).c_str())));
0287       lines.back().second->SetTextFont(42);
0288       lines.back().second->SetTextSize(text_size);
0289       lines.back().second->SetTextAlign(23);
0290       if (eta != 0) {
0291         lines.push_back(std::make_pair<std::shared_ptr<TLine>, std::shared_ptr<TText> >(
0292             std::make_shared<TLine>(
0293                 -deltaR.first / tan(theta), deltaR.first, -deltaR.second / tan(theta), deltaR.second),
0294             std::make_shared<TText>(-deltaR.first / tan(theta), deltaR.first, fmt::sprintf("-%2.1f", eta).c_str())));
0295         lines.back().second->SetTextFont(42);
0296         lines.back().second->SetTextSize(text_size);
0297         lines.back().second->SetTextAlign(23);
0298       }
0299     }
0300   }
0301   return lines;
0302 }
0303 
0304 void ListGroups::produceAndSaveSummaryPlot(const edm::EventSetup &setup) {
0305   const double scale = 10.;
0306   std::vector<TText *> nukem_text;
0307   static int markerStyles[10] = {kFullCircle,
0308                                  kFullSquare,
0309                                  kFullTriangleUp,
0310                                  kFullTriangleDown,
0311                                  kOpenCircle,
0312                                  kOpenSquare,
0313                                  kOpenTriangleUp,
0314                                  kOpenDiamond,
0315                                  kOpenCross,
0316                                  kFullStar};
0317 
0318   edm::ESTransientHandle<DDCompactView> hDdd = setup.getTransientHandle(ddToken_);
0319 
0320   for (const auto &n : m_group_names) {
0321     m_groups.push_back(new MaterialAccountingGroup(n, *hDdd));
0322   };
0323 
0324   std::unique_ptr<TCanvas> canvas(
0325       new TCanvas("Grouping_rz", "Grouping - RZ view", (int)(600 * scale * 1.25), (int)(120 * scale * 1.50)));
0326   canvas->GetFrame()->SetFillColor(kWhite);
0327   gStyle->SetOptStat(0);
0328 
0329   unsigned int color_index = 1;
0330   // Setup the legend
0331   std::unique_ptr<TLegend> leg(new TLegend(0.1, 0.1, 0.23, 0.34));
0332   leg->SetHeader("Tracker Material Grouping");
0333   leg->SetTextFont(42);
0334   leg->SetTextSize(0.008);
0335   leg->SetNColumns(3);
0336   std::unique_ptr<TProfile2D> radlen(
0337       new TProfile2D("OverallRadLen", "OverallRadLen", 600., -300., 300, 120., 0., 120.));
0338   std::unique_ptr<TProfile2D> eneloss(
0339       new TProfile2D("OverallEnergyLoss", "OverallEnergyLoss", 600., -300., 300, 120., 0., 120.));
0340   std::unique_ptr<TProfile2D> radlen_diff(
0341       new TProfile2D("OverallDifferencesRadLen", "OverallDifferencesRadLen", 600., -300., 300, 120., 0., 120.));
0342   std::unique_ptr<TProfile2D> eneloss_diff(
0343       new TProfile2D("OverallDifferencesEnergyLoss", "OverallDifferencesEnergyLoss", 600., -300., 300, 120., 0., 120.));
0344 
0345   for (auto g : m_groups) {
0346     m_plots.push_back(
0347         new TH2F(g->name().c_str(), g->name().c_str(), 6000., -300., 300, 1200., 0., 120.));  // 10x10 points per cm2
0348     TH2F &current = *m_plots.back();
0349     current.SetMarkerColor(m_color[color_index]);
0350     current.SetMarkerStyle(markerStyles[color_index % 10]);
0351     current.SetMarkerSize(0.8);
0352     current.SetLineWidth(1);
0353     for (const auto &element : g->elements()) {
0354       current.Fill(element.z(), element.perp());
0355       radlen->Fill(element.z(), element.perp(), m_values[g->name()].first);
0356       eneloss->Fill(element.z(), element.perp(), m_values[g->name()].second);
0357       radlen_diff->Fill(element.z(), element.perp(), m_diff[g->name()].first);
0358       eneloss_diff->Fill(element.z(), element.perp(), m_diff[g->name()].second);
0359     }
0360 
0361     if (color_index == 1)
0362       current.Draw();
0363     else
0364       current.Draw("SAME");
0365 
0366     leg->AddEntry(&current, g->name().c_str(), "lp")->SetTextColor(m_color[color_index]);
0367     color_index++;
0368 
0369     // Loop over the same chromatic scale in case the number of
0370     // allocated colors is not big enough.
0371     color_index = color_index % m_color.size();
0372   }
0373   leg->Draw();
0374   canvas->SaveAs("Grouping.png");
0375 
0376   std::vector<std::pair<std::shared_ptr<TLine>, std::shared_ptr<TText> > > lines = overlayEtaReferences();
0377 
0378   canvas->Clear();
0379   radlen->SetMinimum(0);
0380   radlen->SetMaximum(0.25);
0381   radlen->Draw("COLZ");
0382   for (const auto &line : lines) {
0383     line.first->SetLineWidth(5);
0384     line.first->Draw();
0385     line.second->Draw();
0386   }
0387   canvas->SaveAs("RadLenValues.png");
0388 
0389   canvas->Clear();
0390   eneloss->SetMinimum(0.00001);
0391   eneloss->SetMaximum(0.0005);
0392   eneloss->Draw("COLZ");
0393   for (const auto &line : lines) {
0394     line.first->SetLineWidth(5);
0395     line.first->Draw();
0396     line.second->Draw();
0397   }
0398   canvas->SaveAs("EnergyLossValues.png");
0399 
0400   canvas->Clear();
0401   gStyle->SetPalette(m_gradient.size(), &m_gradient.front());
0402   gStyle->SetNumberContours(m_gradient.size());
0403   radlen_diff->SetMinimum(-100);
0404   radlen_diff->SetMaximum(100);
0405   radlen_diff->Draw("COLZ");
0406   for (const auto &line : lines) {
0407     line.first->SetLineWidth(5);
0408     line.first->Draw();
0409     line.second->Draw();
0410   }
0411   canvas->SaveAs("RadLenChanges.png");
0412 
0413   canvas->Clear();
0414   eneloss_diff->SetMinimum(-100);
0415   eneloss_diff->SetMaximum(100);
0416   eneloss_diff->Draw("COLZ");
0417   for (const auto &line : lines) {
0418     line.first->SetLineWidth(5);
0419     line.first->Draw();
0420     line.second->Draw();
0421   }
0422   canvas->SaveAs("EnergyLossChanges.png");
0423 
0424   for (auto g : nukem_text)
0425     delete g;
0426 }
0427 
0428 void ListGroups::analyze(const edm::Event &evt, const edm::EventSetup &setup) {
0429   edm::ESTransientHandle<DDCompactView> hDdd = setup.getTransientHandle(ddToken_);
0430 
0431   DDSpecificsHasNamedValueFilter filter{"TrackingMaterialGroup"};
0432   DDFilteredView fv(*hDdd, filter);
0433 
0434   while (fv.next()) {
0435     // print the group name and full hierarchy of all items
0436     std::cout << dddGetString("TrackingMaterialGroup", fv) << '\t';
0437     m_group_names.insert(dddGetString("TrackingMaterialGroup", fv));
0438 
0439     // start from 2 to skip the leading /OCMS[0]/CMSE[1] part
0440     const DDGeoHistory &history = fv.geoHistory();
0441     std::cout << '/';
0442     for (unsigned int h = 2; h < history.size(); ++h)
0443       std::cout << '/' << history[h].logicalPart().name().name() << '[' << history[h].copyno() << ']';
0444 
0445     // DD3Vector and DDTranslation are the same type as math::XYZVector
0446     math::XYZVector position = fv.translation() / 10.;  // mm -> cm
0447     std::cout << "\t(" << position.x() << ", " << position.y() << ", " << position.z() << ") "
0448               << "[rho] " << position.Rho() << std::endl;
0449   };
0450   std::cout << std::endl;
0451 
0452   if (m_saveSummaryPlot)
0453     produceAndSaveSummaryPlot(setup);
0454 }
0455 
0456 void ListGroups::endJob() {}
0457 
0458 //-------------------------------------------------------------------------
0459 // define as a plugin
0460 #include "FWCore/PluginManager/interface/ModuleDef.h"
0461 #include "FWCore/Framework/interface/MakerMacros.h"
0462 DEFINE_FWK_MODULE(ListGroups);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!