Page MenuHomeHEPForge
Files F10664276

No OneTemporary
Actions

View Options

diff --git a/HZTOOL/HZTool.cc b/HZTOOL/HZTool.cc
--- a/HZTOOL/HZTool.cc
+++ b/HZTOOL/HZTool.cc
@@ -1,254 +1,273 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the HZTool class.
//
#include "HZTool.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/ParVector.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Utilities/Throw.h"
#include "ThePEG/CLHEPWrap/HepMCConverter.h"
#include "ThePEG/PDT/StandardMatchers.h"
#include "HepMC/GenEvent.h"
#include "HepMC/IO_HEPEVT.h"
#ifdef ThePEG_TEMPLATES_IN_CC_FILE
// #include "HZTool.tcc"
#endif
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
namespace ThePEG {
template<>
struct HepMCTraits<HepMC::GenEvent>:
public HepMCTraitsBase<HepMC::GenEvent,HepMC::GenParticle,
- HepMC::GenVertex,HepMC::Polarization> {};
+ HepMC::GenVertex,HepMC::Polarization> {
+ /** Create a new particle object with momentum \a p, PDG number \a
+ id and status code \a status. */
+ static ParticleT * newParticle(const Lorentz5Momentum & p,
+ long id, int status) {
+ // Note that according to the documentation the momentum is stored in a
+ // HepLorentzVector in GeV (event though the CLHEP standard is MeV).
+ CLHEP::HepLorentzVector pp(p.x()/GeV, p.y()/GeV, p.z()/GeV, p.t()/GeV);
+ ParticleT * genp = new ParticleT(pp, id, status);
+ genp->setGeneratedMass(p.mass()/GeV);
+ return genp;
+ }
+
+ /** Set the position \a p for the vertex, \a v. */
+ static void setPosition(VertexT & v, const LorentzPoint & p) {
+ // We assume that the position is measured in millimeters.
+ CLHEP::HepLorentzVector pp(p.x()/mm, p.y()/mm, p.z()/mm, p.t()/mm);
+ v.set_position(pp);
+ }
+};
}
extern "C" {
void hztoolinit_(const char *, int);
void hztoolsetxsec_(const double &, const double &);
void hztoolfinish_();
void hztoolredirect_(const char *, int);
void hztoolfixdaughters_();
void hztoolfixbeams_();
void hztoolfixboson_(const double &, const double &, const double &,
const double &, const double &);
void hztoolfixmirdir_();
void hzfilhep_();
void hz95108_(const int &);
void hz96215_(const int &);
void hz98076_(const int &);
void hz98050_(const int &);
void hz98051_(const int &);
void hz98143_(const int &);
void hzf01211e_(const int &);
void hzf89201e_(const int &);
void hzh9807014_(const int &);
void hzh9807018_(const int &);
void hzh9905024_(const int &);
void hzh9907009_(const int &);
void hzh9912022_(const int &);
void hzh0001021_(const int &);
void hzh0010026_(const int &);
void hz0307080_(const int &);
void hzh0412071_(const int &);
}
using namespace ThePEG;
HZTool::HZAnaFn HZTool::getFunctionPointer(string name) {
if ( name == "" ) return 0;
else if ( name == "hz95108" ) return hz95108_;
else if ( name == "hz96215" ) return hz96215_;
else if ( name == "hz98076" ) return hz98076_;
else if ( name == "hz98050" ) return hz98050_;
else if ( name == "hz98051" ) return hz98051_;
else if ( name == "hz98143" ) return hz98143_;
else if ( name == "hzf01211e" ) return hzf01211e_;
else if ( name == "hzf89201e" ) return hzf89201e_;
else if ( name == "hzh9807014" ) return hzh9807014_;
else if ( name == "hzh9807018" ) return hzh9807018_;
else if ( name == "hzh9905024" ) return hzh9905024_;
else if ( name == "hzh9907009" ) return hzh9907009_;
else if ( name == "hzh9912022" ) return hzh9912022_;
else if ( name == "hzh0001021" ) return hzh0001021_;
else if ( name == "hzh0010026" ) return hzh0010026_;
else if ( name == "hz0307080" ) return hz0307080_;
else if ( name == "hzh0412071" ) return hzh0412071_;
else return 0;
}
HZTool::~HZTool() {}
void HZTool::analyze(tEventPtr event, long ieve, int loop, int state) {
sumweight += event->weight();
HepMC::GenEvent * geneve =
HepMCConverter<HepMC::GenEvent>::convert(*event, true);
static HepMC::IO_HEPEVT converter;
converter.set_trust_both_mothers_and_daughters(true);
// static int count = 0;
// if ( count++ < 10 ) geneve->print(cerr);
converter.write_event(geneve);
addDISLepton(*event->primarySubProcess());
hztoolfixbeams_();
hzfilhep_();
hztoolfixdaughters_();
for ( int i = 0, N = functions.size(); i < N; ++i )
if ( functions[i] ) (*functions[i])(2);
delete geneve;
}
LorentzRotation HZTool::transform(tEventPtr event) const {
return LorentzRotation();
// Return the Rotation to the frame in which you want to perform the analysis.
}
void HZTool::analyze(const tPVector & particles) {
AnalysisHandler::analyze(particles);
// Calls analyze() for each particle.
}
void HZTool::analyze(tPPtr) {}
void HZTool::addDISLepton(const SubProcess & sub) {
tcPPtr in;
tcPPtr out;
if ( LeptonMatcher::Check(sub.incoming().first->data()) )
in = sub.incoming().first;
if ( LeptonMatcher::Check(sub.incoming().second->data()) ) {
if ( in ) return;
in = sub.incoming().second;
}
if ( !in ) return;
long idin = in->id();
for ( int i = 0, N = sub.outgoing().size(); i < N; ++i ) {
long idout = sub.outgoing()[i]->id();
if ( !LeptonMatcher::Check(idout) ) continue;
if ( idin == idout ) out = sub.outgoing()[i];
}
if ( !out ) return;
Lorentz5Momentum pgam = in->momentum() - out->momentum();
hztoolfixboson_(pgam.x()/GeV, pgam.y()/GeV, pgam.z()/GeV,
pgam.t()/GeV, pgam.m()/GeV);
if ( in->momentum().z() < 0 )
hztoolfixmirdir_();
}
void HZTool::dofinish() {
AnalysisHandler::dofinish();
hztoolsetxsec_(generator()->integratedXSec()/nanobarn, sumweight);
for ( int i = 0, N = functions.size(); i < N; ++i )
if ( functions[i] ) (*functions[i])(3);
hztoolfinish_();
}
void HZTool::doinitrun() {
sumweight = 0.0;
AnalysisHandler::doinitrun();
string file = ( filename().empty()? generator()->filename(): filename() )
+ ".hzlog";
hztoolredirect_(file.c_str(), file.length());
file = ( filename().empty()? generator()->filename(): filename() )
+ ".rz";
hztoolinit_(file.c_str(), file.length());
for ( int i = 0, N = functions.size(); i < N; ++i )
if ( functions[i] ) (*functions[i])(1);
}
void HZTool::persistentOutput(PersistentOStream & os) const {
os << theFilename << functionNames << sumweight;
}
void HZTool::persistentInput(PersistentIStream & is, int) {
is >> theFilename >> functionNames >> sumweight;
functions.resize(functionNames.size());
for ( int i = 0, N = functionNames.size(); i < N; ++i ) {
functions[i] = getFunctionPointer(functionNames[i]);
if ( !functions[i] ) Throw<MissingFunction>()
<< "While reading '" << fullName() << "': The function '"
<< functionNames[i]
<< "' was not present in the installed HZTool version."
<< Exception::runerror;
}
}
void HZTool::insertFunction(string name, int pos) {
HZAnaFn fp = getFunctionPointer(name);
if ( !fp ) Throw<MissingFunction>()
<< "The function '" << name
<< "' was not present in the installed HZTool version."
<< Exception::eventerror;
pos = max(0, min(int(functionNames.size()), pos));
functionNames.insert(functionNames.begin() + pos, name);
functions.insert(functions.begin() + pos, fp);
}
void HZTool::setFunction(string name, int pos) {
if ( pos < 0 || pos >= int(functionNames.size()) ) return;
HZAnaFn fp = getFunctionPointer(name);
if ( !fp ) Throw<MissingFunction>()
<< "The function '" << name
<< "' was not present in the installed HZTool version."
<< Exception::eventerror;
functionNames[pos] = name;
functions[pos] = fp;
}
void HZTool::delFunction(int pos) {
if ( pos < 0 || pos >= int(functionNames.size()) ) return;
functionNames.erase(functionNames.begin() + pos);
functions.erase(functions.begin() + pos);
}
ClassDescription<HZTool> HZTool::initHZTool;
// Definition of the static class description member.
void HZTool::Init() {
static ClassDocumentation<HZTool> documentation
("This class wraps the HZTool fortran library. The specified "
"<interface>Functions</interface> of in the HZTool library will be "
"called for each event (after the ThePEG::Event is first converted "
"to a HepMC::GenEvent and then translated to the HEPEVT fortran common "
"block). Note that only one HZTool AnalysisHandler can be used for "
"a given EventGenerator, otherwise the result is undefined.");
static Parameter<HZTool,string> interfaceFilename
("Filename",
"The filename to which the HZTool histograms are written. If empty the "
"name of the controlling EventGenerator is used instead. The standard "
"'.rz' suffix is added to the name. Note that since the file is created "
"within the cernlib hbook routines, the actual filename will be in lower "
"case irrespectively of what is specified here. A file with the same "
"name but with the suffix '.hzlog' will also be created containing all "
"standard output from the HZTool library.",
&HZTool::theFilename, "",
true, false);
static ParVector<HZTool,string> interfaceFunctions
("Functions",
"The names of the HZTool analysis routines to be called in this analysis "
"handler. Note that the routine names must be given in lower case "
"letters.",
&HZTool::functionNames, -1, "", "", "",
true, false, Interface::nolimits,
&HZTool::setFunction, &HZTool::insertFunction,
&HZTool::delFunction, (vector<string>(HZTool::*)()const)(0),
(string(HZTool::*)(int)const)(0), (string(HZTool::*)(int)const)(0),
(string(HZTool::*)(int)const)(0), (vector<string>(HZTool::*)()const)(0));
}
diff --git a/HerwigTest/MEqq2gZ2ll.cc b/HerwigTest/MEqq2gZ2ll.cc
--- a/HerwigTest/MEqq2gZ2ll.cc
+++ b/HerwigTest/MEqq2gZ2ll.cc
@@ -1,175 +1,175 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the MEqq2gZ2ll class.
//
#include "MEqq2gZ2ll.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/PDT/EnumParticles.h"
#include "ThePEG/Utilities/Timer.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/StandardModel/StandardModelBase.h"
#include "ThePEG/Handlers/StandardXComb.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
using namespace ThePEG;
using namespace Herwig;
MEqq2gZ2ll::MEqq2gZ2ll()
: coefs(20), mZ2(0.0*GeV2), GZ2(0.0*GeV2), lastCont(0.0), lastBW(0.0) {}
MEqq2gZ2ll::MEqq2gZ2ll(const MEqq2gZ2ll & x)
: ME2to2QCD(x), coefs(x.coefs), mZ2(x.mZ2), GZ2(x.GZ2),
lastCont(x.lastCont), lastBW(x.lastBW) {}
MEqq2gZ2ll::~MEqq2gZ2ll() {}
unsigned int MEqq2gZ2ll::orderInAlphaS() const {
return 0;
}
unsigned int MEqq2gZ2ll::orderInAlphaEW() const {
return 2;
}
void MEqq2gZ2ll::getDiagrams() const {
tcPDPtr gamma = getParticleData(ParticleID::gamma);
tcPDPtr Z0 = getParticleData(ParticleID::Z0);
for(int i = 1; i <= maxFlavour(); ++i) {
tcPDPtr q = getParticleData(i);
tcPDPtr qb = q->CC();
for(long lid = ParticleID::eminus; lid <= ParticleID::tauminus; lid+=2) {
tcPDPtr l = getParticleData(lid);
tcPDPtr lb = l->CC();
add(new_ptr((Tree2toNDiagram(2), q, qb, 1, gamma, 3, l, 3, lb, -1)));
add(new_ptr((Tree2toNDiagram(2), q, qb, 1, Z0, 3, l, 3, lb, -2)));
}
}
}
Energy2 MEqq2gZ2ll::scale() const {
return sHat();
}
double MEqq2gZ2ll::me2() const {
Energy2 m2 = meMomenta()[2].mass2();
// Energy2 p1p3 = 0.5*(m2 - tHat());
// Energy2 p1p2 = 0.5*sHat();
Energy2 p1p3 = meMomenta()[0].dot(meMomenta()[2]);
Energy2 p1p2 = meMomenta()[0].dot(meMomenta()[1]);
Energy4 pt2 = sqr(p1p3);
Energy4 pts = p1p3*p1p2;
Energy4 ps2 = sqr(p1p2);
Energy4 psm = p1p2*m2;
int up = abs(mePartonData()[0]->id() + 1)%2;
lastCont =
(coefs[0 + up]*(pt2 - pts) + coefs[2 + up]*(ps2 + psm))/sqr(sHat());
double intr = 0.25*(coefs[4 + up]*pt2 + coefs[6 + up]*pts +
coefs[8 + up]*ps2 + coefs[10 + up]*psm)*
(sHat() - mZ2)/(sHat()*(sqr(sHat() - mZ2) + mZ2*GZ2));
lastBW = 0.25*(coefs[12 + up]*pt2 + coefs[14 + up]*pts +
coefs[16 + up]*ps2 + coefs[18 + up]*psm)/
(sqr(sHat() - mZ2) + mZ2*GZ2);
double alphaS = SM().alphaS(scale());
//int Nf = SM().Nf(scale());
int Nf = 3;
DVector save;
meInfo(save << lastCont << lastBW);
return (lastCont + intr + lastBW)*sqr(SM().alphaEM(scale()))*
- (1.0 + alphaS/pi + (1.986-0.115*Nf)*sqr(alphaS/pi));
+ (1.0 + alphaS/Constants::pi + (1.986-0.115*Nf)*sqr(alphaS/Constants::pi));
}
Selector<MEqq2gZ2ll::DiagramIndex>
MEqq2gZ2ll::diagrams(const DiagramVector & diags) const {
if ( lastXCombPtr() ) {
lastCont = meInfo()[0];
lastBW = meInfo()[1];
}
Selector<DiagramIndex> sel;
for ( DiagramIndex i = 0; i < diags.size(); ++i ) {
if ( diags[i]->id() == -1 ) sel.insert(lastCont, i);
else if ( diags[i]->id() == -2 ) sel.insert(lastBW, i);
}
return sel;
}
Selector<const ColourLines *>
MEqq2gZ2ll::colourGeometries(tcDiagPtr diag) const {
static ColourLines c("1 -2");
Selector<const ColourLines *> sel;
sel.insert(1.0, &c);
return sel;
}
IBPtr MEqq2gZ2ll::clone() const {
return new_ptr(*this);
}
IBPtr MEqq2gZ2ll::fullclone() const {
return new_ptr(*this);
}
void MEqq2gZ2ll::doinit() throw(InitException) {
double C = sqr(4.0*Constants::pi)/3.0;
double SW2 = SM().sin2ThetaW();
double SW4 = sqr(SW2);
double SW6 = SW2*SW4;
double SW8 = SW2*SW6;
double CW2 = 1.0 - SW2;
coefs[0] = 16.0*C;
coefs[1] = 64.0*C;
coefs[2] = 8.0*C;
coefs[3] = 32.0*C;
C /= (CW2*SW2);
coefs[4] = 4.0*(32.0*SW4 - 32.0*SW2 + 6.0)*C;
coefs[5] = 8.0*(64.0*SW4 - 40.0*SW2 + 6.0)*C;
coefs[6] = -4.0*(32.0*SW4 - 32.0*SW2 + 12.0)*C;
coefs[7] = -8.0*(64.0*SW4 - 40.0*SW2 + 12.0)*C;
coefs[8] = 4.0*(16.0*SW4 - 16.0*SW2 + 6.0)*C;
coefs[9] = 8.0*(32.0*SW4 - 20.0*SW2 + 6.0)*C;
coefs[10] = 4.0*(16.0*SW4 - 16.0*SW2 + 3.0)*C;
coefs[11] = 8.0*(32.0*SW4 - 20.0*SW2 + 3.0)*C;
C /= (CW2*SW2);
coefs[12] = ( 64.0*SW8 - 128.0*SW6 + 128.0*SW4 - 48.0*SW2 + 9.0)*C;
coefs[13] = (256.0*SW8 - 320.0*SW6 + 200.0*SW4 - 60.0*SW2 + 9.0)*C;
coefs[14] = -( 64.0*SW8 - 128.0*SW6 + 176.0*SW4 - 96.0*SW2 + 18.0)*C;
coefs[15] = -(256.0*SW8 - 320.0*SW6 + 296.0*SW4 - 120.0*SW2 + 18.0)*C;
coefs[16] = ( 32.0*SW8 - 64.0*SW6 + 88.0*SW4 - 48.0*SW2 + 9.0)*C;
coefs[17] = (128.0*SW8 - 160.0*SW6 + 148.0*SW4 - 60.0*SW2 + 9.0)*C;
coefs[18] = ( 32.0*SW8 - 64.0*SW6 + 28.0*SW4 - 6.0*SW2)*C;
coefs[19] = (128.0*SW8 - 160.0*SW6 + 64.0*SW4 - 12.0*SW2)*C;
tcPDPtr Z0 = getParticleData(ParticleID::Z0);
mZ2 = sqr(Z0->mass());
GZ2 = sqr(Z0->width());
ME2to2QCD::doinit();
}
void MEqq2gZ2ll::persistentOutput(PersistentOStream & os) const {
os << coefs << ounit(mZ2, GeV2) << ounit(GZ2, GeV2) << lastCont << lastBW;
}
void MEqq2gZ2ll::persistentInput(PersistentIStream & is, int) {
is >> coefs >> iunit(mZ2, GeV2) >> iunit(GZ2, GeV2) >> lastCont >> lastBW;
}
ClassDescription<MEqq2gZ2ll> MEqq2gZ2ll::initMEqq2gZ2ll;
void MEqq2gZ2ll::Init() {
static ClassDocumentation<MEqq2gZ2ll> documentation
("The ThePEG::MEqq2gZ2ll class implements the full"
"\\f$q\\bar{q}\\rightarrow\\gamma/Z^0\\rightarrow\\ell^+\\ell^-\\f$ "
"matrix element including the interference terms.");
}

File Metadata

Mime Type
text/x-diff
Expires
Thu, Apr 24, 6:35 AM (1 d, 15 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4887758
Default Alt Text
(15 KB)

Event Timeline