File indexing completed on 2026-05-06 08:52:10
0001
0002 #ifndef RIVET_RivetHandler_HH
0003 #define RIVET_RivetHandler_HH
0004
0005 #include "Rivet/Config/RivetCommon.hh"
0006 #include "Rivet/Particle.hh"
0007 #include "Rivet/AnalysisLoader.hh"
0008 #include "Rivet/Tools/RivetYODA.hh"
0009 #include "Rivet/Tools/Utils.hh"
0010 #include "Rivet/ProjectionHandler.hh"
0011 #include "YODA/ReaderYODA.h"
0012
0013 #include <fstream>
0014 #include <unordered_map>
0015
0016 namespace Rivet {
0017
0018
0019
0020 class Analysis;
0021 using AnaHandle = std::shared_ptr<Analysis>;
0022
0023
0024
0025
0026
0027
0028
0029 class AnalysisHandler {
0030
0031 using TypeHandlePtr = std::shared_ptr<TypeBaseHandle>;
0032 using TypeRegister = std::unordered_map<string, TypeHandlePtr>;
0033 using TypeRegisterItr = typename TypeRegister::const_iterator;
0034
0035 public:
0036
0037 using Annotations = std::map<std::string, std::string>;
0038
0039
0040 AnalysisHandler();
0041
0042
0043 AnalysisHandler(const AnalysisHandler&) = delete;
0044
0045
0046 AnalysisHandler& operator=(const AnalysisHandler&) = delete;
0047
0048
0049 ~AnalysisHandler();
0050
0051
0052
0053
0054
0055
0056
0057
0058
0059
0060 size_t numEvents() const {
0061 const double N = _eventCounter.get()->persistent(defaultWeightIndex())->numEntries();
0062 return size_t(N + 0.5 - (N<0));
0063 }
0064
0065
0066
0067
0068
0069
0070 double effNumEvents() const {
0071 if ((bool)_eventCounter) { return _eventCounter->effNumEntries(); }
0072 return _eventCounter.get()->persistent(defaultWeightIndex())->effNumEntries();
0073 }
0074
0075
0076
0077
0078
0079 double sumW() const {
0080 if ((bool)_eventCounter) { return _eventCounter->sumW(); }
0081 return _eventCounter.get()->persistent(defaultWeightIndex())->sumW();
0082 }
0083
0084 double sumW2() const {
0085 if ((bool)_eventCounter) { return _eventCounter->sumW2(); }
0086 return _eventCounter.get()->persistent(defaultWeightIndex())->sumW2();
0087 }
0088
0089
0090
0091
0092
0093
0094
0095
0096 const vector<string>& weightNames() const { return _weightNames; }
0097
0098
0099 size_t numWeights() const { return _weightNames.size(); }
0100
0101
0102 bool haveNamedWeights() const;
0103
0104
0105 void setWeightNames(const GenEvent& ge);
0106
0107
0108 void setWeightNames(const vector<string>& weightNames);
0109
0110
0111 size_t defaultWeightIndex() const { return _rivetDefaultWeightIdx; }
0112
0113
0114 vector<double> weightSumWs() const;
0115
0116
0117 void setWeightCap(const double maxWeight) { _weightCap = maxWeight; }
0118
0119
0120 void setNominalWeightName(const std::string& name) { _nominalWeightName = name; }
0121
0122
0123 void skipMultiWeights(bool skip=false) { _skipMultiWeights = skip; }
0124
0125
0126 void matchWeightNames(const std::string& patterns) { _matchWeightNames = patterns; }
0127
0128
0129 void unmatchWeightNames(const std::string& patterns) { _unmatchWeightNames = patterns; }
0130
0131
0132 void setNLOSmearing(double frac) { _NLOSmearing = frac; }
0133
0134
0135
0136
0137
0138
0139
0140
0141 Estimate0DPtr crossSection() const { return _xs; }
0142
0143
0144 void setCrossSection(const vector<pair<double,double>>& xsecs, bool isUserSupplied = false);
0145
0146
0147 void setCrossSection(const pair<double, double>& xsec, bool isUserSupplied=false);
0148
0149
0150 void setCrossSection(double xsec, double xsecerr, bool isUserSupplied=false) {
0151 setCrossSection({xsec, xsecerr}, isUserSupplied);
0152 }
0153
0154
0155
0156
0157 void updateCrossSection();
0158
0159
0160 void notifyEndOfFile() { _isEndOfFile = true; }
0161
0162
0163 double nominalCrossSection() const;
0164
0165
0166 double nominalCrossSectionError() const;
0167
0168
0169
0170
0171
0172
0173
0174
0175 AnalysisHandler& setRunBeams(const ParticlePair& beams);
0176
0177
0178 const ParticlePair& runBeams() const { return _beams; }
0179
0180
0181 PdgIdPair runBeamIDs() const;
0182
0183
0184 pair<double,double> runBeamEnergies() const;
0185
0186
0187 double runSqrtS() const;
0188
0189
0190 void setCheckBeams(bool check=true) { _checkBeams = check; }
0191
0192
0193
0194
0195
0196
0197
0198
0199
0200
0201
0202
0203
0204 std::vector<std::string> annotations() const {
0205 return _beaminfo->annotations();
0206 }
0207
0208
0209 bool hasAnnotation(const std::string& name) const {
0210 return _beaminfo->hasAnnotation(name);
0211 }
0212
0213
0214 const std::string& annotation(const std::string& name) const {
0215 return _beaminfo->annotation(name);
0216 }
0217
0218
0219 const std::string& annotation(const std::string& name, const std::string& defaultreturn) const {
0220 return _beaminfo->annotation(name, defaultreturn);
0221 }
0222
0223
0224
0225
0226 template <typename T>
0227 const T annotation(const std::string& name) const {
0228 return _beaminfo->annotation<T>(name);
0229 }
0230
0231
0232
0233
0234 template <typename T>
0235 const T annotation(const std::string& name, T&& defaultreturn) const {
0236 return _beaminfo->annotation<T>(name, std::forward<T>(defaultreturn));
0237 }
0238
0239
0240
0241
0242 template <typename T>
0243 void setAnnotation(const std::string& name, T&& value) {
0244 _beaminfo->setAnnotation<T>(name, std::forward<T>(value));
0245 }
0246
0247
0248
0249 void setAnnotations(const Annotations& anns) {
0250 _beaminfo->setAnnotations(anns);
0251 }
0252
0253
0254 void rmAnnotation(const std::string& name) {
0255 _beaminfo->rmAnnotation(name);
0256 }
0257
0258
0259
0260 void clearAnnotations() {
0261 _beaminfo->clearAnnotations();
0262 }
0263
0264
0265
0266
0267
0268
0269
0270
0271 template<typename T>
0272 void registerType() {
0273 const std::string name = T().type();
0274 const TypeRegisterItr& res = _register.find(name);
0275 if (res == _register.end()) {
0276 _register[name] = make_shared<TypeHandle<T>>();
0277 }
0278 _reader.registerType<T>();
0279 }
0280
0281
0282
0283 bool copyAO(YODA::AnalysisObjectPtr src, YODA::AnalysisObjectPtr dst, const double scale=1.0);
0284
0285
0286
0287 bool addAO(YODA::AnalysisObjectPtr src, YODA::AnalysisObjectPtr& dst, const double scale);
0288
0289
0290
0291
0292
0293
0294
0295
0296 std::vector<std::string> analysisNames() const;
0297
0298
0299 std::vector<std::string> stdAnalysisNames() const;
0300
0301
0302 const std::map<std::string, AnaHandle>& analysesMap() const {
0303 return _analyses;
0304 }
0305
0306
0307 std::vector<AnaHandle> analyses() const {
0308 std::vector<AnaHandle> rtn;
0309 rtn.reserve(_analyses.size());
0310 for (const auto& apair : _analyses) rtn.push_back(apair.second);
0311 return rtn;
0312 }
0313
0314
0315 AnaHandle analysis(const std::string& analysisname) {
0316 if ( _analyses.find(analysisname) == _analyses.end() )
0317 throw LookupError("No analysis named '" + analysisname + "' registered in AnalysisHandler");
0318 try {
0319 return _analyses[analysisname];
0320 } catch (...) {
0321 throw LookupError("No analysis named '" + analysisname + "' registered in AnalysisHandler");
0322 }
0323 }
0324
0325
0326 AnalysisHandler& addAnalysis(Analysis* analysis);
0327
0328
0329
0330
0331
0332
0333 AnalysisHandler& addAnalysis(const std::string& analysisname);
0334
0335
0336 AnalysisHandler& addAnalysis(const std::string& analysisname, std::map<string, string> pars);
0337
0338
0339
0340
0341
0342
0343
0344 AnalysisHandler& addAnalyses(const std::vector<std::string>& analysisnames);
0345
0346
0347
0348 AnalysisHandler& removeAnalysis(const std::string& analysisname);
0349
0350
0351 AnalysisHandler& removeAnalyses(const std::vector<std::string>& analysisnames);
0352
0353
0354
0355
0356
0357
0358
0359
0360 void init(const GenEvent& event);
0361
0362
0363
0364
0365
0366
0367
0368 void analyze(GenEvent& event);
0369
0370
0371
0372
0373
0374 void analyze(GenEvent* event);
0375
0376
0377
0378 void finalize();
0379
0380
0381
0382
0383
0384
0385
0386
0387
0388 void collapseEventGroup();
0389
0390
0391
0392
0393 void readData(std::istream& istr, const string& fmt, bool preload = true);
0394
0395
0396 void readData(const std::string& filename, bool preload = true);
0397
0398
0399
0400
0401 vector<YODA::AnalysisObjectPtr> getYodaAOs(const bool includeraw=false, const bool mkinert=true) const;
0402
0403
0404 vector<YODA::AnalysisObjectPtr> getRawAOs() const;
0405
0406
0407 vector<std::string> getRawAOPaths() const;
0408
0409
0410
0411 const YODA::AnalysisObjectPtr getPreload(const string& path) const {
0412 auto it = _preloads.find(path);
0413 if ( it == _preloads.end() ) return nullptr;
0414 return it->second;
0415 }
0416
0417
0418
0419
0420 void writeData(std::ostream& ostr, const string& fmt) const;
0421
0422
0423 void writeData(const string& filename) const;
0424
0425
0426
0427
0428
0429
0430 void setFinalizePeriod(const string& dumpfile, int period) {
0431 _dumpPeriod = period;
0432 _dumpFile = dumpfile;
0433 }
0434
0435 void setNoFinalizePeriod() {
0436 setFinalizePeriod("DUMMY", -1);
0437 }
0438
0439
0440 void setBootstrapFilename(const string& filename) {
0441 _bootstrapfilename = filename;
0442 }
0443
0444
0445 vector<pair<string,size_t>> fillLayout() const;
0446
0447
0448 vector<bool> fillOutcomes() const;
0449
0450
0451 vector<double> fillFractions() const;
0452
0453
0454
0455
0456
0457
0458
0459
0460
0461
0462
0463
0464
0465
0466
0467
0468 void mergeYODAs(const vector<string>& aofiles,
0469 const vector<string>& delopts=vector<string>(),
0470 const vector<string>& addopts=vector<string>(),
0471 const vector<string>& matches=vector<string>(),
0472 const vector<string>& unmatches=vector<string>(),
0473 const bool equiv=false,
0474 const bool reentrantOnly = true,
0475 const bool checkbeams = true);
0476
0477
0478 void merge(AnalysisHandler &other);
0479
0480
0481
0482
0483
0484
0485 void loadAOs(const vector<string>& aoPaths, const vector<double>& aoData);
0486
0487
0488
0489
0490
0491
0492 vector<double> serializeContent(bool fixed_length = false) {
0493 if (!_initialised)
0494 throw Error("AnalysisHandler has not been initialised!");
0495
0496 collapseEventGroup();
0497
0498
0499 const vector<YODA::AnalysisObjectPtr> raos = getRawAOs();
0500 size_t total = 0;
0501 for (size_t i = 0; i < raos.size(); ++i) {
0502 total += raos[i]->lengthContent(fixed_length)+1;
0503 }
0504 total += _beaminfo->numBins()+1;
0505
0506
0507 std::vector<double> data;
0508 data.reserve(total);
0509
0510 data.push_back(_beaminfo->numBins());
0511 for (const string& beamID : _beaminfo->xEdges()) {
0512 data.push_back(_beamInfoLabelToID(beamID));
0513 }
0514
0515 for (size_t i = 0; i < raos.size(); ++i) {
0516 vector<double> tmp = raos[i]->serializeContent(fixed_length);
0517 data.push_back(tmp.size());
0518 data.insert(std::end(data),
0519 std::make_move_iterator(std::begin(tmp)),
0520 std::make_move_iterator(std::end(tmp)));
0521 }
0522 return data;
0523 }
0524
0525 void deserializeContent(const vector<double>& data, size_t nprocs = 0) {
0526 if (!_initialised)
0527 throw Error("AnalysisHandler has not been initialised!");
0528
0529 collapseEventGroup();
0530
0531
0532 vector<MultiplexAOPtr> raos = getRivetAOs();
0533
0534
0535
0536 size_t iAO = 0, iW = 0, nBeams = data[0], offset = 1;
0537 if (nprocs) nBeams /= nprocs;
0538 const auto itr = data.cbegin();
0539
0540 vector<int> edges{itr+offset, itr+offset+nBeams};
0541 vector<string> labels; labels.reserve(edges.size());
0542 size_t id = 0;
0543 for (int edge : edges) {
0544 if (nprocs >= 2) edge /= nprocs;
0545 labels.push_back(_mkBeamInfoLabel(++id, edge));
0546 }
0547
0548 _beaminfo = make_shared<YODA::BinnedEstimate<string>>(labels, "/TMP/_BEAMPZ");
0549 offset += nBeams;
0550
0551 size_t beamLen = *(itr + offset); ++offset;
0552 if (nprocs) beamLen /= nprocs;
0553 std::vector<double> energies{itr+offset, itr+offset+beamLen};
0554 if (nprocs >= 2) {
0555 for (double& e : energies) { e /= nprocs; }
0556 }
0557 _beaminfo->deserializeContent(energies);
0558 for (auto& b : _beaminfo->bins(true)) b.rmErrs();
0559 offset += beamLen;
0560
0561
0562 while (offset < data.size()) {
0563 if (iW < numWeights()) raos[iAO].get()->setActiveWeightIdx(iW);
0564 else {
0565 raos[iAO].get()->unsetActiveWeight();
0566 iW = 0; ++iAO;
0567 raos[iAO].get()->setActiveWeightIdx(iW);
0568 }
0569
0570
0571 size_t aoLen = *(itr + offset); ++offset;
0572 if (nprocs) aoLen /= nprocs;
0573 auto first = itr + offset;
0574 auto last = first + aoLen;
0575
0576 raos[iAO].get()->activeAO()->deserializeContent(std::vector<double>{first, last});
0577
0578 ++iW; offset += aoLen;
0579 }
0580 raos[iAO].get()->unsetActiveWeight();
0581
0582 _ntrials = 0.0;
0583 _fileCounter = CounterPtr(weightNames(), Counter("_FILECOUNT"));
0584 _xserr = CounterPtr(weightNames(), Counter("XSECERR"));
0585 if (nprocs >= 2) {
0586 for (size_t iW = 0; iW < numWeights(); ++iW) {
0587 *_fileCounter.get()->persistent(iW) = *_eventCounter.get()->persistent(iW);
0588 _xs.get()->persistent(iW)->scale(1.0/nprocs);
0589 }
0590 }
0591 }
0592
0593
0594
0595
0596
0597
0598
0599
0600
0601 enum class Stage { OTHER, INIT, FINALIZE };
0602
0603
0604 Stage stage() const { return _stage; }
0605
0606
0607
0608 private:
0609
0610
0611
0612
0613
0614 Log& getLog() const;
0615
0616
0617 vector<MultiplexAOPtr> getRivetAOs() const;
0618
0619
0620 void stripOptions(YODA::AnalysisObjectPtr ao, const vector<string>& delopts) const;
0621
0622
0623 void mergeAOS(map<string, YODA::AnalysisObjectPtr> &allaos,
0624 const map<string, YODA::AnalysisObjectPtr> &newaos,
0625 map<string, std::array<double,4>> &allxsecs,
0626 const vector<string>& delopts=vector<string>(),
0627 const vector<string>& optAnas=vector<string>(),
0628 const vector<string>& optKeys=vector<string>(),
0629 const vector<string>& optVals=vector<string>(),
0630 const bool equiv=false,
0631 const bool overwrite_xsec = false,
0632 const double user_xsec = 1.0);
0633
0634
0635
0636
0637
0638
0639 void loadAOs(const map<string, YODA::AnalysisObjectPtr>& allAOs,
0640 const bool unscale = false, const bool reentrantOnly = true);
0641
0642
0643 void _setRunBeamInfo(const ParticlePair& beams);
0644
0645
0646 void _setRunBeamInfo(YODA::AnalysisObjectPtr ao);
0647
0648
0649 string _mkBeamInfoLabel(size_t n, PdgId id) {
0650 return "BEAM"+std::to_string(n)+"("+std::to_string(id)+")";
0651 }
0652
0653
0654 PdgId _beamInfoLabelToID(const string& label) {
0655 size_t pos = label.find("(");
0656 string beamID = label.substr(pos+1, label.size()-pos-2);
0657 return std::stoi(beamID);
0658 }
0659
0660
0661
0662
0663
0664 private:
0665
0666
0667 Stage _stage = Stage::OTHER;
0668
0669
0670 std::map<std::string, AnaHandle> _analyses;
0671
0672
0673
0674
0675 map<string,YODA::AnalysisObjectPtr> _preloads;
0676
0677
0678 vector<YODA::AnalysisObjectPtr> _finalizedAOs;
0679
0680
0681 template<typename... Args>
0682 void registerDefaultTypes();
0683
0684
0685 TypeRegister _register;
0686
0687
0688 YODA::Reader& _reader = YODA::ReaderYODA::create();
0689
0690
0691
0692
0693
0694
0695 std::vector<std::string> _weightNames;
0696 std::vector<std::valarray<double> > _subEventWeights;
0697
0698
0699
0700 std::vector<size_t> _weightIndices;
0701
0702
0703 CounterPtr _eventCounter;
0704
0705
0706 Estimate0DPtr _xs;
0707
0708
0709 CounterPtr _xserr;
0710
0711
0712 YODA::BinnedEstimatePtr<string> _beaminfo;
0713
0714
0715 double _ntrials;
0716
0717
0718 CounterPtr _fileCounter;
0719
0720
0721 bool _isEndOfFile;
0722
0723
0724 bool _subeventWarning;
0725
0726
0727 std::pair<double,double> _userxs;
0728
0729
0730 ParticlePair _beams;
0731
0732
0733 bool _initialised;
0734
0735
0736 bool _checkBeams;
0737
0738
0739 bool _skipMultiWeights;
0740
0741
0742 std::string _matchWeightNames;
0743
0744
0745 std::string _unmatchWeightNames;
0746
0747
0748 std::string _nominalWeightName;
0749
0750
0751 double _weightCap;
0752
0753
0754
0755
0756 double _NLOSmearing;
0757
0758
0759 int _eventNumber;
0760
0761
0762 size_t _defaultWeightIdx;
0763
0764
0765 size_t _rivetDefaultWeightIdx;
0766
0767
0768 int _customDefaultWeightIdx;
0769
0770
0771 int _dumpPeriod;
0772
0773
0774 string _dumpFile;
0775
0776
0777 bool _dumping;
0778
0779
0780 ofstream _fbootstrap;
0781
0782
0783 std::string _bootstrapfilename;
0784
0785
0786 ProjectionHandler _projHandler;
0787
0788
0789
0790 };
0791
0792
0793 }
0794
0795 #endif