HadronSelector.h

// -*- C++ -*-
//
// HadronSelector.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_HadronSelector_H
#define HERWIG_HadronSelector_H
//
// This is the declaration of the HadronSelector class.
//
 
#include "ThePEG/Interface/Interfaced.h"
#include <ThePEG/Persistency/PersistentOStream.h>
#include <ThePEG/Persistency/PersistentIStream.h>
#include <ThePEG/PDT/ParticleData.h>
#include <ThePEG/PDT/StandardMatchers.h>
#include <ThePEG/Repository/EventGenerator.h>
#include "HadronSelector.fh"
#include "HadronInfo.h"
#include "Kupco.h"
#include <tuple>
 
namespace Herwig {
 
using namespace ThePEG;
 
class HadronSelector: public Interfaced {
 
public:
 
  HadronSelector(unsigned int);
 
  virtual tcPDPair chooseHadronPair(const Energy cluMass, 
                                    tcPDPtr par1, tcPDPtr par2) const;
 
  tcPDPtr chooseSingleHadron(tcPDPtr par1, tcPDPtr par2, Energy mass) const;
 
  tcPDPair lightestHadronPair(tcPDPtr ptr1, tcPDPtr ptr2) const;
 
  Energy massLightestHadronPair(tcPDPtr ptr1, tcPDPtr ptr2) const  { 
    map<pair<long,long>,tcPDPair>::const_iterator lightest =
      lightestHadrons_.find(make_pair(abs(ptr1->id()),abs(ptr2->id())));
    if(lightest!=lightestHadrons_.end())
      return lightest->second.first->mass()+lightest->second.second->mass();
    else
      return ZERO;
  }
 
   tcPDPtr lightestHadron(tcPDPtr ptr1, tcPDPtr ptr2) const;
 
  vector<pair<tcPDPtr,double> >
  hadronsBelowThreshold(Energy threshold,
                        tcPDPtr ptr1, tcPDPtr ptr2) const;
 
   Energy massLightestHadron(tcPDPtr ptr1, tcPDPtr ptr2) const {
    // find entry in the table
    pair<long,long> ids(abs(ptr1->id()),abs(ptr2->id()));
    HadronTable::const_iterator tit=_table.find(ids);
    // throw exception if flavours wrong
    if(tit==_table.end()||tit->second.empty())
      throw Exception() <<  "HadronSelector::massLightestHadron "
                        << "failed for particle" << ptr1->id()  << " "
                        << ptr2->id()
                        << Exception::eventerror;
    // return the mass
    return tit->second.begin()->mass;
  }
 
   double pwt(long pid) const {
    map<long,double>::const_iterator it = _pwt.find(abs(pid));
    assert( it != _pwt.end() );
    return it->second;
  }
 
  virtual std::tuple<bool,bool,bool> selectBaryon(const Energy cluMass, tcPDPtr par1, tcPDPtr par2) const;
 
  virtual double strangeWeight(const Energy cluMass, tcPDPtr par1, tcPDPtr par2) const;
 
  virtual PDPtr makeDiquark(tcPDPtr par1, tcPDPtr par2);
  
public:
 
  void persistentOutput(PersistentOStream & os) const;
 
  void persistentInput(PersistentIStream & is, int version);
 
  static void Init();
 
protected:
 
  virtual void doinit();
 
protected:
 
  void insertToHadronTable(tPDPtr &particle, int flav1, int flav2);
 
  virtual void constructHadronTable();
 
  const HadronTable & table() const {
    return _table;
  }
 
  HadronTable & table() {
    return _table;
  }
  
  const vector<PDPtr> & partons() const {
    return _partons;
  }
  
  vector<PDPtr> & partons() {
    return _partons;
  }
 
  map<long,double> & pwt() {
    return _pwt;
  }
 
   double probabilityMixing(const double angleMix,
                                  const int order) const {
    static double convert=Constants::pi/180.0;
    if (order == 1)
      return sqr( cos( angleMix*convert + atan( sqrt(2.0) ) ) );
    else if (order == 2)
      return sqr( sin( angleMix*convert + atan( sqrt(2.0) ) ) );
    else
      return 1.;
  }
 
  virtual double mixingStateWeight(long id) const;
 
  double specialWeight(long id) const {
    const int pspin = id % 10;
    // Only K0L and K0S have pspin == 0, should
    // not get them until Decay step
    assert( pspin != 0 );
    // Baryon : J = 1/2 or 3/2
    if(pspin%2==0)
      return baryonWeight(id);
    // Meson
    else 
      return mesonWeight(id); 
  }
  
  virtual double mesonWeight(long id) const;
 
  virtual double baryonWeight(long id) const = 0;
 
  int signHadron(tcPDPtr ptr1, tcPDPtr ptr2, tcPDPtr hadron) const;
 
  virtual void insertMeson(HadronInfo a, int flav1, int flav2);
 
  virtual void insertOneHalf(HadronInfo a, int flav1, int flav2);
 
  virtual void insertThreeHalf(HadronInfo a, int flav1, int flav2);
  
private:
 
  HadronSelector & operator=(const HadronSelector &) = delete;
 
private:
 
  vector<PDPtr> _partons;
 
  vector<PDPtr> _forbidden;
 
  map<long,double> _pwt;
 
  double _etamix;
 
  double _phimix;
 
  double _h1mix;
 
  double _f0mix;
 
  double _f1mix;
 
  double _f2mix;
 
  double _eta2mix;
 
  double _omhmix;
 
  double _ph3mix;
 
  double _eta2Smix;
 
  double _phi2Smix;
 
  vector<double> _weight1S0;
 
  vector<double> _weight3S1;
 
  vector<double> _weight1P1;
 
  vector<double> _weight3P0;
 
  vector<double> _weight3P1;
 
  vector<double> _weight3P2;
 
  vector<double> _weight1D2;
 
  vector<double> _weight3D1;
 
  vector<double> _weight3D2;
 
  vector<double> _weight3D3;
 
  vector<vector<vector<double> > > _repwt;
 
  HadronTable _table;
 
  enum MesonMultiplets { Lmax = 3, Jmax = 4, Nmax = 4};
 
  unsigned int _topt;
 
  unsigned int _trial;
 
  double _limBottom;
  double _limCharm;
  double _limExotic;
 
  unsigned int belowThreshold_;
  
  map<pair<long,long>,tcPDPair> lightestHadrons_;
};
 
 
}
 
#endif /* HERWIG_HadronSelector_H */