DecayPhaseSpaceMode.h

// -*- C++ -*-
//
// DecayPhaseSpaceMode.h is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2011 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef HERWIG_DecayPhaseSpaceMode_H
#define HERWIG_DecayPhaseSpaceMode_H
//
// This is the declaration of the DecayPhaseSpaceMode class.
//
#include "ThePEG/Interface/Interfaced.h"
#include "DecayPhaseSpaceMode.fh"
#include "DecayPhaseSpaceChannel.h"
#include "Herwig++/PDT/GenericWidthGenerator.h"
#include "Herwig++/PDT/GenericMassGenerator.h"
#include <Herwig++/Decay/DecayVertex.h>
#include "DecayIntegrator.h"

namespace Herwig {
using namespace ThePEG;

class DecayPhaseSpaceMode: public Interfaced {

  friend ostream & operator<<(ostream &, const DecayPhaseSpaceMode &);

  friend class DecayIntegrator;

  friend class DecayPhaseSpaceChannel;

public:

  DecayPhaseSpaceMode() :  _maxweight(0.),_niter(10), _npoint(10000), _ntry(500),
                           _partial(-1), _testOnShell(false), _ichannel(999) {}

  DecayPhaseSpaceMode(tPDVector in, tcDecayIntegratorPtr intin,bool onShell=false) 
    :  _integrator(intin), _maxweight(0.),
       _niter(10), _npoint(10000), _ntry(500),
       _extpart(in),  _partial(-1), _testOnShell(onShell), _ichannel(999) {}

  tcPDPtr externalParticles(int ix) const {return _extpart[ix];}

  unsigned int numberofParticles() const {return _extpart.size();}

  unsigned int numberChannels() const {return _channels.size();}

  void addChannel(DecayPhaseSpaceChannelPtr channel) {
    channel->init();
    _channels.push_back(channel);
  }

  void resetIntermediate(int ichan, tcPDPtr part, Energy mass, Energy width) {
    if(!part) return;
    _channels[ichan]->resetIntermediate(part,mass,width);
  }

  void resetIntermediate(tcPDPtr part, Energy mass, Energy width) {
    for(unsigned int ix=0,N=_channels.size();ix<N;++ix)
      _channels[ix]->resetIntermediate(part,mass,width);
  }

  double maxWeight() const {return _maxweight;}

  void setMaxWeight(double wgt) const {_maxweight=wgt;}

  double channelWeight(unsigned int ichan) const {return _channelwgts[ichan];}

  void setWeights(const vector<double> & in) const {_channelwgts=in;}

  unsigned int selectedChannel() const {return _ichannel;}

  bool testOnShell() const { return _testOnShell; }

protected:

  Energy initializePhaseSpace(bool init, bool onShell=false);

  void setIntegrate(int iter,int points,int ntry) {
    _niter=iter;
    _npoint=points;
    _ntry=ntry;
  }

  void setPartialWidth(int in) {_partial=in;}



  double me2(const int ichan ,const Particle & inpart,
             const ParticleVector &outpart,DecayIntegrator::MEOption opt) const {
    return _integrator->me2(ichan,inpart,outpart,opt);
  }
  
  ParticleVector generate(bool intermediates,bool cc,const Particle & inpart) const;


  int selectChannel(const Particle & inpart, ParticleVector & outpart) const {
    // if using flat phase-space don't need to do this
    if(_channelwgts.empty()) return 0;
    vector<double> mewgts(_channels.size(),0.0);
    double total=0.;
    for(unsigned int ix=0,N=_channels.size();ix<N;++ix) {
      mewgts[ix]=me2(ix,inpart,outpart,DecayIntegrator::Calculate);
      total+=mewgts[ix];
    }
    // randomly pick a channel
    total*=UseRandom::rnd();
    int ichan=-1;
    do {
      ++ichan;
      total-=mewgts[ichan];
    }
    while(ichan<int(_channels.size())&&total>0.);
    return ichan;
  }

  Energy weight(bool cc,int & ichan, const Particle & in,
                ParticleVector & particles,bool first,
                bool onShell=false) const {
    ichan=0;
    Energy phwgt = (_channels.size()==0) ? 
      flatPhaseSpace(cc,in,particles,onShell) : channelPhaseSpace(cc,ichan,in,particles,onShell);
    // generate the matrix element
    return me2(-1,in,particles,
               first ? DecayIntegrator::Initialize : DecayIntegrator::Calculate)*phwgt;
  }
    
  Energy flatPhaseSpace(bool cc,const Particle & inpart, ParticleVector & outpart,
                        bool onShell=false) const;
  
  Energy channelPhaseSpace(bool cc,int & ichan, const Particle & in, 
                           ParticleVector & particles,
                           bool onShell=false) const;

  void constructVertex(const Particle & in, const ParticleVector & out) const;

  vector<Energy> externalMasses(Energy inmass,double & wgt, bool onShell) const;

public:

  void persistentOutput(PersistentOStream & os) const;

  void persistentInput(PersistentIStream & is, int version);

  static void Init();

protected:

  virtual IBPtr clone() const {return new_ptr(*this);}

  virtual IBPtr fullclone() const {return new_ptr(*this);}

protected:

  virtual void doinit();

  virtual void doinitrun();

private:

  static ClassDescription<DecayPhaseSpaceMode> initDecayPhaseSpaceMode;

  DecayPhaseSpaceMode & operator=(const DecayPhaseSpaceMode &);

 private:

  tcDecayIntegratorPtr _integrator;

  vector<DecayPhaseSpaceChannelPtr> _channels;

  mutable vector<double> _channelwgts;

  mutable double _maxweight;

  int _niter;

  int _npoint;

  int _ntry;

  tPDVector _extpart;

  int _partial;

  cGenericWidthGeneratorPtr _widthgen;

  vector<cGenericMassGeneratorPtr> _massgen;

  bool _testOnShell;

  mutable unsigned int _ichannel;

};

ostream & operator<<(ostream &, const DecayPhaseSpaceMode &);

}


#include "ThePEG/Utilities/ClassTraits.h"

namespace ThePEG {

template <>
 struct BaseClassTrait<Herwig::DecayPhaseSpaceMode,1> {
  typedef Interfaced NthBase;
};

template <>
 struct ClassTraits<Herwig::DecayPhaseSpaceMode>
  : public ClassTraitsBase<Herwig::DecayPhaseSpaceMode> {
   static string className() { return "Herwig::DecayPhaseSpaceMode"; }
};

}

#endif /* HERWIG_DecayPhaseSpaceMode_H */