Page MenuHomeHEPForge
Files F7877701

No OneTemporary
Actions

View Options

diff --git a/Analysis/NLOHepMCFile.cc b/Analysis/NLOHepMCFile.cc
--- a/Analysis/NLOHepMCFile.cc
+++ b/Analysis/NLOHepMCFile.cc
@@ -1,288 +1,285 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the NLOHepMCFile class.
//
#include "NLOHepMCFile.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/EventRecord/Event.h"
#include "ThePEG/EventRecord/SubProcess.h"
#include "ThePEG/EventRecord/SubProcessGroup.h"
#include "HepMC/IO_AsciiParticles.h"
#include "HepMC/IO_GenEvent.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
using namespace ThePEG;
NLOHepMCFile::NLOHepMCFile()
: _remnantId(82), _format(1), _filename(), _unitchoice(),
_geneventPrecision(16), _eventNumber(1) {}
// Cannot copy streams.
// Let doinitrun() take care of their initialization.
NLOHepMCFile::NLOHepMCFile(const NLOHepMCFile & x)
: AnalysisHandler(x),
_remnantId(x._remnantId), _format(x._format),
_filename(x._filename), _hepmcio(), _hepmcdump(),
_unitchoice(x._unitchoice),
_geneventPrecision(x._geneventPrecision),
_eventNumber(x._eventNumber) {}
HepMC::GenEvent * NLOHepMCFile::makeEvent(tEventPtr event, tSubProPtr sub, long no,
Energy eUnit, Length lUnit,
CrossSection xsec, CrossSection xsecErr) const {
// generate beam particles
const PPair& beam = event->incoming();
HepMC::GenParticle * b1 =
HepMCTraits<HepMC::GenEvent>::newParticle(beam.first->momentum(),beam.first->id(),
1,eUnit);
HepMC::GenParticle * b2 =
HepMCTraits<HepMC::GenEvent>::newParticle(beam.second->momentum(),beam.second->id(),
1,eUnit);
// generate remnants
HepMC::GenParticle * r1 =
HepMCTraits<HepMC::GenEvent>::newParticle(beam.first->momentum() -
sub->incoming().first->momentum(),
_remnantId,1,eUnit);
HepMC::GenParticle * r2 =
HepMCTraits<HepMC::GenEvent>::newParticle(beam.second->momentum() -
sub->incoming().second->momentum(),
_remnantId,1,eUnit);
// generate outgoing particles
vector<HepMC::GenParticle*> outgoing;
for ( ParticleVector::const_iterator p = sub->outgoing().begin();
p != sub->outgoing().end(); ++p ) {
outgoing.push_back(HepMCTraits<HepMC::GenEvent>::newParticle((**p).momentum(),(**p).id(),
1,eUnit));
}
// generate one blob vertex
HepMC::GenVertex * vertex = HepMCTraits<HepMC::GenEvent>::newVertex();
HepMCTraits<HepMC::GenEvent>::addIncoming(*vertex,b1);
HepMCTraits<HepMC::GenEvent>::addIncoming(*vertex,b2);
HepMCTraits<HepMC::GenEvent>::addOutgoing(*vertex,r1);
HepMCTraits<HepMC::GenEvent>::addOutgoing(*vertex,r2);
for ( vector<HepMC::GenParticle*>::const_iterator p = outgoing.begin();
p != outgoing.end(); ++p )
HepMCTraits<HepMC::GenEvent>::addOutgoing(*vertex,*p);
HepMC::GenEvent * ev =
HepMCTraits<HepMC::GenEvent>::newEvent(no,event->weight()*sub->groupWeight(),
event->optionalWeights());
HepMCTraits<HepMC::GenEvent>::setUnits(*ev,eUnit,lUnit);
HepMCTraits<HepMC::GenEvent>::setBeamParticles(*ev,b1,b2);
HepMCTraits<HepMC::GenEvent>::addVertex(*ev,vertex);
HepMCTraits<HepMC::GenEvent>::setCrossSection(*ev,xsec/picobarn,
xsecErr/picobarn);
return ev;
}
void NLOHepMCFile::analyze(tEventPtr event, long, int, int) {
Energy eUnit;
Length lUnit;
switch (_unitchoice) {
default: eUnit = GeV; lUnit = millimeter; break;
case 1: eUnit = MeV; lUnit = millimeter; break;
case 2: eUnit = GeV; lUnit = centimeter; break;
case 3: eUnit = MeV; lUnit = centimeter; break;
}
tcEHPtr eh = dynamic_ptr_cast<tcEHPtr>(event->primaryCollision()->handler());
assert(eh);
CrossSection xsec = eh->integratedXSec();
CrossSection xsecErr = eh->integratedXSecErr();
tSubProPtr sub = event->primarySubProcess();
Ptr<SubProcessGroup>::tptr grp =
dynamic_ptr_cast<Ptr<SubProcessGroup>::tptr>(sub);
HepMC::GenEvent * hepmc =
makeEvent(event,sub,_eventNumber,eUnit,lUnit,xsec,xsecErr);
if (_hepmcio)
_hepmcio->write_event(hepmc);
else
hepmc->print(_hepmcdump);
delete hepmc;
if ( grp ) {
for ( SubProcessVector::const_iterator s = grp->dependent().begin();
s != grp->dependent().end(); ++s ) {
- // remove this when proper Rivet NLO capabilities are ready
- ++_eventNumber;
-
hepmc = makeEvent(event,*s,_eventNumber,eUnit,lUnit,xsec,xsecErr);
if (_hepmcio)
_hepmcio->write_event(hepmc);
else
hepmc->print(_hepmcdump);
delete hepmc;
}
}
++_eventNumber;
}
IBPtr NLOHepMCFile::clone() const {
return new_ptr(*this);
}
IBPtr NLOHepMCFile::fullclone() const {
return new_ptr(*this);
}
void NLOHepMCFile::doinitrun() {
AnalysisHandler::doinitrun();
// set default filename unless user-specified name exists
if ( _filename.empty() )
_filename = generator()->filename() + ".hepmc";
switch ( _format ) {
default: {
HepMC::IO_GenEvent * tmpio
= new HepMC::IO_GenEvent(_filename.c_str(), ios::out);
tmpio->precision(_geneventPrecision);
_hepmcio = tmpio;
break;
}
case 2:
_hepmcio = new HepMC::IO_AsciiParticles(_filename.c_str(), ios::out);
break;
case 5:
_hepmcio = 0;
_hepmcdump.open(_filename.c_str());
break;
}
}
void NLOHepMCFile::dofinish() {
if (_hepmcio) {
delete _hepmcio;
_hepmcio = 0;
}
else
_hepmcdump.close();
AnalysisHandler::dofinish();
cout << "\nNLOHepMCFile: generated HepMC output.\n";
}
// If needed, insert default implementations of virtual function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void NLOHepMCFile::persistentOutput(PersistentOStream & os) const {
os << _remnantId << _format << _filename
<< _unitchoice << _geneventPrecision << _eventNumber;
}
void NLOHepMCFile::persistentInput(PersistentIStream & is, int) {
is >> _remnantId >> _format >> _filename
>> _unitchoice >> _geneventPrecision >> _eventNumber;
}
// *** Attention *** The following static variable is needed for the type
// description system in ThePEG. Please check that the template arguments
// are correct (the class and its base class), and that the constructor
// arguments are correct (the class name and the name of the dynamically
// loadable library where the class implementation can be found).
DescribeClass<NLOHepMCFile,AnalysisHandler>
describeHerwigNLOHepMCFile("ThePEG::NLOHepMCFile", "HepMCAnalysis.so");
void NLOHepMCFile::Init() {
static ClassDocumentation<NLOHepMCFile> documentation
("Write hard sub processes or sub process groups to HepMC.");
static Parameter<NLOHepMCFile,long> interfaceRemnantId
("RemnantId",
"Set the PDG id to be used for remnants.",
&NLOHepMCFile::_remnantId, 82, 0, 0,
false, false, Interface::nolimits);
static Switch<NLOHepMCFile,int> interfaceFormat
("Format",
"Output format (1 = GenEvent, 2 = AsciiParticles, 5 = HepMC dump)",
&NLOHepMCFile::_format, 1, false, false);
static SwitchOption interfaceFormatGenEvent
(interfaceFormat,
"GenEvent",
"IO_GenEvent format",
1);
static SwitchOption interfaceFormatAsciiParticles
(interfaceFormat,
"AsciiParticles",
"Deprecated (IO_AsciiParticles format)",
2);
static SwitchOption interfaceFormatDump
(interfaceFormat,
"Dump",
"Event dump (human readable)",
5);
static Parameter<NLOHepMCFile,string> interfaceFilename
("Filename", "Name of the output file",
&NLOHepMCFile::_filename, "");
static Parameter<NLOHepMCFile,unsigned int> interfacePrecision
("Precision",
"Choice of output precision for the GenEvent format "
" (as number of digits).",
&NLOHepMCFile::_geneventPrecision, 16, 6, 16,
false, false, Interface::limited);
static Switch<NLOHepMCFile,int> interfaceUnits
("Units",
"Unit choice for energy and length",
&NLOHepMCFile::_unitchoice, 0, false, false);
static SwitchOption interfaceUnitsGeV_mm
(interfaceUnits,
"GeV_mm",
"Use GeV and mm as units.",
0);
static SwitchOption interfaceUnitsMeV_mm
(interfaceUnits,
"MeV_mm",
"Use MeV and mm as units.",
1);
static SwitchOption interfaceUnitsGeV_cm
(interfaceUnits,
"GeV_cm",
"Use GeV and cm as units.",
2);
static SwitchOption interfaceUnitsMeV_cm
(interfaceUnits,
"MeV_cm",
"Use MeV and cm as units.",
3);
}
diff --git a/Analysis/NLORivetAnalysis.cc b/Analysis/NLORivetAnalysis.cc
--- a/Analysis/NLORivetAnalysis.cc
+++ b/Analysis/NLORivetAnalysis.cc
@@ -1,281 +1,276 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the NLORivetAnalysis class.
//
#include "NLORivetAnalysis.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/ParVector.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Vectors/HepMCConverter.h"
#include "ThePEG/Config/HepMCHelper.h"
#include "ThePEG/EventRecord/Event.h"
#include "ThePEG/EventRecord/SubProcess.h"
#include "ThePEG/EventRecord/SubProcessGroup.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Repository/CurrentGenerator.h"
#include "HepMC/GenEvent.h"
#include "Rivet/AnalysisHandler.hh"
#include "Rivet/Tools/Logging.hh"
#include "HepMC/IO_AsciiParticles.h"
#include "HepMC/IO_GenEvent.h"
#include <config.h>
using namespace ThePEG;
NLORivetAnalysis::NLORivetAnalysis()
: _remnantId(82), _format(1),_unitchoice(),
_geneventPrecision(16), debug(false), _rivet(), _nevent(0) {}
HepMC::GenEvent * NLORivetAnalysis::makeEvent(tEventPtr event, tSubProPtr sub, long no,
Energy eUnit, Length lUnit,
CrossSection xsec, CrossSection xsecErr) const {
// generate beam particles
const PPair& beam = event->incoming();
HepMC::GenParticle * b1 =
HepMCTraits<HepMC::GenEvent>::newParticle(beam.first->momentum(),beam.first->id(),
1,eUnit);
HepMC::GenParticle * b2 =
HepMCTraits<HepMC::GenEvent>::newParticle(beam.second->momentum(),beam.second->id(),
1,eUnit);
// generate remnants
HepMC::GenParticle * r1 =
HepMCTraits<HepMC::GenEvent>::newParticle(beam.first->momentum() -
sub->incoming().first->momentum(),
_remnantId,1,eUnit);
HepMC::GenParticle * r2 =
HepMCTraits<HepMC::GenEvent>::newParticle(beam.second->momentum() -
sub->incoming().second->momentum(),
_remnantId,1,eUnit);
// generate outgoing particles
vector<HepMC::GenParticle*> outgoing;
for ( ParticleVector::const_iterator p = sub->outgoing().begin();
p != sub->outgoing().end(); ++p ) {
outgoing.push_back(HepMCTraits<HepMC::GenEvent>::newParticle((**p).momentum(),(**p).id(),
1,eUnit));
}
// generate one blob vertex
HepMC::GenVertex * vertex = HepMCTraits<HepMC::GenEvent>::newVertex();
HepMCTraits<HepMC::GenEvent>::addIncoming(*vertex,b1);
HepMCTraits<HepMC::GenEvent>::addIncoming(*vertex,b2);
HepMCTraits<HepMC::GenEvent>::addOutgoing(*vertex,r1);
HepMCTraits<HepMC::GenEvent>::addOutgoing(*vertex,r2);
for ( vector<HepMC::GenParticle*>::const_iterator p = outgoing.begin();
p != outgoing.end(); ++p )
HepMCTraits<HepMC::GenEvent>::addOutgoing(*vertex,*p);
HepMC::GenEvent * ev =
HepMCTraits<HepMC::GenEvent>::newEvent(no,event->weight()*sub->groupWeight(),
event->optionalWeights());
HepMCTraits<HepMC::GenEvent>::setUnits(*ev,eUnit,lUnit);
HepMCTraits<HepMC::GenEvent>::setBeamParticles(*ev,b1,b2);
HepMCTraits<HepMC::GenEvent>::addVertex(*ev,vertex);
HepMCTraits<HepMC::GenEvent>::setCrossSection(*ev,xsec/picobarn,
xsecErr/picobarn);
return ev;
}
void NLORivetAnalysis::analyze(ThePEG::tEventPtr event, long ieve, int loop, int state) {
Energy eUnit;
Length lUnit;
switch (_unitchoice) {
default: eUnit = GeV; lUnit = millimeter; break;
case 1: eUnit = MeV; lUnit = millimeter; break;
case 2: eUnit = GeV; lUnit = centimeter; break;
case 3: eUnit = MeV; lUnit = centimeter; break;
}
tcEHPtr eh = dynamic_ptr_cast<tcEHPtr>(event->primaryCollision()->handler());
assert(eh);
CrossSection xsec = eh->integratedXSec();
CrossSection xsecErr = eh->integratedXSecErr();
tSubProPtr sub = event->primarySubProcess();
Ptr<SubProcessGroup>::tptr grp =
dynamic_ptr_cast<Ptr<SubProcessGroup>::tptr>(sub);
AnalysisHandler::analyze(event, ieve, loop, state);
// Rotate to CMS, extract final state particles and call analyze(particles).
// convert to hepmc
HepMC::GenEvent * hepmc =
makeEvent(event,sub,_nevent,eUnit,lUnit,xsec,xsecErr);
CurrentGenerator::Redirect stdout(cout);
if ( _rivet ) _rivet->analyze(*hepmc);
// delete hepmc event
delete hepmc;
if ( grp ) {
for ( SubProcessVector::const_iterator s = grp->dependent().begin();
s != grp->dependent().end(); ++s ) {
- // remove this when proper Rivet NLO capabilities are ready
- ++_nevent;
+ hepmc = makeEvent(event,*s,_nevent,eUnit,lUnit,xsec,xsecErr);
- hepmc = makeEvent(event,*s,_nevent,eUnit,lUnit,xsec,xsecErr);
-
-
-
if ( _rivet ) _rivet->analyze(*hepmc);
// delete hepmc event
delete hepmc;
}
}
++_nevent;
}
ThePEG::IBPtr NLORivetAnalysis::clone() const {
return new_ptr(*this);
}
ThePEG::IBPtr NLORivetAnalysis::fullclone() const {
return new_ptr(*this);
}
void NLORivetAnalysis::persistentOutput(ThePEG::PersistentOStream & os) const {
os << _analyses << filename << debug;
}
void NLORivetAnalysis::persistentInput(ThePEG::PersistentIStream & is, int) {
is >> _analyses >> filename >> debug;
}
ThePEG::ClassDescription<NLORivetAnalysis> NLORivetAnalysis::initNLORivetAnalysis;
// Definition of the static class description member.
void NLORivetAnalysis::Init() {
static ThePEG::ClassDocumentation<NLORivetAnalysis> documentation
("The NLORivetAnalysis class is a simple class to allow analyses"
" from the Rivet library to be called from ThePEG");
static ThePEG::ParVector<NLORivetAnalysis,string> interfaceAnalyses
("Analyses",
"The names of the Rivet analyses to use",
&NLORivetAnalysis::_analyses, -1, "", "","" "",
false, false, ThePEG::Interface::nolimits);
static Parameter<NLORivetAnalysis,long> interfaceRemnantId
("RemnantId",
"Set the PDG id to be used for remnants.",
&NLORivetAnalysis::_remnantId, 82, 0, 0,
false, false, Interface::nolimits);
static Parameter<NLORivetAnalysis,string> interfaceFilename
("Filename",
#if ThePEG_RIVET_VERSION == 1
"The name of the file where the AIDA histograms are put. If empty, "
"the run name will be used instead. '.aida' will in any case be "
"appended to the file name.",
#elif ThePEG_RIVET_VERSION > 1
"The name of the file where the YODA histograms are put. If empty, "
"the run name will be used instead. '.yoda' will in any case be "
"appended to the file name.",
#else
#error "Unknown ThePEG_RIVET_VERSION"
#endif
&NLORivetAnalysis::filename, "", true, false);
static Switch<NLORivetAnalysis,bool> interfaceDebug
("Debug",
"Enable debug information from Rivet",
&NLORivetAnalysis::debug, false, true, false);
static SwitchOption interfaceDebugNo
(interfaceDebug,
"No",
"Disable debug information.",
false);
static SwitchOption interfaceDebugYes
(interfaceDebug,
"Yes",
"Enable debug information from Rivet.",
true);
interfaceAnalyses.rank(10);
}
void NLORivetAnalysis::dofinish() {
AnalysisHandler::dofinish();
if( _nevent > 0 && _rivet ) {
CurrentGenerator::Redirect stdout(cout);
_rivet->setCrossSection(generator()->integratedXSec()/picobarn);
_rivet->finalize();
string fname = filename;
#if ThePEG_RIVET_VERSION == 1
if ( fname.empty() ) fname = generator()->runName() + ".aida";
#elif ThePEG_RIVET_VERSION > 1
if ( fname.empty() ) fname = generator()->runName() + ".yoda";
#else
#error "Unknown ThePEG_RIVET_VERSION"
#endif
_rivet->writeData(fname);
}
delete _rivet;
_rivet = 0;
}
void NLORivetAnalysis::doinit() {
AnalysisHandler::doinit();
if(_analyses.empty())
throw ThePEG::Exception() << "Must have at least one analysis loaded in "
<< "NLORivetAnalysis::doinitrun()"
<< ThePEG::Exception::runerror;
// check that analysis list is available
_rivet = new Rivet::AnalysisHandler; //(fname);
_rivet->addAnalyses(_analyses);
if ( _rivet->analysisNames().size() != _analyses.size() ) {
throw ThePEG::Exception()
<< "Rivet could not find all requested analyses.\n"
<< "Use 'rivet --list-analyses' to check availability.\n"
<< ThePEG::Exception::runerror;
}
delete _rivet;
_rivet = 0;
}
void NLORivetAnalysis::doinitrun() {
AnalysisHandler::doinitrun();
// create NLORivet analysis handler
CurrentGenerator::Redirect stdout(cout);
_rivet = new Rivet::AnalysisHandler; //(fname);
_rivet->addAnalyses(_analyses);
// check that analysis list is still available
if ( _rivet->analysisNames().size() != _analyses.size() ) {
throw ThePEG::Exception()
<< "Rivet could not find all requested analyses.\n"
<< "Use 'rivet --list-analyses' to check availability.\n"
<< ThePEG::Exception::runerror;
}
if ( debug )
Rivet::Log::setLevel("Rivet",Rivet::Log::DEBUG);
}

File Metadata

Mime Type
text/x-diff
Expires
Tue, Nov 19, 4:16 PM (1 d, 15 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3805102
Default Alt Text
(18 KB)

Event Timeline