|
[Rivet-svn] r2187 - in trunk: . include/Rivet include/Rivet/Projections pyext src/Analyses src/Core src/Projectionsblackhole at projects.hepforge.org blackhole at projects.hepforge.orgMon Dec 14 14:42:56 GMT 2009
Author: buckley Date: Mon Dec 14 14:42:56 2009 New Revision: 2187 Log: More overhauls of the Run/AnalysisHandler API, and related improvements to various beam property functionality Modified: trunk/ChangeLog trunk/include/Rivet/Analysis.hh trunk/include/Rivet/AnalysisHandler.hh trunk/include/Rivet/Particle.fhh trunk/include/Rivet/Particle.hh trunk/include/Rivet/ParticleName.hh trunk/include/Rivet/Projections/Beam.hh trunk/pyext/Makefile.am trunk/pyext/rivet.i trunk/src/Analyses/UA5_1982_S875503.cc trunk/src/Core/Analysis.cc trunk/src/Core/AnalysisHandler.cc trunk/src/Core/Run.cc trunk/src/Projections/Beam.cc Modified: trunk/ChangeLog ============================================================================== --- trunk/ChangeLog Sun Dec 13 20:39:43 2009 (r2186) +++ trunk/ChangeLog Mon Dec 14 14:42:56 2009 (r2187) @@ -1,3 +1,10 @@ +2009-12-14 Andy Buckley <andy at insectnation.org> + + * Adding more beam configuration features to Beam and + AnalysisHandler: the setRunBeams(...) methods on the latter now + allows a beam configuration for the run to be specified without + using the Run class. + 2009-12-11 Andy Buckley <andy at insectnation.org> * Removing use of PVertex from few remaining analyses. Still used @@ -21,7 +28,7 @@ 2009-12-10 Hendrik Hoeth <hendrik.hoeth at cern.ch> - * propagate SPECIAL and HISTOGRAM sections from .plot files + * Propagate SPECIAL and HISTOGRAM sections from .plot files through compare-histos * STAR_2006_S6860818: <pT> vs particle mass, validate analysis Modified: trunk/include/Rivet/Analysis.hh ============================================================================== --- trunk/include/Rivet/Analysis.hh Sun Dec 13 20:39:43 2009 (r2186) +++ trunk/include/Rivet/Analysis.hh Mon Dec 14 14:42:56 2009 (r2187) @@ -170,7 +170,10 @@ /// @name Run conditions /// Incoming beams for this run - const BeamPair& beams() const; + const ParticlePair& beams() const; + + /// Incoming beam IDs for this run + const BeamPair beamIds() const; /// Centre of mass energy for this run double sqrtS() const; Modified: trunk/include/Rivet/AnalysisHandler.hh ============================================================================== --- trunk/include/Rivet/AnalysisHandler.hh Sun Dec 13 20:39:43 2009 (r2186) +++ trunk/include/Rivet/AnalysisHandler.hh Mon Dec 14 14:42:56 2009 (r2187) @@ -8,6 +8,7 @@ #include "Rivet/Event.hh" #include "Rivet/Analysis.hh" #include "Rivet/AnalysisLoader.hh" +#include "Rivet/Projections/Beam.hh" namespace Rivet { @@ -96,32 +97,33 @@ bool hasCrossSection() const; - /// Get beam IDs for this run, determined from first event - const BeamPair& beams() const { - return _beams; + /// Set beams for this run (as determined from first event) + AnalysisHandler& setRunBeams(const Event& event) { + return setRunBeams(Rivet::beams(event)); } - /// Set beam IDs for this run (as determined from first event) - AnalysisHandler& setBeams(const BeamPair& beams) { + /// Set beams for this run (as determined from beam particles) + AnalysisHandler& setRunBeams(const ParticlePair& beams) { getLog() << Log::DEBUG << "Setting run beams = " << beams << endl; _beams = beams; return *this; } - - /// Get energy for this run, determined from first event - double sqrtS() const { - return _sqrts; + /// Get beam IDs for this run, determined from first event + const ParticlePair& beams() const { + return _beams; } - /// Set energy for this run (as determined from first event) - AnalysisHandler& setSqrtS(double sqrts) { - getLog() << Log::DEBUG << "Setting run sqrt(s) = " << sqrts << endl; - _sqrts = sqrts; - return *this; + /// Get beam IDs for this run, determined from first event + BeamPair beamIds() const { + return Rivet::beamIds(beams()); } - + /// Get energy for this run, determined from first event + double sqrtS() const { + return Rivet::sqrtS(beams()); + } + //@} @@ -157,7 +159,7 @@ AnalysisHandler& removeAnalyses(const std::vector<std::string>& analysisnames); - /// Add an analysis to the run list explicitely + /// Add an analysis to the run list by object AnalysisHandler& addAnalysis(Analysis* analysis) { analysis->_analysishandler = this; _analyses.insert(analysis); @@ -173,12 +175,15 @@ /// @name Main init/execute/finalise //@{ - /// Initialize a run. If this run is to be joined together with other - /// runs, \a N should be set to the total number of runs to be - /// combined, and \a i should be the index of this run. This function - /// will initialize the histogram factory and then call the - /// AnalysisBase::init() function of all included analysis objects. - void init(int i=0, int N=0); + /// Initialize a run (beam configuration should be specified first, via @c{setRunBeams}). + void init(); + + + /// Initialize a run, with the run beams taken from the example event. + void init(const GenEvent& event) { + setRunBeams(event); + init(); + } /// Analyze the given \a event. This function will call the @@ -232,13 +237,6 @@ /// Run name std::string _runname; - /// If non-zero the number of runs to be combined into one analysis. - int _nRun; - - /// If non-zero, the index of this run, if a part of several runs to - /// be combined into one. - int _iRun; - /// Number of events seen. size_t _numEvents; @@ -249,11 +247,8 @@ double _xs; /// Beams used by this run. - BeamPair _beams; - - /// Centre of mass energy for this run - double _sqrts; - + ParticlePair _beams; + //@} Modified: trunk/include/Rivet/Particle.fhh ============================================================================== --- trunk/include/Rivet/Particle.fhh Sun Dec 13 20:39:43 2009 (r2186) +++ trunk/include/Rivet/Particle.fhh Mon Dec 14 14:42:56 2009 (r2187) @@ -4,9 +4,9 @@ #include "Rivet/Rivet.hh" - namespace Rivet { + class Particle; /** Typedef a vector of Particle objects. */ @@ -15,6 +15,15 @@ /** Typedef a pair of Particle objects. */ typedef std::pair<Particle, Particle> ParticlePair; + /// Typedef for a PDG ID code. + typedef int PdgId; + + /// Typedef for a pair of particle names. + typedef std::pair<PdgId, PdgId> PidPair; + + /// Typedef for a pair of beam particle names. + typedef std::pair<PdgId, PdgId> BeamPair; + } #endif Modified: trunk/include/Rivet/Particle.hh ============================================================================== --- trunk/include/Rivet/Particle.hh Sun Dec 13 20:39:43 2009 (r2186) +++ trunk/include/Rivet/Particle.hh Mon Dec 14 14:42:56 2009 (r2187) @@ -101,51 +101,88 @@ }; + /// @name String representation + //@{ + /// Print a ParticlePair as a string. + inline std::string toString(const ParticlePair& pair) { + stringstream out; + out << "[" + << toParticleName(pair.first.pdgId()) << " @ " + << pair.first.momentum().E()/GeV << " GeV, " + << toParticleName(pair.second.pdgId()) << " @ " + << pair.second.momentum().E()/GeV << " GeV]"; + return out.str(); + } + + /// Allow ParticlePair to be passed to an ostream. + inline std::ostream& operator<<(std::ostream& os, const ParticlePair& pp) { + os << toString(pp); + return os; + } + + //@} + + + /// @name Comparison functors + //@{ + /// Sort by descending transverse momentum, \f$ p_\perp \f$ inline bool cmpParticleByPt(const Particle& a, const Particle& b) { return a.momentum().pT() > b.momentum().pT(); } + /// Sort by ascending transverse momentum, \f$ p_\perp \f$ inline bool cmpParticleByAscPt(const Particle& a, const Particle& b) { return a.momentum().pT() < b.momentum().pT(); } + /// Sort by descending transverse energy, \f$ E_\perp \f$ inline bool cmpParticleByEt(const Particle& a, const Particle& b) { return a.momentum().Et() > b.momentum().Et(); } + /// Sort by ascending transverse energy, \f$ E_\perp \f$ inline bool cmpParticleByAscEt(const Particle& a, const Particle& b) { return a.momentum().Et() < b.momentum().Et(); } + /// Sort by descending energy, \f$ E \f$ inline bool cmpParticleByE(const Particle& a, const Particle& b) { return a.momentum().E() > b.momentum().E(); } + /// Sort by ascending energy, \f$ E \f$ inline bool cmpParticleByAscE(const Particle& a, const Particle& b) { return a.momentum().E() < b.momentum().E(); } + /// Sort by descending pseudorapidity, \f$ \eta \f$ inline bool cmpParticleByDescPseudorapidity(const Particle& a, const Particle& b) { return a.momentum().pseudorapidity() > b.momentum().pseudorapidity(); } + /// Sort by ascending pseudorapidity, \f$ \eta \f$ inline bool cmpParticleByAscPseudorapidity(const Particle& a, const Particle& b) { return a.momentum().pseudorapidity() < b.momentum().pseudorapidity(); } + /// Sort by descending absolute pseudorapidity, \f$ |\eta| \f$ inline bool cmpParticleByDescAbsPseudorapidity(const Particle& a, const Particle& b) { return fabs(a.momentum().pseudorapidity()) > fabs(b.momentum().pseudorapidity()); } + /// Sort by ascending absolute pseudorapidity, \f$ |\eta| \f$ inline bool cmpParticleByAscAbsPseudorapidity(const Particle& a, const Particle& b) { return fabs(a.momentum().pseudorapidity()) < fabs(b.momentum().pseudorapidity()); } + /// Sort by descending rapidity, \f$ y \f$ inline bool cmpParticleByDescRapidity(const Particle& a, const Particle& b) { return a.momentum().rapidity() > b.momentum().rapidity(); } + /// Sort by ascending rapidity, \f$ y \f$ inline bool cmpParticleByAscRapidity(const Particle& a, const Particle& b) { return a.momentum().rapidity() < b.momentum().rapidity(); } + /// Sort by descending absolute rapidity, \f$ |y| \f$ inline bool cmpParticleByDescAbsRapidity(const Particle& a, const Particle& b) { return fabs(a.momentum().rapidity()) > fabs(b.momentum().rapidity()); } + /// Sort by ascending absolute rapidity, \f$ |y| \f$ inline bool cmpParticleByAscAbsRapidity(const Particle& a, const Particle& b) { return fabs(a.momentum().rapidity()) < fabs(b.momentum().rapidity()); } - - + //@} } Modified: trunk/include/Rivet/ParticleName.hh ============================================================================== --- trunk/include/Rivet/ParticleName.hh Sun Dec 13 20:39:43 2009 (r2186) +++ trunk/include/Rivet/ParticleName.hh Mon Dec 14 14:42:56 2009 (r2187) @@ -2,9 +2,11 @@ #define RIVET_PARTICLENAME_HH #include "Rivet/Rivet.hh" +#include "Rivet/Particle.fhh" namespace Rivet { + /// Enumeration of available beam particles (using PDG IDs where available) enum ParticleName { ELECTRON = 11, @@ -61,8 +63,6 @@ PHOTOANTITAU }; - /// Typedef for a PDG ID code. - typedef int PdgId; /// Convenience maker of particle ID pairs. inline std::pair<PdgId,PdgId> make_pdgid_pair(PdgId a, PdgId b) { @@ -192,14 +192,10 @@ return os; } + ///////////////////////////////////////////////// // Beams - - /// Typedef for a pair of beam particle names. - typedef std::pair<PdgId, PdgId> BeamPair; - - /// Print a BeamPair as a string. inline std::string toString(const BeamPair& pair) { string out = "[" + Modified: trunk/include/Rivet/Projections/Beam.hh ============================================================================== --- trunk/include/Rivet/Projections/Beam.hh Sun Dec 13 20:39:43 2009 (r2186) +++ trunk/include/Rivet/Projections/Beam.hh Mon Dec 14 14:42:56 2009 (r2187) @@ -9,9 +9,34 @@ namespace Rivet { + /// @name Stand-alone functions + //@{ + + /// Function to get beam particles from an event + ParticlePair beams(const Event& e); + + /// Function to get beam particle IDs from an event + BeamPair beamIds(const Event& e); + + /// Function to get beam particle IDs from a pair of particles + BeamPair beamIds(const ParticlePair& beams); + + /// Function to get beam centre of mass energy from an event + double sqrtS(const Event& e); + + /// Function to get beam centre of mass energy from a pair of particles + double sqrtS(const ParticlePair& beams); + + /// Function to get beam centre of mass energy from a pair of beam momenta + double sqrtS(const FourMomentum& pa, const FourMomentum& pb); + + //@} + + + + /// Project out the incoming beams - class Beam : public Projection { - + class Beam : public Projection { public: /// The default constructor. @@ -33,13 +58,12 @@ } /// The pair of beam particle PDG codes in the current collision. - const BeamPair beamIDs() const { - return make_pair(beams().first.pdgId(), - beams().second.pdgId()); + const BeamPair beamIds() const { + return Rivet::beamIds(beams()); } /// Get centre of mass energy, \f$ \sqrt{s} \f$. - const double sqrtS() const; + double sqrtS() const; public: @@ -57,29 +81,13 @@ private: - /// The beam particles in the current collision in GenEvent + + /// The beam particles in the current collision ParticlePair _theBeams; }; - /////////////////////////////////////////////////////// - - - /// @name Stand-alone functions - //@{ - - /// Function to get beam particles from an event - ParticlePair beams(const Event& e); - - /// Function to get beam particle IDs from an event - BeamPair beamIds(const Event& e); - - /// Function to get beam centre of mass energy from an event - double sqrtS(const Event& e); - - //@} - } #endif Modified: trunk/pyext/Makefile.am ============================================================================== --- trunk/pyext/Makefile.am Sun Dec 13 20:39:43 2009 (r2186) +++ trunk/pyext/Makefile.am Mon Dec 14 14:42:56 2009 (r2187) @@ -1,10 +1,7 @@ EXTRA_DIST = rivet.i ez_setup.py rivet_wrap.cc rivet.py: rivet.i - if test ! -e rivet_wrap.cc -o ! -e rivet.py; then \ - rm -f rivet_wrap.cc rivet.py; \ - swig -c++ -python -I$(top_srcdir)/include -o rivet_wrap.cc $<; \ - fi + swig -c++ -python -I$(top_srcdir)/include -o rivet_wrap.cc $<; if ENABLE_PYEXT Modified: trunk/pyext/rivet.i ============================================================================== --- trunk/pyext/rivet.i Sun Dec 13 20:39:43 2009 (r2186) +++ trunk/pyext/rivet.i Mon Dec 14 14:42:56 2009 (r2187) @@ -26,8 +26,9 @@ // Particle ID stuff +%include "Rivet/Particle.fhh" %include "Rivet/ParticleName.hh" -%template(BeamPair) std::pair<long,long>; +//%template(BeamPair) std::pair<PdgId,PdgId>; // Logging interface @@ -70,6 +71,7 @@ double sqrtS(const Event& e); + // Mapping of just the metadata parts of the Analysis API class Analysis { public: virtual std::string name() const; @@ -84,7 +86,9 @@ virtual const std::vector<std::pair<double,double> >& energies() const; virtual std::vector<std::string> authors() const; virtual std::vector<std::string> references() const; - virtual const BeamPair& beams() const; + // virtual const ParticlePair& beams() const; + // virtual const BeamPair& beamIds() const; + // virtual double sqrtS() const; virtual const BeamPair& requiredBeams() const; virtual const bool isCompatible(const ParticleName& beam1, const ParticleName& beam2) const; @@ -98,20 +102,23 @@ class AnalysisHandler { public: - AnalysisHandler(std::string basefilename="Rivet", std::string runname="", + AnalysisHandler(std::string basefilename="Rivet", + std::string runname="", HistoFormat storetype=AIDAML); std::string runName() const; size_t numEvents() const; double sumOfWeights() const; double sqrtS() const; - const BeamPair& beams() const; + const ParticlePair& beams() const; + const BeamPair& beamIds() const; std::vector<std::string> analysisNames(); AnalysisHandler& addAnalysis(const std::string& analysisname); AnalysisHandler& addAnalyses(const std::vector<std::string>& analysisnames); AnalysisHandler& removeAnalysis(const std::string& analysisname); AnalysisHandler& removeAnalyses(const std::vector<std::string>& analysisnames); AnalysisHandler& removeIncompatibleAnalyses(const BeamPair& beams); - void init(int i=0, int N=0); + void init(); + void init(const HepMC::GenEvent& event); void analyze(const HepMC::GenEvent& event); void finalize(); bool needCrossSection(); @@ -126,6 +133,7 @@ static Analysis* getAnalysis(const std::string& analysisname); }; + } %include "Rivet/Run.hh" Modified: trunk/src/Analyses/UA5_1982_S875503.cc ============================================================================== --- trunk/src/Analyses/UA5_1982_S875503.cc Sun Dec 13 20:39:43 2009 (r2186) +++ trunk/src/Analyses/UA5_1982_S875503.cc Mon Dec 14 14:42:56 2009 (r2187) @@ -25,7 +25,7 @@ addProjection(ChargedFinalState(-3.5, 3.5), "CFS"); // Book histos based on pp or ppbar beams - if (beams().first == beams().second) { + if (beamIds().first == beamIds().second) { _hist_nch = bookHistogram1D(2,1,1); _hist_eta = bookHistogram1D(3,1,1); } else { @@ -59,7 +59,7 @@ void finalize() { /// @todo Why the factor of 2 on Nch for ppbar? - if (beams().first == beams().second) { + if (beamIds().first == beamIds().second) { scale(_hist_nch, 1.0/_sumWTrig); } else { scale(_hist_nch, 0.5/_sumWTrig); Modified: trunk/src/Core/Analysis.cc ============================================================================== --- trunk/src/Core/Analysis.cc Sun Dec 13 20:39:43 2009 (r2186) +++ trunk/src/Core/Analysis.cc Mon Dec 14 14:42:56 2009 (r2187) @@ -76,10 +76,14 @@ return handler().sqrtS(); } - const BeamPair& Analysis::beams() const { + const ParticlePair& Analysis::beams() const { return handler().beams(); } + const BeamPair Analysis::beamIds() const { + return handler().beamIds(); + } + const string Analysis::histoDir() const { string path = "/" + name(); Modified: trunk/src/Core/AnalysisHandler.cc ============================================================================== --- trunk/src/Core/AnalysisHandler.cc Sun Dec 13 20:39:43 2009 (r2186) +++ trunk/src/Core/AnalysisHandler.cc Mon Dec 14 14:42:56 2009 (r2187) @@ -12,9 +12,8 @@ AnalysisHandler::AnalysisHandler(string basefilename, string runname, HistoFormat storetype) - : _runname(runname), _nRun(0), _iRun(0), _numEvents(0), - _sumOfWeights(0.0), _xs(-1.0), - _beams(BeamPair(ANY,ANY)), _sqrts(-1.0) + : _runname(runname), _numEvents(0), + _sumOfWeights(0.0), _xs(-1.0) { _theAnalysisFactory = createAnalysisFactory(); _setupFactories(basefilename, storetype); @@ -23,9 +22,8 @@ AnalysisHandler::AnalysisHandler(IAnalysisFactory& afac, string basefilename, string runname, HistoFormat storetype) - : _runname(runname), _nRun(0), _iRun(0), _numEvents(0), + : _runname(runname), _numEvents(0), _sumOfWeights(0.0), _xs(-1.0), - _beams(BeamPair(ANY,ANY)), _sqrts(-1.0), _theAnalysisFactory(&afac) { _setupFactories(basefilename, storetype); @@ -33,8 +31,7 @@ AnalysisHandler::~AnalysisHandler() - { - } + { } Log& AnalysisHandler::getLog() { @@ -42,10 +39,8 @@ } - void AnalysisHandler::init(int i, int N) { + void AnalysisHandler::init() { getLog() << Log::DEBUG << "Initialising the analysis handler" << endl; - _nRun = N; - _iRun = i; _numEvents = 0; _sumOfWeights = 0.0; foreach (Analysis* a, _analyses) { Modified: trunk/src/Core/Run.cc ============================================================================== --- trunk/src/Core/Run.cc Sun Dec 13 20:39:43 2009 (r2186) +++ trunk/src/Core/Run.cc Mon Dec 14 14:42:56 2009 (r2187) @@ -71,8 +71,7 @@ Log::getLog("Rivet.Run") << Log::INFO << "First event beams: " << this->beams() << " @ " << this->sqrtS()/GeV << " GeV" << endl; // Pass to analysis handler - _ah.setBeams(_beams); - _ah.setSqrtS(_sqrts); + _ah.setRunBeams(*_evt); // Set cross-section from command line if (_xs >= 0.0) { Modified: trunk/src/Projections/Beam.cc ============================================================================== --- trunk/src/Projections/Beam.cc Sun Dec 13 20:39:43 2009 (r2186) +++ trunk/src/Projections/Beam.cc Mon Dec 14 14:42:56 2009 (r2187) @@ -6,44 +6,6 @@ namespace Rivet { - void Beam::project(const Event& e) { - assert(e.genEvent().particles_size() >= 2); - std::pair<HepMC::GenParticle*, HepMC::GenParticle*> beams = - e.genEvent().beam_particles(); - assert(beams.first); - _theBeams.first = *(beams.first); - assert(beams.second); - _theBeams.second = *(beams.second); - - getLog() << Log::DEBUG << "Beam particle IDs = " - << _theBeams.first.pdgId() << ", " - << _theBeams.second.pdgId() << endl; - } - - - const double Beam::sqrtS() const { - const double mom1 = beams().first.momentum().pz(); - const double mom2 = beams().second.momentum().pz(); - assert(sign(mom1) != sign(mom2)); - double sqrts = 0.0; - if (fuzzyEquals(fabs(mom1), fabs(mom2))) { - sqrts = fabs(mom1) + fabs(mom2); - } else { - /// @todo Implement general sqrt(s) for asymmetric beams... requires particle masses. - getLog() << Log::WARN << "Asymmetric beams: mass treatment still to be implemented!" << endl; - const double E1 = beams().first.momentum().E(); - const double E2 = beams().second.momentum().E(); - sqrts = (E1+E2)*(E1+E2) - (mom1+mom2)*(mom1+mom2); - sqrts = sqrt(sqrts); - } - getLog() << Log::DEBUG << "sqrt(s) = " << sqrts/GeV << " GeV" << endl; - return sqrts; - } - - - ///////////////////////////////////////////////// - - ParticlePair beams(const Event& e) { Beam beamproj; beamproj.project(e); @@ -53,7 +15,11 @@ BeamPair beamIds(const Event& e) { Beam beamproj; beamproj.project(e); - return beamproj.beamIDs(); + return beamproj.beamIds(); + } + + BeamPair beamIds(const ParticlePair& beams) { + return make_pair(beams.first.pdgId(), beams.second.pdgId()); } double sqrtS(const Event& e) { @@ -62,5 +28,41 @@ return beamproj.sqrtS(); } + double sqrtS(const ParticlePair& beams) { + return sqrtS(beams.first.momentum(), beams.second.momentum()); + } + + double sqrtS(const FourMomentum& pa, const FourMomentum& pb) { + const double mom1 = pa.pz(); + const double e1 = pa.E(); + const double mom2 = pb.pz(); + const double e2 = pb.E(); + double sqrts = sqrt( sqr(e1+e2) - sqr(mom1+mom2) ); + return sqrts; + } + + + + ///////////////////////////////////////////// + + + + void Beam::project(const Event& e) { + assert(e.genEvent().particles_size() >= 2); + pair<HepMC::GenParticle*, HepMC::GenParticle*> beams = e.genEvent().beam_particles(); + assert(beams.first && beams.second); + _theBeams.first = *(beams.first); + _theBeams.second = *(beams.second); + getLog() << Log::DEBUG << "Beam particle IDs = " << beamIds() << endl; + } + + + double Beam::sqrtS() const { + double sqrts = Rivet::sqrtS(beams()); + getLog() << Log::DEBUG << "sqrt(s) = " << sqrts/GeV << " GeV" << endl; + return sqrts; + } + + }
More information about the Rivet-svn mailing list |