File indexing completed on 2025-09-17 09:13:44
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->annotation<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, const bool reentrantOnly = true);
0474
0475
0476 void merge(AnalysisHandler &other);
0477
0478
0479
0480
0481
0482
0483 void loadAOs(const vector<string>& aoPaths, const vector<double>& aoData);
0484
0485
0486
0487
0488
0489
0490 vector<double> serializeContent(bool fixed_length = false) {
0491 if (!_initialised)
0492 throw Error("AnalysisHandler has not been initialised!");
0493
0494 collapseEventGroup();
0495
0496
0497 const vector<YODA::AnalysisObjectPtr> raos = getRawAOs();
0498 size_t total = 0;
0499 for (size_t i = 0; i < raos.size(); ++i) {
0500 total += raos[i]->lengthContent(fixed_length)+1;
0501 }
0502 total += _beaminfo->numBins()+1;
0503
0504
0505 std::vector<double> data;
0506 data.reserve(total);
0507
0508 data.push_back(_beaminfo->numBins());
0509 for (const string& beamID : _beaminfo->xEdges()) {
0510 data.push_back(_beamInfoLabelToID(beamID));
0511 }
0512
0513 for (size_t i = 0; i < raos.size(); ++i) {
0514 vector<double> tmp = raos[i]->serializeContent(fixed_length);
0515 data.push_back(tmp.size());
0516 data.insert(std::end(data),
0517 std::make_move_iterator(std::begin(tmp)),
0518 std::make_move_iterator(std::end(tmp)));
0519 }
0520 return data;
0521 }
0522
0523 void deserializeContent(const vector<double>& data, size_t nprocs = 0) {
0524 if (!_initialised)
0525 throw Error("AnalysisHandler has not been initialised!");
0526
0527 collapseEventGroup();
0528
0529
0530 vector<MultiplexAOPtr> raos = getRivetAOs();
0531
0532
0533
0534 size_t iAO = 0, iW = 0, nBeams = data[0], offset = 1;
0535 if (nprocs) nBeams /= nprocs;
0536 const auto itr = data.cbegin();
0537
0538 vector<int> edges{itr+offset, itr+offset+nBeams};
0539 vector<string> labels; labels.reserve(edges.size());
0540 size_t id = 0;
0541 for (int edge : edges) {
0542 if (nprocs >= 2) edge /= nprocs;
0543 labels.push_back(_mkBeamInfoLabel(++id, edge));
0544 }
0545
0546 _beaminfo = make_shared<YODA::BinnedEstimate<string>>(labels, "/TMP/_BEAMPZ");
0547 offset += nBeams;
0548
0549 size_t beamLen = *(itr + offset); ++offset;
0550 if (nprocs) beamLen /= nprocs;
0551 std::vector<double> energies{itr+offset, itr+offset+beamLen};
0552 if (nprocs >= 2) {
0553 for (double& e : energies) { e /= nprocs; }
0554 }
0555 _beaminfo->deserializeContent(energies);
0556 for (auto& b : _beaminfo->bins(true)) b.rmErrs();
0557 offset += beamLen;
0558
0559
0560 while (offset < data.size()) {
0561 if (iW < numWeights()) raos[iAO].get()->setActiveWeightIdx(iW);
0562 else {
0563 raos[iAO].get()->unsetActiveWeight();
0564 iW = 0; ++iAO;
0565 raos[iAO].get()->setActiveWeightIdx(iW);
0566 }
0567
0568
0569 size_t aoLen = *(itr + offset); ++offset;
0570 if (nprocs) aoLen /= nprocs;
0571 auto first = itr + offset;
0572 auto last = first + aoLen;
0573
0574 raos[iAO].get()->activeAO()->deserializeContent(std::vector<double>{first, last});
0575
0576 ++iW; offset += aoLen;
0577 }
0578 raos[iAO].get()->unsetActiveWeight();
0579
0580 _ntrials = 0.0;
0581 _fileCounter = CounterPtr(weightNames(), Counter("_FILECOUNT"));
0582 _xserr = CounterPtr(weightNames(), Counter("XSECERR"));
0583 if (nprocs >= 2) {
0584 for (size_t iW = 0; iW < numWeights(); ++iW) {
0585 *_fileCounter.get()->persistent(iW) = *_eventCounter.get()->persistent(iW);
0586 _xs.get()->persistent(iW)->scale(1.0/nprocs);
0587 }
0588 }
0589 }
0590
0591
0592
0593
0594
0595
0596
0597
0598
0599 enum class Stage { OTHER, INIT, FINALIZE };
0600
0601
0602 Stage stage() const { return _stage; }
0603
0604
0605
0606 private:
0607
0608
0609
0610
0611
0612 Log& getLog() const;
0613
0614
0615 vector<MultiplexAOPtr> getRivetAOs() const;
0616
0617
0618 void stripOptions(YODA::AnalysisObjectPtr ao, const vector<string>& delopts) const;
0619
0620
0621 void mergeAOS(map<string, YODA::AnalysisObjectPtr> &allaos,
0622 const map<string, YODA::AnalysisObjectPtr> &newaos,
0623 map<string, std::array<double,4>> &allxsecs,
0624 const vector<string>& delopts=vector<string>(),
0625 const vector<string>& optAnas=vector<string>(),
0626 const vector<string>& optKeys=vector<string>(),
0627 const vector<string>& optVals=vector<string>(),
0628 const bool equiv=false,
0629 const bool overwrite_xsec = false,
0630 const double user_xsec = 1.0);
0631
0632
0633
0634
0635
0636
0637 void loadAOs(const map<string, YODA::AnalysisObjectPtr>& allAOs,
0638 const bool unscale = false, const bool reentrantOnly = true);
0639
0640
0641 void _setRunBeamInfo(const ParticlePair& beams);
0642
0643
0644 void _setRunBeamInfo(YODA::AnalysisObjectPtr ao);
0645
0646
0647 string _mkBeamInfoLabel(size_t n, PdgId id) {
0648 return "BEAM"+std::to_string(n)+"("+std::to_string(id)+")";
0649 }
0650
0651
0652 PdgId _beamInfoLabelToID(const string& label) {
0653 size_t pos = label.find("(");
0654 string beamID = label.substr(pos+1, label.size()-pos-2);
0655 return std::stoi(beamID);
0656 }
0657
0658
0659
0660
0661
0662 private:
0663
0664
0665 Stage _stage = Stage::OTHER;
0666
0667
0668 std::map<std::string, AnaHandle> _analyses;
0669
0670
0671
0672
0673 map<string,YODA::AnalysisObjectPtr> _preloads;
0674
0675
0676 vector<YODA::AnalysisObjectPtr> _finalizedAOs;
0677
0678
0679 template<typename... Args>
0680 void registerDefaultTypes();
0681
0682
0683 TypeRegister _register;
0684
0685
0686 YODA::Reader& _reader = YODA::ReaderYODA::create();
0687
0688
0689
0690
0691
0692
0693 std::vector<std::string> _weightNames;
0694 std::vector<std::valarray<double> > _subEventWeights;
0695
0696
0697
0698 std::vector<size_t> _weightIndices;
0699
0700
0701 CounterPtr _eventCounter;
0702
0703
0704 Estimate0DPtr _xs;
0705
0706
0707 CounterPtr _xserr;
0708
0709
0710 YODA::BinnedEstimatePtr<string> _beaminfo;
0711
0712
0713 double _ntrials;
0714
0715
0716 CounterPtr _fileCounter;
0717
0718
0719 bool _isEndOfFile;
0720
0721
0722 std::pair<double,double> _userxs;
0723
0724
0725 ParticlePair _beams;
0726
0727
0728 bool _initialised;
0729
0730
0731 bool _checkBeams;
0732
0733
0734 bool _skipMultiWeights;
0735
0736
0737 std::string _matchWeightNames;
0738
0739
0740 std::string _unmatchWeightNames;
0741
0742
0743 std::string _nominalWeightName;
0744
0745
0746 double _weightCap;
0747
0748
0749
0750
0751 double _NLOSmearing;
0752
0753
0754 int _eventNumber;
0755
0756
0757 size_t _defaultWeightIdx;
0758
0759
0760 size_t _rivetDefaultWeightIdx;
0761
0762
0763 int _customDefaultWeightIdx;
0764
0765
0766 int _dumpPeriod;
0767
0768
0769 string _dumpFile;
0770
0771
0772 bool _dumping;
0773
0774
0775 ofstream _fbootstrap;
0776
0777
0778 std::string _bootstrapfilename;
0779
0780
0781 ProjectionHandler _projHandler;
0782
0783
0784
0785 };
0786
0787
0788 }
0789
0790 #endif