|
|
[Rivet-svn] r2186 - in trunk: bin data/plotinfo include/Rivet include/Rivet/Projections src/Analyses src/Core src/Projections
blackhole at projects.hepforge.org
blackhole at projects.hepforge.org
Sun Dec 13 20:39:44 GMT 2009
Author: buckley
Date: Sun Dec 13 20:39:43 2009
New Revision: 2186
Log:
Various API improvements, provoked by trying to get the D0 2008 W asymmetry analysis in a working state. I can't get any asymmetry with Pythia, though: this needs to be investigated before any release. Aargh.
Modified:
trunk/bin/compare-histos
trunk/data/plotinfo/D0_2008_S7837160.plot
trunk/include/Rivet/Analysis.hh
trunk/include/Rivet/Projections/InvMassFinalState.hh
trunk/include/Rivet/Projections/WFinder.hh
trunk/src/Analyses/CDF_2006_S6653332.cc
trunk/src/Analyses/CDF_2008_S7541902.cc
trunk/src/Analyses/CDF_2008_S8095620.cc
trunk/src/Analyses/D0_2008_S7837160.cc
trunk/src/Analyses/MC_WANALYSIS.cc
trunk/src/Core/Analysis.cc
trunk/src/Projections/InvMassFinalState.cc
trunk/src/Projections/WFinder.cc
Modified: trunk/bin/compare-histos
==============================================================================
--- trunk/bin/compare-histos Sun Dec 13 16:39:36 2009 (r2185)
+++ trunk/bin/compare-histos Sun Dec 13 20:39:43 2009 (r2186)
@@ -355,4 +355,4 @@
f.write(headstr + "\n" + "\n".join(histstrs))
f.close()
num_written += 1
- logging.info("Wrote %d histo files to %s" % (num_written, outdir))
+ logging.info("Wrote %d histo files" % num_written)
Modified: trunk/data/plotinfo/D0_2008_S7837160.plot
==============================================================================
--- trunk/data/plotinfo/D0_2008_S7837160.plot Sun Dec 13 16:39:36 2009 (r2185)
+++ trunk/data/plotinfo/D0_2008_S7837160.plot Sun Dec 13 20:39:43 2009 (r2186)
@@ -1,23 +1,20 @@
# BEGIN PLOT /D0_2008_S7837160/d01-x01-y01
-Title=W+ charge asymmetry for $25 > E_\perp > 35$ GeV
-XLabel=$|\eta|$ of leading electron
-#YLabel=A = $(\frac{\mathrm{d}{\sigma^+}}{\mathrm{d}{|\eta|}} - \frac{\mathrm{d}{\sigma^-}}{\mathrm{d}{|\eta|}}) / (\frac{\mathrm{d}{\sigma^+}}{\mathrm{d}{|\eta|}} + \frac{\mathrm{d}{\sigma^-}}{\mathrm{d}{|\eta|}})$
+Title=W charge asymmetry for $25 > E_\perp > 35$ GeV
+XLabel=$|\eta|$ of charged lepton
YLabel=$\mathcal{A}(\mathrm{d}{\sigma^+}\!/\!\mathrm{d}{|\eta|},\, \mathrm{d}{\sigma^-}\!/\!\mathrm{d}{|\eta|})$
LogY=0
# END PLOT
-# BEGIN PLOT /D0_2008_S7837160/d02-x01-y01
-Title=W- charge asymmetry for $25 > E_\perp > 35$ GeV
-XLabel=$|\eta|$ of leading positron
-#YLabel=A = $(\frac{\mathrm{d}{\sigma^+}}{\mathrm{d}{|\eta|}} - \frac{\mathrm{d}{\sigma^-}}{\mathrm{d}{|\eta|}}) / (\frac{\mathrm{d}{\sigma^+}}{\mathrm{d}{|\eta|}} + \frac{\mathrm{d}{\sigma^-}}{\mathrm{d}{|\eta|}})$
+# BEGIN PLOT /D0_2008_S7837160/d01-x01-y02
+Title=W charge asymmetry for $E_\perp > 35$ GeV
+XLabel=$|\eta|$ of charged lepton
YLabel=$\mathcal{A}(\mathrm{d}{\sigma^+}\!/\!\mathrm{d}{|\eta|},\, \mathrm{d}{\sigma^-}\!/\!\mathrm{d}{|\eta|})$
LogY=0
# END PLOT
-# BEGIN PLOT /D0_2008_S7837160/d03-x01-y01
-Title=W$\pm$ charge asymmetry for $E_\perp > 35$ GeV
-XLabel=$|\eta|$ of leading electron
-#YLabel=A = $(\frac{\mathrm{d}{\sigma^+}}{\mathrm{d}{|\eta|}} - \frac{\mathrm{d}{\sigma^-}}{\mathrm{d}{|\eta|}}) / (\frac{\mathrm{d}{\sigma^+}}{\mathrm{d}{|\eta|}} + \frac{\mathrm{d}{\sigma^-}}{\mathrm{d}{|\eta|}})$
+# BEGIN PLOT /D0_2008_S7837160/d01-x01-y03
+Title=W charge asymmetry for $E_\perp > 25$ GeV
+XLabel=$|\eta|$ of charged lepton
YLabel=$\mathcal{A}(\mathrm{d}{\sigma^+}\!/\!\mathrm{d}{|\eta|},\, \mathrm{d}{\sigma^-}\!/\!\mathrm{d}{|\eta|})$
LogY=0
# END PLOT
Modified: trunk/include/Rivet/Analysis.hh
==============================================================================
--- trunk/include/Rivet/Analysis.hh Sun Dec 13 16:39:36 2009 (r2185)
+++ trunk/include/Rivet/Analysis.hh Sun Dec 13 20:39:43 2009 (r2186)
@@ -196,17 +196,17 @@
/// @param histo The histogram to be normalised.
/// @param norm The new area of the histogram.
/// @warning The old histogram will be deleted, and its pointer set to zero.
- void normalize(AIDA::IHistogram1D*& histo, const double norm=1.0);
+ void normalize(AIDA::IHistogram1D*& histo, double norm=1.0);
/// Multiplicatively scale the given histogram, @a histo. After this call the
/// histogram will have been transformed to a DataPointSet with the same name and path.
/// @param histo The histogram to be scaled.
/// @param scale The factor used to multiply the histogram bin heights.
/// @warning The old histogram will be deleted, and its pointer set to zero.
- void scale(AIDA::IHistogram1D*& histo, const double scale);
+ void scale(AIDA::IHistogram1D*& histo, double scale);
/// Set the cross section from the generator
- Analysis& setCrossSection(const double& xs);
+ Analysis& setCrossSection(double xs);
/// Return true if this analysis needs to know the process cross-section.
bool needsCrossSection() const;
@@ -262,12 +262,18 @@
/// @name Internal histogram booking (for use by Analysis sub-classes).
//@{
+ /// Get bin edges for a named histo (using ref AIDA caching)
+ const BinEdges& binEdges(const std::string& hname) const;
+
+ /// Get bin edges for a numbered histo (using ref AIDA caching)
+ const BinEdges& binEdges(size_t datasetId, size_t xAxisId, size_t yAxisId) const;
+
/// Book a 1D histogram with @a nbins uniformly distributed across the range @a lower - @a upper .
/// (NB. this returns a pointer rather than a reference since it will
/// have to be stored in the analysis class - there's no point in forcing users to explicitly
/// get the pointer from a reference before they can use it!)
AIDA::IHistogram1D* bookHistogram1D(const std::string& name,
- const size_t nbins, const double lower, const double upper,
+ size_t nbins, double lower, double upper,
const std::string& title="",
const std::string& xtitle="", const std::string& ytitle="");
@@ -288,8 +294,8 @@
/// Book a 1D histogram based on the paper, dataset and x/y-axis IDs in the corresponding
/// HepData record. The binnings will be obtained by reading the bundled AIDA data record file
/// of the same filename as the analysis' name() property.
- AIDA::IHistogram1D* bookHistogram1D(const size_t datasetId, const size_t xAxisId,
- const size_t yAxisId, const std::string& title="",
+ AIDA::IHistogram1D* bookHistogram1D(size_t datasetId, size_t xAxisId, size_t yAxisId,
+ const std::string& title="",
const std::string& xtitle="", const std::string& ytitle="");
//@}
@@ -303,7 +309,7 @@
/// have to be stored in the analysis class - there's no point in forcing users to explicitly
/// get the pointer from a reference before they can use it!)
AIDA::IProfile1D* bookProfile1D(const std::string& name,
- const size_t nbins, const double lower, const double upper,
+ size_t nbins, double lower, double upper,
const std::string& title="",
const std::string& xtitle="", const std::string& ytitle="");
@@ -325,8 +331,8 @@
/// Book a 1D profile histogram based on the paper, dataset and x/y-axis IDs in the corresponding
/// HepData record. The binnings will be obtained by reading the bundled AIDA data record file
/// of the same filename as the analysis' name() property.
- AIDA::IProfile1D* bookProfile1D(const size_t datasetId, const size_t xAxisId,
- const size_t yAxisId, const std::string& title="",
+ AIDA::IProfile1D* bookProfile1D(size_t datasetId, size_t xAxisId, size_t yAxisId,
+ const std::string& title="",
const std::string& xtitle="", const std::string& ytitle="");
//@}
@@ -347,7 +353,7 @@
/// have to be stored in the analysis class - there's no point in forcing users to explicitly
/// get the pointer from a reference before they can use it!)
AIDA::IDataPointSet* bookDataPointSet(const std::string& name,
- const size_t npts, const double lower, const double upper,
+ size_t npts, double lower, double upper,
const std::string& title="",
const std::string& xtitle="", const std::string& ytitle="");
@@ -360,9 +366,10 @@
/// Book a 2-dimensional data point set based on the paper, dataset and x/y-axis IDs in the corresponding
/// HepData record. The binnings (x-errors) will be obtained by reading the bundled AIDA data record file
/// of the same filename as the analysis' name() property.
- AIDA::IDataPointSet* bookDataPointSet(const size_t datasetId, const size_t xAxisId,
- const size_t yAxisId, const std::string& title="",
+ AIDA::IDataPointSet* bookDataPointSet(size_t datasetId, size_t xAxisId, size_t yAxisId,
+ const std::string& title="",
const std::string& xtitle="", const std::string& ytitle="");
+
//@}
@@ -375,13 +382,10 @@
void _makeHistoDir();
/// Get the bin edges for this paper from the reference AIDA file, and cache them.
- void _cacheBinEdges();
+ void _cacheBinEdges() const;
/// Get the x-axis points for this paper from the reference AIDA file, and cache them.
- void _cacheXAxisData();
-
- /// Make the axis code string (dsDD-xXX-yYY)
- string _makeAxisCode(const size_t datasetId, const size_t xAxisId, const size_t yAxisId) const;
+ void _cacheXAxisData() const;
//@}
@@ -418,17 +422,17 @@
AnalysisHandler* _analysishandler;
/// Flag to indicate whether the histogram directory is present
- bool _madeHistoDir;
+ mutable bool _madeHistoDir;
/// Collection of x-axis point data to speed up many autobookings: the
/// reference data file should only be read once.
/// @todo Reduce memory occupancy, or clear after initialisation?
- map<string, vector<DPSXPoint> > _dpsData;
+ mutable map<string, vector<DPSXPoint> > _dpsData;
/// Collection of cached bin edges to speed up many autobookings: the
/// reference data file should only be read once.
/// @todo Reduce memory occupancy, or clear after initialisation?
- map<string, BinEdges> _histBinEdges;
+ mutable map<string, BinEdges> _histBinEdges;
private:
Modified: trunk/include/Rivet/Projections/InvMassFinalState.hh
==============================================================================
--- trunk/include/Rivet/Projections/InvMassFinalState.hh Sun Dec 13 16:39:36 2009 (r2185)
+++ trunk/include/Rivet/Projections/InvMassFinalState.hh Sun Dec 13 20:39:43 2009 (r2186)
@@ -9,18 +9,17 @@
/// Identify particles which can be paired to make an invariant mass within a given mass window
class InvMassFinalState : public FinalState {
-
public:
// Constructor for a single inv-mass pair
InvMassFinalState(const FinalState& fsp,
- const std::pair<long, long>& idpair, // pair of decay products
+ const std::pair<PdgId, PdgId>& idpair, // pair of decay products
double minmass, // min inv mass
double maxmass); // max inv mass
InvMassFinalState(const FinalState& fsp,
- const std::vector<std::pair<long, long> >& idpairs, // vector of pairs of decay products
+ const std::vector<std::pair<PdgId, PdgId> >& idpairs, // vector of pairs of decay products
double minmass, // min inv mass
double maxmass); // max inv mass
@@ -29,6 +28,12 @@
virtual const Projection* clone() const {
return new InvMassFinalState(*this);
}
+
+
+ public:
+
+ /// Constituent pairs
+ const std::vector<std::pair<Particle, Particle> >& particlePairs() const;
protected:
@@ -42,13 +47,19 @@
private:
- /// ids of the decay products
- std::vector<std::pair<long, long> > _decayids;
+ /// Handy typedef for a pair of PID codes
+ typedef pair<PdgId,PdgId> PidPair;
+
+ /// IDs of the decay products
+ std::vector<PidPair> _decayids;
+
+ /// Constituent pairs
+ std::vector<std::pair<Particle, Particle> > _particlePairs;
- /// min inv mass
+ /// Min inv mass
double _minmass;
- /// max inv mass
+ /// Max inv mass
double _maxmass;
};
Modified: trunk/include/Rivet/Projections/WFinder.hh
==============================================================================
--- trunk/include/Rivet/Projections/WFinder.hh Sun Dec 13 16:39:36 2009 (r2185)
+++ trunk/include/Rivet/Projections/WFinder.hh Sun Dec 13 20:39:43 2009 (r2186)
@@ -12,10 +12,9 @@
namespace Rivet {
- /// Chain together different projections as convenience for finding Z's
+ /// Chain together different projections as convenience for finding W's
/// from two leptons in the final state
class WFinder : public FinalState {
-
public:
/// @name Constructors
@@ -57,15 +56,16 @@
//@}
- /// Access to the remaining particles, after the Z and clustered photons
+ /// Access to the remaining particles, after the W and clustered photons
/// have been removed from the full final state
/// (e.g. for running a jet finder on it)
const FinalState& remainingFinalState() const;
- /// Access to the Z constituent leptons final state
+ /// Access to the W constituent leptons final state
/// (e.g. for more fine-grained cuts on the leptons)
const FinalState& constituentsFinalState() const;
+
protected:
/// Apply the projection on the supplied event.
@@ -89,6 +89,12 @@
double m2_min, double m2_max,
double dRmax);
+
+ private:
+
+ // Mass range
+ double _m2_min, _m2_max;
+
};
Modified: trunk/src/Analyses/CDF_2006_S6653332.cc
==============================================================================
--- trunk/src/Analyses/CDF_2006_S6653332.cc Sun Dec 13 16:39:36 2009 (r2185)
+++ trunk/src/Analyses/CDF_2006_S6653332.cc Sun Dec 13 20:39:43 2009 (r2186)
@@ -43,7 +43,7 @@
// Create a final state with any e+e- or mu+mu- pair with
// invariant mass 76 -> 106 GeV and ET > 20 (Z decay products)
- vector<pair<long,long> > vids;
+ vector<pair<PdgId,PdgId> > vids;
vids.push_back(make_pair(ELECTRON, POSITRON));
vids.push_back(make_pair(MUON, ANTIMUON));
FinalState fs2(-3.6, 3.6);
Modified: trunk/src/Analyses/CDF_2008_S7541902.cc
==============================================================================
--- trunk/src/Analyses/CDF_2008_S7541902.cc Sun Dec 13 16:39:36 2009 (r2185)
+++ trunk/src/Analyses/CDF_2008_S7541902.cc Sun Dec 13 20:39:43 2009 (r2186)
@@ -43,7 +43,7 @@
addProjection(fs, "FS");
// Create a final state with any e-nu pair with invariant mass 65 -> 95 GeV and ET > 20 (W decay products)
- vector<pair<long,long> > vids;
+ vector<pair<PdgId,PdgId> > vids;
vids += make_pair(ELECTRON, NU_EBAR);
vids += make_pair(POSITRON, NU_E);
FinalState fs2(-3.6, 3.6, 20*GeV);
Modified: trunk/src/Analyses/CDF_2008_S8095620.cc
==============================================================================
--- trunk/src/Analyses/CDF_2008_S8095620.cc Sun Dec 13 16:39:36 2009 (r2185)
+++ trunk/src/Analyses/CDF_2008_S8095620.cc Sun Dec 13 20:39:43 2009 (r2186)
@@ -35,7 +35,7 @@
addProjection(fs, "FS");
// Create a final state with any e+e- or mu+mu- pair with
// invariant mass 76 -> 106 GeV and ET > 20 (Z decay products)
- vector<pair<long,long> > vids;
+ vector<pair<PdgId,PdgId> > vids;
vids.push_back(make_pair(ELECTRON, POSITRON));
vids.push_back(make_pair(MUON, ANTIMUON));
FinalState fs2(-3.6, 3.6);
Modified: trunk/src/Analyses/D0_2008_S7837160.cc
==============================================================================
--- trunk/src/Analyses/D0_2008_S7837160.cc Sun Dec 13 16:39:36 2009 (r2185)
+++ trunk/src/Analyses/D0_2008_S7837160.cc Sun Dec 13 20:39:43 2009 (r2186)
@@ -3,10 +3,14 @@
#include "Rivet/Tools/Logging.hh"
#include "Rivet/Tools/ParticleIdUtils.hh"
#include "Rivet/Projections/FinalState.hh"
+#include "Rivet/Projections/WFinder.hh"
#include "Rivet/Projections/LeadingParticlesFinalState.hh"
#include "Rivet/Projections/IdentifiedFinalState.hh"
#include "Rivet/RivetAIDA.hh"
+#include "LWH/Histogram1D.h"
+#include "LWH/HistogramFactory.h"
+
namespace Rivet {
@@ -31,110 +35,72 @@
// Book histograms and set up projections
void init() {
- // Leading electrons
- FinalState fs(-5.0, 5.0);
-
- LeadingParticlesFinalState efs(fs);
- efs.addParticleId(ELECTRON).addParticleId(POSITRON);
- addProjection(efs, "WDecayE");
-
- LeadingParticlesFinalState nufs(fs);
- nufs.addParticleId(NU_E).addParticleId(NU_EBAR);
- addProjection(nufs, "WDecayNu");
-
- // Final state w/o electron
- IdentifiedFinalState ifs(fs);
- ifs.acceptId(PHOTON);
- addProjection(ifs, "PhotonFS");
+ // Projections
+ const WFinder wfe(-5, 5, 0.0*GeV, ELECTRON, 60.0*GeV, 100.0*GeV, 0.2);
+ addProjection(wfe, "WFe");
// Cross-section histograms
- /// @todo Take binnings from data (since division, subtraction etc. must have compatible binnings)
- _h_dsigplus_deta_25_35 = bookHistogram1D("/dsigplus_deta_25_35", 10, 0.0, 3.2);
- _h_dsigminus_deta_25_35 = bookHistogram1D("/dsigminus_deta_25_35", 10, 0.0, 3.2);
- _h_dsigplus_deta_35 = bookHistogram1D("/dsigplus_deta_35", 10, 0.0, 3.2);
- _h_dsigminus_deta_35 = bookHistogram1D("/dsigminus_deta_35", 10, 0.0, 3.2);
- _h_dsigplus_deta_25 = bookHistogram1D("/dsigplus_deta_25", 10, 0.0, 3.2);
- _h_dsigminus_deta_25 = bookHistogram1D("/dsigminus_deta_25", 10, 0.0, 3.2);
+ const BinEdges& edges = binEdges(1,1,1);
+ _h_dsigplus_deta_25_35 = bookHistogram1D("/dsigplus_deta_25_35", edges);
+ _h_dsigminus_deta_25_35 = bookHistogram1D("/dsigminus_deta_25_35", edges);
+ _h_dsigplus_deta_35 = bookHistogram1D("/dsigplus_deta_35", edges);
+ _h_dsigminus_deta_35 = bookHistogram1D("/dsigminus_deta_35", edges);
+ _h_dsigplus_deta_25 = bookHistogram1D("/dsigplus_deta_25", edges);
+ _h_dsigminus_deta_25 = bookHistogram1D("/dsigminus_deta_25", edges);
}
-
+
/// Do the analysis
void analyze(const Event & event) {
- const double weight = event.weight();
-
- /// @todo Use WFinder projection (includes photon summing)
-
- // Find the W decay products
- const FinalState& efs = applyProjection<FinalState>(event, "WDecayE");
- const FinalState& nufs = applyProjection<FinalState>(event, "WDecayNu");
-
- // If there is no e/nu_e pair in the FinalState, skip the event
- if (efs.particles().size() < 1 || nufs.particles().size() < 1) {
- getLog() << Log::DEBUG << "No e/nu_e pair found " << endl;
- vetoEvent;
- }
-
- // Identify leading nu and electron
- ParticleVector es = efs.particles();
- sort(es.begin(), es.end(), cmpParticleByEt);
- Particle leading_e = es[0];
- //
- ParticleVector nus = nufs.particles();
- sort(nus.begin(), nus.end(), cmpParticleByEt);
- Particle leading_nu = nus[0];
-
- // Require that the neutrino has Et > 25 GeV
- const FourMomentum nu = leading_nu.momentum();
- if (nu.Et() < 25*GeV) {
- getLog() << Log::DEBUG << "Neutrino fails Et cut" << endl;
+ const WFinder& wf = applyProjection<WFinder>(event, "WFe");
+ if (wf.size() == 0) {
+ getLog() << Log::DEBUG << "No W candidates found: vetoing" << endl;
vetoEvent;
}
-
- // Get "raw" electron 4-momentum and add back in photons that could have radiated from the electron
- FourMomentum e = leading_e.momentum();
- const ParticleVector photons = applyProjection<FinalState>(event, "PhotonFS").particles();
- const double HALO_RADIUS = 0.2;
- foreach (const Particle& p, photons) {
- const double pho_eta = p.momentum().pseudorapidity();
- const double pho_phi = p.momentum().azimuthalAngle();
- if (deltaR(e.pseudorapidity(), e.azimuthalAngle(), pho_eta, pho_phi) < HALO_RADIUS) {
- e += p.momentum();
+
+ // Require that leptons have Et >= 25 GeV
+ FourMomentum p_e;
+ int chg_e = 0;
+ foreach (const Particle& l, wf.constituentsFinalState().particles()) {
+ const FourMomentum pl = l.momentum();
+ if (abs(l.pdgId()) == ELECTRON) {
+ chg_e = PID::threeCharge(l.pdgId());
+ p_e = pl;
+ }
+ if (pl.Et()/GeV < 25.0) {
+ getLog() << Log::DEBUG << l.pdgId() << " fails Et cut" << endl;
+ vetoEvent;
}
}
-
- // Require that the electron has Et > 25 GeV
- if (e.Et() < 25*GeV) {
- getLog() << Log::DEBUG << "Electron fails Et cut" << endl;
- vetoEvent;
- }
-
- const double eta_e = fabs(e.pseudorapidity());
- const double et_e = e.Et();
- const int chg_e = PID::threeCharge(leading_e.pdgId());
+ assert(chg_e != 0);
+
+ const double weight = event.weight();
+ const double eta_e = fabs(p_e.pseudorapidity());
+ const double et_e = p_e.Et();
if (et_e < 35*GeV) {
- // 25 < ET < 35
+ // 25 <= ET < 35
if (chg_e < 0) {
_h_dsigminus_deta_25_35->fill(eta_e, weight);
} else {
_h_dsigplus_deta_25_35->fill(eta_e, weight);
}
} else {
- // ET > 35
+ // ET >= 35
if (chg_e < 0) {
_h_dsigminus_deta_35->fill(eta_e, weight);
} else {
_h_dsigplus_deta_35->fill(eta_e, weight);
}
}
- // Inclusive: ET > 25
+ // Inclusive: ET >= 25
if (chg_e < 0) {
_h_dsigminus_deta_25->fill(eta_e, weight);
} else {
_h_dsigplus_deta_25->fill(eta_e, weight);
}
}
-
-
+
+
/// Finalize
void finalize() {
// Construct asymmetry: (dsig+/deta - dsig-/deta) / (dsig+/deta + dsig-/deta) for each Et region
@@ -150,14 +116,14 @@
IHistogram1D* num35 = hf.subtract("/num35", *_h_dsigplus_deta_35, *_h_dsigminus_deta_35);
IHistogram1D* denom35 = hf.add("/denom35", *_h_dsigplus_deta_35, *_h_dsigminus_deta_35);
assert(num35 && denom35);
- hf.divide(histoDir() + "/d02-x01-y01", *num35, *denom35);
+ hf.divide(histoDir() + "/d01-x01-y02", *num35, *denom35);
hf.destroy(num35);
hf.destroy(denom35);
//
IHistogram1D* num25 = hf.subtract("/num25", *_h_dsigplus_deta_25, *_h_dsigminus_deta_25);
IHistogram1D* denom25 = hf.add("/denom25", *_h_dsigplus_deta_25, *_h_dsigminus_deta_25);
assert(num25 && denom25);
- hf.divide(histoDir() + "/d03-x01-y01", *num25, *denom25);
+ hf.divide(histoDir() + "/d01-x01-y03", *num25, *denom25);
hf.destroy(num25);
hf.destroy(denom25);
@@ -183,7 +149,6 @@
//@}
};
-
// This global object acts as a hook for the plugin system
Modified: trunk/src/Analyses/MC_WANALYSIS.cc
==============================================================================
--- trunk/src/Analyses/MC_WANALYSIS.cc Sun Dec 13 16:39:36 2009 (r2185)
+++ trunk/src/Analyses/MC_WANALYSIS.cc Sun Dec 13 20:39:43 2009 (r2186)
@@ -38,7 +38,7 @@
addProjection(fastjets, "Jets");
// Histos
- _hist_chargemulti = bookHistogram1D("track-n", 50, -0.5, 999.5);
+ _hist_chargemulti = bookHistogram1D("track-n", 50, -0.5, 199.5);
_hist_chargept = bookHistogram1D("track-pt", 20, 0, 20);
/// @todo Use profile plots instead:
_hist_chargemeanpt = bookHistogram1D("track-ptavg", 20, 0, 10);
Modified: trunk/src/Core/Analysis.cc
==============================================================================
--- trunk/src/Core/Analysis.cc Sun Dec 13 16:39:36 2009 (r2185)
+++ trunk/src/Core/Analysis.cc Sun Dec 13 20:39:43 2009 (r2186)
@@ -11,6 +11,26 @@
namespace Rivet {
+ namespace {
+ string makeAxisCode(const size_t datasetId, const size_t xAxisId, const size_t yAxisId) {
+ 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();
+ }
+ }
+
+
+ ////////////////////////
+
+
Analysis::Analysis(const string& name)
: _crossSection(-1.0),
_gotCrossSection(false),
@@ -187,7 +207,7 @@
}
- Analysis& Analysis::setCrossSection(const double& xs) {
+ Analysis& Analysis::setCrossSection(double xs) {
_crossSection = xs;
_gotCrossSection = true;
return *this;
@@ -227,7 +247,7 @@
// Histogramming
- void Analysis::_cacheBinEdges() {
+ void Analysis::_cacheBinEdges() const {
_cacheXAxisData();
if (_histBinEdges.empty()) {
getLog() << Log::TRACE << "Getting histo bin edges from AIDA for paper " << name() << endl;
@@ -236,7 +256,7 @@
}
- void Analysis::_cacheXAxisData() {
+ void Analysis::_cacheXAxisData() const {
if (_dpsData.empty()) {
getLog() << Log::TRACE << "Getting DPS x-axis data from AIDA for paper " << name() << endl;
_dpsData = getDPSXValsErrs(name());
@@ -244,26 +264,32 @@
}
- string Analysis::_makeAxisCode(const size_t datasetId, const size_t xAxisId, const size_t 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();
+ const BinEdges& Analysis::binEdges(const string& hname) const {
+ _cacheBinEdges();
+ getLog() << Log::TRACE << "Using histo bin edges for " << name() << ":" << hname << endl;
+ const BinEdges& edges = _histBinEdges.find(hname)->second;
+ if (getLog().isActive(Log::TRACE)) {
+ stringstream edges_ss;
+ foreach (const double be, edges) {
+ edges_ss << " " << be;
+ }
+ getLog() << Log::TRACE << "Edges:" << edges_ss.str() << endl;
+ }
+ return edges;
+ }
+
+
+ const BinEdges& Analysis::binEdges(size_t datasetId, size_t xAxisId, size_t yAxisId) const {
+ const string hname = makeAxisCode(datasetId, xAxisId, yAxisId);
+ return binEdges(hname);
}
- IHistogram1D* Analysis::bookHistogram1D(const size_t datasetId, const size_t xAxisId,
- const size_t yAxisId, const string& title,
+ IHistogram1D* Analysis::bookHistogram1D(size_t datasetId, size_t xAxisId,
+ size_t yAxisId, const string& title,
const string& xtitle, const string& ytitle)
{
- const string axisCode = _makeAxisCode(datasetId, xAxisId, yAxisId);
+ const string axisCode = makeAxisCode(datasetId, xAxisId, yAxisId);
return bookHistogram1D(axisCode, title, xtitle, ytitle);
}
@@ -272,9 +298,7 @@
const string& xtitle, const string& ytitle)
{
// Get the bin edges (only read the AIDA file once)
- _cacheBinEdges();
- getLog() << Log::TRACE << "Using histo bin edges for " << name() << ":" << hname << endl;
- const BinEdges edges = _histBinEdges.find(hname)->second;
+ const BinEdges edges = binEdges(hname);
_makeHistoDir();
const string path = histoPath(hname);
IHistogram1D* hist = histogramFactory().createHistogram1D(path, title, edges);
@@ -286,7 +310,7 @@
IHistogram1D* Analysis::bookHistogram1D(const string& hname,
- const size_t nbins, const double lower, const double upper,
+ size_t nbins, double lower, double upper,
const string& title,
const string& xtitle, const string& ytitle) {
_makeHistoDir();
@@ -316,10 +340,10 @@
/////////////////
- IProfile1D* Analysis::bookProfile1D(const size_t datasetId, const size_t xAxisId,
- const size_t yAxisId, const string& title,
+ IProfile1D* Analysis::bookProfile1D(size_t datasetId, size_t xAxisId,
+ size_t yAxisId, const string& title,
const string& xtitle, const string& ytitle) {
- const string axisCode = _makeAxisCode(datasetId, xAxisId, yAxisId);
+ const string axisCode = makeAxisCode(datasetId, xAxisId, yAxisId);
return bookProfile1D(axisCode, title, xtitle, ytitle);
}
@@ -328,16 +352,7 @@
const string& xtitle, const string& ytitle)
{
// Get the bin edges (only read the AIDA file once)
- _cacheBinEdges();
- getLog() << Log::TRACE << "Using profile histo bin edges for " << name() << ":" << hname << endl;
- const BinEdges edges = _histBinEdges.find(hname)->second;
- if (getLog().isActive(Log::TRACE)) {
- stringstream edges_ss;
- foreach (const double be, edges) {
- edges_ss << " " << be;
- }
- getLog() << Log::TRACE << "Edges:" << edges_ss.str() << endl;
- }
+ const BinEdges edges = binEdges(hname);
_makeHistoDir();
const string path = histoPath(hname);
IProfile1D* prof = histogramFactory().createProfile1D(path, title, edges);
@@ -349,7 +364,7 @@
IProfile1D* Analysis::bookProfile1D(const string& hname,
- const size_t nbins, const double lower, const double upper,
+ size_t nbins, double lower, double upper,
const string& title,
const string& xtitle, const string& ytitle) {
_makeHistoDir();
@@ -393,7 +408,7 @@
IDataPointSet* Analysis::bookDataPointSet(const string& hname,
- const size_t npts, const double lower, const double upper,
+ size_t npts, double lower, double upper,
const string& title,
const string& xtitle, const string& ytitle) {
IDataPointSet* dps = bookDataPointSet(hname, title, xtitle, ytitle);
@@ -410,13 +425,13 @@
}
- IDataPointSet* Analysis::bookDataPointSet(const size_t datasetId, const size_t xAxisId,
- const size_t yAxisId, const string& title,
+ IDataPointSet* Analysis::bookDataPointSet(size_t datasetId, size_t xAxisId,
+ size_t yAxisId, const string& title,
const string& xtitle, const string& ytitle) {
// Get the bin edges (only read the AIDA file once)
_cacheXAxisData();
// Build the axis code
- const string axisCode = _makeAxisCode(datasetId, xAxisId, yAxisId);
+ const string axisCode = makeAxisCode(datasetId, xAxisId, yAxisId);
//const map<string, vector<DPSXPoint> > xpoints = getDPSXValsErrs(papername);
getLog() << Log::TRACE << "Using DPS x-positions for " << name() << ":" << axisCode << endl;
IDataPointSet* dps = bookDataPointSet(axisCode, title, xtitle, ytitle);
@@ -452,7 +467,7 @@
}
- void Analysis::normalize(AIDA::IHistogram1D*& histo, const double norm) {
+ void Analysis::normalize(AIDA::IHistogram1D*& histo, double norm) {
if (!histo) {
getLog() << Log::ERROR << "Failed to normalise histo=NULL in analysis "
<< name() << "(norm=" << norm << ")" << endl;
@@ -477,7 +492,7 @@
}
- void Analysis::scale(AIDA::IHistogram1D*& histo, const double scale) {
+ void Analysis::scale(AIDA::IHistogram1D*& histo, double scale) {
if (!histo) {
getLog() << Log::ERROR << "Failed to scale histo=NULL in analysis "
<< name() << "(scale=" << scale << ")" << endl;
Modified: trunk/src/Projections/InvMassFinalState.cc
==============================================================================
--- trunk/src/Projections/InvMassFinalState.cc Sun Dec 13 16:39:36 2009 (r2185)
+++ trunk/src/Projections/InvMassFinalState.cc Sun Dec 13 20:39:43 2009 (r2186)
@@ -7,7 +7,7 @@
InvMassFinalState::InvMassFinalState(const FinalState& fsp,
- const std::pair<long, long>& idpair, // pair of decay products
+ const pair<PdgId, PdgId>& idpair, // pair of decay products
double minmass, // min inv mass
double maxmass) // max inv mass
: _minmass(minmass), _maxmass(maxmass)
@@ -19,7 +19,7 @@
InvMassFinalState::InvMassFinalState(const FinalState& fsp,
- const std::vector<std::pair<long, long> >& idpairs, // vector of pairs of decay products
+ const vector<pair<PdgId, PdgId> >& idpairs, // vector of pairs of decay products
double minmass, // min inv mass
double maxmass) // max inv mass
: _decayids(idpairs), _minmass(minmass), _maxmass(maxmass)
@@ -57,6 +57,7 @@
void InvMassFinalState::project(const Event& e) {
_theParticles.clear();
+ _particlePairs.clear();
const FinalState& fs = applyProjection<FinalState>(e, "FS");
// Containers for the particles of type specified in the pair
@@ -64,9 +65,8 @@
vector<const Particle*> type2;
// Get all the particles of the type specified in the pair from the particle list
foreach (const Particle& ipart, fs.particles()) {
- // Loop around possible particle pairs
- typedef pair<long,long> longpair;
- foreach (const longpair& ipair, _decayids) {
+ // Loop around possible particle pairs (typedef needed to keep BOOST_FOREACH happy)
+ foreach (const PidPair& ipair, _decayids) {
if (ipart.pdgId() == ipair.first) {
if (accept(ipart.genParticle())) {
type1 += &ipart;
@@ -89,23 +89,26 @@
foreach (const Particle* i2, type2) {
FourMomentum v4 = i1->momentum() + i2->momentum();
if (v4.mass() > _minmass && v4.mass() < _maxmass) {
- // Avoid duplicates
+ getLog() << Log::DEBUG << "Selecting particles with IDs "
+ << i1->pdgId() << " & " << i2->pdgId()
+ << " and mass = " << v4.mass()/GeV << " GeV" << endl;
+ // Store accepted particles, avoiding duplicates
if (find(tmp.begin(), tmp.end(), i1) == tmp.end()) {
tmp.push_back(i1);
- _theParticles.push_back(*i1);
+ _theParticles += *i1;
}
if (find(tmp.begin(), tmp.end(), i2) == tmp.end()) {
tmp.push_back(i2);
- _theParticles.push_back(*i2);
+ _theParticles += *i2;
}
- getLog() << Log::DEBUG << "Selecting particles with IDs "
- << i1->pdgId() << " & " << i2->pdgId()
- << " and mass = " << v4.mass()/GeV << " GeV" << endl;
+ // Store accepted particle pairs
+ _particlePairs += make_pair(*i1, *i2);
}
}
}
- getLog() << Log::DEBUG << "Selected " << _theParticles.size() << " particles." << endl;
+ getLog() << Log::DEBUG << "Selected " << _theParticles.size() << " particles "
+ << "(" << _particlePairs.size() << " pairs)" << endl;
if (getLog().isActive(Log::TRACE)) {
foreach (const Particle& p, _theParticles) {
getLog() << Log::TRACE << "ID: " << p.pdgId()
@@ -115,4 +118,10 @@
}
+ /// Constituent pairs
+ const std::vector<std::pair<Particle, Particle> >& InvMassFinalState::particlePairs() const {
+ return _particlePairs;
+ }
+
+
}
Modified: trunk/src/Projections/WFinder.cc
==============================================================================
--- trunk/src/Projections/WFinder.cc Sun Dec 13 16:39:36 2009 (r2185)
+++ trunk/src/Projections/WFinder.cc Sun Dec 13 20:39:43 2009 (r2186)
@@ -55,12 +55,14 @@
{
setName("WFinder");
- //addProjection(fs, "FS");
+ // Mass range
+ _m2_min = m2_min;
+ _m2_max = m2_max;
assert(abs(pid) == ELECTRON || abs(pid) == MUON || abs(pid) == TAU);
PdgId nu_pid = abs(pid) + 1;
assert(abs(nu_pid) == NU_E || abs(nu_pid) == NU_MU || abs(nu_pid) == NU_TAU);
- vector<pair<long, long> > l_nu_ids;
+ vector<pair<PdgId, PdgId> > l_nu_ids;
l_nu_ids += make_pair(abs(pid), -abs(nu_pid));
l_nu_ids += make_pair(-abs(pid), abs(nu_pid));
InvMassFinalState imfs(fs, l_nu_ids, m2_min, m2_max);
@@ -128,8 +130,11 @@
pW += photon.momentum();
}
msg << " = " << pW;
- getLog() << Log::DEBUG << msg.str() << endl;
+ // Check mass range again
+ if (!inRange(pW.mass()/GeV, _m2_min, _m2_max)) return;
+ getLog() << Log::DEBUG << msg.str() << endl;
+
// Make W Particle and insert into particles list
const PdgId wpid = (wcharge == 1) ? WPLUSBOSON : WMINUSBOSON;
_theParticles.push_back(Particle(wpid, pW));
More information about the Rivet-svn
mailing list
|