herwig is hosted by Hepforge, IPPP Durham
Herwig  7.2.1
DecayPhaseSpaceMode.h
1 // -*- C++ -*-
2 //
3 // DecayPhaseSpaceMode.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_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 
43 class DecayPhaseSpaceMode: public Base {
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 
177  Energy epsilonPS() const {return _eps;}
178 
179 protected:
180 
187  Energy initializePhaseSpace(bool init, bool onShell=false);
188 
195  void setIntegrate(int iter,int points,int ntry) {
196  _niter=iter;
197  _npoint=points;
198  _ntry=ntry;
199  }
200 
206  void setPartialWidth(int in) {_partial=in;}
208 
209 
212 
221  double me2(const int ichan ,const Particle & inpart,
222  const ParticleVector &outpart,DecayIntegrator::MEOption opt) const {
223  return _integrator->me2(ichan,inpart,outpart,opt);
224  }
225 
235  ParticleVector generate(bool intermediates,bool cc,const Particle & inpart) const;
236 
237 
243  int selectChannel(const Particle & inpart, ParticleVector & outpart) const {
244  // if using flat phase-space don't need to do this
245  if(_channelwgts.empty()) return 0;
246  vector<double> mewgts(_channels.size(),0.0);
247  double total=0.;
248  for(unsigned int ix=0,N=_channels.size();ix<N;++ix) {
249  mewgts[ix]=me2(ix,inpart,outpart,DecayIntegrator::Calculate);
250  total+=mewgts[ix];
251  }
252  // randomly pick a channel
253  total*=UseRandom::rnd();
254  int ichan=-1;
255  do {
256  ++ichan;
257  total-=mewgts[ichan];
258  }
259  while(ichan<int(_channels.size())&&total>0.);
260  return ichan;
261  }
262 
273  Energy weight(bool cc,int & ichan, const Particle & in,
274  ParticleVector & particles,bool first,
275  bool onShell=false) const {
276  ichan=0;
277  Energy phwgt = (_channels.size()==0) ?
278  flatPhaseSpace(cc,in,particles,onShell) : channelPhaseSpace(cc,ichan,in,particles,onShell);
279  // generate the matrix element
280  return me2(-1,in,particles,
281  first ? DecayIntegrator::Initialize : DecayIntegrator::Calculate)*phwgt;
282  }
283 
292  Energy flatPhaseSpace(bool cc,const Particle & inpart, ParticleVector & outpart,
293  bool onShell=false) const;
294 
304  Energy channelPhaseSpace(bool cc,int & ichan, const Particle & in,
305  ParticleVector & particles,
306  bool onShell=false) const;
307 
313  void constructVertex(const Particle & in, const ParticleVector & out) const;
314 
321  vector<Energy> externalMasses(Energy inmass,double & wgt, bool onShell) const;
323 
324 public:
325 
332  void persistentOutput(PersistentOStream & os) const;
333 
339  void persistentInput(PersistentIStream & is, int version);
341 
345  static void Init();
346 
347 protected:
348 
356  virtual void init();
357 
361  virtual void initrun();
363 
364 private:
365 
369  DecayPhaseSpaceMode & operator=(const DecayPhaseSpaceMode &) = delete;
370 
371  private:
372 
376  tcDecayIntegratorPtr _integrator;
377 
381  vector<DecayPhaseSpaceChannelPtr> _channels;
382 
386  mutable vector<double> _channelwgts;
387 
391  mutable double _maxweight;
392 
396  int _niter;
397 
401  int _npoint;
402 
406  int _ntry;
407 
412 
416  int _partial;
417 
421  cGenericWidthGeneratorPtr _widthgen;
422 
426  vector<cGenericMassGeneratorPtr> _massgen;
427 
433 
437  mutable unsigned int _ichannel;
438 
442  Energy _eps;
443 
444 };
445 
449 ostream & operator<<(ostream &, const DecayPhaseSpaceMode &);
450 
451 }
452 
453 
454 #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.
tcPDPtr externalParticles(int ix) const
Access to the external particles.
double channelWeight(unsigned int ichan) const
Get the weight for a channel.
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
Energy _eps
Epsilon parameter for phase-space integration.
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.
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.
vector< PPtr > ParticleVector
vector< cGenericMassGeneratorPtr > _massgen
The mass generators for the outgoing particles.
unsigned int numberofParticles() const
Number of external particles.
Energy epsilonPS() const
Access to the epsilon parameter.
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.