diff --git a/analyses/pluginALICE/ALICE_2013_I1225979.cc b/analyses/pluginALICE/ALICE_2013_I1225979.cc new file mode 100644 --- /dev/null +++ b/analyses/pluginALICE/ALICE_2013_I1225979.cc @@ -0,0 +1,106 @@ +// -*- C++ -*- +#include "Rivet/Analysis.hh" +#include "Rivet/Projections/ChargedFinalState.hh" +#include "Rivet/Tools/AliceCommon.hh" +#include "Rivet/Projections/AliceCommon.hh" +namespace Rivet { + + + /// @brief ALICE PbPb at 2.76 TeV eta distributions. + class ALICE_2013_I1225979 : public Analysis { + public: + + /// Constructor + DEFAULT_RIVET_ANALYSIS_CTOR(ALICE_2013_I1225979); + + + /// @name Analysis methods + //@{ + + /// Book histograms and initialise projections before the run + void init() { + // Initialise and register projections + // Centrality projection. + declareCentrality(ALICE::V0MMultiplicity(), "ALICE_2015_PBPBCentrality", + "V0M","V0M"); + // Projections for the 2-out-of-3 trigger. + declare(ChargedFinalState( (Cuts::eta > 2.8 && Cuts::eta < 5.1) && + Cuts::pT > 0.1*GeV), "VZERO1"); + declare(ChargedFinalState( (Cuts::eta > -3.7 && Cuts::eta < -1.7) && + Cuts::pT > 0.1*GeV), "VZERO2"); + declare(ChargedFinalState(Cuts::abseta < 1. && Cuts::pT > 0.15*GeV), + "SPD"); + + // Primary particles. + declare(ALICE::PrimaryParticles(Cuts::abseta < 5),"APRIM"); + + // The centrality bins upper bin edges. + centralityBins = { 5., 10., 20., 30. }; + // Centrality histograms and corresponding sow counters. + for (int i = 0; i < 4; ++i) { + histEta[centralityBins[i]] = bookHisto1D(1, 1, i + 1); + sow[centralityBins[i]] = bookCounter("sow_" + toString(i)); + } + } + + + /// Perform the per-event analysis + void analyze(const Event& event) { + const double weight = event.weight(); + // Trigger projections. + const ChargedFinalState& vz1 = + applyProjection(event,"VZERO1"); + const ChargedFinalState& vz2 = + applyProjection(event,"VZERO2"); + const ChargedFinalState& spd = + applyProjection(event,"SPD"); + int fwdTrig = (vz1.particles().size() > 0 ? 1 : 0); + int bwdTrig = (vz2.particles().size() > 0 ? 1 : 0); + int cTrig = (spd.particles().size() > 0 ? 1 : 0); + + if (fwdTrig + bwdTrig + cTrig < 2) vetoEvent; + // We must have direct acces to the centrality projection. + const CentralityProjection& cent = apply(event,"V0M"); + double c = cent(); + // Find the correct centrality histogram + auto hItr = histEta.upper_bound(c); + if (hItr == histEta.end()) return; + // Find the correct sow. + auto sItr = sow.upper_bound(c); + if (sItr == sow.end()) return; + sItr->second->fill(weight); + + // Fill the histograms. + for ( const auto& p : + applyProjection(event,"APRIM").particles() ) + if(p.abscharge() > 0) hItr->second->fill(p.eta(), weight); + + } + + + /// Normalise histograms etc., after the run + void finalize() { + for (int i = 0; i < 4; ++i) + histEta[centralityBins[i]]->scaleW(1./sow[centralityBins[i]]->sumW()); + + } + + //@} + + + /// @name Histograms + //@{ + vector centralityBins; + map histEta; + map sow; + //@} + + + }; + + + // The hook for the plugin system + DECLARE_RIVET_PLUGIN(ALICE_2013_I1225979); + + +} diff --git a/analyses/pluginALICE/ALICE_2013_I1225979.info b/analyses/pluginALICE/ALICE_2013_I1225979.info new file mode 100644 --- /dev/null +++ b/analyses/pluginALICE/ALICE_2013_I1225979.info @@ -0,0 +1,41 @@ +Name: ALICE_2013_I1225979 +Year: 2013 +Summary: Charged particle production as function of centrality, central events only, in PbPb collisions at 2.76 TeV. +Experiment: ALICE +Collider: LHC +InspireID: 1225979 +Status: UNVALIDATED +Authors: + - Christian Bierlich +References: +- Phys.Lett.B726(2013)610-622 +- DOI:10.1016/j.physletb.2013.09.022 +- arXiv:1304.0347 +RunInfo: PbPb minimum bias events. The analysis holds the Primary Particle definition, so don't limit decays on generator level. +NeedCrossSection: no +Beams: [1000822080, 1000822080] +# This is _total_ energy of beams, so this becomes 208*2760=574080 +Energies: [574080] +Description: + 'Charged particle pseudorapidity density in centrality classes 0-5,5-10,10-20,20-30. Measurements cover a wide $\eta$ range from -3.5 to 5. Centrality classes refer to forward V0 spectrum, as also measured by ALICE, can be modified to use a user definition instead.' +BibKey: Abbas:2013bpa +BibTeX: '@article{Abbas:2013bpa, + author = "Abbas, Ehab and others", + title = "{Centrality dependence of the pseudorapidity density + distribution for charged particles in Pb-Pb collisions at + $\sqrt{s_{\rm NN}}$ = 2.76 TeV}", + collaboration = "ALICE", + journal = "Phys. Lett.", + volume = "B726", + year = "2013", + pages = "610-622", + doi = "10.1016/j.physletb.2013.09.022", + eprint = "1304.0347", + archivePrefix = "arXiv", + primaryClass = "nucl-ex", + reportNumber = "CERN-PH-EP-2013-045", + SLACcitation = "%%CITATION = ARXIV:1304.0347;%%" +}' +ToDo: + - The corresponding centrality analysis. + diff --git a/analyses/pluginALICE/ALICE_2013_I1225979.plot b/analyses/pluginALICE/ALICE_2013_I1225979.plot new file mode 100644 --- /dev/null +++ b/analyses/pluginALICE/ALICE_2013_I1225979.plot @@ -0,0 +1,8 @@ +# BEGIN PLOT /ALICE_2013_I1225979/d01-x01-y01 +#Title=[Uncomment and insert title for histogram d01-x01-y01 here] +#XLabel=[Uncomment and insert x-axis label for histogram d01-x01-y01 here] +#YLabel=[Uncomment and insert y-axis label for histogram d01-x01-y01 here] +# + any additional plot settings you might like, see make-plots documentation +# END PLOT + +# ... add more histograms as you need them ... diff --git a/analyses/pluginALICE/ALICE_2013_I1225979.yoda b/analyses/pluginALICE/ALICE_2013_I1225979.yoda new file mode 100644 --- /dev/null +++ b/analyses/pluginALICE/ALICE_2013_I1225979.yoda @@ -0,0 +1,349 @@ +BEGIN YODA_SCATTER2D_V2 /REF/ALICE_2013_I1225979/d01-x01-y01 +IsRef: 1 +Path: /REF/ALICE_2013_I1225979/d01-x01-y01 +Title: doi:10.17182/hepdata.68753.v1/t1 +Type: Scatter2D +--- +# xval xerr- xerr+ yval yerr- yerr+ +-4.875000e+00 1.250000e-01 1.250000e-01 8.880000e+02 3.900000e+01 3.900000e+01 +-4.625000e+00 1.250000e-01 1.250000e-01 9.740000e+02 4.100000e+01 4.100000e+01 +-4.375000e+00 1.250000e-01 1.250000e-01 1.046000e+03 4.400000e+01 4.400000e+01 +-4.125000e+00 1.250000e-01 1.250000e-01 1.128000e+03 4.600000e+01 4.600000e+01 +-3.875000e+00 1.250000e-01 1.250000e-01 1.209000e+03 5.000000e+01 5.000000e+01 +-3.625000e+00 1.250000e-01 1.250000e-01 1.304000e+03 5.400000e+01 5.400000e+01 +-3.375000e+00 1.250000e-01 1.250000e-01 1.388000e+03 5.700000e+01 5.700000e+01 +-3.125000e+00 1.250000e-01 1.250000e-01 1.449000e+03 5.900000e+01 5.900000e+01 +-2.875000e+00 1.250000e-01 1.250000e-01 1.504000e+03 5.800000e+01 5.800000e+01 +-2.625000e+00 1.250000e-01 1.250000e-01 1.578000e+03 6.100000e+01 6.100000e+01 +-2.375000e+00 1.250000e-01 1.250000e-01 1.627000e+03 6.200000e+01 6.200000e+01 +-2.125000e+00 1.250000e-01 1.250000e-01 1.674000e+03 6.200000e+01 6.200000e+01 +-1.875000e+00 1.250000e-01 1.250000e-01 1.709000e+03 6.600000e+01 6.600000e+01 +-1.625000e+00 1.250000e-01 1.250000e-01 1.736000e+03 6.400000e+01 6.400000e+01 +-1.375000e+00 1.250000e-01 1.250000e-01 1.739000e+03 6.500000e+01 6.500000e+01 +-1.125000e+00 1.250000e-01 1.250000e-01 1.699000e+03 6.800000e+01 6.800000e+01 +-8.750000e-01 1.250000e-01 1.250000e-01 1.674000e+03 4.600000e+01 4.600000e+01 +-6.250000e-01 1.250000e-01 1.250000e-01 1.658000e+03 4.000000e+01 4.000000e+01 +-3.750000e-01 1.250000e-01 1.250000e-01 1.627000e+03 3.900000e+01 3.900000e+01 +-1.250000e-01 1.250000e-01 1.250000e-01 1.615000e+03 3.900000e+01 3.900000e+01 +1.250000e-01 1.250000e-01 1.250000e-01 1.615000e+03 3.900000e+01 3.900000e+01 +3.750000e-01 1.250000e-01 1.250000e-01 1.627000e+03 3.900000e+01 3.900000e+01 +6.250000e-01 1.250000e-01 1.250000e-01 1.658000e+03 4.000000e+01 4.000000e+01 +8.750000e-01 1.250000e-01 1.250000e-01 1.674000e+03 4.600000e+01 4.600000e+01 +1.125000e+00 1.250000e-01 1.250000e-01 1.699000e+03 6.800000e+01 6.800000e+01 +1.375000e+00 1.250000e-01 1.250000e-01 1.739000e+03 6.500000e+01 6.500000e+01 +1.625000e+00 1.250000e-01 1.250000e-01 1.736000e+03 6.400000e+01 6.400000e+01 +1.875000e+00 1.250000e-01 1.250000e-01 1.709000e+03 6.600000e+01 6.600000e+01 +2.125000e+00 1.250000e-01 1.250000e-01 1.674000e+03 6.200000e+01 6.200000e+01 +2.375000e+00 1.250000e-01 1.250000e-01 1.627000e+03 6.200000e+01 6.200000e+01 +2.625000e+00 1.250000e-01 1.250000e-01 1.578000e+03 6.100000e+01 6.100000e+01 +2.875000e+00 1.250000e-01 1.250000e-01 1.504000e+03 5.800000e+01 5.800000e+01 +3.125000e+00 1.250000e-01 1.250000e-01 1.449000e+03 5.900000e+01 5.900000e+01 +3.375000e+00 1.250000e-01 1.250000e-01 1.388000e+03 5.700000e+01 5.700000e+01 +3.625000e+00 1.250000e-01 1.250000e-01 1.304000e+03 5.400000e+01 5.400000e+01 +3.875000e+00 1.250000e-01 1.250000e-01 1.209000e+03 5.000000e+01 5.000000e+01 +4.125000e+00 1.250000e-01 1.250000e-01 1.128000e+03 4.600000e+01 4.600000e+01 +4.375000e+00 1.250000e-01 1.250000e-01 1.046000e+03 4.400000e+01 4.400000e+01 +4.625000e+00 1.250000e-01 1.250000e-01 9.740000e+02 4.100000e+01 4.100000e+01 +4.875000e+00 1.250000e-01 1.250000e-01 8.880000e+02 3.900000e+01 3.900000e+01 +5.125000e+00 1.250000e-01 1.250000e-01 7.700000e+02 4.200000e+01 4.200000e+01 +5.375000e+00 1.250000e-01 1.250000e-01 6.670000e+02 3.900000e+01 3.900000e+01 +END YODA_SCATTER2D_V2 +BEGIN YODA_SCATTER2D_V2 /REF/ALICE_2013_I1225979/d01-x01-y02 +IsRef: 1 +Path: /REF/ALICE_2013_I1225979/d01-x01-y02 +Title: doi:10.17182/hepdata.68753.v1/t1 +Type: Scatter2D +--- +# xval xerr- xerr+ yval yerr- yerr+ +-4.875000e+00 1.250000e-01 1.250000e-01 7.440000e+02 3.400000e+01 3.400000e+01 +-4.625000e+00 1.250000e-01 1.250000e-01 8.070000e+02 3.500000e+01 3.500000e+01 +-4.375000e+00 1.250000e-01 1.250000e-01 8.610000e+02 3.700000e+01 3.700000e+01 +-4.125000e+00 1.250000e-01 1.250000e-01 9.290000e+02 3.900000e+01 3.900000e+01 +-3.875000e+00 1.250000e-01 1.250000e-01 1.001000e+03 4.300000e+01 4.300000e+01 +-3.625000e+00 1.250000e-01 1.250000e-01 1.072000e+03 4.600000e+01 4.600000e+01 +-3.375000e+00 1.250000e-01 1.250000e-01 1.144000e+03 4.900000e+01 4.900000e+01 +-3.125000e+00 1.250000e-01 1.250000e-01 1.196000e+03 5.100000e+01 5.100000e+01 +-2.875000e+00 1.250000e-01 1.250000e-01 1.243000e+03 5.100000e+01 5.100000e+01 +-2.625000e+00 1.250000e-01 1.250000e-01 1.298000e+03 5.200000e+01 5.200000e+01 +-2.375000e+00 1.250000e-01 1.250000e-01 1.333000e+03 5.400000e+01 5.400000e+01 +-2.125000e+00 1.250000e-01 1.250000e-01 1.361000e+03 5.400000e+01 5.400000e+01 +-1.875000e+00 1.250000e-01 1.250000e-01 1.393000e+03 5.800000e+01 5.800000e+01 +-1.625000e+00 1.250000e-01 1.250000e-01 1.422000e+03 5.600000e+01 5.600000e+01 +-1.375000e+00 1.250000e-01 1.250000e-01 1.418000e+03 5.600000e+01 5.600000e+01 +-1.125000e+00 1.250000e-01 1.250000e-01 1.380000e+03 5.800000e+01 5.800000e+01 +-8.750000e-01 1.250000e-01 1.250000e-01 1.362000e+03 3.700000e+01 3.700000e+01 +-6.250000e-01 1.250000e-01 1.250000e-01 1.338000e+03 3.200000e+01 3.200000e+01 +-3.750000e-01 1.250000e-01 1.250000e-01 1.330000e+03 3.200000e+01 3.200000e+01 +-1.250000e-01 1.250000e-01 1.250000e-01 1.318000e+03 3.200000e+01 3.200000e+01 +1.250000e-01 1.250000e-01 1.250000e-01 1.318000e+03 3.200000e+01 3.200000e+01 +3.750000e-01 1.250000e-01 1.250000e-01 1.330000e+03 3.200000e+01 3.200000e+01 +6.250000e-01 1.250000e-01 1.250000e-01 1.338000e+03 3.200000e+01 3.200000e+01 +8.750000e-01 1.250000e-01 1.250000e-01 1.362000e+03 3.700000e+01 3.700000e+01 +1.125000e+00 1.250000e-01 1.250000e-01 1.380000e+03 5.800000e+01 5.800000e+01 +1.375000e+00 1.250000e-01 1.250000e-01 1.418000e+03 5.600000e+01 5.600000e+01 +1.625000e+00 1.250000e-01 1.250000e-01 1.422000e+03 5.600000e+01 5.600000e+01 +1.875000e+00 1.250000e-01 1.250000e-01 1.393000e+03 5.800000e+01 5.800000e+01 +2.125000e+00 1.250000e-01 1.250000e-01 1.361000e+03 5.400000e+01 5.400000e+01 +2.375000e+00 1.250000e-01 1.250000e-01 1.333000e+03 5.400000e+01 5.400000e+01 +2.625000e+00 1.250000e-01 1.250000e-01 1.298000e+03 5.200000e+01 5.200000e+01 +2.875000e+00 1.250000e-01 1.250000e-01 1.243000e+03 5.100000e+01 5.100000e+01 +3.125000e+00 1.250000e-01 1.250000e-01 1.196000e+03 5.100000e+01 5.100000e+01 +3.375000e+00 1.250000e-01 1.250000e-01 1.144000e+03 4.900000e+01 4.900000e+01 +3.625000e+00 1.250000e-01 1.250000e-01 1.072000e+03 4.600000e+01 4.600000e+01 +3.875000e+00 1.250000e-01 1.250000e-01 1.001000e+03 4.300000e+01 4.300000e+01 +4.125000e+00 1.250000e-01 1.250000e-01 9.290000e+02 3.900000e+01 3.900000e+01 +4.375000e+00 1.250000e-01 1.250000e-01 8.610000e+02 3.700000e+01 3.700000e+01 +4.625000e+00 1.250000e-01 1.250000e-01 8.070000e+02 3.500000e+01 3.500000e+01 +4.875000e+00 1.250000e-01 1.250000e-01 7.440000e+02 3.400000e+01 3.400000e+01 +5.125000e+00 1.250000e-01 1.250000e-01 6.510000e+02 3.900000e+01 3.900000e+01 +5.375000e+00 1.250000e-01 1.250000e-01 5.770000e+02 3.500000e+01 3.500000e+01 +END YODA_SCATTER2D_V2 +BEGIN YODA_SCATTER2D_V2 /REF/ALICE_2013_I1225979/d01-x01-y03 +IsRef: 1 +Path: /REF/ALICE_2013_I1225979/d01-x01-y03 +Title: doi:10.17182/hepdata.68753.v1/t1 +Type: Scatter2D +--- +# xval xerr- xerr+ yval yerr- yerr+ +-4.875000e+00 1.250000e-01 1.250000e-01 5.590000e+02 2.900000e+01 2.900000e+01 +-4.625000e+00 1.250000e-01 1.250000e-01 6.150000e+02 3.100000e+01 3.100000e+01 +-4.375000e+00 1.250000e-01 1.250000e-01 6.560000e+02 3.200000e+01 3.200000e+01 +-4.125000e+00 1.250000e-01 1.250000e-01 7.050000e+02 3.300000e+01 3.300000e+01 +-3.875000e+00 1.250000e-01 1.250000e-01 7.540000e+02 3.600000e+01 3.600000e+01 +-3.625000e+00 1.250000e-01 1.250000e-01 8.120000e+02 3.900000e+01 3.900000e+01 +-3.375000e+00 1.250000e-01 1.250000e-01 8.620000e+02 4.200000e+01 4.200000e+01 +-3.125000e+00 1.250000e-01 1.250000e-01 8.990000e+02 4.300000e+01 4.300000e+01 +-2.875000e+00 1.250000e-01 1.250000e-01 9.280000e+02 4.400000e+01 4.400000e+01 +-2.625000e+00 1.250000e-01 1.250000e-01 9.690000e+02 4.500000e+01 4.500000e+01 +-2.375000e+00 1.250000e-01 1.250000e-01 9.980000e+02 4.500000e+01 4.500000e+01 +-2.125000e+00 1.250000e-01 1.250000e-01 1.022000e+03 4.400000e+01 4.400000e+01 +-1.875000e+00 1.250000e-01 1.250000e-01 1.048000e+03 4.600000e+01 4.600000e+01 +-1.625000e+00 1.250000e-01 1.250000e-01 1.062000e+03 4.700000e+01 4.700000e+01 +-1.375000e+00 1.250000e-01 1.250000e-01 1.064000e+03 4.700000e+01 4.700000e+01 +-1.125000e+00 1.250000e-01 1.250000e-01 1.026000e+03 4.600000e+01 4.600000e+01 +-8.750000e-01 1.250000e-01 1.250000e-01 1.014000e+03 2.800000e+01 2.800000e+01 +-6.250000e-01 1.250000e-01 1.250000e-01 1.006000e+03 2.400000e+01 2.400000e+01 +-3.750000e-01 1.250000e-01 1.250000e-01 9.910000e+02 2.400000e+01 2.400000e+01 +-1.250000e-01 1.250000e-01 1.250000e-01 9.820000e+02 2.400000e+01 2.400000e+01 +1.250000e-01 1.250000e-01 1.250000e-01 9.820000e+02 2.400000e+01 2.400000e+01 +3.750000e-01 1.250000e-01 1.250000e-01 9.910000e+02 2.400000e+01 2.400000e+01 +6.250000e-01 1.250000e-01 1.250000e-01 1.006000e+03 2.400000e+01 2.400000e+01 +8.750000e-01 1.250000e-01 1.250000e-01 1.014000e+03 2.800000e+01 2.800000e+01 +1.125000e+00 1.250000e-01 1.250000e-01 1.026000e+03 4.600000e+01 4.600000e+01 +1.375000e+00 1.250000e-01 1.250000e-01 1.064000e+03 4.700000e+01 4.700000e+01 +1.625000e+00 1.250000e-01 1.250000e-01 1.062000e+03 4.700000e+01 4.700000e+01 +1.875000e+00 1.250000e-01 1.250000e-01 1.048000e+03 4.600000e+01 4.600000e+01 +2.125000e+00 1.250000e-01 1.250000e-01 1.022000e+03 4.400000e+01 4.400000e+01 +2.375000e+00 1.250000e-01 1.250000e-01 9.980000e+02 4.500000e+01 4.500000e+01 +2.625000e+00 1.250000e-01 1.250000e-01 9.690000e+02 4.500000e+01 4.500000e+01 +2.875000e+00 1.250000e-01 1.250000e-01 9.280000e+02 4.400000e+01 4.400000e+01 +3.125000e+00 1.250000e-01 1.250000e-01 8.990000e+02 4.300000e+01 4.300000e+01 +3.375000e+00 1.250000e-01 1.250000e-01 8.620000e+02 4.200000e+01 4.200000e+01 +3.625000e+00 1.250000e-01 1.250000e-01 8.120000e+02 3.900000e+01 3.900000e+01 +3.875000e+00 1.250000e-01 1.250000e-01 7.540000e+02 3.600000e+01 3.600000e+01 +4.125000e+00 1.250000e-01 1.250000e-01 7.050000e+02 3.300000e+01 3.300000e+01 +4.375000e+00 1.250000e-01 1.250000e-01 6.560000e+02 3.200000e+01 3.200000e+01 +4.625000e+00 1.250000e-01 1.250000e-01 6.150000e+02 3.100000e+01 3.100000e+01 +4.875000e+00 1.250000e-01 1.250000e-01 5.590000e+02 2.900000e+01 2.900000e+01 +5.125000e+00 1.250000e-01 1.250000e-01 4.910000e+02 3.700000e+01 3.700000e+01 +5.375000e+00 1.250000e-01 1.250000e-01 4.510000e+02 2.900000e+01 2.900000e+01 +END YODA_SCATTER2D_V2 +BEGIN YODA_SCATTER2D_V2 /REF/ALICE_2013_I1225979/d01-x01-y04 +IsRef: 1 +Path: /REF/ALICE_2013_I1225979/d01-x01-y04 +Title: doi:10.17182/hepdata.68753.v1/t1 +Type: Scatter2D +--- +# xval xerr- xerr+ yval yerr- yerr+ +-4.875000e+00 1.250000e-01 1.250000e-01 3.890000e+02 2.100000e+01 2.100000e+01 +-4.625000e+00 1.250000e-01 1.250000e-01 4.290000e+02 2.200000e+01 2.200000e+01 +-4.375000e+00 1.250000e-01 1.250000e-01 4.590000e+02 2.400000e+01 2.400000e+01 +-4.125000e+00 1.250000e-01 1.250000e-01 4.910000e+02 2.500000e+01 2.500000e+01 +-3.875000e+00 1.250000e-01 1.250000e-01 5.260000e+02 2.700000e+01 2.700000e+01 +-3.625000e+00 1.250000e-01 1.250000e-01 5.630000e+02 3.000000e+01 3.000000e+01 +-3.375000e+00 1.250000e-01 1.250000e-01 5.960000e+02 3.100000e+01 3.100000e+01 +-3.125000e+00 1.250000e-01 1.250000e-01 6.190000e+02 3.200000e+01 3.200000e+01 +-2.875000e+00 1.250000e-01 1.250000e-01 6.410000e+02 3.400000e+01 3.400000e+01 +-2.625000e+00 1.250000e-01 1.250000e-01 6.690000e+02 3.500000e+01 3.500000e+01 +-2.375000e+00 1.250000e-01 1.250000e-01 6.830000e+02 3.300000e+01 3.300000e+01 +-2.125000e+00 1.250000e-01 1.250000e-01 7.000000e+02 3.400000e+01 3.400000e+01 +-1.875000e+00 1.250000e-01 1.250000e-01 7.160000e+02 3.500000e+01 3.500000e+01 +-1.625000e+00 1.250000e-01 1.250000e-01 7.250000e+02 3.600000e+01 3.600000e+01 +-1.375000e+00 1.250000e-01 1.250000e-01 7.240000e+02 3.700000e+01 3.700000e+01 +-1.125000e+00 1.250000e-01 1.250000e-01 7.010000e+02 3.400000e+01 3.400000e+01 +-8.750000e-01 1.250000e-01 1.250000e-01 6.920000e+02 1.900000e+01 1.900000e+01 +-6.250000e-01 1.250000e-01 1.250000e-01 6.820000e+02 1.600000e+01 1.600000e+01 +-3.750000e-01 1.250000e-01 1.250000e-01 6.730000e+02 1.600000e+01 1.600000e+01 +-1.250000e-01 1.250000e-01 1.250000e-01 6.660000e+02 1.600000e+01 1.600000e+01 +1.250000e-01 1.250000e-01 1.250000e-01 6.660000e+02 1.600000e+01 1.600000e+01 +3.750000e-01 1.250000e-01 1.250000e-01 6.730000e+02 1.600000e+01 1.600000e+01 +6.250000e-01 1.250000e-01 1.250000e-01 6.820000e+02 1.600000e+01 1.600000e+01 +8.750000e-01 1.250000e-01 1.250000e-01 6.920000e+02 1.900000e+01 1.900000e+01 +1.125000e+00 1.250000e-01 1.250000e-01 7.010000e+02 3.400000e+01 3.400000e+01 +1.375000e+00 1.250000e-01 1.250000e-01 7.240000e+02 3.700000e+01 3.700000e+01 +1.625000e+00 1.250000e-01 1.250000e-01 7.250000e+02 3.600000e+01 3.600000e+01 +1.875000e+00 1.250000e-01 1.250000e-01 7.160000e+02 3.500000e+01 3.500000e+01 +2.125000e+00 1.250000e-01 1.250000e-01 7.000000e+02 3.400000e+01 3.400000e+01 +2.375000e+00 1.250000e-01 1.250000e-01 6.830000e+02 3.300000e+01 3.300000e+01 +2.625000e+00 1.250000e-01 1.250000e-01 6.690000e+02 3.500000e+01 3.500000e+01 +2.875000e+00 1.250000e-01 1.250000e-01 6.410000e+02 3.400000e+01 3.400000e+01 +3.125000e+00 1.250000e-01 1.250000e-01 6.190000e+02 3.200000e+01 3.200000e+01 +3.375000e+00 1.250000e-01 1.250000e-01 5.960000e+02 3.100000e+01 3.100000e+01 +3.625000e+00 1.250000e-01 1.250000e-01 5.630000e+02 3.000000e+01 3.000000e+01 +3.875000e+00 1.250000e-01 1.250000e-01 5.260000e+02 2.700000e+01 2.700000e+01 +4.125000e+00 1.250000e-01 1.250000e-01 4.910000e+02 2.500000e+01 2.500000e+01 +4.375000e+00 1.250000e-01 1.250000e-01 4.590000e+02 2.400000e+01 2.400000e+01 +4.625000e+00 1.250000e-01 1.250000e-01 4.290000e+02 2.200000e+01 2.200000e+01 +4.875000e+00 1.250000e-01 1.250000e-01 3.890000e+02 2.100000e+01 2.100000e+01 +5.125000e+00 1.250000e-01 1.250000e-01 3.430000e+02 2.600000e+01 2.600000e+01 +5.375000e+00 1.250000e-01 1.250000e-01 3.130000e+02 2.200000e+01 2.200000e+01 +END YODA_SCATTER2D_V2 +BEGIN YODA_SCATTER3D_V2 /REF/ALICE_2013_I1225979/d02-x01-y01 +IsRef: 1 +Path: /REF/ALICE_2013_I1225979/d02-x01-y01 +Title: doi:10.17182/hepdata.68753.v1/t2 +Type: Scatter3D +--- +# xval xerr- xerr+ yval yerr- yerr+ zval zerr- zerr+ +2.500000e+00 2.500000e+00 2.500000e+00 3.828000e+02 3.100000e+00 3.100000e+00 1.716500e+04 7.722072e+02 7.722072e+02 +7.500000e+00 2.500000e+00 2.500000e+00 3.297000e+02 4.600000e+00 4.600000e+00 1.409900e+04 6.548191e+02 6.548191e+02 +1.500000e+01 5.000000e+00 5.000000e+00 2.605000e+02 4.400000e+00 4.400000e+00 1.058100e+04 5.353177e+02 5.353177e+02 +2.500000e+01 5.000000e+00 5.000000e+00 1.864000e+02 3.900000e+00 3.900000e+00 7.278000e+03 3.868036e+02 3.868036e+02 +END YODA_SCATTER3D_V2 +BEGIN YODA_SCATTER2D_V2 /REF/ALICE_2013_I1225979/d03-x01-y01 +IsRef: 1 +Path: /REF/ALICE_2013_I1225979/d03-x01-y01 +Title: doi:10.17182/hepdata.68753.v1/t3 +Type: Scatter2D +--- +# xval xerr- xerr+ yval yerr- yerr+ +1.864000e+02 3.900000e+00 3.900000e+00 7.510000e+00 1.500000e-01 1.500000e-01 +2.605000e+02 4.400000e+00 4.400000e+00 7.890000e+00 1.400000e-01 1.400000e-01 +3.297000e+02 4.600000e+00 4.600000e+00 8.340000e+00 1.400000e-01 1.400000e-01 +3.828000e+02 3.100000e+00 3.100000e+00 8.840000e+00 1.500000e-01 1.500000e-01 +END YODA_SCATTER2D_V2 +BEGIN YODA_SCATTER2D_V2 /REF/ALICE_2013_I1225979/d03-x01-y02 +IsRef: 1 +Path: /REF/ALICE_2013_I1225979/d03-x01-y02 +Title: doi:10.17182/hepdata.68753.v1/t3 +Type: Scatter2D +--- +# xval xerr- xerr+ yval yerr- yerr+ +1.864000e+02 3.900000e+00 3.900000e+00 7.580000e+00 1.800000e-01 1.800000e-01 +2.605000e+02 4.400000e+00 4.400000e+00 7.930000e+00 1.700000e-01 1.700000e-01 +3.297000e+02 4.600000e+00 4.600000e+00 8.350000e+00 1.700000e-01 1.700000e-01 +3.828000e+02 3.100000e+00 3.100000e+00 8.810000e+00 1.700000e-01 1.700000e-01 +END YODA_SCATTER2D_V2 +BEGIN YODA_SCATTER2D_V2 /REF/ALICE_2013_I1225979/d03-x01-y03 +IsRef: 1 +Path: /REF/ALICE_2013_I1225979/d03-x01-y03 +Title: doi:10.17182/hepdata.68753.v1/t3 +Type: Scatter2D +--- +# xval xerr- xerr+ yval yerr- yerr+ +1.864000e+02 3.900000e+00 3.900000e+00 6.770000e+00 1.800000e-01 1.800000e-01 +2.605000e+02 4.400000e+00 4.400000e+00 7.020000e+00 1.700000e-01 1.700000e-01 +3.297000e+02 4.600000e+00 4.600000e+00 7.400000e+00 1.500000e-01 1.500000e-01 +3.828000e+02 3.100000e+00 3.100000e+00 7.730000e+00 1.500000e-01 1.500000e-01 +END YODA_SCATTER2D_V2 +BEGIN YODA_SCATTER2D_V2 /REF/ALICE_2013_I1225979/d03-x01-y04 +IsRef: 1 +Path: /REF/ALICE_2013_I1225979/d03-x01-y04 +Title: doi:10.17182/hepdata.68753.v1/t3 +Type: Scatter2D +--- +# xval xerr- xerr+ yval yerr- yerr+ +1.864000e+02 3.900000e+00 3.900000e+00 5.470000e+00 1.400000e-01 1.400000e-01 +2.605000e+02 4.400000e+00 4.400000e+00 5.620000e+00 1.400000e-01 1.400000e-01 +3.297000e+02 4.600000e+00 4.600000e+00 5.860000e+00 1.300000e-01 1.300000e-01 +3.828000e+02 3.100000e+00 3.100000e+00 6.120000e+00 1.300000e-01 1.300000e-01 +END YODA_SCATTER2D_V2 +BEGIN YODA_SCATTER2D_V2 /REF/ALICE_2013_I1225979/d03-x01-y05 +IsRef: 1 +Path: /REF/ALICE_2013_I1225979/d03-x01-y05 +Title: doi:10.17182/hepdata.68753.v1/t3 +Type: Scatter2D +--- +# xval xerr- xerr+ yval yerr- yerr+ +1.864000e+02 3.900000e+00 3.900000e+00 3.950000e+00 1.200000e-01 1.200000e-01 +2.605000e+02 4.400000e+00 4.400000e+00 4.060000e+00 1.200000e-01 1.200000e-01 +3.297000e+02 4.600000e+00 4.600000e+00 4.210000e+00 1.100000e-01 1.100000e-01 +3.828000e+02 3.100000e+00 3.100000e+00 4.310000e+00 1.100000e-01 1.100000e-01 +END YODA_SCATTER2D_V2 +BEGIN YODA_SCATTER2D_V2 /REF/ALICE_2013_I1225979/d04-x01-y01 +IsRef: 1 +Path: /REF/ALICE_2013_I1225979/d04-x01-y01 +Title: doi:10.17182/hepdata.68753.v1/t4 +Type: Scatter2D +--- +# xval xerr- xerr+ yval yerr- yerr+ +-4.694000e+00 2.020000e-01 2.020000e-01 8.870000e+02 3.900000e+01 3.900000e+01 +-4.444000e+00 1.940000e-01 1.940000e-01 9.780000e+02 4.100000e+01 4.100000e+01 +-4.194000e+00 1.870000e-01 1.870000e-01 1.047000e+03 4.400000e+01 4.400000e+01 +-3.945000e+00 1.810000e-01 1.810000e-01 1.132000e+03 4.600000e+01 4.600000e+01 +-3.695000e+00 1.740000e-01 1.740000e-01 1.212000e+03 5.000000e+01 5.000000e+01 +-3.445000e+00 1.680000e-01 1.680000e-01 1.307000e+03 5.400000e+01 5.400000e+01 +-3.196000e+00 1.630000e-01 1.630000e-01 1.393000e+03 5.800000e+01 5.800000e+01 +-2.946000e+00 1.570000e-01 1.570000e-01 1.455000e+03 6.000000e+01 6.000000e+01 +-2.698000e+00 1.520000e-01 1.520000e-01 1.514000e+03 5.900000e+01 5.900000e+01 +-2.450000e+00 1.470000e-01 1.470000e-01 1.594000e+03 6.200000e+01 6.200000e+01 +-2.202000e+00 1.400000e-01 1.400000e-01 1.653000e+03 6.400000e+01 6.400000e+01 +-1.956000e+00 1.330000e-01 1.330000e-01 1.708000e+03 6.500000e+01 6.500000e+01 +-1.712000e+00 1.270000e-01 1.270000e-01 1.761000e+03 7.100000e+01 7.100000e+01 +-1.471000e+00 1.220000e-01 1.220000e-01 1.813000e+03 7.100000e+01 7.100000e+01 +-1.232000e+00 1.170000e-01 1.170000e-01 1.843000e+03 7.500000e+01 7.500000e+01 +-9.980000e-01 1.140000e-01 1.140000e-01 1.831000e+03 8.000000e+01 8.000000e+01 +-7.690000e-01 1.110000e-01 1.110000e-01 1.844000e+03 6.100000e+01 6.100000e+01 +-5.440000e-01 1.090000e-01 1.090000e-01 1.867000e+03 5.700000e+01 5.700000e+01 +-3.240000e-01 1.080000e-01 1.080000e-01 1.866000e+03 5.600000e+01 5.600000e+01 +-1.080000e-01 1.080000e-01 1.080000e-01 1.873000e+03 5.400000e+01 5.400000e+01 +1.080000e-01 1.080000e-01 1.080000e-01 1.873000e+03 5.400000e+01 5.400000e+01 +3.240000e-01 1.080000e-01 1.080000e-01 1.866000e+03 5.600000e+01 5.600000e+01 +5.440000e-01 1.090000e-01 1.090000e-01 1.867000e+03 5.700000e+01 5.700000e+01 +7.690000e-01 1.110000e-01 1.110000e-01 1.844000e+03 6.100000e+01 6.100000e+01 +9.980000e-01 1.140000e-01 1.140000e-01 1.831000e+03 8.000000e+01 8.000000e+01 +1.232000e+00 1.170000e-01 1.170000e-01 1.843000e+03 7.500000e+01 7.500000e+01 +1.471000e+00 1.220000e-01 1.220000e-01 1.813000e+03 7.100000e+01 7.100000e+01 +1.712000e+00 1.270000e-01 1.270000e-01 1.761000e+03 7.100000e+01 7.100000e+01 +1.956000e+00 1.330000e-01 1.330000e-01 1.708000e+03 6.500000e+01 6.500000e+01 +2.202000e+00 1.400000e-01 1.400000e-01 1.653000e+03 6.400000e+01 6.400000e+01 +2.450000e+00 1.470000e-01 1.470000e-01 1.594000e+03 6.200000e+01 6.200000e+01 +2.698000e+00 1.520000e-01 1.520000e-01 1.514000e+03 5.900000e+01 5.900000e+01 +2.946000e+00 1.570000e-01 1.570000e-01 1.455000e+03 6.000000e+01 6.000000e+01 +3.196000e+00 1.630000e-01 1.630000e-01 1.393000e+03 5.800000e+01 5.800000e+01 +3.445000e+00 1.680000e-01 1.680000e-01 1.307000e+03 5.400000e+01 5.400000e+01 +3.695000e+00 1.740000e-01 1.740000e-01 1.212000e+03 5.000000e+01 5.000000e+01 +3.945000e+00 1.810000e-01 1.810000e-01 1.132000e+03 4.600000e+01 4.600000e+01 +4.194000e+00 1.870000e-01 1.870000e-01 1.047000e+03 4.400000e+01 4.400000e+01 +4.444000e+00 1.940000e-01 1.940000e-01 9.780000e+02 4.100000e+01 4.100000e+01 +4.694000e+00 2.020000e-01 2.020000e-01 8.870000e+02 3.900000e+01 3.900000e+01 +4.944000e+00 2.100000e-01 2.100000e-01 7.700000e+02 4.200000e+01 4.200000e+01 +5.194000e+00 2.180000e-01 2.180000e-01 6.670000e+02 3.900000e+01 3.900000e+01 +END YODA_SCATTER2D_V2 +BEGIN YODA_SCATTER2D_V2 /REF/ALICE_2013_I1225979/d05-x01-y01 +IsRef: 1 +Path: /REF/ALICE_2013_I1225979/d05-x01-y01 +Title: doi:10.17182/hepdata.68753.v1/t5 +Type: Scatter2D +--- +# xval xerr- xerr+ yval yerr- yerr+ +1.000000e+00 5.000000e-01 5.000000e-01 1.250000e+00 2.000000e-02 2.000000e-02 +END YODA_SCATTER2D_V2 +BEGIN YODA_SCATTER2D_V2 /REF/ALICE_2013_I1225979/d05-x01-y02 +IsRef: 1 +Path: /REF/ALICE_2013_I1225979/d05-x01-y02 +Title: doi:10.17182/hepdata.68753.v1/t5 +Type: Scatter2D +--- +# xval xerr- xerr+ yval yerr- yerr+ +1.000000e+00 5.000000e-01 5.000000e-01 1.180000e+00 2.000000e-02 2.000000e-02 +END YODA_SCATTER2D_V2 +BEGIN YODA_SCATTER2D_V2 /REF/ALICE_2013_I1225979/d05-x01-y03 +IsRef: 1 +Path: /REF/ALICE_2013_I1225979/d05-x01-y03 +Title: doi:10.17182/hepdata.68753.v1/t5 +Type: Scatter2D +--- +# xval xerr- xerr+ yval yerr- yerr+ +1.000000e+00 5.000000e-01 5.000000e-01 1.430000e+00 2.000000e-02 2.000000e-02 +END YODA_SCATTER2D_V2 diff --git a/analyses/pluginALICE/ALICE_2015_PBPBCentrality.cc b/analyses/pluginALICE/ALICE_2015_PBPBCentrality.cc new file mode 100644 --- /dev/null +++ b/analyses/pluginALICE/ALICE_2015_PBPBCentrality.cc @@ -0,0 +1,82 @@ +/** + * @file ALICE_2015_PBPBCentrality.cc + * @author Christian Holm Christensen + * @date Wed Aug 22 15:56:23 2018 + * + * @brief Dummy analysis for centrality calibration in Pb-Pb at 5.02TeV + */ +#include +#include + +namespace Rivet +{ + /** + * Dummy analysis for centrality calibration in Pb-Pb at 5.02TeV + */ + class ALICE_2015_PBPBCentrality : public Analysis + { + public: + /** + * Constructor + */ + ALICE_2015_PBPBCentrality() + : Analysis("ALICE_2015_PBPBCentrality") + { + } + /** + * Initialize this analysis. + */ + void init() + { + ALICE::V0AndTrigger v0and; + declare(v0and,"V0-AND"); + + ALICE::V0MMultiplicity v0m; + declare(v0m,"V0M"); + + _v0m = bookHisto1D("V0M","Forward multiplicity","V0M","Events"); + _imp = bookHisto1D("V0M_IMP",100,0,20, + "Impact parameter","b (fm)","Events"); + } + /** + * Analyse a single event. + * + * @param event The event + */ + void analyze(const Event& event) + { + // Get and fill in the impact parameter value if the information + // is valid. + const HepMC::GenEvent* ge = event.genEvent(); + const HepMC::HeavyIon* hi = ge->heavy_ion(); + if (hi && hi->is_valid()) + _imp->fill(hi->impact_parameter(), event.weight()); + + + // Check if we have any hit in either V0-A or -C. If not, the + // event is not selected and we get out. + if (!apply(event,"V0-AND")()) return; + + // Fill in the V0 multiplicity for this event + _v0m->fill(apply(event,"V0M")(), event.weight()); + } + /** + * Finalize this analysis + */ + void finalize() + { + _v0m->normalize(); + _imp->normalize(); + } + + /** The distribution of V0M multiplicity */ + Histo1DPtr _v0m; + /** The distribution of impact parameters */ + Histo1DPtr _imp; + }; + // The hook for the plugin system + DECLARE_RIVET_PLUGIN(ALICE_2015_PBPBCentrality); +} +// +// EOF +// diff --git a/analyses/pluginALICE/ALICE_2015_PBPBCentrality.info b/analyses/pluginALICE/ALICE_2015_PBPBCentrality.info new file mode 100644 --- /dev/null +++ b/analyses/pluginALICE/ALICE_2015_PBPBCentrality.info @@ -0,0 +1,20 @@ +# -*- mode: conf-colon-mode -*- +Name: ALICE_2015_PBPBCentrality +Year: 2018 +Summary: ALICE centrality calibration for Pb-Pb at 5.02TeV +Experiment: ALICE +Collider: LHC +Status: NEW +Authors: + - Christian Holm Christensen +RunInfo: Pb-Pb at sqrt(sNN)=5.02TeV +NumEvents: 1000000 +NeedCrossSection: no +Beams: [1000822080, 1000822080] +# This _total_ energy of beams, so this becomes 208*5023=1044784 +Energies: [1044784] +Description: + 'Dummy analysis to bring in the centrality calibration for Pb-Pb at 5.02TeV' +ToDo: + 'Validate' + diff --git a/analyses/pluginALICE/ALICE_2015_PBPBCentrality.plot b/analyses/pluginALICE/ALICE_2015_PBPBCentrality.plot new file mode 100644 --- /dev/null +++ b/analyses/pluginALICE/ALICE_2015_PBPBCentrality.plot @@ -0,0 +1,11 @@ +BEGIN PLOT /ALICE_2015_PBPBCentrality/V0M + Title=Forward multiplicity + XLabel=V0M + YLabel=Events +END PLOT + +BEGIN PLOT /ALICE_2015_PBPBCentrality/V0M_IMP + Title=Impact parameter + XLabel=$b$ (fm) + YLabel=Events +END PLOT diff --git a/analyses/pluginALICE/ALICE_2015_PBPBCentrality.yoda b/analyses/pluginALICE/ALICE_2015_PBPBCentrality.yoda new file mode 100644 --- /dev/null +++ b/analyses/pluginALICE/ALICE_2015_PBPBCentrality.yoda @@ -0,0 +1,508 @@ +BEGIN YODA_SCATTER2D_V2 /REF/ALICE_2015_PBPBCentrality/V0M +Variations: [""] +Path: /REF/ALICE_2015_PBPBCentrality/V0M +Title: V0M +Type: Scatter2D +--- +# xval xerr- xerr+ yval yerr- yerr- +3.500000e+01 3.500000e+01 4.000000e+01 7.685643e-01 1.918994e-03 1.918994e-03 +1.150000e+02 4.000000e+01 4.000000e+01 1.565472e-02 4.189411e-04 4.189411e-04 +1.950000e+02 4.000000e+01 4.000000e+01 9.106134e-03 3.211281e-04 3.211281e-04 +2.750000e+02 4.000000e+01 4.000000e+01 7.007011e-03 2.817957e-04 2.817957e-04 +3.550000e+02 4.000000e+01 4.000000e+01 5.706946e-03 2.543899e-04 2.543899e-04 +4.350000e+02 4.000000e+01 4.000000e+01 4.885180e-03 2.353170e-04 2.353170e-04 +5.150000e+02 4.000000e+01 4.000000e+01 4.101734e-03 2.157543e-04 2.157543e-04 +5.950000e+02 4.000000e+01 4.000000e+01 3.826393e-03 2.083850e-04 2.083850e-04 +6.750000e+02 4.000000e+01 4.000000e+01 3.333901e-03 1.945097e-04 1.945097e-04 +7.550000e+02 4.000000e+01 4.000000e+01 2.993272e-03 1.842787e-04 1.842787e-04 +8.350000e+02 4.000000e+01 4.000000e+01 2.786057e-03 1.777550e-04 1.777550e-04 +9.150000e+02 4.000000e+01 4.000000e+01 2.598711e-03 1.716658e-04 1.716658e-04 +9.950000e+02 4.000000e+01 4.000000e+01 2.391496e-03 1.646705e-04 1.646705e-04 +1.075000e+03 4.000000e+01 4.000000e+01 2.336143e-03 1.627642e-04 1.627642e-04 +1.155000e+03 4.000000e+01 4.000000e+01 2.155894e-03 1.563456e-04 1.563456e-04 +1.235000e+03 4.000000e+01 4.000000e+01 2.028159e-03 1.517183e-04 1.517183e-04 +1.315000e+03 4.000000e+01 4.000000e+01 1.996934e-03 1.504815e-04 1.504815e-04 +1.395000e+03 4.000000e+01 4.000000e+01 1.891907e-03 1.464720e-04 1.464720e-04 +1.475000e+03 4.000000e+01 4.000000e+01 1.747140e-03 1.408029e-04 1.408029e-04 +1.555000e+03 4.000000e+01 4.000000e+01 1.717335e-03 1.394954e-04 1.394954e-04 +1.635000e+03 4.000000e+01 4.000000e+01 1.568310e-03 1.333549e-04 1.333549e-04 +1.715000e+03 4.000000e+01 4.000000e+01 1.527151e-03 1.316365e-04 1.316365e-04 +1.795000e+03 4.000000e+01 4.000000e+01 1.497346e-03 1.302743e-04 1.302743e-04 +1.875000e+03 4.000000e+01 4.000000e+01 1.510119e-03 1.309109e-04 1.309109e-04 +1.955000e+03 4.000000e+01 4.000000e+01 1.415027e-03 1.264448e-04 1.264448e-04 +2.035000e+03 4.000000e+01 4.000000e+01 1.382384e-03 1.252666e-04 1.252666e-04 +2.115000e+03 4.000000e+01 4.000000e+01 1.317097e-03 1.221327e-04 1.221327e-04 +2.195000e+03 4.000000e+01 4.000000e+01 1.300065e-03 1.214522e-04 1.214522e-04 +2.275000e+03 4.000000e+01 4.000000e+01 1.253229e-03 1.190653e-04 1.190653e-04 +2.355000e+03 4.000000e+01 4.000000e+01 1.247552e-03 1.189358e-04 1.189358e-04 +2.435000e+03 4.000000e+01 4.000000e+01 1.219166e-03 1.175022e-04 1.175022e-04 +2.515000e+03 4.000000e+01 4.000000e+01 1.158137e-03 1.144672e-04 1.144672e-04 +2.595000e+03 4.000000e+01 4.000000e+01 1.179426e-03 1.154510e-04 1.154510e-04 +2.675000e+03 4.000000e+01 4.000000e+01 1.172330e-03 1.152779e-04 1.152779e-04 +2.755000e+03 4.000000e+01 4.000000e+01 1.074399e-03 1.103092e-04 1.103092e-04 +2.835000e+03 4.000000e+01 4.000000e+01 9.991768e-04 1.062779e-04 1.062779e-04 +2.915000e+03 4.000000e+01 4.000000e+01 1.031820e-03 1.081240e-04 1.081240e-04 +2.995000e+03 4.000000e+01 4.000000e+01 1.060206e-03 1.096694e-04 1.096694e-04 +3.075000e+03 4.000000e+01 4.000000e+01 9.934996e-04 1.059795e-04 1.059795e-04 +3.155000e+03 4.000000e+01 4.000000e+01 9.807261e-04 1.053959e-04 1.053959e-04 +3.235000e+03 4.000000e+01 4.000000e+01 1.053110e-03 1.091852e-04 1.091852e-04 +3.315000e+03 4.000000e+01 4.000000e+01 9.991768e-04 1.063075e-04 1.063075e-04 +3.395000e+03 4.000000e+01 4.000000e+01 9.523403e-04 1.039303e-04 1.039303e-04 +3.475000e+03 4.000000e+01 4.000000e+01 9.636946e-04 1.045466e-04 1.045466e-04 +3.555000e+03 4.000000e+01 4.000000e+01 9.736296e-04 1.050719e-04 1.050719e-04 +3.635000e+03 4.000000e+01 4.000000e+01 9.182775e-04 1.019979e-04 1.019979e-04 +3.715000e+03 4.000000e+01 4.000000e+01 8.870532e-04 1.000890e-04 1.000890e-04 +3.795000e+03 4.000000e+01 4.000000e+01 9.097618e-04 1.014717e-04 1.014717e-04 +3.875000e+03 4.000000e+01 4.000000e+01 9.267932e-04 1.022912e-04 1.022912e-04 +3.955000e+03 4.000000e+01 4.000000e+01 8.586675e-04 9.863270e-05 9.863270e-05 +4.035000e+03 4.000000e+01 4.000000e+01 8.274433e-04 9.685953e-05 9.685953e-05 +4.115000e+03 4.000000e+01 4.000000e+01 8.288625e-04 9.683860e-05 9.683860e-05 +4.195000e+03 4.000000e+01 4.000000e+01 8.345397e-04 9.726063e-05 9.726063e-05 +4.275000e+03 4.000000e+01 4.000000e+01 8.033154e-04 9.543075e-05 9.543075e-05 +4.355000e+03 4.000000e+01 4.000000e+01 7.493826e-04 9.215666e-05 9.215666e-05 +4.435000e+03 4.000000e+01 4.000000e+01 7.791876e-04 9.384564e-05 9.384564e-05 +4.515000e+03 4.000000e+01 4.000000e+01 7.408669e-04 9.159535e-05 9.159535e-05 +4.595000e+03 4.000000e+01 4.000000e+01 7.962190e-04 9.499222e-05 9.499222e-05 +4.675000e+03 4.000000e+01 4.000000e+01 7.337704e-04 9.121345e-05 9.121345e-05 +4.755000e+03 4.000000e+01 4.000000e+01 7.877033e-04 9.429941e-05 9.429941e-05 +4.835000e+03 4.000000e+01 4.000000e+01 7.011269e-04 8.904765e-05 8.904765e-05 +4.915000e+03 4.000000e+01 4.000000e+01 7.351897e-04 9.131073e-05 9.131073e-05 +4.995000e+03 4.000000e+01 4.000000e+01 7.578983e-04 9.261763e-05 9.261763e-05 +5.075000e+03 4.000000e+01 4.000000e+01 7.295126e-04 9.095511e-05 9.095511e-05 +5.155000e+03 4.000000e+01 4.000000e+01 7.933804e-04 9.467316e-05 9.467316e-05 +5.235000e+03 4.000000e+01 4.000000e+01 7.578983e-04 9.246926e-05 9.246926e-05 +5.315000e+03 4.000000e+01 4.000000e+01 6.514519e-04 8.573586e-05 8.573586e-05 +5.395000e+03 4.000000e+01 4.000000e+01 6.997076e-04 8.887232e-05 8.887232e-05 +5.475000e+03 4.000000e+01 4.000000e+01 6.997076e-04 8.903380e-05 8.903380e-05 +5.555000e+03 4.000000e+01 4.000000e+01 6.642255e-04 8.671529e-05 8.671529e-05 +5.635000e+03 4.000000e+01 4.000000e+01 6.855147e-04 8.819625e-05 8.819625e-05 +5.715000e+03 4.000000e+01 4.000000e+01 6.315819e-04 8.448270e-05 8.448270e-05 +5.795000e+03 4.000000e+01 4.000000e+01 5.989383e-04 8.221161e-05 8.221161e-05 +5.875000e+03 4.000000e+01 4.000000e+01 6.670640e-04 8.671868e-05 8.671868e-05 +5.955000e+03 4.000000e+01 4.000000e+01 7.451247e-04 9.184278e-05 9.184278e-05 +6.035000e+03 4.000000e+01 4.000000e+01 6.443555e-04 8.533151e-05 8.533151e-05 +6.115000e+03 4.000000e+01 4.000000e+01 6.046155e-04 8.280223e-05 8.280223e-05 +6.195000e+03 4.000000e+01 4.000000e+01 6.159698e-04 8.354831e-05 8.354831e-05 +6.275000e+03 4.000000e+01 4.000000e+01 6.031962e-04 8.259599e-05 8.259599e-05 +6.355000e+03 4.000000e+01 4.000000e+01 6.486133e-04 8.567241e-05 8.567241e-05 +6.435000e+03 4.000000e+01 4.000000e+01 6.670640e-04 8.691943e-05 8.691943e-05 +6.515000e+03 4.000000e+01 4.000000e+01 6.003576e-04 8.244357e-05 8.244357e-05 +6.595000e+03 4.000000e+01 4.000000e+01 5.989383e-04 8.234529e-05 8.234529e-05 +6.675000e+03 4.000000e+01 4.000000e+01 5.819069e-04 8.109768e-05 8.109768e-05 +6.755000e+03 4.000000e+01 4.000000e+01 6.443555e-04 8.536594e-05 8.536594e-05 +6.835000e+03 4.000000e+01 4.000000e+01 6.188083e-04 8.354397e-05 8.354397e-05 +6.915000e+03 4.000000e+01 4.000000e+01 6.386783e-04 8.493479e-05 8.493479e-05 +6.995000e+03 4.000000e+01 4.000000e+01 5.989383e-04 8.238760e-05 8.238760e-05 +7.075000e+03 4.000000e+01 4.000000e+01 6.003576e-04 8.246909e-05 8.246909e-05 +7.155000e+03 4.000000e+01 4.000000e+01 6.017769e-04 8.242978e-05 8.242978e-05 +7.235000e+03 4.000000e+01 4.000000e+01 5.890033e-04 8.153480e-05 8.153480e-05 +7.315000e+03 4.000000e+01 4.000000e+01 5.577791e-04 7.950831e-05 7.950831e-05 +7.395000e+03 4.000000e+01 4.000000e+01 5.719719e-04 8.045131e-05 8.045131e-05 +7.475000e+03 4.000000e+01 4.000000e+01 5.379091e-04 7.808159e-05 7.808159e-05 +7.555000e+03 4.000000e+01 4.000000e+01 5.393284e-04 7.800198e-05 7.800198e-05 +7.635000e+03 4.000000e+01 4.000000e+01 5.450055e-04 7.843271e-05 7.843271e-05 +7.715000e+03 4.000000e+01 4.000000e+01 5.506826e-04 7.890894e-05 7.890894e-05 +7.795000e+03 4.000000e+01 4.000000e+01 5.364898e-04 7.779588e-05 7.779588e-05 +7.875000e+03 4.000000e+01 4.000000e+01 5.166198e-04 7.644255e-05 7.644255e-05 +7.955000e+03 4.000000e+01 4.000000e+01 5.024269e-04 7.540603e-05 7.540603e-05 +8.035000e+03 4.000000e+01 4.000000e+01 5.379091e-04 7.793945e-05 7.793945e-05 +8.115000e+03 4.000000e+01 4.000000e+01 4.825570e-04 7.376986e-05 7.376986e-05 +8.195000e+03 4.000000e+01 4.000000e+01 5.308127e-04 7.743073e-05 7.743073e-05 +8.275000e+03 4.000000e+01 4.000000e+01 5.308127e-04 7.740862e-05 7.740862e-05 +8.355000e+03 4.000000e+01 4.000000e+01 5.350705e-04 7.785405e-05 7.785405e-05 +8.435000e+03 4.000000e+01 4.000000e+01 4.754605e-04 7.333535e-05 7.333535e-05 +8.515000e+03 4.000000e+01 4.000000e+01 4.683641e-04 7.272903e-05 7.272903e-05 +8.595000e+03 4.000000e+01 4.000000e+01 4.782991e-04 7.358791e-05 7.358791e-05 +8.675000e+03 4.000000e+01 4.000000e+01 4.967498e-04 7.480694e-05 7.480694e-05 +8.755000e+03 4.000000e+01 4.000000e+01 5.251355e-04 7.700405e-05 7.700405e-05 +8.835000e+03 4.000000e+01 4.000000e+01 5.137812e-04 7.594366e-05 7.594366e-05 +8.915000e+03 4.000000e+01 4.000000e+01 4.669448e-04 7.259666e-05 7.259666e-05 +8.995000e+03 4.000000e+01 4.000000e+01 4.953305e-04 7.468629e-05 7.468629e-05 +9.075000e+03 4.000000e+01 4.000000e+01 4.385591e-04 7.027335e-05 7.027335e-05 +9.155000e+03 4.000000e+01 4.000000e+01 4.924920e-04 7.469025e-05 7.469025e-05 +9.235000e+03 4.000000e+01 4.000000e+01 4.413977e-04 7.060201e-05 7.060201e-05 +9.315000e+03 4.000000e+01 4.000000e+01 4.726220e-04 7.295060e-05 7.295060e-05 +9.395000e+03 4.000000e+01 4.000000e+01 4.924920e-04 7.450889e-05 7.450889e-05 +9.475000e+03 4.000000e+01 4.000000e+01 4.272048e-04 6.944204e-05 6.944204e-05 +9.555000e+03 4.000000e+01 4.000000e+01 4.584291e-04 7.199381e-05 7.199381e-05 +9.635000e+03 4.000000e+01 4.000000e+01 4.456555e-04 7.073106e-05 7.073106e-05 +9.715000e+03 4.000000e+01 4.000000e+01 4.910727e-04 7.463518e-05 7.463518e-05 +9.795000e+03 4.000000e+01 4.000000e+01 4.882341e-04 7.437876e-05 7.437876e-05 +9.875000e+03 4.000000e+01 4.000000e+01 4.413977e-04 7.017109e-05 7.017109e-05 +9.955000e+03 4.000000e+01 4.000000e+01 4.484941e-04 7.125967e-05 7.125967e-05 +1.003500e+04 4.000000e+01 4.000000e+01 3.888841e-04 6.622162e-05 6.622162e-05 +1.011500e+04 4.000000e+01 4.000000e+01 4.087541e-04 6.783051e-05 6.783051e-05 +1.019500e+04 4.000000e+01 4.000000e+01 4.229470e-04 6.915721e-05 6.915721e-05 +1.027500e+04 4.000000e+01 4.000000e+01 4.655255e-04 7.242896e-05 7.242896e-05 +1.035500e+04 4.000000e+01 4.000000e+01 4.385591e-04 7.036485e-05 7.036485e-05 +1.043500e+04 4.000000e+01 4.000000e+01 4.456555e-04 7.098055e-05 7.098055e-05 +1.051500e+04 4.000000e+01 4.000000e+01 4.243663e-04 6.933172e-05 6.933172e-05 +1.059500e+04 4.000000e+01 4.000000e+01 3.903034e-04 6.645454e-05 6.645454e-05 +1.067500e+04 4.000000e+01 4.000000e+01 4.101734e-04 6.805903e-05 6.805903e-05 +1.075500e+04 4.000000e+01 4.000000e+01 4.470748e-04 7.113272e-05 7.113272e-05 +1.083500e+04 4.000000e+01 4.000000e+01 4.385591e-04 7.043288e-05 7.043288e-05 +1.091500e+04 4.000000e+01 4.000000e+01 4.172698e-04 6.844703e-05 6.844703e-05 +1.099500e+04 4.000000e+01 4.000000e+01 4.059156e-04 6.766259e-05 6.766259e-05 +1.107500e+04 4.000000e+01 4.000000e+01 4.130120e-04 6.827501e-05 6.827501e-05 +1.115500e+04 4.000000e+01 4.000000e+01 3.761106e-04 6.503908e-05 6.503908e-05 +1.123500e+04 4.000000e+01 4.000000e+01 3.903034e-04 6.632548e-05 6.632548e-05 +1.131500e+04 4.000000e+01 4.000000e+01 4.328820e-04 6.992365e-05 6.992365e-05 +1.139500e+04 4.000000e+01 4.000000e+01 4.002384e-04 6.719360e-05 6.719360e-05 +1.147500e+04 4.000000e+01 4.000000e+01 4.016577e-04 6.744339e-05 6.744339e-05 +1.155500e+04 4.000000e+01 4.000000e+01 4.115927e-04 6.797943e-05 6.797943e-05 +1.163500e+04 4.000000e+01 4.000000e+01 4.115927e-04 6.782516e-05 6.782516e-05 +1.171500e+04 4.000000e+01 4.000000e+01 4.201084e-04 6.897960e-05 6.897960e-05 +1.179500e+04 4.000000e+01 4.000000e+01 4.385591e-04 7.038483e-05 7.038483e-05 +1.187500e+04 4.000000e+01 4.000000e+01 4.172698e-04 6.854046e-05 6.854046e-05 +1.195500e+04 4.000000e+01 4.000000e+01 3.604984e-04 6.379710e-05 6.379710e-05 +1.203500e+04 4.000000e+01 4.000000e+01 3.817877e-04 6.546254e-05 6.546254e-05 +1.211500e+04 4.000000e+01 4.000000e+01 3.846263e-04 6.569903e-05 6.569903e-05 +1.219500e+04 4.000000e+01 4.000000e+01 4.030770e-04 6.732386e-05 6.732386e-05 +1.227500e+04 4.000000e+01 4.000000e+01 3.718527e-04 6.493605e-05 6.493605e-05 +1.235500e+04 4.000000e+01 4.000000e+01 3.846263e-04 6.592789e-05 6.592789e-05 +1.243500e+04 4.000000e+01 4.000000e+01 3.534020e-04 6.321706e-05 6.321706e-05 +1.251500e+04 4.000000e+01 4.000000e+01 3.633370e-04 6.398286e-05 6.398286e-05 +1.259500e+04 4.000000e+01 4.000000e+01 3.505634e-04 6.294789e-05 6.294789e-05 +1.267500e+04 4.000000e+01 4.000000e+01 3.917227e-04 6.658293e-05 6.658293e-05 +1.275500e+04 4.000000e+01 4.000000e+01 3.448863e-04 6.251679e-05 6.251679e-05 +1.283500e+04 4.000000e+01 4.000000e+01 3.448863e-04 6.247813e-05 6.247813e-05 +1.291500e+04 4.000000e+01 4.000000e+01 3.647563e-04 6.416131e-05 6.416131e-05 +1.299500e+04 4.000000e+01 4.000000e+01 3.789491e-04 6.543936e-05 6.543936e-05 +1.307500e+04 4.000000e+01 4.000000e+01 3.420477e-04 6.217087e-05 6.217087e-05 +1.315500e+04 4.000000e+01 4.000000e+01 3.789491e-04 6.527152e-05 6.527152e-05 +1.323500e+04 4.000000e+01 4.000000e+01 3.746913e-04 6.475468e-05 6.475468e-05 +1.331500e+04 4.000000e+01 4.000000e+01 3.860456e-04 6.591232e-05 6.591232e-05 +1.339500e+04 4.000000e+01 4.000000e+01 3.874649e-04 6.616330e-05 6.616330e-05 +1.347500e+04 4.000000e+01 4.000000e+01 3.704334e-04 6.474444e-05 6.474444e-05 +1.355500e+04 4.000000e+01 4.000000e+01 3.633370e-04 6.415344e-05 6.415344e-05 +1.363500e+04 4.000000e+01 4.000000e+01 3.860456e-04 6.561053e-05 6.561053e-05 +1.371500e+04 4.000000e+01 4.000000e+01 3.604984e-04 6.382324e-05 6.382324e-05 +1.379500e+04 4.000000e+01 4.000000e+01 3.477249e-04 6.263385e-05 6.263385e-05 +1.387500e+04 4.000000e+01 4.000000e+01 3.675949e-04 6.413710e-05 6.413710e-05 +1.395500e+04 4.000000e+01 4.000000e+01 3.534020e-04 6.318940e-05 6.318940e-05 +1.403500e+04 4.000000e+01 4.000000e+01 3.448863e-04 6.249412e-05 6.249412e-05 +1.411500e+04 4.000000e+01 4.000000e+01 3.420477e-04 6.221383e-05 6.221383e-05 +1.419500e+04 4.000000e+01 4.000000e+01 3.278549e-04 6.076007e-05 6.076007e-05 +1.427500e+04 4.000000e+01 4.000000e+01 3.207584e-04 6.004780e-05 6.004780e-05 +1.435500e+04 4.000000e+01 4.000000e+01 3.079849e-04 5.901938e-05 5.901938e-05 +1.443500e+04 4.000000e+01 4.000000e+01 3.150813e-04 5.970140e-05 5.970140e-05 +1.451500e+04 4.000000e+01 4.000000e+01 3.278549e-04 6.089889e-05 6.089889e-05 +1.459500e+04 4.000000e+01 4.000000e+01 3.321127e-04 6.122016e-05 6.122016e-05 +1.467500e+04 4.000000e+01 4.000000e+01 3.661756e-04 6.435808e-05 6.435808e-05 +1.475500e+04 4.000000e+01 4.000000e+01 3.335320e-04 6.132461e-05 6.132461e-05 +1.483500e+04 4.000000e+01 4.000000e+01 2.895342e-04 5.679261e-05 5.679261e-05 +1.491500e+04 4.000000e+01 4.000000e+01 3.321127e-04 6.091433e-05 6.091433e-05 +1.499500e+04 4.000000e+01 4.000000e+01 2.994692e-04 5.816034e-05 5.816034e-05 +1.507500e+04 4.000000e+01 4.000000e+01 3.108235e-04 5.887701e-05 5.887701e-05 +1.515500e+04 4.000000e+01 4.000000e+01 3.065656e-04 5.878974e-05 5.878974e-05 +1.523500e+04 4.000000e+01 4.000000e+01 3.122427e-04 5.934299e-05 5.934299e-05 +1.531500e+04 4.000000e+01 4.000000e+01 3.746913e-04 6.512309e-05 6.512309e-05 +1.539500e+04 4.000000e+01 4.000000e+01 3.675949e-04 6.439024e-05 6.439024e-05 +1.547500e+04 4.000000e+01 4.000000e+01 3.264356e-04 6.064579e-05 6.064579e-05 +1.555500e+04 4.000000e+01 4.000000e+01 2.909535e-04 5.714572e-05 5.714572e-05 +1.563500e+04 4.000000e+01 4.000000e+01 3.235970e-04 6.044865e-05 6.044865e-05 +1.571500e+04 4.000000e+01 4.000000e+01 3.094042e-04 5.900923e-05 5.900923e-05 +1.579500e+04 4.000000e+01 4.000000e+01 3.150813e-04 5.961247e-05 5.961247e-05 +1.587500e+04 4.000000e+01 4.000000e+01 3.065656e-04 5.851284e-05 5.851284e-05 +1.595500e+04 4.000000e+01 4.000000e+01 3.023077e-04 5.805852e-05 5.805852e-05 +1.603500e+04 4.000000e+01 4.000000e+01 3.122427e-04 5.942215e-05 5.942215e-05 +1.611500e+04 4.000000e+01 4.000000e+01 2.540520e-04 5.337551e-05 5.337551e-05 +1.619500e+04 4.000000e+01 4.000000e+01 2.810185e-04 5.634906e-05 5.634906e-05 +1.627500e+04 4.000000e+01 4.000000e+01 3.165006e-04 5.972942e-05 5.972942e-05 +1.635500e+04 4.000000e+01 4.000000e+01 3.278549e-04 6.058417e-05 6.058417e-05 +1.643500e+04 4.000000e+01 4.000000e+01 2.583099e-04 5.392980e-05 5.392980e-05 +1.651500e+04 4.000000e+01 4.000000e+01 3.023077e-04 5.819505e-05 5.819505e-05 +1.659500e+04 4.000000e+01 4.000000e+01 3.051463e-04 5.876303e-05 5.876303e-05 +1.667500e+04 4.000000e+01 4.000000e+01 2.966306e-04 5.789410e-05 5.789410e-05 +1.675500e+04 4.000000e+01 4.000000e+01 2.966306e-04 5.761280e-05 5.761280e-05 +1.683500e+04 4.000000e+01 4.000000e+01 2.923727e-04 5.726081e-05 5.726081e-05 +1.691500e+04 4.000000e+01 4.000000e+01 3.008885e-04 5.830696e-05 5.830696e-05 +1.699500e+04 4.000000e+01 4.000000e+01 2.937920e-04 5.767984e-05 5.767984e-05 +1.707500e+04 4.000000e+01 4.000000e+01 2.895342e-04 5.709018e-05 5.709018e-05 +1.715500e+04 4.000000e+01 4.000000e+01 3.023077e-04 5.767562e-05 5.767562e-05 +1.723500e+04 4.000000e+01 4.000000e+01 3.136620e-04 5.952522e-05 5.952522e-05 +1.731500e+04 4.000000e+01 4.000000e+01 2.994692e-04 5.815755e-05 5.815755e-05 +1.739500e+04 4.000000e+01 4.000000e+01 2.682449e-04 5.490307e-05 5.490307e-05 +1.747500e+04 4.000000e+01 4.000000e+01 2.654063e-04 5.469327e-05 5.469327e-05 +1.755500e+04 4.000000e+01 4.000000e+01 3.094042e-04 5.898838e-05 5.898838e-05 +1.763500e+04 4.000000e+01 4.000000e+01 2.583099e-04 5.391735e-05 5.391735e-05 +1.771500e+04 4.000000e+01 4.000000e+01 2.639870e-04 5.465870e-05 5.465870e-05 +1.779500e+04 4.000000e+01 4.000000e+01 2.725028e-04 5.536254e-05 5.536254e-05 +1.787500e+04 4.000000e+01 4.000000e+01 2.682449e-04 5.498942e-05 5.498942e-05 +1.795500e+04 4.000000e+01 4.000000e+01 2.540520e-04 5.343271e-05 5.343271e-05 +1.803500e+04 4.000000e+01 4.000000e+01 2.852763e-04 5.637301e-05 5.637301e-05 +1.811500e+04 4.000000e+01 4.000000e+01 2.356013e-04 5.137383e-05 5.137383e-05 +1.819500e+04 4.000000e+01 4.000000e+01 2.710835e-04 5.535635e-05 5.535635e-05 +1.827500e+04 4.000000e+01 4.000000e+01 2.952113e-04 5.767442e-05 5.767442e-05 +1.835500e+04 4.000000e+01 4.000000e+01 2.725028e-04 5.547180e-05 5.547180e-05 +1.843500e+04 4.000000e+01 4.000000e+01 2.725028e-04 5.517967e-05 5.517967e-05 +1.851500e+04 4.000000e+01 4.000000e+01 3.023077e-04 5.838547e-05 5.838547e-05 +1.859500e+04 4.000000e+01 4.000000e+01 2.512135e-04 5.314622e-05 5.314622e-05 +1.867500e+04 4.000000e+01 4.000000e+01 2.881149e-04 5.695395e-05 5.695395e-05 +1.875500e+04 4.000000e+01 4.000000e+01 2.327628e-04 5.112999e-05 5.112999e-05 +1.883500e+04 4.000000e+01 4.000000e+01 2.469556e-04 5.277762e-05 5.277762e-05 +1.891500e+04 4.000000e+01 4.000000e+01 2.568906e-04 5.363349e-05 5.363349e-05 +1.899500e+04 4.000000e+01 4.000000e+01 2.469556e-04 5.274202e-05 5.274202e-05 +1.907500e+04 4.000000e+01 4.000000e+01 2.980499e-04 5.782645e-05 5.782645e-05 +1.915500e+04 4.000000e+01 4.000000e+01 2.824377e-04 5.614573e-05 5.614573e-05 +1.923500e+04 4.000000e+01 4.000000e+01 2.327628e-04 5.120079e-05 5.120079e-05 +1.931500e+04 4.000000e+01 4.000000e+01 2.654063e-04 5.437152e-05 5.437152e-05 +1.939500e+04 4.000000e+01 4.000000e+01 2.497942e-04 5.315352e-05 5.315352e-05 +1.947500e+04 4.000000e+01 4.000000e+01 2.881149e-04 5.704033e-05 5.704033e-05 +1.955500e+04 4.000000e+01 4.000000e+01 2.795992e-04 5.579061e-05 5.579061e-05 +1.963500e+04 4.000000e+01 4.000000e+01 2.497942e-04 5.308180e-05 5.308180e-05 +1.971500e+04 4.000000e+01 4.000000e+01 2.356013e-04 5.148588e-05 5.148588e-05 +1.979500e+04 4.000000e+01 4.000000e+01 2.341821e-04 5.149227e-05 5.149227e-05 +1.987500e+04 4.000000e+01 4.000000e+01 2.540520e-04 5.352094e-05 5.352094e-05 +1.995500e+04 4.000000e+01 4.000000e+01 2.341821e-04 5.118034e-05 5.118034e-05 +2.003500e+04 4.000000e+01 4.000000e+01 2.639870e-04 5.445681e-05 5.445681e-05 +2.011500e+04 4.000000e+01 4.000000e+01 2.795992e-04 5.606828e-05 5.606828e-05 +2.019500e+04 4.000000e+01 4.000000e+01 2.654063e-04 5.484531e-05 5.484531e-05 +2.027500e+04 4.000000e+01 4.000000e+01 2.597292e-04 5.410610e-05 5.410610e-05 +2.035500e+04 4.000000e+01 4.000000e+01 2.256663e-04 5.031816e-05 5.031816e-05 +2.043500e+04 4.000000e+01 4.000000e+01 2.398592e-04 5.210676e-05 5.210676e-05 +2.051500e+04 4.000000e+01 4.000000e+01 2.441170e-04 5.245287e-05 5.245287e-05 +2.059500e+04 4.000000e+01 4.000000e+01 2.256663e-04 5.038660e-05 5.038660e-05 +2.067500e+04 4.000000e+01 4.000000e+01 2.568906e-04 5.391805e-05 5.391805e-05 +2.075500e+04 4.000000e+01 4.000000e+01 2.469556e-04 5.260272e-05 5.260272e-05 +2.083500e+04 4.000000e+01 4.000000e+01 2.384399e-04 5.185052e-05 5.185052e-05 +2.091500e+04 4.000000e+01 4.000000e+01 2.753413e-04 5.571849e-05 5.571849e-05 +2.099500e+04 4.000000e+01 4.000000e+01 2.583099e-04 5.401698e-05 5.401698e-05 +2.107500e+04 4.000000e+01 4.000000e+01 2.639870e-04 5.434508e-05 5.434508e-05 +2.115500e+04 4.000000e+01 4.000000e+01 2.199892e-04 4.967608e-05 4.967608e-05 +2.123500e+04 4.000000e+01 4.000000e+01 2.185699e-04 4.925870e-05 4.925870e-05 +2.131500e+04 4.000000e+01 4.000000e+01 2.441170e-04 5.254370e-05 5.254370e-05 +2.139500e+04 4.000000e+01 4.000000e+01 2.398592e-04 5.187614e-05 5.187614e-05 +2.147500e+04 4.000000e+01 4.000000e+01 2.228278e-04 4.986479e-05 4.986479e-05 +2.155500e+04 4.000000e+01 4.000000e+01 2.668256e-04 5.420680e-05 5.420680e-05 +2.163500e+04 4.000000e+01 4.000000e+01 2.412785e-04 5.213500e-05 5.213500e-05 +2.171500e+04 4.000000e+01 4.000000e+01 2.242471e-04 5.011550e-05 5.011550e-05 +2.179500e+04 4.000000e+01 4.000000e+01 2.185699e-04 4.965863e-05 4.965863e-05 +2.187500e+04 4.000000e+01 4.000000e+01 2.455363e-04 5.221552e-05 5.221552e-05 +2.195500e+04 4.000000e+01 4.000000e+01 2.143121e-04 4.891220e-05 4.891220e-05 +2.203500e+04 4.000000e+01 4.000000e+01 2.327628e-04 5.080642e-05 5.080642e-05 +2.211500e+04 4.000000e+01 4.000000e+01 2.313435e-04 5.083872e-05 5.083872e-05 +2.219500e+04 4.000000e+01 4.000000e+01 2.128928e-04 4.903816e-05 4.903816e-05 +2.227500e+04 4.000000e+01 4.000000e+01 2.412785e-04 5.220440e-05 5.220440e-05 +2.235500e+04 4.000000e+01 4.000000e+01 2.270856e-04 5.062004e-05 5.062004e-05 +2.243500e+04 4.000000e+01 4.000000e+01 2.285049e-04 5.060711e-05 5.060711e-05 +2.251500e+04 4.000000e+01 4.000000e+01 2.157313e-04 4.935963e-05 4.935963e-05 +2.259500e+04 4.000000e+01 4.000000e+01 2.100542e-04 4.835760e-05 4.835760e-05 +2.267500e+04 4.000000e+01 4.000000e+01 2.057963e-04 4.827165e-05 4.827165e-05 +2.275500e+04 4.000000e+01 4.000000e+01 2.114735e-04 4.874811e-05 4.874811e-05 +2.283500e+04 4.000000e+01 4.000000e+01 2.114735e-04 4.880105e-05 4.880105e-05 +2.291500e+04 4.000000e+01 4.000000e+01 2.313435e-04 5.116281e-05 5.116281e-05 +2.299500e+04 4.000000e+01 4.000000e+01 2.214085e-04 4.982117e-05 4.982117e-05 +2.307500e+04 4.000000e+01 4.000000e+01 2.469556e-04 5.261702e-05 5.261702e-05 +2.315500e+04 4.000000e+01 4.000000e+01 2.171506e-04 4.941320e-05 4.941320e-05 +2.323500e+04 4.000000e+01 4.000000e+01 2.057963e-04 4.814668e-05 4.814668e-05 +2.331500e+04 4.000000e+01 4.000000e+01 2.100542e-04 4.854065e-05 4.854065e-05 +2.339500e+04 4.000000e+01 4.000000e+01 2.143121e-04 4.916825e-05 4.916825e-05 +2.347500e+04 4.000000e+01 4.000000e+01 1.930228e-04 4.644337e-05 4.644337e-05 +2.355500e+04 4.000000e+01 4.000000e+01 2.015385e-04 4.755326e-05 4.755326e-05 +2.363500e+04 4.000000e+01 4.000000e+01 2.072156e-04 4.838871e-05 4.838871e-05 +2.371500e+04 4.000000e+01 4.000000e+01 1.958614e-04 4.661463e-05 4.661463e-05 +2.379500e+04 4.000000e+01 4.000000e+01 2.299242e-04 5.078741e-05 5.078741e-05 +2.387500e+04 4.000000e+01 4.000000e+01 1.958614e-04 4.667592e-05 4.667592e-05 +2.395500e+04 4.000000e+01 4.000000e+01 2.057963e-04 4.816618e-05 4.816618e-05 +2.403500e+04 4.000000e+01 4.000000e+01 2.015385e-04 4.756282e-05 4.756282e-05 +2.411500e+04 4.000000e+01 4.000000e+01 1.774106e-04 4.464010e-05 4.464010e-05 +2.419500e+04 4.000000e+01 4.000000e+01 2.256663e-04 5.032353e-05 5.032353e-05 +2.427500e+04 4.000000e+01 4.000000e+01 1.930228e-04 4.629368e-05 4.629368e-05 +2.435500e+04 4.000000e+01 4.000000e+01 1.575407e-04 4.176292e-05 4.176292e-05 +2.443500e+04 4.000000e+01 4.000000e+01 2.015385e-04 4.750595e-05 4.750595e-05 +2.451500e+04 4.000000e+01 4.000000e+01 2.057963e-04 4.828137e-05 4.828137e-05 +2.459500e+04 4.000000e+01 4.000000e+01 2.128928e-04 4.902944e-05 4.902944e-05 +2.467500e+04 4.000000e+01 4.000000e+01 2.341821e-04 5.135554e-05 5.135554e-05 +2.475500e+04 4.000000e+01 4.000000e+01 1.916035e-04 4.584117e-05 4.584117e-05 +2.483500e+04 4.000000e+01 4.000000e+01 1.958614e-04 4.679823e-05 4.679823e-05 +2.491500e+04 4.000000e+01 4.000000e+01 2.128928e-04 4.903658e-05 4.903658e-05 +2.499500e+04 4.000000e+01 4.000000e+01 2.072156e-04 4.812362e-05 4.812362e-05 +2.507500e+04 4.000000e+01 4.000000e+01 1.887649e-04 4.613457e-05 4.613457e-05 +2.515500e+04 4.000000e+01 4.000000e+01 2.057963e-04 4.808616e-05 4.808616e-05 +2.523500e+04 4.000000e+01 4.000000e+01 1.873456e-04 4.567130e-05 4.567130e-05 +2.531500e+04 4.000000e+01 4.000000e+01 2.185699e-04 4.965761e-05 4.965761e-05 +2.539500e+04 4.000000e+01 4.000000e+01 1.688949e-04 4.355975e-05 4.355975e-05 +2.547500e+04 4.000000e+01 4.000000e+01 1.802492e-04 4.501540e-05 4.501540e-05 +2.555500e+04 4.000000e+01 4.000000e+01 2.384399e-04 5.192490e-05 5.192490e-05 +2.563500e+04 4.000000e+01 4.000000e+01 1.972806e-04 4.662594e-05 4.662594e-05 +2.571500e+04 4.000000e+01 4.000000e+01 1.802492e-04 4.484941e-05 4.484941e-05 +2.579500e+04 4.000000e+01 4.000000e+01 1.972806e-04 4.698231e-05 4.698231e-05 +2.587500e+04 4.000000e+01 4.000000e+01 1.816685e-04 4.528034e-05 4.528034e-05 +2.595500e+04 4.000000e+01 4.000000e+01 2.256663e-04 5.047599e-05 5.047599e-05 +2.603500e+04 4.000000e+01 4.000000e+01 2.185699e-04 4.948139e-05 4.948139e-05 +2.611500e+04 4.000000e+01 4.000000e+01 1.774106e-04 4.475650e-05 4.475650e-05 +2.619500e+04 4.000000e+01 4.000000e+01 1.603792e-04 4.232364e-05 4.232364e-05 +2.627500e+04 4.000000e+01 4.000000e+01 1.859264e-04 4.546758e-05 4.546758e-05 +2.635500e+04 4.000000e+01 4.000000e+01 1.887649e-04 4.617734e-05 4.617734e-05 +2.643500e+04 4.000000e+01 4.000000e+01 2.001192e-04 4.726598e-05 4.726598e-05 +2.651500e+04 4.000000e+01 4.000000e+01 1.944421e-04 4.656725e-05 4.656725e-05 +2.659500e+04 4.000000e+01 4.000000e+01 2.015385e-04 4.777960e-05 4.777960e-05 +2.667500e+04 4.000000e+01 4.000000e+01 2.001192e-04 4.743971e-05 4.743971e-05 +2.675500e+04 4.000000e+01 4.000000e+01 2.043771e-04 4.789295e-05 4.789295e-05 +2.683500e+04 4.000000e+01 4.000000e+01 1.759914e-04 4.410870e-05 4.410870e-05 +2.691500e+04 4.000000e+01 4.000000e+01 2.199892e-04 4.977607e-05 4.977607e-05 +2.699500e+04 4.000000e+01 4.000000e+01 1.575407e-04 4.175121e-05 4.175121e-05 +2.707500e+04 4.000000e+01 4.000000e+01 1.788299e-04 4.468757e-05 4.468757e-05 +2.715500e+04 4.000000e+01 4.000000e+01 1.930228e-04 4.680477e-05 4.680477e-05 +2.723500e+04 4.000000e+01 4.000000e+01 1.802492e-04 4.500445e-05 4.500445e-05 +2.731500e+04 4.000000e+01 4.000000e+01 1.688949e-04 4.360841e-05 4.360841e-05 +2.739500e+04 4.000000e+01 4.000000e+01 1.603792e-04 4.222390e-05 4.222390e-05 +2.747500e+04 4.000000e+01 4.000000e+01 1.632178e-04 4.275044e-05 4.275044e-05 +2.755500e+04 4.000000e+01 4.000000e+01 1.944421e-04 4.674978e-05 4.674978e-05 +2.763500e+04 4.000000e+01 4.000000e+01 1.731528e-04 4.375995e-05 4.375995e-05 +2.771500e+04 4.000000e+01 4.000000e+01 1.731528e-04 4.366122e-05 4.366122e-05 +2.779500e+04 4.000000e+01 4.000000e+01 1.745721e-04 4.431599e-05 4.431599e-05 +2.787500e+04 4.000000e+01 4.000000e+01 1.830878e-04 4.509609e-05 4.509609e-05 +2.795500e+04 4.000000e+01 4.000000e+01 1.788299e-04 4.451405e-05 4.451405e-05 +2.803500e+04 4.000000e+01 4.000000e+01 1.816685e-04 4.526342e-05 4.526342e-05 +2.811500e+04 4.000000e+01 4.000000e+01 1.717335e-04 4.382221e-05 4.382221e-05 +2.819500e+04 4.000000e+01 4.000000e+01 1.617985e-04 4.237143e-05 4.237143e-05 +2.827500e+04 4.000000e+01 4.000000e+01 1.788299e-04 4.469869e-05 4.469869e-05 +2.835500e+04 4.000000e+01 4.000000e+01 1.703142e-04 4.355423e-05 4.355423e-05 +2.843500e+04 4.000000e+01 4.000000e+01 1.717335e-04 4.404040e-05 4.404040e-05 +2.851500e+04 4.000000e+01 4.000000e+01 1.532828e-04 4.145775e-05 4.145775e-05 +2.859500e+04 4.000000e+01 4.000000e+01 1.646371e-04 4.298562e-05 4.298562e-05 +2.867500e+04 4.000000e+01 4.000000e+01 1.703142e-04 4.367195e-05 4.367195e-05 +2.875500e+04 4.000000e+01 4.000000e+01 1.376707e-04 3.941178e-05 3.941178e-05 +2.883500e+04 4.000000e+01 4.000000e+01 1.646371e-04 4.286606e-05 4.286606e-05 +2.891500e+04 4.000000e+01 4.000000e+01 1.646371e-04 4.292074e-05 4.292074e-05 +2.899500e+04 4.000000e+01 4.000000e+01 1.518635e-04 4.113820e-05 4.113820e-05 +2.907500e+04 4.000000e+01 4.000000e+01 1.788299e-04 4.492151e-05 4.492151e-05 +2.915500e+04 4.000000e+01 4.000000e+01 1.504442e-04 4.096957e-05 4.096957e-05 +2.923500e+04 4.000000e+01 4.000000e+01 1.603792e-04 4.232051e-05 4.232051e-05 +2.931500e+04 4.000000e+01 4.000000e+01 1.589599e-04 4.204809e-05 4.204809e-05 +2.939500e+04 4.000000e+01 4.000000e+01 1.518635e-04 4.111326e-05 4.111326e-05 +2.947500e+04 4.000000e+01 4.000000e+01 1.447671e-04 4.037690e-05 4.037690e-05 +2.955500e+04 4.000000e+01 4.000000e+01 1.745721e-04 4.391775e-05 4.391775e-05 +2.963500e+04 4.000000e+01 4.000000e+01 1.561214e-04 4.152373e-05 4.152373e-05 +2.971500e+04 4.000000e+01 4.000000e+01 1.617985e-04 4.250165e-05 4.250165e-05 +2.979500e+04 4.000000e+01 4.000000e+01 1.901842e-04 4.627518e-05 4.627518e-05 +2.987500e+04 4.000000e+01 4.000000e+01 1.660564e-04 4.315896e-05 4.315896e-05 +2.995500e+04 4.000000e+01 4.000000e+01 1.788299e-04 4.471751e-05 4.471751e-05 +3.003500e+04 4.000000e+01 4.000000e+01 1.419285e-04 3.970952e-05 3.970952e-05 +3.011500e+04 4.000000e+01 4.000000e+01 1.348321e-04 3.890900e-05 3.890900e-05 +3.019500e+04 4.000000e+01 4.000000e+01 1.873456e-04 4.595201e-05 4.595201e-05 +3.027500e+04 4.000000e+01 4.000000e+01 1.376707e-04 3.889911e-05 3.889911e-05 +3.035500e+04 4.000000e+01 4.000000e+01 1.532828e-04 4.138807e-05 4.138807e-05 +3.043500e+04 4.000000e+01 4.000000e+01 1.759914e-04 4.421429e-05 4.421429e-05 +3.051500e+04 4.000000e+01 4.000000e+01 1.447671e-04 4.018666e-05 4.018666e-05 +3.059500e+04 4.000000e+01 4.000000e+01 1.476057e-04 4.073679e-05 4.073679e-05 +3.067500e+04 4.000000e+01 4.000000e+01 1.405092e-04 3.964898e-05 3.964898e-05 +3.075500e+04 4.000000e+01 4.000000e+01 1.646371e-04 4.302566e-05 4.302566e-05 +3.083500e+04 4.000000e+01 4.000000e+01 1.589599e-04 4.207328e-05 4.207328e-05 +3.091500e+04 4.000000e+01 4.000000e+01 1.348321e-04 3.832361e-05 3.832361e-05 +3.099500e+04 4.000000e+01 4.000000e+01 1.788299e-04 4.471903e-05 4.471903e-05 +3.107500e+04 4.000000e+01 4.000000e+01 1.376707e-04 3.931292e-05 3.931292e-05 +3.115500e+04 4.000000e+01 4.000000e+01 1.263164e-04 3.755523e-05 3.755523e-05 +3.123500e+04 4.000000e+01 4.000000e+01 1.476057e-04 4.043961e-05 4.043961e-05 +3.131500e+04 4.000000e+01 4.000000e+01 1.390899e-04 3.892625e-05 3.892625e-05 +3.139500e+04 4.000000e+01 4.000000e+01 1.291550e-04 3.800613e-05 3.800613e-05 +3.147500e+04 4.000000e+01 4.000000e+01 1.277357e-04 3.785455e-05 3.785455e-05 +3.155500e+04 4.000000e+01 4.000000e+01 1.490249e-04 4.058002e-05 4.058002e-05 +3.163500e+04 4.000000e+01 4.000000e+01 1.348321e-04 3.862771e-05 3.862771e-05 +3.171500e+04 4.000000e+01 4.000000e+01 1.263164e-04 3.761028e-05 3.761028e-05 +3.179500e+04 4.000000e+01 4.000000e+01 1.348321e-04 3.895302e-05 3.895302e-05 +3.187500e+04 4.000000e+01 4.000000e+01 1.504442e-04 4.120451e-05 4.120451e-05 +3.195500e+04 4.000000e+01 4.000000e+01 1.220585e-04 3.692279e-05 3.692279e-05 +3.203500e+04 4.000000e+01 4.000000e+01 1.433478e-04 4.019322e-05 4.019322e-05 +3.211500e+04 4.000000e+01 4.000000e+01 1.532828e-04 4.124929e-05 4.124929e-05 +3.219500e+04 4.000000e+01 4.000000e+01 1.121235e-04 3.529266e-05 3.529266e-05 +3.227500e+04 4.000000e+01 4.000000e+01 1.291550e-04 3.791230e-05 3.791230e-05 +3.235500e+04 4.000000e+01 4.000000e+01 1.220585e-04 3.674997e-05 3.674997e-05 +3.243500e+04 4.000000e+01 4.000000e+01 1.192200e-04 3.640876e-05 3.640876e-05 +3.251500e+04 4.000000e+01 4.000000e+01 1.036078e-04 3.415217e-05 3.415217e-05 +3.259500e+04 4.000000e+01 4.000000e+01 1.248971e-04 3.730474e-05 3.730474e-05 +3.267500e+04 4.000000e+01 4.000000e+01 1.433478e-04 3.985952e-05 3.985952e-05 +3.275500e+04 4.000000e+01 4.000000e+01 1.192200e-04 3.644155e-05 3.644155e-05 +3.283500e+04 4.000000e+01 4.000000e+01 1.348321e-04 3.877586e-05 3.877586e-05 +3.291500e+04 4.000000e+01 4.000000e+01 9.509211e-05 3.217821e-05 3.217821e-05 +3.299500e+04 4.000000e+01 4.000000e+01 1.007692e-04 3.324930e-05 3.324930e-05 +3.307500e+04 4.000000e+01 4.000000e+01 8.373782e-05 3.029251e-05 3.029251e-05 +3.315500e+04 4.000000e+01 4.000000e+01 1.206392e-04 3.680709e-05 3.680709e-05 +3.323500e+04 4.000000e+01 4.000000e+01 9.367282e-05 3.214416e-05 3.214416e-05 +3.331500e+04 4.000000e+01 4.000000e+01 1.036078e-04 3.390299e-05 3.390299e-05 +3.339500e+04 4.000000e+01 4.000000e+01 8.799568e-05 3.140884e-05 3.140884e-05 +3.347500e+04 4.000000e+01 4.000000e+01 8.515711e-05 3.068063e-05 3.068063e-05 +3.355500e+04 4.000000e+01 4.000000e+01 8.941497e-05 3.144366e-05 3.144366e-05 +3.363500e+04 4.000000e+01 4.000000e+01 6.954497e-05 2.743959e-05 2.743959e-05 +3.371500e+04 4.000000e+01 4.000000e+01 8.231854e-05 3.028948e-05 3.028948e-05 +3.379500e+04 4.000000e+01 4.000000e+01 7.380283e-05 2.849747e-05 2.849747e-05 +3.387500e+04 4.000000e+01 4.000000e+01 7.096426e-05 2.757165e-05 2.757165e-05 +3.395500e+04 4.000000e+01 4.000000e+01 6.954497e-05 2.684988e-05 2.684988e-05 +3.403500e+04 4.000000e+01 4.000000e+01 7.096426e-05 2.800789e-05 2.800789e-05 +3.411500e+04 4.000000e+01 4.000000e+01 5.819069e-05 2.489078e-05 2.489078e-05 +3.419500e+04 4.000000e+01 4.000000e+01 3.690141e-05 1.998709e-05 1.998709e-05 +3.427500e+04 4.000000e+01 4.000000e+01 6.954497e-05 2.782947e-05 2.782947e-05 +3.435500e+04 4.000000e+01 4.000000e+01 4.967498e-05 2.305776e-05 2.305776e-05 +3.443500e+04 4.000000e+01 4.000000e+01 4.541713e-05 2.184607e-05 2.184607e-05 +3.451500e+04 4.000000e+01 4.000000e+01 3.690141e-05 1.963130e-05 1.963130e-05 +3.459500e+04 4.000000e+01 4.000000e+01 3.690141e-05 2.035598e-05 2.035598e-05 +3.467500e+04 4.000000e+01 4.000000e+01 2.980499e-05 1.699551e-05 1.699551e-05 +3.475500e+04 4.000000e+01 4.000000e+01 3.264356e-05 1.856779e-05 1.856779e-05 +3.483500e+04 4.000000e+01 4.000000e+01 2.412785e-05 1.494309e-05 1.494309e-05 +3.491500e+04 4.000000e+01 4.000000e+01 2.980499e-05 1.803450e-05 1.803450e-05 +3.499500e+04 4.000000e+01 4.000000e+01 1.703142e-05 1.273765e-05 1.273765e-05 +3.507500e+04 4.000000e+01 4.000000e+01 2.554713e-05 1.636237e-05 1.636237e-05 +3.515500e+04 4.000000e+01 4.000000e+01 2.270856e-05 1.303123e-05 1.303123e-05 +3.523500e+04 4.000000e+01 4.000000e+01 1.703142e-05 1.073047e-05 1.073047e-05 +3.531500e+04 4.000000e+01 4.000000e+01 1.277357e-05 9.174407e-06 9.174407e-06 +3.539500e+04 4.000000e+01 4.000000e+01 1.561214e-05 1.014259e-05 1.014259e-05 +3.547500e+04 4.000000e+01 4.000000e+01 1.703142e-05 1.176946e-05 1.176946e-05 +3.555500e+04 4.000000e+01 4.000000e+01 8.515711e-06 6.852916e-06 6.852916e-06 +3.563500e+04 4.000000e+01 4.000000e+01 1.135428e-05 8.723301e-06 8.723301e-06 +3.571500e+04 4.000000e+01 4.000000e+01 9.934996e-06 7.304018e-06 7.304018e-06 +3.579500e+04 4.000000e+01 4.000000e+01 4.257856e-06 4.257856e-06 4.257856e-06 +3.587500e+04 4.000000e+01 4.000000e+01 5.677141e-06 5.677141e-06 5.677141e-06 +3.595500e+04 4.000000e+01 4.000000e+01 8.515711e-06 6.852916e-06 6.852916e-06 +3.603500e+04 4.000000e+01 4.000000e+01 4.257856e-06 3.426459e-06 3.426459e-06 +3.611500e+04 4.000000e+01 4.000000e+01 1.419285e-06 1.419285e-06 1.419285e-06 +3.619500e+04 4.000000e+01 4.000000e+01 1.419285e-06 1.419285e-06 1.419285e-06 +3.627500e+04 4.000000e+01 4.000000e+01 0.000000e+00 0.000000e+00 0.000000e+00 +3.635500e+04 4.000000e+01 4.000000e+01 0.000000e+00 0.000000e+00 0.000000e+00 +3.643500e+04 4.000000e+01 4.000000e+01 0.000000e+00 0.000000e+00 0.000000e+00 +3.651500e+04 4.000000e+01 4.000000e+01 1.419285e-06 1.419285e-06 1.419285e-06 +3.659500e+04 4.000000e+01 4.000000e+01 0.000000e+00 0.000000e+00 0.000000e+00 +3.667500e+04 4.000000e+01 4.000000e+01 2.838570e-06 2.838570e-06 2.838570e-06 +3.675500e+04 4.000000e+01 4.000000e+01 1.419285e-06 1.419285e-06 1.419285e-06 +3.683500e+04 4.000000e+01 4.000000e+01 0.000000e+00 0.000000e+00 0.000000e+00 +3.691500e+04 4.000000e+01 4.000000e+01 1.419285e-06 1.419285e-06 1.419285e-06 +3.699500e+04 4.000000e+01 4.000000e+01 0.000000e+00 0.000000e+00 0.000000e+00 +3.707500e+04 4.000000e+01 4.000000e+01 1.419285e-06 1.419285e-06 1.419285e-06 +3.715500e+04 4.000000e+01 4.000000e+01 0.000000e+00 0.000000e+00 0.000000e+00 +3.723500e+04 4.000000e+01 4.000000e+01 0.000000e+00 0.000000e+00 0.000000e+00 +3.731500e+04 4.000000e+01 4.000000e+01 0.000000e+00 0.000000e+00 0.000000e+00 +3.739500e+04 4.000000e+01 4.000000e+01 0.000000e+00 0.000000e+00 0.000000e+00 +3.747500e+04 4.000000e+01 4.000000e+01 1.419285e-06 1.419285e-06 1.419285e-06 +3.755500e+04 4.000000e+01 4.000000e+01 0.000000e+00 0.000000e+00 0.000000e+00 +3.763500e+04 4.000000e+01 4.000000e+01 0.000000e+00 0.000000e+00 0.000000e+00 +3.771500e+04 4.000000e+01 4.000000e+01 0.000000e+00 0.000000e+00 0.000000e+00 +3.779500e+04 4.000000e+01 4.000000e+01 0.000000e+00 0.000000e+00 0.000000e+00 +3.787500e+04 4.000000e+01 4.000000e+01 0.000000e+00 0.000000e+00 0.000000e+00 +3.795500e+04 4.000000e+01 4.000000e+01 0.000000e+00 0.000000e+00 0.000000e+00 +3.803500e+04 4.000000e+01 4.000000e+01 0.000000e+00 0.000000e+00 0.000000e+00 +3.811500e+04 4.000000e+01 4.000000e+01 0.000000e+00 0.000000e+00 0.000000e+00 +3.819500e+04 4.000000e+01 4.000000e+01 0.000000e+00 0.000000e+00 0.000000e+00 +3.827500e+04 4.000000e+01 4.000000e+01 0.000000e+00 0.000000e+00 0.000000e+00 +3.835500e+04 4.000000e+01 4.000000e+01 0.000000e+00 0.000000e+00 0.000000e+00 +3.843500e+04 4.000000e+01 4.000000e+01 0.000000e+00 0.000000e+00 0.000000e+00 +3.851500e+04 4.000000e+01 4.000000e+01 0.000000e+00 0.000000e+00 0.000000e+00 +3.859500e+04 4.000000e+01 4.000000e+01 0.000000e+00 0.000000e+00 0.000000e+00 +3.867500e+04 4.000000e+01 4.000000e+01 0.000000e+00 0.000000e+00 0.000000e+00 +3.875500e+04 4.000000e+01 4.000000e+01 0.000000e+00 0.000000e+00 0.000000e+00 +3.883500e+04 4.000000e+01 4.000000e+01 0.000000e+00 0.000000e+00 0.000000e+00 +3.891500e+04 4.000000e+01 4.000000e+01 0.000000e+00 0.000000e+00 0.000000e+00 +3.899500e+04 4.000000e+01 4.000000e+01 0.000000e+00 0.000000e+00 0.000000e+00 +3.907500e+04 4.000000e+01 4.000000e+01 0.000000e+00 0.000000e+00 0.000000e+00 +3.915500e+04 4.000000e+01 4.000000e+01 0.000000e+00 0.000000e+00 0.000000e+00 +3.923500e+04 4.000000e+01 4.000000e+01 0.000000e+00 0.000000e+00 0.000000e+00 +3.931500e+04 4.000000e+01 4.000000e+01 0.000000e+00 0.000000e+00 0.000000e+00 +3.939500e+04 4.000000e+01 4.000000e+01 0.000000e+00 0.000000e+00 0.000000e+00 +3.947500e+04 4.000000e+01 4.000000e+01 0.000000e+00 0.000000e+00 0.000000e+00 +3.955500e+04 4.000000e+01 4.000000e+01 0.000000e+00 0.000000e+00 0.000000e+00 +3.963500e+04 4.000000e+01 4.000000e+01 0.000000e+00 0.000000e+00 0.000000e+00 +3.971500e+04 4.000000e+01 4.000000e+01 0.000000e+00 0.000000e+00 0.000000e+00 +3.979500e+04 4.000000e+01 4.000000e+01 0.000000e+00 0.000000e+00 0.000000e+00 +3.987500e+04 4.000000e+01 4.000000e+01 0.000000e+00 0.000000e+00 0.000000e+00 +3.995500e+04 4.000000e+01 4.000000e+01 0.000000e+00 0.000000e+00 0.000000e+00 +END YODA_SCATTER2D_V2 diff --git a/analyses/pluginMC/MC_OPTIONS.cc b/analyses/pluginMC/MC_OPTIONS.cc new file mode 100644 --- /dev/null +++ b/analyses/pluginMC/MC_OPTIONS.cc @@ -0,0 +1,80 @@ +// -*- C++ -*- +#include "Rivet/Analysis.hh" + +namespace Rivet { + // Example analysis to show how to use options in an analysis + // Example of a custom class to be read in from an option. + class A { + public: + A() : a(-1.0) {} + private: + double a; + // Custom class must be streamable. + friend istream& operator>> (istream& is, A& a); + friend ostream& operator<< (ostream& os, const A& a); + }; + + // Custom class must be streamable. + istream& operator>> (istream& is, A& a) { + is >> a.a; + return is; + } + ostream& operator<< (ostream& os, const A& a) { + os << a.a; + return os; + } + + class MC_OPTIONS : public Analysis { + public: + + /// Constructor + DEFAULT_RIVET_ANALYSIS_CTOR(MC_OPTIONS); + + + /// @name Analysis methods + //@{ + + /// Book histograms and initialise projections before the run + void init() { + + // Parameters read in. + // A double. + double f = getOption("foo", 1.0); + // A string. + string s = getOption("bar", ""); + // A custom object. + A a = getOption("baz", A()); + + cout << "foo = " << f << endl; + cout << "bar = " << s << endl; + cout << "baz = " << a << endl; + } + + + /// Perform the per-event analysis + void analyze(const Event& event) { + + } + + + /// Finalize + void finalize() { + } + + //@} + + + private: + + /// @name Histograms + //@{ + //@} + + + }; + + + // The hook for the plugin system + DECLARE_RIVET_PLUGIN(MC_OPTIONS); + +} diff --git a/include/Rivet/Analysis.hh b/include/Rivet/Analysis.hh --- a/include/Rivet/Analysis.hh +++ b/include/Rivet/Analysis.hh @@ -1,1124 +1,1162 @@ // -*- C++ -*- #ifndef RIVET_Analysis_HH #define RIVET_Analysis_HH #include "Rivet/Config/RivetCommon.hh" #include "Rivet/AnalysisInfo.hh" #include "Rivet/Event.hh" #include "Rivet/Projection.hh" #include "Rivet/ProjectionApplier.hh" #include "Rivet/ProjectionHandler.hh" #include "Rivet/AnalysisLoader.hh" #include "Rivet/Tools/Cuts.hh" #include "Rivet/Tools/Logging.hh" #include "Rivet/Tools/ParticleUtils.hh" #include "Rivet/Tools/BinnedHistogram.hh" #include "Rivet/Tools/RivetMT2.hh" #include "Rivet/Tools/RivetYODA.hh" +#include "Rivet/Projections/CentralityProjection.hh" /// @def vetoEvent /// Preprocessor define for vetoing events, including the log message and return. #define vetoEvent \ do { MSG_DEBUG("Vetoing event on line " << __LINE__ << " of " << __FILE__); return; } while(0) namespace Rivet { // Forward declaration class AnalysisHandler; /// @brief This is the base class of all analysis classes in Rivet. /// /// There are /// three virtual functions which should be implemented in base classes: /// /// void init() is called by Rivet before a run is started. Here the /// analysis class should book necessary histograms. The needed /// projections should probably rather be constructed in the /// constructor. /// /// void analyze(const Event&) is called once for each event. Here the /// analysis class should apply the necessary Projections and fill the /// histograms. /// /// void finalize() is called after a run is finished. Here the analysis /// class should do whatever manipulations are necessary on the /// histograms. Writing the histograms to a file is, however, done by /// the Rivet class. class Analysis : public ProjectionApplier { /// The AnalysisHandler is a friend. friend class AnalysisHandler; public: /// @name Standard constructors and destructors. //@{ // /// The default constructor. // Analysis(); /// Constructor Analysis(const std::string& name); /// The destructor. virtual ~Analysis() {} //@} public: /// @name Main analysis methods //@{ /// Initialize this analysis object. A concrete class should here /// book all necessary histograms. An overridden function must make /// sure it first calls the base class function. virtual void init() { } /// Analyze one event. A concrete class should here apply the /// necessary projections on the \a event and fill the relevant /// histograms. An overridden function must make sure it first calls /// the base class function. virtual void analyze(const Event& event) = 0; /// Finalize this analysis object. A concrete class should here make /// all necessary operations on the histograms. Writing the /// histograms to a file is, however, done by the Rivet class. An /// overridden function must make sure it first calls the base class /// function. virtual void finalize() { } //@} public: /// @name Metadata /// Metadata is used for querying from the command line and also for /// building web pages and the analysis pages in the Rivet manual. //@{ /// Get the actual AnalysisInfo object in which all this metadata is stored. const AnalysisInfo& info() const { assert(_info && "No AnalysisInfo object :O"); return *_info; } /// @brief Get the name of the analysis. /// /// By default this is computed by combining the results of the /// experiment, year and Spires ID metadata methods and you should /// only override it if there's a good reason why those won't /// work. If options has been set for this instance, a /// corresponding string is appended at the end. virtual std::string name() const { return ( (info().name().empty()) ? _defaultname : info().name() ) + _optstring; } // get name of reference data file, which could be different from plugin name virtual std::string getRefDataName() const { return (info().getRefDataName().empty()) ? _defaultname : info().getRefDataName(); } // set name of reference data file, which could be different from plugin name virtual void setRefDataName(const std::string& ref_data="") { info().setRefDataName(!ref_data.empty() ? ref_data : name()); } /// Get the Inspire ID code for this analysis. virtual std::string inspireId() const { return info().inspireId(); } /// Get the SPIRES ID code for this analysis (~deprecated). virtual std::string spiresId() const { return info().spiresId(); } /// @brief Names & emails of paper/analysis authors. /// /// Names and email of authors in 'NAME \' format. The first /// name in the list should be the primary contact person. virtual std::vector authors() const { return info().authors(); } /// @brief Get a short description of the analysis. /// /// Short (one sentence) description used as an index entry. /// Use @a description() to provide full descriptive paragraphs /// of analysis details. virtual std::string summary() const { return info().summary(); } /// @brief Get a full description of the analysis. /// /// Full textual description of this analysis, what it is useful for, /// what experimental techniques are applied, etc. Should be treated /// as a chunk of restructuredText (http://docutils.sourceforge.net/rst.html), /// with equations to be rendered as LaTeX with amsmath operators. virtual std::string description() const { return info().description(); } /// @brief Information about the events needed as input for this analysis. /// /// Event types, energies, kinematic cuts, particles to be considered /// stable, etc. etc. Should be treated as a restructuredText bullet list /// (http://docutils.sourceforge.net/rst.html) virtual std::string runInfo() const { return info().runInfo(); } /// Experiment which performed and published this analysis. virtual std::string experiment() const { return info().experiment(); } /// Collider on which the experiment ran. virtual std::string collider() const { return info().collider(); } /// When the original experimental analysis was published. virtual std::string year() const { return info().year(); } /// The luminosity in inverse femtobarn virtual std::string luminosityfb() const { return info().luminosityfb(); } /// Journal, and preprint references. virtual std::vector references() const { return info().references(); } /// BibTeX citation key for this article. virtual std::string bibKey() const { return info().bibKey(); } /// BibTeX citation entry for this article. virtual std::string bibTeX() const { return info().bibTeX(); } /// Whether this analysis is trusted (in any way!) virtual std::string status() const { return (info().status().empty()) ? "UNVALIDATED" : info().status(); } /// Any work to be done on this analysis. virtual std::vector todos() const { return info().todos(); } /// Return the allowed pairs of incoming beams required by this analysis. virtual const std::vector& requiredBeams() const { return info().beams(); } /// Declare the allowed pairs of incoming beams required by this analysis. virtual Analysis& setRequiredBeams(const std::vector& requiredBeams) { info().setBeams(requiredBeams); return *this; } /// Sets of valid beam energy pairs, in GeV virtual const std::vector >& requiredEnergies() const { return info().energies(); } /// Get vector of analysis keywords virtual const std::vector & keywords() const { return info().keywords(); } /// Declare the list of valid beam energy pairs, in GeV virtual Analysis& setRequiredEnergies(const std::vector >& requiredEnergies) { info().setEnergies(requiredEnergies); return *this; } /// Return true if this analysis needs to know the process cross-section. /// @todo Remove this and require HepMC >= 2.06 bool needsCrossSection() const { return info().needsCrossSection(); } /// Declare whether this analysis needs to know the process cross-section from the generator. /// @todo Remove this and require HepMC >= 2.06 Analysis& setNeedsCrossSection(bool needed=true) { info().setNeedsCrossSection(needed); return *this; } //@} /// @name Internal metadata modifying methods //@{ /// Get the actual AnalysisInfo object in which all this metadata is stored (non-const). AnalysisInfo& info() { assert(_info && "No AnalysisInfo object :O"); return *_info; } //@} /// @name Run conditions //@{ /// Incoming beams for this run const ParticlePair& beams() const; /// Incoming beam IDs for this run const PdgIdPair beamIds() const; /// Centre of mass energy for this run double sqrtS() const; //@} /// @name Analysis / beam compatibility testing //@{ /// Check if analysis is compatible with the provided beam particle IDs and energies bool isCompatible(const ParticlePair& beams) const; /// Check if analysis is compatible with the provided beam particle IDs and energies bool isCompatible(PdgId beam1, PdgId beam2, double e1, double e2) const; /// Check if analysis is compatible with the provided beam particle IDs and energies bool isCompatible(const PdgIdPair& beams, const std::pair& energies) const; //@} /// Set the cross section from the generator Analysis& setCrossSection(double xs); //, double xserr=0.0); /// Access the controlling AnalysisHandler object. AnalysisHandler& handler() const { return *_analysishandler; } // get name of reference data file, which could be different from plugin name - virtual std::string getRefDataName() const { + /* virtual std::string getRefDataName() const { return (info().getRefDataName().empty()) ? _defaultname : info().getRefDataName(); } // set name of reference data file, which could be different from plugin name virtual void setRefDataName(const std::string& ref_data="") { info().setRefDataName(!ref_data.empty() ? ref_data : name()); } - +*/ protected: /// Get a Log object based on the name() property of the calling analysis object. Log& getLog() const; /// Get the process cross-section in pb. Throws if this hasn't been set. double crossSection() const; /// Get the process cross-section per generated event in pb. Throws if this /// hasn't been set. double crossSectionPerEvent() const; /// @brief Get the number of events seen (via the analysis handler). /// /// @note Use in the finalize phase only. size_t numEvents() const; /// @brief Get the sum of event weights seen (via the analysis handler). /// /// @note Use in the finalize phase only. double sumW() const; /// Alias double sumOfWeights() const { return sumW(); } /// @brief Get the sum of squared event weights seen (via the analysis handler). /// /// @note Use in the finalize phase only. double sumW2() const; protected: /// @name Histogram paths //@{ /// Get the canonical histogram "directory" path for this analysis. const std::string histoDir() const; /// Get the canonical histogram path for the named histogram in this analysis. const std::string histoPath(const std::string& hname) const; /// Get the canonical histogram path for the numbered histogram in this analysis. const std::string histoPath(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const; /// Get the internal histogram name for given d, x and y (cf. HepData) const std::string mkAxisCode(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const; /// Alias /// @deprecated Prefer the "mk" form, consistent with other "making function" names const std::string makeAxisCode(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const { return mkAxisCode(datasetId, xAxisId, yAxisId); } //@} /// @name Histogram reference data //@{ /// Get reference data for a named histo /// @todo SFINAE to ensure that the type inherits from YODA::AnalysisObject? template const T& refData(const string& hname) const { _cacheRefData(); MSG_TRACE("Using histo bin edges for " << name() << ":" << hname); if (!_refdata[hname]) { MSG_ERROR("Can't find reference histogram " << hname); throw Exception("Reference data " + hname + " not found."); } return dynamic_cast(*_refdata[hname]); } /// Get reference data for a numbered histo /// @todo SFINAE to ensure that the type inherits from YODA::AnalysisObject? template const T& refData(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const { const string hname = makeAxisCode(datasetId, xAxisId, yAxisId); return refData(hname); } //@} /// @name Counter booking //@{ /// Book a counter. CounterPtr bookCounter(const std::string& name, const std::string& title=""); // const std::string& valtitle="" /// Book a counter, using a path generated from the dataset and axis ID codes /// /// The paper, dataset and x/y-axis IDs will be used to build the histo name in the HepData standard way. CounterPtr bookCounter(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId, const std::string& title=""); // const std::string& valtitle="" //@} /// @name 1D histogram booking //@{ /// Book a 1D histogram with @a nbins uniformly distributed across the range @a lower - @a upper . Histo1DPtr bookHisto1D(const std::string& name, size_t nbins, double lower, double upper, const std::string& title="", const std::string& xtitle="", const std::string& ytitle=""); /// Book a 1D histogram with non-uniform bins defined by the vector of bin edges @a binedges . Histo1DPtr bookHisto1D(const std::string& name, const std::vector& binedges, const std::string& title="", const std::string& xtitle="", const std::string& ytitle=""); /// Book a 1D histogram with non-uniform bins defined by the vector of bin edges @a binedges . Histo1DPtr bookHisto1D(const std::string& name, const std::initializer_list& binedges, const std::string& title="", const std::string& xtitle="", const std::string& ytitle=""); /// Book a 1D histogram with binning from a reference scatter. Histo1DPtr bookHisto1D(const std::string& name, const Scatter2D& refscatter, const std::string& title="", const std::string& xtitle="", const std::string& ytitle=""); /// Book a 1D histogram, using the binnings in the reference data histogram. Histo1DPtr bookHisto1D(const std::string& name, const std::string& title="", const std::string& xtitle="", const std::string& ytitle=""); /// Book a 1D histogram, using the binnings in the reference data histogram. /// /// The paper, dataset and x/y-axis IDs will be used to build the histo name in the HepData standard way. Histo1DPtr bookHisto1D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId, const std::string& title="", const std::string& xtitle="", const std::string& ytitle=""); //@} /// @name 2D histogram booking //@{ /// Book a 2D histogram with @a nxbins and @a nybins uniformly /// distributed across the ranges @a xlower - @a xupper and @a /// ylower - @a yupper respectively along the x- and y-axis. Histo2DPtr bookHisto2D(const std::string& name, size_t nxbins, double xlower, double xupper, size_t nybins, double ylower, double yupper, const std::string& title="", const std::string& xtitle="", const std::string& ytitle="", const std::string& ztitle=""); /// Book a 2D histogram with non-uniform bins defined by the /// vectors of bin edges @a xbinedges and @a ybinedges. Histo2DPtr bookHisto2D(const std::string& name, const std::vector& xbinedges, const std::vector& ybinedges, const std::string& title="", const std::string& xtitle="", const std::string& ytitle="", const std::string& ztitle=""); /// Book a 2D histogram with non-uniform bins defined by the /// vectors of bin edges @a xbinedges and @a ybinedges. Histo2DPtr bookHisto2D(const std::string& name, const std::initializer_list& xbinedges, const std::initializer_list& ybinedges, const std::string& title="", const std::string& xtitle="", const std::string& ytitle="", const std::string& ztitle=""); /// Book a 2D histogram with binning from a reference scatter. Histo2DPtr bookHisto2D(const std::string& name, const Scatter3D& refscatter, const std::string& title="", const std::string& xtitle="", const std::string& ytitle="", const std::string& ztitle=""); /// Book a 2D histogram, using the binnings in the reference data histogram. Histo2DPtr bookHisto2D(const std::string& name, const std::string& title="", const std::string& xtitle="", const std::string& ytitle="", const std::string& ztitle=""); /// Book a 2D histogram, using the binnings in the reference data histogram. /// /// The paper, dataset and x/y-axis IDs will be used to build the histo name in the HepData standard way. Histo2DPtr bookHisto2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId, const std::string& title="", const std::string& xtitle="", const std::string& ytitle="", const std::string& ztitle=""); //@} /// @name 1D profile histogram booking //@{ /// Book a 1D profile histogram with @a nbins uniformly distributed across the range @a lower - @a upper . Profile1DPtr bookProfile1D(const std::string& name, size_t nbins, double lower, double upper, const std::string& title="", const std::string& xtitle="", const std::string& ytitle=""); /// Book a 1D profile histogram with non-uniform bins defined by the vector of bin edges @a binedges . Profile1DPtr bookProfile1D(const std::string& name, const std::vector& binedges, const std::string& title="", const std::string& xtitle="", const std::string& ytitle=""); /// Book a 1D profile histogram with non-uniform bins defined by the vector of bin edges @a binedges . Profile1DPtr bookProfile1D(const std::string& name, const std::initializer_list& binedges, const std::string& title="", const std::string& xtitle="", const std::string& ytitle=""); /// Book a 1D profile histogram with binning from a reference scatter. Profile1DPtr bookProfile1D(const std::string& name, const Scatter2D& refscatter, const std::string& title="", const std::string& xtitle="", const std::string& ytitle=""); /// Book a 1D profile histogram, using the binnings in the reference data histogram. Profile1DPtr bookProfile1D(const std::string& name, const std::string& title="", const std::string& xtitle="", const std::string& ytitle=""); /// Book a 1D profile histogram, using the binnings in the reference data histogram. /// /// The paper, dataset and x/y-axis IDs will be used to build the histo name in the HepData standard way. Profile1DPtr bookProfile1D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId, const std::string& title="", const std::string& xtitle="", const std::string& ytitle=""); //@} /// @name 2D profile histogram booking //@{ /// Book a 2D profile histogram with @a nxbins and @a nybins uniformly /// distributed across the ranges @a xlower - @a xupper and @a ylower - @a /// yupper respectively along the x- and y-axis. Profile2DPtr bookProfile2D(const std::string& name, size_t nxbins, double xlower, double xupper, size_t nybins, double ylower, double yupper, const std::string& title="", const std::string& xtitle="", const std::string& ytitle="", const std::string& ztitle=""); /// Book a 2D profile histogram with non-uniform bins defined by the vectorx /// of bin edges @a xbinedges and @a ybinedges. Profile2DPtr bookProfile2D(const std::string& name, const std::vector& xbinedges, const std::vector& ybinedges, const std::string& title="", const std::string& xtitle="", const std::string& ytitle="", const std::string& ztitle=""); /// Book a 2D profile histogram with non-uniform bins defined by the vectorx /// of bin edges @a xbinedges and @a ybinedges. Profile2DPtr bookProfile2D(const std::string& name, const std::initializer_list& xbinedges, const std::initializer_list& ybinedges, const std::string& title="", const std::string& xtitle="", const std::string& ytitle="", const std::string& ztitle=""); /// Book a 2D profile histogram with binning from a reference scatter. Profile2DPtr bookProfile2D(const std::string& name, const Scatter3D& refscatter, const std::string& title="", const std::string& xtitle="", const std::string& ytitle="", const std::string& ztitle=""); /// Book a 2D profile histogram, using the binnings in the reference data histogram. Profile2DPtr bookProfile2D(const std::string& name, const std::string& title="", const std::string& xtitle="", const std::string& ytitle="", const std::string& ztitle=""); /// Book a 2D profile histogram, using the binnings in the reference data histogram. /// /// The paper, dataset and x/y-axis IDs will be used to build the histo name in the HepData standard way. Profile2DPtr bookProfile2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId, const std::string& title="", const std::string& xtitle="", const std::string& ytitle="", const std::string& ztitle=""); //@} /// @name 2D scatter booking //@{ /// @brief Book a 2-dimensional data point set with the given name. /// /// @note Unlike histogram booking, scatter booking by default makes no /// attempt to use reference data to pre-fill the data object. If you want /// this, which is sometimes useful e.g. when the x-position is not really /// meaningful and can't be extracted from the data, then set the @a /// copy_pts parameter to true. This creates points to match the reference /// data's x values and errors, but with the y values and errors zeroed... /// assuming that there is a reference histo with the same name: if there /// isn't, an exception will be thrown. Scatter2DPtr bookScatter2D(const std::string& name, bool copy_pts=false, const std::string& title="", const std::string& xtitle="", const std::string& ytitle=""); /// @brief Book a 2-dimensional data point set, using the binnings in the reference data histogram. /// /// The paper, dataset and x/y-axis IDs will be used to build the histo name in the HepData standard way. /// /// @note Unlike histogram booking, scatter booking by default makes no /// attempt to use reference data to pre-fill the data object. If you want /// this, which is sometimes useful e.g. when the x-position is not really /// meaningful and can't be extracted from the data, then set the @a /// copy_pts parameter to true. This creates points to match the reference /// data's x values and errors, but with the y values and errors zeroed. Scatter2DPtr bookScatter2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId, bool copy_pts=false, const std::string& title="", const std::string& xtitle="", const std::string& ytitle=""); /// @brief Book a 2-dimensional data point set with equally spaced x-points in a range. /// /// The y values and errors will be set to 0. Scatter2DPtr bookScatter2D(const std::string& name, size_t npts, double lower, double upper, const std::string& title="", const std::string& xtitle="", const std::string& ytitle=""); /// @brief Book a 2-dimensional data point set based on provided contiguous "bin edges". /// /// The y values and errors will be set to 0. Scatter2DPtr bookScatter2D(const std::string& hname, const std::vector& binedges, const std::string& title, const std::string& xtitle, const std::string& ytitle); //@} public: - /// @name accessing options for this Anslysis instance. + /// @name accessing options for this Analysis instance. //@{ /// Get an option for this analysis instance as a string. std::string getOption(std::string optname) { if ( _options.find(optname) != _options.end() ) return _options.find(optname)->second; return ""; } /// Get an option for this analysis instance converted to a /// specific type (given by the specified @a def value). template T getOption(std::string optname, T def) { if (_options.find(optname) == _options.end()) return def; std::stringstream ss; ss << _options.find(optname)->second; T ret; ss >> ret; return ret; } //@} + /// @brief Book a CentralityProjection + /// + /// Using a SingleValueProjection, @a proj, giving the value of an + /// experimental observable to be used as a centrality estimator, + /// book a CentralityProjection based on the experimentally + /// measured pecentiles of this observable (as given by the + /// reference data for the @a calHistName histogram in the @a + /// calAnaName analysis. If a preloaded file with the output of a + /// run using the @a calAnaName analysis contains a valid + /// generated @a calHistName histogram, it will be used as an + /// optional percentile binning. Also if this preloaded file + /// contains a histogram with the name @a calHistName with an + /// appended "_IMP" This histogram will be used to add an optional + /// centrality percentile based on the generated impact + /// parameter. If @increasing is true, a low (high) value of @proj + /// is assumed to correspond to a more peripheral (central) event. + const CentralityProjection& + declareCentrality(const SingleValueProjection &proj, + string calAnaName, string calHistName, + const string projName, bool increasing = false); + /// @name Analysis object manipulation /// @todo Should really be protected: only public to keep BinnedHistogram happy for now... //@{ /// Multiplicatively scale the given counter, @a cnt, by factor @s factor. void scale(CounterPtr cnt, double factor); /// Multiplicatively scale the given counters, @a cnts, by factor @s factor. /// @note Constness intentional, if weird, to allow passing rvalue refs of smart ptrs (argh) /// @todo Use SFINAE for a generic iterable of CounterPtrs void scale(const std::vector& cnts, double factor) { for (auto& c : cnts) scale(c, factor); } /// @todo YUCK! template void scale(const CounterPtr (&cnts)[array_size], double factor) { // for (size_t i = 0; i < std::extent::value; ++i) scale(cnts[i], factor); for (auto& c : cnts) scale(c, factor); } /// Normalize the given histogram, @a histo, to area = @a norm. void normalize(Histo1DPtr histo, double norm=1.0, bool includeoverflows=true); /// Normalize the given histograms, @a histos, to area = @a norm. /// @note Constness intentional, if weird, to allow passing rvalue refs of smart ptrs (argh) /// @todo Use SFINAE for a generic iterable of Histo1DPtrs void normalize(const std::vector& histos, double norm=1.0, bool includeoverflows=true) { for (auto& h : histos) normalize(h, norm, includeoverflows); } /// @todo YUCK! template void normalize(const Histo1DPtr (&histos)[array_size], double norm=1.0, bool includeoverflows=true) { for (auto& h : histos) normalize(h, norm, includeoverflows); } /// Multiplicatively scale the given histogram, @a histo, by factor @s factor. void scale(Histo1DPtr histo, double factor); /// Multiplicatively scale the given histograms, @a histos, by factor @s factor. /// @note Constness intentional, if weird, to allow passing rvalue refs of smart ptrs (argh) /// @todo Use SFINAE for a generic iterable of Histo1DPtrs void scale(const std::vector& histos, double factor) { for (auto& h : histos) scale(h, factor); } /// @todo YUCK! template void scale(const Histo1DPtr (&histos)[array_size], double factor) { for (auto& h : histos) scale(h, factor); } /// Normalize the given histogram, @a histo, to area = @a norm. void normalize(Histo2DPtr histo, double norm=1.0, bool includeoverflows=true); /// Normalize the given histograms, @a histos, to area = @a norm. /// @note Constness intentional, if weird, to allow passing rvalue refs of smart ptrs (argh) /// @todo Use SFINAE for a generic iterable of Histo2DPtrs void normalize(const std::vector& histos, double norm=1.0, bool includeoverflows=true) { for (auto& h : histos) normalize(h, norm, includeoverflows); } /// @todo YUCK! template void normalize(const Histo2DPtr (&histos)[array_size], double norm=1.0, bool includeoverflows=true) { for (auto& h : histos) normalize(h, norm, includeoverflows); } /// Multiplicatively scale the given histogram, @a histo, by factor @s factor. void scale(Histo2DPtr histo, double factor); /// Multiplicatively scale the given histograms, @a histos, by factor @s factor. /// @note Constness intentional, if weird, to allow passing rvalue refs of smart ptrs (argh) /// @todo Use SFINAE for a generic iterable of Histo2DPtrs void scale(const std::vector& histos, double factor) { for (auto& h : histos) scale(h, factor); } /// @todo YUCK! template void scale(const Histo2DPtr (&histos)[array_size], double factor) { for (auto& h : histos) scale(h, factor); } /// Helper for counter division. /// /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target. void divide(CounterPtr c1, CounterPtr c2, Scatter1DPtr s) const; /// Helper for histogram division with raw YODA objects. /// /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target. void divide(const YODA::Counter& c1, const YODA::Counter& c2, Scatter1DPtr s) const; /// Helper for histogram division. /// /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target. void divide(Histo1DPtr h1, Histo1DPtr h2, Scatter2DPtr s) const; /// Helper for histogram division with raw YODA objects. /// /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target. void divide(const YODA::Histo1D& h1, const YODA::Histo1D& h2, Scatter2DPtr s) const; /// Helper for profile histogram division. /// /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target. void divide(Profile1DPtr p1, Profile1DPtr p2, Scatter2DPtr s) const; /// Helper for profile histogram division with raw YODA objects. /// /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target. void divide(const YODA::Profile1D& p1, const YODA::Profile1D& p2, Scatter2DPtr s) const; /// Helper for 2D histogram division. /// /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target. void divide(Histo2DPtr h1, Histo2DPtr h2, Scatter3DPtr s) const; /// Helper for 2D histogram division with raw YODA objects. /// /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target. void divide(const YODA::Histo2D& h1, const YODA::Histo2D& h2, Scatter3DPtr s) const; /// Helper for 2D profile histogram division. /// /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target. void divide(Profile2DPtr p1, Profile2DPtr p2, Scatter3DPtr s) const; /// Helper for 2D profile histogram division with raw YODA objects /// /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target. void divide(const YODA::Profile2D& p1, const YODA::Profile2D& p2, Scatter3DPtr s) const; /// Helper for histogram efficiency calculation. /// /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target. void efficiency(Histo1DPtr h1, Histo1DPtr h2, Scatter2DPtr s) const; /// Helper for histogram efficiency calculation. /// /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target. void efficiency(const YODA::Histo1D& h1, const YODA::Histo1D& h2, Scatter2DPtr s) const; /// Helper for histogram asymmetry calculation. /// /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target. void asymm(Histo1DPtr h1, Histo1DPtr h2, Scatter2DPtr s) const; /// Helper for histogram asymmetry calculation. /// /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target. void asymm(const YODA::Histo1D& h1, const YODA::Histo1D& h2, Scatter2DPtr s) const; /// Helper for converting a differential histo to an integral one. /// /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target. void integrate(Histo1DPtr h, Scatter2DPtr s) const; /// Helper for converting a differential histo to an integral one. /// /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target. void integrate(const Histo1D& h, Scatter2DPtr s) const; //@} public: /// List of registered analysis data objects const vector& analysisObjects() const { return _analysisobjects; } protected: /// @name Data object registration, retrieval, and removal //@{ /// Register a data object in the histogram system void addAnalysisObject(AnalysisObjectPtr ao); /// Get a data object from the histogram system template const std::shared_ptr getAnalysisObject(const std::string& name) const { foreach (const AnalysisObjectPtr& ao, analysisObjects()) { if (ao->path() == histoPath(name)) return dynamic_pointer_cast(ao); } throw Exception("Data object " + histoPath(name) + " not found"); } /// Get a data object from the histogram system (non-const) template std::shared_ptr getAnalysisObject(const std::string& name) { foreach (const AnalysisObjectPtr& ao, analysisObjects()) { if (ao->path() == histoPath(name)) return dynamic_pointer_cast(ao); } throw Exception("Data object " + histoPath(name) + " not found"); } /// Unregister a data object from the histogram system (by name) void removeAnalysisObject(const std::string& path); /// Unregister a data object from the histogram system (by pointer) void removeAnalysisObject(AnalysisObjectPtr ao); + /// Get all data object from the AnalysisHandler. + vector getAllData(bool includeorphans) const; + + /// Get a data object from another analysis (e.g. preloaded + /// calibration histogram). + /// Get a data object from the histogram system (non-const) + template + std::shared_ptr getAnalysisObject(const std::string & ananame, + const std::string& name) { + std::string path = "/" + ananame + "/" + name; + for ( AnalysisObjectPtr ao : getAllData(true) ) { + if ( ao->path() == path ) + return dynamic_pointer_cast(ao); + } + return std::shared_ptr(); + } /// Get a named Histo1D object from the histogram system const Histo1DPtr getHisto1D(const std::string& name) const { return getAnalysisObject(name); } /// Get a named Histo1D object from the histogram system (non-const) Histo1DPtr getHisto1D(const std::string& name) { return getAnalysisObject(name); } /// Get a Histo1D object from the histogram system by axis ID codes (non-const) const Histo1DPtr getHisto1D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const { return getAnalysisObject(makeAxisCode(datasetId, xAxisId, yAxisId)); } /// Get a Histo1D object from the histogram system by axis ID codes (non-const) Histo1DPtr getHisto1D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) { return getAnalysisObject(makeAxisCode(datasetId, xAxisId, yAxisId)); } /// Get a named Histo2D object from the histogram system const Histo2DPtr getHisto2D(const std::string& name) const { return getAnalysisObject(name); } /// Get a named Histo2D object from the histogram system (non-const) Histo2DPtr getHisto2D(const std::string& name) { return getAnalysisObject(name); } /// Get a Histo2D object from the histogram system by axis ID codes (non-const) const Histo2DPtr getHisto2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const { return getAnalysisObject(makeAxisCode(datasetId, xAxisId, yAxisId)); } /// Get a Histo2D object from the histogram system by axis ID codes (non-const) Histo2DPtr getHisto2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) { return getAnalysisObject(makeAxisCode(datasetId, xAxisId, yAxisId)); } /// Get a named Profile1D object from the histogram system const Profile1DPtr getProfile1D(const std::string& name) const { return getAnalysisObject(name); } /// Get a named Profile1D object from the histogram system (non-const) Profile1DPtr getProfile1D(const std::string& name) { return getAnalysisObject(name); } /// Get a Profile1D object from the histogram system by axis ID codes (non-const) const Profile1DPtr getProfile1D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const { return getAnalysisObject(makeAxisCode(datasetId, xAxisId, yAxisId)); } /// Get a Profile1D object from the histogram system by axis ID codes (non-const) Profile1DPtr getProfile1D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) { return getAnalysisObject(makeAxisCode(datasetId, xAxisId, yAxisId)); } /// Get a named Profile2D object from the histogram system const Profile2DPtr getProfile2D(const std::string& name) const { return getAnalysisObject(name); } /// Get a named Profile2D object from the histogram system (non-const) Profile2DPtr getProfile2D(const std::string& name) { return getAnalysisObject(name); } /// Get a Profile2D object from the histogram system by axis ID codes (non-const) const Profile2DPtr getProfile2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const { return getAnalysisObject(makeAxisCode(datasetId, xAxisId, yAxisId)); } /// Get a Profile2D object from the histogram system by axis ID codes (non-const) Profile2DPtr getProfile2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) { return getAnalysisObject(makeAxisCode(datasetId, xAxisId, yAxisId)); } /// Get a named Scatter2D object from the histogram system const Scatter2DPtr getScatter2D(const std::string& name) const { return getAnalysisObject(name); } /// Get a named Scatter2D object from the histogram system (non-const) Scatter2DPtr getScatter2D(const std::string& name) { return getAnalysisObject(name); } /// Get a Scatter2D object from the histogram system by axis ID codes (non-const) const Scatter2DPtr getScatter2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const { return getAnalysisObject(makeAxisCode(datasetId, xAxisId, yAxisId)); } /// Get a Scatter2D object from the histogram system by axis ID codes (non-const) Scatter2DPtr getScatter2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) { return getAnalysisObject(makeAxisCode(datasetId, xAxisId, yAxisId)); } //@} private: /// Name passed to constructor (used to find .info analysis data file, and as a fallback) string _defaultname; /// Pointer to analysis metadata object unique_ptr _info; /// Storage of all plot objects /// @todo Make this a map for fast lookup by path? vector _analysisobjects; - + /// @name Cross-section variables //@{ double _crossSection; bool _gotCrossSection; //@} /// The controlling AnalysisHandler object. AnalysisHandler* _analysishandler; /// Collection of cached refdata to speed up many autobookings: the /// reference data file should only be read once. mutable std::map _refdata; /// Options the (this instance of) the analysis map _options; /// The string of options. string _optstring; private: /// @name Utility functions //@{ /// Get the reference data for this paper and cache it. void _cacheRefData() const; //@} /// The assignment operator is private and must never be called. /// In fact, it should not even be implemented. Analysis& operator=(const Analysis&); }; } // Include definition of analysis plugin system so that analyses automatically see it when including Analysis.hh #include "Rivet/AnalysisBuilder.hh" /// @def DECLARE_RIVET_PLUGIN /// Preprocessor define to prettify the global-object plugin hook mechanism. #define DECLARE_RIVET_PLUGIN(clsname) Rivet::AnalysisBuilder plugin_ ## clsname /// @def DECLARE_ALIASED_RIVET_PLUGIN /// Preprocessor define to prettify the global-object plugin hook mechanism, with an extra alias name for this analysis. // #define DECLARE_ALIASED_RIVET_PLUGIN(clsname, alias) Rivet::AnalysisBuilder plugin_ ## clsname ## ( ## #alias ## ) #define DECLARE_ALIASED_RIVET_PLUGIN(clsname, alias) DECLARE_RIVET_PLUGIN(clsname)( #alias ) /// @def DEFAULT_RIVET_ANALYSIS_CONSTRUCTOR /// Preprocessor define to prettify the manky constructor with name string argument #define DEFAULT_RIVET_ANALYSIS_CONSTRUCTOR(clsname) clsname() : Analysis(# clsname) {} /// @def DEFAULT_RIVET_ANALYSIS_CTOR /// Slight abbreviation for DEFAULT_RIVET_ANALYSIS_CONSTRUCTOR #define DEFAULT_RIVET_ANALYSIS_CTOR(clsname) DEFAULT_RIVET_ANALYSIS_CONSTRUCTOR(clsname) #endif diff --git a/include/Rivet/AnalysisHandler.hh b/include/Rivet/AnalysisHandler.hh --- a/include/Rivet/AnalysisHandler.hh +++ b/include/Rivet/AnalysisHandler.hh @@ -1,265 +1,269 @@ // -*- C++ -*- #ifndef RIVET_RivetHandler_HH #define RIVET_RivetHandler_HH #include "Rivet/Config/RivetCommon.hh" #include "Rivet/Particle.hh" #include "Rivet/AnalysisLoader.hh" #include "Rivet/Tools/RivetYODA.hh" namespace Rivet { // Forward declaration and smart pointer for Analysis class Analysis; typedef std::shared_ptr AnaHandle; // Needed to make smart pointers compare equivalent in the STL set struct CmpAnaHandle { bool operator() (const AnaHandle& a, const AnaHandle& b) const { return a.get() < b.get(); } }; /// A class which handles a number of analysis objects to be applied to /// generated events. An {@link Analysis}' AnalysisHandler is also responsible /// for handling the final writing-out of histograms. class AnalysisHandler { public: /// @name Constructors and destructors. */ //@{ /// Preferred constructor, with optional run name. AnalysisHandler(const string& runname=""); /// @brief Destructor /// The destructor is not virtual, as this class should not be inherited from. ~AnalysisHandler(); //@} private: /// Get a logger object. Log& getLog() const; public: /// @name Run properties //@{ /// Get the name of this run. string runName() const { return _runname; } /// Get the number of events seen. Should only really be used by external /// steering code or analyses in the finalize phase. size_t numEvents() const { return _eventcounter.numEntries(); } /// @brief Access the sum of the event weights seen /// /// This is the weighted equivalent of the number of events. It should only /// be used by external steering code or analyses in the finalize phase. double sumW() const { return _eventcounter.sumW(); } /// Access to the sum of squared-weights double sumW2() const { return _eventcounter.sumW2(); } /// @brief Compatibility alias for sumOfWeights /// /// @deprecated Prefer sumW double sumOfWeights() const { return sumW(); } /// @brief Set the sum of weights /// /// This is useful if Rivet is steered externally and /// the analyses are run for a sub-contribution of the events /// (but of course have to be normalised to the total sum of weights) /// /// @todo What about the sumW2 term? That needs to be set coherently. Need a /// new version, with all three N,sumW,sumW2 numbers (or a counter) /// supplied. /// /// @deprecated Weight sums are no longer tracked this way... void setSumOfWeights(const double& sum) { //_sumOfWeights = sum; throw Error("Can't set sum of weights independently, since it's now tracked by a Counter. " "Please contact the Rivet authors if you need this."); } /// Is cross-section information required by at least one child analysis? /// @deprecated Should no-longer be an issue: does any generator not write the cross-section? bool needCrossSection() const; /// Whether the handler knows about a cross-section. /// @deprecated Should no-longer be an issue: does any generator not write the cross-section? bool hasCrossSection() const; /// Get the cross-section known to the handler. double crossSection() const { return _xs; } /// Set the cross-section for the process being generated. /// @todo What about the xsec uncertainty? Add a second, optional argument? AnalysisHandler& setCrossSection(double xs); /// Set the beam particles for this run AnalysisHandler& setRunBeams(const ParticlePair& beams) { _beams = beams; MSG_DEBUG("Setting run beams = " << beams << " @ " << sqrtS()/GeV << " GeV"); return *this; } /// Get the beam particles for this run, usually determined from the first event. const ParticlePair& beams() const { return _beams; } /// Get beam IDs for this run, usually determined from the first event. /// @deprecated Use standalone beamIds(ah.beams()), to clean AH interface PdgIdPair beamIds() const; /// Get energy for this run, usually determined from the first event. /// @deprecated Use standalone sqrtS(ah.beams()), to clean AH interface double sqrtS() const; /// Setter for _ignoreBeams void setIgnoreBeams(bool ignore=true); //@} /// @name Handle analyses //@{ /// Get a list of the currently registered analyses' names. std::vector analysisNames() const; /// Get the collection of currently registered analyses. const std::set& analyses() const { return _analyses; } /// Get a registered analysis by name. const AnaHandle analysis(const std::string& analysisname) const; /// Add an analysis to the run list by object AnalysisHandler& addAnalysis(Analysis* analysis); /// @brief Add an analysis to the run list using its name. /// /// The actual Analysis to be used will be obtained via /// AnalysisLoader::getAnalysis(string). If no matching analysis is found, /// no analysis is added (i.e. the null pointer is checked and discarded. AnalysisHandler& addAnalysis(const std::string& analysisname); /// @brief Add analyses to the run list using their names. /// /// The actual {@link Analysis}' to be used will be obtained via /// AnalysisHandler::addAnalysis(string), which in turn uses /// AnalysisLoader::getAnalysis(string). If no matching analysis is found /// for a given name, no analysis is added, but also no error is thrown. AnalysisHandler& addAnalyses(const std::vector& analysisnames); /// Remove an analysis from the run list using its name. AnalysisHandler& removeAnalysis(const std::string& analysisname); /// Remove analyses from the run list using their names. AnalysisHandler& removeAnalyses(const std::vector& analysisnames); //@} /// @name Main init/execute/finalise //@{ /// Initialize a run, with the run beams taken from the example event. void init(const GenEvent& event); /// @brief Analyze the given \a event by reference. /// /// This function will call the AnalysisBase::analyze() function of all /// included analysis objects. void analyze(const GenEvent& event); /// @brief Analyze the given \a event by pointer. /// /// This function will call the AnalysisBase::analyze() function of all /// included analysis objects, after checking the event pointer validity. void analyze(const GenEvent* event); /// Finalize a run. This function calls the AnalysisBase::finalize() /// functions of all included analysis objects. void finalize(); //@} /// @name Histogram / data object access //@{ /// Add a vector of analysis objects to the current state. void addData(const std::vector& aos); /// Read analysis plots into the histo collection (via addData) from the named file. void readData(const std::string& filename); /// Get all analyses' plots as a vector of analysis objects. - std::vector getData() const; + std::vector getData(bool includeorphans = false) const; /// Write all analyses' plots (via getData) to the named file. void writeData(const std::string& filename) const; //@} private: /// The collection of Analysis objects to be used. set _analyses; + /// A vector of pre-loaded object which do not have a valid + /// Analysis plugged in. + vector _orphanedPreloads; + /// @name Run properties //@{ /// Run name std::string _runname; /// Event counter Counter _eventcounter; /// Cross-section known to AH double _xs, _xserr; /// Beams used by this run. ParticlePair _beams; /// Flag to check if init has been called bool _initialised; /// Flag whether input event beams should be ignored in compatibility check bool _ignoreBeams; //@} private: /// The assignment operator is private and must never be called. /// In fact, it should not even be implemented. AnalysisHandler& operator=(const AnalysisHandler&); /// The copy constructor is private and must never be called. In /// fact, it should not even be implemented. AnalysisHandler(const AnalysisHandler&); }; } #endif diff --git a/include/Rivet/Makefile.am b/include/Rivet/Makefile.am --- a/include/Rivet/Makefile.am +++ b/include/Rivet/Makefile.am @@ -1,171 +1,179 @@ ## Internal headers - not to be installed nobase_dist_noinst_HEADERS = ## Public headers - to be installed nobase_pkginclude_HEADERS = ## Rivet interface nobase_pkginclude_HEADERS += \ Rivet.hh \ Run.hh \ Event.hh \ ParticleBase.hh \ Particle.fhh Particle.hh \ Jet.fhh Jet.hh \ Projection.fhh Projection.hh \ ProjectionApplier.hh \ ProjectionHandler.hh \ Analysis.hh \ AnalysisHandler.hh \ AnalysisInfo.hh \ AnalysisBuilder.hh \ AnalysisLoader.hh ## Build config stuff nobase_pkginclude_HEADERS += \ Config/RivetCommon.hh \ Config/RivetConfig.hh \ Config/BuildOptions.hh ## Projections nobase_pkginclude_HEADERS += \ + Projections/AliceCommon.hh \ Projections/AxesDefinition.hh \ Projections/Beam.hh \ Projections/BeamThrust.hh \ Projections/CentralEtHCM.hh \ + Projections/CentralityProjection.hh \ Projections/ChargedFinalState.hh \ Projections/ChargedLeptons.hh \ Projections/ConstLossyFinalState.hh \ Projections/DirectFinalState.hh \ Projections/DISFinalState.hh \ Projections/DISKinematics.hh \ Projections/DISLepton.hh \ Projections/DressedLeptons.hh \ Projections/FastJets.hh \ Projections/FinalPartons.hh \ Projections/FinalState.hh \ Projections/FoxWolframMoments.hh \ Projections/FParameter.hh \ + Projections/GeneratedPercentileProjection.hh \ Projections/HadronicFinalState.hh \ Projections/HeavyHadrons.hh \ Projections/Hemispheres.hh \ Projections/IdentifiedFinalState.hh \ + Projections/ImpactParameterProjection.hh \ Projections/IndirectFinalState.hh \ Projections/InitialQuarks.hh \ Projections/InvMassFinalState.hh \ Projections/JetAlg.hh \ Projections/JetShape.hh \ Projections/LeadingParticlesFinalState.hh \ Projections/LossyFinalState.hh \ Projections/MergedFinalState.hh \ Projections/MissingMomentum.hh \ Projections/NeutralFinalState.hh \ Projections/NonHadronicFinalState.hh \ Projections/NonPromptFinalState.hh \ Projections/ParisiTensor.hh \ Projections/ParticleFinder.hh \ Projections/PartonicTops.hh \ Projections/PrimaryHadrons.hh \ + Projections/PrimaryParticles.hh \ Projections/PromptFinalState.hh \ + Projections/SingleValueProjection.hh \ Projections/SmearedParticles.hh \ Projections/SmearedJets.hh \ Projections/SmearedMET.hh \ Projections/Sphericity.hh \ Projections/Spherocity.hh \ Projections/TauFinder.hh \ Projections/Thrust.hh \ Projections/TriggerCDFRun0Run1.hh \ Projections/TriggerCDFRun2.hh \ Projections/TriggerUA5.hh \ Projections/UnstableFinalState.hh \ + Projections/UserCentEstimate.hh \ Projections/VetoedFinalState.hh \ Projections/VisibleFinalState.hh \ Projections/WFinder.hh \ Projections/ZFinder.hh ## Meta-projection convenience headers nobase_pkginclude_HEADERS += \ Projections/FinalStates.hh \ Projections/Smearing.hh ## Analysis base class headers # TODO: Move to Rivet/AnalysisTools header dir? nobase_pkginclude_HEADERS += \ Analyses/MC_ParticleAnalysis.hh \ Analyses/MC_JetAnalysis.hh \ Analyses/MC_JetSplittings.hh ## Tools nobase_pkginclude_HEADERS += \ + Tools/AliceCommon.hh \ Tools/BeamConstraint.hh \ Tools/BinnedHistogram.hh \ Tools/CentralityBinner.hh \ Tools/Cmp.fhh \ Tools/Cmp.hh \ Tools/Cutflow.hh \ Tools/Cuts.fhh \ Tools/Cuts.hh \ Tools/Exceptions.hh \ Tools/JetUtils.hh \ Tools/Logging.hh \ Tools/Random.hh \ Tools/ParticleBaseUtils.hh \ Tools/ParticleIdUtils.hh \ Tools/ParticleUtils.hh \ Tools/ParticleName.hh \ Tools/PrettyPrint.hh \ Tools/RivetPaths.hh \ Tools/RivetSTL.hh \ Tools/RivetFastJet.hh \ Tools/RivetHepMC.hh \ Tools/RivetYODA.hh \ Tools/RivetMT2.hh \ Tools/SmearingFunctions.hh \ Tools/MomentumSmearingFunctions.hh \ Tools/ParticleSmearingFunctions.hh \ Tools/JetSmearingFunctions.hh \ Tools/TypeTraits.hh \ Tools/Utils.hh \ Tools/Nsubjettiness/AxesDefinition.hh \ Tools/Nsubjettiness/ExtraRecombiners.hh \ Tools/Nsubjettiness/MeasureDefinition.hh \ Tools/Nsubjettiness/Njettiness.hh \ Tools/Nsubjettiness/NjettinessPlugin.hh \ Tools/Nsubjettiness/Nsubjettiness.hh \ Tools/Nsubjettiness/TauComponents.hh \ Tools/Nsubjettiness/XConePlugin.hh nobase_dist_noinst_HEADERS += \ Tools/osdir.hh ## Maths nobase_pkginclude_HEADERS += \ Math/Matrices.hh \ Math/Vector3.hh \ Math/VectorN.hh \ Math/MatrixN.hh \ Math/MatrixDiag.hh \ Math/MathHeader.hh \ Math/Vectors.hh \ Math/LorentzTrans.hh \ Math/Matrix3.hh \ Math/MathUtils.hh \ Math/Vector4.hh \ Math/Math.hh \ Math/Units.hh \ Math/Constants.hh \ Math/eigen/util.h \ Math/eigen/regressioninternal.h \ Math/eigen/regression.h \ Math/eigen/vector.h \ Math/eigen/ludecompositionbase.h \ Math/eigen/ludecomposition.h \ Math/eigen/linearsolver.h \ Math/eigen/linearsolverbase.h \ Math/eigen/matrix.h \ Math/eigen/vectorbase.h \ Math/eigen/projective.h \ Math/eigen/matrixbase.h diff --git a/include/Rivet/Projections/AliceCommon.hh b/include/Rivet/Projections/AliceCommon.hh new file mode 100644 --- /dev/null +++ b/include/Rivet/Projections/AliceCommon.hh @@ -0,0 +1,428 @@ +/** + * @file AliceCommon.hh + * @author Christian Holm Christensen + * @date Fri Aug 24 09:50:33 2018 + * + * @brief Common AlICE specific projections + * @ingroup alice_rivet + */ +#ifndef PROJECTIONS_ALICECOMMON_HH +#define PROJECTIONS_ALICECOMMON_HH +#include +#include +#include +#include +#include + +namespace Rivet +{ + namespace ALICE + { + /** + * Template for ALICE V0 multiplicity projection. Which + * acceptance to look in depends on the template argument @a MODE: + * + * - @c MODE=-1 Check the V0-C acceptance (@f$-3.7<\eta<-1.7@f$) + * - @c MODE=+1 Check the V0-A acceptance (@f$+2.8<\eta<+5.1@f$) + * - @c MODE=0 Check both V0-A and -C acceptances (sum) + * + * + * @ingroup alice_rivet + */ + template + class V0Multiplicity : public SingleValueProjection + { + public: + /** + * Constructor + */ + V0Multiplicity() + : SingleValueProjection() + { + setName("ALICE::V0Multiplicity"); + Cut cut; + if (MODE < 0) cut = V0Cacceptance; + else if (MODE > 0) cut = V0Aacceptance; + else cut = (V0Aacceptance || V0Cacceptance); + // Declare our projection. Note, the cuts stipulate charged + // particles, so we just use a final state (rather than + // charged-final state) projection here. + const FinalState fs(cut); + this->declare(fs, "FinalState"); + } + /** + * Destructor + */ + virtual ~V0Multiplicity() {} + /** + * Do the projection. Sums number of charged final state + * particles within the acceptances of the specified V0 + * sub-detectors. + * + * @param e Event to project from + */ + virtual void project(const Event& e) + { + clear(); + set(apply(e,"FinalState").particles().size()); + } + /** + * Clone this projection + * + * @return New wrapped pointer to object of this class + */ + virtual std::unique_ptr clone() const + { + return std::unique_ptr(new V0Multiplicity(*this)); + } + /** + * Compare to another projection + * + * @param p Projection to compare against + */ + virtual int compare(const Projection& p) const + { + return dynamic_cast*>(&p) ? + EQUIVALENT : UNDEFINED; + } + }; + //---------------------------------------------------------------- + /** + * Convinience type definition + * + * @ingroup alice_rivet + */ + typedef V0Multiplicity<-1> V0AMultiplicity; + /** + * Convinience type definition + * + * @ingroup alice_rivet + */ + typedef V0Multiplicity<+1> V0CMultiplicity; + /** + * Convinience type definition + * + * @ingroup alice_rivet + */ + typedef V0Multiplicity<0> V0MMultiplicity; + //---------------------------------------------------------------- + /** + * Template for ALICE CL multiplicity projection. Which + * acceptance to look in depends on the template argument @a INNER: + * + * - @c INNER=true Check the inner SPD layer + * - @c INNER=false Check the outer SPD layer + * + * @ingroup alice_rivet + * + */ + template + class CLMultiplicity : public SingleValueProjection + { + public: + /** + * Constructor + */ + CLMultiplicity() + : SingleValueProjection() + { + setName("ALICE::CLMultiplicity"); + Cut cut; + if (INNER) cut = CL0acceptance; + else cut = CL1acceptance; + // Declare our projection. Note, the cuts stipulate charged + // particles, so we just use a final state (rather than + // charged-final state) projection here. + const FinalState fs(cut); + this->declare(fs, "FinalState"); + } + /** + * Destructor + */ + virtual ~CLMultiplicity() {} + /** + * Do the projection. Sums number of charged final state + * particles within the acceptances of the specified CL + * sub-detectors. + * + * @param e Event to project from + */ + virtual void project(const Event& e) + { + clear(); + set(apply(e,"FinalState").particles().size()); + } + /** + * Clone this projection + * + * @return New wrapped pointer to object of this class + */ + virtual std::unique_ptr clone() const + { + return std::unique_ptr(new CLMultiplicity(*this)); + } + /** + * Compare to another projection + * + * @param p Projection to compare against + */ + virtual int compare(const Projection& p) const + { + return dynamic_cast*>(&p) ? + EQUIVALENT : UNDEFINED; + } + }; + //---------------------------------------------------------------- + /** + * Convinience type definition + * + * @ingroup alice_rivet + */ + typedef CLMultiplicity CL0Multiplicity; + /** + * Convinience type definition + * + * @ingroup alice_rivet + */ + typedef CLMultiplicity CL1Multiplicity; + //================================================================ + /** + * A template of ALICE V0-based triggers. + * + * - @c MODE=-1 Check in the V0-C acceptance (@f$-3.7<\eta<-1.7@f$) + * - @c MODE=+1 Check in the V0-A acceptance (@f$+2.8<\eta<+5.1@f$) + * - @c MODE=0 Check in both V0-A and -C acceptances (V0-OR) + * + * @ingroup alice_rivet + */ + template + class V0Trigger : public TriggerProjection + { + public: + /** + * Constructor + */ + V0Trigger() + : TriggerProjection() + { + setName("ALICE::V0Trigger"); + // Declare our projection. Note, the cuts stipulate charged + // particles, so we just use a final state (rather than + // charged-final state) projection here. + const V0Multiplicity fs; + this->declare(fs, "FinalState"); + } + /** + * Destructor + */ + virtual ~V0Trigger() {} + /** + * Do the projection. Checks if the number of projected + * particles is larger than 0 + * + * @param e Event to project from + */ + virtual void project(const Event& e) + { + fail(); // Assume failure + if (apply>(e,"FinalState")() > 0) pass(); + } + /** + * Clone this projection + * + * @return New wrapped pointer to object of this class + */ + virtual std::unique_ptr clone() const + { + return std::unique_ptr(new V0Trigger(*this)); + } + /** + * Compare to projections. + * + * @param p Projection to compare to. + * + * @return true (EQUIVALENT) if the projection @a p is of the same + * type as this. + */ + virtual int compare(const Projection& p) const + { + return dynamic_cast*>(&p) ? + EQUIVALENT : UNDEFINED; + } + }; + //---------------------------------------------------------------- + /** + * Convinience type definition + * + * @ingroup alice_rivet + */ + // typedef V0Trigger<-1> V0ATrigger; + using V0ATrigger = V0Trigger<-1>; + /** + * Convinience type definition + * + * @ingroup alice_rivet + */ + // typedef V0Trigger<+1> V0CTrigger; + using V0CTrigger = V0Trigger<+1>; + /** + * Convinience type definition + * + * @ingroup alice_rivet + */ + // typedef V0Trigger<0> V0OrTrigger; + using V0OrTrigger = V0Trigger<0>; + //---------------------------------------------------------------- + /** + * A trigger projetion for the ALICE V0-AND (a.k.a. CINT7) + * requirement. + */ + class V0AndTrigger : public TriggerProjection + { + public: + /** + * Constructor + */ + V0AndTrigger() + : TriggerProjection() + { + const V0ATrigger v0a; + const V0CTrigger v0c; + this->declare(v0a, "V0A"); + this->declare(v0c, "V0C"); + } + /** + * Destructor + */ + virtual ~V0AndTrigger() {} + /** + * Do the projection. Checks if the numbers of projected + * particles on both sides, are larger than 0 + * + * @param e Event to project from + */ + virtual void project(const Event& e) + { + fail(); // Assume failure + if (apply(e,"V0A")() && + apply(e,"V0C")()) pass(); + } + /** + * Compare to projections. + * + * @param p Projection to compare to. + * + * @return true (EQUIVALENT) if the projection @a p is of the same + * type as this. + */ + virtual int compare(const Projection& p) const + { + return dynamic_cast(&p) ? + EQUIVALENT : UNDEFINED; + } + /** + * Clone this projection + * + * @return New wrapped pointer to object of this class + */ + virtual std::unique_ptr clone() const + { + return std::unique_ptr(new V0AndTrigger(*this)); + } + }; + + + //================================================================ + /** + * Standard ALICE primary particle definition. + * also + * Primary particle definition according to public note + * ALICE-PUBLIC-2017-005 + * + * @ingroup alice_rivet + */ + class PrimaryParticles : public Rivet::PrimaryParticles + { + public: + PrimaryParticles(const Cut& c=Cuts::open()) + : Rivet::PrimaryParticles({},c) + {} + /** + * Compare to projections. + * + * @param p Projection to compare to. + * + * @return true (EQUIVALENT) if the projection @a p is of the same + * type as this, if the cuts are equal, and that the list of PDG + * IDs are the same. + */ + virtual int compare(const Projection& p) const + { + const PrimaryParticles* o = dynamic_cast(&p); + if (_cuts != o->_cuts) return UNDEFINED; + return mkPCmp(*o,"PrimaryParticles"); + } + /** + * Clone this projection + * + * @return New wrapped pointer to object of this class + */ + virtual std::unique_ptr clone() const + { + return std::unique_ptr(new PrimaryParticles(*this)); + } + protected: + /** + * Check PDG ID of particle @a p is in the list of accepted + * primaries. + * + * @param p Particle to investigate. + * + * @return true if the particle PDG ID is in the list of known + * primary PDG IDs. + * + * @note We explicitly override this to allow for nuclei, and we + * explicitly check for a specific set of particles (and + * anti-particles). This means we do not use the base class + * list of particles. Therefore, we also need to override the + * compare method. + */ + bool isPrimaryPID(const HepMC::GenParticle* p) const + { + int pdg = PID::abspid(p->pdg_id()); + // Check for nuclus + if (pdg > 1000000000) return true; + + switch (pdg) { + case Rivet::PID::MUON: + case Rivet::PID::ELECTRON: + case Rivet::PID::GAMMA: + case Rivet::PID::PIPLUS: + case Rivet::PID::KPLUS: + case Rivet::PID::K0S: + case Rivet::PID::K0L: + case Rivet::PID::PROTON: + case Rivet::PID::NEUTRON: + case Rivet::PID::LAMBDA: + case Rivet::PID::SIGMAMINUS: + case Rivet::PID::SIGMAPLUS: + case Rivet::PID::XIMINUS: + case Rivet::PID::XI0: + case Rivet::PID::OMEGAMINUS: + case Rivet::PID::NU_E: + case Rivet::PID::NU_MU: + case Rivet::PID::NU_TAU: + return true; + } + return false; + } + }; + } +} +#endif +// +// EOF +// + + + diff --git a/include/Rivet/Projections/CentralityProjection.hh b/include/Rivet/Projections/CentralityProjection.hh new file mode 100644 --- /dev/null +++ b/include/Rivet/Projections/CentralityProjection.hh @@ -0,0 +1,103 @@ +// -*- C++ -*- +#ifndef RIVET_CENTRALITYPROJECTION_HH +#define RIVET_CENTRALITYPROJECTION_HH + +#include "Rivet/Projections/PercentileProjection.hh" +#include "Rivet/Tools/RivetYODA.hh" +#include + +namespace Rivet { + +/** + @brief CentralityProjection is used together with the + percentile-based analysis objects Percentile and PercentileXaxsis. + + The interior actually defines several different centrality + estimates: the centrality observable used in the experiment with a + reference calibration ("REF"); the same but using a user-defined + calibration done with the corresponding minimum bias analysis + ("GEN"); a centrality based on the impact parameter reported in + HepMC::HeavyIon::impact_parameter, using a calibration histogram + generated with the same minimum bias analysis ("IMP"). For HepMC3 + it may optionally also include a direct report from the generator + about the centrality, if available in HepMC::HeavyIon::centrality + ("RAW"), and a user-defined generated centrality estimate + communicated via the HepMC::HeavyIon::user_cent_estimate ("USR"). + + @author Leif Lönnblad + +*/ + +class CentralityProjection: public SingleValueProjection { + +public: + + /// Default constructor. + CentralityProjection() {} + + + DEFAULT_RIVET_PROJ_CLONE(CentralityProjection); + + /// @BRIEF Add a new centality estimate. + /// + /// The SingelValueProjection, @a p, should return a value between 0 + /// and 100, and the @a pname should be one of "REF", "GEN", "IMP", + /// "USR", or "RAW", as described above. + void add(const SingleValueProjection & p, string pname) { + _projNames.push_back(pname); + declare(p, pname); + } + + /// Perform all internal projections. + void project(const Event& e) { + _values.clear(); + for ( string pname : _projNames ) + _values.push_back(apply(e, pname)()); + if ( !_values.empty() ) set(_values[0]); + } + + /// Cheek if no internal projections have been added. + bool empty() const { + return _projNames.empty(); + } + + /// Return the percentile of the @a i'th projection. + /// + /// Note that operator() will return the zero'th projection. + double operator[](int i) const { + return _values[i]; + } + + // Standard comparison function. + int compare(const Projection& p) const { + const CentralityProjection* other = + dynamic_cast(&p); + if (other->_projNames.size() == 0) return UNDEFINED; + for (string pname : _projNames) { + bool hasPname = true; + for (string p2name : other->_projNames){ + if (pname != p2name) hasPname = false; + } + if (!hasPname) return UNDEFINED; + } + return EQUIVALENT; + } + + /// THe list of names of the internal projections. + vector projections() const { + return _projNames; + } + +private: + + /// THe list of names of the internal projections. + vector _projNames; + + /// The list of percentiles resulting from the last projection. + vector _values; + +}; + +} + +#endif diff --git a/include/Rivet/Projections/GeneratedPercentileProjection.hh b/include/Rivet/Projections/GeneratedPercentileProjection.hh new file mode 100644 --- /dev/null +++ b/include/Rivet/Projections/GeneratedPercentileProjection.hh @@ -0,0 +1,41 @@ +// -*- C++ -*- +#ifndef RIVET_GENERATEDPERCENTILEPROJECTION_HH +#define RIVET_GENERATEDPERCENTILEPROJECTION_HH + +#include "Rivet/Projections/SingleValueProjection.hh" + + +namespace Rivet { + +class GeneratedPercentileProjection: public SingleValueProjection { + +public: + + GeneratedPercentileProjection() { + setName("GeneratedPercentileProjection"); + } + + /// Clone on the heap. + DEFAULT_RIVET_PROJ_CLONE(GeneratedPercentileProjection); + +protected: + + void project(const Event& e) { + clear(); +#if HEPMC_VERSION_CODE >= 3000000 + const HepMC::HeavyIon * hi = e.genEvent()->heavy_ion(); + if ( hi && hi->centrality >= 0.0 ) + set(hi->centrality*100.0); +#endif + } + + int compare(const Projection& p) const { + return 0; + } + +}; + +} + +#endif + diff --git a/include/Rivet/Projections/ImpactParameterProjection.hh b/include/Rivet/Projections/ImpactParameterProjection.hh new file mode 100644 --- /dev/null +++ b/include/Rivet/Projections/ImpactParameterProjection.hh @@ -0,0 +1,39 @@ +// -*- C++ -*- +#ifndef RIVET_IMPACTPARAMETERPROJECTION_HH +#define RIVET_IMPACTPARAMETERPROJECTION_HH + +#include "Rivet/Projections/SingleValueProjection.hh" + + +namespace Rivet { + +class ImpactParameterProjection: public SingleValueProjection { + +public: + + ImpactParameterProjection() { + setName("ImpactParameterProjection"); + } + + /// Clone on the heap. + DEFAULT_RIVET_PROJ_CLONE(ImpactParameterProjection); + +protected: + + void project(const Event& e) { + clear(); + const HepMC::HeavyIon * hi = e.genEvent()->heavy_ion(); + if ( hi && hi->impact_parameter() >= 0.0 ) + set(hi->impact_parameter()); + } + + int compare(const Projection& p) const { + return 0; + } + +}; + +} + +#endif + diff --git a/include/Rivet/Projections/PrimaryParticles.hh b/include/Rivet/Projections/PrimaryParticles.hh new file mode 100644 --- /dev/null +++ b/include/Rivet/Projections/PrimaryParticles.hh @@ -0,0 +1,142 @@ +// -*- C++ -*- +/** + * @file PrimaryParticles.hh + * @author Christian Holm Christensen + * @date Mon Aug 27 08:46:34 2018 + * @brief Primary particle definition based on PDG code. + */ +#ifndef RIVET_PrimaryParticles_HH +#define RIVET_PrimaryParticles_HH + +#include +#include + +namespace Rivet +{ + /// @brief Project out primary particles according to definition. + /// A Rivet projection that mimics an experimental primary partcile + /// definition by projecting out according to particle id. + /// The projection can be further specialized to accomodate + /// specific experimental definitions. + // + class PrimaryParticles : public ParticleFinder + { + public: + /** + * Constructor + * + * @param cuts Normal particle cuts + * @param pdgIds List of PDG IDs which are considered primary. + * + * @todo Instead of using a real vector use an initializer list - + * more flexible. Also, do not provide a default for the PDG IDs + * (even if empty) as this will force the user to actively select + * which codes are primary. + */ + PrimaryParticles(std::initializer_list pdgIds, + const Cut& c=Cuts::open()) : + ParticleFinder(c), _pdgIds(pdgIds) { + setName("PrimaryParticles"); + } + // Clone on the heap. + DEFAULT_RIVET_PROJ_CLONE(PrimaryParticles); + + /** + * Copy constructor + * + * @param other Other object to copy from + */ + PrimaryParticles(const PrimaryParticles& other) : + ParticleFinder(other), _pdgIds(other._pdgIds) { + } + /** + * Compare to projections. + * + * @param p Projection to compare to. + * + * @return true (EQUIVALENT) if the projection @a p is of the same + * type as this, if the cuts are equal, and that the list of PDG + * IDs are the same. + */ + virtual int compare(const Projection& p) const + { + const PrimaryParticles* other = dynamic_cast(&p); + if (!other) return UNDEFINED; + if (_cuts != other->_cuts || _pdgIds != other->_pdgIds) return UNDEFINED; + return EQUIVALENT; + + } + protected: + /** + * Do the projection. + * + * @param e Event to project from + */ + virtual void project(const Event& e); + + /** + * @{ + * @name Internally used member functions + */ + /** + * Check if the particle is a priamry. + * + * @param p Pointer to a HepMC particle + * + * @return true if the particle @a p is considered primary + */ + virtual bool isPrimary(const HepMC::GenParticle* p) const; + /** + * Check if the particle should be ignored by the status code of + * the particle. + */ + virtual bool isIgnored(const HepMC::GenParticle* p) const; + /** + * Check PDG ID of particle @a p is in the list of accepted + * primaries. + * + * @param p Particle to investigate. + * + * @return true if the particle PDG ID is in the list of known + * primary PDG IDs. + */ + virtual bool isPrimaryPID(const HepMC::GenParticle* p) const; + /* + * Check if a particle @a p has decayed. + * + * @param p Pointer to HepMC particle + * + * @return true if the particle has decayed according to the + * status flag of the particle @a p + */ + virtual bool hasDecayed(const HepMC::GenParticle* p) const; + /** + * Check if a particle is a beam (remnant) particle. + * + * @param p Particle to check + * + * @return true if the particle @a p is a (remnant) beam particle + */ + virtual bool isBeam(const HepMC::GenParticle* p) const; + /* + * Get the immediate ancestor of a particle. + * + * @param p Particle for which to get the immediate ancestor + * + * @return Pointer to immediate ancestor or null if there's no ancestor. + */ + const HepMC::GenParticle* ancestor(const HepMC::GenParticle* p) const; + /* + * Get the immediate ancestor of a particle, which is @e not an + * ignored particle. + * + * @param p Particle for which to get the immediate ancestor + * + * @return Pointer to immediate ancestor or null if there's no ancestor. + */ + const HepMC::GenParticle* ancestor(const HepMC::GenParticle* p, bool) const; + /** Particle types to test for */ + std::vector _pdgIds; + }; +} +#endif diff --git a/include/Rivet/Projections/SingleValueProjection.hh b/include/Rivet/Projections/SingleValueProjection.hh new file mode 100644 --- /dev/null +++ b/include/Rivet/Projections/SingleValueProjection.hh @@ -0,0 +1,63 @@ +// -*- C++ -*- +#ifndef RIVET_SINGLEVALUEPROJECTION_HH +#define RIVET_SINGLEVALUEPROJECTION_HH + +#include "Rivet/Projection.hh" + +namespace Rivet { + +/** @brief Base class for projections returning a single floating + point value. + + @author Leif Lönnblad + + Project an event down to a single floating point value accessible + through the operator() function. + +*/ +class SingleValueProjection: public Projection { + +public: + + /// The default constructor. + SingleValueProjection() : _value(-1.0), _isSet(false) { + setName("SingleValueProjection"); + } + + /// Returns true if the value has been set. + bool isSet() const { + return _isSet; + } + + /// Return the single value. + double operator()() const { + return _value; + } + +protected: + + /// Set the value. + void set(double v) { + _value = v; + _isSet = true; + } + + + /// Unset the value. + void clear() { + _value = -1.0; + _isSet = false; + } + +private: + + double _value; + + bool _isSet; + +}; + +} + +#endif + diff --git a/include/Rivet/Projections/UserCentEstimate.hh b/include/Rivet/Projections/UserCentEstimate.hh new file mode 100644 --- /dev/null +++ b/include/Rivet/Projections/UserCentEstimate.hh @@ -0,0 +1,41 @@ +// -*- C++ -*- +#ifndef RIVET_USERCENTESTIMATE_HH +#define RIVET_USERCENTESTIMATE_HH + +#include "Rivet/Projections/SingleValueProjection.hh" + + +namespace Rivet { + +class UserCentEstimate: public SingleValueProjection { + +public: + + UserCentEstimate() { + setName("UserCentEstimate"); + } + + /// Clone on the heap. + DEFAULT_RIVET_PROJ_CLONE(UserCentEstimate); + +protected: + + void project(const Event& e) { + clear(); +#if HEPMC_VERSION_CODE >= 3000000 + const HepMC::HeavyIon * hi = e.genEvent()->heavy_ion(); + if ( hi && hi->user_cent_estimate >= 0.0 ) + set(hi->centrality*100.0); +#endif + } + + int compare(const Projection& p) const { + return 0; + } + +}; + +} + +#endif + diff --git a/include/Rivet/Tools/AliceCommon.hh b/include/Rivet/Tools/AliceCommon.hh new file mode 100644 --- /dev/null +++ b/include/Rivet/Tools/AliceCommon.hh @@ -0,0 +1,85 @@ +/** + * @file AliceCommon.hh + * @author Christian Holm Christensen + * @date Thu Apr 27 17:07:35 2017 + * + * @brief Constants, etc. used in the ALICE Rivet code + * @copyright GNU Lesser Public License. + * + * @ingroup alice_rivet + */ +#ifndef TOOLS_ALICECOMMON_HH +#define TOOLS_ALICECOMMON_HH +#include +#include + +namespace Rivet +{ + /** + * @defgroup alice_rivet ALICE specific code in Rivet + * + * This include projections to emulate trigger conditions, + * centrality, and selection of primary particles. + */ + /** + * Namespace for ALICE specific core code + * + * @ingroup alice_rivet + */ + namespace ALICE + { + /** + * The acceptance cut for the V0A + * + * @ingroup alice_rivet + */ + const Cut V0Aacceptance = (Cuts::etaIn(+2.8,+5.1)&&(Cuts::abscharge3 > 0)); + /** + * The acceptance cut for the V0C + * + * @ingroup alice_rivet + */ + const Cut V0Cacceptance = (Cuts::etaIn(-3.7,-1.7)&&(Cuts::abscharge3 > 0)); + /** + * The acceptance cut for clusters on layer 0 of the SPD + * + * @ingroup alice_rivet + */ + const Cut CL0acceptance = (Cuts::etaIn(-2.0,2.0) && (Cuts::abscharge3 > 0)); + /** + * The acceptance cut for clusters on layer 1 of the SPD + * + * @ingroup alice_rivet + */ + const Cut CL1acceptance = (Cuts::etaIn(-1.4,1.4) && (Cuts::abscharge3 > 0)); + /** + * The acceptance cut for mid-rapidity + * + * @ingroup alice_rivet + */ + const Cut Eta1acceptance = (Cuts::etaIn(-1,1) && (Cuts::abscharge3 > 0)); + /** + * The acceptance cut for SPD FASTOR + * + * @ingroup alice_rivet + */ + const Cut FASTORacceptance = CL0acceptance; +#if 0 + /** + * Pb identification + * + * @ingroup alice_rivet + */ + const int PbId = (1000000000 + // ION identifier + 0*10000000 + // # strange quarks + 82*10000 + // atomic number + 208*10 + // atomic weight + 0*1); // Isomer number +#endif + } +} + +#endif +// +// EOF +// diff --git a/src/Core/Analysis.cc b/src/Core/Analysis.cc --- a/src/Core/Analysis.cc +++ b/src/Core/Analysis.cc @@ -1,899 +1,1006 @@ // -*- C++ -*- #include "Rivet/Config/RivetCommon.hh" #include "Rivet/Analysis.hh" #include "Rivet/AnalysisHandler.hh" #include "Rivet/AnalysisInfo.hh" #include "Rivet/Tools/BeamConstraint.hh" +#include "Rivet/Projections/ImpactParameterProjection.hh" +#include "Rivet/Projections/GeneratedPercentileProjection.hh" +#include "Rivet/Projections/UserCentEstimate.hh" +#include "Rivet/Projections/CentralityProjection.hh" namespace Rivet { Analysis::Analysis(const string& name) : _crossSection(-1.0), _gotCrossSection(false), _analysishandler(NULL) { ProjectionApplier::_allowProjReg = false; _defaultname = name; unique_ptr ai = AnalysisInfo::make(name); assert(ai); _info = move(ai); assert(_info); } double Analysis::sqrtS() const { return handler().sqrtS(); } const ParticlePair& Analysis::beams() const { return handler().beams(); } const PdgIdPair Analysis::beamIds() const { return handler().beamIds(); } const string Analysis::histoDir() const { /// @todo Cache in a member variable string _histoDir; if (_histoDir.empty()) { _histoDir = "/" + name(); if (handler().runName().length() > 0) { _histoDir = "/" + handler().runName() + _histoDir; } replace_all(_histoDir, "//", "/"); //< iterates until none } return _histoDir; } const string Analysis::histoPath(const string& hname) const { const string path = histoDir() + "/" + hname; return path; } const string Analysis::histoPath(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const { return histoDir() + "/" + mkAxisCode(datasetId, xAxisId, yAxisId); } const string Analysis::mkAxisCode(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const { stringstream axisCode; axisCode << "d"; if (datasetId < 10) axisCode << 0; axisCode << datasetId; axisCode << "-x"; if (xAxisId < 10) axisCode << 0; axisCode << xAxisId; axisCode << "-y"; if (yAxisId < 10) axisCode << 0; axisCode << yAxisId; return axisCode.str(); } Log& Analysis::getLog() const { string logname = "Rivet.Analysis." + name(); return Log::getLog(logname); } /////////////////////////////////////////// size_t Analysis::numEvents() const { return handler().numEvents(); } double Analysis::sumW() const { return handler().sumW(); } double Analysis::sumW2() const { return handler().sumW2(); } /////////////////////////////////////////// bool Analysis::isCompatible(const ParticlePair& beams) const { return isCompatible(beams.first.pid(), beams.second.pid(), beams.first.energy(), beams.second.energy()); } bool Analysis::isCompatible(PdgId beam1, PdgId beam2, double e1, double e2) const { PdgIdPair beams(beam1, beam2); pair energies(e1, e2); return isCompatible(beams, energies); } bool Analysis::isCompatible(const PdgIdPair& beams, const pair& energies) const { // First check the beam IDs bool beamIdsOk = false; for (const PdgIdPair& bp : requiredBeams()) { if (compatible(beams, bp)) { beamIdsOk = true; break; } } if (!beamIdsOk) return false; // Next check that the energies are compatible (within 1% or 1 GeV, whichever is larger, for a bit of UI forgiveness) /// @todo Use some sort of standard ordering to improve comparisons, esp. when the two beams are different particles bool beamEnergiesOk = requiredEnergies().size() > 0 ? false : true; typedef pair DoublePair; for (const DoublePair& ep : requiredEnergies()) { if ((fuzzyEquals(ep.first, energies.first, 0.01) && fuzzyEquals(ep.second, energies.second, 0.01)) || (fuzzyEquals(ep.first, energies.second, 0.01) && fuzzyEquals(ep.second, energies.first, 0.01)) || (abs(ep.first - energies.first) < 1*GeV && abs(ep.second - energies.second) < 1*GeV) || (abs(ep.first - energies.second) < 1*GeV && abs(ep.second - energies.first) < 1*GeV)) { beamEnergiesOk = true; break; } } return beamEnergiesOk; /// @todo Need to also check internal consistency of the analysis' /// beam requirements with those of the projections it uses. } /////////////////////////////////////////// Analysis& Analysis::setCrossSection(double xs) { _crossSection = xs; _gotCrossSection = true; return *this; } double Analysis::crossSection() const { if (!_gotCrossSection || std::isnan(_crossSection)) { string errMsg = "You did not set the cross section for the analysis " + name(); throw Error(errMsg); } return _crossSection; } double Analysis::crossSectionPerEvent() const { const double sumW = sumOfWeights(); assert(sumW != 0.0); return _crossSection / sumW; } //////////////////////////////////////////////////////////// // Histogramming void Analysis::_cacheRefData() const { if (_refdata.empty()) { MSG_TRACE("Getting refdata cache for paper " << name()); _refdata = getRefData(getRefDataName()); } } + vector Analysis::getAllData(bool includeorphans) const{ + return handler().getData(includeorphans); + } CounterPtr Analysis::bookCounter(const string& cname, const string& title) { // const string& xtitle, // const string& ytitle) { const string path = histoPath(cname); CounterPtr ctr = make_shared(path, title); addAnalysisObject(ctr); MSG_TRACE("Made counter " << cname << " for " << name()); // hist->setAnnotation("XLabel", xtitle); // hist->setAnnotation("YLabel", ytitle); return ctr; } CounterPtr Analysis::bookCounter(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId, const string& title) { // const string& xtitle, // const string& ytitle) { const string axisCode = mkAxisCode(datasetId, xAxisId, yAxisId); return bookCounter(axisCode, title); } Histo1DPtr Analysis::bookHisto1D(const string& hname, size_t nbins, double lower, double upper, const string& title, const string& xtitle, const string& ytitle) { Histo1DPtr hist; try { // try to bind to pre-existing // AnalysisObjectPtr ao = getAnalysisObject(path); // hist = dynamic_pointer_cast(ao); hist = getHisto1D(hname); /// @todo Test that cast worked /// @todo Also test that binning is as expected? MSG_TRACE("Bound pre-existing histogram " << hname << " for " << name()); } catch (...) { // binding failed; make it from scratch hist = make_shared(nbins, lower, upper, histoPath(hname), title); addAnalysisObject(hist); MSG_TRACE("Made histogram " << hname << " for " << name()); } hist->setTitle(title); hist->setAnnotation("XLabel", xtitle); hist->setAnnotation("YLabel", ytitle); return hist; } Histo1DPtr Analysis::bookHisto1D(const string& hname, const vector& binedges, const string& title, const string& xtitle, const string& ytitle) { Histo1DPtr hist; try { // try to bind to pre-existing // AnalysisObjectPtr ao = getAnalysisObject(path); // hist = dynamic_pointer_cast(ao); hist = getHisto1D(hname); /// @todo Test that cast worked /// @todo Also test that binning is as expected? MSG_TRACE("Bound pre-existing histogram " << hname << " for " << name()); } catch (...) { // binding failed; make it from scratch hist = make_shared(binedges, histoPath(hname), title); addAnalysisObject(hist); MSG_TRACE("Made histogram " << hname << " for " << name()); } hist->setTitle(title); hist->setAnnotation("XLabel", xtitle); hist->setAnnotation("YLabel", ytitle); return hist; } Histo1DPtr Analysis::bookHisto1D(const string& hname, const initializer_list& binedges, const string& title, const string& xtitle, const string& ytitle) { return bookHisto1D(hname, vector{binedges}, title, xtitle, ytitle); } Histo1DPtr Analysis::bookHisto1D(const string& hname, const Scatter2D& refscatter, const string& title, const string& xtitle, const string& ytitle) { Histo1DPtr hist; try { // try to bind to pre-existing // AnalysisObjectPtr ao = getAnalysisObject(path); // hist = dynamic_pointer_cast(ao); hist = getHisto1D(hname); /// @todo Test that cast worked /// @todo Also test that binning is as expected? MSG_TRACE("Bound pre-existing histogram " << hname << " for " << name()); } catch (...) { // binding failed; make it from scratch hist = make_shared(refscatter, histoPath(hname)); if (hist->hasAnnotation("IsRef")) hist->rmAnnotation("IsRef"); addAnalysisObject(hist); MSG_TRACE("Made histogram " << hname << " for " << name()); } hist->setTitle(title); hist->setAnnotation("XLabel", xtitle); hist->setAnnotation("YLabel", ytitle); return hist; } Histo1DPtr Analysis::bookHisto1D(const string& hname, const string& title, const string& xtitle, const string& ytitle) { const Scatter2D& refdata = refData(hname); return bookHisto1D(hname, refdata, title, xtitle, ytitle); } Histo1DPtr Analysis::bookHisto1D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId, const string& title, const string& xtitle, const string& ytitle) { const string axisCode = mkAxisCode(datasetId, xAxisId, yAxisId); return bookHisto1D(axisCode, title, xtitle, ytitle); } /// @todo Add booking methods which take a path, titles and *a reference Scatter from which to book* ///////////////// Histo2DPtr Analysis::bookHisto2D(const string& hname, size_t nxbins, double xlower, double xupper, size_t nybins, double ylower, double yupper, const string& title, const string& xtitle, const string& ytitle, const string& ztitle) { const string path = histoPath(hname); Histo2DPtr hist = make_shared(nxbins, xlower, xupper, nybins, ylower, yupper, path, title); addAnalysisObject(hist); MSG_TRACE("Made 2D histogram " << hname << " for " << name()); hist->setAnnotation("XLabel", xtitle); hist->setAnnotation("YLabel", ytitle); hist->setAnnotation("ZLabel", ztitle); return hist; } Histo2DPtr Analysis::bookHisto2D(const string& hname, const vector& xbinedges, const vector& ybinedges, const string& title, const string& xtitle, const string& ytitle, const string& ztitle) { const string path = histoPath(hname); Histo2DPtr hist = make_shared(xbinedges, ybinedges, path, title); addAnalysisObject(hist); MSG_TRACE("Made 2D histogram " << hname << " for " << name()); hist->setAnnotation("XLabel", xtitle); hist->setAnnotation("YLabel", ytitle); hist->setAnnotation("ZLabel", ztitle); return hist; } Histo2DPtr Analysis::bookHisto2D(const string& hname, const initializer_list& xbinedges, const initializer_list& ybinedges, const string& title, const string& xtitle, const string& ytitle, const string& ztitle) { return bookHisto2D(hname, vector{xbinedges}, vector{ybinedges}, title, xtitle, ytitle, ztitle); } Histo2DPtr Analysis::bookHisto2D(const string& hname, const Scatter3D& refscatter, const string& title, const string& xtitle, const string& ytitle, const string& ztitle) { const string path = histoPath(hname); Histo2DPtr hist( new Histo2D(refscatter, path) ); addAnalysisObject(hist); MSG_TRACE("Made 2D histogram " << hname << " for " << name()); if (hist->hasAnnotation("IsRef")) hist->rmAnnotation("IsRef"); hist->setTitle(title); hist->setAnnotation("XLabel", xtitle); hist->setAnnotation("YLabel", ytitle); hist->setAnnotation("ZLabel", ztitle); return hist; } Histo2DPtr Analysis::bookHisto2D(const string& hname, const string& title, const string& xtitle, const string& ytitle, const string& ztitle) { const Scatter3D& refdata = refData(hname); return bookHisto2D(hname, refdata, title, xtitle, ytitle, ztitle); } Histo2DPtr Analysis::bookHisto2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId, const string& title, const string& xtitle, const string& ytitle, const string& ztitle) { const string axisCode = mkAxisCode(datasetId, xAxisId, yAxisId); return bookHisto2D(axisCode, title, xtitle, ytitle, ztitle); } ///////////////// Profile1DPtr Analysis::bookProfile1D(const string& hname, size_t nbins, double lower, double upper, const string& title, const string& xtitle, const string& ytitle) { const string path = histoPath(hname); Profile1DPtr prof = make_shared(nbins, lower, upper, path, title); addAnalysisObject(prof); MSG_TRACE("Made profile histogram " << hname << " for " << name()); prof->setAnnotation("XLabel", xtitle); prof->setAnnotation("YLabel", ytitle); return prof; } Profile1DPtr Analysis::bookProfile1D(const string& hname, const vector& binedges, const string& title, const string& xtitle, const string& ytitle) { const string path = histoPath(hname); Profile1DPtr prof = make_shared(binedges, path, title); addAnalysisObject(prof); MSG_TRACE("Made profile histogram " << hname << " for " << name()); prof->setAnnotation("XLabel", xtitle); prof->setAnnotation("YLabel", ytitle); return prof; } Profile1DPtr Analysis::bookProfile1D(const string& hname, const initializer_list& binedges, const string& title, const string& xtitle, const string& ytitle) { return bookProfile1D(hname, vector{binedges}, title, xtitle, ytitle); } Profile1DPtr Analysis::bookProfile1D(const string& hname, const Scatter2D& refscatter, const string& title, const string& xtitle, const string& ytitle) { const string path = histoPath(hname); Profile1DPtr prof = make_shared(refscatter, path); addAnalysisObject(prof); MSG_TRACE("Made profile histogram " << hname << " for " << name()); if (prof->hasAnnotation("IsRef")) prof->rmAnnotation("IsRef"); prof->setTitle(title); prof->setAnnotation("XLabel", xtitle); prof->setAnnotation("YLabel", ytitle); return prof; } Profile1DPtr Analysis::bookProfile1D(const string& hname, const string& title, const string& xtitle, const string& ytitle) { const Scatter2D& refdata = refData(hname); return bookProfile1D(hname, refdata, title, xtitle, ytitle); } Profile1DPtr Analysis::bookProfile1D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId, const string& title, const string& xtitle, const string& ytitle) { const string axisCode = mkAxisCode(datasetId, xAxisId, yAxisId); return bookProfile1D(axisCode, title, xtitle, ytitle); } /////////////////// Profile2DPtr Analysis::bookProfile2D(const string& hname, size_t nxbins, double xlower, double xupper, size_t nybins, double ylower, double yupper, const string& title, const string& xtitle, const string& ytitle, const string& ztitle) { const string path = histoPath(hname); Profile2DPtr prof = make_shared(nxbins, xlower, xupper, nybins, ylower, yupper, path, title); addAnalysisObject(prof); MSG_TRACE("Made 2D profile histogram " << hname << " for " << name()); prof->setAnnotation("XLabel", xtitle); prof->setAnnotation("YLabel", ytitle); prof->setAnnotation("ZLabel", ztitle); return prof; } Profile2DPtr Analysis::bookProfile2D(const string& hname, const vector& xbinedges, const vector& ybinedges, const string& title, const string& xtitle, const string& ytitle, const string& ztitle) { const string path = histoPath(hname); Profile2DPtr prof = make_shared(xbinedges, ybinedges, path, title); addAnalysisObject(prof); MSG_TRACE("Made 2D profile histogram " << hname << " for " << name()); prof->setAnnotation("XLabel", xtitle); prof->setAnnotation("YLabel", ytitle); prof->setAnnotation("ZLabel", ztitle); return prof; } Profile2DPtr Analysis::bookProfile2D(const string& hname, const initializer_list& xbinedges, const initializer_list& ybinedges, const string& title, const string& xtitle, const string& ytitle, const string& ztitle) { return bookProfile2D(hname, vector{xbinedges}, vector{ybinedges}, title, xtitle, ytitle, ztitle); } Profile2DPtr Analysis::bookProfile2D(const string& hname, const Scatter3D& refscatter, const string& title, const string& xtitle, const string& ytitle, const string& ztitle) { const string path = histoPath(hname); Profile2DPtr prof( new Profile2D(refscatter, path) ); addAnalysisObject(prof); MSG_TRACE("Made 2D profile histogram " << hname << " for " << name()); if (prof->hasAnnotation("IsRef")) prof->rmAnnotation("IsRef"); prof->setTitle(title); prof->setAnnotation("XLabel", xtitle); prof->setAnnotation("YLabel", ytitle); prof->setAnnotation("ZLabel", ztitle); return prof; } Profile2DPtr Analysis::bookProfile2D(const string& hname, const string& title, const string& xtitle, const string& ytitle, const string& ztitle) { const Scatter3D& refdata = refData(hname); return bookProfile2D(hname, refdata, title, xtitle, ytitle, ztitle); } Profile2DPtr Analysis::bookProfile2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId, const string& title, const string& xtitle, const string& ytitle, const string& ztitle) { const string axisCode = mkAxisCode(datasetId, xAxisId, yAxisId); return bookProfile2D(axisCode, title, xtitle, ytitle, ztitle); } ///////////////// Scatter2DPtr Analysis::bookScatter2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId, bool copy_pts, const string& title, const string& xtitle, const string& ytitle) { const string axisCode = mkAxisCode(datasetId, xAxisId, yAxisId); return bookScatter2D(axisCode, copy_pts, title, xtitle, ytitle); } Scatter2DPtr Analysis::bookScatter2D(const string& hname, bool copy_pts, const string& title, const string& xtitle, const string& ytitle) { Scatter2DPtr s; const string path = histoPath(hname); if (copy_pts) { const Scatter2D& refdata = refData(hname); s = make_shared(refdata, path); for (Point2D& p : s->points()) p.setY(0, 0); } else { s = make_shared(path); } addAnalysisObject(s); MSG_TRACE("Made scatter " << hname << " for " << name()); if (s->hasAnnotation("IsRef")) s->rmAnnotation("IsRef"); s->setTitle(title); s->setAnnotation("XLabel", xtitle); s->setAnnotation("YLabel", ytitle); return s; } Scatter2DPtr Analysis::bookScatter2D(const string& hname, size_t npts, double lower, double upper, const string& title, const string& xtitle, const string& ytitle) { Scatter2DPtr s; const string path = histoPath(hname); try { // try to bind to pre-existing s = getAnalysisObject(hname); /// @todo Also test that binning is as expected? MSG_TRACE("Bound pre-existing scatter " << path << " for " << name()); } catch (...) { // binding failed; make it from scratch s = make_shared(path); const double binwidth = (upper-lower)/npts; for (size_t pt = 0; pt < npts; ++pt) { const double bincentre = lower + (pt + 0.5) * binwidth; s->addPoint(bincentre, 0, binwidth/2.0, 0); } addAnalysisObject(s); MSG_TRACE("Made scatter " << hname << " for " << name()); } s->setTitle(title); s->setAnnotation("XLabel", xtitle); s->setAnnotation("YLabel", ytitle); return s; } Scatter2DPtr Analysis::bookScatter2D(const string& hname, const vector& binedges, const string& title, const string& xtitle, const string& ytitle) { const string path = histoPath(hname); Scatter2DPtr s = make_shared(path); for (size_t pt = 0; pt < binedges.size()-1; ++pt) { const double bincentre = (binedges[pt] + binedges[pt+1]) / 2.0; const double binwidth = binedges[pt+1] - binedges[pt]; s->addPoint(bincentre, 0, binwidth/2.0, 0); } addAnalysisObject(s); MSG_TRACE("Made scatter " << hname << " for " << name()); s->setTitle(title); s->setAnnotation("XLabel", xtitle); s->setAnnotation("YLabel", ytitle); return s; } ///////////////////// void Analysis::divide(CounterPtr c1, CounterPtr c2, Scatter1DPtr s) const { const string path = s->path(); *s = *c1 / *c2; s->setPath(path); } void Analysis::divide(const Counter& c1, const Counter& c2, Scatter1DPtr s) const { const string path = s->path(); *s = c1 / c2; s->setPath(path); } void Analysis::divide(Histo1DPtr h1, Histo1DPtr h2, Scatter2DPtr s) const { const string path = s->path(); *s = *h1 / *h2; s->setPath(path); } void Analysis::divide(const Histo1D& h1, const Histo1D& h2, Scatter2DPtr s) const { const string path = s->path(); *s = h1 / h2; s->setPath(path); } void Analysis::divide(Profile1DPtr p1, Profile1DPtr p2, Scatter2DPtr s) const { const string path = s->path(); *s = *p1 / *p2; s->setPath(path); } void Analysis::divide(const Profile1D& p1, const Profile1D& p2, Scatter2DPtr s) const { const string path = s->path(); *s = p1 / p2; s->setPath(path); } void Analysis::divide(Histo2DPtr h1, Histo2DPtr h2, Scatter3DPtr s) const { const string path = s->path(); *s = *h1 / *h2; s->setPath(path); } void Analysis::divide(const Histo2D& h1, const Histo2D& h2, Scatter3DPtr s) const { const string path = s->path(); *s = h1 / h2; s->setPath(path); } void Analysis::divide(Profile2DPtr p1, Profile2DPtr p2, Scatter3DPtr s) const { const string path = s->path(); *s = *p1 / *p2; s->setPath(path); } void Analysis::divide(const Profile2D& p1, const Profile2D& p2, Scatter3DPtr s) const { const string path = s->path(); *s = p1 / p2; s->setPath(path); } /// @todo Counter and Histo2D efficiencies and asymms void Analysis::efficiency(Histo1DPtr h1, Histo1DPtr h2, Scatter2DPtr s) const { const string path = s->path(); *s = YODA::efficiency(*h1, *h2); s->setPath(path); } void Analysis::efficiency(const Histo1D& h1, const Histo1D& h2, Scatter2DPtr s) const { const string path = s->path(); *s = YODA::efficiency(h1, h2); s->setPath(path); } void Analysis::asymm(Histo1DPtr h1, Histo1DPtr h2, Scatter2DPtr s) const { const string path = s->path(); *s = YODA::asymm(*h1, *h2); s->setPath(path); } void Analysis::asymm(const Histo1D& h1, const Histo1D& h2, Scatter2DPtr s) const { const string path = s->path(); *s = YODA::asymm(h1, h2); s->setPath(path); } void Analysis::scale(CounterPtr cnt, double factor) { if (!cnt) { MSG_WARNING("Failed to scale counter=NULL in analysis " << name() << " (scale=" << factor << ")"); return; } if (std::isnan(factor) || std::isinf(factor)) { MSG_WARNING("Failed to scale counter=" << cnt->path() << " in analysis: " << name() << " (invalid scale factor = " << factor << ")"); factor = 0; } MSG_TRACE("Scaling counter " << cnt->path() << " by factor " << factor); try { cnt->scaleW(factor); } catch (YODA::Exception& we) { MSG_WARNING("Could not scale counter " << cnt->path()); return; } } void Analysis::normalize(Histo1DPtr histo, double norm, bool includeoverflows) { if (!histo) { MSG_WARNING("Failed to normalize histo=NULL in analysis " << name() << " (norm=" << norm << ")"); return; } MSG_TRACE("Normalizing histo " << histo->path() << " to " << norm); try { const double hint = histo->integral(includeoverflows); if (hint == 0) MSG_WARNING("Skipping histo with null area " << histo->path()); else histo->normalize(norm, includeoverflows); } catch (YODA::Exception& we) { MSG_WARNING("Could not normalize histo " << histo->path()); return; } } void Analysis::scale(Histo1DPtr histo, double factor) { if (!histo) { MSG_WARNING("Failed to scale histo=NULL in analysis " << name() << " (scale=" << factor << ")"); return; } if (std::isnan(factor) || std::isinf(factor)) { MSG_WARNING("Failed to scale histo=" << histo->path() << " in analysis: " << name() << " (invalid scale factor = " << factor << ")"); factor = 0; } MSG_TRACE("Scaling histo " << histo->path() << " by factor " << factor); try { histo->scaleW(factor); } catch (YODA::Exception& we) { MSG_WARNING("Could not scale histo " << histo->path()); return; } } void Analysis::normalize(Histo2DPtr histo, double norm, bool includeoverflows) { if (!histo) { MSG_ERROR("Failed to normalize histo=NULL in analysis " << name() << " (norm=" << norm << ")"); return; } MSG_TRACE("Normalizing histo " << histo->path() << " to " << norm); try { const double hint = histo->integral(includeoverflows); if (hint == 0) MSG_WARNING("Skipping histo with null area " << histo->path()); else histo->normalize(norm, includeoverflows); } catch (YODA::Exception& we) { MSG_WARNING("Could not normalize histo " << histo->path()); return; } } void Analysis::scale(Histo2DPtr histo, double factor) { if (!histo) { MSG_ERROR("Failed to scale histo=NULL in analysis " << name() << " (scale=" << factor << ")"); return; } if (std::isnan(factor) || std::isinf(factor)) { MSG_ERROR("Failed to scale histo=" << histo->path() << " in analysis: " << name() << " (invalid scale factor = " << factor << ")"); factor = 0; } MSG_TRACE("Scaling histo " << histo->path() << " by factor " << factor); try { histo->scaleW(factor); } catch (YODA::Exception& we) { MSG_WARNING("Could not scale histo " << histo->path()); return; } } void Analysis::integrate(Histo1DPtr h, Scatter2DPtr s) const { // preserve the path info const string path = s->path(); *s = toIntegralHisto(*h); s->setPath(path); } void Analysis::integrate(const Histo1D& h, Scatter2DPtr s) const { // preserve the path info const string path = s->path(); *s = toIntegralHisto(h); s->setPath(path); } /// @todo 2D versions of integrate... defined how, exactly?!? ////////////////////////////////// void Analysis::addAnalysisObject(AnalysisObjectPtr ao) { _analysisobjects.push_back(ao); } void Analysis::removeAnalysisObject(const string& path) { for (vector::iterator it = _analysisobjects.begin(); it != _analysisobjects.end(); ++it) { if ((*it)->path() == path) { _analysisobjects.erase(it); break; } } } void Analysis::removeAnalysisObject(AnalysisObjectPtr ao) { for (vector::iterator it = _analysisobjects.begin(); it != _analysisobjects.end(); ++it) { if (*it == ao) { _analysisobjects.erase(it); break; } } } +const CentralityProjection & +Analysis::declareCentrality(const SingleValueProjection &proj, + string calAnaName, string calHistName, + const string projName, bool increasing) { + + CentralityProjection cproj; + + // Select the centrality variable from option. Use REF as default. + // Other selections are "GEN", "IMP" and "USR" (USR only in HEPMC 3). + string sel = getOption("cent","REF"); + set done; + + if ( sel == "REF" ) { + Scatter2DPtr refscat; + auto refmap = getRefData(calAnaName); + if ( refmap.find(calHistName) != refmap.end() ) + refscat = + dynamic_pointer_cast(refmap.find(calHistName)->second); + + if ( !refscat ) { + MSG_WARNING("No reference calibration histogram for " << + "CentralityProjection " << projName << " found " << + "(requested histogram " << calHistName << " in " << + calAnaName << ")"); + } + else { + MSG_INFO("Found calibration histogram " << sel << " " << refscat->path()); + cproj.add(PercentileProjection(proj, refscat, increasing), sel); + } + } + else if ( sel == "GEN" ) { + Histo1DPtr genhist; + string histpath = "/" + calAnaName + "/" + calHistName; + for ( AnalysisObjectPtr ao : handler().getData(true) ) { + if ( ao->path() == histpath ) + genhist = dynamic_pointer_cast(ao); + } + + if ( !genhist || genhist->numEntries() <= 1 ) { + MSG_WARNING("No generated calibration histogram for " << + "CentralityProjection " << projName << " found " << + "(requested histogram " << calHistName << " in " << + calAnaName << ")"); + } + else { + MSG_INFO("Found calibration histogram " << sel << " " << genhist->path()); + cproj.add(PercentileProjection(proj, genhist, increasing), sel); + } + } + else if ( sel == "IMP" ) { + Histo1DPtr imphist = + getAnalysisObject(calAnaName, calHistName + "_IMP"); + if ( !imphist || imphist->numEntries() <= 1 ) { + MSG_WARNING("No impact parameter calibration histogram for " << + "CentralityProjection " << projName << " found " << + "(requested histogram " << calHistName << "_IMP in " << + calAnaName << ")"); + } + else { + MSG_INFO("Found calibration histogram " << sel << " " << imphist->path()); + cproj.add(PercentileProjection(ImpactParameterProjection(), + imphist, true), sel); + } + } + else if ( sel == "USR" ) { +#if HEPMC_VERSION_CODE >= 3000000 + Histo1DPtr usrhist = + getAnalysisObject(calAnaName, calHistName + "_USR"); + if ( !usrhist || usrhist->numEntries() <= 1 ) { + MSG_WARNING("No user-defined calibration histogram for " << + "CentralityProjection " << projName << " found " << + "(requested histogram " << calHistName << "_USR in " << + calAnaName << ")"); + continue; + } + else { + MSG_INFO("Found calibration histogram " << sel << " " << usrhist->path()); + cproj.add((UserCentEstimate(), usrhist, true), sel); + } +#else + MSG_WARNING("UserCentEstimate is only available with HepMC3."); +#endif + } + else if ( sel == "RAW" ) { +#if HEPMC_VERSION_CODE >= 3000000 + cproj.add(GeneratedCentrality(), sel); +#else + MSG_WARNING("GeneratedCentrality is only available with HepMC3."); +#endif + } + else + MSG_WARNING("'" << sel << "' is not a valid PercentileProjection tag."); + + if ( cproj.empty() ) + MSG_WARNING("CentralityProjection " << projName + << " did not contain any valid PercentileProjections."); + + return declare(cproj, projName); + +} } diff --git a/src/Core/AnalysisHandler.cc b/src/Core/AnalysisHandler.cc --- a/src/Core/AnalysisHandler.cc +++ b/src/Core/AnalysisHandler.cc @@ -1,379 +1,381 @@ // -*- C++ -*- #include "Rivet/Config/RivetCommon.hh" #include "Rivet/AnalysisHandler.hh" #include "Rivet/Analysis.hh" #include "Rivet/Tools/ParticleName.hh" #include "Rivet/Tools/BeamConstraint.hh" #include "Rivet/Tools/Logging.hh" #include "Rivet/Projections/Beam.hh" #include "YODA/IO.h" namespace Rivet { AnalysisHandler::AnalysisHandler(const string& runname) : _runname(runname), _eventcounter("/_EVTCOUNT"), _xs(NAN), _xserr(NAN), _initialised(false), _ignoreBeams(false) { } AnalysisHandler::~AnalysisHandler() { } Log& AnalysisHandler::getLog() const { return Log::getLog("Rivet.Analysis.Handler"); } void AnalysisHandler::init(const GenEvent& ge) { if (_initialised) throw UserError("AnalysisHandler::init has already been called: cannot re-initialize!"); setRunBeams(Rivet::beams(ge)); MSG_DEBUG("Initialising the analysis handler"); _eventcounter.reset(); // Check that analyses are beam-compatible, and remove those that aren't const size_t num_anas_requested = analysisNames().size(); vector anamestodelete; for (const AnaHandle a : _analyses) { if (!_ignoreBeams && !a->isCompatible(beams())) { //MSG_DEBUG(a->name() << " requires beams " << a->requiredBeams() << " @ " << a->requiredEnergies() << " GeV"); anamestodelete.push_back(a->name()); } } for (const string& aname : anamestodelete) { MSG_WARNING("Analysis '" << aname << "' is incompatible with the provided beams: removing"); removeAnalysis(aname); } if (num_anas_requested > 0 && analysisNames().empty()) { cerr << "All analyses were incompatible with the first event's beams\n" << "Exiting, since this probably wasn't intentional!" << endl; exit(1); } // Warn if any analysis' status is not unblemished for (const AnaHandle a : analyses()) { if (toUpper(a->status()) == "PRELIMINARY") { MSG_WARNING("Analysis '" << a->name() << "' is preliminary: be careful, it may change and/or be renamed!"); } else if (toUpper(a->status()) == "OBSOLETE") { MSG_WARNING("Analysis '" << a->name() << "' is obsolete: please update!"); } else if (toUpper(a->status()).find("UNVALIDATED") != string::npos) { MSG_WARNING("Analysis '" << a->name() << "' is unvalidated: be careful, it may be broken!"); } } // Initialize the remaining analyses for (AnaHandle a : _analyses) { MSG_DEBUG("Initialising analysis: " << a->name()); try { // Allow projection registration in the init phase onwards a->_allowProjReg = true; a->init(); //MSG_DEBUG("Checking consistency of analysis: " << a->name()); //a->checkConsistency(); } catch (const Error& err) { cerr << "Error in " << a->name() << "::init method: " << err.what() << endl; exit(1); } MSG_DEBUG("Done initialising analysis: " << a->name()); } _initialised = true; MSG_DEBUG("Analysis handler initialised"); } void AnalysisHandler::analyze(const GenEvent& ge) { // Call init with event as template if not already initialised if (!_initialised) init(ge); assert(_initialised); // Ensure that beam details match those from the first event (if we're checking beams) if ( !_ignoreBeams ) { const PdgIdPair beams = Rivet::beamIds(ge); const double sqrts = Rivet::sqrtS(ge); if (!compatible(beams, _beams) || !fuzzyEquals(sqrts, sqrtS())) { cerr << "Event beams mismatch: " << PID::toBeamsString(beams) << " @ " << sqrts/GeV << " GeV" << " vs. first beams " << this->beams() << " @ " << this->sqrtS()/GeV << " GeV" << endl; exit(1); } } // Create the Rivet event wrapper /// @todo Filter/normalize the event here Event event(ge); // Weights /// @todo Drop this / just report first weight when we support multiweight events _eventcounter.fill(event.weight()); MSG_DEBUG("Event #" << _eventcounter.numEntries() << " weight = " << event.weight()); // Cross-section #ifdef HEPMC_HAS_CROSS_SECTION if (ge.cross_section()) { _xs = ge.cross_section()->cross_section(); _xserr = ge.cross_section()->cross_section_error(); } #endif // Run the analyses for (AnaHandle a : _analyses) { MSG_TRACE("About to run analysis " << a->name()); try { a->analyze(event); } catch (const Error& err) { cerr << "Error in " << a->name() << "::analyze method: " << err.what() << endl; exit(1); } MSG_TRACE("Finished running analysis " << a->name()); } } void AnalysisHandler::analyze(const GenEvent* ge) { if (ge == nullptr) { MSG_ERROR("AnalysisHandler received null pointer to GenEvent"); //throw Error("AnalysisHandler received null pointer to GenEvent"); } analyze(*ge); } void AnalysisHandler::finalize() { if (!_initialised) return; MSG_INFO("Finalising analyses"); for (AnaHandle a : _analyses) { a->setCrossSection(_xs); try { a->finalize(); } catch (const Error& err) { cerr << "Error in " << a->name() << "::finalize method: " << err.what() << endl; exit(1); } } // Print out number of events processed const int nevts = _eventcounter.numEntries(); MSG_INFO("Processed " << nevts << " event" << (nevts != 1 ? "s" : "")); // // Delete analyses // MSG_DEBUG("Deleting analyses"); // _analyses.clear(); // Print out MCnet boilerplate cout << endl; cout << "The MCnet usage guidelines apply to Rivet: see http://www.montecarlonet.org/GUIDELINES" << endl; cout << "Please acknowledge plots made with Rivet analyses, and cite arXiv:1003.0694 (http://arxiv.org/abs/1003.0694)" << endl; } AnalysisHandler& AnalysisHandler::addAnalysis(const string& analysisname) { // Check for a duplicate analysis /// @todo Might we want to be able to run an analysis twice, with different params? /// Requires avoiding histo tree clashes, i.e. storing the histos on the analysis objects. string ananame = analysisname; vector anaopt = split(analysisname, ":"); if ( anaopt.size() > 1 ) ananame = anaopt[0]; AnaHandle analysis( AnalysisLoader::getAnalysis(ananame) ); if (analysis.get() != 0) { // < Check for null analysis. MSG_DEBUG("Adding analysis '" << analysisname << "'"); map opts; for ( int i = 1, N = anaopt.size(); i < N; ++i ) { vector opt = split(anaopt[i], "="); if ( opt.size() != 2 ) { MSG_WARNING("Error in option specification. Skipping analysis " << analysisname); return *this; } if ( !analysis->info().validOption(opt[0], opt[1]) ) { MSG_WARNING("Cannot set option '" << opt[0] << "' to '" << opt[1] << "'. Skipping analysis " << analysisname); return *this; } opts[opt[0]] = opt[1]; } for ( auto opt: opts) { analysis->_options[opt.first] = opt.second; analysis->_optstring += ":" + opt.first + "=" + opt.second; } for (const AnaHandle& a : _analyses) { if (a->name() == analysis->name() ) { MSG_WARNING("Analysis '" << analysisname << "' already registered: skipping duplicate"); return *this; } } analysis->_analysishandler = this; _analyses.insert(analysis); } else { MSG_WARNING("Analysis '" << analysisname << "' not found."); } // MSG_WARNING(_analyses.size()); // for (const AnaHandle& a : _analyses) MSG_WARNING(a->name()); return *this; } AnalysisHandler& AnalysisHandler::removeAnalysis(const string& analysisname) { std::shared_ptr toremove; for (const AnaHandle a : _analyses) { if (a->name() == analysisname) { toremove = a; break; } } if (toremove.get() != 0) { MSG_DEBUG("Removing analysis '" << analysisname << "'"); _analyses.erase(toremove); } return *this; } ///////////////////////////// void AnalysisHandler::addData(const std::vector& aos) { for (const AnalysisObjectPtr ao : aos) { const string path = ao->path(); if (path.size() > 1) { // path > "/" try { const string ananame = split(path, "/")[0]; AnaHandle a = analysis(ananame); a->addAnalysisObject(ao); /// @todo Need to statistically merge... } catch (const Error& e) { MSG_WARNING(e.what()); } } } } void AnalysisHandler::readData(const string& filename) { vector aos; try { /// @todo Use new YODA SFINAE to fill the smart ptr vector directly vector aos_raw; YODA::read(filename, aos_raw); for (AnalysisObject* aor : aos_raw) aos.push_back(AnalysisObjectPtr(aor)); } catch (...) { //< YODA::ReadError& throw UserError("Unexpected error in reading file: " + filename); } if (!aos.empty()) addData(aos); } - vector AnalysisHandler::getData() const { + vector AnalysisHandler::getData(bool includeorphans) const { vector rtn; // Event counter rtn.push_back( make_shared(_eventcounter) ); // Cross-section + err as scatter YODA::Scatter1D::Points pts; pts.insert(YODA::Point1D(_xs, _xserr)); rtn.push_back( make_shared(pts, "/_XSEC") ); // Analysis histograms for (const AnaHandle a : analyses()) { vector aos = a->analysisObjects(); // MSG_WARNING(a->name() << " " << aos.size()); for (const AnalysisObjectPtr ao : aos) { // Exclude paths from final write-out if they contain a "TMP" layer (i.e. matching "/TMP/") /// @todo This needs to be much more nuanced for re-entrant histogramming if (ao->path().find("/TMP/") != string::npos) continue; rtn.push_back(ao); } } // Sort histograms alphanumerically by path before write-out sort(rtn.begin(), rtn.end(), [](AnalysisObjectPtr a, AnalysisObjectPtr b) {return a->path() < b->path();}); + if ( includeorphans ) + rtn.insert(rtn.end(), _orphanedPreloads.begin(), _orphanedPreloads.end()); return rtn; } void AnalysisHandler::writeData(const string& filename) const { const vector aos = getData(); try { YODA::write(filename, aos.begin(), aos.end()); } catch (...) { //< YODA::WriteError& throw UserError("Unexpected error in writing file: " + filename); } } std::vector AnalysisHandler::analysisNames() const { std::vector rtn; for (AnaHandle a : _analyses) { rtn.push_back(a->name()); } return rtn; } const AnaHandle AnalysisHandler::analysis(const std::string& analysisname) const { for (const AnaHandle a : analyses()) if (a->name() == analysisname) return a; throw Error("No analysis named '" + analysisname + "' registered in AnalysisHandler"); } AnalysisHandler& AnalysisHandler::addAnalyses(const std::vector& analysisnames) { for (const string& aname : analysisnames) { //MSG_DEBUG("Adding analysis '" << aname << "'"); addAnalysis(aname); } return *this; } AnalysisHandler& AnalysisHandler::removeAnalyses(const std::vector& analysisnames) { for (const string& aname : analysisnames) { removeAnalysis(aname); } return *this; } bool AnalysisHandler::needCrossSection() const { bool rtn = false; for (const AnaHandle a : _analyses) { if (!rtn) rtn = a->needsCrossSection(); if (rtn) break; } return rtn; } AnalysisHandler& AnalysisHandler::setCrossSection(double xs) { _xs = xs; return *this; } bool AnalysisHandler::hasCrossSection() const { return (!std::isnan(crossSection())); } AnalysisHandler& AnalysisHandler::addAnalysis(Analysis* analysis) { analysis->_analysishandler = this; _analyses.insert(AnaHandle(analysis)); return *this; } PdgIdPair AnalysisHandler::beamIds() const { return Rivet::beamIds(beams()); } double AnalysisHandler::sqrtS() const { return Rivet::sqrtS(beams()); } void AnalysisHandler::setIgnoreBeams(bool ignore) { _ignoreBeams=ignore; } } diff --git a/src/Projections/Makefile.am b/src/Projections/Makefile.am --- a/src/Projections/Makefile.am +++ b/src/Projections/Makefile.am @@ -1,45 +1,46 @@ noinst_LTLIBRARIES = libRivetProjections.la libRivetProjections_la_SOURCES = \ Beam.cc \ BeamThrust.cc \ ChargedFinalState.cc \ ChargedLeptons.cc \ CentralEtHCM.cc \ DISFinalState.cc \ DISKinematics.cc \ DISLepton.cc \ DressedLeptons.cc \ FastJets.cc \ FinalPartons.cc \ FinalState.cc \ FParameter.cc \ HadronicFinalState.cc \ HeavyHadrons.cc \ Hemispheres.cc \ IdentifiedFinalState.cc \ InitialQuarks.cc \ InvMassFinalState.cc \ JetAlg.cc \ JetShape.cc \ LeadingParticlesFinalState.cc \ MergedFinalState.cc \ MissingMomentum.cc \ NeutralFinalState.cc \ NonHadronicFinalState.cc \ NonPromptFinalState.cc \ ParisiTensor.cc \ ParticleFinder.cc \ PrimaryHadrons.cc \ + PrimaryParticles.cc \ PromptFinalState.cc \ Sphericity.cc \ Spherocity.cc \ TauFinder.cc \ Thrust.cc \ TriggerCDFRun0Run1.cc \ TriggerCDFRun2.cc \ TriggerUA5.cc \ UnstableFinalState.cc \ VetoedFinalState.cc \ VisibleFinalState.cc \ WFinder.cc \ ZFinder.cc diff --git a/src/Projections/PrimaryParticles.cc b/src/Projections/PrimaryParticles.cc new file mode 100644 --- /dev/null +++ b/src/Projections/PrimaryParticles.cc @@ -0,0 +1,91 @@ +/** + * @file PrimaryParticles.cc + * @author Christian Holm Christensen + * @date Mon Aug 27 09:06:03 2018 + * + * @brief Primary particle definition based on PDG IDs. + */ + +#include +#include +#include +#include +#include +#include +#include + +namespace Rivet { + + void PrimaryParticles::project(const Event& e) + { + _theParticles.clear(); // Clear cache + bool open = _cuts == Cuts::open(); + for (auto p : Rivet::particles(e.genEvent())) { + if (isPrimary(p) && (open || _cuts->accept(Particle(p)))) + _theParticles.push_back(Particle(*p)); + } + } + + bool PrimaryParticles::isPrimary(const HepMC::GenParticle* p) const + { + if(isIgnored(p)) return false; + if(!isPrimaryPID(p)) return false; + + // Loop back over ancestors that are _not_ ignored + const HepMC::GenParticle* m = p; + while ((m = ancestor(m,true))) { + if (isBeam(m)) return true; + if (isPrimaryPID(m)) return false; + if (!hasDecayed(m)) return false; + } + return true; + } + bool PrimaryParticles::isIgnored(const HepMC::GenParticle* p) const + { + return (p->status()==0 || (p->status()>=11 && p->status()<=200)); + } + + bool PrimaryParticles::isPrimaryPID(const HepMC::GenParticle* p) const + { + int thisPID = PID::abspid(p->pdg_id()); + for(const auto pid : _pdgIds) + if (thisPID == pid) return true; + return false; + } + + bool PrimaryParticles::isBeam(const HepMC::GenParticle* p) const + { + // Pythia6 uses 3 for initial state + return p && (p->status() == 3 || p->status() == 4); + } + + bool PrimaryParticles::hasDecayed(const HepMC::GenParticle* p) const + { + return p && p->status() == 2; + } + + const HepMC::GenParticle* + PrimaryParticles::ancestor(const HepMC::GenParticle* p) const + { + const HepMC::GenVertex* vtx = p->production_vertex(); + if (!vtx) return 0; + + HepMC::GenVertex::particles_in_const_iterator i = + vtx->particles_in_const_begin(); + if (i == vtx->particles_in_const_end()) return 0; + return *i; + } + + const HepMC::GenParticle* + PrimaryParticles::ancestor(const HepMC::GenParticle* p, bool) const + { + const HepMC::GenParticle* m = p; + do { + m = ancestor(m); + } while (m && isIgnored(m)); + return m; + } +} +// +// EOF +//