herwig is hosted by Hepforge, IPPP Durham
Herwig  7.1.5
DecayPhaseSpaceMode.h
1 // -*- C++ -*-
2 //
3 // DecayPhaseSpaceMode.h is a part of Herwig - A multi-purpose Monte Carlo event generator
4 // Copyright (C) 2002-2017 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_DecayPhaseSpaceMode_H
10 #define HERWIG_DecayPhaseSpaceMode_H
11 //
12 // This is the declaration of the DecayPhaseSpaceMode class.
13 //
14 #include "ThePEG/Interface/Interfaced.h"
15 #include "DecayPhaseSpaceMode.fh"
16 #include "DecayPhaseSpaceChannel.h"
17 #include "Herwig/PDT/GenericWidthGenerator.h"
18 #include "Herwig/PDT/GenericMassGenerator.h"
19 #include <Herwig/Decay/DecayVertex.h>
20 #include "DecayIntegrator.h"
21 
22 namespace Herwig {
23 using namespace ThePEG;
24 
44 
48  friend ostream & operator<<(ostream &, const DecayPhaseSpaceMode &);
49 
54  friend class DecayIntegrator;
55 
60  friend class DecayPhaseSpaceChannel;
61 
62 public:
63 
69  DecayPhaseSpaceMode() : _maxweight(0.),_niter(10), _npoint(10000), _ntry(500),
70  _partial(-1), _testOnShell(false), _ichannel(999) {}
71 
80  DecayPhaseSpaceMode(tPDVector in, tcDecayIntegratorPtr intin,bool onShell=false)
81  : _integrator(intin), _maxweight(0.),
82  _niter(10), _npoint(10000), _ntry(500),
83  _extpart(in), _partial(-1), _testOnShell(onShell), _ichannel(999) {}
85 
91  tcPDPtr externalParticles(int ix) const {return _extpart[ix];}
92 
97  unsigned int numberofParticles() const {return _extpart.size();}
98 
103  unsigned int numberChannels() const {return _channels.size();}
104 
109  void addChannel(DecayPhaseSpaceChannelPtr channel) {
110  channel->init();
111  _channels.push_back(channel);
112  }
113 
122  void resetIntermediate(int ichan, tcPDPtr part, Energy mass, Energy width) {
123  if(!part) return;
124  _channels[ichan]->resetIntermediate(part,mass,width);
125  }
126 
134  void resetIntermediate(tcPDPtr part, Energy mass, Energy width) {
135  for(unsigned int ix=0,N=_channels.size();ix<N;++ix)
136  _channels[ix]->resetIntermediate(part,mass,width);
137  }
138 
143  double maxWeight() const {return _maxweight;}
144 
149  void setMaxWeight(double wgt) const {_maxweight=wgt;}
150 
156  double channelWeight(unsigned int ichan) const {return _channelwgts[ichan];}
157 
162  void setWeights(const vector<double> & in) const {_channelwgts=in;}
163 
167  unsigned int selectedChannel() const {return _ichannel;}
168 
172  bool testOnShell() const { return _testOnShell; }
173 
174 protected:
175 
182  Energy initializePhaseSpace(bool init, bool onShell=false);
183 
190  void setIntegrate(int iter,int points,int ntry) {
191  _niter=iter;
192  _npoint=points;
193  _ntry=ntry;
194  }
195 
201  void setPartialWidth(int in) {_partial=in;}
203 
204 
207 
216  double me2(const int ichan ,const Particle & inpart,
217  const ParticleVector &outpart,DecayIntegrator::MEOption opt) const {
218  return _integrator->me2(ichan,inpart,outpart,opt);
219  }
220 
230  ParticleVector generate(bool intermediates,bool cc,const Particle & inpart) const;
231 
232 
238  int selectChannel(const Particle & inpart, ParticleVector & outpart) const {
239  // if using flat phase-space don't need to do this
240  if(_channelwgts.empty()) return 0;
241  vector<double> mewgts(_channels.size(),0.0);
242  double total=0.;
243  for(unsigned int ix=0,N=_channels.size();ix<N;++ix) {
244  mewgts[ix]=me2(ix,inpart,outpart,DecayIntegrator::Calculate);
245  total+=mewgts[ix];
246  }
247  // randomly pick a channel
248  total*=UseRandom::rnd();
249  int ichan=-1;
250  do {
251  ++ichan;
252  total-=mewgts[ichan];
253  }
254  while(ichan<int(_channels.size())&&total>0.);
255  return ichan;
256  }
257 
268  Energy weight(bool cc,int & ichan, const Particle & in,
269  ParticleVector & particles,bool first,
270  bool onShell=false) const {
271  ichan=0;
272  Energy phwgt = (_channels.size()==0) ?
273  flatPhaseSpace(cc,in,particles,onShell) : channelPhaseSpace(cc,ichan,in,particles,onShell);
274  // generate the matrix element
275  return me2(-1,in,particles,
276  first ? DecayIntegrator::Initialize : DecayIntegrator::Calculate)*phwgt;
277  }
278 
287  Energy flatPhaseSpace(bool cc,const Particle & inpart, ParticleVector & outpart,
288  bool onShell=false) const;
289 
299  Energy channelPhaseSpace(bool cc,int & ichan, const Particle & in,
300  ParticleVector & particles,
301  bool onShell=false) const;
302 
308  void constructVertex(const Particle & in, const ParticleVector & out) const;
309 
316  vector<Energy> externalMasses(Energy inmass,double & wgt, bool onShell) const;
318 
319 public:
320 
327  void persistentOutput(PersistentOStream & os) const;
328 
334  void persistentInput(PersistentIStream & is, int version);
336 
340  static void Init();
341 
342 protected:
343 
350  virtual IBPtr clone() const {return new_ptr(*this);}
351 
356  virtual IBPtr fullclone() const {return new_ptr(*this);}
358 
359 protected:
360 
368  virtual void doinit();
369 
373  virtual void doinitrun();
375 
376 private:
377 
382 
386  DecayPhaseSpaceMode & operator=(const DecayPhaseSpaceMode &) = delete;
387 
388  private:
389 
393  tcDecayIntegratorPtr _integrator;
394 
398  vector<DecayPhaseSpaceChannelPtr> _channels;
399 
403  mutable vector<double> _channelwgts;
404 
408  mutable double _maxweight;
409 
413  int _niter;
414 
418  int _npoint;
419 
423  int _ntry;
424 
429 
433  int _partial;
434 
438  cGenericWidthGeneratorPtr _widthgen;
439 
443  vector<cGenericMassGeneratorPtr> _massgen;
444 
450 
454  mutable unsigned int _ichannel;
455 
456 };
457 
461 ostream & operator<<(ostream &, const DecayPhaseSpaceMode &);
462 
463 }
464 
465 
466 #include "ThePEG/Utilities/ClassTraits.h"
467 
468 namespace ThePEG {
469 
472 template <>
479  typedef Interfaced NthBase;
480 };
481 
482 template <>
488  : public ClassTraitsBase<Herwig::DecayPhaseSpaceMode> {
490  static string className() { return "Herwig::DecayPhaseSpaceMode"; }
491 };
492 
495 }
496 
497 #endif /* HERWIG_DecayPhaseSpaceMode_H */
DecayPhaseSpaceMode(tPDVector in, tcDecayIntegratorPtr intin, bool onShell=false)
Constructor with a pointer to a DecayPhaseIntegrator and a vector of particle data objects for the ex...
int _npoint
Number of weights for each iteration of the initialization.
int _partial
Which of the partial widths of the incoming particle to use.
bool testOnShell() const
test on/off-shell kinematics
unsigned int selectedChannel() const
Access to the selected channel.
cGenericWidthGeneratorPtr _widthgen
The width generator for the incoming particle.
The DecayPhaseSpaceMode class is designed to store a group of phase-space channels for use by the Dec...
void setIntegrate(int iter, int points, int ntry)
Set the integration parameters.
ThePEG::Ptr< InterfacedBase >::pointer IBPtr
tcPDPtr externalParticles(int ix) const
Access to the external particles.
double channelWeight(unsigned int ichan) const
Get the weight for a channel.
virtual IBPtr fullclone() const
Make a clone of this object, possibly modifying the cloned object to make it sane.
void setWeights(const vector< double > &in) const
Set the weights for the different channels.
tPDVector _extpart
External particles.
static double rnd()
int _ntry
Number of attempts to generate the decay.
vector< double > _channelwgts
the weights for the different channels
int selectChannel(const Particle &inpart, ParticleVector &outpart) const
Select which channel to use to output the particles.
tcDecayIntegratorPtr _integrator
pointer to the decayer
DecayPhaseSpaceMode()
Default constructor.
MEOption
Enum for the matrix element option.
ostream & operator<<(ostream &, const DecayIntegrator &)
Output information on the DecayIntegrator for debugging purposes.
void resetIntermediate(int ichan, tcPDPtr part, Energy mass, Energy width)
Reset the properties of one of the intermediate particles.
vector< tPDPtr > tPDVector
double me2(const int ichan, const Particle &inpart, const ParticleVector &outpart, DecayIntegrator::MEOption opt) const
Access to the matrix element from the decayer.
void setPartialWidth(int in)
Set the partial width to use for normalization.
Main class for Decayers implementing multi-channel phase space integration.
ThePEG::Ptr< ParticleData >::transient_const_pointer tcPDPtr
void resetIntermediate(tcPDPtr part, Energy mass, Energy width)
Reset the properties of one of the intermediate particles.
double _maxweight
the maximum weight for the decay
vector< DecayPhaseSpaceChannelPtr > _channels
pointers to the phase-space channels
Energy weight(bool cc, int &ichan, const Particle &in, ParticleVector &particles, bool first, bool onShell=false) const
Return the weight for a given phase-space point.
static ClassDescription< DecayPhaseSpaceMode > initDecayPhaseSpaceMode
Describe a concrete class with persistent data.
pair< double, double > generate(const Generator< Density > &gen, double r)
Generate a random variable and return its weight.
void addChannel(DecayPhaseSpaceChannelPtr channel)
Add a new channel.
This class is designed to store the information needed for a given phase-space channel for use by the...
bool _testOnShell
Whether to check on-shell or off-shell kinematics in doinit, if on-shell off-shell is tested in initr...
-*- C++ -*-
double maxWeight() const
Get the maximum weight for the decay.
int _niter
Number of iterations for the initialization.
virtual IBPtr clone() const
Make a simple clone of this object.
vector< PPtr > ParticleVector
vector< cGenericMassGeneratorPtr > _massgen
The mass generators for the outgoing particles.
unsigned int numberofParticles() const
Number of external particles.
unsigned int _ichannel
The selected channel.
unsigned int numberChannels() const
Number of channels.
void setMaxWeight(double wgt) const
Set the maximum weight for the decay.