rivet is hosted by Hepforge, IPPP Durham
Rivet 4.0.0
AnalysisHandler.hh
1// -*- C++ -*-
2#ifndef RIVET_RivetHandler_HH
3#define RIVET_RivetHandler_HH
4
5#include "Rivet/Config/RivetCommon.hh"
6#include "Rivet/Particle.hh"
7#include "Rivet/AnalysisLoader.hh"
8#include "Rivet/Tools/RivetYODA.hh"
9#include "Rivet/Tools/Utils.hh"
10#include "Rivet/ProjectionHandler.hh"
11#include "YODA/ReaderYODA.h"
12
13#include <fstream>
14#include <unordered_map>
15
16namespace Rivet {
17
18
19 // Forward declaration and smart pointer for Analysis
20 class Analysis;
21 using AnaHandle = std::shared_ptr<Analysis>;
22
23
30
31 using TypeHandlePtr = std::shared_ptr<TypeBaseHandle>;
32 using TypeRegister = std::unordered_map<string, TypeHandlePtr>;
33 using TypeRegisterItr = typename TypeRegister::const_iterator;
34
35 public:
36
38 AnalysisHandler(const string& runname="");
39
42
45
48
49
52
54 string runName() const;
55
61 size_t numEvents() const {
62 const double N = _eventCounter.get()->persistent(defaultWeightIndex())->numEntries();
63 return size_t(N + 0.5 - (N<0)); // round to nearest integer
64 }
65
71 double effNumEvents() const {
72 if ((bool)_eventCounter) { return _eventCounter->effNumEntries(); }
73 return _eventCounter.get()->persistent(defaultWeightIndex())->effNumEntries();
74 }
75
80 double sumW() const {
81 if ((bool)_eventCounter) { return _eventCounter->sumW(); }
82 return _eventCounter.get()->persistent(defaultWeightIndex())->sumW();
83 }
85 double sumW2() const {
86 if ((bool)_eventCounter) { return _eventCounter->sumW2(); }
87 return _eventCounter.get()->persistent(defaultWeightIndex())->sumW2();
88 }
89
91
92
95
97 const vector<string>& weightNames() const { return _weightNames; }
98
100 size_t numWeights() const { return _weightNames.size(); }
101
103 bool haveNamedWeights() const;
104
106 void setWeightNames(const GenEvent& ge);
107
109 void setWeightNames(const vector<string>& weightNames);
110
112 size_t defaultWeightIndex() const { return _rivetDefaultWeightIdx; }
113
115 vector<double> weightSumWs() const;
116
118 void setWeightCap(const double maxWeight) { _weightCap = maxWeight; }
119
121 void setNominalWeightName(const std::string& name) { _nominalWeightName = name; }
122
124 void skipMultiWeights(bool skip=false) { _skipMultiWeights = skip; }
125
127 void matchWeightNames(const std::string& patterns) { _matchWeightNames = patterns; }
128
130 void unmatchWeightNames(const std::string& patterns) { _unmatchWeightNames = patterns; }
131
133 void setNLOSmearing(double frac) { _NLOSmearing = frac; }
134
136
137
140
142 Estimate0DPtr crossSection() const { return _xs; }
143
145 void setCrossSection(const vector<pair<double,double>>& xsecs, bool isUserSupplied = false);
146
148 void setCrossSection(const pair<double, double>& xsec, bool isUserSupplied=false);
149
151 void setCrossSection(double xsec, double xsecerr, bool isUserSupplied=false) {
152 setCrossSection({xsec, xsecerr}, isUserSupplied);
153 }
154
159
161 void notifyEndOfFile() { _isEndOfFile = true; }
162
164 double nominalCrossSection() const;
165
168
170
171
174
177
179 const ParticlePair& runBeams() const { return _beams; }
180
182 PdgIdPair runBeamIDs() const;
183
185 pair<double,double> runBeamEnergies() const;
186
188 double runSqrtS() const;
189
191 void setCheckBeams(bool check=true) { _checkBeams = check; }
192
194 // void setCheckConsistency(bool check=true) { _checkConsistency = check; }
195 // Check event consistency with the run, usually determined from the first event
196 // bool consistentWithRun(Event& event) {
197
199
200
203
205 template<typename T>
207 const std::string name = T().type();
208 const TypeRegisterItr& res = _register.find(name);
209 if (res == _register.end()) {
210 _register[name] = make_shared<TypeHandle<T>>();
211 }
212 _reader.registerType<T>(); // also let YODA know
213 }
214
217 bool copyAO(YODA::AnalysisObjectPtr src, YODA::AnalysisObjectPtr dst, const double scale=1.0);
218
221 bool addAO(YODA::AnalysisObjectPtr src, YODA::AnalysisObjectPtr& dst, const double scale);
222
224
225
228
230 std::vector<std::string> analysisNames() const;
231
233 std::vector<std::string> stdAnalysisNames() const;
234
236 const std::map<std::string, AnaHandle>& analysesMap() const {
237 return _analyses;
238 }
239
241 std::vector<AnaHandle> analyses() const {
242 std::vector<AnaHandle> rtn;
243 rtn.reserve(_analyses.size());
244 for (const auto& apair : _analyses) rtn.push_back(apair.second);
245 return rtn;
246 }
247
249 AnaHandle analysis(const std::string& analysisname) {
250 if ( _analyses.find(analysisname) == _analyses.end() )
251 throw LookupError("No analysis named '" + analysisname + "' registered in AnalysisHandler");
252 try {
253 return _analyses[analysisname];
254 } catch (...) {
255 throw LookupError("No analysis named '" + analysisname + "' registered in AnalysisHandler");
256 }
257 }
258
261
267 AnalysisHandler& addAnalysis(const std::string& analysisname);
268
270 AnalysisHandler& addAnalysis(const std::string& analysisname, std::map<string, string> pars);
271
278 AnalysisHandler& addAnalyses(const std::vector<std::string>& analysisnames);
279
280
282 AnalysisHandler& removeAnalysis(const std::string& analysisname);
283
285 AnalysisHandler& removeAnalyses(const std::vector<std::string>& analysisnames);
286
288
289
292
294 void init(const GenEvent& event);
295
302 void analyze(GenEvent& event);
303
308 void analyze(GenEvent* event);
309
312 void finalize();
313
315
316
319
323
327 void readData(std::istream& istr, const string& fmt, bool preload = true);
328
330 void readData(const std::string& filename, bool preload = true);
331
335 vector<YODA::AnalysisObjectPtr> getYodaAOs(const bool includeraw=false, const bool mkinert=true) const;
336
338 vector<YODA::AnalysisObjectPtr> getRawAOs() const;
339
341 vector<std::string> getRawAOPaths() const;
342
345 const YODA::AnalysisObjectPtr getPreload(const string& path) const {
346 auto it = _preloads.find(path);
347 if ( it == _preloads.end() ) return nullptr;
348 return it->second;
349 }
350
354 void writeData(std::ostream& ostr, const string& fmt) const;
355
357 void writeData(const string& filename) const;
358
364 void setFinalizePeriod(const string& dumpfile, int period) {
365 _dumpPeriod = period;
366 _dumpFile = dumpfile;
367 }
370 setFinalizePeriod("DUMMY", -1);
371 }
372
374 void setBootstrapFilename(const string& filename) {
375 _bootstrapfilename = filename;
376 }
377
379 vector<pair<string,size_t>> fillLayout() const;
380
382 vector<bool> fillOutcomes() const;
383
385 vector<double> fillFractions() const;
386
401
402 void mergeYODAs(const vector<string>& aofiles,
403 const vector<string>& delopts=vector<string>(),
404 const vector<string>& addopts=vector<string>(),
405 const vector<string>& matches=vector<string>(),
406 const vector<string>& unmatches=vector<string>(),
407 const bool equiv=false, const bool reentrantOnly = true);
408
410 void merge(AnalysisHandler &other);
411
417 void loadAOs(const vector<string>& aoPaths, const vector<double>& aoData);
418
420
423
424 vector<double> serializeContent(bool fixed_length = false) {
425 if (!_initialised)
426 throw Error("AnalysisHandler has not been initialised!");
427
429
430 // Loop over raw AOs and work out the size of the content data
431 const vector<YODA::AnalysisObjectPtr> raos = getRawAOs();
432 size_t total = 0;
433 for (size_t i = 0; i < raos.size(); ++i) {
434 total += raos[i]->lengthContent(fixed_length)+1;
435 }
436 total += _beaminfo->numBins()+1;
437
438 // Loop over raw AOs and retrieve the content data
439 std::vector<double> data; // serialized data vector
440 data.reserve(total); // pre-allocate enough memory
441 // Add beam IDs
442 data.push_back(_beaminfo->numBins());
443 for (const string& beamID : _beaminfo->xEdges()) {
444 data.push_back(_beamInfoLabelToID(beamID));
445 }
446 // Add raw YODA AO content
447 for (size_t i = 0; i < raos.size(); ++i) {
448 vector<double> tmp = raos[i]->serializeContent(fixed_length);
449 data.insert(std::end(data),
450 std::make_move_iterator(std::begin(tmp)),
451 std::make_move_iterator(std::end(tmp)));
452 }
453 return data;
454 }
455
456 void deserializeContent(const vector<double>& data, size_t nprocs = 0) {
457 if (!_initialised)
458 throw Error("AnalysisHandler has not been initialised!");
459
461
462 // get Rivet AOs for access to raw AO pointers
463 vector<MultiplexAOPtr> raos = getRivetAOs();
464
465
466 // beam info first
467 size_t iAO = 0, iW = 0, nBeams = data[0], offset = 1;
468 if (nprocs) nBeams /= nprocs;
469 const auto itr = data.cbegin();
470 // set beam IDs
471 vector<int> edges{itr+offset, itr+offset+nBeams};
472 vector<string> labels; labels.reserve(edges.size());
473 size_t id = 0;
474 for (int edge : edges) {
475 if (nprocs >= 2) edge /= nprocs;
476 labels.push_back(_mkBeamInfoLabel(++id, edge));
477 }
478
479 _beaminfo = make_shared<YODA::BinnedEstimate<string>>(labels, "/TMP/_BEAMPZ");
480 offset+= nBeams;
481 // set beam momenta
482 size_t beamLen = *(itr + offset); ++offset;
483 if (nprocs) beamLen /= nprocs;
484 std::vector<double> energies{itr+offset, itr+offset+beamLen};
485 if (nprocs >= 2) {
486 for (double& e : energies) { e /= nprocs; }
487 }
488 _beaminfo->deserializeContent(energies);
489 offset += beamLen;
490
491 // then the multiweighted AOs
492 while (offset < data.size()) {
493 if (iW < numWeights()) raos[iAO].get()->setActiveWeightIdx(iW);
494 else {
495 raos[iAO].get()->unsetActiveWeight();
496 iW = 0; ++iAO; // move on to next AO
497 raos[iAO].get()->setActiveWeightIdx(iW);
498 }
499
500 // obtain content length and set content iterators
501 size_t aoLen = *(itr + offset); ++offset;
502 if (nprocs) aoLen /= nprocs;
503 auto first = itr + offset;
504 auto last = first + aoLen;
505 // load data into AO
506 raos[iAO].get()->activeAO()->deserializeContent(std::vector<double>{first, last});
507
508 ++iW; offset += aoLen; // increment offset
509 }
510 raos[iAO].get()->unsetActiveWeight();
511 }
512
514
515
518
521 enum class Stage { OTHER, INIT, FINALIZE };
522
524 Stage stage() const { return _stage; }
525
527
528 private:
529
532
534 Log& getLog() const;
535
537 vector<MultiplexAOPtr> getRivetAOs() const;
538
540 void stripOptions(YODA::AnalysisObjectPtr ao, const vector<string>& delopts) const;
541
543 void mergeAOS(map<string, YODA::AnalysisObjectPtr> &allaos,
544 const map<string, YODA::AnalysisObjectPtr> &newaos,
545 map<string, pair<double, double>> &allxsecs,
546 const vector<string>& delopts=vector<string>(),
547 const vector<string>& optAnas=vector<string>(),
548 const vector<string>& optKeys=vector<string>(),
549 const vector<string>& optVals=vector<string>(),
550 const bool equiv=false,
551 const bool overwrite_xsec = false,
552 const double user_xsec = 1.0);
553
554
559 void loadAOs(const map<string, YODA::AnalysisObjectPtr>& allAOs,
560 const bool unscale = false, const bool reentrantOnly = true);
561
563 void _setRunBeamInfo(const ParticlePair& beams);
564
566 void _setRunBeamInfo(YODA::AnalysisObjectPtr ao);
567
569 string _mkBeamInfoLabel(size_t n, PdgId id) {
570 return "BEAM"+std::to_string(n)+"("+std::to_string(id)+")";
571 }
572
574 PdgId _beamInfoLabelToID(const string& label) {
575 size_t pos = label.find("(");
576 string beamID = label.substr(pos+1, label.size()-pos-2);
577 return std::stoi(beamID);
578 }
579
580
582
583
584 private:
585
587 Stage _stage = Stage::OTHER;
588
590 std::map<std::string, AnaHandle> _analyses;
591
595 map<string,YODA::AnalysisObjectPtr> _preloads;
596
598 vector<YODA::AnalysisObjectPtr> _finalizedAOs;
599
601 template<typename... Args>
602 void registerDefaultTypes();
603
605 TypeRegister _register;
606
608 YODA::Reader& _reader = YODA::ReaderYODA::create();
609
610
613
615 std::vector<std::string> _weightNames;
616 std::vector<std::valarray<double> > _subEventWeights;
617 //size_t _numWeightTypes; // always == WeightVector.size()
618
620 std::vector<size_t> _weightIndices;
621
623 std::string _runname;
624
626 CounterPtr _eventCounter;
627
629 Estimate0DPtr _xs;
630
632 YODA::BinnedEstimatePtr<string> _beaminfo;
633
635 vector<Estimate0D> _xsAvg;
636
638 double _numEntriesAggregate;
639
641 bool _isEndOfFile;
642
644 std::pair<double,double> _userxs;
645
647 ParticlePair _beams;
648
650 bool _initialised;
651
653 bool _checkBeams;
654
656 bool _skipMultiWeights;
657
659 std::string _matchWeightNames;
660
662 std::string _unmatchWeightNames;
663
665 std::string _nominalWeightName;
666
668 double _weightCap;
669
673 double _NLOSmearing;
674
676 int _eventNumber;
677
679 size_t _defaultWeightIdx;
680
682 size_t _rivetDefaultWeightIdx;
683
685 int _customDefaultWeightIdx;
686
688 int _dumpPeriod;
689
691 string _dumpFile;
692
694 bool _dumping;
695
697 ofstream _fbootstrap;
698
700 std::string _bootstrapfilename;
701
703 ProjectionHandler _projHandler;
704
706
707 };
708
709
710}
711
712#endif
The key class for coordination of Analysis objects and the event loop.
Definition AnalysisHandler.hh:29
Estimate0DPtr crossSection() const
Get the cross-section known to the handler.
Definition AnalysisHandler.hh:142
AnalysisHandler & addAnalysis(const std::string &analysisname, std::map< string, string > pars)
Add an analysis with a map of analysis options.
void setNLOSmearing(double frac)
Set the relative width of the NLO smearing window.
Definition AnalysisHandler.hh:133
void matchWeightNames(const std::string &patterns)
Specify weight-name patterns to accept.
Definition AnalysisHandler.hh:127
std::vector< std::string > stdAnalysisNames() const
Get a list of the official analysis names for this release.
vector< bool > fillOutcomes() const
Return a vector of the binary fill outcome (was/wasn't filled) at each fill position.
size_t defaultWeightIndex() const
Get the index of the nominal weight-stream.
Definition AnalysisHandler.hh:112
void mergeYODAs(const vector< string > &aofiles, const vector< string > &delopts=vector< string >(), const vector< string > &addopts=vector< string >(), const vector< string > &matches=vector< string >(), const vector< string > &unmatches=vector< string >(), const bool equiv=false, const bool reentrantOnly=true)
Merge the vector of YODA files, using the cross-section and weight information provided in each.
void setNominalWeightName(const std::string &name)
Set the name of the nominal weight stream.
Definition AnalysisHandler.hh:121
void writeData(const string &filename) const
Write all analyses' plots (via getData) to the named file.
AnalysisHandler & setRunBeams(const ParticlePair &beams)
Set the beam particles for this run.
void setCrossSection(double xsec, double xsecerr, bool isUserSupplied=false)
Set the cross-section for the process being generated (alternative signature)
Definition AnalysisHandler.hh:151
void init(const GenEvent &event)
Initialize a run, with the run beams taken from the example event.
AnalysisHandler & addAnalyses(const std::vector< std::string > &analysisnames)
Add analyses to the run list using their names.
AnalysisHandler & removeAnalyses(const std::vector< std::string > &analysisnames)
Remove analyses from the run list using their names.
size_t numWeights() const
Are any of the weights non-numeric?
Definition AnalysisHandler.hh:100
void analyze(GenEvent &event)
Analyze the given event by reference.
vector< YODA::AnalysisObjectPtr > getYodaAOs(const bool includeraw=false, const bool mkinert=true) const
vector< double > fillFractions() const
Return a vector of the fill fraction at each fill position.
void setCrossSection(const pair< double, double > &xsec, bool isUserSupplied=false)
Set all cross-sections for the process being generated, based on nominal weight.
string runName() const
Get the name of this run.
void registerType()
Register an AO type handle into type map and YODA reader.
Definition AnalysisHandler.hh:206
AnalysisHandler(const string &runname="")
Preferred constructor, with optional run name.
pair< double, double > runBeamEnergies() const
Get beam IDs for this run, usually determined from the first event.
const YODA::AnalysisObjectPtr getPreload(const string &path) const
Definition AnalysisHandler.hh:345
void analyze(GenEvent *event)
Analyze the given event by pointer.
~AnalysisHandler()
The destructor is not virtual, as this class should not be inherited from.
void merge(AnalysisHandler &other)
A method to merge another AnalysisHandler into the current one.
void readData(const std::string &filename, bool preload=true)
Read analysis plots into the histo collection (via addData) from the named file.
bool copyAO(YODA::AnalysisObjectPtr src, YODA::AnalysisObjectPtr dst, const double scale=1.0)
const vector< string > & weightNames() const
Names of event weight categories.
Definition AnalysisHandler.hh:97
void setFinalizePeriod(const string &dumpfile, int period)
Configure the AnalysisObject dump rate and destination.
Definition AnalysisHandler.hh:364
void setCrossSection(const vector< pair< double, double > > &xsecs, bool isUserSupplied=false)
Set all cross-sections for the process being generated specifically (preferred)
double nominalCrossSectionError() const
Get the nominal cross-section.
void readData(std::istream &istr, const string &fmt, bool preload=true)
Read analysis plots into the histo collection from the given stream.
void notifyEndOfFile()
Toggle to signal a change in HepMC input file.
Definition AnalysisHandler.hh:161
Stage stage() const
Return the current processing stage.
Definition AnalysisHandler.hh:524
double sumW2() const
Access to the sum of squared-weights.
Definition AnalysisHandler.hh:85
AnalysisHandler & addAnalysis(Analysis *analysis)
Add an analysis to the run list by object.
PdgIdPair runBeamIDs() const
Get beam IDs for this run, usually determined from the first event.
bool addAO(YODA::AnalysisObjectPtr src, YODA::AnalysisObjectPtr &dst, const double scale)
double runSqrtS() const
Get energy for this run, usually determined from the first event.
bool haveNamedWeights() const
Are any of the weights non-numeric?
void setCheckBeams(bool check=true)
Option to disable analysis-compatibility checks.
Definition AnalysisHandler.hh:191
void setWeightNames(const GenEvent &ge)
Set the weight names from a GenEvent.
size_t numEvents() const
Definition AnalysisHandler.hh:61
void loadAOs(const vector< string > &aoPaths, const vector< double > &aoData)
A method to prepare a re-entrant run for a given set of AO paths and serialized AO data.
vector< double > weightSumWs() const
Access the array of sum of the event weights seen.
std::vector< AnaHandle > analyses() const
Get the collection of currently registered analyses.
Definition AnalysisHandler.hh:241
void setWeightCap(const double maxWeight)
Set the weight cap.
Definition AnalysisHandler.hh:118
AnaHandle analysis(const std::string &analysisname)
Get a registered analysis by name.
Definition AnalysisHandler.hh:249
void writeData(std::ostream &ostr, const string &fmt) const
Write all analyses' plots (via getData) to the given stream.
const std::map< std::string, AnaHandle > & analysesMap() const
Get the collection of currently registered analyses.
Definition AnalysisHandler.hh:236
AnalysisHandler & addAnalysis(const std::string &analysisname)
Add an analysis to the run list using its name.
Stage
Definition AnalysisHandler.hh:521
AnalysisHandler & removeAnalysis(const std::string &analysisname)
Remove an analysis from the run list using its name.
vector< pair< string, size_t > > fillLayout() const
Return a vector of (AO path, AO numBins) pairs to decode the fills layout.
AnalysisHandler(const AnalysisHandler &)=delete
The copy constructor is deleted, so it can never be called.
void skipMultiWeights(bool skip=false)
Ignore all weight streams other than the nominal.
Definition AnalysisHandler.hh:124
const ParticlePair & runBeams() const
Get the beam particles for this run, usually determined from the first event.
Definition AnalysisHandler.hh:179
vector< std::string > getRawAOPaths() const
Get all raw YODA analysis object paths (across all weights)
double effNumEvents() const
Definition AnalysisHandler.hh:71
vector< YODA::AnalysisObjectPtr > getRawAOs() const
Get all raw YODA analysis objects (across all weights)
void unmatchWeightNames(const std::string &patterns)
Specify weight-name patterns to reject.
Definition AnalysisHandler.hh:130
void setWeightNames(const vector< string > &weightNames)
Set the weight names from a vector<string>
void setBootstrapFilename(const string &filename)
Set filename of the bootstrap file.
Definition AnalysisHandler.hh:374
double nominalCrossSection() const
Get the nominal cross-section.
void setNoFinalizePeriod()
Configure the AnalysisObject dump rate and destination.
Definition AnalysisHandler.hh:369
double sumW() const
Access the sum of the event weights seen.
Definition AnalysisHandler.hh:80
AnalysisHandler & operator=(const AnalysisHandler &)=delete
The assignment operator is deleted, so it can never be called.
std::vector< std::string > analysisNames() const
Get a list of the currently registered analyses' names.
This is the base class of all analysis classes in Rivet.
Definition Analysis.hh:67
Logging system for controlled & formatted writing to stdout.
Definition Logging.hh:10
Definition RivetYODA.hh:1330
ParticlePair beams(const Event &e)
Get beam particles from an event.
Definition MC_CENT_PPB_Projections.hh:10
std::pair< Particle, Particle > ParticlePair
Typedef for a pair of Particle objects.
Definition Particle.hh:38
Generic runtime Rivet error.
Definition Exceptions.hh:12
Error relating to looking up analysis objects in the register.
Definition Exceptions.hh:61