|
[Rivet-svn] r2130 - in trunk: . bin include/Rivet pyext src/Coreblackhole at projects.hepforge.org blackhole at projects.hepforge.orgWed Dec 2 16:37:20 GMT 2009
Author: buckley Date: Wed Dec 2 16:37:20 2009 New Revision: 2130 Log: Restructuring of Run, and making first event sqrt(s) and beams accessible from later in the run via the AnalysisHandler. Requires a Python build dir deletion since Python-facing Run and AnalysisHandler APIs have changed. Modified: trunk/ChangeLog trunk/bin/rivet trunk/include/Rivet/AnalysisHandler.hh trunk/include/Rivet/Run.hh trunk/pyext/rivet.i trunk/src/Core/AnalysisHandler.cc trunk/src/Core/Run.cc Modified: trunk/ChangeLog ============================================================================== --- trunk/ChangeLog Wed Dec 2 15:37:56 2009 (r2129) +++ trunk/ChangeLog Wed Dec 2 16:37:20 2009 (r2130) @@ -1,3 +1,12 @@ +2009-12-02 Andy Buckley <andy at insectnation.org> + + * Adding passing of first event sqrt(s) and beams to analysis handler. + + * Restructuring running to only use one HepMC input file (no-one + was using multiple ones, right?) and to break down the Run class + to cleanly separate the init and event loop phases. End of file is + now neater. + 2009-12-01 Andy Buckley <andy at insectnation.org> * Adding parsing of beam types and pairs of energies from YAML. Modified: trunk/bin/rivet ============================================================================== --- trunk/bin/rivet Wed Dec 2 15:37:56 2009 (r2129) +++ trunk/bin/rivet Wed Dec 2 16:37:20 2009 (r2130) @@ -56,7 +56,7 @@ usage="""Run Rivet analyses on inputted events from file or Unix pipe Examples: - %prog [options] <hepmcfile> [<hepmcfile> ...] + %prog [options] <hepmcfile> my_generator -o myfifo & \ %prog [options] myfifo agile-runmc <genname> -n 100k | %prog [options] @@ -377,42 +377,37 @@ ## Read and process events run = rivet.Run(ah) if opts.CROSS_SECTION is not None: + logging.info("User-supplied cross-section = %e pb" % opts.CROSS_SECTION) run.setCrossSection(opts.CROSS_SECTION) if opts.LIST_USED_ANALYSES is not None: run.setListAnalyses(opts.LIST_USED_ANALYSES) -EVTNUM = 0 -starttime = time.time() +## Print platform type +import platform +logging.info("Rivet running on machine %s (%s)" % (platform.node(), platform.machine())) + +## Init run based on one event +evtfile = HEPMCFILES[0] +logging.info("Reading events from '%s'" % evtfile) +if not run.init(evtfile): + log.error("Failed to initialise on event file %s" % evtfile) + sys.exit(2) -for evtfile in HEPMCFILES: - import platform - logging.info("Rivet running on machine %s (%s)" % (platform.node(), platform.machine())) - if opts.CROSS_SECTION is not None: - logging.info("User-supplied cross-section = %e pb" % opts.CROSS_SECTION) - logging.info("Reading events from '%s'" % evtfile) - if not run.prepareFile(evtfile): - break - - while opts.MAXEVTNUM is None or EVTNUM < opts.MAXEVTNUM: - EVTNUM += 1 - logNEvt(EVTNUM, starttime, opts.MAXEVTNUM) - if not run.processEvent(EVTNUM==1): - #logging.warn("run.processEvent failed for evt #%i!" % (EVTNUM)) - break - if RECVD_KILL_SIGNAL is not None: - break - - if not run.finalizeFile(): +## Event loop +starttime = time.time() +EVTNUM = 0 +while opts.MAXEVTNUM is None or EVTNUM < opts.MAXEVTNUM: + EVTNUM += 1 + logNEvt(EVTNUM, starttime, opts.MAXEVTNUM) + if not run.processEvent(): + logging.warn("Event processing failed for evt #%i!" % (EVTNUM)) break - if RECVD_KILL_SIGNAL is not None: break - + if not run.readEvent(): + break logging.info("Finished event loop") - -## Override any generator cross-section to that supplied on the command line, if there was one -if opts.CROSS_SECTION is not None: - run.setCrossSection(opts.CROSS_SECTION) +run.finalize() ## Finalize and write out data file ah.finalize() Modified: trunk/include/Rivet/AnalysisHandler.hh ============================================================================== --- trunk/include/Rivet/AnalysisHandler.hh Wed Dec 2 15:37:56 2009 (r2129) +++ trunk/include/Rivet/AnalysisHandler.hh Wed Dec 2 16:37:20 2009 (r2130) @@ -59,9 +59,13 @@ public: + /// @name Run properties + //@{ + /// Get the name of this run. string runName() const; + /// 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; @@ -74,8 +78,53 @@ /// Set 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) - void setSumOfWeights(const double& sum); - + void setSumOfWeights(const double& sum); + + + /// Is cross-section information required by at least one child analysis? + bool needCrossSection() const; + + /// Set the cross-section for the process being generated. + AnalysisHandler& setCrossSection(double xs); + + /// Get the cross-section known to the handler. + double crossSection() const { + return _xs; + } + + /// Whether the handler knows about a cross-section. + bool hasCrossSection() const; + + + /// Get beam IDs for this run, determined from first event + const BeamPair& beams() const { + return _beams; + } + + /// Set beam IDs for this run (as determined from first event) + AnalysisHandler& setBeams(const BeamPair& beams) { + _beams = beams; + return *this; + } + + + /// Get energy for this run, determined from first event + double sqrtS() const { + return _sqrts; + } + + /// Set energy for this run (as determined from first event) + AnalysisHandler& setSqrtS(double sqrts) { + _sqrts = sqrts; + return *this; + } + + + //@} + + + /// @name handle analyses + //@{ /// Get a list of the currently registered analyses' names. std::vector<std::string> analysisNames(); @@ -111,6 +160,11 @@ /// Remove beam-incompatible analyses from the run list. AnalysisHandler& removeIncompatibleAnalyses(const BeamPair& beams); + //@} + + + /// @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 @@ -131,6 +185,11 @@ /// for further analysis or writing to file is left to the API user. void finalize(); + //@} + + + /// @name AIDA factories etc. + //@{ /// The AIDA analysis factory. AIDA::IAnalysisFactory& analysisFactory(); @@ -151,13 +210,7 @@ /// The AIDA histogram factory. AIDA::IDataPointSetFactory& datapointsetFactory(); - - /// Is cross-section information required by at least one child analysis? - bool needCrossSection() const; - - - /// Set the cross-section for the process being generated. - AnalysisHandler& setCrossSection(double xs); + //@} private: @@ -165,6 +218,10 @@ /// The collection of Analysis objects to be used. set<Analysis*> _analyses; + + /// @name Run properties + //@{ + /// Run name std::string _runname; @@ -181,6 +238,21 @@ /// Sum of event weights seen. double _sumOfWeights; + /// Cross-section known to AH + double _xs; + + /// Beams used by this run. + BeamPair _beams; + + /// Centre of mass energy for this run + double _sqrts; + + //@} + + + /// @name AIDA factory handles + //@{ + /// The AIDA analysis factory. AIDA::IAnalysisFactory* _theAnalysisFactory; @@ -193,6 +265,8 @@ /// The AIDA data point set factory. AIDA::IDataPointSetFactory* _theDataPointSetFactory; + //@} + private: Modified: trunk/include/Rivet/Run.hh ============================================================================== --- trunk/include/Rivet/Run.hh Wed Dec 2 15:37:56 2009 (r2129) +++ trunk/include/Rivet/Run.hh Wed Dec 2 16:37:20 2009 (r2130) @@ -55,22 +55,26 @@ //@{ /// Set up HepMC file readers - bool prepareFile(const std::string& evtfile); + bool init(const std::string& evtfile); + + /// Read the next HepMC event + bool readEvent(); /// Handle next event - bool processEvent(bool firstEvent); + bool processEvent(); /// Close up HepMC I/O - bool finalizeFile(); + bool finalize(); //@} - + private: /// AnalysisHandler object AnalysisHandler& _ah; + /// @name Run variables obtained from events or command line //@{ @@ -93,6 +97,9 @@ /// @name HepMC I/O members //@{ + /// Current event + shared_ptr<GenEvent> _evt; + /// Output stream for HepMC writer shared_ptr<std::istream> _istr; Modified: trunk/pyext/rivet.i ============================================================================== --- trunk/pyext/rivet.i Wed Dec 2 15:37:56 2009 (r2129) +++ trunk/pyext/rivet.i Wed Dec 2 16:37:20 2009 (r2130) @@ -103,6 +103,8 @@ std::string runName() const; size_t numEvents() const; double sumOfWeights() const; + double sqrtS() const; + const BeamPair& beams() const; std::vector<std::string> analysisNames(); AnalysisHandler& addAnalysis(const std::string& analysisname); AnalysisHandler& addAnalyses(const std::vector<std::string>& analysisnames); Modified: trunk/src/Core/AnalysisHandler.cc ============================================================================== --- trunk/src/Core/AnalysisHandler.cc Wed Dec 2 15:37:56 2009 (r2129) +++ trunk/src/Core/AnalysisHandler.cc Wed Dec 2 16:37:20 2009 (r2130) @@ -12,7 +12,9 @@ AnalysisHandler::AnalysisHandler(string basefilename, string runname, HistoFormat storetype) - : _runname(runname), _nRun(0), _iRun(0), _numEvents(0), _sumOfWeights(0.0) { + : _runname(runname), _nRun(0), _iRun(0), _numEvents(0), + _sumOfWeights(0.0), _xs(-1.0) + { _theAnalysisFactory = createAnalysisFactory(); _setupFactories(basefilename, storetype); } @@ -20,8 +22,10 @@ AnalysisHandler::AnalysisHandler(IAnalysisFactory& afac, string basefilename, string runname, HistoFormat storetype) - : _runname(runname), _nRun(0), _iRun(0), _numEvents(0), _sumOfWeights(0.0), - _theAnalysisFactory(&afac) { + : _runname(runname), _nRun(0), _iRun(0), _numEvents(0), + _sumOfWeights(0.0), _xs(-1.0), + _theAnalysisFactory(&afac) + { _setupFactories(basefilename, storetype); } @@ -288,10 +292,17 @@ AnalysisHandler& AnalysisHandler::setCrossSection(double xs) { + _xs = xs; foreach (Analysis* a, _analyses) { a->setCrossSection(xs); } return *this; } + + bool AnalysisHandler::hasCrossSection() const { + return (crossSection() >= 0); + } + + } Modified: trunk/src/Core/Run.cc ============================================================================== --- trunk/src/Core/Run.cc Wed Dec 2 15:37:56 2009 (r2129) +++ trunk/src/Core/Run.cc Wed Dec 2 16:37:20 2009 (r2130) @@ -9,7 +9,7 @@ Run::Run(AnalysisHandler& ah) - : _ah(ah), _xs(-1.0) + : _ah(ah), _xs(-1.0), _sqrts(-1.0) { } @@ -28,7 +28,20 @@ } - bool Run::prepareFile(const std::string& evtfile) { + // Fill event and check for a bad read state + bool Run::readEvent() { + /// @todo Clear rather than new the GenEvent object per-event? + _evt.reset(new GenEvent()); + if (_io->rdstate() != 0 || !_io->fill_next_event(_evt.get()) ) { + //Log::getLog("Rivet.Run") << Log::DEBUG << "Read failed. End of file?" << endl; + return false; + } + return true; + } + + + bool Run::init(const std::string& evtfile) { + // Set up HepMC input reader objects if (evtfile == "-") { _io.reset(new HepMC::IO_GenEvent(std::cin)); } else { @@ -41,103 +54,94 @@ return false; } - return true; - } - + // Read first event to define run conditions + bool ok = readEvent(); + if (!ok) return false; + if (_evt->particles_size() == 0) { + Log::getLog("Rivet.Run") << Log::ERROR << "Empty first event." << endl; + return false; + } - bool Run::processEvent(bool firstEvent) { - // Fill event and check for a bad read state - shared_ptr<GenEvent> evt; - evt.reset(new GenEvent()); - if (_io->rdstate() != 0 || !_io->fill_next_event(evt.get()) ) { - Log::getLog("Rivet.Run") << Log::DEBUG << "Read failed. End of file?" << endl; + // Set required beams for run based on first beams + const BeamPair beams = beamIds(*_evt); + const double sqrts = Rivet::sqrtS(*_evt); + _beams = beams; + _sqrts = sqrts; + 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); + + // Set cross-section from command line + if (_xs >= 0.0) { + Log::getLog("Rivet.Run") + << Log::DEBUG << "Setting user cross-section = " << _xs << " pb" << endl; + _ah.setCrossSection(_xs); + } + + // Check that analyses are beam-compatible + const size_t num_anas_requested = _ah.analysisNames().size(); + _ah.removeIncompatibleAnalyses(beams); + if (num_anas_requested > 0 && _ah.analysisNames().size() == 0) { + Log::getLog("Rivet.Run") << Log::ERROR + << "All analyses were incompatible with the first event's beams\n" + << "Exiting, since this probably isn't intentional!" << endl; return false; } - - // Get beam details from first event, and ensure they match for all following events - if (evt->particles_size() != 0) { - const BeamPair beams = beamIds(*evt); - const double sqrts = Rivet::sqrtS(*evt); - Log::getLog("Rivet.Run") << Log::DEBUG << "Beams: " - << beams << " @ " << sqrts/GeV << " GeV" << endl; - if (firstEvent) { - _beams = beams; - _sqrts = sqrts; - Log::getLog("Rivet.Run") << Log::INFO << "First event beams: " - << this->beams() << " @ " << this->sqrtS()/GeV << " GeV" << endl; - } else { - if (beams != _beams || !fuzzyEquals(sqrts, sqrtS())) { - Log::getLog("Rivet.Run") << Log::ERROR << "Event beams mismatch: " - << beams << " @ " << sqrts/GeV << " GeV" << " vs. first beams " - << this->beams() << " @ " << this->sqrtS()/GeV << " GeV" << endl; - return false; - } + + // List the chosen & compatible analyses if requested + if (_listAnalyses) { + foreach (const std::string& ana, _ah.analysisNames()) { + cout << ana << endl; } } - // Set up system based on properties of first event - if (firstEvent) { - // If empty - if (evt->particles_size() == 0) { - Log::getLog("Rivet.Run") << Log::ERROR << "Empty first event." << endl; - return false; - } + return true; + } - const size_t num_anas_requested = _ah.analysisNames().size(); - _ah.removeIncompatibleAnalyses(beams()); - if (num_anas_requested > 0 && _ah.analysisNames().size() == 0) { - Log::getLog("Rivet.Run") << Log::ERROR - << "All analyses were incompatible with the first event's beams\n" - << "Exiting, since this probably isn't intentional!" << endl; - return false; - } - - if (_listAnalyses) { - foreach (const std::string& ana, _ah.analysisNames()) { - cout << ana << endl; - } - } + bool Run::processEvent() { + // Ensure that beam details match those from first event + const BeamPair beams = beamIds(*_evt); + const double sqrts = Rivet::sqrtS(*_evt); + if (beams != _beams || !fuzzyEquals(sqrts, sqrtS())) { + Log::getLog("Rivet.Run") + << Log::ERROR << "Event beams mismatch: " + << beams << " @ " << sqrts/GeV << " GeV" << " vs. first beams " + << this->beams() << " @ " << this->sqrtS()/GeV << " GeV" << endl; + return false; } - - // Set cross-section if specified from command line - if (_xs > 0.0) { - if (firstEvent) { - Log::getLog("Rivet.Run") << Log::DEBUG - << "Setting user cross-section = " << _xs << " pb" << endl; - _ah.setCrossSection(_xs); - } - } - // Set cross-section if found in event + // Set cross-section if found in event and not from command line #ifdef HEPMC_HAS_CROSS_SECTION - else if (evt->cross_section()) { - const double xs = evt->cross_section()->cross_section(); //< in pb - Log::getLog("Rivet.Run") << Log::DEBUG - << "Setting cross-section = " << xs << " pb" << endl; + if (_xs < 0.0 && _evt->cross_section()) { + const double xs = _evt->cross_section()->cross_section(); //< in pb + Log::getLog("Rivet.Run") + << Log::DEBUG << "Setting cross-section = " << xs << " pb" << endl; _ah.setCrossSection(xs); } #endif // Complain about absence of cross-section if required! - else { - if (_ah.needCrossSection()) { - Log::getLog("Rivet.Run") << Log::ERROR - << "Total cross-section needed for at least one of the analyses. " - << "Please set it (on the command line with '-x' if using the 'rivet' program)" << endl; - return false; - } + if (_ah.needCrossSection() && !_ah.hasCrossSection()) { + Log::getLog("Rivet.Run") + << Log::ERROR + << "Total cross-section needed for at least one of the analyses. " + << "Please set it (on the command line with '-x' if using the 'rivet' program)" << endl; + return false; } - - /// @todo If NOT first event, check that beams aren't changed - + // Analyze event - _ah.analyze(*evt); + _ah.analyze(*_evt); return true; } - bool Run::finalizeFile() { + bool Run::finalize() { + _evt.reset(); + _istr.reset(); + _io.reset(); return true; }
More information about the Rivet-svn mailing list |