ClusterFissioner.h

// -*- C++ -*-
//
// ClusterFissioner.h is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2019 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef HERWIG_ClusterFissioner_H
#define HERWIG_ClusterFissioner_H
 
#include <ThePEG/Interface/Interfaced.h>
#include "CluHadConfig.h"
#include "HadronSelector.h"
#include "ClusterFissioner.fh"
#include "CheckId.h"
 
namespace Herwig {
using namespace ThePEG;
 
  //class Cluster;          // forward declaration
 
class ClusterFissioner: public Interfaced {
 
public:
 
   ClusterFissioner();
 
 
  tPVector fission(ClusterVector & clusters, bool softUEisOn);
 
public:
 
  void persistentOutput(PersistentOStream & os) const;
 
  void persistentInput(PersistentIStream & is, int version);
 
  static void Init();
 
protected:
 
  virtual IBPtr clone() const;
 
  virtual IBPtr fullclone() const;
 
private:
 
  ClusterFissioner & operator=(const ClusterFissioner &) = delete;
 
  void cut(stack<ClusterPtr> &,
           ClusterVector&, tPVector & finalhadrons, bool softUEisOn);
 
public:
 
  typedef pair<PPtr,PPtr> PPair;
 
  typedef pair<PPair,PPair> cutType;
 
  virtual cutType cutTwo(ClusterPtr &, tPVector & finalhadrons, bool softUEisOn);
 
  virtual cutType cutThree(ClusterPtr &, tPVector & finalhadrons, bool softUEisOn);
public:
 
  PPair produceHadron(tcPDPtr hadron, tPPtr newPtr, const Lorentz5Momentum &a,
                      const LorentzPoint &b) const;
 
protected:
 
  Energy2 clustermass(const ClusterPtr & cluster) const;
 
  virtual void drawNewFlavour(PPtr& newPtr1, PPtr& newPtr2, const ClusterPtr & cluster) const {
    if (_enhanceSProb == 0){
      drawNewFlavour(newPtr1,newPtr2);
    }
    else {
      drawNewFlavourEnhanced(newPtr1,newPtr2,clustermass(cluster));
    }
  }
 
  virtual pair<Energy,Energy> drawNewMasses(Energy Mc, bool soft1, bool soft2,
                                            Lorentz5Momentum& pClu1, Lorentz5Momentum& pClu2,
                                            tPPtr ptrQ1, Lorentz5Momentum& pQ1,
                                            tPPtr, Lorentz5Momentum& pQone,
                                            tPPtr, Lorentz5Momentum& pQtwo,
                                            tPPtr ptrQ2,  Lorentz5Momentum& pQ2) const {
 
    pair<Energy,Energy> result;
 
    double exp1=_pSplitLight;
    double exp2=_pSplitLight;
 
    if     (CheckId::isExotic(ptrQ1->dataPtr())) exp1 = _pSplitExotic;
    else if(CheckId::hasBottom(ptrQ1->dataPtr()))exp1 = _pSplitBottom;
    else if(CheckId::hasCharm(ptrQ1->dataPtr())) exp1 = _pSplitCharm;
 
    if     (CheckId::isExotic(ptrQ2->dataPtr()))  exp2 = _pSplitExotic;
    else if(CheckId::hasBottom(ptrQ2->dataPtr())) exp2 = _pSplitBottom;
    else if(CheckId::hasCharm(ptrQ2->dataPtr()))  exp2 = _pSplitCharm;
 
    result.first = drawChildMass(Mc,pQ1.mass(),pQ2.mass(),pQone.mass(),exp1,soft1);
    result.second = drawChildMass(Mc,pQ2.mass(),pQ1.mass(),pQtwo.mass(),exp2,soft2);
 
    pClu1.setMass(result.first);
    pClu2.setMass(result.second);
 
    return result;
 
  }
 
  virtual void calculateKinematics(const Lorentz5Momentum &pClu,
                                   const Lorentz5Momentum &p0Q1,
                                   const bool toHadron1, const bool toHadron2,
                                   Lorentz5Momentum &pClu1, Lorentz5Momentum &pClu2,
                                   Lorentz5Momentum &pQ1, Lorentz5Momentum &pQb,
                                   Lorentz5Momentum &pQ2, Lorentz5Momentum &pQ2b) const;
 
protected:
 
  PPair produceCluster(tPPtr ptrQ, tPPtr newPtr, const Lorentz5Momentum &a,
                       const LorentzPoint &b, const Lorentz5Momentum &c,
                       const Lorentz5Momentum &d, const bool rem,
                       tPPtr spect=tPPtr(), bool remSpect=false) const;
 
  void drawNewFlavour(PPtr& newPtrPos, PPtr& newPtrNeg) const;
 
  void drawNewFlavourEnhanced(PPtr& newPtrPos,PPtr& newPtrNeg, Energy2 mass2) const;
 
 
  Energy drawChildMass(const Energy M, const Energy m1, const Energy m2,
                       const Energy m, const double exp, const bool soft) const;
 
  void calculatePositions(const Lorentz5Momentum &pClu,
                          const LorentzPoint & positionClu,
                          const Lorentz5Momentum & pClu1,
                          const Lorentz5Momentum & pClu2,
                          LorentzPoint & positionClu1,
                          LorentzPoint & positionClu2 ) const;
 
protected:
 
  HadronSelectorPtr hadronSelector() const {return _hadronSelector;}
 
  Energy btClM() const {return _btClM;}
 
  double pSplitLight() const {return _pSplitLight;}
 
  double pSplitBottom() const {return _pSplitBottom;}
 
  double pSplitCharm() const {return _pSplitCharm;}
 
  double pSplitExotic() const {return _pSplitExotic;}
 
private:
 
  bool ProbablityFunction(double scale, double threshold);
 
  bool isHeavy(tcClusterPtr );
 
  HadronSelectorPtr _hadronSelector;
 
  Energy _clMaxLight;
  Energy _clMaxBottom;
  Energy _clMaxCharm;
  Energy _clMaxExotic;
 
  double _clPowLight;
  double _clPowBottom;
  double _clPowCharm;
  double _clPowExotic;
 
  double _pSplitLight;
  double _pSplitBottom;
  double _pSplitCharm;
  double _pSplitExotic;
 
 
  // weights for alternaive cluster fission
  double _fissionPwtUquark;
  double _fissionPwtDquark;
  double _fissionPwtSquark;
 
  int _fissionCluster;
 
  int _kinematicThresholdChoice;
 
 
  Energy _btClM;
 
  int _iopRem;
 
  Tension _kappa;
 
  int _enhanceSProb;
 
  Energy _m0Fission;
 
  int _massMeasure;
 
  const double _maxScale = 20.;
 
  double _probPowFactor;
 
  double _probShift;
 
  Energy2 _kinThresholdShift;
 
};
 
}
 
#endif /* HERWIG_ClusterFissioner_H */