File indexing completed on 2025-04-19 09:06:56
0001
0002 #ifndef RIVET_Analysis_HH
0003 #define RIVET_Analysis_HH
0004
0005 #include "Rivet/Config/RivetCommon.hh"
0006 #include "Rivet/AnalysisInfo.hh"
0007 #include "Rivet/Event.hh"
0008 #include "Rivet/Projection.hh"
0009 #include "Rivet/ProjectionApplier.hh"
0010 #include "Rivet/ProjectionHandler.hh"
0011 #include "Rivet/AnalysisLoader.hh"
0012 #include "Rivet/Tools/Cuts.hh"
0013 #include "Rivet/Tools/Logging.hh"
0014 #include "Rivet/Tools/ParticleUtils.hh"
0015 #ifdef HAVE_H5
0016 #include "Rivet/Tools/RivetHDF5.hh"
0017 #endif
0018 #include "Rivet/Tools/HistoGroup.hh"
0019 #include "Rivet/Tools/RivetMT2.hh"
0020 #include "Rivet/Tools/RivetPaths.hh"
0021 #include "Rivet/Tools/RivetYODA.hh"
0022 #include "Rivet/Tools/Percentile.hh"
0023 #include "Rivet/Tools/Cutflow.hh"
0024 #include "Rivet/Projections/CentralityProjection.hh"
0025 #include <tuple>
0026
0027
0028
0029
0030 #define vetoEvent \
0031 do { MSG_DEBUG("Vetoing event on line " << __LINE__ << " of " << __FILE__); return; } while(0)
0032
0033
0034 namespace Rivet {
0035
0036
0037
0038 using std::cout;
0039 using std::cerr;
0040 using std::endl;
0041 using std::tuple;
0042 using std::stringstream;
0043 using std::swap;
0044 using std::numeric_limits;
0045
0046
0047
0048 class AnalysisHandler;
0049
0050
0051
0052
0053
0054
0055
0056
0057
0058
0059
0060
0061
0062
0063
0064
0065
0066
0067
0068
0069 class Analysis : public ProjectionApplier {
0070 public:
0071
0072
0073 friend class AnalysisHandler;
0074
0075
0076
0077 Analysis(const std::string& name);
0078
0079
0080 virtual ~Analysis() {}
0081
0082
0083 Analysis& operator = (const Analysis&) = delete;
0084
0085
0086 public:
0087
0088
0089
0090
0091
0092
0093
0094 virtual void init() { }
0095
0096
0097
0098
0099
0100 virtual void analyze(const Event& event) = 0;
0101
0102
0103
0104
0105
0106
0107 virtual void finalize() { }
0108
0109
0110
0111
0112 void syncDeclQueue() {
0113 this->_syncDeclQueue();
0114 this->markAsOwned();
0115 }
0116
0117
0118
0119
0120 public:
0121
0122
0123
0124
0125
0126
0127
0128
0129 void loadInfo() { info().parseInfoFile(); }
0130
0131
0132 const AnalysisInfo& info() const {
0133 if (!_info) throw Error("No AnalysisInfo object :-O");
0134 return *_info;
0135 }
0136
0137
0138
0139
0140
0141
0142
0143
0144 virtual std::string name() const {
0145 return ( (info().name().empty()) ? _defaultname : info().name() ) + _optstring;
0146 }
0147
0148
0149
0150
0151
0152 std::string analysisDataPath(const std::string& extn, const std::string& suffix="") {
0153 string filename = name() + (suffix.empty() ? "" : "-") + suffix + "." + extn;
0154 return findAnalysisDataFile(filename);
0155 }
0156
0157
0158 virtual std::string inspireID() const {
0159 return info().inspireID();
0160 }
0161
0162
0163 virtual std::string spiresID() const {
0164 return info().spiresID();
0165 }
0166
0167
0168
0169
0170
0171 virtual std::vector<std::string> authors() const {
0172 return info().authors();
0173 }
0174
0175
0176
0177
0178
0179
0180 virtual std::string summary() const {
0181 return info().summary();
0182 }
0183
0184
0185
0186
0187
0188
0189
0190 virtual std::string description() const {
0191 return info().description();
0192 }
0193
0194
0195
0196
0197
0198
0199 virtual std::string runInfo() const {
0200 return info().runInfo();
0201 }
0202
0203
0204 virtual std::string experiment() const {
0205 return info().experiment();
0206 }
0207
0208
0209 virtual std::string collider() const {
0210 return info().collider();
0211 }
0212
0213
0214 virtual std::string year() const {
0215 return info().year();
0216 }
0217
0218
0219 virtual double luminosityfb() const {
0220 return info().luminosityfb();
0221 }
0222
0223 virtual double luminosity() const {
0224 return info().luminosity();
0225 }
0226
0227 double luminositypb() const { return luminosity(); }
0228
0229
0230 virtual std::vector<std::string> references() const {
0231 return info().references();
0232 }
0233
0234
0235 virtual std::string bibKey() const {
0236 return info().bibKey();
0237 }
0238
0239
0240 virtual std::string bibTeX() const {
0241 return info().bibTeX();
0242 }
0243
0244
0245 virtual std::string status() const {
0246 return (info().status().empty()) ? "UNVALIDATED" : info().status();
0247 }
0248
0249
0250 virtual std::string warning() const {
0251 return info().warning();
0252 }
0253
0254
0255 virtual std::vector<std::string> todos() const {
0256 return info().todos();
0257 }
0258
0259
0260 virtual std::vector<std::string> validation() const {
0261 return info().validation();
0262 }
0263
0264
0265 virtual bool reentrant() const {
0266 return info().reentrant();
0267 }
0268
0269
0270 virtual const std::vector<std::string>& keywords() const {
0271 return info().keywords();
0272 }
0273
0274
0275 virtual std::string refMatch() const {
0276 return info().refMatch();
0277 }
0278
0279
0280 virtual std::string refUnmatch() const {
0281 return info().refUnmatch();
0282 }
0283
0284
0285 virtual std::string writerDoublePrecision() const {
0286 return info().writerDoublePrecision();
0287 }
0288
0289
0290 virtual const std::vector<PdgIdPair>& requiredBeamIDs() const {
0291 return info().beamIDs();
0292 }
0293
0294 virtual Analysis& setRequiredBeamIDs(const std::vector<PdgIdPair>& beamids) {
0295 info().setBeamIDs(beamids);
0296 return *this;
0297 }
0298
0299
0300 virtual const std::vector<std::pair<double, double> >& requiredBeamEnergies() const {
0301 return info().energies();
0302 }
0303
0304 virtual Analysis& setRequiredBeamEnergies(const std::vector<std::pair<double, double> >& energies) {
0305 info().setEnergies(energies);
0306 return *this;
0307 }
0308
0309
0310
0311 virtual std::string refFile() const {
0312 return info().refFile();
0313 }
0314
0315 virtual std::string refDataName() const {
0316 return (info().getRefDataName().empty()) ? _defaultname : info().getRefDataName();
0317 }
0318
0319 virtual void setRefDataName(const std::string& ref_data="") {
0320 info().setRefDataName(!ref_data.empty() ? ref_data : name());
0321 }
0322
0323
0324
0325
0326
0327 AnalysisInfo& info() {
0328 if (!_info) throw Error("No AnalysisInfo object :-O");
0329 return *_info;
0330 }
0331
0332
0333
0334
0335
0336
0337
0338
0339 bool merging() const {
0340 return sqrtS() <= 0.0;
0341 }
0342
0343
0344 bool compatibleWithRun() const;
0345
0346
0347
0348
0349
0350
0351
0352
0353
0354 const ParticlePair& beams() const;
0355
0356
0357 PdgIdPair beamIDs() const;
0358
0359
0360 pair<double,double> beamEnergies() const;
0361
0362
0363 vector<double> allowedEnergies() const;
0364
0365
0366 double sqrtS() const;
0367
0368
0369 bool beamsMatch(const ParticlePair& beams) const;
0370
0371
0372 bool beamsMatch(PdgId beam1, PdgId beam2, double e1, double e2) const;
0373
0374
0375 bool beamsMatch(const PdgIdPair& beams, const std::pair<double,double>& energies) const;
0376
0377
0378 bool beamIDsMatch(PdgId beam1, PdgId beam2) const;
0379
0380
0381 bool beamIDsMatch(const PdgIdPair& beamids) const;
0382
0383
0384 bool beamEnergiesMatch(double e1, double e2) const;
0385
0386
0387 bool beamEnergiesMatch(const std::pair<double,double>& energies) const;
0388
0389
0390 bool beamEnergyMatch(const std::pair<double,double>& energies) const;
0391
0392
0393 bool beamEnergyMatch(double sqrts) const;
0394
0395
0396 bool isCompatibleWithSqrtS(double energy, double tolerance=1e-5) const;
0397
0398
0399
0400
0401
0402 AnalysisHandler& handler() const { return *_analysishandler; }
0403
0404
0405 protected:
0406
0407
0408 Log& getLog() const;
0409
0410
0411 double crossSection() const;
0412
0413
0414
0415 double crossSectionPerEvent() const;
0416
0417
0418 double crossSectionError() const;
0419
0420
0421
0422 double crossSectionErrorPerEvent() const;
0423
0424
0425
0426
0427 size_t numEvents() const;
0428
0429
0430
0431
0432 double sumW() const;
0433
0434 double sumOfWeights() const { return sumW(); }
0435
0436
0437
0438
0439 double sumW2() const;
0440
0441
0442 protected:
0443
0444
0445
0446
0447
0448
0449
0450
0451 const std::string histoDir() const;
0452
0453
0454 const std::string histoPath(const std::string& hname) const;
0455
0456
0457 const std::string histoPath(unsigned int datasetID, unsigned int xAxisID, unsigned int yAxisID) const;
0458
0459
0460 const std::string mkAxisCode(unsigned int datasetID, unsigned int xAxisID, unsigned int yAxisID) const;
0461
0462
0463
0464
0465 #ifdef HAVE_H5
0466
0467
0468
0469
0470 H5::File auxFile() const {
0471 return H5::readFile(name()+".h5");
0472 }
0473
0474
0475 template <typename T>
0476 bool auxData(const string& dsname, T& rtndata) {
0477
0478 return H5::readData(name()+".h5", dsname, rtndata);
0479 }
0480
0481 template <typename T>
0482 T auxData(const string& dsname) {
0483
0484 return H5::readData<T>(name()+".h5", dsname);
0485 }
0486
0487
0488 #endif
0489
0490
0491
0492
0493
0494
0495 const std::map<std::string, YODA::AnalysisObjectPtr>& refData() const {
0496 _cacheRefData();
0497 return _refdata;
0498 }
0499
0500
0501
0502
0503 template <typename T=YODA::Estimate1D>
0504 const T& refData(const string& hname) const {
0505 _cacheRefData();
0506 MSG_TRACE("Using histo bin edges for " << name() << ":" << hname);
0507 if (!_refdata[hname]) {
0508 MSG_ERROR("Can't find reference histogram " << hname);
0509 throw Exception("Reference data " + hname + " not found.");
0510 }
0511 return dynamic_cast<T&>(*_refdata[hname]);
0512 }
0513
0514
0515
0516
0517 template <typename T=YODA::Estimate1D>
0518 const T& refData(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const {
0519 const string hname = mkAxisCode(datasetId, xAxisId, yAxisId);
0520 return refData<T>(hname);
0521 }
0522
0523
0524
0525
0526
0527
0528
0529
0530
0531
0532
0533 CounterPtr& book(CounterPtr&, const std::string& name);
0534
0535
0536
0537
0538 CounterPtr& book(CounterPtr&, unsigned int datasetID, unsigned int xAxisID, unsigned int yAxisID);
0539
0540
0541
0542
0543
0544
0545
0546
0547 Estimate0DPtr& book(Estimate0DPtr&, const std::string& name);
0548
0549
0550
0551
0552 Estimate0DPtr& book(Estimate0DPtr&, unsigned int datasetID, unsigned int xAxisID, unsigned int yAxisID);
0553
0554
0555
0556
0557
0558
0559
0560
0561 template<size_t DbnN, typename... AxisT, typename = YODA::enable_if_all_CAxisT<AxisT...>>
0562 BinnedDbnPtr<DbnN, AxisT...>& book(BinnedDbnPtr<DbnN, AxisT...>& ao,
0563 const std::string& name, const std::vector<size_t>& nbins,
0564 const std::vector<std::pair<double,double>>& loUpPairs) {
0565 if (nbins.size() != loUpPairs.size()) {
0566 throw RangeError("Vectors should have the same size!");
0567 }
0568 const string path = histoPath(name);
0569
0570 YODA::BinnedDbn<DbnN, AxisT...> yao(nbins, loUpPairs, path);
0571 _setWriterPrecision(path, yao);
0572
0573 return ao = registerAO(yao);
0574 }
0575
0576
0577 Histo1DPtr& book(Histo1DPtr& ao, const std::string& name,
0578 const size_t nbins, const double lower, const double upper) {
0579 return book(ao, name, vector<size_t>{nbins},
0580 vector<pair<double,double>>{{lower,upper}});
0581 }
0582
0583 Profile1DPtr& book(Profile1DPtr& ao, const std::string& name,
0584 const size_t nbins, const double lower, const double upper) {
0585 return book(ao, name, vector<size_t>{nbins},
0586 vector<pair<double,double>>{{lower,upper}});
0587 }
0588
0589
0590 Histo2DPtr& book(Histo2DPtr& ao, const std::string& name,
0591 const size_t nbinsX, const double lowerX, const double upperX,
0592 const size_t nbinsY, const double lowerY, const double upperY) {
0593 return book(ao, name, vector<size_t>{nbinsX,nbinsY},
0594 vector<pair<double,double>>{{lowerX,upperX}, {lowerY,upperY}});
0595 }
0596
0597 Profile2DPtr& book(Profile2DPtr& ao, const std::string& name,
0598 const size_t nbinsX, const double lowerX, const double upperX,
0599 const size_t nbinsY, const double lowerY, const double upperY) {
0600 return book(ao, name, vector<size_t>{nbinsX,nbinsY},
0601 vector<pair<double,double>>{{lowerX,upperX}, {lowerY,upperY}});
0602 }
0603
0604
0605 Histo3DPtr& book(Histo3DPtr& ao, const std::string& name,
0606 const size_t nbinsX, const double lowerX, const double upperX,
0607 const size_t nbinsY, const double lowerY, const double upperY,
0608 const size_t nbinsZ, const double lowerZ, const double upperZ) {
0609 return book(ao, name, vector<size_t>{nbinsX,nbinsY,nbinsZ},
0610 vector<pair<double,double>>{{lowerX,upperX}, {lowerY,upperY}, {lowerZ,upperZ}});
0611 }
0612
0613 Profile3DPtr& book(Profile3DPtr& ao, const std::string& name,
0614 const size_t nbinsX, const double lowerX, const double upperX,
0615 const size_t nbinsY, const double lowerY, const double upperY,
0616 const size_t nbinsZ, const double lowerZ, const double upperZ) {
0617 return book(ao, name, vector<size_t>{nbinsX,nbinsY,nbinsZ},
0618 vector<pair<double,double>>{{lowerX,upperX}, {lowerY,upperY}, {lowerZ,upperZ}});
0619 }
0620
0621
0622 template<size_t DbnN, typename... AxisT>
0623 BinnedDbnPtr<DbnN, AxisT...>& book(BinnedDbnPtr<DbnN, AxisT...>& ao, const std::string& name,
0624 const std::vector<AxisT>&... binedges) {
0625 const string path = histoPath(name);
0626 YODA::BinnedDbn<DbnN, AxisT...> yao(binedges..., path);
0627 _setWriterPrecision(path, yao);
0628
0629 return ao = registerAO(yao);
0630 }
0631
0632
0633 template<size_t DbnN, typename... AxisT>
0634 BinnedDbnPtr<DbnN, AxisT...>& book(BinnedDbnPtr<DbnN, AxisT...>& ao, const std::string& name,
0635 const std::initializer_list<AxisT>&... binedges) {
0636 return book(ao, name, vector<AxisT>{binedges} ...);
0637 }
0638
0639
0640 template<size_t DbnN, typename... AxisT>
0641 BinnedDbnPtr<DbnN, AxisT...>& book(BinnedDbnPtr<DbnN, AxisT...>& ao, const std::string& name,
0642 const YODA::BinnedEstimate<AxisT...>& refest) {
0643 const string path = histoPath(name);
0644
0645 YODA::BinnedDbn<DbnN, AxisT...> yao(refest.binning(), path);
0646 for (const string& a : yao.annotations()) {
0647 if (a != "Path") yao.rmAnnotation(a);
0648 }
0649 _setWriterPrecision(path, yao);
0650 return ao = registerAO(yao);
0651 }
0652
0653
0654 template<size_t DbnN, typename... AxisT>
0655 BinnedDbnPtr<DbnN, AxisT...>& book(BinnedDbnPtr<DbnN, AxisT...>& ao, const std::string& name) {
0656 return book(ao, name, refData<YODA::BinnedEstimate<AxisT...>>(name));
0657 }
0658
0659
0660
0661
0662 template<size_t DbnN, typename... AxisT>
0663 BinnedDbnPtr<DbnN, AxisT...>& book(BinnedDbnPtr<DbnN, AxisT...>& ao, const unsigned int datasetID,
0664 const unsigned int xAxisID, const unsigned int yAxisID) {
0665 const string name = mkAxisCode(datasetID, xAxisID, yAxisID);
0666 return book(ao, name);
0667 }
0668
0669
0670
0671
0672
0673
0674 template <typename GroupAxisT, typename... AxisT>
0675 HistoGroupPtr<GroupAxisT, AxisT...>& book(HistoGroupPtr<GroupAxisT, AxisT...>& ao,
0676 const std::vector<GroupAxisT>& edges,
0677 const std::vector<std::string>& names) {
0678 ao = make_shared<HistoGroup<GroupAxisT, AxisT...>>(edges);
0679 if (ao->numBins() != names.size()) {
0680 throw RangeError("Binning and reference-data names don't match!");
0681 }
0682 for (auto& b : ao->bins()) {
0683 const string& refname = names[b.index()-1];
0684 book(b, refname, refData<YODA::BinnedEstimate<AxisT...>>(refname));
0685 }
0686 return ao;
0687 }
0688
0689 template <typename GroupAxisT, typename... AxisT>
0690 HistoGroupPtr<GroupAxisT, AxisT...>& book(HistoGroupPtr<GroupAxisT, AxisT...>& ao,
0691 const std::vector<GroupAxisT>& edges) {
0692 return ao = make_shared<HistoGroup<GroupAxisT, AxisT...>>(edges);
0693 }
0694
0695 template <typename GroupAxisT, typename... AxisT>
0696 HistoGroupPtr<GroupAxisT, AxisT...>& book(HistoGroupPtr<GroupAxisT, AxisT...>& ao,
0697 std::initializer_list<GroupAxisT>&& edges) {
0698 return ao = make_shared<HistoGroup<GroupAxisT, AxisT...>>(std::move(edges));
0699 }
0700
0701
0702
0703
0704
0705
0706
0707 template<typename... AxisT, typename = YODA::enable_if_all_CAxisT<AxisT...>>
0708 BinnedEstimatePtr<AxisT...>& book(BinnedEstimatePtr<AxisT...>& ao,
0709 const std::string& name, const std::vector<size_t>& nbins,
0710 const std::vector<std::pair<double,double>>& loUpPairs) {
0711 if (nbins.size() != loUpPairs.size()) {
0712 throw RangeError("Vectors should have the same size!");
0713 }
0714 const string path = histoPath(name);
0715
0716 YODA::BinnedEstimate<AxisT...> yao(nbins, loUpPairs, path);
0717 _setWriterPrecision(path, yao);
0718
0719 return ao = registerAO(yao);
0720 }
0721
0722
0723
0724 Estimate1DPtr& book(Estimate1DPtr& ao, const std::string& name,
0725 const size_t nbins, const double lower, const double upper) {
0726 return book(ao, name, vector<size_t>{nbins},
0727 vector<pair<double,double>>{{lower,upper}});
0728 }
0729
0730
0731 Estimate2DPtr& book(Estimate2DPtr& ao, const std::string& name,
0732 const size_t nbinsX, const double lowerX, const double upperX,
0733 const size_t nbinsY, const double lowerY, const double upperY) {
0734 return book(ao, name, vector<size_t>{nbinsX,nbinsY},
0735 vector<pair<double,double>>{{lowerX,upperX}, {lowerY,upperY}});
0736 }
0737
0738
0739 Estimate3DPtr& book(Estimate3DPtr& ao, const std::string& name,
0740 const size_t nbinsX, const double lowerX, const double upperX,
0741 const size_t nbinsY, const double lowerY, const double upperY,
0742 const size_t nbinsZ, const double lowerZ, const double upperZ) {
0743 return book(ao, name, vector<size_t>{nbinsX,nbinsY,nbinsZ},
0744 vector<pair<double,double>>{{lowerX,upperX}, {lowerY,upperY}, {lowerZ,upperZ}});
0745 }
0746
0747
0748 template<typename... AxisT>
0749 BinnedEstimatePtr<AxisT...>& book(BinnedEstimatePtr<AxisT...>& ao, const std::string& name,
0750 const std::vector<AxisT>&... binedges) {
0751 const string path = histoPath(name);
0752 YODA::BinnedEstimate<AxisT...> yao(binedges..., path);
0753 _setWriterPrecision(path, yao);
0754
0755 return ao = registerAO(yao);
0756 }
0757
0758
0759 template<typename... AxisT>
0760 BinnedEstimatePtr<AxisT...>& book(BinnedEstimatePtr<AxisT...>& ao, const std::string& name,
0761 const std::initializer_list<AxisT>&... binedges) {
0762 return book(ao, name, vector<AxisT>{binedges} ...);
0763 }
0764
0765
0766 template<typename... AxisT>
0767 BinnedEstimatePtr<AxisT...>& book(BinnedEstimatePtr<AxisT...>& ao, const std::string& name) {
0768
0769 const string path = histoPath(name);
0770 YODA::BinnedEstimate<AxisT...> yao;
0771 try {
0772 yao = YODA::BinnedEstimate<AxisT...>(refData<YODA::BinnedEstimate<AxisT...>>(name).binning());
0773 } catch (...) {
0774 MSG_DEBUG("Couldn't retrieve reference binning, continue with nullary AO constructor.");
0775 }
0776 yao.setPath(path);
0777 _setWriterPrecision(path, yao);
0778 return ao = registerAO(yao);
0779 }
0780
0781
0782
0783
0784 template<typename... AxisT>
0785 BinnedEstimatePtr<AxisT...>& book(BinnedEstimatePtr<AxisT...>& ao, const unsigned int datasetID,
0786 const unsigned int xAxisID, const unsigned int yAxisID) {
0787 const string name = mkAxisCode(datasetID, xAxisID, yAxisID);
0788 return book(ao, name);
0789 }
0790
0791
0792
0793
0794
0795
0796
0797
0798
0799
0800
0801
0802
0803
0804
0805
0806
0807 template<size_t N>
0808 ScatterNDPtr<N>& book(ScatterNDPtr<N>& snd, const string& name, const bool copy_pts = false) {
0809 const string path = histoPath(name);
0810 YODA::ScatterND<N> scat(path);
0811 if (copy_pts) {
0812 const YODA::ScatterND<N> ref = refData<YODA::EstimateND<N-1>>(name).mkScatter();
0813 for (YODA::PointND<N> p : ref.points()) {
0814 p.setVal(N-1, 0.0);
0815 p.setErr(N-1, 0.0);
0816 scat.addPoint(std::move(p));
0817 }
0818 }
0819 _setWriterPrecision(path, scat);
0820 return snd = registerAO(scat);
0821 }
0822
0823
0824
0825
0826
0827
0828
0829
0830
0831
0832
0833 template<size_t N>
0834 ScatterNDPtr<N>& book(ScatterNDPtr<N>& snd, const unsigned int datasetID, const unsigned int xAxisID,
0835 const unsigned int yAxisID, const bool copy_pts = false) {
0836 const string axisCode = mkAxisCode(datasetID, xAxisID, yAxisID);
0837 return book(snd, axisCode, copy_pts);
0838 }
0839
0840
0841
0842
0843
0844
0845 Scatter2DPtr& book(Scatter2DPtr& snd, const string& name,
0846 const size_t npts, const double lower, const double upper) {
0847 const string path = histoPath(name);
0848
0849 Scatter2D scat(path);
0850 const double binwidth = (upper-lower)/npts;
0851 for (size_t pt = 0; pt < npts; ++pt) {
0852 const double bincentre = lower + (pt + 0.5) * binwidth;
0853 scat.addPoint(bincentre, 0, binwidth/2.0, 0);
0854 }
0855 _setWriterPrecision(path, scat);
0856
0857 return snd = registerAO(scat);
0858 }
0859
0860 Scatter3DPtr& book(Scatter3DPtr& snd, const string& name,
0861 const size_t nptsX, const double lowerX, const double upperX,
0862 const size_t nptsY, const double lowerY, const double upperY) {
0863 const string path = histoPath(name);
0864
0865 Scatter3D scat(path);
0866 const double xbinwidth = (upperX-lowerX)/nptsX;
0867 const double ybinwidth = (upperY-lowerY)/nptsY;
0868 for (size_t xpt = 0; xpt < nptsX; ++xpt) {
0869 const double xbincentre = lowerX + (xpt + 0.5) * xbinwidth;
0870 for (size_t ypt = 0; ypt < nptsY; ++ypt) {
0871 const double ybincentre = lowerY + (ypt + 0.5) * ybinwidth;
0872 scat.addPoint(xbincentre, ybincentre, 0, 0.5*xbinwidth, 0.5*ybinwidth, 0);
0873 }
0874 }
0875 _setWriterPrecision(path, scat);
0876
0877 return snd = registerAO(scat);
0878 }
0879
0880
0881
0882
0883
0884
0885 Scatter2DPtr& book(Scatter2DPtr& snd, const string& name,
0886 const std::vector<double>& binedges) {
0887 const string path = histoPath(name);
0888
0889 Scatter2D scat(path);
0890 for (size_t pt = 0; pt < binedges.size()-1; ++pt) {
0891 const double bincentre = (binedges[pt] + binedges[pt+1]) / 2.0;
0892 const double binwidth = binedges[pt+1] - binedges[pt];
0893 scat.addPoint(bincentre, 0, binwidth/2.0, 0);
0894 }
0895 _setWriterPrecision(path, scat);
0896
0897 return snd = registerAO(scat);
0898 }
0899
0900 Scatter3DPtr& book(Scatter3DPtr& snd, const string& name,
0901 const std::vector<double>& binedgesX,
0902 const std::vector<double>& binedgesY) {
0903 const string path = histoPath(name);
0904
0905 Scatter3D scat(path);
0906 for (size_t xpt = 0; xpt < binedgesX.size()-1; ++xpt) {
0907 const double xbincentre = (binedgesX[xpt] + binedgesX[xpt+1]) / 2.0;
0908 const double xbinwidth = binedgesX[xpt+1] - binedgesX[xpt];
0909 for (size_t ypt = 0; ypt < binedgesY.size()-1; ++ypt) {
0910 const double ybincentre = (binedgesY[ypt] + binedgesY[ypt+1]) / 2.0;
0911 const double ybinwidth = binedgesY[ypt+1] - binedgesY[ypt];
0912 scat.addPoint(xbincentre, ybincentre, 0, 0.5*xbinwidth, 0.5*ybinwidth, 0);
0913 }
0914 }
0915 _setWriterPrecision(path, scat);
0916
0917 return snd = registerAO(scat);
0918 }
0919
0920
0921 template<size_t N>
0922 ScatterNDPtr<N>& book(ScatterNDPtr<N>& snd, const string& name, const YODA::ScatterND<N>& refscatter) {
0923 const string path = histoPath(name);
0924
0925 YODA::ScatterND<N> scat(refscatter, path);
0926 for (const string& a : scat.annotations()) {
0927 if (a != "Path") scat.rmAnnotation(a);
0928 }
0929 _setWriterPrecision(path, scat);
0930
0931 return snd = registerAO(scat);
0932 }
0933
0934
0935
0936
0937
0938
0939
0940 CutflowPtr& book(CutflowPtr& ao, const string& name, const std::vector<std::string>& edges) {
0941 const string path = histoPath(name);
0942 Cutflow yao(edges, path);
0943 _setWriterPrecision(path, yao);
0944 return ao = registerAO(yao);
0945 }
0946
0947
0948 CutflowPtr& book(CutflowPtr& ao, const string& name, const std::initializer_list<std::string>& edges) {
0949 return book(ao, name, vector<std::string>{edges});
0950 }
0951
0952
0953
0954
0955
0956
0957 CutflowsPtr& book(CutflowsPtr& ao, const std::vector<std::string>& edges,
0958 const std::vector<std::vector<std::string>>& innerEdges) {
0959 ao = make_shared<Cutflows>(edges);
0960 if (ao->numBins() !=innerEdges.size()) {
0961 throw RangeError("Outer and Inner edges don't match");
0962 }
0963 for (auto& b : ao->bins()) {
0964 book(b, b.xEdge(), innerEdges[b.index()-1]);
0965 }
0966 return ao;
0967 }
0968
0969 CutflowsPtr& book(CutflowsPtr& ao, const std::vector<std::string>& edges) {
0970 return ao = make_shared<Cutflows>(edges);
0971 }
0972
0973 CutflowsPtr& book(CutflowsPtr& ao, std::initializer_list<std::string>&& edges) {
0974 return ao = make_shared<Cutflows>(std::move(edges));
0975 }
0976
0977
0978
0979
0980
0981
0982
0983 virtual void rawHookIn(YODA::AnalysisObjectPtr yao) {
0984 (void) yao;
0985 }
0986
0987
0988
0989
0990 virtual void rawHookOut(const vector<MultiplexAOPtr>& raos, size_t iW) {
0991 (void) raos;
0992 (void) iW;
0993 }
0994
0995
0996 public:
0997
0998
0999
1000
1001
1002 const std::map<std::string,std::string>& options() const {
1003 return _options;
1004 }
1005
1006
1007 std::string getOption(std::string optname, string def="") const {
1008 if ( _options.find(optname) != _options.end() )
1009 return _options.find(optname)->second;
1010 return def;
1011 }
1012
1013
1014
1015
1016
1017 std::string getOption(std::string optname, const char* def) {
1018 return getOption<std::string>(optname, def);
1019 }
1020
1021
1022
1023
1024
1025
1026
1027
1028 template<typename T>
1029 T getOption(std::string optname, T def) const {
1030 if (_options.find(optname) == _options.end()) return def;
1031 std::stringstream ss;
1032 ss.exceptions(std::ios::failbit);
1033 T ret;
1034 ss << _options.find(optname)->second;
1035 try {
1036 ss >> ret;
1037 } catch (...) {
1038 throw ReadError("Could not read user-provided option into requested type");
1039 }
1040 return ret;
1041 }
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055 bool getOption(std::string optname, bool def) const {
1056 if (_options.find(optname) == _options.end()) return def;
1057 const std::string val = getOption(optname);
1058 const std::string lval = toLower(val);
1059 if (lval.empty()) return false;
1060 if (lval == "true" || lval == "yes" || lval == "on") return true;
1061 if (lval == "false" || lval == "no" || lval == "off") return false;
1062 return bool(getOption<int>(optname, 0));
1063 }
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087 const CentralityProjection&
1088 declareCentrality(const SingleValueProjection &proj,
1089 string calAnaName, string calHistName,
1090 const string projName,
1091 PercentileOrder pctorder=PercentileOrder::DECREASING);
1092
1093
1094
1095
1096
1097
1098
1099
1100 template<typename T>
1101 Percentile<T> book(const string& projName,
1102 const vector<pair<double, double>>& centralityBins,
1103 const vector<tuple<size_t, size_t, size_t>>& ref) {
1104
1105 using RefT = typename ReferenceTraits<T>::RefT;
1106 using WrapT = MultiplexPtr<Multiplexer<T>>;
1107
1108 Percentile<T> pctl(this, projName);
1109
1110 const size_t nCent = centralityBins.size();
1111 for (size_t iCent = 0; iCent < nCent; ++iCent) {
1112 const string axisCode = mkAxisCode(std::get<0>(ref[iCent]),
1113 std::get<1>(ref[iCent]),
1114 std::get<2>(ref[iCent]));
1115 const RefT& refscatter = refData<RefT>(axisCode);
1116
1117 WrapT wtf(_weightNames(), T(refscatter, histoPath(axisCode)));
1118 wtf = addAnalysisObject(wtf);
1119
1120 CounterPtr cnt(_weightNames(), Counter(histoPath("TMP/COUNTER/" + axisCode)));
1121 cnt = addAnalysisObject(cnt);
1122
1123 pctl.add(wtf, cnt, centralityBins[iCent]);
1124 }
1125 return pctl;
1126 }
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
1164 private:
1165
1166
1167
1168
1169 vector<string> _weightNames() const;
1170
1171
1172 YODA::AnalysisObjectPtr _getPreload (const string& name) const;
1173
1174
1175 MultiplexAOPtr _getOtherAnalysisObject(const std::string & ananame, const std::string& name);
1176
1177
1178 void _checkBookInit() const;
1179
1180
1181 bool _inInit() const;
1182
1183
1184 bool _inFinalize() const;
1185
1186
1187 template <typename YODAT>
1188 void _setWriterPrecision(const string& path, YODAT& yao) {
1189 const string re = _info->writerDoublePrecision();
1190 if (re != "") {
1191 std::smatch match;
1192 const bool needsDP = std::regex_search(path, match, std::regex(re));
1193 if (needsDP) yao.template setAnnotation("WriterDoublePrecision", "1");
1194 }
1195 }
1196
1197
1198 private:
1199
1200
1201 class CounterAdapter {
1202 public:
1203
1204 CounterAdapter(double x) : x_(x) {}
1205
1206 CounterAdapter(const YODA::Counter& c) : x_(c.val()) {}
1207
1208 CounterAdapter(const YODA::Estimate& e) : x_(e.val()) {}
1209
1210 CounterAdapter(const YODA::Scatter1D& s) : x_(s.points()[0].x()) {
1211 if (s.numPoints() != 1) throw RangeError("Can only scale by a single value.");
1212 }
1213
1214 operator double() const { return x_; }
1215
1216 private:
1217 double x_;
1218
1219 };
1220
1221
1222 public:
1223
1224 double dbl(double x) { return x; }
1225 double dbl(const YODA::Counter& c) { return c.val(); }
1226 double dbl(const YODA::Estimate0D& e) { return e.val(); }
1227 double dbl(const YODA::Scatter1D& s) {
1228 if ( s.numPoints() != 1 ) throw RangeError("Only scatter with single value supported.");
1229 return s.points()[0].x();
1230 }
1231
1232
1233 protected:
1234
1235
1236
1237
1238
1239
1240 template<typename T>
1241 void scale(MultiplexPtr<Multiplexer<T>>& ao, CounterAdapter factor) {
1242 if (!ao) {
1243 MSG_WARNING("Failed to scale AnalysisObject=NULL in analysis "
1244 << name() << " (scale=" << double(factor) << ")");
1245 return;
1246 }
1247 if (std::isnan(double(factor)) || std::isinf(double(factor))) {
1248 MSG_WARNING("Failed to scale AnalysisObject=" << ao->path() << " in analysis: "
1249 << name() << " (invalid scale factor = " << double(factor) << ")");
1250 factor = 0;
1251 }
1252 MSG_TRACE("Scaling AnalysisObject " << ao->path() << " by factor " << double(factor));
1253 try {
1254 if constexpr( isFillable<T>::value ) {
1255 ao->scaleW(factor);
1256 }
1257 else {
1258 ao->scale(factor);
1259 }
1260 }
1261 catch (YODA::Exception& we) {
1262 MSG_WARNING("Could not scale AnalysisObject " << ao->path());
1263 return;
1264 }
1265 }
1266
1267
1268 template<typename GroupAxisT, typename... AxisT>
1269 void scale(HistoGroupPtr<GroupAxisT, AxisT...>& group, CounterAdapter factor) {
1270 if (!group) {
1271 MSG_WARNING("Failed to scale AnalysisObject=NULL in analysis "
1272 << name() << " (scale=" << double(factor) << ")");
1273 return;
1274 }
1275 if (std::isnan(double(factor)) || std::isinf(double(factor))) {
1276 MSG_WARNING("Failed to scale histo group in analysis: "
1277 << name() << " (invalid scale factor = " << double(factor) << ")");
1278 factor = 0;
1279 }
1280 MSG_TRACE("Scaling histo group by factor " << double(factor));
1281 try {
1282 group->scaleW(factor);
1283 }
1284 catch (YODA::Exception& we) {
1285 MSG_WARNING("Could not scale histo group.");
1286 return;
1287 }
1288 }
1289
1290
1291 template<typename GroupAxisT, typename... AxisT>
1292 void scale(HistoGroupPtr<GroupAxisT, AxisT...>& group, const vector<double>& factors) {
1293 if (!group) {
1294 MSG_WARNING("Failed to scale AnalysisObject=NULL in analysis " << name());
1295 return;
1296 }
1297 if (group->numBins(true) != factors.size()) {
1298 throw RangeError(name() + ": Number of scale factors does not match group binning");
1299 return;
1300 }
1301 for (auto& b : group->bins(true)) {
1302 if (!b.get()) continue;
1303 double factor = factors[b.index()];
1304 if (std::isnan(factor) || std::isinf(factor)) {
1305 MSG_WARNING("Failed to scale componment of histo group in analysis: "
1306 << name() << " (invalid scale factor = " << factor << ")");
1307 factor = 0;
1308 }
1309 MSG_TRACE("Scaling histo group element by factor " << factor);
1310 try {
1311 b->scaleW(factor);
1312 }
1313 catch (YODA::Exception& we) {
1314 MSG_WARNING("Could not scale component of histo group.");
1315 }
1316 }
1317 }
1318
1319
1320 void scale(CutflowsPtr& group, CounterAdapter factor) {
1321 if (!group) {
1322 MSG_WARNING("Failed to scale AnalysisObject=NULL in analysis "
1323 << name() << " (scale=" << double(factor) << ")");
1324 return;
1325 }
1326 if (std::isnan(double(factor)) || std::isinf(double(factor))) {
1327 MSG_WARNING("Failed to scale histo group in analysis: "
1328 << name() << " (invalid scale factor = " << double(factor) << ")");
1329 factor = 0;
1330 }
1331 MSG_TRACE("Scaling histo group by factor " << double(factor));
1332 try {
1333 group->scale(factor);
1334 }
1335 catch (YODA::Exception& we) {
1336 MSG_WARNING("Could not scale histo group.");
1337 return;
1338 }
1339 }
1340
1341
1342 template<typename T, typename U>
1343 void scale(std::map<T, U>& aos, CounterAdapter factor) {
1344 for (auto& item : aos) scale(item.second, factor);
1345 }
1346
1347
1348 template <typename AORange, typename = std::enable_if_t<YODA::isIterable<AORange>>>
1349 void scale(AORange& aos, CounterAdapter factor) {
1350 for (auto& ao : aos) scale(ao, factor);
1351 }
1352
1353
1354 template <typename T>
1355 void scale(std::initializer_list<T> aos, CounterAdapter factor) {
1356 for (auto& ao : std::vector<T>{aos}) scale(ao, factor);
1357 }
1358
1359
1360 template<typename T, typename U>
1361 void scale(std::map<T, U>& aos, const vector<double>& factors) {
1362 for (auto& item : aos) scale(item.second, factors);
1363 }
1364
1365
1366 template <typename AORange, typename = std::enable_if_t<YODA::isIterable<AORange>>>
1367 void scale(AORange& aos, const vector<double>& factors) {
1368 for (auto& ao : aos) scale(ao, factors);
1369 }
1370
1371
1372 template <typename T>
1373 void scale(std::initializer_list<T> aos, const vector<double>& factors) {
1374 for (auto& ao : std::vector<T>{aos}) scale(ao, factors);
1375 }
1376
1377
1378 template<typename GroupAxisT, typename... AxisT>
1379 void divByGroupWidth(HistoGroupPtr<GroupAxisT, AxisT...>& group) {
1380 if (!group) {
1381 MSG_WARNING("Failed to scale HistoGroup=NULL in analysis "
1382 << name() << " by group axis width");
1383 return;
1384 }
1385 group->divByGroupWidth();
1386 }
1387
1388
1389 template<typename T, typename U>
1390 void divByGroupWidth(std::map<T, U>& aos) {
1391 for (auto& item : aos) divByGroupWidth(item.second);
1392 }
1393
1394
1395 template <typename AORange, typename = std::enable_if_t<YODA::isIterable<AORange>>>
1396 void divByGroupWidth(AORange& aos) {
1397 for (auto& ao : aos) divByGroupWidth(ao);
1398 }
1399
1400
1401 template <typename T>
1402 void divByGroupWidth(std::initializer_list<T> aos) {
1403 for (auto& ao : std::vector<T>{aos}) divByGroupWidth(ao);
1404 }
1405
1406
1407
1408 template <size_t DbnN, typename... AxisT>
1409 void normalize(BinnedDbnPtr<DbnN, AxisT...> ao, const CounterAdapter norm=1.0, const bool includeoverflows=true) {
1410 if (!ao) {
1411 MSG_WARNING("Failed to normalize histo=NULL in analysis " << name() << " (norm=" << double(norm) << ")");
1412 return;
1413 }
1414 MSG_TRACE("Normalizing histo " << ao->path() << " to " << double(norm));
1415 try {
1416 const double hint = ao->integral(includeoverflows);
1417 if (hint == 0) MSG_DEBUG("Skipping histo with null area " << ao->path());
1418 else ao->normalize(norm, includeoverflows);
1419 }
1420 catch (YODA::Exception& we) {
1421 MSG_WARNING("Could not normalize histo " << ao->path());
1422 return;
1423 }
1424 }
1425
1426
1427 template <typename GroupAxisT, typename... AxisT>
1428 void normalize(HistoGroupPtr<GroupAxisT, AxisT...> group, const CounterAdapter norm=1.0, const bool includeoverflows=true) {
1429 if (!group) {
1430 MSG_WARNING("Failed to normalize histo=NULL in analysis " << name() << " (norm=" << double(norm) << ")");
1431 return;
1432 }
1433 MSG_TRACE("Normalizing histo group to " << double(norm));
1434 try {
1435 const double hint = group->integral(includeoverflows);
1436 if (hint == 0) MSG_DEBUG("Skipping histo group with null area.");
1437 else group->normalize(norm, includeoverflows);
1438 }
1439 catch (YODA::Exception& we) {
1440 MSG_WARNING("Could not normalize histo group.");
1441 return;
1442 }
1443 }
1444
1445
1446
1447 template <typename AORange, typename = std::enable_if_t<YODA::isIterable<AORange>>>
1448 void normalize(AORange& aos, const CounterAdapter norm=1.0, const bool includeoverflows=true) {
1449 for (auto& ao : aos) normalize(ao, norm, includeoverflows);
1450 }
1451
1452
1453 template<typename T>
1454 void normalize(std::initializer_list<T>&& aos, const CounterAdapter norm=1.0, const bool includeoverflows=true) {
1455 for (auto& ao : aos) normalize(ao, norm, includeoverflows);
1456 }
1457
1458
1459 template<typename T, typename U>
1460 void normalize(std::map<T, U>& aos,
1461 const CounterAdapter norm=1.0, const bool includeoverflows=true) {
1462 for (auto& item : aos) normalize(item.second, norm, includeoverflows);
1463 }
1464
1465
1466 template <typename GroupAxisT, typename... AxisT>
1467 void normalizeGroup(HistoGroupPtr<GroupAxisT, AxisT...> group, const CounterAdapter norm=1.0, const bool includeoverflows=true) {
1468 if (!group) {
1469 MSG_WARNING("Failed to normalize histo=NULL in analysis " << name() << " (norm=" << double(norm) << ")");
1470 return;
1471 }
1472 MSG_TRACE("Normalizing histo group to " << double(norm));
1473 try {
1474 const double hint = group->integral(includeoverflows);
1475 if (hint == 0) MSG_DEBUG("Skipping histo group with null area.");
1476 else group->normalizeGroup(norm, includeoverflows);
1477 }
1478 catch (YODA::Exception& we) {
1479 MSG_WARNING("Could not normalize histo group.");
1480 return;
1481 }
1482 }
1483
1484
1485 template <typename AORange, typename = std::enable_if_t<YODA::isIterable<AORange>>>
1486 void normalizeGroup(AORange& aos, const CounterAdapter norm=1.0, const bool includeoverflows=true) {
1487 for (auto& ao : aos) normalizeGroup(ao, norm, includeoverflows);
1488 }
1489
1490
1491 template<typename T>
1492 void normalizeGroup(std::initializer_list<T>&& aos, const CounterAdapter norm=1.0, const bool includeoverflows=true) {
1493 for (auto& ao : aos) normalizeGroup(ao, norm, includeoverflows);
1494 }
1495
1496
1497 template<typename T, typename U>
1498 void normalizeGroup(std::map<T, U>& aos,
1499 const CounterAdapter norm=1.0, const bool includeoverflows=true) {
1500 for (auto& item : aos) normalizeGroup(item.second, norm, includeoverflows);
1501 }
1502
1503
1504
1505
1506
1507
1508 template<size_t DbnN, typename... AxisT>
1509 void barchart(BinnedDbnPtr<DbnN, AxisT...> ao, BinnedEstimatePtr<AxisT...> est) const {
1510 const string path = est->path();
1511 *est = ao->mkEstimate(path, "stats", false);
1512 }
1513
1514
1515
1516
1517 void divide(CounterPtr c1, CounterPtr c2, Estimate0DPtr est) const;
1518
1519
1520
1521
1522 void divide(const YODA::Counter& c1, const YODA::Counter& c2, Estimate0DPtr est) const;
1523
1524
1525
1526
1527 void divide(Estimate0DPtr e1, Estimate0DPtr e2, Estimate0DPtr est) const;
1528
1529
1530
1531
1532 void divide(const YODA::Estimate0D& e1, const YODA::Estimate0D& e2, Estimate0DPtr est) const;
1533
1534
1535
1536
1537
1538 template<size_t DbnN, typename... AxisT>
1539 void divide(const YODA::BinnedDbn<DbnN, AxisT...>& h1, const YODA::BinnedDbn<DbnN, AxisT...>& h2,
1540 BinnedEstimatePtr<AxisT...> est) const {
1541 const string path = est->path();
1542 *est = h1 / h2;
1543 est->setPath(path);
1544 }
1545
1546 template<size_t DbnN, typename... AxisT>
1547 void divide(BinnedDbnPtr<DbnN, AxisT...> h1, BinnedDbnPtr<DbnN, AxisT...> h2,
1548 BinnedEstimatePtr<AxisT...> est) const {
1549 return divide(*h1, *h2, est);
1550 }
1551
1552
1553
1554
1555 template<typename... AxisT>
1556 void divide(const YODA::BinnedEstimate<AxisT...>& e1, const YODA::BinnedEstimate<AxisT...>& e2,
1557 BinnedEstimatePtr<AxisT...> est) const {
1558 const string path = est->path();
1559 *est = e1 / e2;
1560 est->setPath(path);
1561 }
1562
1563 template<typename... AxisT>
1564 void divide(BinnedEstimatePtr<AxisT...> e1, BinnedEstimatePtr<AxisT...> e2,
1565 BinnedEstimatePtr<AxisT...> est) const {
1566 return divide(*e1, *e2, est);
1567 }
1568
1569
1570
1571
1572
1573
1574 void efficiency(CounterPtr c1, CounterPtr c2, Estimate0DPtr est) const {
1575 efficiency(*c1, *c2, est);
1576 }
1577
1578
1579
1580
1581 void efficiency(const YODA::Counter& c1, const YODA::Counter& c2, Estimate0DPtr est) const {
1582 const string path = est->path();
1583 *est = YODA::efficiency(c1, c2);
1584 est->setPath(path);
1585 }
1586
1587
1588
1589
1590
1591
1592 template<size_t DbnN, typename... AxisT>
1593 void efficiency(const YODA::BinnedDbn<DbnN, AxisT...>& h1, const YODA::BinnedDbn<DbnN, AxisT...>& h2,
1594 BinnedEstimatePtr<AxisT...> est) const {
1595 const string path = est->path();
1596 *est = YODA::efficiency(h1, h2);
1597 est->setPath(path);
1598 }
1599
1600 template<size_t DbnN, typename... AxisT>
1601 void efficiency(BinnedDbnPtr<DbnN, AxisT...> h1, BinnedDbnPtr<DbnN, AxisT...> h2,
1602 BinnedEstimatePtr<AxisT...> est) const {
1603 efficiency(*h1, *h2, est);
1604 }
1605
1606
1607
1608
1609
1610 template<typename... AxisT>
1611 void efficiency(const YODA::BinnedEstimate<AxisT...>& e1, const YODA::BinnedEstimate<AxisT...>& e2,
1612 BinnedEstimatePtr<AxisT...> est) const {
1613 const string path = est->path();
1614 *est = YODA::efficiency(e1, e2);
1615 est->setPath(path);
1616 }
1617
1618 template<typename... AxisT>
1619 void efficiency(BinnedEstimatePtr<AxisT...> e1, BinnedEstimatePtr<AxisT...> e2,
1620 BinnedEstimatePtr<AxisT...> est) const {
1621 efficiency(*e1, *e2, est);
1622 }
1623
1624
1625
1626
1627
1628 template<size_t DbnN, typename... AxisT>
1629 void asymm(const YODA::BinnedDbn<DbnN, AxisT...>& h1, const YODA::BinnedDbn<DbnN, AxisT...>& h2,
1630 BinnedEstimatePtr<AxisT...> est) const {
1631 const string path = est->path();
1632 *est = YODA::asymm(h1, h2);
1633 est->setPath(path);
1634 }
1635
1636 template<size_t DbnN, typename... AxisT>
1637 void asymm(BinnedDbnPtr<DbnN, AxisT...> h1, BinnedDbnPtr<DbnN, AxisT...> h2,
1638 BinnedEstimatePtr<AxisT...> est) const {
1639 asymm(*h1, *h2, est);
1640 }
1641
1642
1643
1644
1645 template<typename... AxisT>
1646 void asymm(const YODA::BinnedEstimate<AxisT...>& e1, const YODA::BinnedEstimate<AxisT...>& e2,
1647 BinnedEstimatePtr<AxisT...> est) const {
1648 const string path = est->path();
1649 *est = YODA::asymm(e1, e2);
1650 est->setPath(path);
1651 }
1652
1653 template<typename... AxisT>
1654 void asymm(BinnedEstimatePtr<AxisT...> e1, BinnedEstimatePtr<AxisT...> e2,
1655 BinnedEstimatePtr<AxisT...> est) const {
1656 asymm(*e1, *e2, est);
1657 }
1658
1659
1660
1661
1662 template<size_t DbnN, typename... AxisT>
1663 void integrate(const YODA::BinnedDbn<DbnN, AxisT...>& h, BinnedEstimatePtr<AxisT...> est) const {
1664 const string path = est->path();
1665 *est = mkIntegral(h);
1666 est->setPath(path);
1667 }
1668
1669 template<size_t DbnN, typename... AxisT>
1670 void integrate(BinnedDbnPtr<DbnN, AxisT...>& h, BinnedEstimatePtr<AxisT...> est) const {
1671 integrate(*h, est);
1672 }
1673
1674
1675
1676
1677 public:
1678
1679
1680 const vector<MultiplexAOPtr>& analysisObjects() const {
1681 return _analysisobjects;
1682 }
1683
1684
1685 protected:
1686
1687
1688
1689
1690
1691 size_t defaultWeightIndex() const;
1692
1693
1694 template <typename YODAT>
1695 shared_ptr<YODAT> getPreload(const string& path) const {
1696 return dynamic_pointer_cast<YODAT>(_getPreload(path));
1697 }
1698
1699
1700
1701 template <typename YODAT>
1702 MultiplexPtr< Multiplexer<YODAT> > registerAO(const YODAT& yao) {
1703 using MultiplexerT = Multiplexer<YODAT>;
1704 using YODAPtrT = shared_ptr<YODAT>;
1705 using RAOT = MultiplexPtr<MultiplexerT>;
1706
1707 if ( !_inInit() && !_inFinalize() ) {
1708 MSG_ERROR("Can't book objects outside of init() or finalize()");
1709 throw UserError(name() + ": Can't book objects outside of init() or finalize().");
1710 }
1711
1712
1713
1714
1715 for (auto& waold : analysisObjects()) {
1716 if ( yao.path() == waold.get()->basePath() ) {
1717 const string msg = "Found double-booking of " + yao.path() + " in " + name();
1718 if ( _inInit() ) {
1719 MSG_ERROR(msg);
1720 throw LookupError(msg);
1721 } else {
1722 MSG_WARNING(msg + ". Keeping previous booking");
1723 }
1724 return RAOT(dynamic_pointer_cast<MultiplexerT>(waold.get()));
1725 }
1726 }
1727
1728 shared_ptr<MultiplexerT> wao = make_shared<MultiplexerT>();
1729 wao->_basePath = yao.path();
1730 YODAPtrT yaop = make_shared<YODAT>(yao);
1731
1732 for (const string& weightname : _weightNames()) {
1733
1734
1735 string finalpath = yao.path();
1736 if ( weightname != "" ) finalpath += "[" + weightname + "]";
1737 YODAPtrT preload = getPreload<YODAT>(finalpath);
1738 if ( preload ) {
1739 if ( !bookingCompatible(preload, yaop) ) {
1740
1741 MSG_WARNING("Found incompatible pre-existing data object with same base path "
1742 << finalpath << " for " << name());
1743 preload = nullptr;
1744 } else {
1745 MSG_TRACE("Using preloaded " << finalpath << " in " <<name());
1746 wao->_final.push_back(make_shared<YODAT>(*preload));
1747 }
1748 }
1749 else {
1750 wao->_final.push_back(make_shared<YODAT>(yao));
1751 wao->_final.back()->setPath(finalpath);
1752 }
1753
1754
1755 string rawpath = "/RAW" + finalpath;
1756 preload = getPreload<YODAT>(rawpath);
1757 if ( preload ) {
1758 if ( !bookingCompatible(preload, yaop) ) {
1759 MSG_WARNING("Found incompatible pre-existing data object with same base path "
1760 << rawpath << " for " << name());
1761 preload = nullptr;
1762 } else {
1763 MSG_TRACE("Using preloaded " << rawpath << " in " <<name());
1764 wao->_persistent.push_back(make_shared<YODAT>(*preload));
1765 }
1766 }
1767 else {
1768 wao->_persistent.push_back(make_shared<YODAT>(yao));
1769 wao->_persistent.back()->setPath(rawpath);
1770 }
1771 }
1772 MultiplexPtr<MultiplexerT> ret(wao);
1773
1774 ret.get()->unsetActiveWeight();
1775 if ( _inFinalize() ) {
1776
1777
1778 ret.get()->pushToFinal();
1779 ret.get()->setActiveFinalWeightIdx(0);
1780 }
1781 _analysisobjects.push_back(ret);
1782
1783 return ret;
1784 }
1785
1786
1787
1788 template <typename AO=MultiplexAOPtr>
1789 AO addAnalysisObject(const AO& aonew) {
1790 _checkBookInit();
1791
1792 for (const MultiplexAOPtr& ao : analysisObjects()) {
1793
1794
1795 ao.get()->setActiveWeightIdx(defaultWeightIndex());
1796 aonew.get()->setActiveWeightIdx(defaultWeightIndex());
1797 if (ao->path() != aonew->path()) continue;
1798
1799
1800
1801 AO aoold = AO(dynamic_pointer_cast<typename AO::value_type>(ao.get()));
1802 if ( !aoold || !bookingCompatible(aonew, aoold) ) {
1803 MSG_WARNING("Found incompatible pre-existing data object with same base path "
1804 << aonew->path() << " for " << name());
1805 throw LookupError("Found incompatible pre-existing data object with same base path during AO booking");
1806 }
1807
1808
1809 for (size_t weightIdx = 0; weightIdx < _weightNames().size(); ++weightIdx) {
1810 aoold.get()->setActiveWeightIdx(weightIdx);
1811 aonew.get()->setActiveWeightIdx(weightIdx);
1812 if (aoold->path() != aonew->path()) {
1813 MSG_WARNING("Found incompatible pre-existing data object with different weight-path "
1814 << aonew->path() << " for " << name());
1815 throw LookupError("Found incompatible pre-existing data object with same weight-path during AO booking");
1816 }
1817 }
1818
1819
1820 aoold.get()->unsetActiveWeight();
1821 MSG_TRACE("Bound pre-existing data object " << aoold->path() << " for " << name());
1822 return aoold;
1823 }
1824
1825
1826 MSG_TRACE("Registered " << aonew->annotation("Type") << " " << aonew->path() << " for " << name());
1827 aonew.get()->unsetActiveWeight();
1828
1829 _analysisobjects.push_back(aonew);
1830 return aonew;
1831 }
1832
1833
1834 void removeAnalysisObject(const std::string& path);
1835
1836
1837 void removeAnalysisObject(const MultiplexAOPtr& ao);
1838
1839
1840 template <typename AO=MultiplexAOPtr>
1841 const AO getAnalysisObject(const std::string& aoname) const {
1842 for (const MultiplexAOPtr& ao : analysisObjects()) {
1843 ao.get()->setActiveWeightIdx(defaultWeightIndex());
1844 if (ao->path() == histoPath(aoname)) {
1845
1846 return AO(dynamic_pointer_cast<typename AO::value_type>(ao.get()));
1847 }
1848 }
1849 throw LookupError("Data object " + histoPath(aoname) + " not found");
1850 }
1851
1852
1853
1854
1855
1856
1857
1858
1859
1860
1861
1862
1863
1864
1865
1866
1867
1868
1869
1870
1871
1872
1873
1874 template <typename AO=MultiplexAOPtr>
1875 AO getAnalysisObject(const std::string& ananame,
1876 const std::string& aoname) {
1877 MultiplexAOPtr ao = _getOtherAnalysisObject(ananame, aoname);
1878
1879 return AO(dynamic_pointer_cast<typename AO::value_type>(ao.get()));
1880 }
1881
1882
1883
1884
1885 private:
1886
1887
1888 string _defaultname;
1889
1890
1891 unique_ptr<AnalysisInfo> _info;
1892
1893
1894
1895 vector<MultiplexAOPtr> _analysisobjects;
1896
1897
1898
1899 double _crossSection;
1900 bool _gotCrossSection;
1901
1902
1903
1904 AnalysisHandler* _analysishandler;
1905
1906
1907
1908 mutable std::map<std::string, YODA::AnalysisObjectPtr> _refdata;
1909
1910
1911 map<string, string> _options;
1912
1913
1914 string _optstring;
1915
1916
1917 private:
1918
1919
1920
1921
1922
1923 void _cacheRefData() const;
1924
1925
1926
1927 };
1928
1929
1930
1931
1932
1933
1934
1935
1936
1937 }
1938
1939
1940
1941 #include "Rivet/AnalysisBuilder.hh"
1942
1943
1944
1945
1946
1947
1948
1949 #define RIVET_DECLARE_PLUGIN(clsname) ::Rivet::AnalysisBuilder<clsname> plugin_ ## clsname
1950
1951
1952
1953 #define RIVET_DECLARE_ALIASED_PLUGIN(clsname, alias) RIVET_DECLARE_PLUGIN(clsname)( #alias )
1954
1955
1956
1957 #define RIVET_DEFAULT_ANALYSIS_CTOR(clsname) clsname() : Analysis(# clsname) {}
1958
1959
1960
1961 #define RIVET_REGISTER_TYPE(...) handler().registerType<__VA_ARGS__>()
1962
1963
1964
1965 #define RIVET_REGISTER_BINNED_SET(...) { \
1966 RIVET_REGISTER_TYPE(YODA::BinnedHisto<__VA_ARGS__>); \
1967 RIVET_REGISTER_TYPE(YODA::BinnedProfile<__VA_ARGS__>); \
1968 RIVET_REGISTER_TYPE(YODA::BinnedEstimate<__VA_ARGS__>); }
1969
1970
1971
1972
1973 #endif