herwig is hosted by Hepforge, IPPP Durham
Herwig 7.3.0
PhaseSpaceMode.h
1// -*- C++ -*-
2//
3// PhaseSpaceMode.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_PhaseSpaceMode_H
10#define Herwig_PhaseSpaceMode_H
11//
12// This is the declaration of the PhaseSpaceMode class.
13//
14
16#include "PhaseSpaceMode.fh"
17#include "PhaseSpaceChannel.h"
18#include "Herwig/PDT/GenericWidthGenerator.h"
19#include "Herwig/PDT/GenericMassGenerator.h"
20
21namespace Herwig {
22using namespace ThePEG;
23
42class PhaseSpaceMode: public Base {
43
44 friend class PhaseSpaceChannel;
45
46public:
47
54 testOnShell_(false), eMax_(ZERO),
55 eps_(ZERO), nRand_(0)
56 {};
57
62 double wMax, tPDPtr in2=tPDPtr(),
63 Energy eMax=ZERO) : incoming_(make_pair(in1,in2)),
64 maxWeight_(wMax),
65 outgoing_(out), partial_(-1),
66 testOnShell_(false), eMax_(eMax),
67 eps_(ZERO), nRand_(0)
68 {};
70
71public:
72
83 tcDecayIntegratorPtr decayer,
84 bool intermediates,bool cc);
85
91 channel.init(this);
92 channels_.push_back(channel);
93 }
94
103 void resetIntermediate(int ichan, tcPDPtr part, Energy mass, Energy width) {
104 if(!part) return;
105 channels_[ichan].resetIntermediate(part,mass,width);
106 }
107
115 void resetIntermediate(tcPDPtr part, Energy mass, Energy width) {
116 for(PhaseSpaceChannel & channel : channels_)
117 channel.resetIntermediate(part,mass,width);
118 }
119
123 const vector<PhaseSpaceChannel> & channels() const {return channels_;}
124
128 void setWeights(const vector<double> & wgts) {
129 assert(wgts.size()==channels_.size());
130 for(unsigned int ix=0;ix<channels_.size();++ix)
131 channels_[ix].weight(wgts[ix]);
132 }
133
137 unsigned int selectedChannel() const {return iChannel_;}
138
142 unsigned int nRand() const {return nRand_;}
143
147 void checkOnShell(bool in) {testOnShell_=in;}
148
149public:
150
158
164 void persistentInput(PersistentIStream & is, int version);
166
173 static void Init();
174
175public :
176
180 void init();
181
185 void initrun();
186
191 double maxWeight() const {return maxWeight_;}
192
197 void maxWeight(double wgt) const {maxWeight_=wgt;}
198
203 Energy initializePhaseSpace(bool init, tcDecayIntegratorPtr decayer,
204 bool onShell=false);
205
209 pair<PDPtr,PDPtr> incoming() const {return incoming_;}
210
215 tPDVector outgoing() const {return outgoing_;}
216
221 tPDVector outgoingCC() const {return outgoingCC_;}
222
227 unsigned int numberOfParticles() const {return outgoing_.size();}
228
234 void setPartialWidth(int in) {partial_=in;}
235
239 Energy epsilonPS() const {return eps_;}
240
244 void fillStack(const double * r) {
245 assert(rStack_.empty());
246 for(unsigned int ix=nRand_;ix>0;--ix)
247 rStack_.push(r[nRand_-1]);
248 }
249
253 void fillStack() {
254 assert(rStack_.empty());
255 for(unsigned int ix=0;ix<nRand_;++ix) rStack_.push(UseRandom::rnd());
256 }
257
265 Energy weight(int & ichan, const Lorentz5Momentum & in,
266 vector<Lorentz5Momentum> & momenta,
267 bool onShell=false) const {
268 ichan=0;
269 // flat phase-space
270 if(channels_.empty())
271 return flatPhaseSpace(in,momenta,onShell);
272 // multi-channel
273 else
274 return channelPhaseSpace(ichan,in,momenta,onShell);
275 }
276
277public :
278
282 friend ostream & operator<<(ostream & os, const PhaseSpaceMode & mode) {
283 os << "The mode has " << mode.channels_.size() << " channels\n";
284 if(mode.incoming_.second==PDPtr())
285 os << "This is a mode for the decay of " << mode.incoming_.first->PDGName() << " to ";
286 else
287 os << "This is a mode for " << mode.incoming_.first->PDGName() << ", "
288 << mode.incoming_.second->PDGName() << " to ";
289 for(tPDPtr out : mode.outgoing_) os << out->PDGName() << " ";
290 os << "\n";
291 for(const PhaseSpaceChannel & channel : mode.channels_) os << channel;
292 return os;
293 }
294
295private:
296
304 Energy flatPhaseSpace(const Lorentz5Momentum & in,
305 vector<Lorentz5Momentum> & momenta,
306 bool onShell=false) const;
307
316 Energy channelPhaseSpace(int & ichan, const Lorentz5Momentum & in,
317 vector<Lorentz5Momentum> & momenta,
318 bool onShell=false) const;
319
326 vector<Energy> externalMasses(Energy inmass,double & wgt, bool onShell) const;
327
333 void constructVertex(const Particle & in, const ParticleVector & out,
334 tcDecayIntegratorPtr decayer) const;
335
336private:
337
343
344private:
345
349 pair<PDPtr,PDPtr> incoming_;
350
354 vector<PhaseSpaceChannel> channels_;
355
359 mutable double maxWeight_;
360
365 tPDVector outgoingCC_;
366
371
375 cGenericWidthGeneratorPtr widthGen_;
376
380 vector<cGenericMassGeneratorPtr> massGen_;
381 vector<double> BRsum_;
387
391 Energy eMax_;
392
396 mutable unsigned int iChannel_;
397
401 Energy eps_;
402
406 unsigned int nRand_;
407
411 mutable stack<double> rStack_;
412};
413
414}
415
416#endif /* Herwig_PhaseSpaceMode_H */
This class is designed to store the information needed for a given phase-space channel for use by the...
void init(tPhaseSpaceModePtr mode)
Initialize the channel.
The PhaseSpaceMode class is designed to store a group of phase-space channels for use by the DecayInt...
PhaseSpaceMode & operator=(const PhaseSpaceMode &)=delete
The assignment operator is private and must never be called.
unsigned int nRand_
Number of random numbers required.
void persistentOutput(PersistentOStream &os) const
Function used to write out object persistently.
bool testOnShell_
Whether to check on-shell or off-shell kinematics in doinit, if on-shell off-shell is tested in initr...
pair< PDPtr, PDPtr > incoming() const
The incoming particles.
tPDVector outgoing() const
Access to the outging particles.
int partial_
Which of the partial widths of the incoming particle to use.
void persistentInput(PersistentIStream &is, int version)
Function used to read in object persistently.
void initrun()
Initialisation before the run stage.
double maxWeight_
The maximum weight.
void setPartialWidth(int in)
Set the partial width to use for normalization.
Energy initializePhaseSpace(bool init, tcDecayIntegratorPtr decayer, bool onShell=false)
Initialise the phase space.
Energy weight(int &ichan, const Lorentz5Momentum &in, vector< Lorentz5Momentum > &momenta, bool onShell=false) const
Return the weight for a given phase-space point.
PhaseSpaceMode(tPDPtr in1, tPDVector out, double wMax, tPDPtr in2=tPDPtr(), Energy eMax=ZERO)
The default constructor.
vector< cGenericMassGeneratorPtr > massGen_
The mass generators for the outgoing particles.
stack< double > rStack_
Stack for the random numbers.
void checkOnShell(bool in)
Set whether of not decays are generated on-shell.
tPDVector outgoing_
The external particles.
friend ostream & operator<<(ostream &os, const PhaseSpaceMode &mode)
A friend operator to allow the mode to be outputted for debugging purposes.
unsigned int numberOfParticles() const
Number of outgoing particles.
static void Init()
The standard Init function used to initialize the interfaces.
void init()
Initialize the phase space.
double maxWeight() const
Get the maximum weight for the decay.
Energy flatPhaseSpace(const Lorentz5Momentum &in, vector< Lorentz5Momentum > &momenta, bool onShell=false) const
Return the weight and momenta for a flat phase-space decay.
void resetIntermediate(tcPDPtr part, Energy mass, Energy width)
Reset the properties of one of the intermediate particles.
unsigned int iChannel_
The selected channel.
ParticleVector generateDecay(const Particle &inpart, tcDecayIntegratorPtr decayer, bool intermediates, bool cc)
Generate the decay.
vector< Energy > externalMasses(Energy inmass, double &wgt, bool onShell) const
Generate the masses of the external particles.
Energy epsilonPS() const
Access to the epsilon parameter.
Energy channelPhaseSpace(int &ichan, const Lorentz5Momentum &in, vector< Lorentz5Momentum > &momenta, bool onShell=false) const
Generate a phase-space point using multichannel phase space.
void resetIntermediate(int ichan, tcPDPtr part, Energy mass, Energy width)
Reset the properties of one of the intermediate particles.
void setWeights(const vector< double > &wgts)
Set the weights.
Energy eMax_
The maximum energy for the mode.
cGenericWidthGeneratorPtr widthGen_
The width generator for the incoming particle.
void addChannel(PhaseSpaceChannel channel)
Add a new channel.
unsigned int selectedChannel() const
Access to the selected channel.
pair< PDPtr, PDPtr > incoming_
The incoming particles.
void fillStack(const double *r)
Fill the stack.
PhaseSpaceMode()
The default constructor.
tPDVector outgoingCC() const
Access to the outging particles.
const vector< PhaseSpaceChannel > & channels() const
The phase-space channels.
void fillStack()
Fill the stack.
void constructVertex(const Particle &in, const ParticleVector &out, tcDecayIntegratorPtr decayer) const
Construct the vertex for spin corrections.
unsigned int nRand() const
Number of rnadom numbers needed.
void maxWeight(double wgt) const
Set the maximum weight for the decay.
vector< PhaseSpaceChannel > channels_
The phase-space channels.
static double rnd()
-*- C++ -*-
ThePEG::Ptr< ParticleData >::transient_pointer tPDPtr
vector< tPDPtr > tPDVector
vector< PPtr > ParticleVector
constexpr ZeroUnit ZERO
ThePEG::Ptr< ParticleData >::transient_const_pointer tcPDPtr
ThePEG::Ptr< ParticleData >::pointer PDPtr