PYTHIA8/pythia8175/include/MergingHooks.h

   1 // MergingHooks.h is a part of the PYTHIA event generator.
   2 // Copyright (C) 2013 Torbjorn Sjostrand.
   3 // PYTHIA is licenced under the GNU GPL version 2, see COPYING for details.
   4 // Please respect the MCnet Guidelines, see GUIDELINES for details.
   5
   6 // This file is written by Stefan Prestel.
   7 // Header file to allow user access to program at different stages.
   8 // HardProcess: Container class for the hard process to be merged. Holds the
   9 //              bookkeeping of particles not be be reclustered
  10 // MergingHooks: Steering class for matrix element merging. Some functions can
  11 //               be redefined in a derived class to have access to the merging
  12
  13 #ifndef Pythia8_MergingHooks_H
  14 #define Pythia8_MergingHooks_H
  15
  16 #include "Basics.h"
  17 #include "BeamParticle.h"
  18 #include "Event.h"
  19 #include "Info.h"
  20 #include "ParticleData.h"
  21 #include "PartonSystems.h"
  22 #include "PythiaStdlib.h"
  23 #include "Settings.h"
  24
  25 namespace Pythia8 {
  26
  27 //==========================================================================
  28
  29 // Declaration of hard process class
  30 // This class holds information on the desired hard 2->2 process
  31 // for the merging.
  32 // This class is a container class for History class use.
  33
  34 class HardProcess {
  35
  36 public:
  37
  38    // Flavour of the first incoming particle
  39   int hardIncoming1;
  40   // Flavour of the second incoming particle
  41   int hardIncoming2;
  42   // Flavours of the outgoing particles
  43   vector<int> hardOutgoing1;
  44   vector<int> hardOutgoing2;
  45   // Flavour of intermediate bosons in the hard 2->2
  46   vector<int> hardIntermediate;
  47
  48   // Current reference event
  49   Event state;
  50   // Potential positions of outgoing particles in reference event
  51   vector<int> PosOutgoing1;
  52   vector<int> PosOutgoing2;
  53   // Potential positions of intermediate bosons in reference event
  54   vector<int> PosIntermediate;
  55
  56   // Information on merging scale read from LHE file
  57   double tms;
  58
  59   // Default constructor
  60   HardProcess(){}
  61   // Default destructor
  62   ~HardProcess(){}
  63
  64   // Copy constructor
  65   HardProcess( const HardProcess& hardProcessIn )
  66     : state(hardProcessIn.state),
  67       tms(hardProcessIn.tms) {
  68     hardIncoming1 = hardProcessIn.hardIncoming1;
  69     hardIncoming2 = hardProcessIn.hardIncoming2;
  70     for(int i =0; i < int(hardProcessIn.hardOutgoing1.size());++i)
  71       hardOutgoing1.push_back( hardProcessIn.hardOutgoing1[i]);
  72     for(int i =0; i < int(hardProcessIn.hardOutgoing2.size());++i)
  73       hardOutgoing2.push_back( hardProcessIn.hardOutgoing2[i]);
  74     for(int i =0; i < int(hardProcessIn.hardIntermediate.size());++i)
  75       hardIntermediate.push_back( hardProcessIn.hardIntermediate[i]);
  76     for(int i =0; i < int(hardProcessIn.PosOutgoing1.size());++i)
  77       PosOutgoing1.push_back( hardProcessIn.PosOutgoing1[i]);
  78     for(int i =0; i < int(hardProcessIn.PosOutgoing2.size());++i)
  79       PosOutgoing2.push_back( hardProcessIn.PosOutgoing2[i]);
  80     for(int i =0; i < int(hardProcessIn.PosIntermediate.size());++i)
  81       PosIntermediate.push_back( hardProcessIn.PosIntermediate[i]);
  82   }
  83
  84   // Constructor with path to LHE file
  85   HardProcess( string LHEfile, ParticleData* particleData) {
  86     state = Event();
  87     state.init("(hard process)", particleData);
  88     translateLHEFString(LHEfile);
  89   }
  90
  91   // Constructor with core process input
  92   void initOnProcess( string process, ParticleData* particleData);
  93
  94   // Constructor with path to LHE file input
  95   void initOnLHEF( string LHEfile, ParticleData* particleData);
  96
  97   // Function to access the LHE file and read relevant information
  98   void translateLHEFString( string LHEpath);
  99
 100   // Function to translate the process string (in MG/ME notation)
 101   void translateProcessString( string process);
 102
 103   // Function to clear hard process information
 104   void clear();
 105
 106   // Function to check whether the sets of candidates Pos1, Pos2, together
 107   // with the proposed candidate iPos give an allowed hard process state
 108   bool allowCandidates(int iPos, vector<int> Pos1, vector<int> Pos2,
 109     const Event& event);
 110   // Function to identify the hard subprocess in the current event
 111   void storeCandidates( const Event& event, string process);
 112   // Function to check if the particle event[iPos] matches any of
 113   // the stored outgoing particles of the hard subprocess
 114   bool matchesAnyOutgoing(int iPos, const Event& event);
 115
 116   // Function to get the number of coloured final state partons in the
 117   // hard process
 118   int nQuarksOut();
 119   // Function to get the number of uncoloured final state particles in the
 120   // hard process
 121   int nLeptonOut();
 122   // Function to get the number of electroweak final state bosons in the
 123   // hard process
 124   int nBosonsOut();
 125
 126   // Function to get the number of coloured initial state partons in the
 127   // hard process
 128   int nQuarksIn();
 129   // Function to get the number of uncoloured initial state particles in the
 130   // hard process
 131   int nLeptonIn();
 132   // Function to report if a resonace decay was found in the 2->2 sub-process
 133   // of the  current state
 134   int hasResInCurrent();
 135   // Function to report the number of resonace decays in the 2->2 sub-process
 136   // of the  current state
 137   int nResInCurrent();
 138   // Function to report if a resonace decay was found in the 2->2 hard process
 139   bool hasResInProc();
 140   // Function to print the hard process (for debug)
 141   void list() const;
 142   // Function to print the hard process candidates in the
 143   // Matrix element state (for debug)
 144   void listCandidates() const;
 145
 146 };
 147
 148 //==========================================================================
 149
 150 // MergingHooks is base class for user input to the merging procedure.
 151
 152 class MergingHooks {
 153
 154 public:
 155
 156   // Constructor.
 157   MergingHooks() :
 158     doUserMergingSave(false),
 159     doMGMergingSave(false),
 160     doKTMergingSave(false),
 161     doPTLundMergingSave(false),
 162     doCutBasedMergingSave(false),
 163     doNL3TreeSave(false),
 164     doNL3LoopSave(false),
 165     doNL3SubtSave(false),
 166     doUNLOPSTreeSave(false),
 167     doUNLOPSLoopSave(false),
 168     doUNLOPSSubtSave(false),
 169     doUNLOPSSubtNLOSave(false),
 170     doUMEPSTreeSave(false),
 171     doUMEPSSubtSave(false),
 172     doEstimateXSection(false),
 173     doRemoveDecayProducts(false),
 174     doOrderHistoriesSave(true),
 175     doCutOnRecStateSave(false),
 176     doWClusteringSave(false),
 177     doSQCDClusteringSave(false),
 178     doIgnoreEmissionsSave(true),
 179     doIgnoreStepSave(true) {
 180       inputEvent = Event(); resonances.resize(0); infoPtr = 0;
 181       particleDataPtr = 0; partonSystemsPtr = 0;}
 182
 183   // Make History class friend to allow access to advanced switches
 184   friend class History;
 185   // Make Pythia class friend
 186   friend class Pythia;
 187   // Make PartonLevel class friend
 188   friend class PartonLevel;
 189   // Make SpaceShower class friend
 190   friend class SpaceShower;
 191   // Make TimeShower class friend
 192   friend class TimeShower;
 193   // Make Merging class friend
 194   friend class Merging;
 195
 196   //----------------------------------------------------------------------//
 197   // Functions that allow user interference
 198   //----------------------------------------------------------------------//
 199
 200   // Destructor.
 201   virtual ~MergingHooks(){}
 202   // Function encoding the functional definition of the merging scale
 203   virtual double tmsDefinition( const Event& event){ return event[0].e();}
 204   // Function to dampen weights calculated from histories with lowest
 205   // multiplicity reclustered events that do not pass the ME cuts
 206   virtual double dampenIfFailCuts( const Event& inEvent ) {
 207     // Dummy statement to avoid compiler warnings
 208     if(false) cout << inEvent[0].e();
 209     return 1.;
 210   }
 211   // Hooks to disallow states in the construction of all histories, e.g.
 212   // because jets are below the merging scale or fail the matrix element cuts
 213   // Function to allow interference in the construction of histories
 214   virtual bool canCutOnRecState() { return doCutOnRecStateSave; }
 215   // Function to check reclustered state while generating all possible
 216   // histories
 217   // Function implementing check of reclustered events while constructing
 218   // all possible histories
 219   virtual bool doCutOnRecState( const Event& event ) {
 220     // Dummy statement to avoid compiler warnings.
 221     if(false) cout << event[0].e();
 222     // Count number of final state partons.
 223     int nPartons = 0;
 224     for( int i=0; i < int(event.size()); ++i)
 225       if(  event[i].isFinal()
 226         && (event[i].isGluon() || event[i].isQuark()) )
 227         nPartons++;
 228     // For gg -> h, allow only histories with gluons in initial state
 229     if( hasEffectiveG2EW() && nPartons < 2){
 230       if(event[3].id() != 21 && event[4].id() != 21)
 231         return true;
 232     }
 233     return false;
 234   }
 235   // Function to allow not counting a trial emission.
 236   virtual bool canVetoTrialEmission() { return false;}
 237   // Function to check if trial emission should be rejected.
 238   virtual bool doVetoTrialEmission( const Event&, const Event& ) {
 239     return false; }
 240
 241   // Function to calculate the hard process matrix element.
 242   virtual double hardProcessME( const Event& inEvent ) {
 243     // Dummy statement to avoid compiler warnings.
 244     if ( false ) cout << inEvent[0].e(); return 1.; }
 245
 246   //----------------------------------------------------------------------//
 247   // Simple output functions
 248   //----------------------------------------------------------------------//
 249
 250   // Function returning the value of the merging scale.
 251   double tms() {
 252     if(doCutBasedMergingSave) return 0.;
 253     else return tmsValueSave;
 254   }
 255   // Function returning the value of the Delta R_{ij} cut for
 256   // cut based merging scale definition.
 257   double dRijMS() {
 258     return ((tmsListSave.size() == 3) ? tmsListSave[0] : 0.);
 259   }
 260   // Function returning the value of the pT_{i} cut for
 261   // cut based merging scale definition.
 262   double pTiMS() {
 263     return ((tmsListSave.size() == 3) ? tmsListSave[1] : 0.);
 264   }
 265   // Function returning the value of the pT_{i} cut for
 266   // cut based merging scale definition.
 267   double QijMS() {
 268     return ((tmsListSave.size() == 3) ? tmsListSave[2] : 0.);
 269   }
 270   // Function returning the value of the maximal number of merged jets.
 271   int nMaxJets() { return nJetMaxSave;}
 272   // Function returning the value of the maximal number of merged jets,
 273   // for which NLO corrections are available.
 274   int nMaxJetsNLO() { return nJetMaxNLOSave;}
 275   // Function to return hard process string.
 276   string getProcessString() { return processSave;}
 277   // Function to return the number of outgoing partons in the core process
 278   int nHardOutPartons(){ return hardProcess.nQuarksOut();}
 279   // Function to return the number of outgoing leptons in the core process
 280   int nHardOutLeptons(){ return hardProcess.nLeptonOut();}
 281   // Function to return the number of outgoing electroweak bosons in the core
 282   // process.
 283   int nHardOutBosons(){ return hardProcess.nBosonsOut();}
 284   // Function to return the number of incoming partons (hadrons) in the core
 285   // process.
 286   int nHardInPartons(){ return hardProcess.nQuarksIn();}
 287   // Function to return the number of incoming leptons in the core process.
 288   int nHardInLeptons(){ return hardProcess.nLeptonIn();}
 289   // Function to report the number of resonace decays in the 2->2 sub-process
 290   // of the  current state.
 291   int nResInCurrent(){ return hardProcess.nResInCurrent();}
 292   // Function to determine if user defined merging should be applied.
 293   bool doUserMerging(){ return doUserMergingSave;}
 294   // Function to determine if automated MG/ME merging should be applied.
 295   bool doMGMerging() { return doMGMergingSave;}
 296   // Function to determine if KT merging should be applied.
 297   bool doKTMerging() { return doKTMergingSave;}
 298   // Function to determine if PTLund merging should be applied.
 299   bool doPTLundMerging() { return doPTLundMergingSave;}
 300   // Function to determine if cut based merging should be applied.
 301   bool doCutBasedMerging() { return doCutBasedMergingSave;}
 302   bool doCKKWLMerging() { return (doUserMergingSave || doMGMergingSave
 303     || doKTMergingSave || doPTLundMergingSave || doCutBasedMergingSave); }
 304   // Functions to determine if and which part of  UMEPS merging
 305   // should be applied
 306   bool doUMEPSTree() { return doUMEPSTreeSave;}
 307   bool doUMEPSSubt() { return doUMEPSSubtSave;}
 308   bool doUMEPSMerging() { return (doUMEPSTreeSave || doUMEPSSubtSave);}
 309   // Functions to determine if and which part of  NL3 merging
 310   // should be applied
 311   bool doNL3Tree() { return doNL3TreeSave;}
 312   bool doNL3Loop() { return doNL3LoopSave;}
 313   bool doNL3Subt() { return doNL3SubtSave;}
 314   bool doNL3Merging() { return (doNL3TreeSave || doNL3LoopSave
 315                              || doNL3SubtSave); }
 316   // Functions to determine if and which part of UNLOPS merging
 317   // should be applied
 318   bool doUNLOPSTree() { return doUNLOPSTreeSave;}
 319   bool doUNLOPSLoop() { return doUNLOPSLoopSave;}
 320   bool doUNLOPSSubt() { return doUNLOPSSubtSave;}
 321   bool doUNLOPSSubtNLO() { return doUNLOPSSubtNLOSave;}
 322   bool doUNLOPSMerging() { return (doUNLOPSTreeSave || doUNLOPSLoopSave
 323                              || doUNLOPSSubtSave || doUNLOPSSubtNLOSave); }
 324   // Return the number clustering steps that have actually been done.
 325   int nRecluster() { return nReclusterSave;}
 326
 327   //----------------------------------------------------------------------//
 328   // Output functions to analyse/prepare event for merging
 329   //----------------------------------------------------------------------//
 330
 331   // Function to check if event contains an emission not present in the hard
 332   // process.
 333   bool isFirstEmission(const Event& event);
 334
 335   // Function to allow effective gg -> EW boson couplings.
 336   bool hasEffectiveG2EW() {
 337     if ( getProcessString().compare("pp>h") == 0 ) return true;
 338     return false; }
 339
 340   // Return event stripped from decay products.
 341   Event bareEvent( const Event& inputEventIn, bool storeInputEvent );
 342   // Write event with decay products attached to argument.
 343   bool reattachResonanceDecays( Event& process );
 344
 345   // Check if particle at position iPos is part of the hard sub-system
 346   bool isInHard( int iPos, const Event& event);
 347
 348   // Function to return the number of clustering steps for the current event
 349   int getNumberOfClusteringSteps(const Event& event);
 350
 351   //----------------------------------------------------------------------//
 352   // Functions to steer contruction of histories
 353   //----------------------------------------------------------------------//
 354
 355   // Function to force preferred picking of ordered histories. By default,
 356   // unordered histories will only be considered if no ordered histories
 357   // were found.
 358   void orderHistories( bool doOrderHistoriesIn) {
 359     doOrderHistoriesSave = doOrderHistoriesIn; }
 360   // Function to force cut on reconstructed states internally, as needed
 361   // for gg -> Higgs to ensure that e.g. uu~ -> Higgs is not constructed.
 362   void allowCutOnRecState( bool doCutOnRecStateIn) {
 363     doCutOnRecStateSave = doCutOnRecStateIn; }
 364
 365   // Function to allow final state clusterings of W-bosons
 366   void doWClustering( bool doWClusteringIn ) {
 367     doWClusteringSave = doWClusteringIn; }
 368
 369   //----------------------------------------------------------------------//
 370   // Functions used as default merging scales
 371   //----------------------------------------------------------------------//
 372
 373   // Function to return the value of the merging scale function in the
 374   // current event.
 375   double tmsNow( const Event& event );
 376   // Find the minimal Lund pT between coloured partons in the event
 377   double rhoms( const Event& event, bool withColour);
 378   // Function to calculate the minimal kT in the event
 379   double kTms(const Event & event);
 380   // Find the if the event passes the Delta R_{ij}, pT_{i} and Q_{ij} cuts on
 381   // the matrix element, as needed for cut-based merging scale definition
 382   double cutbasedms( const Event& event );
 383
 384 protected:
 385
 386   //----------------------------------------------------------------------//
 387   // The members, switches etc.
 388   //----------------------------------------------------------------------//
 389
 390   // Helper class doing all the core process book-keeping
 391   HardProcess hardProcess;
 392
 393   // Pointer to various information on the generation.
 394   Info*          infoPtr;
 395
 396   // Pointer to the particle data table.
 397   ParticleData*  particleDataPtr;
 398
 399   // Pointer to the particle systems.
 400   PartonSystems* partonSystemsPtr;
 401
 402   // AlphaS objects for alphaS reweighting use
 403   AlphaStrong AlphaS_FSRSave;
 404   AlphaStrong AlphaS_ISRSave;
 405   AlphaEM AlphaEM_FSRSave;
 406
 407   // Saved path to LHE file for more automated merging
 408   string lheInputFile;
 409
 410   // Flags for merging procedure definition.
 411   bool   doUserMergingSave, doMGMergingSave, doKTMergingSave,
 412          doPTLundMergingSave, doCutBasedMergingSave,
 413          includeMassiveSave, enforceStrongOrderingSave, orderInRapiditySave,
 414          pickByFullPSave, pickByPoPT2Save, includeRedundantSave,
 415          pickBySumPTSave, allowColourShufflingSave, resetHardQRenSave,
 416          resetHardQFacSave;
 417   int    unorderedScalePrescipSave, unorderedASscalePrescipSave,
 418          unorderedPDFscalePrescipSave, incompleteScalePrescipSave,
 419          ktTypeSave, nReclusterSave;
 420   double scaleSeparationFactorSave, nonJoinedNormSave,
 421          fsrInRecNormSave, herwigAcollFSRSave, herwigAcollISRSave,
 422          pT0ISRSave, pTcutSave;
 423   bool   doNL3TreeSave, doNL3LoopSave, doNL3SubtSave;
 424   bool   doUNLOPSTreeSave, doUNLOPSLoopSave, doUNLOPSSubtSave,
 425          doUNLOPSSubtNLOSave;
 426   bool   doUMEPSTreeSave, doUMEPSSubtSave;
 427
 428   // Flag to only do phase space cut, rejecting events below the tms cut.
 429   bool   doEstimateXSection;
 430
 431   // Save input event in case decay products need to be detached.
 432   Event inputEvent;
 433   vector< pair<int,int> > resonances;
 434   bool doRemoveDecayProducts;
 435
 436   // Starting scale for attaching MPI.
 437   double muMISave;
 438
 439   // Precalculated K-factors.
 440   double kFactor0jSave;
 441   double kFactor1jSave;
 442   double kFactor2jSave;
 443
 444   // Saved members.
 445   double tmsValueSave;
 446   int nJetMaxSave;
 447   int nJetMaxNLOSave;
 448   string processSave;
 449
 450   // List of cut values to used to define a merging scale. Ordering:
 451   // 0: DeltaR_{jet_i,jet_j,min}
 452   // 1: p_{T,jet_i,min}
 453   // 2: Q_{jet_i,jet_j,min}
 454   vector<double> tmsListSave;
 455
 456   // INTERNAL Hooks to allow construction of all histories,
 457   // including un-ordered ones
 458   bool doOrderHistoriesSave;
 459
 460   // INTERNAL Hooks to disallow states in the construction of all histories,
 461   // e.g. because jets are below the merging scale, of to avoid the
 462   // construction of uu~ -> Higgs histories.
 463   bool doCutOnRecStateSave;
 464
 465   // INTERNAL Hooks to allow clustering W bosons.
 466   bool doWClusteringSave, doSQCDClusteringSave;
 467
 468   // Store / get first scale in PDF's that Pythia should have used
 469   double muFSave;
 470   double muRSave;
 471
 472   // Store / get factorisation scale used in matrix element calculation.
 473   double muFinMESave;
 474   double muRinMESave;
 475
 476   // Flag to indicate trial shower usage.
 477   bool doIgnoreEmissionsSave;
 478   // Flag to indicate if events should be vetoed.
 479   bool doIgnoreStepSave;
 480   // Stored weights in case veot needs to be revoked
 481   double pTsave;
 482   double weightCKKWL1Save, weightCKKWL2Save;
 483   int nMinMPISave;
 484   // Save CKKW-L weight / O(\alpha_s) weight.
 485   double weightCKKWLSave, weightFIRSTSave;
 486
 487   //----------------------------------------------------------------------//
 488   // Generic setup functions
 489   //----------------------------------------------------------------------//
 490
 491   // Functions for internal use inside Pythia source code
 492   // Initialize.
 493   void init( Settings* settings, Info* infoPtrIn,
 494     ParticleData* particleDataPtrIn, PartonSystems* partonSystemsPtrIn,
 495     ostream& os = cout);
 496
 497   // Function storing candidates for the hard process in the current event
 498   // Needed in order not to cluster members of the core process
 499   void storeHardProcessCandidates(const Event& event){
 500     hardProcess.storeCandidates(event,getProcessString());
 501   }
 502   // Function to set the path to the LHE file, so that more automated merging
 503   // can be used.
 504   // Remove "_1.lhe" suffix from LHE file name.
 505   // This is done so that the HarsProcess class can access both the +0 and +1
 506   // LHE files to find both the merging scale and the core process string
 507   // Store.
 508   void setLHEInputFile( string lheFile) {
 509     lheInputFile = lheFile.substr(0,lheFile.size()-6); }
 510
 511   //----------------------------------------------------------------------//
 512   // Functions for output of class members.
 513   //----------------------------------------------------------------------//
 514
 515   // Return AlphaStrong objects
 516   AlphaStrong* AlphaS_FSR() { return &AlphaS_FSRSave;}
 517   AlphaStrong* AlphaS_ISR() { return &AlphaS_ISRSave;}
 518   AlphaEM* AlphaEM_FSR() { return &AlphaEM_FSRSave;}
 519
 520   // Functions to return advanced merging switches
 521   // Include masses in definition of evolution pT and splitting kernels
 522   bool includeMassive() { return includeMassiveSave;}
 523   // Prefer strongly ordered histories
 524   bool enforceStrongOrdering() { return enforceStrongOrderingSave;}
 525   // Prefer histories ordered in rapidity and evolution pT
 526   bool orderInRapidity() { return orderInRapiditySave;}
 527   // Pick history probabilistically by full (correct) splitting probabilities
 528   bool pickByFull() { return pickByFullPSave;}
 529   // Pick history probabilistically, with easier form of probabilities
 530   bool pickByPoPT2() { return pickByPoPT2Save;}
 531   // Include redundant terms (e.g. PDF ratios) in the splitting probabilities
 532   bool includeRedundant(){ return includeRedundantSave;}
 533   // Pick by winner-takes-it-all, with minimum sum of scalar evolution pT
 534   bool pickBySumPT(){ return pickBySumPTSave;}
 535
 536   // Prescription for combined scale of unordered emissions
 537   // 0 : use larger scale
 538   // 1 : use smaller scale
 539   int unorderedScalePrescip() { return unorderedScalePrescipSave;}
 540   // Prescription for combined scale used in alpha_s for unordered emissions
 541   // 0 : use combined emission scale in alpha_s weight for both (!) splittings
 542   // 1 : use original reclustered scales of each emission in alpha_s weight
 543   int unorderedASscalePrescip() { return unorderedASscalePrescipSave;}
 544   // Prescription for combined scale used in PDF ratios for unordered
 545   // emissions
 546   // 0 : use combined emission scale in PDFs for both (!) splittings
 547   // 1 : use original reclustered scales of each emission in PDF ratiost
 548   int unorderedPDFscalePrescip() { return unorderedPDFscalePrescipSave;}
 549   // Prescription for starting scale of incomplete histories
 550   // 0: use factorization scale
 551   // 1: use sHat
 552   // 2: use s
 553   int incompleteScalePrescip() { return incompleteScalePrescipSave;}
 554
 555   // Allow swapping one colour index while reclustering
 556   bool allowColourShuffling() { return allowColourShufflingSave;}
 557
 558   // Allow use of dynamical renormalisation scale of the core 2-> 2 process.
 559   bool resetHardQRen() { return resetHardQRenSave; }
 560   // Allow use of dynamical factorisation scale of the core 2-> 2 process.
 561   bool resetHardQFac() { return resetHardQFacSave; }
 562
 563   // Factor by which two scales should differ to be classified strongly
 564   // ordered.
 565   double scaleSeparationFactor() { return scaleSeparationFactorSave;}
 566   // Absolute normalization of splitting probability for non-joined
 567   // evolution.
 568   double nonJoinedNorm() { return nonJoinedNormSave;}
 569   // Absolute normalization of splitting probability for final state
 570   // splittings with initial state recoiler
 571   double fsrInRecNorm() { return fsrInRecNormSave;}
 572   // Factor multiplying scalar evolution pT for FSR splitting, when picking
 573   // history by minimum scalar pT (see Jonathan Tully's thesis)
 574   double herwigAcollFSR() { return herwigAcollFSRSave;}
 575   // Factor multiplying scalar evolution pT for ISR splitting, when picking
 576   // history by minimum scalar pT (see Jonathan Tully's thesis)
 577   double herwigAcollISR() { return herwigAcollISRSave;}
 578   // ISR regularisation scale
 579   double pT0ISR() { return pT0ISRSave;}
 580   // Shower cut-off scale
 581   double pTcut() { return pTcutSave;}
 582
 583   // MI starting scale.
 584   void muMI( double mu) { muMISave = mu; }
 585   double muMI() { return muMISave;}
 586
 587   // Full k-Factor for current event
 588   double kFactor(int njet = 0) {
 589     return (njet == 0) ? kFactor0jSave
 590           :(njet == 1) ? kFactor1jSave
 591           : kFactor2jSave;
 592   }
 593   // O(\alhpa_s)-term of the k-Factor for current event
 594   double k1Factor( int njet = 0) {
 595     return (kFactor(njet) - 1)/infoPtr->alphaS();
 596   }
 597
 598   // Function to return if construction of histories is biased towards ordered
 599   // histories.
 600   bool orderHistories() { return doOrderHistoriesSave;}
 601
 602   // INTERNAL Hooks to disallow states in the construction of all histories,
 603   // e.g. because jets are below the merging scale, of to avoid the
 604   // construction of uu~ -> Higgs histories.
 605   bool allowCutOnRecState() { return doCutOnRecStateSave;}
 606
 607   // INTERNAL Hooks to allow clustering W bosons.
 608   bool doWClustering() { return doWClusteringSave;}
 609   // INTERNAL Hooks to allow clustering clustering of gluons to squarks.
 610   bool doSQCDClustering() { return doSQCDClusteringSave;}
 611
 612   // Store / get first scale in PDF's that Pythia should have used
 613   double muF() { return (muFSave > 0.) ? muFSave : infoPtr->QFac();}
 614   double muR() { return (muRSave > 0.) ? muRSave : infoPtr->QRen();}
 615   // Store / get factorisation scale used in matrix element calculation.
 616   double muFinME() { return (muFinMESave > 0.) ? muFinMESave
 617                        : infoPtr->QFac();}
 618   double muRinME() { return (muRinMESave > 0.) ? muRinMESave
 619                        : infoPtr->QRen();}
 620
 621   //----------------------------------------------------------------------//
 622   // Functions to steer shower evolution
 623   //----------------------------------------------------------------------//
 624
 625   // Flag to indicate trial shower usage.
 626   void doIgnoreEmissions( bool doIgnoreIn ) {
 627     doIgnoreEmissionsSave = doIgnoreIn;
 628   }
 629   // Function to allow not counting a trial emission.
 630   bool canVetoEmission() { return !doIgnoreEmissionsSave; }
 631   // Function to check if emission should be rejected.
 632   bool doVetoEmission( const Event& );
 633
 634   // Flag to indicate if events should be vetoed.
 635   void doIgnoreStep( bool doIgnoreIn ) { doIgnoreStepSave = doIgnoreIn; }
 636   // Function to allow event veto.
 637   bool canVetoStep() { return !doIgnoreStepSave; }
 638   // Function to check event veto.
 639   bool doVetoStep( const Event& process, const Event& event,
 640     bool doResonance = false );
 641
 642   // Stored weights in case veot needs to be revoked
 643   void storeWeights( double weight ){ weightCKKWL1Save = weightCKKWL2Save
 644      = weight; }
 645
 646   // Set starting scales
 647   bool setShowerStartingScales( bool isTrial, bool doMergeFirstEmm,
 648     double& pTscaleIn, const Event& event,
 649     double& pTmaxFSRIn, bool& limitPTmaxFSRin,
 650     double& pTmaxISRIn, bool& limitPTmaxISRin,
 651     double& pTmaxMPIIn, bool& limitPTmaxMPIin );
 652   void nMinMPI( int nMinMPIIn ) { nMinMPISave = nMinMPIIn; }
 653   int nMinMPI() { return nMinMPISave;}
 654
 655   //----------------------------------------------------------------------//
 656   // Functions for internal merging scale definions
 657   //----------------------------------------------------------------------//
 658
 659   // Function to calculate the kT separation between two particles
 660   double kTdurham(const Particle& RadAfterBranch,
 661     const Particle& EmtAfterBranch, int Type, double D );
 662   // Function to compute "pythia pT separation" from Particle input
 663   double rhoPythia(const Particle& RadAfterBranch,
 664     const Particle& EmtAfterBranch, const Particle& RecAfterBranch,
 665     int ShowerType);
 666   // Function to find a colour (anticolour) index in the input event,
 667   // used to find colour-connected recoilers
 668   int findColour(int col, int iExclude1, int iExclude2,
 669     const Event& event, int type, bool isHardIn);
 670   // Function to check if the properties of the input particle should be
 671   // checked against the cut-based merging scale defintion.
 672   bool checkAgainstCut( const Particle& particle);
 673   // Function to compute Delta R separation from 4-vector input
 674   double deltaRij(Vec4 jet1, Vec4 jet2);
 675
 676   //----------------------------------------------------------------------//
 677   // Functions for weight management
 678   //----------------------------------------------------------------------//
 679
 680   // Function to get the CKKW-L weight for the current event
 681   double getWeightNLO() { return (weightCKKWLSave - weightFIRSTSave);}
 682   // Return CKKW-L weight.
 683   double getWeightCKKWL() { return weightCKKWLSave; }
 684   // Return O(\alpha_s) weight.
 685   double getWeightFIRST() { return weightFIRSTSave; }
 686   // Set CKKW-L weight.
 687   void setWeightCKKWL( double weightIn){
 688     weightCKKWLSave = weightIn;
 689     infoPtr->setWeightCKKWL(weightIn); }
 690   // Set O(\alpha_s) weight.
 691   void setWeightFIRST( double weightIn){
 692     weightFIRSTSave = weightIn;
 693     infoPtr->setWeightFIRST(weightIn); }
 694
 695 };
 696
 697 //==========================================================================
 698
 699 } // end namespace Pythia8
 700
 701 #endif // Pythia8_MergingHooks_H