herwig is hosted by Hepforge, IPPP Durham
Herwig 7.3.0
PhaseSpaceChannel.h
1// -*- C++ -*-
2//
3// PhaseSpaceChannel.h is a part of Herwig - A multi-purpose Monte Carlo event generator
4// Copyright (C) 2002-2019 The Herwig Collaboration
5//
6// Herwig is licenced under version 3 of the GPL, see COPYING for details.
7// Please respect the MCnet academic guidelines, see GUIDELINES for details.
8//
9#ifndef Herwig_PhaseSpaceChannel_H
10#define Herwig_PhaseSpaceChannel_H
11//
12// This is the declaration of the PhaseSpaceChannel class.
13//
14
16#include "ThePEG/PDT/ParticleData.h"
17#include "ThePEG/Persistency/PersistentOStream.h"
18#include "ThePEG/Persistency/PersistentIStream.h"
19#include "ThePEG/Utilities/EnumIO.h"
20#include "PhaseSpaceMode.fh"
21#include "ThePEG/MatrixElement/Tree2toNDiagram.h"
22
23namespace Herwig {
24
25using namespace ThePEG;
26
53
54public:
55
56 friend PersistentOStream;
57 friend PersistentIStream;
58
63
67 enum Jacobian {BreitWigner,Power,OnShell};
68
76PhaseSpaceResonance(cPDPtr part) :particle(part), mass2(sqr(part->mass())), mWidth(part->mass()*part->width()),
77 jacobian(BreitWigner), power(0.), children(make_pair(0,0))
78{};
83
87 Energy2 mass2;
88
92 Energy2 mWidth;
93
98
102 double power;
103
107 pair<int,int> children;
108
112 vector<int> descendents;
113
114};
115
116public:
117
122
126 PhaseSpaceChannel(tPhaseSpaceModePtr inm, bool skip=false);
127
133 if(intermediates_.size()==1&&skipFirst_) {
134 skipFirst_=false;
135 }
136 else
137 intermediates_.push_back(PhaseSpaceResonance(res));
138 if(iAdd_<0) return *this;
139 if(intermediates_[iAdd_].children.first==0)
140 intermediates_[iAdd_].children.first = 1-int(intermediates_.size());
141 else
142 intermediates_[iAdd_].children.second = 1-int(intermediates_.size());
143 iAdd_=-1;
144 return *this;
145 }
146
152 if(iAdd_<0) iAdd_ = o;
153 else if(o>=0) {
154 if(intermediates_[iAdd_].children.first==0)
155 intermediates_[iAdd_].children.first = o;
156 else
157 intermediates_[iAdd_].children.second = o;
158 iAdd_=-1;
159 }
160 else if(o<0) {
161 assert(false);
162 }
163 return *this;
164 }
165
169 void setJacobian(unsigned int ires, PhaseSpaceResonance::Jacobian jac, double power) {
170 intermediates_[ires].jacobian = jac;
171 intermediates_[ires].power = power;
172 }
173
174public:
175
179 void init(tPhaseSpaceModePtr mode);
180
184 void initrun(tPhaseSpaceModePtr mode);
185
190
194 const double & weight() const {return weight_;}
195
199 void weight(double in) {weight_=in;}
200
210 void resetIntermediate(tcPDPtr part,Energy mass,Energy width) {
211 if(!part) return;
213 if(res.particle!=part) continue;
214 res.mass2 = sqr(mass);
215 res.mWidth = mass*width;
216 }
217 }
218
227 vector<Lorentz5Momentum> generateMomenta(const Lorentz5Momentum & pin,
228 const vector<Energy> & massext) const;
229
230
238 double generateWeight(const vector<Lorentz5Momentum> & output) const;
239
251 void generateIntermediates(bool cc,const Particle & in, ParticleVector & out);
252
256 ThePEG::Ptr<ThePEG::Tree2toNDiagram>::pointer createDiagram() const;
257
258public:
259
266 const PhaseSpaceChannel & x) {
267 os << x.weight_ << x.initialized_ << x.intermediates_;
268 return os;
269 }
270
277 PhaseSpaceChannel & x) {
278 is >> x.weight_ >> x.initialized_ >> x.intermediates_;
279 return is;
280 }
281
286 friend ostream & operator<<(ostream & os, const PhaseSpaceChannel & channel);
287
288private:
289
294 vector<int> & children) {
295 if(res.children.first>0)
296 children.push_back(res.children.first);
297 else
298 findChildren(intermediates_[abs(res.children.first)],children);
299 if(!res.particle) return;
300 if(res.children.second>0)
301 children.push_back(res.children.second);
302 else
303 findChildren(intermediates_[abs(res.children.second)],children);
304 }
305
315 void twoBodyDecay(const Lorentz5Momentum & p,
316 const Energy m1, const Energy m2,
317 Lorentz5Momentum & p1, Lorentz5Momentum & p2) const;
318
328 Energy lower,Energy upper,
329 const double & rnd) const;
330
338 InvEnergy2 massWeight(const PhaseSpaceResonance & res,
339 Energy moff,Energy lower,Energy upper) const;
341
347 double atanhelper(const PhaseSpaceResonance & res, Energy limit) const;
348
349private:
350
354 tPhaseSpaceModePtr mode_;
355
359 vector<PhaseSpaceResonance> intermediates_;
360
364 int iAdd_ = -1;
365
369 double weight_;
370
375
380
381};
382
383
391 os << x.particle << ounit(x.mass2,GeV2) << ounit(x.mWidth,GeV2)
392 << oenum(x.jacobian) << x.power << x.children << x.descendents;
393 return os;
394}
395
403 is >> x.particle >> iunit(x.mass2,GeV2) >> iunit(x.mWidth,GeV2)
404 >> ienum(x.jacobian) >> x.power >> x.children >> x.descendents;
405 return is;
406}
407
412ostream & operator<<(ostream & os, const PhaseSpaceChannel & channel);
413
417class PhaseSpaceError: public Exception {};
418
419}
420
421#endif /* Herwig_PhaseSpaceChannel_H */
This class is designed to store the information needed for a given phase-space channel for use by the...
vector< Lorentz5Momentum > generateMomenta(const Lorentz5Momentum &pin, const vector< Energy > &massext) const
Generate the momenta of the external particles.
vector< PhaseSpaceResonance > intermediates_
The intermediates.
bool skipFirst_
Wheter or not to skiip the first resonance.
int iAdd_
Integer to keep track of what we are adding.
const double & weight() const
The weight.
void init(tPhaseSpaceModePtr mode)
Initialize the channel.
bool checkKinematics()
Check the kinematics.
void resetIntermediate(tcPDPtr part, Energy mass, Energy width)
Reset the properties of an intermediate particle.
PhaseSpaceChannel(tPhaseSpaceModePtr inm, bool skip=false)
Constructor with incoming particles.
ThePEG::Ptr< ThePEG::Tree2toNDiagram >::pointer createDiagram() const
Create a diagrams for the channel.
PhaseSpaceChannel & operator,(tPDPtr res)
If less than zero indicate that this channel is competed.
double atanhelper(const PhaseSpaceResonance &res, Energy limit) const
Helper function for the weight calculation.
Energy generateMass(const PhaseSpaceResonance &res, Energy lower, Energy upper, const double &rnd) const
Generate the mass of a resonance.
void initrun(tPhaseSpaceModePtr mode)
Initialize the channel.
double generateWeight(const vector< Lorentz5Momentum > &output) const
Generate the weight for this channel given a phase space configuration.
void generateIntermediates(bool cc, const Particle &in, ParticleVector &out)
Generate the final-state particles including the intermediate resonances.
bool initialized_
Whether or not its been initialized.
tPhaseSpaceModePtr mode_
Pointer to the phase-space mode.
void twoBodyDecay(const Lorentz5Momentum &p, const Energy m1, const Energy m2, Lorentz5Momentum &p1, Lorentz5Momentum &p2) const
Calculate the momenta for a two body decay The return value indicates success or failure.
friend PersistentIStream & operator>>(PersistentIStream &is, PhaseSpaceChannel &x)
Input operator to allow persistently written data to be read in.
void weight(double in)
Set the weight.
InvEnergy2 massWeight(const PhaseSpaceResonance &res, Energy moff, Energy lower, Energy upper) const
Return the weight for a given resonance.
friend PersistentOStream & operator<<(PersistentOStream &os, const PhaseSpaceChannel &x)
Output operator to allow the structure to be persistently written.
PhaseSpaceChannel()
Default constructior.
friend ostream & operator<<(ostream &os, const PhaseSpaceChannel &channel)
A friend output operator to allow the channel to be outputted for debugging purposes.
void setJacobian(unsigned int ires, PhaseSpaceResonance::Jacobian jac, double power)
Set the jacobian for a given resonance.
void findChildren(const PhaseSpaceResonance &res, vector< int > &children)
Find the external particles which are the children of a given resonance.
exception for this class and those inheriting from it
-*- C++ -*-
ThePEG::Ptr< ParticleData >::transient_pointer tPDPtr
IUnit< T, UT > iunit(T &t, const UT &ut)
PersistentIStream & operator>>(PersistentIStream &is, HandlerGroup< HDLR > &hg)
IEnum< T > ienum(T &t)
OEnum< T > oenum(const T &t)
vector< PPtr > ParticleVector
ThePEG::Ptr< ParticleData >::transient_const_pointer tcPDPtr
ostream & operator<<(ostream &, const Collision &)
constexpr auto sqr(const T &x) -> decltype(x *x)
ThePEG::Ptr< ParticleData >::const_pointer cPDPtr
OUnit< T, UT > ounit(const T &t, const UT &ut)
Struct for the intermediates in the phase-space channel.
PhaseSpaceResonance(cPDPtr part)
The constructor.
vector< int > descendents
The final descendents.