QTilde/Kinematics/KinematicsReconstructor.h

0001 // -*- C++ -*-
0002 //
0003 // KinematicsReconstructor.h is a part of Herwig - A multi-purpose Monte Carlo event generator
0004 // Copyright (C) 2002-2019 The Herwig Collaboration
0005 //
0006 // Herwig is licenced under version 3 of the GPL, see COPYING for details.
0007 // Please respect the MCnet academic guidelines, see GUIDELINES for details.
0008 //
0009 #ifndef HERWIG_KinematicsReconstructor_H
0010 #define HERWIG_KinematicsReconstructor_H
0011 //
0012 // This is the declaration of the KinematicsReconstructor class.
0013 //
0014
0015 #include "ThePEG/Interface/Interfaced.h"
0016 #include "Herwig/Shower/QTilde/Base/ShowerParticle.h"
0017 #include "Herwig/Shower/QTilde/Base/ShowerProgenitor.h"
0018 #include "Herwig/Shower/QTilde/Base/ShowerTree.h"
0019 #include "Herwig/Shower/QTilde/Base/HardTree.h"
0020 #include "KinematicsReconstructor.fh"
0021 #include <cassert>
0022
0023 namespace Herwig {
0024
0025 using namespace ThePEG;
0026
0027   /**\ingroup Shower
0028    * Exception class
0029    * used to communicate failure of kinematics
0030    * reconstruction.
0031    */
0032   struct KinematicsReconstructionVeto {};
0033
0034
0035 /** \ingroup Shower
0036  *  A simple struct to store the information we need on the
0037  *  showering
0038  */
0039 struct JetKinStruct {
0040
0041   /**
0042    *  Parent particle of the jet
0043    */
0044   tShowerParticlePtr parent;
0045
0046   /**
0047    *  Momentum of the particle before reconstruction
0048    */
0049   Lorentz5Momentum p;
0050
0051   /**
0052    *  Momentum of the particle after reconstruction
0053    */
0054   Lorentz5Momentum q;
0055 };
0056
0057 /**
0058  * typedef for a vector of JetKinStruct
0059  */
0060 typedef vector<JetKinStruct> JetKinVect;
0061
0062 /**
0063  *  Enum to identify types of colour singlet systems
0064  */
0065 enum SystemType { UNDEFINED=-1, II, IF, F ,I };
0066
0067 /**
0068  *  Struct to store colour singlets
0069  */
0070 template<typename Value> struct ColourSinglet {
0071
0072   typedef vector<ColourSinglet<Value> > VecType;
0073
0074   ColourSinglet() : type(UNDEFINED) {}
0075
0076   ColourSinglet(SystemType intype,Value inpart)
0077     : type(intype),jets(1,inpart) {}
0078
0079   /**
0080    * The type of system
0081    */
0082   SystemType type;
0083
0084   /**
0085    *  The jets in the system
0086    */
0087   vector<Value> jets;
0088
0089 };
0090
0091 /**
0092  *  Struct to store a colour singlet system of particles
0093  */
0094 typedef ColourSinglet<ShowerProgenitorPtr> ColourSingletSystem;
0095
0096 /**
0097  * Struct to store a colour singlet shower
0098  */
0099 typedef ColourSinglet<HardBranchingPtr> ColourSingletShower;
0100
0101 /** \ingroup Shower
0102  *
0103  * This class is responsible for the kinematical reconstruction
0104  * after each showering step, and also for the necessary Lorentz boosts
0105  * in order to preserve energy-momentum conservation in the overall collision,
0106  * and also the invariant mass and the rapidity of the hard subprocess system.
0107  * In the case of multi-step showering, there will be not unnecessary
0108  * kinematical reconstructions.
0109  *
0110  * There is also the option of taking a set of momenta for the particles
0111  * and inverting the reconstruction to give the evolution variables for the
0112  * shower.
0113  *
0114  * Notice:
0115  * - although we often use the term "jet" in either methods or variables names,
0116  *   or in comments, which could appear applicable only for QCD showering,
0117  *   there is indeed no "dynamics" represented in this class: only kinematics
0118  *   is involved, as the name of this class remainds. Therefore it can be used
0119  *   for any kind of showers (QCD-,QED-,EWK-,... bremsstrahlung).
0120  *
0121  * @see ShowerParticle
0122  * @see ShowerKinematics
0123  * @see \ref KinematicsReconstructorInterfaces "The interfaces"
0124  * defined for KinematicsReconstructor.
0125  */
0126 class KinematicsReconstructor: public Interfaced {
0127
0128 public:
0129
0130   /**
0131    *  Default constructor
0132    */
0133   KinematicsReconstructor() : _reconopt(0), _initialBoost(0),
0134                   _finalStateReconOption(0),
0135                   _initialStateReconOption(0),
0136                   _finalFinalWeight(false), _minQ(MeV) {};
0137
0138   /**
0139    *  Methods to reconstruct the kinematics of a scattering or decay process
0140    */
0141   //@{
0142   /**
0143    * Given in input a vector of the particles which initiated the showers
0144    * the method does the reconstruction of such jets,
0145    * including the appropriate boosts (kinematics reshufflings)
0146    * needed to conserve the total energy-momentum of the collision
0147    * and preserving the invariant mass and the rapidity of the
0148    * hard subprocess system.
0149    */
0150   virtual bool reconstructHardJets(ShowerTreePtr hard,
0151                    const map<tShowerProgenitorPtr,
0152                    pair<Energy,double> > & pt,
0153                    ShowerInteraction type,
0154                    bool switchRecon) const;
0155
0156   /**
0157    * Given in input a vector of the particles which initiated the showers
0158    * the method does the reconstruction of such jets,
0159    * including the appropriate boosts (kinematics reshufflings)
0160    * needed to conserve the total energy-momentum of the collision
0161    * and preserving the invariant mass and the rapidity of the
0162    * hard subprocess system.
0163    */
0164   virtual bool reconstructDecayJets(ShowerTreePtr decay,
0165                     ShowerInteraction type) const;
0166   //@}
0167
0168   /**
0169    *  Methods to invert the reconstruction of the shower for
0170    *  a scattering or decay process and calculate
0171    *  the variables used to generate the
0172    *  shower given the particles produced.
0173    *  This is needed for the CKKW and POWHEG approaches
0174    */
0175   //@{
0176   /**
0177    *  Given the particles, with a history which we wish to interpret
0178    *  as a shower reconstruct the variables used to generate the
0179    * shower
0180    */
0181   virtual bool deconstructDecayJets(HardTreePtr,ShowerInteraction) const;
0182
0183   /**
0184    *  Given the particles, with a history which we wish to interpret
0185    *  as a shower reconstruct the variables used to generate the shower
0186    *  for a hard process
0187    */
0188   virtual bool deconstructHardJets(HardTreePtr,ShowerInteraction) const;
0189   //@}
0190
0191 public:
0192
0193   /** @name Functions used by the persistent I/O system. */
0194   //@{
0195   /**
0196    * Function used to write out object persistently.
0197    * @param os the persistent output stream written to.
0198    */
0199   void persistentOutput(PersistentOStream & os) const;
0200
0201   /**
0202    * Function used to read in object persistently.
0203    * @param is the persistent input stream read from.
0204    * @param version the version number of the object when written.
0205    */
0206   void persistentInput(PersistentIStream & is, int version);
0207   //@}
0208
0209   /**
0210    * The standard Init function used to initialize the interfaces.
0211    * Called exactly once for each class by the class description system
0212    * before the main function starts or
0213    * when this class is dynamically loaded.
0214    */
0215   static void Init();
0216
0217 public:
0218
0219   /**
0220    * Given the particle (ShowerParticle object) that
0221    * originates a forward (time-like) jet, this method reconstructs the kinematics
0222    * of the jet. That is, by starting from the final grand-children (which
0223    * originates directly or indirectly from particleJetParent,
0224    * and which don't have children), and moving "backwards" (in a physical
0225    * time picture), towards the particleJetParent, the
0226    * ShowerKinematics objects associated with the various particles,
0227    * which have been created during the showering, are now completed.
0228    * In particular, at the end, we get the mass of the jet, which is the
0229    * main information we want.
0230    * This methods returns false if there was no radiation or rescaling required
0231    */
0232   virtual bool reconstructTimeLikeJet(const tShowerParticlePtr particleJetParent,
0233                       const tShowerParticlePtr progenitor) const;
0234
0235   /**
0236    * Given a vector of 5-momenta of jets, where the 3-momenta are the initial
0237    * ones before showering and the masses are reconstructed after the showering,
0238    * this method returns the overall scaling factor for the 3-momenta of the
0239    * vector of particles, vec{P}_i -> k * vec{P}_i, such to preserve energy-
0240    * momentum conservation, i.e. after the rescaling the center of mass 5-momentum
0241    * is equal to the one specified in input, cmMomentum.
0242    * The method returns 0 if such factor cannot be found.
0243    * @param root_s Centre-of-mass energy
0244    * @param jets The jets
0245    */
0246   double solveKfactor( const Energy & root_s, const JetKinVect & jets ) const;
0247
0248   /**
0249    * Compute the boost to get from the the old momentum to the new
0250    */
0251   LorentzRotation solveBoost(const double k,
0252                  const Lorentz5Momentum & newq,
0253                  const Lorentz5Momentum & oldp) const;
0254
0255   /**
0256    *  Apply a transform to the particle and any child, including child ShowerTree
0257    *  objects
0258    * @param particle The particle
0259    * @param r The Lorentz transformation
0260    * @param match Whether or not to look at children etc
0261    * @param original The original particle
0262    */
0263   void deepTransform(PPtr particle,const LorentzRotation & r,
0264              bool match=true,PPtr original=PPtr()) const;
0265
0266 protected:
0267
0268   /**
0269    *  Methods to reconstruct the kinematics of individual jets
0270    */
0271   //@{
0272
0273   /**
0274    * Exactly similar to the previous one, but for a space-like jet.
0275    * Also in this case we start from the final grand-children (which
0276    * are childless) of the particle which originates the jet, but in
0277    * this case we proceed "forward" (in the physical time picture)
0278    * towards the particleJetParent.
0279    * This methods returns false if there was no radiation or rescaling required
0280    */
0281   bool reconstructSpaceLikeJet(const tShowerParticlePtr particleJetParent) const;
0282
0283   /**
0284    * Exactly similar to the previous one, but for a decay jet
0285    * This methods returns false if there was no radiation or rescaling required
0286    */
0287   bool reconstructDecayJet(const tShowerParticlePtr particleJetParent) const;
0288   //@}
0289
0290   /**
0291    *  Methods to perform the reconstruction of various types of colour
0292    *  singlet systems
0293    */
0294   //@{
0295   /**
0296    *  Perform the reconstruction of a system with one incoming and at least one
0297    *  outgoing particle
0298    */
0299   void reconstructInitialFinalSystem(vector<ShowerProgenitorPtr>) const;
0300
0301   /**
0302    *  Perform the reconstruction of a system with only final-state
0303    *  particles
0304    */
0305   void reconstructFinalStateSystem(bool applyBoost,
0306                    const LorentzRotation & toRest,
0307                    const LorentzRotation & fromRest,
0308                    vector<ShowerProgenitorPtr>) const;
0309
0310   /**
0311    *  Reconstruction of a general coloured system
0312    */
0313   void reconstructGeneralSystem(vector<ShowerProgenitorPtr> & ShowerHardJets) const;
0314
0315   /**
0316    * Reconstruction of a general coloured system doing
0317    * final-final, then initial-final and then initial-initial
0318    */
0319   void reconstructFinalFirst(vector<ShowerProgenitorPtr> & ShowerHardJets) const;
0320
0321   /**
0322    *  Reconstruction of a general coloured system doing
0323    *  colour parners
0324    */
0325   void reconstructColourPartner(vector<ShowerProgenitorPtr> & ShowerHardJets) const;
0326
0327   /**
0328    *  Reconstruction based on colour singlet systems
0329    */
0330   void reconstructColourSinglets(vector<ShowerProgenitorPtr> & ShowerHardJets,
0331                  ShowerInteraction type) const;
0332
0333   /**
0334    *  Perform the reconstruction of a system with only final-state
0335    *  particles
0336    */
0337   void reconstructInitialInitialSystem(bool & applyBoost,
0338                        LorentzRotation &   toRest,
0339                        LorentzRotation & fromRest,
0340                        vector<ShowerProgenitorPtr>) const;
0341   //@}
0342
0343   /**
0344    *  Methods to perform the inverse reconstruction of various types of
0345    *  colour singlet systems
0346    */
0347   //@{
0348   /**
0349    *  Perform the inverse reconstruction of a system with only final-state
0350    *  particles
0351    */
0352   void deconstructFinalStateSystem(const LorentzRotation &   toRest,
0353                    const LorentzRotation & fromRest,
0354                    HardTreePtr,
0355                    vector<HardBranchingPtr>,
0356                    ShowerInteraction) const;
0357
0358   /**
0359    *  Perform the inverse reconstruction of a system with only initial-state
0360    *  particles
0361    */
0362   void deconstructInitialInitialSystem(bool & applyBoost,
0363                        LorentzRotation &   toRest,
0364                        LorentzRotation & fromRest,
0365                        HardTreePtr,
0366                        vector<HardBranchingPtr>,
0367                        ShowerInteraction ) const;
0368
0369   /**
0370    *  Perform the inverse reconstruction of a system with only initial-state
0371    *  particles
0372    */
0373   void deconstructInitialFinalSystem(HardTreePtr,
0374                      vector<HardBranchingPtr>,
0375                      ShowerInteraction ) const;
0376
0377   bool deconstructGeneralSystem(HardTreePtr,
0378                 ShowerInteraction) const;
0379
0380   bool deconstructColourSinglets(HardTreePtr,
0381                  ShowerInteraction) const;
0382
0383   bool deconstructColourPartner(HardTreePtr,
0384                 ShowerInteraction) const;
0385   //@}
0386
0387   /**
0388    *   Recursively treat the most off-shell paricle seperately
0389    * for final-final reconstruction
0390    */
0391   void reconstructFinalFinalOffShell(JetKinVect orderedJets, Energy2 s,
0392                      bool recursive) const;
0393
0394   /**
0395    *  Various methods for the Lorentz transforms needed to do the
0396    *  rescalings
0397    */
0398   //@{
0399   /**
0400    * Compute the boost to get from the the old momentum to the new
0401    */
0402   LorentzRotation solveBoost(const Lorentz5Momentum & newq,
0403                  const Lorentz5Momentum & oldq) const;
0404
0405   /**
0406    * Compute the boost to get from the the old momentum to the new
0407    */
0408   LorentzRotation solveBoostZ(const Lorentz5Momentum & newq,
0409                   const Lorentz5Momentum & oldq) const;
0410
0411   /**
0412    *  Recursively boost the initial-state shower
0413    * @param p The particle
0414    * @param bv The boost
0415    * @param parent The parent of the chain
0416    */
0417   void boostChain(tPPtr p, const LorentzRotation & bv, tPPtr & parent) const;
0418
0419   /**
0420    * Given a 5-momentum and a scale factor, the method returns the
0421    * Lorentz boost that transforms the 3-vector vec{momentum} --->
0422    * k*vec{momentum}. The method returns the null boost in the case no
0423    * solution exists. This will only work in the case where the
0424    * outgoing jet-momenta are parallel to the momenta of the particles
0425    * leaving the hard subprocess.
0426    */
0427   Boost solveBoostBeta( const double k, const Lorentz5Momentum & newq,
0428               const Lorentz5Momentum & oldp);
0429
0430   /**
0431    * Compute boost parameter along z axis to get (Ep, any perp, qp)
0432    * from (E, same perp, q).
0433    */
0434   double getBeta(const double E, const double q,
0435          const double Ep, const double qp) const
0436   {return (q*E-qp*Ep)/(sqr(qp)+sqr(E));}
0437   //@}
0438
0439   /**
0440    *  Methods to calculate the various scaling factors
0441    */
0442   //@{
0443
0444   /**
0445    *  Calculate the rescaling factors for the jets in a particle decay where
0446    *  there was initial-state radiation
0447    * @param mb The mass of the decaying particle
0448    * @param n  The reference vector for the initial state radiation
0449    * @param pjet The momentum of the initial-state jet
0450    * @param jetKinematics The JetKinStruct objects for the jets
0451    * @param partner The colour partner
0452    * @param ppartner The momentum of the colour partner of the decaying particle
0453    * before and after radiation
0454    * @param k1 The rescaling parameter for the partner
0455    * @param k2 The rescaling parameter for the outgoing singlet
0456    * @param qt The transverse momentum vector
0457    */
0458   bool solveDecayKFactor(Energy mb,
0459              const Lorentz5Momentum & n,
0460              const Lorentz5Momentum & pjet,
0461              const JetKinVect & jetKinematics,
0462              ShowerParticlePtr partner,
0463              Lorentz5Momentum ppartner[2],
0464              double & k1,
0465              double & k2,
0466              Lorentz5Momentum & qt) const;
0467
0468   /**
0469    * Compute the momentum rescaling factor needed to invert the shower
0470    * @param pout The momenta of the outgoing particles
0471    * @param mon  The on-shell masses
0472    * @param roots The mass of the decaying particle
0473    */
0474   double inverseRescalingFactor(vector<Lorentz5Momentum> pout,
0475                 vector<Energy> mon,Energy roots) const;
0476
0477   /**
0478    * Compute the momentum rescaling factor needed to invert the shower
0479    * @param pout The momenta of the outgoing particles
0480    * @param mon  The on-shell masses
0481    * @param roots The mass of the decaying particle
0482    * @param ppartner The momentum of the colour partner
0483    * @param mbar The mass of the decaying particle
0484    * @param k1 The first scaling factor
0485    * @param k2 The second scaling factor
0486    */
0487   bool inverseDecayRescalingFactor(vector<Lorentz5Momentum> pout,
0488                    vector<Energy> mon,Energy roots,
0489                    Lorentz5Momentum ppartner, Energy mbar,
0490                    double & k1, double & k2) const;
0491
0492   /**
0493    * Check the rescaling conserves momentum
0494    * @param k The rescaling
0495    * @param root_s The centre-of-mass energy
0496    * @param jets The jets
0497    */
0498   Energy momConsEq(double k, const Energy & root_s,
0499            const JetKinVect & jets) const;
0500
0501
0502   void findInitialBoost(const Lorentz5Momentum & pold, const Lorentz5Momentum & pnew,
0503             LorentzRotation & toRest, LorentzRotation & fromRest) const;
0504   //@}
0505
0506   /**
0507    *  Find the colour partners of a particle to identify the colour singlet
0508    *  systems for the reconstruction.
0509    */
0510   template<typename Value> void findPartners(Value branch,set<Value> & done,
0511                          const set<Value> & branchings,
0512                          vector<Value> & jets) const;
0513
0514   /**
0515    *  Add the intrinsic \f$p_T\f$ to the system if needed
0516    */
0517   bool addIntrinsicPt(vector<ShowerProgenitorPtr>) const;
0518
0519   /**
0520    *  Find the mass of a particle in the hard branching
0521    */
0522   Energy findMass(HardBranchingPtr) const;
0523
0524   /**
0525    *  Calculate the initial-state rescaling factors
0526    */
0527   vector<double> initialStateRescaling(double x1, double x2,
0528                        const Lorentz5Momentum & pold,
0529                        const vector<Lorentz5Momentum> & p,
0530                        const vector<Lorentz5Momentum> & pq,
0531                        const vector<Energy>& highespts) const;
0532
0533   /**
0534    *  Calculate the inverse of the initial-state rescaling factor
0535    */
0536   vector<double> inverseInitialStateRescaling(double & x1, double & x2,
0537                           const Lorentz5Momentum & pold,
0538                           const vector<Lorentz5Momentum> & p,
0539                           const vector<Lorentz5Momentum> & pq) const;
0540
0541   /**
0542    *  Find the colour singlet systems
0543    */
0544   template<typename Value >
0545   typename ColourSinglet<Value>::VecType identifySystems(set<Value> jets,
0546                              unsigned int & nnun,unsigned int & nnii,
0547                              unsigned int & nnif,unsigned int & nnf,
0548                              unsigned int & nni) const;
0549
0550   /**
0551    *  Combine final-state colour systems
0552    */
0553   template<typename Value>
0554   void combineFinalState(vector<ColourSinglet<Value> > & systems) const;
0555
0556 protected:
0557
0558   /** @name Clone Methods. */
0559   //@{
0560   /**
0561    * Make a simple clone of this object.
0562    * @return a pointer to the new object.
0563    */
0564   virtual IBPtr clone() const {return new_ptr(*this);}
0565
0566   /** Make a clone of this object, possibly modifying the cloned object
0567    * to make it sane.
0568    * @return a pointer to the new object.
0569    */
0570   virtual IBPtr fullclone() const {return new_ptr(*this);}
0571   //@}
0572
0573 protected:
0574
0575   /** @name Standard Interfaced functions. */
0576   //@{
0577   /**
0578    * Initialize this object after the setup phase before saving an
0579    * EventGenerator to disk.
0580    * @throws InitException if object could not be initialized properly.
0581    */
0582   virtual void doinit();
0583   //@}
0584
0585 private:
0586
0587   /**
0588    * The assignment operator is private and must never be called.
0589    * In fact, it should not even be implemented.
0590    */
0591   KinematicsReconstructor & operator=(const KinematicsReconstructor &) = delete;
0592
0593 private:
0594
0595   /**
0596    *  Option for handling the reconstruction
0597    */
0598   unsigned int _reconopt;
0599
0600   /**
0601    *  Option for the boost for initial-initial reconstruction
0602    */
0603   unsigned int _initialBoost;
0604
0605   /**
0606    * Option for the reconstruction of final state systems
0607    */
0608   unsigned int _finalStateReconOption;
0609
0610   /**
0611    * Option for the initial state reconstruction
0612    */
0613   unsigned int _initialStateReconOption;
0614
0615   /**
0616    *  Option for FF kinematic factor
0617    */
0618   bool _finalFinalWeight;
0619
0620   /**
0621    * Minimum invariant mass for initial-final dipoles to allow the
0622    * reconstruction
0623    */
0624   Energy _minQ;
0625
0626   /**
0627    * Storage of the intrinsic \f$p_T\f$
0628    */
0629   mutable map<tShowerProgenitorPtr,pair<Energy,double> > _intrinsic;
0630
0631   /**
0632    *  Current ShowerTree
0633    */
0634   mutable tShowerTreePtr _currentTree;
0635
0636   /**
0637    * Particles which shouldn't have their masses rescaled as
0638    * vector for the interface
0639    */
0640   PDVector _noRescaleVector;
0641
0642   /**
0643    * Particles which shouldn't have their masses rescaled as
0644    * set for quick access
0645    */
0646   set<cPDPtr> _noRescale;
0647
0648   /**
0649    * Storage of the boosts applied to enable resetting after failure
0650    */
0651   mutable map<tPPtr,vector<LorentzRotation> > _boosts;
0652
0653   /**
0654    * Storage of the boosts applied to enable resetting after failure
0655    */
0656   mutable map<tShowerTreePtr,vector<LorentzRotation> > _treeBoosts;
0657 };
0658
0659 }
0660
0661 #endif /* HERWIG_KinematicsReconstructor_H */