Hcal/test/visualizeHFPhase1PMTParams.cc

0001 #include <iostream>
0002 #include <cassert>
0003 #include <sstream>
0004 #include <cfloat>
0005 #include <algorithm>
0006 #include <memory>
0007
0008 #include "CondTools/Hcal/interface/visualizeHFPhase1PMTParams.h"
0009 #include "CondTools/Hcal/interface/parseHcalDetId.h"
0010
0011 #include "TCanvas.h"
0012 #include "TPad.h"
0013 #include "TGraph.h"
0014 #include "TAxis.h"
0015 #include "TH1F.h"
0016
0017 typedef std::vector<std::shared_ptr<TGraph> > Garbage;
0018
0019 //
0020 // Return a usable root name (without spaces) for the given detector id
0021 //
0022 static std::string rootDetectorName(const HcalDetId& id) {
0023   using namespace std;
0024   ostringstream os;
0025   os << hcalSubdetectorName(id.subdet()) << '_' << id.ieta() << '_' << id.iphi() << '_' << id.depth();
0026   return os.str();
0027 }
0028
0029 static TH1F* plotFunctors(TVirtualPad* pad,
0030                           const char* plotTitle,
0031                           const char* yAxisTitle,
0032                           const double cut,
0033                           const double ymin,
0034                           const double ymax,
0035                           const AbsHcalFunctor& f1,
0036                           const AbsHcalFunctor& f2,
0037                           const VisualizationOptions& opts,
0038                           Garbage& garbage,
0039                           const bool isReference,
0040                           TH1F* frame) {
0041   const int lineWidth = 2;
0042   const int refLineWidth = 4;
0043
0044   assert(opts.plotPoints > 1);
0045   const unsigned nDraw = opts.plotPoints + 1;
0046   std::vector<double> xvec(nDraw);
0047   std::vector<double> y1(nDraw);
0048   std::vector<double> y2(nDraw);
0049   const double xrange = opts.maxCharge - cut;
0050   const double yrange = ymax - ymin;
0051   xvec[0] = cut;
0052   y1[0] = ymin + 0.05 * yrange;
0053   y2[0] = ymax - 0.05 * yrange;
0054   const double h = xrange / (opts.plotPoints - 1);
0055   for (unsigned i = 1; i <= opts.plotPoints; ++i) {
0056     const double x = cut + h * (i - 1);
0057     xvec[i] = x;
0058     y1[i] = f1(x);
0059     y2[i] = f2(x);
0060   }
0061
0062   // Root was designed by (insert your favorite explective here),
0063   // so the plot frame on the pad is of type TH1F (a histogram!)
0064   //
0065   if (frame == nullptr) {
0066     frame = pad->DrawFrame(opts.minCharge, ymin, opts.maxCharge + 0.1 * xrange, ymax);
0067     frame->GetXaxis()->SetTitle("Q [fC]");
0068     frame->GetYaxis()->SetTitle(yAxisTitle);
0069     frame->SetTitle(plotTitle);
0070   }
0071
0072   std::shared_ptr<TGraph> gr1(new TGraph(nDraw, &xvec[0], &y1[0]));
0073   gr1->SetLineWidth(isReference ? refLineWidth : lineWidth);
0074   gr1->SetLineColor(isReference ? kGreen : kBlue);
0075   gr1->Draw("L");
0076   garbage.push_back(gr1);
0077
0078   std::shared_ptr<TGraph> gr2(new TGraph(nDraw, &xvec[0], &y2[0]));
0079   gr2->SetLineWidth(isReference ? refLineWidth : lineWidth);
0080   gr2->SetLineColor(isReference ? kGreen : kRed);
0081   gr2->Draw("L");
0082   garbage.push_back(gr2);
0083
0084   return frame;
0085 }
0086
0087 static void plotConfig(const HcalDetId& id,
0088                        const unsigned tableIndex,
0089                        const HFPhase1PMTData& cut,
0090                        const VisualizationOptions& opts,
0091                        const HFPhase1PMTData* ref) {
0092   using namespace std;
0093
0094   const double canvasWidth = 1000;
0095   const double canvasHeight = 300;
0096
0097   // Root was designed by idiots oblivious to smart pointers and unable
0098   // to carry their objects around. Because of this, most plottable
0099   // object in "Root" must have a unique name, even if it makes no sense
0100   // whatsoever to name such an object. Keep it in mind that idiots are
0101   // incapable of handling spaces inside the names.
0102   //
0103   const bool isDefault = tableIndex == HcalIndexLookup::InvalidIndex;
0104   const std::string& detName = rootDetectorName(id);
0105   const std::string& canvName = isDefault ? std::string("Default") : detName;
0106   ostringstream ostitle;
0107   if (isDefault)
0108     ostitle << "Default Cuts";
0109   else
0110     ostitle << "Cuts for " << detName;
0111   const std::string& canvTitle = ostitle.str();
0112
0113   TCanvas* canv = new TCanvas(canvName.c_str(), canvTitle.c_str(), canvasWidth, canvasHeight);
0114
0115   // Root was designed by idiots oblivious to the difference between
0116   // window size and canvas size. Because of this, the canvas size
0117   // has to be set explicitly again, despite the fact that the
0118   // constructor already had these arguments.
0119   //
0120   canv->SetCanvasSize(canvasWidth, canvasHeight);
0121   canv->Divide(3, 1);
0122
0123   // Root was designed by idiots changing their minds too often.
0124   // Graphs are an exception from the naming scheme and ownership
0125   // convention. They are not managed by root, so graphs on all
0126   // pads must persist until the whole canvas is written out.
0127   // This is why we need a garbage bin -- to manage these things.
0128   //
0129   Garbage bin;
0130
0131   TVirtualPad* pad = canv->cd(1);
0132   {
0133     TH1F* frame = nullptr;
0134     if (ref)
0135       frame = plotFunctors(pad,
0136                            "PMT Anode 0",
0137                            "Rise Time Cuts [ns]",
0138                            ref->minCharge0(),
0139                            opts.minTDC,
0140                            opts.maxTDC,
0141                            ref->cut(HFPhase1PMTData::T_0_MIN),
0142                            ref->cut(HFPhase1PMTData::T_0_MAX),
0143                            opts,
0144                            bin,
0145                            true,
0146                            frame);
0147     plotFunctors(pad,
0148                  "PMT Anode 0",
0149                  "Rise Time Cuts [ns]",
0150                  cut.minCharge0(),
0151                  opts.minTDC,
0152                  opts.maxTDC,
0153                  cut.cut(HFPhase1PMTData::T_0_MIN),
0154                  cut.cut(HFPhase1PMTData::T_0_MAX),
0155                  opts,
0156                  bin,
0157                  false,
0158                  frame);
0159   }
0160   pad->Draw();
0161
0162   ostringstream midtitle;
0163   midtitle << "PMT Anode 1 ";
0164   if (isDefault)
0165     midtitle << "(Default Cuts)";
0166   else
0167     midtitle << '(' << detName << ", Index " << tableIndex << ')';
0168   const string& mtstr = midtitle.str();
0169   pad = canv->cd(2);
0170   {
0171     TH1F* frame = nullptr;
0172     if (ref)
0173       frame = plotFunctors(pad,
0174                            mtstr.c_str(),
0175                            "Rise Time Cuts [ns]",
0176                            ref->minCharge1(),
0177                            opts.minTDC,
0178                            opts.maxTDC,
0179                            ref->cut(HFPhase1PMTData::T_1_MIN),
0180                            ref->cut(HFPhase1PMTData::T_1_MAX),
0181                            opts,
0182                            bin,
0183                            true,
0184                            frame);
0185     plotFunctors(pad,
0186                  mtstr.c_str(),
0187                  "Rise Time Cuts [ns]",
0188                  cut.minCharge1(),
0189                  opts.minTDC,
0190                  opts.maxTDC,
0191                  cut.cut(HFPhase1PMTData::T_1_MIN),
0192                  cut.cut(HFPhase1PMTData::T_1_MAX),
0193                  opts,
0194                  bin,
0195                  false,
0196                  frame);
0197   }
0198   pad->Draw();
0199
0200   pad = canv->cd(3);
0201   {
0202     TH1F* frame = nullptr;
0203     if (ref)
0204       frame = plotFunctors(pad,
0205                            "Charge Asymmetry",
0206                            "Asymmetry Cuts",
0207                            ref->minChargeAsymm(),
0208                            opts.minAsymm,
0209                            opts.maxAsymm,
0210                            ref->cut(HFPhase1PMTData::ASYMM_MIN),
0211                            ref->cut(HFPhase1PMTData::ASYMM_MAX),
0212                            opts,
0213                            bin,
0214                            true,
0215                            frame);
0216     plotFunctors(pad,
0217                  "Charge Asymmetry",
0218                  "Asymmetry Cuts",
0219                  cut.minChargeAsymm(),
0220                  opts.minAsymm,
0221                  opts.maxAsymm,
0222                  cut.cut(HFPhase1PMTData::ASYMM_MIN),
0223                  cut.cut(HFPhase1PMTData::ASYMM_MAX),
0224                  opts,
0225                  bin,
0226                  false,
0227                  frame);
0228   }
0229   pad->Draw();
0230
0231   canv->Write();
0232 }
0233
0234 void visualizeHFPhase1PMTParams(const std::vector<HcalDetId>& idVec,
0235                                 const HFPhase1PMTParams& cuts,
0236                                 const VisualizationOptions& options,
0237                                 const HFPhase1PMTParams* reference) {
0238   using namespace std;
0239
0240   const unsigned n = idVec.size();
0241   if (options.verbose) {
0242     cout << "Cut visualization is requested for the following PMTs:\n";
0243     for (unsigned i = 0; i < n; ++i)
0244       cout << idVec[i] << '\n';
0245     cout.flush();
0246   }
0247
0248   // Plot the default cut
0249   const HFPhase1PMTData* defaultCut = cuts.getDefault();
0250   const HFPhase1PMTData* defaultRef = nullptr;
0251   if (reference)
0252     defaultRef = reference->getDefault();
0253   if (defaultCut)
0254     plotConfig(HcalDetId(), HcalIndexLookup::InvalidIndex, *defaultCut, options, defaultRef);
0255
0256   // Plot all other cuts
0257   for (unsigned i = 0; i < n; ++i) {
0258     const HcalDetId& id = idVec[i];
0259     const unsigned tableIndex = cuts.getIndex(id);
0260     if (tableIndex == HcalIndexLookup::InvalidIndex) {
0261       if (defaultCut)
0262         cout << "PMT " << id << " is using default config" << endl;
0263       else
0264         cout << "ERROR! No configuration found for PMT " << id << endl;
0265     } else {
0266       cout << "PMT " << id << " is using config " << tableIndex << endl;
0267       const HFPhase1PMTData* cut = cuts.getByIndex(tableIndex);
0268       assert(cut);
0269       const HFPhase1PMTData* ref = nullptr;
0270       if (reference)
0271         ref = reference->getByIndex(tableIndex);
0272       plotConfig(id, tableIndex, *cut, options, ref);
0273     }
0274   }
0275 }