Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-09-13 09:05:22

0001 // MergingHooks.h is a part of the PYTHIA event generator.
0002 // Copyright (C) 2025 Torbjorn Sjostrand.
0003 // PYTHIA is licenced under the GNU GPL v2 or later, see COPYING for details.
0004 // Please respect the MCnet Guidelines, see GUIDELINES for details.
0005 
0006 // This file is written by Stefan Prestel.
0007 // Header file to allow user access to program at different stages.
0008 // HardProcess: Container class for the hard process to be merged. Holds the
0009 //              bookkeeping of particles not be be reclustered
0010 // MergingHooks: Steering class for matrix element merging. Some functions can
0011 //               be redefined in a derived class to have access to the merging
0012 
0013 #ifndef Pythia8_MergingHooks_H
0014 #define Pythia8_MergingHooks_H
0015 
0016 #include "Pythia8/Basics.h"
0017 #include "Pythia8/BeamParticle.h"
0018 #include "Pythia8/Event.h"
0019 #include "Pythia8/Info.h"
0020 #include "Pythia8/ParticleData.h"
0021 #include "Pythia8/PartonSystems.h"
0022 #include "Pythia8/PhysicsBase.h"
0023 #include "Pythia8/PythiaStdlib.h"
0024 #include "Pythia8/Settings.h"
0025 
0026 
0027 namespace Pythia8 {
0028 
0029 class PartonLevel;
0030 
0031 //==========================================================================
0032 
0033 // Declaration of hard process class
0034 // This class holds information on the desired hard 2->2 process
0035 // for the merging.
0036 // This class is a container class for History class use.
0037 
0038 class HardProcess {
0039 
0040 public:
0041 
0042    // Flavour of the first incoming particle
0043   int hardIncoming1;
0044   // Flavour of the second incoming particle
0045   int hardIncoming2;
0046   // Flavours of the outgoing particles
0047   vector<int> hardOutgoing1;
0048   vector<int> hardOutgoing2;
0049   // Flavour of intermediate bosons in the hard 2->2
0050   vector<int> hardIntermediate;
0051 
0052   // Current reference event
0053   Event state;
0054   // Potential positions of outgoing particles in reference event
0055   vector<int> PosOutgoing1;
0056   vector<int> PosOutgoing2;
0057   // Potential positions of intermediate bosons in reference event
0058   vector<int> PosIntermediate;
0059 
0060   // Information on merging scale read from LHE file
0061   double tms;
0062 
0063   // Default constructor
0064   HardProcess() : hardIncoming1(), hardIncoming2(), tms(){}
0065   // Default destructor
0066   virtual ~HardProcess(){}
0067 
0068   // Copy constructor
0069   HardProcess( const HardProcess& hardProcessIn )
0070     : state(hardProcessIn.state),
0071       tms(hardProcessIn.tms) {
0072     hardIncoming1 = hardProcessIn.hardIncoming1;
0073     hardIncoming2 = hardProcessIn.hardIncoming2;
0074     for(int i =0; i < int(hardProcessIn.hardOutgoing1.size());++i)
0075       hardOutgoing1.push_back( hardProcessIn.hardOutgoing1[i]);
0076     for(int i =0; i < int(hardProcessIn.hardOutgoing2.size());++i)
0077       hardOutgoing2.push_back( hardProcessIn.hardOutgoing2[i]);
0078     for(int i =0; i < int(hardProcessIn.hardIntermediate.size());++i)
0079       hardIntermediate.push_back( hardProcessIn.hardIntermediate[i]);
0080     for(int i =0; i < int(hardProcessIn.PosOutgoing1.size());++i)
0081       PosOutgoing1.push_back( hardProcessIn.PosOutgoing1[i]);
0082     for(int i =0; i < int(hardProcessIn.PosOutgoing2.size());++i)
0083       PosOutgoing2.push_back( hardProcessIn.PosOutgoing2[i]);
0084     for(int i =0; i < int(hardProcessIn.PosIntermediate.size());++i)
0085       PosIntermediate.push_back( hardProcessIn.PosIntermediate[i]);
0086   }
0087 
0088   // Constructor with path to LHE file
0089   HardProcess( string LHEfile, ParticleData* particleData) : hardIncoming1(),
0090     hardIncoming2(), tms() {
0091     state = Event();
0092     state.init("(hard process)", particleData);
0093     translateLHEFString(LHEfile);
0094   }
0095 
0096   // Constructor with core process input
0097   virtual void initOnProcess( string process, ParticleData* particleData);
0098 
0099   // Constructor with path to LHE file input
0100   void initOnLHEF( string LHEfile, ParticleData* particleData);
0101 
0102   // Function to access the LHE file and read relevant information
0103   void translateLHEFString( string LHEpath);
0104 
0105   // Function to translate the process string (in MG/ME notation)
0106   virtual void translateProcessString( string process);
0107 
0108   // Function to clear hard process information
0109   void clear();
0110 
0111   // Function to check whether the sets of candidates Pos1, Pos2, together
0112   // with the proposed candidate iPos give an allowed hard process state
0113   virtual bool allowCandidates(int iPos, vector<int> Pos1, vector<int> Pos2,
0114     const Event& event);
0115   // Function to identify the hard subprocess in the current event
0116   virtual void storeCandidates( const Event& event, string process);
0117   // Function to check if the particle event[iPos] matches any of
0118   // the stored outgoing particles of the hard subprocess
0119   virtual bool matchesAnyOutgoing(int iPos, const Event& event);
0120   // Function to check if instead of the particle event[iCandidate], another
0121   // particle could serve as part of the hard process. Assumes that iCandidate
0122   // is already stored as part of the hard process.
0123   virtual bool findOtherCandidates(int iPos, const Event& event,
0124     bool doReplace);
0125   // Function to exchange a stored hard process candidate with another choice.
0126   virtual bool exchangeCandidates( vector<int> candidates1,
0127     vector<int> candidates2,
0128     unordered_map<int,int> further1, unordered_map<int,int> further2);
0129 
0130   // Function to get the number of coloured final state partons in the
0131   // hard process
0132   int nQuarksOut();
0133   // Function to get the number of uncoloured final state particles in the
0134   // hard process
0135   int nLeptonOut();
0136   // Function to get the number of electroweak final state bosons in the
0137   // hard process
0138   int nBosonsOut();
0139 
0140   // Function to get the number of coloured initial state partons in the
0141   // hard process
0142   int nQuarksIn();
0143   // Function to get the number of uncoloured initial state particles in the
0144   // hard process
0145   int nLeptonIn();
0146   // Function to report if a resonace decay was found in the 2->2 sub-process
0147   // of the  current state
0148   int hasResInCurrent();
0149   // Function to report the number of resonace decays in the 2->2 sub-process
0150   // of the  current state
0151   int nResInCurrent();
0152   // Function to report if a resonace decay was found in the 2->2 hard process
0153   bool hasResInProc();
0154   // Function to print the hard process (for debug)
0155   void list() const;
0156   // Function to print the hard process candidates in the
0157   // Matrix element state (for debug)
0158   void listCandidates() const;
0159 
0160 };
0161 
0162 //==========================================================================
0163 
0164 // MergingHooks is base class for user input to the merging procedure.
0165 
0166 class MergingHooks : public PhysicsBase {
0167 
0168 public:
0169 
0170   // Constructor.
0171   MergingHooks() : useShowerPluginSave(), showers(),
0172     doUserMergingSave(false),
0173     doMGMergingSave(false),
0174     doKTMergingSave(false),
0175     doPTLundMergingSave(false),
0176     doCutBasedMergingSave(false),
0177     doDynamicMergingSave(false), includeMassiveSave(),
0178     enforceStrongOrderingSave(),
0179     orderInRapiditySave(), pickByFullPSave(), pickByPoPT2Save(),
0180     includeRedundantSave(), pickBySumPTSave(), allowColourShufflingSave(),
0181     resetHardQRenSave(), resetHardQFacSave(), unorderedScalePrescipSave(),
0182     unorderedASscalePrescipSave(), unorderedPDFscalePrescipSave(),
0183     incompleteScalePrescipSave(), ktTypeSave(), nReclusterSave(),
0184     nQuarksMergeSave(), nRequestedSave(), scaleSeparationFactorSave(),
0185     nonJoinedNormSave(), fsrInRecNormSave(), herwigAcollFSRSave(),
0186     herwigAcollISRSave(), pT0ISRSave(), pTcutSave(),
0187     doNL3TreeSave(false),
0188     doNL3LoopSave(false),
0189     doNL3SubtSave(false),
0190     doUNLOPSTreeSave(false),
0191     doUNLOPSLoopSave(false),
0192     doUNLOPSSubtSave(false),
0193     doUNLOPSSubtNLOSave(false),
0194     doUMEPSTreeSave(false),
0195     doUMEPSSubtSave(false),
0196     doEstimateXSection(false),
0197     doRuntimeAMCATNLOInterfaceSave(false),
0198     applyVeto(),
0199     doRemoveDecayProducts(false), muMISave(), kFactor0jSave(), kFactor1jSave(),
0200     kFactor2jSave(), tmsValueSave(), tmsValueNow(),
0201     DparameterSave(), SparameterSave(),
0202     nJetMaxSave(), nJetMaxNLOSave(),
0203     doOrderHistoriesSave(true),
0204     doCutOnRecStateSave(false),
0205     doWeakClusteringSave(false),
0206     doSQCDClusteringSave(false), muFSave(), muRSave(), muFinMESave(),
0207     muRinMESave(),
0208     doIgnoreEmissionsSave(true),
0209     doIgnoreStepSave(true), pTsave(), weightCKKWL1Save(), weightCKKWL2Save(),
0210     nMinMPISave(), weightCKKWLSave(), weightFIRSTSave(),
0211     doVariations(false), nWgts(0), nJetMaxLocal(), nJetMaxNLOLocal(),
0212     hasJetMaxLocal(false),
0213     includeWGTinXSECSave(false), nHardNowSave(), nJetNowSave(),
0214     tmsHardNowSave(), tmsNowSave() {
0215       inputEvent = Event(); resonances.resize(0);
0216       useOwnHardProcess = false; hardProcess = 0; stopScaleSave= 0.0;
0217       nVetoedInMainShower = 0;}
0218 
0219   // Make History class friend to allow access to advanced switches
0220   friend class History;
0221   // Make Pythia class friend
0222   friend class Pythia;
0223   // Make PartonLevel class friend
0224   friend class PartonLevel;
0225   // Make SpaceShower class friend
0226   friend class SpaceShower;
0227   // Make TimeShower class friend
0228   friend class TimeShower;
0229   // Make Merging class friend
0230   friend class Merging;
0231 
0232   // Destructor.
0233   virtual ~MergingHooks();
0234   // Function encoding the functional definition of the merging scale
0235   virtual double tmsDefinition( const Event& event){ return event[0].e();}
0236   // Function to dampen weights calculated from histories with lowest
0237   // multiplicity reclustered events that do not pass the ME cuts
0238   virtual double dampenIfFailCuts( const Event& inEvent ) {
0239     // Dummy statement to avoid compiler warnings
0240     if(false) cout << inEvent[0].e();
0241     return 1.;
0242   }
0243   // Hooks to disallow states in the construction of all histories, e.g.
0244   // because jets are below the merging scale or fail the matrix element cuts
0245   // Function to allow interference in the construction of histories
0246   virtual bool canCutOnRecState() { return doCutOnRecStateSave; }
0247   // Function to check reclustered state while generating all possible
0248   // histories
0249   // Function implementing check of reclustered events while constructing
0250   // all possible histories
0251   virtual bool doCutOnRecState( const Event& event ) {
0252     // Dummy statement to avoid compiler warnings.
0253     if(false) cout << event[0].e();
0254     // Count number of final state partons.
0255     int nPartons = 0;
0256     for( int i=0; i < int(event.size()); ++i)
0257       if(  event[i].isFinal()
0258         && (event[i].isGluon() || event[i].isQuark()) )
0259         nPartons++;
0260     // For gg -> h, allow only histories with gluons in initial state
0261     if( hasEffectiveG2EW() && nPartons < 2){
0262       if(event[3].id() != 21 && event[4].id() != 21)
0263         return true;
0264     }
0265     return false;
0266   }
0267   // Function to allow not counting a trial emission.
0268   virtual bool canVetoTrialEmission() { return false;}
0269   // Function to check if trial emission should be rejected.
0270   virtual bool doVetoTrialEmission( const Event&, const Event& ) {
0271     return false; }
0272 
0273   // Function to calculate the hard process matrix element.
0274   virtual double hardProcessME( const Event& inEvent ) {
0275     // Dummy statement to avoid compiler warnings.
0276     if ( false ) cout << inEvent[0].e();
0277     return 1.; }
0278 
0279   // Functions for internal use inside Pythia source code
0280   // Initialize.
0281   virtual void init();
0282 
0283   // Function returning the value of the merging scale.
0284   double tms() {
0285     if(doCutBasedMergingSave) return 0.;
0286     //else return tmsValueSave;
0287     else return tmsValueNow;
0288   }
0289   double tmsCut() {
0290     if (doCutBasedMergingSave) return 0.;
0291     if (doDynamicMergingSave) {
0292       // Fetch fixed merging cut and factorisation scale.
0293       const double QbarCut = tmsValueSave;
0294       const double Q2      = infoPtr->Q2DIS();
0295       // Calculate dynamic merging scale, cf. eq. (2) in arXiv:0912.3715.
0296       return QbarCut/sqrt(1. + pow2(QbarCut/SparameterSave)/Q2);
0297     }
0298     return tmsValueSave;
0299   }
0300   void tms( double tmsIn ) { tmsValueNow = tmsIn; }
0301 
0302   // Function returning the value of the Delta R_{ij} cut for
0303   // cut based merging scale definition.
0304   double dRijMS() {
0305     return ((tmsListSave.size() == 3) ? tmsListSave[0] : 0.);
0306   }
0307   // Function returning the value of the pT_{i} cut for
0308   // cut based merging scale definition.
0309   double pTiMS() {
0310     return ((tmsListSave.size() == 3) ? tmsListSave[1] : 0.);
0311   }
0312   // Function returning the value of the pT_{i} cut for
0313   // cut based merging scale definition.
0314   double QijMS() {
0315     return ((tmsListSave.size() == 3) ? tmsListSave[2] : 0.);
0316   }
0317   // Function returning the value of the maximal number of merged jets.
0318   int nMaxJets() { return (hasJetMaxLocal) ? nJetMaxLocal : nJetMaxSave;}
0319   // Function returning the value of the maximal number of merged jets,
0320   // for which NLO corrections are available.
0321   int nMaxJetsNLO()
0322     { return (hasJetMaxLocal) ? nJetMaxNLOLocal : nJetMaxNLOSave;}
0323   // Function to return hard process string.
0324   string getProcessString() { return processNow;}
0325   // Function to return the number of outgoing partons in the core process
0326   int nHardOutPartons(){ return hardProcess->nQuarksOut();}
0327   // Function to return the number of outgoing leptons in the core process
0328   int nHardOutLeptons(){ return hardProcess->nLeptonOut();}
0329   // Function to return the number of outgoing electroweak bosons in the core
0330   // process.
0331   int nHardOutBosons(){ return hardProcess->nBosonsOut();}
0332   // Function to return the number of incoming partons (hadrons) in the core
0333   // process.
0334   int nHardInPartons(){ return hardProcess->nQuarksIn();}
0335   // Function to return the number of incoming leptons in the core process.
0336   int nHardInLeptons(){ return hardProcess->nLeptonIn();}
0337   // Function to report the number of resonace decays in the 2->2 sub-process
0338   // of the  current state.
0339   int nResInCurrent(){ return hardProcess->nResInCurrent();}
0340   // Function to determine if user defined merging should be applied.
0341   bool doUserMerging(){ return doUserMergingSave;}
0342   // Function to determine if automated MG/ME merging should be applied.
0343   bool doMGMerging() { return doMGMergingSave;}
0344   // Function to determine if KT merging should be applied.
0345   bool doKTMerging() { return doKTMergingSave;}
0346   // Function to determine if PTLund merging should be applied.
0347   bool doPTLundMerging() { return doPTLundMergingSave;}
0348   // Function to determine if cut based merging should be applied.
0349   bool doCutBasedMerging() { return doCutBasedMergingSave;}
0350   // Function to determine if merging with dynamic scale should be applied.
0351   bool doDynamicMerging() { return doDynamicMergingSave;}
0352   bool doCKKWLMerging() { return (doUserMergingSave || doMGMergingSave
0353       || doKTMergingSave || doPTLundMergingSave || doCutBasedMergingSave
0354       || doDynamicMergingSave); }
0355   // Functions to determine if and which part of  UMEPS merging
0356   // should be applied
0357   bool doUMEPSTree() { return doUMEPSTreeSave;}
0358   bool doUMEPSSubt() { return doUMEPSSubtSave;}
0359   bool doUMEPSMerging() { return (doUMEPSTreeSave || doUMEPSSubtSave);}
0360   // Functions to determine if and which part of  NL3 merging
0361   // should be applied
0362   bool doNL3Tree() { return doNL3TreeSave;}
0363   bool doNL3Loop() { return doNL3LoopSave;}
0364   bool doNL3Subt() { return doNL3SubtSave;}
0365   bool doNL3Merging() { return (doNL3TreeSave || doNL3LoopSave
0366                              || doNL3SubtSave); }
0367   // Functions to determine if and which part of UNLOPS merging
0368   // should be applied
0369   bool doUNLOPSTree() { return doUNLOPSTreeSave;}
0370   bool doUNLOPSLoop() { return doUNLOPSLoopSave;}
0371   bool doUNLOPSSubt() { return doUNLOPSSubtSave;}
0372   bool doUNLOPSSubtNLO() { return doUNLOPSSubtNLOSave;}
0373   bool doUNLOPSMerging() { return (doUNLOPSTreeSave || doUNLOPSLoopSave
0374                              || doUNLOPSSubtSave || doUNLOPSSubtNLOSave); }
0375 
0376   // Function to determine if we have a runtime interface to aMC@NLO.
0377   bool doRuntimeAMCATNLOInterface() { return doRuntimeAMCATNLOInterfaceSave;}
0378 
0379   // Return the number clustering steps that have actually been done.
0380   int nRecluster() { return nReclusterSave;}
0381 
0382   // Return number of requested additional jets on top of the Born process.
0383   int nRequested() { return nRequestedSave;}
0384 
0385   //----------------------------------------------------------------------//
0386   // Output functions to analyse/prepare event for merging
0387   //----------------------------------------------------------------------//
0388 
0389   // Function to check if event contains an emission not present in the hard
0390   // process.
0391   bool isFirstEmission(const Event& event);
0392 
0393   // Function to allow effective gg -> EW boson couplings.
0394   bool hasEffectiveG2EW() {
0395     if ( getProcessString().compare("pp>h") == 0 ) return true;
0396     return false; }
0397 
0398   // Function to allow effective gg -> EW boson couplings.
0399   bool allowEffectiveVertex( vector<int> in, vector<int> out) {
0400     if ( getProcessString().compare("ta+ta->jj") == 0
0401       || getProcessString().compare("ta-ta+>jj") == 0 ) {
0402       int nInFermions(0), nOutFermions(0);
0403       for (int i=0; i < int(in.size()); ++i)
0404         if (abs(in[i])<20) nInFermions++;
0405       for (int i=0; i < int(out.size()); ++i)
0406         if (abs(out[i])<20) nOutFermions++;
0407       return (nInFermions%2==0 && nOutFermions%2==0);
0408     }
0409     return false;
0410   }
0411 
0412   // Return event stripped from decay products.
0413   Event bareEvent( const Event& inputEventIn, bool storeInputEvent );
0414   // Write event with decay products attached to argument.
0415   bool reattachResonanceDecays( Event& process );
0416 
0417   // Check if particle at position iPos is part of the hard sub-system
0418   bool isInHard( int iPos, const Event& event);
0419 
0420   // Function to return the number of clustering steps for the current event
0421   virtual int getNumberOfClusteringSteps(const Event& event,
0422     bool resetNjetMax = false);
0423 
0424   //----------------------------------------------------------------------//
0425   // Functions to steer construction of histories
0426   //----------------------------------------------------------------------//
0427 
0428   // Function to force preferred picking of ordered histories. By default,
0429   // unordered histories will only be considered if no ordered histories
0430   // were found.
0431   void orderHistories( bool doOrderHistoriesIn) {
0432     doOrderHistoriesSave = doOrderHistoriesIn; }
0433   // Function to force cut on reconstructed states internally, as needed
0434   // for gg -> Higgs to ensure that e.g. uu~ -> Higgs is not constructed.
0435   void allowCutOnRecState( bool doCutOnRecStateIn) {
0436     doCutOnRecStateSave = doCutOnRecStateIn; }
0437 
0438   // Function to allow final state clusterings of weak bosons
0439   void doWeakClustering( bool doWeakClusteringIn ) {
0440     doWeakClusteringSave = doWeakClusteringIn; }
0441 
0442   //----------------------------------------------------------------------//
0443   // Functions used as default merging scales
0444   //----------------------------------------------------------------------//
0445 
0446   // Function to check if the input particle is a light jet, i.e. should be
0447   // checked against the merging scale defintion.
0448   bool checkAgainstCut( const Particle& particle);
0449   // Function to return the value of the merging scale function in the
0450   // current event.
0451   virtual double tmsNow( const Event& event );
0452   // Find the minimal Lund pT between coloured partons in the event
0453   double rhoms( const Event& event, bool withColour);
0454   // Function to calculate the minimal kT in the event
0455   double kTms(const Event & event);
0456   // Find the if the event passes the Delta R_{ij}, pT_{i} and Q_{ij} cuts on
0457   // the matrix element, as needed for cut-based merging scale definition
0458   double cutbasedms( const Event& event );
0459 
0460   //----------------------------------------------------------------------//
0461   // Functions to steer shower evolution (public to allow for PS plugin)
0462   //----------------------------------------------------------------------//
0463 
0464   // Flag to indicate trial shower usage.
0465   void doIgnoreEmissions( bool doIgnoreIn ) {
0466     doIgnoreEmissionsSave = doIgnoreIn;}
0467   // Function to allow not counting a trial emission.
0468   virtual bool canVetoEmission() { return !doIgnoreEmissionsSave; }
0469   // Function to check if emission should be rejected.
0470   virtual bool doVetoEmission( const Event& );
0471 
0472   //----------------------------------------------------------------------//
0473   // Functions used as clusterings / probabilities
0474   //----------------------------------------------------------------------//
0475 
0476   bool useShowerPluginSave;
0477   virtual bool useShowerPlugin() { return useShowerPluginSave; }
0478   virtual bool usesVincia() {return false;}
0479 
0480   //----------------------------------------------------------------------//
0481   // Functions to retrieve if merging weight should countin the internal
0482   // cross section and the event weight.
0483   //----------------------------------------------------------------------//
0484 
0485   bool includeWGTinXSEC() { return includeWGTinXSECSave;}
0486 
0487   //----------------------------------------------------------------------//
0488   // Functions to retrieve event veto information
0489   //----------------------------------------------------------------------//
0490 
0491   int nHardNow()      { return nHardNowSave; }
0492   double tmsHardNow() { return tmsHardNowSave; }
0493   int nJetsNow()      { return nJetNowSave; }
0494   double tmsNow()     { return tmsNowSave;}
0495 
0496   void setHardProcessPtr(HardProcess* hardProcIn) { hardProcess = hardProcIn; }
0497 
0498   //----------------------------------------------------------------------//
0499   // Functions related to renormalization scale variations
0500   //----------------------------------------------------------------------//
0501 
0502   int nMuRVar() { return muRVarFactors.size(); }
0503   void printIndividualWeights();
0504 
0505   //----------------------------------------------------------------------//
0506   // The members, switches etc.
0507   //----------------------------------------------------------------------//
0508 
0509   // Helper class doing all the core process book-keeping
0510   bool useOwnHardProcess;
0511   HardProcess* hardProcess;
0512 
0513   PartonLevel* showers;
0514   void setShowerPointer(PartonLevel* psIn ) {showers = psIn;}
0515 
0516   // AlphaS objects for alphaS reweighting use
0517   AlphaStrong AlphaS_FSRSave;
0518   AlphaStrong AlphaS_ISRSave;
0519   AlphaEM AlphaEM_FSRSave;
0520   AlphaEM AlphaEM_ISRSave;
0521 
0522   // Saved path to LHE file for more automated merging
0523   string lheInputFile;
0524 
0525   // Flags for merging procedure definition.
0526   bool   doUserMergingSave, doMGMergingSave, doKTMergingSave,
0527          doPTLundMergingSave, doCutBasedMergingSave, doDynamicMergingSave,
0528          includeMassiveSave, enforceStrongOrderingSave, orderInRapiditySave,
0529          pickByFullPSave, pickByPoPT2Save, includeRedundantSave,
0530          pickBySumPTSave, allowColourShufflingSave, resetHardQRenSave,
0531          resetHardQFacSave;
0532   int    unorderedScalePrescipSave, unorderedASscalePrescipSave,
0533          unorderedPDFscalePrescipSave, incompleteScalePrescipSave,
0534          ktTypeSave, nReclusterSave, nQuarksMergeSave, nRequestedSave;
0535 
0536   double scaleSeparationFactorSave, nonJoinedNormSave,
0537          fsrInRecNormSave, herwigAcollFSRSave, herwigAcollISRSave,
0538          pT0ISRSave, pTcutSave, pTminISRSave, pTminFSRSave;
0539   bool   doNL3TreeSave, doNL3LoopSave, doNL3SubtSave;
0540   bool   doUNLOPSTreeSave, doUNLOPSLoopSave, doUNLOPSSubtSave,
0541          doUNLOPSSubtNLOSave;
0542   bool   doUMEPSTreeSave, doUMEPSSubtSave;
0543 
0544   // Flag to only do phase space cut, rejecting events below the tms cut.
0545   bool   doEstimateXSection;
0546 
0547   // Flag for runtime aMC@NLO interface. Needed for aMC@NLO-Delta.
0548   bool   doRuntimeAMCATNLOInterfaceSave;
0549 
0550   // Flag for postponing event vetos. If false, events are not vetoed
0551   // based on the merging scale in CKKW-L merging, but have to be
0552   // vetoed by the user in the main program.
0553   bool   applyVeto;
0554 
0555   // Save input event in case decay products need to be detached.
0556   Event inputEvent;
0557   vector< pair<int,int> > resonances;
0558   bool doRemoveDecayProducts;
0559 
0560   // Starting scale for attaching MPI.
0561   double muMISave;
0562 
0563   // Precalculated K-factors.
0564   double kFactor0jSave;
0565   double kFactor1jSave;
0566   double kFactor2jSave;
0567 
0568   // Saved members.
0569   double tmsValueSave, tmsValueNow, DparameterSave, SparameterSave;
0570   int nJetMaxSave;
0571   int nJetMaxNLOSave;
0572 
0573   string processSave, processNow;
0574 
0575   // List of cut values to used to define a merging scale. Ordering:
0576   // 0: DeltaR_{jet_i,jet_j,min}
0577   // 1: p_{T,jet_i,min}
0578   // 2: Q_{jet_i,jet_j,min}
0579   vector<double> tmsListSave;
0580 
0581   // INTERNAL Hooks to allow construction of all histories,
0582   // including un-ordered ones
0583   bool doOrderHistoriesSave;
0584 
0585   // INTERNAL Hooks to disallow states in the construction of all histories,
0586   // e.g. because jets are below the merging scale, of to avoid the
0587   // construction of uu~ -> Higgs histories.
0588   bool doCutOnRecStateSave;
0589 
0590   // INTERNAL Hooks to allow clustering W bosons.
0591   bool doWeakClusteringSave, doSQCDClusteringSave;
0592 
0593   // Store / get first scale in PDF's that Pythia should have used
0594   double muFSave;
0595   double muRSave;
0596 
0597   // Store / get factorisation scale used in matrix element calculation.
0598   double muFinMESave;
0599   double muRinMESave;
0600 
0601   // Flag to indicate trial shower usage.
0602   bool doIgnoreEmissionsSave;
0603   // Flag to indicate if events should be vetoed.
0604   bool doIgnoreStepSave;
0605   // Stored weights in case veot needs to be revoked
0606   double pTsave;
0607   vector<double> weightCKKWL1Save, weightCKKWL2Save;
0608   int nMinMPISave;
0609   // Save CKKW-L weight / O(\alpha_s) weight.
0610   vector<double> weightCKKWLSave, weightFIRSTSave;
0611 
0612   // Struct to save individual weights
0613   struct IndividualWeights {
0614     vector<double> wtSave;
0615     vector<double> pdfWeightSave;
0616     vector<double> mpiWeightSave;
0617     vector<double> asWeightSave;
0618     vector<double> aemWeightSave;
0619     vector<double> bornAsVarFac;
0620   };
0621 
0622   IndividualWeights individualWeights;
0623 
0624   // Flag to indicate whether renormalization scale variations are performed
0625   bool doVariations;
0626   // Vector of variation factors applied to renormalization scale
0627   vector<double> muRVarFactors;
0628   // Number of weights, nominal + variations
0629   int nWgts;
0630 
0631   // Local copies of nJetMax inputs, if recalculation is necessary.
0632   int nJetMaxLocal;
0633   int nJetMaxNLOLocal;
0634   bool hasJetMaxLocal;
0635 
0636   // Event veto and hard process information, if veto should not applied be
0637   // directly, but is up to the user.
0638   bool includeWGTinXSECSave;
0639   int nHardNowSave, nJetNowSave;
0640   double tmsHardNowSave, tmsNowSave;
0641 
0642   //----------------------------------------------------------------------//
0643   // Generic setup functions
0644   //----------------------------------------------------------------------//
0645 
0646   // Function storing candidates for the hard process in the current event
0647   // Needed in order not to cluster members of the core process
0648   void storeHardProcessCandidates(const Event& event){
0649     hardProcess->storeCandidates(event,getProcessString());
0650   }
0651 
0652   // Function to set the path to the LHE file, so that more automated merging
0653   // can be used.
0654   // Remove "_1.lhe" suffix from LHE file name.
0655   // This is done so that the HarsProcess class can access both the +0 and +1
0656   // LHE files to find both the merging scale and the core process string
0657   // Store.
0658   void setLHEInputFile( string lheFile) {
0659     lheInputFile = lheFile.substr(0,lheFile.size()-6); }
0660 
0661   //----------------------------------------------------------------------//
0662   // Functions for output of class members.
0663   //----------------------------------------------------------------------//
0664 
0665   // Return AlphaStrong objects
0666   AlphaStrong* AlphaS_FSR() { return &AlphaS_FSRSave;}
0667   AlphaStrong* AlphaS_ISR() { return &AlphaS_ISRSave;}
0668   AlphaEM* AlphaEM_FSR() { return &AlphaEM_FSRSave;}
0669   AlphaEM* AlphaEM_ISR() { return &AlphaEM_ISRSave;}
0670 
0671   // Functions to return advanced merging switches
0672   // Include masses in definition of evolution pT and splitting kernels
0673   bool includeMassive() { return includeMassiveSave;}
0674   // Prefer strongly ordered histories
0675   bool enforceStrongOrdering() { return enforceStrongOrderingSave;}
0676   // Prefer histories ordered in rapidity and evolution pT
0677   bool orderInRapidity() { return orderInRapiditySave;}
0678   // Pick history probabilistically by full (correct) splitting probabilities
0679   bool pickByFull() { return pickByFullPSave;}
0680   // Pick history probabilistically, with easier form of probabilities
0681   bool pickByPoPT2() { return pickByPoPT2Save;}
0682   // Include redundant terms (e.g. PDF ratios) in the splitting probabilities
0683   bool includeRedundant(){ return includeRedundantSave;}
0684   // Pick by winner-takes-it-all, with minimum sum of scalar evolution pT
0685   bool pickBySumPT(){ return pickBySumPTSave;}
0686 
0687   // Prescription for combined scale of unordered emissions
0688   // 0 : use larger scale
0689   // 1 : use smaller scale
0690   int unorderedScalePrescip() { return unorderedScalePrescipSave;}
0691   // Prescription for combined scale used in alpha_s for unordered emissions
0692   // 0 : use combined emission scale in alpha_s weight for both (!) splittings
0693   // 1 : use original reclustered scales of each emission in alpha_s weight
0694   int unorderedASscalePrescip() { return unorderedASscalePrescipSave;}
0695   // Prescription for combined scale used in PDF ratios for unordered
0696   // emissions
0697   // 0 : use combined emission scale in PDFs for both (!) splittings
0698   // 1 : use original reclustered scales of each emission in PDF ratiost
0699   int unorderedPDFscalePrescip() { return unorderedPDFscalePrescipSave;}
0700   // Prescription for starting scale of incomplete histories
0701   // 0: use factorization scale
0702   // 1: use sHat
0703   // 2: use s
0704   int incompleteScalePrescip() { return incompleteScalePrescipSave;}
0705 
0706   // Allow swapping one colour index while reclustering
0707   bool allowColourShuffling() { return allowColourShufflingSave;}
0708 
0709   // Allow use of dynamical renormalisation scale of the core 2-> 2 process.
0710   bool resetHardQRen() { return resetHardQRenSave; }
0711   // Allow use of dynamical factorisation scale of the core 2-> 2 process.
0712   bool resetHardQFac() { return resetHardQFacSave; }
0713 
0714   // Factor by which two scales should differ to be classified strongly
0715   // ordered.
0716   double scaleSeparationFactor() { return scaleSeparationFactorSave;}
0717   // Absolute normalization of splitting probability for non-joined
0718   // evolution.
0719   double nonJoinedNorm() { return nonJoinedNormSave;}
0720   // Absolute normalization of splitting probability for final state
0721   // splittings with initial state recoiler
0722   double fsrInRecNorm() { return fsrInRecNormSave;}
0723   // Factor multiplying scalar evolution pT for FSR splitting, when picking
0724   // history by minimum scalar pT (see Jonathan Tully's thesis)
0725   double herwigAcollFSR() { return herwigAcollFSRSave;}
0726   // Factor multiplying scalar evolution pT for ISR splitting, when picking
0727   // history by minimum scalar pT (see Jonathan Tully's thesis)
0728   double herwigAcollISR() { return herwigAcollISRSave;}
0729   // ISR regularisation scale
0730   double pT0ISR() { return pT0ISRSave;}
0731   // Shower cut-off scale
0732   double pTcut() { return pTcutSave;}
0733 
0734   // MI starting scale.
0735   void muMI( double mu) { muMISave = mu; }
0736   double muMI() { return muMISave;}
0737 
0738   // Full k-Factor for current event
0739   double kFactor(int njet = 0) {
0740     return (njet == 0) ? kFactor0jSave
0741           :(njet == 1) ? kFactor1jSave
0742           : kFactor2jSave;
0743   }
0744   // O(\alhpa_s)-term of the k-Factor for current event
0745   double k1Factor( int njet = 0) {
0746     return (kFactor(njet) - 1)/infoPtr->alphaS();
0747   }
0748 
0749   // Function to return if construction of histories is biased towards ordered
0750   // histories.
0751   bool orderHistories() { return doOrderHistoriesSave;}
0752 
0753   // INTERNAL Hooks to disallow states in the construction of all histories,
0754   // e.g. because jets are below the merging scale, of to avoid the
0755   // construction of uu~ -> Higgs histories.
0756   bool allowCutOnRecState() { return doCutOnRecStateSave;}
0757 
0758   // INTERNAL Hooks to allow clustering W bosons.
0759   bool doWeakClustering() { return doWeakClusteringSave;}
0760   // INTERNAL Hooks to allow clustering clustering of gluons to squarks.
0761   bool doSQCDClustering() { return doSQCDClusteringSave;}
0762 
0763   // Store / get first scale in PDF's that Pythia should have used
0764   double muF() { return (muFSave > 0.) ? muFSave : infoPtr->QFac();}
0765   double muR() { return (muRSave > 0.) ? muRSave : infoPtr->QRen();}
0766   // Store / get factorisation scale used in matrix element calculation.
0767   double muFinME() {
0768     // Start with checking the event attribute called "muf".
0769     string mus = infoPtr->getEventAttribute("muf2",true);
0770     double mu  = (mus.empty()) ? 0. : atof((char*)mus.c_str());
0771     mu = sqrt(mu);
0772     // Check the scales tag of the event.
0773     if (infoPtr->scales) mu = infoPtr->getScalesAttribute("muf");
0774     return (mu > 0.) ? mu : (muFinMESave > 0.) ? muFinMESave
0775       : infoPtr->QFac();
0776   }
0777   double muRinME() {
0778     // Start with checking the event attribute called "mur2".
0779     string mus = infoPtr->getEventAttribute("mur2",true);
0780     double mu  = (mus.empty()) ? 0. : atof((char*)mus.c_str());
0781     mu = sqrt(mu);
0782     // Check the scales tag of the event.
0783     if (infoPtr->scales) mu = infoPtr->getScalesAttribute("mur");
0784     return (mu > 0.) ? mu : (muRinMESave > 0.) ? muRinMESave
0785       : infoPtr->QRen();
0786   }
0787 
0788 
0789   //----------------------------------------------------------------------//
0790   // Functions to steer merging code
0791   //----------------------------------------------------------------------//
0792 
0793   // Flag to indicate if events should be vetoed.
0794   void doIgnoreStep(bool doIgnoreIn) {doIgnoreStepSave = doIgnoreIn;}
0795   // Function to allow event veto.
0796   virtual bool canVetoStep() {return !doIgnoreStepSave;}
0797 
0798   // Stored weights in case veto needs to be revoked
0799   void storeWeights(vector<double> weight) {
0800     weightCKKWL1Save = weightCKKWL2Save = weight;}
0801 
0802   // Function to check event veto.
0803   virtual bool doVetoStep(const Event& process, const Event& event,
0804     bool doResonance = false);
0805 
0806   // Set starting scales
0807   virtual bool setShowerStartingScales( bool isTrial, bool doMergeFirstEmm,
0808     double& pTscaleIn, const Event& event,
0809     double& pTmaxFSRIn, bool& limitPTmaxFSRin,
0810     double& pTmaxISRIn, bool& limitPTmaxISRin,
0811     double& pTmaxMPIIn, bool& limitPTmaxMPIin );
0812 
0813   // Set shower stopping scale. Necessary to e.g. avoid accumulation of
0814   // incorrect (low-pT) shower weights through trial showering.
0815   double stopScaleSave;
0816   void setShowerStoppingScale(double scale = 0.) {stopScaleSave = scale;}
0817   double getShowerStoppingScale() {return stopScaleSave;}
0818 
0819   void nMinMPI(int nMinMPIIn) {nMinMPISave = nMinMPIIn; }
0820   int nMinMPI() {return nMinMPISave;}
0821 
0822   //----------------------------------------------------------------------//
0823   // Functions for internal merging scale definions
0824   //----------------------------------------------------------------------//
0825 
0826   // Function to calculate the kT separation between two particles
0827   double kTdurham(const Particle& RadAfterBranch,
0828     const Particle& EmtAfterBranch, int Type, double D );
0829   // Function to compute "pythia pT separation" from Particle input
0830   double rhoPythia(const Event& event, int rad, int emt, int rec,
0831     int ShowerType);
0832 
0833   // Function to find a colour (anticolour) index in the input event,
0834   // used to find colour-connected recoilers
0835   int findColour(int col, int iExclude1, int iExclude2,
0836     const Event& event, int type, bool isHardIn);
0837   // Function to compute Delta R separation from 4-vector input
0838   double deltaRij(Vec4 jet1, Vec4 jet2);
0839 
0840   //----------------------------------------------------------------------//
0841   // Functions for weight management
0842   //----------------------------------------------------------------------//
0843 
0844   // Function to get the CKKW-L weight for the current event
0845   double getWeightNLO(int i=0) { return (weightCKKWLSave[i]
0846                                          - weightFIRSTSave[i]);}
0847   // Return CKKW-L weight.
0848   vector<double> getWeightCKKWL() { return weightCKKWLSave; }
0849   // Return O(\alpha_s) weight.
0850   vector<double> getWeightFIRST() { return weightFIRSTSave; }
0851   // Set CKKW-L weight.
0852   void setWeightCKKWL( vector<double> weightIn){
0853     weightCKKWLSave = weightIn;
0854     infoPtr->weightContainerPtr
0855       ->weightsMerging.setValueVector(weightIn); }
0856   // Set O(\alpha_s) weight.
0857   void setWeightFIRST( vector<double> weightIn){
0858     weightFIRSTSave = weightIn;
0859     infoPtr->weightContainerPtr->weightsMerging
0860       .setValueFirstVector(weightIn); }
0861   // Function to return Sudakov weight as calculated before, also include MPI
0862   // weight. Only call after regular weight functions, since it is calculated
0863   // there.
0864   vector<double> getSudakovWeight() {
0865     vector<double> ret = individualWeights.wtSave;
0866     for (int i = 0; i < nWgts; ++i) {
0867      ret[i] *= individualWeights.pdfWeightSave[i] *
0868                individualWeights.mpiWeightSave[i];
0869     }
0870    return ret;
0871   }
0872   // Function to return coupling weight.
0873   vector<double> getCouplingWeight() {
0874     vector<double> ret = individualWeights.asWeightSave;
0875     for (int i = 0; i < nWgts; ++i) {
0876       ret[i] *= individualWeights.aemWeightSave[i];
0877     }
0878     return ret;
0879   }
0880 
0881 //--------------------------------------------------------------------------
0882 
0883 
0884 
0885   //----------------------------------------------------------------------//
0886   // Functions and members to store the event veto information
0887   //----------------------------------------------------------------------//
0888 
0889   // Set CKKWL veto information.
0890   void setEventVetoInfo(int nJetNowIn, double tmsNowIn) {
0891     nJetNowSave = nJetNowIn; tmsNowSave = tmsNowIn; }
0892 
0893   // Set the hard process information.
0894   void setHardProcessInfo(int nHardNowIn, double tmsHardNowIn) {
0895     nHardNowSave = nHardNowIn; tmsHardNowSave = tmsHardNowIn; }
0896 
0897   // Statistics.
0898   int nVetoedInMainShower;
0899   void addVetoInMainShower() {++nVetoedInMainShower;}
0900   int getNumberVetoedInMainShower() {return nVetoedInMainShower;}
0901 
0902 };
0903 
0904 //==========================================================================
0905 
0906 } // end namespace Pythia8
0907 
0908 #endif // Pythia8_MergingHooks_H