DecayPhaseSpaceChannel.h

// -*- C++ -*-
//
// DecayPhaseSpaceChannel.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_DecayPhaseSpaceChannel_H
#define HERWIG_DecayPhaseSpaceChannel_H
//
// This is the declaration of the DecayPhaseSpaceChannel class.
//
#include <ThePEG/Interface/Interfaced.h>
#include <ThePEG/PDT/ParticleData.h>
#include <ThePEG/EventRecord/Particle.h>
#include "DecayPhaseSpaceChannel.fh"
#include "DecayIntegrator.h"
#include "Herwig++/Utilities/Kinematics.h"
#include "ThePEG/Repository/UseRandom.h"
#include "DecayPhaseSpaceMode.fh"

namespace Herwig {
using namespace ThePEG;

class DecayPhaseSpaceChannel: public Interfaced {

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

  friend class DecayPhaseSpaceMode;

public:
    
  DecayPhaseSpaceChannel() {}

  DecayPhaseSpaceChannel(tcDecayPhaseSpaceModePtr);
  
public:
  
  void addIntermediate(PDPtr part,int jac,double power,int dau1,int dau2) {
    _intpart.push_back(part);
    _jactype.push_back(jac);
    _intpower.push_back(power);
    _intdau1.push_back(dau1); 
    _intdau2.push_back(dau2);
  }
  
  void resetIntermediate(tcPDPtr part,Energy mass,Energy width) {
    if(!part) return;
    int idin=part->id();
    for(unsigned int ix=0;ix<_intpart.size();++ix) {
      if(_intpart[ix] && _intpart[ix]->id()==idin) {
        _intmass[ix] =mass;_intwidth[ix]=width;
        _intmass2[ix]=mass*mass;_intmwidth[ix]=mass*width;
      }
    }
  }

  /*
   * Reset the one of the daughters
   * @param oldp The id of the particle being reset
   * @param newp The id of the particle replacing it
   */
  void resetDaughter(int oldp, int newp) {
    for(unsigned int ix=0;ix<_intdau1.size();++ix) {
      if(_intdau1[ix]==oldp) _intdau1[ix]=newp;
    }
    for(unsigned int ix=0;ix<_intdau2.size();++ix) {
      if(_intdau2[ix]==oldp) _intdau2[ix]=newp;
    }
  }

protected:

  vector<Lorentz5Momentum> generateMomenta(const Lorentz5Momentum & pin,
                                           const vector<Energy> & massext);
  
  double generateWeight(const vector<Lorentz5Momentum> & output);

  void generateIntermediates(bool cc,const Particle & in, ParticleVector & out);

  void twoBodyDecay(const Lorentz5Momentum & p, 
                    const Energy m1, const Energy m2,
                    Lorentz5Momentum & p1, Lorentz5Momentum & p2);
  
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<DecayPhaseSpaceChannel> initDecayPhaseSpaceChannel;
  
  DecayPhaseSpaceChannel & operator=(const DecayPhaseSpaceChannel &);
  
private:
  
  Energy generateMass(int ires,Energy lower,Energy upper);
  
  InvEnergy2 massWeight(int ires, Energy moff,Energy lower,Energy upper);

private:
  
  tcDecayPhaseSpaceModePtr _mode;

  vector <PDPtr> _intpart;

  vector <int> _jactype;

  vector<Energy> _intmass;

  vector<Energy> _intwidth;

  vector<Energy2> _intmass2;

  vector<Energy2>_intmwidth;

  vector<double> _intpower;

  vector<int> _intdau1; 

  vector<int> _intdau2;

  vector< vector<int> > _intext;
  
  double atanhelper_(int ires, Energy limit);
};


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

class DecayPhaseSpaceError: public Exception {};
}


namespace ThePEG {

template <>

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

template <>

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

}

#endif /* HERWIG_DecayPhaseSpaceChannel_H */