Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-10-25 09:43:10

0001 // -*- C++ -*-
0002 //
0003 // Package:     SiPixelPhase1Common
0004 // Class  :     HistogramManager
0005 //
0006 #include "DQM/SiPixelPhase1Common/interface/HistogramManager.h"
0007 
0008 #include <sstream>
0009 #include <boost/algorithm/string.hpp>
0010 
0011 // Geometry stuff
0012 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0013 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0014 
0015 // Logger
0016 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0017 
0018 HistogramManager::HistogramManager(const edm::ParameterSet& iconfig, GeometryInterface& geo)
0019     : geometryInterface(geo),
0020       enabled(iconfig.getParameter<bool>("enabled")),
0021       perLumiHarvesting(iconfig.getParameter<bool>("perLumiHarvesting")),
0022       bookUndefined(iconfig.getParameter<bool>("bookUndefined")),
0023       top_folder_name(iconfig.getParameter<std::string>("topFolderName")),
0024       name(iconfig.getParameter<std::string>("name")),
0025       title(iconfig.getParameter<std::string>("title")),
0026       xlabel(iconfig.getParameter<std::string>("xlabel")),
0027       ylabel(iconfig.getParameter<std::string>("ylabel")),
0028       dimensions(iconfig.getParameter<int>("dimensions")),
0029       range_x_nbins(iconfig.getParameter<int>("range_nbins")),
0030       range_x_min(iconfig.getParameter<double>("range_min")),
0031       range_x_max(iconfig.getParameter<double>("range_max")),
0032       range_y_nbins(iconfig.getParameter<int>("range_y_nbins")),
0033       range_y_min(iconfig.getParameter<double>("range_y_min")),
0034       range_y_max(iconfig.getParameter<double>("range_y_max")),
0035       statsOverflows(iconfig.getParameter<bool>("statsOverflows")) {
0036   auto spec_configs = iconfig.getParameter<edm::VParameterSet>("specs");
0037   for (const auto& spec : spec_configs) {
0038     // this would fit better in SummationSpecification(...), but it has to
0039     // happen here.
0040     auto conf = spec.getParameter<edm::ParameterSet>("conf");
0041     if (!conf.getParameter<bool>("enabled"))
0042       continue;
0043     addSpec(SummationSpecification(spec, geometryInterface));
0044   }
0045 }
0046 
0047 void HistogramManager::addSpec(SummationSpecification spec) {
0048   specs.push_back(spec);
0049   tables.push_back(Table());
0050   counters.push_back(Table());
0051   significantvalues.push_back(GeometryInterface::Values());
0052   fastpath.push_back(nullptr);
0053 }
0054 
0055 // This is the hottest function in the HistogramManager. Make sure it does not
0056 // allocate memory or other expensive things.
0057 // Currently the GeometryInterface::extract (some virtual calls) and the map
0058 // lookups should be the most expensive things, but they should only happen if
0059 // the module changed from the last call; an optimization that fails when row/
0060 // col are used.
0061 // fillInternal (called from here) does more lookups on the geometry, if EXTEND
0062 // is used; we do not attempt the optimization there since in most cases row/
0063 // col are involved.
0064 void HistogramManager::fill(double x, double y, DetId sourceModule, const edm::Event* sourceEvent, int col, int row) {
0065   if (!enabled)
0066     return;
0067 
0068   // We only need to check the module here, since the fastpath is only used to
0069   // determine which plot is filled (not which bin inside) and fillInternal
0070   // re-extracts whatever it needs for the axis.
0071   // Since do not support booking based on ROC or time, the plot can only depend
0072   // on the module.
0073   // The sourceEvent check is not really effective (pointer is always the same)
0074   // but needed for initialisation.
0075   // NOTE that this might change if we want to support per-ROC booking
0076   // NOTE that we could even cache the bin to fill if iq and spec are identical,
0077   // also across HistogramManagers.
0078   bool cached = false;
0079   if (sourceModule == this->iq.sourceModule && sourceEvent == this->iq.sourceEvent) {
0080     cached = true;
0081   }
0082   iq = GeometryInterface::InterestingQuantities{sourceEvent, sourceModule, int16_t(col), int16_t(row)};
0083   for (unsigned int i = 0; i < specs.size(); i++) {
0084     auto& s = specs[i];
0085     auto& t = s.steps[0].type == SummationStep::COUNT ? counters[i] : tables[i];
0086     if (!cached) {  // recompute ME to fill (aka fastpath)
0087       significantvalues[i].clear();
0088       geometryInterface.extractColumns(s.steps[0].columns, iq, significantvalues[i]);
0089       auto histo = t.find(significantvalues[i]);
0090       if (histo == t.end()) {
0091         if (bookUndefined) {
0092           edm::LogError("HistogramManager") << "Missing Histogram!\n"
0093                                             << "name " << tables[i].begin()->second.th1->GetName() << "\n";
0094           assert(!"Histogram not booked! Probably inconsistent geometry description.");
0095         }
0096         fastpath[i] = nullptr;
0097       } else {
0098         fastpath[i] = &(histo->second);
0099       }
0100     }
0101     // A fastpath of nullptr means there is no ME for this iq, which can be
0102     // a valid cached result.
0103     if (fastpath[i]) {
0104       if (s.steps[0].type == SummationStep::COUNT) {
0105         fastpath[i]->count++;
0106       } else {
0107         fillInternal(x, y, this->dimensions, iq, s.steps.begin() + 1, s.steps.end(), *(fastpath[i]));
0108       }
0109     }
0110   }
0111 }
0112 void HistogramManager::fill(double x, DetId sourceModule, const edm::Event* sourceEvent, int col, int row) {
0113   assert(this->dimensions == 1);
0114   fill(x, 0.0, sourceModule, sourceEvent, col, row);
0115 }
0116 void HistogramManager::fill(DetId sourceModule, const edm::Event* sourceEvent, int col, int row) {
0117   assert(this->dimensions == 0);
0118   fill(0.0, 0.0, sourceModule, sourceEvent, col, row);
0119 }
0120 
0121 void HistogramManager::fillInternal(double x,
0122                                     double y,
0123                                     int n_parameters,
0124                                     GeometryInterface::InterestingQuantities const& iq,
0125                                     std::vector<SummationStep>::iterator first,
0126                                     std::vector<SummationStep>::iterator last,
0127                                     AbstractHistogram& dest) {
0128   double fx = 0, fy = 0, fz = 0;
0129   int tot_parameters = n_parameters;
0130   for (auto it = first; it != last; ++it) {
0131     if (it->stage != SummationStep::STAGE1)
0132       break;
0133     // The specification builder precomputes where x and y go, this loop will
0134     // always do 3 iterations to set x, y, z. The builder does not know how
0135     // many parameters we have, so we have to check that and count the total.
0136     switch (it->type) {
0137       case SummationStep::USE_X:
0138         if (it->arg[0] == '1' && n_parameters >= 1)
0139           fx = x;
0140         if (it->arg[0] == '2' && n_parameters >= 2)
0141           fx = y;
0142         break;
0143       case SummationStep::USE_Y:
0144         if (it->arg[0] == '1' && n_parameters >= 1)
0145           fy = x;
0146         if (it->arg[0] == '2' && n_parameters >= 2)
0147           fy = y;
0148         break;
0149       case SummationStep::USE_Z:
0150         if (it->arg[0] == '1' && n_parameters >= 1)
0151           fz = x;
0152         if (it->arg[0] == '2' && n_parameters >= 2)
0153           fz = y;
0154         break;
0155       case SummationStep::EXTEND_X:
0156         fx = geometryInterface.extract(it->columns[0], iq).second;
0157         tot_parameters++;
0158         break;
0159       case SummationStep::EXTEND_Y:
0160         fy = geometryInterface.extract(it->columns[0], iq).second;
0161         tot_parameters++;
0162         break;
0163       case SummationStep::PROFILE:
0164         break;  // profile does not make a difference here, only in booking
0165       default:
0166         assert(!"illegal step in STAGE1!");
0167     }
0168   }
0169 
0170   switch (tot_parameters) {
0171     case 1:
0172       dest.me->Fill(fx);
0173       break;
0174     case 2:
0175       dest.me->Fill(fx, fy);
0176       break;
0177     case 3:
0178       dest.me->Fill(fx, fy, fz);
0179       break;
0180     default:
0181       edm::LogError("HistogramManager") << "got " << tot_parameters << " dimensions\n"
0182                                         << "name " << dest.th1->GetName() << "\n";
0183       assert(!"More than 3 dimensions should never occur.");
0184   }
0185 }
0186 
0187 // This is only used for ndigis-like counting. It could be more optimized, but
0188 // is probably fine for a per-event thing.
0189 void HistogramManager::executePerEventHarvesting(const edm::Event* sourceEvent) {
0190   if (!enabled)
0191     return;
0192   for (unsigned int i = 0; i < specs.size(); i++) {
0193     auto& s = specs[i];
0194     auto& t = tables[i];
0195     auto& c = counters[i];
0196     if (s.steps[0].type != SummationStep::COUNT)
0197       continue;  // no counting, done
0198     assert((s.steps.size() >= 2 && s.steps[1].type == SummationStep::GROUPBY) ||
0199            !"Incomplete spec (but this cannot be caught in Python)");
0200     for (auto& e : c) {
0201       // the iq on the counter can only be a _sample_, since many modules
0202       // could be grouped on one counter. Even worse, the sourceEvent ptr
0203       // could be dangling if the counter was not touched in this event, so
0204       // we replace it. row/col are most likely useless as well.
0205       auto iq = e.second.iq_sample;
0206       iq.sourceEvent = sourceEvent;
0207 
0208       significantvalues[i].clear();
0209       geometryInterface.extractColumns(s.steps[1].columns, iq, significantvalues[i]);
0210       auto histo = t.find(significantvalues[i]);
0211       if (histo == t.end()) {
0212         if (!bookUndefined)
0213           continue;
0214         edm::LogError("HistogramManager") << "Histogram Missing!\n"
0215                                           << "name " << t.begin()->second.th1->GetName() << "\n"
0216                                           << "ctr "
0217                                           << " detid " << iq.sourceModule << "\n";
0218         assert(!"Histogram not booked! (per-event) Probably inconsistent geometry description.");
0219       }
0220       fillInternal(e.second.count, 0, 1, iq, s.steps.begin() + 2, s.steps.end(), histo->second);
0221       e.second.count = 0;
0222     }
0223   }
0224 }
0225 
0226 std::pair<std::string, std::string> HistogramManager::makePathName(SummationSpecification const& s,
0227                                                                    GeometryInterface::Values const& significantvalues,
0228                                                                    SummationStep const* upto) {
0229   std::ostringstream dir("");
0230   std::string suffix = "";
0231 
0232   // we omit the last value here, to get all disks next to each other etc.
0233   if (!significantvalues.empty()) {
0234     for (auto it = significantvalues.begin(); it != (significantvalues.end() - 1); ++it) {
0235       auto name = geometryInterface.formatValue(it->first, it->second);
0236       if (name.empty())
0237         continue;
0238       dir << name << "/";
0239     }
0240     auto e = significantvalues[significantvalues.size() - 1];
0241     suffix = "_" + geometryInterface.formatValue(e.first, e.second);
0242   }
0243 
0244   // PERF: this is actually independent of significantvalues and iq
0245   std::string name = this->name;
0246   for (SummationStep const& step : s.steps) {
0247     if (&step == upto)
0248       break;
0249     switch (step.type) {
0250       case SummationStep::COUNT:
0251         name = "num_" + name;
0252         break;
0253       case SummationStep::EXTEND_X:
0254       case SummationStep::EXTEND_Y: {
0255         if (step.stage != SummationStep::STAGE1)
0256           break;
0257         GeometryInterface::Column col0 = step.columns[0];
0258         std::string colname = geometryInterface.pretty(col0);
0259         name = name + "_per_" + colname;
0260         break;
0261       }
0262       case SummationStep::REDUCE: {
0263         auto red = step.arg;
0264         boost::algorithm::to_lower(red);
0265         name = red + "_" + name;
0266         break;
0267       }
0268       default:
0269         // Maybe PROFILE is worth showing.
0270         break;  // not visible in name
0271     }
0272   }
0273   return std::make_pair(top_folder_name + "/" + dir.str(), name + suffix);
0274 }
0275 
0276 void HistogramManager::book(DQMStore::IBooker& iBooker, edm::EventSetup const& iSetup) {
0277   if (!geometryInterface.loaded()) {
0278     geometryInterface.load(iSetup);
0279   }
0280   if (!enabled)
0281     return;
0282 
0283   struct MEInfo {
0284     int dimensions = 1;
0285     double range_x_min = 1e12;
0286     double range_x_max = -1e12;
0287     double range_y_min = 1e12;
0288     double range_y_max = -1e12;
0289     double range_z_min = 1e12;  // z range carried around but unused
0290     double range_z_max = -1e12;
0291     int range_x_nbins = 0;
0292     int range_y_nbins = 0;
0293     int range_z_nbins = 0;
0294     GeometryInterface::Value binwidth_x = 0;  // override nbins for geom-things
0295     GeometryInterface::Value binwidth_y = 0;
0296     std::string title, xlabel, ylabel, zlabel;
0297     bool do_profile = false;
0298     bool statsOverflows = true;
0299     ;
0300   };
0301   std::map<GeometryInterface::Values, MEInfo> toBeBooked;
0302 
0303   for (unsigned int i = 0; i < specs.size(); i++) {
0304     auto& s = specs[i];
0305     auto& t = tables[i];
0306     auto& c = counters[i];
0307     toBeBooked.clear();
0308     bool bookCounters = false;
0309 
0310     auto firststep = s.steps.begin();
0311     int n_parameters = this->dimensions;
0312 
0313     if (firststep->type != SummationStep::GROUPBY) {
0314       ++firststep;
0315       n_parameters = 1;
0316       bookCounters = true;
0317     }
0318 
0319     auto laststep = std::find_if(
0320         s.steps.begin(), s.steps.end(), [](SummationStep const& step) { return step.stage == SummationStep::STAGE2; });
0321 
0322     GeometryInterface::Values significantvalues;
0323 
0324     for (auto iq : geometryInterface.allModules()) {
0325       if (bookCounters) {
0326         // add an entry for the counting step if present
0327         geometryInterface.extractColumns(s.steps[0].columns, iq, significantvalues);
0328         c[significantvalues].iq_sample = iq;
0329       }
0330 
0331       geometryInterface.extractColumns(firststep->columns, iq, significantvalues);
0332       if (!bookUndefined) {
0333         // skip if any column is UNDEFINED
0334         bool ok = true;
0335         for (auto e : significantvalues)
0336           if (e.second == GeometryInterface::UNDEFINED)
0337             ok = false;
0338         if (!ok)
0339           continue;
0340       }
0341 
0342       auto histo = toBeBooked.find(significantvalues);
0343       if (histo == toBeBooked.end()) {
0344         // create new histo
0345         MEInfo& mei = toBeBooked[significantvalues];
0346         mei.title = this->title;
0347         mei.statsOverflows = this->statsOverflows;
0348         if (bookCounters)
0349           mei.title =
0350               "Number of " + mei.title + " per Event and " + geometryInterface.pretty(*(s.steps[0].columns.end() - 1));
0351         std::string xlabel = bookCounters ? "#" + this->xlabel : this->xlabel;
0352 
0353         // refer to fillInternal() for the actual execution
0354         // compute labels, title, type, user-set ranges here
0355         int tot_parameters = n_parameters;
0356 
0357 #define SET_AXIS(to, from)                                                          \
0358   mei.to##label = from##label;                                                      \
0359   mei.range_##to##_min = ((it->nbins == -1) ? this->range_##from##_min : it->xmin); \
0360   mei.range_##to##_max = ((it->nbins == -1) ? this->range_##from##_max : it->xmax); \
0361   mei.range_##to##_nbins = ((it->nbins == -1) ? this->range_##from##_nbins : it->nbins)
0362 
0363         for (auto it = firststep + 1; it != laststep; ++it) {
0364           switch (it->type) {
0365             case SummationStep::USE_X:
0366               if (it->arg[0] == '1' && n_parameters >= 1) {
0367                 SET_AXIS(x, x);
0368               }  // TODO: make use of current nbins, xmin, xmax if set
0369               if (it->arg[0] == '2' && n_parameters >= 2) {
0370                 SET_AXIS(x, y);
0371               }
0372               break;
0373             case SummationStep::USE_Y:
0374               if (it->arg[0] == '1' && n_parameters >= 1) {
0375                 SET_AXIS(y, x);
0376               }
0377               if (it->arg[0] == '2' && n_parameters >= 2) {
0378                 SET_AXIS(y, y);
0379               }
0380               break;
0381             case SummationStep::USE_Z:
0382               if (it->arg[0] == '1' && n_parameters >= 1) {
0383                 SET_AXIS(z, x);
0384               }
0385               if (it->arg[0] == '2' && n_parameters >= 2) {
0386                 SET_AXIS(z, y);
0387               }
0388               break;
0389             case SummationStep::EXTEND_X: {
0390               assert(mei.range_x_nbins == 0);
0391               auto col = it->columns[0];
0392               mei.xlabel = geometryInterface.pretty(col);
0393               mei.title = mei.title + " by " + mei.xlabel;
0394               if (geometryInterface.minValue(col) != GeometryInterface::UNDEFINED)
0395                 mei.range_x_min = geometryInterface.minValue(col);
0396               if (geometryInterface.maxValue(col) != GeometryInterface::UNDEFINED)
0397                 mei.range_x_max = geometryInterface.maxValue(col);
0398               mei.binwidth_x = geometryInterface.binWidth(col);
0399               tot_parameters++;
0400             } break;
0401             case SummationStep::EXTEND_Y: {
0402               auto col = it->columns[0];
0403               mei.ylabel = geometryInterface.pretty(col);
0404               mei.title = mei.title + " by " + mei.ylabel;
0405               if (geometryInterface.minValue(col) != GeometryInterface::UNDEFINED)
0406                 mei.range_y_min = geometryInterface.minValue(col);
0407               if (geometryInterface.maxValue(col) != GeometryInterface::UNDEFINED)
0408                 mei.range_y_max = geometryInterface.maxValue(col);
0409               mei.binwidth_y = geometryInterface.binWidth(col);
0410               tot_parameters++;
0411             } break;
0412             case SummationStep::PROFILE:
0413               mei.do_profile = true;
0414               break;
0415             default:
0416               assert(!"illegal step in STAGE1! (booking)");
0417           }
0418         }
0419         mei.dimensions = tot_parameters;
0420         if (mei.do_profile)
0421           mei.title = "Profile of " + mei.title;
0422         if (!mei.zlabel.empty())
0423           mei.title = mei.title + " (Z: " + mei.zlabel + ")";
0424       }
0425       // only update range
0426       MEInfo& mei = toBeBooked[significantvalues];
0427       double val;
0428 
0429       for (auto it = firststep + 1; it != laststep; ++it) {
0430         switch (it->type) {
0431           case SummationStep::EXTEND_X:
0432             val = geometryInterface.extract(it->columns[0], iq).second;
0433             if (val == GeometryInterface::UNDEFINED)
0434               break;
0435             mei.range_x_min = std::min(mei.range_x_min, val);
0436             mei.range_x_max = std::max(mei.range_x_max, val);
0437             break;
0438           case SummationStep::EXTEND_Y:
0439             val = geometryInterface.extract(it->columns[0], iq).second;
0440             if (val == GeometryInterface::UNDEFINED)
0441               break;
0442             mei.range_y_min = std::min(mei.range_y_min, val);
0443             mei.range_y_max = std::max(mei.range_y_max, val);
0444             break;
0445           default:  // ignore the rest, code above will catch bugs
0446             break;
0447         }
0448       }
0449     }
0450 
0451     // Now do the actual booking.
0452     for (auto& e : toBeBooked) {
0453       AbstractHistogram& h = t[e.first];
0454       MEInfo& mei = e.second;
0455       auto name = makePathName(s, e.first, &(*laststep));
0456       iBooker.setCurrentFolder(name.first);
0457 
0458       // determine nbins for geometry derived quantities
0459       // due to how we counted above, we need to include lower and upper bound
0460       // For Coord-values, which are not precisely aligned with the bins, we force
0461       // alignment.
0462       if (mei.binwidth_x != 0) {
0463         double range = (mei.range_x_max - mei.range_x_min) / mei.binwidth_x;
0464         if ((range - int(range)) == 0.0) {
0465           mei.range_x_min -= mei.binwidth_x / 2;
0466           mei.range_x_max += mei.binwidth_x / 2;
0467         } else {
0468           mei.range_x_min = std::floor(mei.range_x_min / mei.binwidth_x) * mei.binwidth_x;
0469           mei.range_x_max = std::ceil(mei.range_x_max / mei.binwidth_x) * mei.binwidth_x;
0470         }
0471         mei.range_x_nbins = int((mei.range_x_max - mei.range_x_min) / mei.binwidth_x);
0472       }
0473       if (mei.binwidth_y != 0) {
0474         double range = (mei.range_y_max - mei.range_y_min) / mei.binwidth_y;
0475         if ((range - int(range)) == 0.0) {
0476           mei.range_y_min -= mei.binwidth_y / 2;
0477           mei.range_y_max += mei.binwidth_y / 2;
0478         } else {
0479           mei.range_y_min = std::floor(mei.range_y_min / mei.binwidth_y) * mei.binwidth_y;
0480           mei.range_y_max = std::ceil(mei.range_y_max / mei.binwidth_y) * mei.binwidth_y;
0481         }
0482         mei.range_y_nbins = int((mei.range_y_max - mei.range_y_min) / mei.binwidth_y);
0483       }
0484 
0485       if (mei.dimensions == 1) {
0486         h.me = iBooker.book1D(
0487             name.second, (mei.title + ";" + mei.xlabel).c_str(), mei.range_x_nbins, mei.range_x_min, mei.range_x_max);
0488       } else if (mei.dimensions == 2 && !mei.do_profile) {
0489         h.me = iBooker.book2D(name.second,
0490                               (mei.title + ";" + mei.xlabel + ";" + mei.ylabel).c_str(),
0491                               mei.range_x_nbins,
0492                               mei.range_x_min,
0493                               mei.range_x_max,
0494                               mei.range_y_nbins,
0495                               mei.range_y_min,
0496                               mei.range_y_max);
0497       } else if (mei.dimensions == 2 && mei.do_profile) {
0498         h.me = iBooker.bookProfile(name.second,
0499                                    (mei.title + ";" + mei.xlabel + ";" + mei.ylabel).c_str(),
0500                                    mei.range_x_nbins,
0501                                    mei.range_x_min,
0502                                    mei.range_x_max,
0503                                    0.0,
0504                                    0.0,
0505                                    "");
0506       } else if (mei.dimensions == 3 && mei.do_profile) {
0507         h.me = iBooker.bookProfile2D(name.second,
0508                                      (mei.title + ";" + mei.xlabel + ";" + mei.ylabel).c_str(),
0509                                      mei.range_x_nbins,
0510                                      mei.range_x_min,
0511                                      mei.range_x_max,
0512                                      mei.range_y_nbins,
0513                                      mei.range_y_min,
0514                                      mei.range_y_max,
0515                                      0.0,
0516                                      0.0);  // Z range is ignored if min==max
0517       } else {
0518         edm::LogError("HistogramManager") << "Illegal Histogram!\n"
0519                                           << "name " << name.second << "\n"
0520                                           << "dim " << mei.dimensions << " profile " << mei.do_profile << "\n";
0521         assert(!"Illegal Histogram kind.");
0522       }
0523       h.th1 = h.me->getTH1();
0524       h.me->setStatOverflows(mei.statsOverflows);
0525     }
0526   }
0527 }
0528 
0529 void HistogramManager::executePerLumiHarvesting(DQMStore::IBooker& iBooker,
0530                                                 DQMStore::IGetter& iGetter,
0531                                                 edm::LuminosityBlock const& lumiBlock,
0532                                                 edm::EventSetup const& iSetup) {
0533   if (!enabled)
0534     return;
0535   // this should also give us the GeometryInterface for offline, though it is a
0536   // bit dirty and might explode.
0537   if (!geometryInterface.loaded()) {
0538     geometryInterface.load(iSetup);
0539   }
0540   if (perLumiHarvesting) {
0541     this->lumisection = &lumiBlock;  // "custom" steps can use this
0542     executeHarvesting(iBooker, iGetter);
0543     this->lumisection = nullptr;
0544   }
0545 }
0546 
0547 void HistogramManager::loadFromDQMStore(SummationSpecification& s, Table& t, DQMStore::IGetter& iGetter) {
0548   t.clear();
0549   GeometryInterface::Values significantvalues;
0550   auto firststep = s.steps.begin();
0551   if (firststep->type != SummationStep::GROUPBY)
0552     ++firststep;
0553   auto laststep = std::find_if(
0554       s.steps.begin(), s.steps.end(), [](SummationStep const& step) { return step.stage == SummationStep::STAGE2; });
0555 
0556   for (auto iq : geometryInterface.allModules()) {
0557     geometryInterface.extractColumns(firststep->columns, iq, significantvalues);
0558 
0559     auto histo = t.find(significantvalues);
0560     if (histo == t.end()) {
0561       auto name = makePathName(s, significantvalues, &(*laststep));
0562       std::string path = name.first + name.second;
0563       MonitorElement* me = iGetter.get(path);
0564       if (!me) {
0565         if (bookUndefined)
0566           edm::LogError("HistogramManager") << "ME " << path << " not found\n";
0567         // else this will happen quite often
0568       } else {
0569         AbstractHistogram& histo = t[significantvalues];
0570         histo.me = me;
0571         histo.th1 = histo.me->getTH1();
0572         histo.iq_sample = iq;
0573       }
0574     }
0575   }
0576 }
0577 
0578 void HistogramManager::executeGroupBy(SummationStep const& step,
0579                                       Table& t,
0580                                       DQMStore::IBooker& iBooker,
0581                                       SummationSpecification const& s) {
0582   // Simple regrouping, sum histos if they end up in the same place.
0583   Table out;
0584   GeometryInterface::Values significantvalues;
0585   for (auto& e : t) {
0586     TH1* th1 = e.second.th1;
0587     geometryInterface.extractColumns(step.columns, e.second.iq_sample, significantvalues);
0588     AbstractHistogram& new_histo = out[significantvalues];
0589     if (!new_histo.me) {
0590       auto name = makePathName(s, significantvalues, &step);
0591       iBooker.setCurrentFolder(name.first);
0592       if (dynamic_cast<TH1F*>(th1))
0593         new_histo.me = iBooker.book1D(name.second, (TH1F*)th1);
0594       else if (dynamic_cast<TH2F*>(th1))
0595         new_histo.me = iBooker.book2D(name.second, (TH2F*)th1);
0596       else if (dynamic_cast<TProfile*>(th1))
0597         new_histo.me = iBooker.bookProfile(name.second, (TProfile*)th1);
0598       else if (dynamic_cast<TProfile2D*>(th1))
0599         new_histo.me = iBooker.bookProfile2D(name.second, (TProfile2D*)th1);
0600       else
0601         assert(!"No idea how to book this.");
0602       new_histo.th1 = new_histo.me->getTH1();
0603       new_histo.iq_sample = e.second.iq_sample;
0604     } else {
0605       new_histo.th1->Add(th1);
0606     }
0607     new_histo.me->setStatOverflows(e.second.me->getStatOverflows());
0608   }
0609   t.swap(out);
0610 }
0611 
0612 void HistogramManager::executeExtend(SummationStep const& step,
0613                                      Table& t,
0614                                      std::string const& reduce_type,
0615                                      DQMStore::IBooker& iBooker,
0616                                      SummationSpecification const& s) {
0617   // For the moment only X.
0618   // first pass determines the range.
0619   std::map<GeometryInterface::Values, int> nbins;
0620   // separators collects meta info for the render plugin about the boundaries.
0621   // for each EXTEND, this is added to the axis label. In total this is not
0622   // fully correct since we have to assume the the substructure of each sub-
0623   // histogram is the same, which is e.g. not true for layers. It still works
0624   // since layers only contain leaves (ladders).
0625   std::map<GeometryInterface::Values, std::string> separators;
0626 
0627   GeometryInterface::Values significantvalues;
0628 
0629   for (auto& e : t) {
0630     geometryInterface.extractColumns(step.columns, e.second.iq_sample, significantvalues);
0631     int& n = nbins[significantvalues];
0632     assert(e.second.th1 || !"invalid histogram");
0633     // if we reduce, every histogram only needs one bin
0634     int bins = !reduce_type.empty() ? 1 : e.second.th1->GetXaxis()->GetNbins();
0635     if (bins > 1)
0636       separators[significantvalues] += std::to_string(n) + ",";
0637     n += bins;
0638   }
0639   for (auto& e : separators)
0640     e.second = "(" + e.second + ")/";
0641 
0642   Table out;
0643   for (auto& e : t) {
0644     geometryInterface.extractColumns(step.columns, e.second.iq_sample, significantvalues);
0645     TH1* th1 = e.second.th1;
0646     assert(th1);
0647 
0648     AbstractHistogram& new_histo = out[significantvalues];
0649     if (!new_histo.me) {
0650       // we put the name of the actual, last column of a input histo there.
0651       std::string colname = geometryInterface.pretty((e.first.end() - 1)->first);
0652 
0653       auto separator = separators[significantvalues];
0654 
0655       auto name = makePathName(s, significantvalues, &step);
0656       auto title =
0657           std::string("") + th1->GetTitle() + " per " + colname + ";" + colname + separator +
0658           (!reduce_type.empty() ? th1->GetYaxis()->GetTitle() : th1->GetXaxis()->GetTitle()) + ";" +
0659           (!reduce_type.empty() ? reduce_type + " of " + th1->GetXaxis()->GetTitle() : th1->GetYaxis()->GetTitle());
0660       iBooker.setCurrentFolder(name.first);
0661 
0662       if (th1->GetDimension() == 1) {
0663         new_histo.me =
0664             iBooker.book1D(name.second, title, nbins[significantvalues], 0.5, nbins[significantvalues] + 0.5);
0665       } else {
0666         assert(!"2D extend not implemented in harvesting.");
0667       }
0668       new_histo.th1 = new_histo.me->getTH1();
0669       new_histo.iq_sample = e.second.iq_sample;
0670       new_histo.count = 1;  // used as a fill pointer. Assumes histograms are
0671                             // ordered correctly (map should provide that)
0672     }
0673 
0674     // now add data.
0675     if (new_histo.th1->GetDimension() == 1) {
0676       if (reduce_type.empty()) {  // no reduction, concatenate
0677         for (int i = 1; i <= th1->GetXaxis()->GetNbins(); i++) {
0678           new_histo.th1->SetBinContent(new_histo.count, th1->GetBinContent(i));
0679           new_histo.th1->SetBinError(new_histo.count, th1->GetBinError(i));
0680           new_histo.count += 1;
0681         }
0682       } else if (reduce_type == "MEAN") {
0683         new_histo.th1->SetBinContent(new_histo.count, th1->GetMean());
0684         new_histo.th1->SetBinError(new_histo.count, th1->GetMeanError());
0685         new_histo.count += 1;
0686       } else {
0687         assert(!"Reduction type not supported");
0688       }
0689       new_histo.me->setStatOverflows(e.second.me->getStatOverflows());
0690     } else {
0691       assert(!"2D extend not implemented in harvesting.");
0692     }
0693   }
0694   t.swap(out);
0695 }
0696 
0697 void HistogramManager::executeHarvesting(DQMStore::IBooker& iBooker, DQMStore::IGetter& iGetter) {
0698   if (!enabled)
0699     return;
0700   // Debug output
0701   for (auto& s : specs) {
0702     edm::LogInfo log("HistogramManager");
0703     log << "Specs for " << name << " ";
0704     s.dump(log, geometryInterface);
0705   }
0706 
0707   for (unsigned int i = 0; i < specs.size(); i++) {
0708     auto& s = specs[i];
0709     auto& t = tables[i];
0710     loadFromDQMStore(s, t, iGetter);
0711 
0712     std::string reduce_type = "";
0713 
0714     // now execute step2.
0715     for (SummationStep const& step : s.steps) {
0716       if (step.stage == SummationStep::STAGE2) {
0717         switch (step.type) {
0718           case SummationStep::SAVE:
0719             // no explicit implementation atm.
0720             break;
0721           case SummationStep::GROUPBY:
0722             executeGroupBy(step, t, iBooker, s);
0723             break;
0724           case SummationStep::REDUCE:
0725             // reduction is done in the following EXTEND
0726             reduce_type = step.arg;
0727             break;
0728           case SummationStep::EXTEND_X:
0729             executeExtend(step, t, reduce_type, iBooker, s);
0730             reduce_type = "";
0731             break;
0732           case SummationStep::EXTEND_Y:
0733             assert(!"EXTEND_Y currently not supported in harvesting.");
0734             break;
0735           default:
0736             assert(!"Operation not supported in harvesting.");
0737         }
0738       }
0739     }
0740   }
0741 }