KinematicsReconstructor.h

// -*- C++ -*-
//
// KinematicsReconstructor.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_KinematicsReconstructor_H
#define HERWIG_KinematicsReconstructor_H
//
// This is the declaration of the KinematicsReconstructor class.
//
 
#include "ThePEG/Interface/Interfaced.h"
#include "Herwig/Shower/QTilde/Base/ShowerParticle.h"
#include "Herwig/Shower/QTilde/Base/ShowerProgenitor.h"
#include "Herwig/Shower/QTilde/Base/ShowerTree.h"
#include "Herwig/Shower/QTilde/Base/HardTree.h"
#include "KinematicsReconstructor.fh"
#include <cassert>
 
namespace Herwig {
 
using namespace ThePEG;
 
  struct KinematicsReconstructionVeto {};
 
 
struct JetKinStruct {
 
  tShowerParticlePtr parent;
 
  Lorentz5Momentum p;
 
  Lorentz5Momentum q;
};
 
typedef vector<JetKinStruct> JetKinVect;
 
enum SystemType { UNDEFINED=-1, II, IF, F ,I };
 
template<typename Value> struct ColourSinglet {
 
  typedef vector<ColourSinglet<Value> > VecType;
  
  ColourSinglet() : type(UNDEFINED) {}
  
  ColourSinglet(SystemType intype,Value inpart) 
    : type(intype),jets(1,inpart) {}
  
  SystemType type;
  
  vector<Value> jets;
 
};
 
typedef ColourSinglet<ShowerProgenitorPtr> ColourSingletSystem;
 
typedef ColourSinglet<HardBranchingPtr> ColourSingletShower;
 
class KinematicsReconstructor: public Interfaced {
 
public:
 
  KinematicsReconstructor() : _reconopt(0), _initialBoost(0),
                              _finalStateReconOption(0),
                              _initialStateReconOption(0),
                              _finalFinalWeight(false), _minQ(MeV) {};
 
  virtual bool reconstructHardJets(ShowerTreePtr hard,
                                   const map<tShowerProgenitorPtr,
                                   pair<Energy,double> > & pt,
                                   ShowerInteraction type,
                                   bool switchRecon) const;
 
  virtual bool reconstructDecayJets(ShowerTreePtr decay,
                                    ShowerInteraction type) const;
 
  virtual bool deconstructDecayJets(HardTreePtr,ShowerInteraction) const;
 
  virtual bool deconstructHardJets(HardTreePtr,ShowerInteraction) const;
 
public:
 
  void persistentOutput(PersistentOStream & os) const;
 
  void persistentInput(PersistentIStream & is, int version);
 
  static void Init();
 
public:
  
  virtual bool reconstructTimeLikeJet(const tShowerParticlePtr particleJetParent,
                                      const tShowerParticlePtr progenitor) const;
  
  double solveKfactor( const Energy & root_s, const JetKinVect & jets ) const;
  
  LorentzRotation solveBoost(const double k, 
                             const Lorentz5Momentum & newq, 
                             const Lorentz5Momentum & oldp) const;
 
  void deepTransform(PPtr particle,const LorentzRotation & r,
                     bool match=true,PPtr original=PPtr()) const;
 
protected:
 
 
  bool reconstructSpaceLikeJet(const tShowerParticlePtr particleJetParent) const;
 
  bool reconstructDecayJet(const tShowerParticlePtr particleJetParent) const;
 
  void reconstructInitialFinalSystem(vector<ShowerProgenitorPtr>) const;
 
  void reconstructFinalStateSystem(bool applyBoost, 
                                   const LorentzRotation & toRest,
                                   const LorentzRotation & fromRest, 
                                   vector<ShowerProgenitorPtr>) const;
 
  void reconstructGeneralSystem(vector<ShowerProgenitorPtr> & ShowerHardJets) const;
 
  void reconstructFinalFirst(vector<ShowerProgenitorPtr> & ShowerHardJets) const;
 
  void reconstructColourPartner(vector<ShowerProgenitorPtr> & ShowerHardJets) const;
 
  void reconstructColourSinglets(vector<ShowerProgenitorPtr> & ShowerHardJets,
                                 ShowerInteraction type) const;
 
  void reconstructInitialInitialSystem(bool & applyBoost,
                                       LorentzRotation &   toRest,
                                       LorentzRotation & fromRest,
                                       vector<ShowerProgenitorPtr>) const;
 
  void deconstructFinalStateSystem(const LorentzRotation &   toRest,
                                   const LorentzRotation & fromRest,
                                   HardTreePtr,
                                   vector<HardBranchingPtr>,
                                   ShowerInteraction) const;
  
  void deconstructInitialInitialSystem(bool & applyBoost,
                                       LorentzRotation &   toRest,
                                       LorentzRotation & fromRest,
                                       HardTreePtr,
                                       vector<HardBranchingPtr>,
                                       ShowerInteraction ) const;
 
  void deconstructInitialFinalSystem(HardTreePtr,
                                     vector<HardBranchingPtr>,
                                     ShowerInteraction ) const;
 
  bool deconstructGeneralSystem(HardTreePtr,
                                ShowerInteraction) const;
 
  bool deconstructColourSinglets(HardTreePtr,
                                 ShowerInteraction) const;
 
  bool deconstructColourPartner(HardTreePtr,
                                ShowerInteraction) const;
 
  void reconstructFinalFinalOffShell(JetKinVect orderedJets, Energy2 s,
                                     bool recursive) const;
 
  LorentzRotation solveBoost(const Lorentz5Momentum & newq, 
                             const Lorentz5Momentum & oldq) const;
  
  LorentzRotation solveBoostZ(const Lorentz5Momentum & newq, 
                              const Lorentz5Momentum & oldq) const;
  
  void boostChain(tPPtr p, const LorentzRotation & bv, tPPtr & parent) const;
 
  Boost solveBoostBeta( const double k, const Lorentz5Momentum & newq, 
                          const Lorentz5Momentum & oldp);
 
  double getBeta(const double E, const double q, 
                 const double Ep, const double qp) const
  {return (q*E-qp*Ep)/(sqr(qp)+sqr(E));}
 
 
  bool solveDecayKFactor(Energy mb, 
                         const Lorentz5Momentum & n, 
                         const Lorentz5Momentum & pjet, 
                         const JetKinVect & jetKinematics, 
                         ShowerParticlePtr partner,
                         Lorentz5Momentum ppartner[2],
                         double & k1, 
                         double & k2,
                         Lorentz5Momentum & qt) const;
 
  double inverseRescalingFactor(vector<Lorentz5Momentum> pout,
                                vector<Energy> mon,Energy roots) const;
 
  bool inverseDecayRescalingFactor(vector<Lorentz5Momentum> pout,
                                   vector<Energy> mon,Energy roots,
                                   Lorentz5Momentum ppartner, Energy mbar,
                                   double & k1, double & k2) const;
 
  Energy momConsEq(double k, const Energy & root_s,
                   const JetKinVect & jets) const;
 
 
  void findInitialBoost(const Lorentz5Momentum & pold, const Lorentz5Momentum & pnew,
                        LorentzRotation & toRest, LorentzRotation & fromRest) const;
 
  template<typename Value> void findPartners(Value branch,set<Value> & done,
                                             const set<Value> & branchings,
                                             vector<Value> & jets) const;
 
  bool addIntrinsicPt(vector<ShowerProgenitorPtr>) const;
 
  Energy findMass(HardBranchingPtr) const;
 
  vector<double> initialStateRescaling(double x1, double x2, 
                                       const Lorentz5Momentum & pold,
                                       const vector<Lorentz5Momentum> & p,
                                       const vector<Lorentz5Momentum> & pq,
                                       const vector<Energy>& highespts) const;
 
  vector<double> inverseInitialStateRescaling(double & x1, double & x2,
                                              const Lorentz5Momentum & pold,
                                              const vector<Lorentz5Momentum> & p,
                                              const vector<Lorentz5Momentum> & pq) const;
 
  template<typename Value >
  typename ColourSinglet<Value>::VecType identifySystems(set<Value> jets,
                                                         unsigned int & nnun,unsigned int & nnii,
                                                         unsigned int & nnif,unsigned int & nnf,
                                                         unsigned int & nni) const;
 
  template<typename Value>
  void combineFinalState(vector<ColourSinglet<Value> > & systems) const;
 
protected:
 
  virtual IBPtr clone() const {return new_ptr(*this);}
 
  virtual IBPtr fullclone() const {return new_ptr(*this);}
 
protected:
 
  virtual void doinit();
 
private:
 
  KinematicsReconstructor & operator=(const KinematicsReconstructor &) = delete;
 
private:
 
  unsigned int _reconopt;
 
  unsigned int _initialBoost;
 
  unsigned int _finalStateReconOption;
 
  unsigned int _initialStateReconOption;
 
  bool _finalFinalWeight;
 
  Energy _minQ;
 
  mutable map<tShowerProgenitorPtr,pair<Energy,double> > _intrinsic;
 
  mutable tShowerTreePtr _currentTree;
 
  PDVector _noRescaleVector;
 
  set<cPDPtr> _noRescale;
 
  mutable map<tPPtr,vector<LorentzRotation> > _boosts;
 
  mutable map<tShowerTreePtr,vector<LorentzRotation> > _treeBoosts;
};
 
}
 
#endif /* HERWIG_KinematicsReconstructor_H */