Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-01-18 10:06:22 

0001 // DireTimes.h is a part of the PYTHIA event generator.
0002 // Copyright (C) 2024 Stefan Prestel, 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 // Header file for the timelike final-state showers.
0007 // DireTimesEnd: data on a radiating dipole end.
0008 // DireTimes: handles the showering description.
0009 
0010 #ifndef Pythia8_DireTimes_H
0011 #define Pythia8_DireTimes_H
0012 
0013 #define DIRE_TIMES_VERSION "2.002"
0014 
0015 #include "Pythia8/Basics.h"
0016 #include "Pythia8/TimeShower.h"
0017 #include "Pythia8/BeamParticle.h"
0018 #include "Pythia8/ProcessLevel.h"
0019 #include "Pythia8/Event.h"
0020 #include "Pythia8/Info.h"
0021 #include "Pythia8/ParticleData.h"
0022 #include "Pythia8/PartonSystems.h"
0023 #include "Pythia8/PythiaStdlib.h"
0024 #include "Pythia8/Settings.h"
0025 #include "Pythia8/StandardModel.h"
0026 #include "Pythia8/UserHooks.h"
0027 #include "Pythia8/MergingHooks.h"
0028 
0029 #include "Pythia8/DireBasics.h"
0030 #include "Pythia8/DireSplittingLibrary.h"
0031 #include "Pythia8/DireWeightContainer.h"
0032 
0033 namespace Pythia8 {
0034 
0035 //==========================================================================
0036 
0037 // Data on radiating dipole ends; only used inside DireTimes class.
0038 
0039 class DireTimesEnd {
0040 
0041 public:
0042 
0043   // Constructors.
0044   DireTimesEnd() : iRadiator(-1), iRecoiler(-1), pTmax(0.), colType(0),
0045     chgType(0), gamType(0), weakType(0), isrType(0), system(0), systemRec(0),
0046     MEtype(0), iMEpartner(-1), weakPol(0), isOctetOnium(false),
0047     isHiddenValley(false), colvType(0), MEmix(0.), MEorder(true),
0048     MEsplit(true), MEgluinoRec(false), isFlexible(false), flavour(0), iAunt(0),
0049     idRadAft(0), idEmtAft(0) {
0050     mRad = m2Rad = mRec = m2Rec = mDip = m2Dip = m2DipCorr = pT2 = m2 = z
0051          = mFlavour = asymPol = flexFactor = phi = 0.;
0052     sa1  = xa = phia1 = pT2start = pT2stop = pT2Old = 0.;
0053   }
0054 
0055   DireTimesEnd(int iRadiatorIn, int iRecoilerIn, double pTmaxIn = 0.,
0056     int colIn = 0, int chgIn = 0, int gamIn = 0, int weakTypeIn = 0,
0057     int isrIn = 0, int systemIn = 0, int MEtypeIn = 0, int iMEpartnerIn = -1,
0058     int weakPolIn = 0, bool isOctetOniumIn = false,
0059     DireSingleColChain iSiblingsIn = DireSingleColChain(),
0060     bool isHiddenValleyIn = false,
0061     int colvTypeIn = 0, double MEmixIn = 0.,
0062     bool MEorderIn = true, bool MEsplitIn = true, bool MEgluinoRecIn = false,
0063     bool isFlexibleIn = false, int idRadAftIn = 0, int idEmtAftIn = 0,
0064     vector<int> iSpectatorIn = vector<int>(),
0065     vector<double> massIn = vector<double>(),
0066     vector<int> allowedIn = vector<int>() ) :
0067     iRadiator(iRadiatorIn), iRecoiler(iRecoilerIn), pTmax(pTmaxIn),
0068     colType(colIn), chgType(chgIn), gamType(gamIn), weakType(weakTypeIn),
0069     isrType(isrIn), system(systemIn), systemRec(systemIn), MEtype(MEtypeIn),
0070     iMEpartner(iMEpartnerIn), weakPol(weakPolIn), isOctetOnium(isOctetOniumIn),
0071     isHiddenValley(isHiddenValleyIn), colvType(colvTypeIn), MEmix(MEmixIn),
0072     MEorder (MEorderIn), MEsplit(MEsplitIn), MEgluinoRec(MEgluinoRecIn),
0073     isFlexible(isFlexibleIn), flavour(0), iAunt(0), mass(massIn),
0074     idRadAft(idRadAftIn), idEmtAft(idEmtAftIn), iSpectator(iSpectatorIn),
0075     allowedEmissions(allowedIn), iSiblings(iSiblingsIn) {
0076     mRad = m2Rad = mRec = m2Rec = mDip = m2Dip = m2DipCorr = pT2 = m2 = z
0077          = mFlavour = asymPol = flexFactor = phi = 0.;
0078     sa1  = xa = phia1 = 0.;
0079     pT2start = pT2stop = pT2Old = 0.;
0080   }
0081 
0082   DireTimesEnd( const DireTimesEnd& dip )
0083     : iRadiator(dip.iRadiator), iRecoiler(dip.iRecoiler), pTmax(dip.pTmax),
0084     colType(dip.colType), chgType(dip.chgType), gamType(dip.gamType),
0085     weakType(dip.weakType), isrType(dip.isrType), system(dip.system),
0086     systemRec(dip.systemRec), MEtype(dip.MEtype), iMEpartner(dip.iMEpartner),
0087     weakPol(dip.weakPol), isOctetOnium(dip.isOctetOnium),
0088     isHiddenValley(dip.isHiddenValley), colvType(dip.colvType),
0089     MEmix(dip.MEmix), MEorder(dip.MEorder), MEsplit(dip.MEsplit),
0090     MEgluinoRec(dip.MEgluinoRec), isFlexible(dip.isFlexible),
0091     flavour(dip.flavour), iAunt(dip.iAunt),
0092     mRad(dip.mRad), m2Rad(dip.m2Rad), mRec(dip.mRec), m2Rec(dip.m2Rec),
0093     mDip(dip.mDip), m2Dip(dip.m2Dip), m2DipCorr(dip.m2DipCorr), pT2(dip.pT2),
0094     m2(dip.m2), z(dip.z), mFlavour(dip.mFlavour), asymPol(dip.asymPol),
0095     flexFactor(dip.flexFactor), phi(dip.phi), pT2start(dip.pT2start),
0096     pT2stop(dip.pT2stop), pT2Old(dip.pT2Old), sa1(dip.sa1), xa(dip.xa),
0097     phia1(dip.phia1), mass(dip.mass), idRadAft(dip.idRadAft),
0098     idEmtAft(dip.idEmtAft), iSpectator(dip.iSpectator),
0099     allowedEmissions(dip.allowedEmissions), iSiblings(dip.iSiblings) {}
0100 
0101   DireTimesEnd & operator=(const DireTimesEnd& t) { if (this != &t)
0102     { iRadiator = t.iRadiator; iRecoiler = t.iRecoiler; pTmax = t.pTmax;
0103       colType = t.colType; chgType = t.chgType; gamType = t.gamType;
0104       weakType = t.weakType; isrType = t.isrType; system = t.system;
0105       systemRec = t.systemRec; MEtype = t.MEtype; iMEpartner = t.iMEpartner;
0106       weakPol = t.weakPol; isOctetOnium = t.isOctetOnium;
0107       isHiddenValley = t.isHiddenValley; colvType = t.colvType;
0108       MEmix = t.MEmix; MEorder = t.MEorder; MEsplit = t.MEsplit;
0109       MEgluinoRec = t.MEgluinoRec; isFlexible = t.isFlexible;
0110       flavour = t.flavour; iAunt = t.iAunt;
0111       mRad = t.mRad; m2Rad = t.m2Rad; mRec = t.mRec; m2Rec = t.m2Rec;
0112       mDip = t.mDip; m2Dip = t.m2Dip; m2DipCorr = t.m2DipCorr; pT2 = t.pT2;
0113       m2 = t.m2; z = t.z; mFlavour = t.mFlavour; asymPol = t.asymPol;
0114       flexFactor = t.flexFactor; phi = t.phi; pT2start = t.pT2start;
0115       pT2stop = t.pT2stop; pT2Old = t.pT2Old; sa1 = t.sa1; xa = t.xa;
0116       phia1 = t.phia1; mass = t.mass; idRadAft = t.idRadAft;
0117       idEmtAft = t.idEmtAft; iSpectator = t.iSpectator;
0118       allowedEmissions = t.allowedEmissions; iSiblings = t.iSiblings;}
0119     return *this; }
0120 
0121   // Basic properties related to dipole and matrix element corrections.
0122   int    iRadiator, iRecoiler;
0123   double pTmax;
0124   int    colType, chgType, gamType, weakType, isrType, system, systemRec,
0125          MEtype, iMEpartner, weakPol;
0126   bool   isOctetOnium, isHiddenValley;
0127   int    colvType;
0128   double MEmix;
0129   bool   MEorder, MEsplit, MEgluinoRec, isFlexible;
0130 
0131   // Properties specific to current trial emission.
0132   int    flavour, iAunt;
0133   double mRad, m2Rad, mRec, m2Rec, mDip, m2Dip, m2DipCorr,
0134          pT2, m2, z, mFlavour, asymPol, flexFactor, phi, pT2start, pT2stop,
0135          pT2Old;
0136 
0137   // Properties of 1->3 splitting.
0138   double sa1, xa,  phia1;
0139 
0140   // Stored masses.
0141   vector<double> mass;
0142 
0143   int idRadAft, idEmtAft;
0144 
0145   // Extended list of recoilers.
0146   vector<int> iSpectator;
0147   // List of allowed emissions (to avoid double-counting, since one
0148   // particle can be part of many different dipoles.
0149   void appendAllowedEmt( int id) {
0150     if ( find(allowedEmissions.begin(), allowedEmissions.end(), id)
0151       == allowedEmissions.end() ) allowedEmissions.push_back(id);
0152   }
0153   void removeAllowedEmt( int id) {
0154     if ( find(allowedEmissions.begin(), allowedEmissions.end(), id)
0155       != allowedEmissions.end() ) allowedEmissions.erase (
0156       remove(allowedEmissions.begin(), allowedEmissions.end(), id),
0157       allowedEmissions.end());
0158   }
0159   void clearAllowedEmt() { allowedEmissions.resize(0); }
0160   vector<int> allowedEmissions;
0161   bool canEmit() { return int(allowedEmissions.size() > 0); }
0162 
0163   void init(const Event& state) {
0164     mRad   = state[iRadiator].m();
0165     mRec   = state[iRecoiler].m();
0166     mDip   = sqrt( abs(2. * state[iRadiator].p() * state[iRecoiler].p()));
0167     m2Rad  = pow2(mRad);
0168     m2Rec  = pow2(mRec);
0169     m2Dip  = pow2(mDip);
0170   }
0171 
0172   void list() const {
0173     // Header.
0174     cout << "\n --------  Begin DireTimesEnd Listing  ----------------"
0175          << "------------------------------------------------------- \n \n  "
0176          << "  rad    rec       pTmax  col  parent1  parent2   isr"
0177          << "  sys sysR type  MErec    pol    m2      allowedIds\n"
0178          << fixed << setprecision(3);
0179     cout << scientific << setprecision(4)
0180       << setw(7) << iRadiator
0181       << setw(7) << iRecoiler
0182       << setw(12)<< pTmax
0183       << setw(5) << colType
0184       << setw(5) << isrType
0185       << setw(5) << system      << setw(5) << systemRec
0186       << setw(5) << MEtype      << setw(7) << iMEpartner
0187       << setw(5) << weakPol
0188       << setw(12) << m2Dip;
0189     for (int j = 0; j < int(allowedEmissions.size()); ++j)
0190       cout << setw(5) << allowedEmissions[j] << " ";
0191     cout << endl;
0192    // Done.
0193     cout << "\n --------  End DireTimesEnd Listing  ------------"
0194          << "-------------------------------------------------------" << endl;
0195   }
0196 
0197   DireSingleColChain iSiblings;
0198   void setSiblings(DireSingleColChain s) { clearSiblings(); iSiblings = s; }
0199   void clearSiblings() { iSiblings.clear(); }
0200 
0201   friend bool operator==(const DireTimesEnd& dip1, const DireTimesEnd& dip2);
0202 
0203 };
0204 
0205 //==========================================================================
0206 
0207 // The DireTimes class does timelike showers.
0208 
0209 class DireTimes : public TimeShower {
0210 
0211 public:
0212 
0213   // Constructor.
0214   DireTimes() {}
0215 
0216   DireTimes( MergingHooksPtr mergingHooksPtrIn, PartonVertexPtr ) {
0217     mergingHooksPtr   = mergingHooksPtrIn;
0218     beamOffset        = 0;
0219     userHooksPtr      = nullptr;
0220     splittingsPtr     = nullptr;
0221     weights           = 0;
0222     direInfoPtr       = nullptr;
0223     printBanner       = true;
0224     isInitSave        = false;
0225     usePDFalphas      = false;
0226     usePDF            = true;
0227     useSystems        = true;
0228     suppressLargeMECs = false;
0229   }
0230 
0231   // Destructor.
0232   virtual ~DireTimes() {}
0233 
0234   // Initialize alphaStrong and related pTmin parameters.
0235   virtual void init( BeamParticle* beamAPtrIn = nullptr,
0236     BeamParticle* beamBPtrIn = nullptr);
0237 
0238   bool initSplits() {
0239     if (splittingsPtr) splits = splittingsPtr->getSplittings();
0240     return (splits.size() > 0);
0241   }
0242 
0243   // Initialize various pointers.
0244   // (Separated from rest of init since not virtual.)
0245   void reinitPtr(Info* infoPtrIn, MergingHooksPtr mergingHooksPtrIn,
0246     DireSplittingLibrary* splittingsPtrIn, DireInfo* direInfoPtrIn) {
0247        infoPtr          = infoPtrIn;
0248        settingsPtr      = infoPtr->settingsPtr;
0249        particleDataPtr  = infoPtr->particleDataPtr;
0250        rndmPtr          = infoPtr->rndmPtr;
0251        partonSystemsPtr = infoPtr->partonSystemsPtr;
0252        userHooksPtr     = infoPtr->userHooksPtr;
0253        mergingHooksPtr  = mergingHooksPtrIn;
0254        splittingsPtr    = splittingsPtrIn;
0255        direInfoPtr      = direInfoPtrIn;
0256   }
0257 
0258   void initVariations();
0259 
0260   // Reset parton shower.
0261   void clear();
0262 
0263   void setWeightContainerPtr(DireWeightContainer* weightsIn) {
0264     weights = weightsIn;}
0265 
0266   // Find whether to limit maximum scale of emissions, and whether to dampen.
0267   virtual bool limitPTmax( Event& event, double Q2Fac = 0.,
0268     double Q2Ren = 0.);
0269 
0270   // Potential enhancement factor of pTmax scale for hardest emission.
0271   virtual double enhancePTmax() { return pTmaxFudge;}
0272 
0273   // Top-level routine to do a full time-like shower in resonance decay.
0274   virtual int shower( int iBeg, int iEnd, Event& event, double pTmax,
0275     int nBranchMax = 0);
0276 
0277   // Top-level routine for QED radiation in hadronic decay to two leptons.
0278   virtual int showerQED( int i1, int i2, Event& event, double pTmax);
0279 
0280   // Global recoil: reset counters and store locations of outgoing partons.
0281   virtual void prepareGlobal( Event&);
0282 
0283   // Prepare system for evolution after each new interaction; identify ME.
0284   virtual void prepare( int iSys, Event& event, bool limitPTmaxIn = true);
0285 
0286   // Finalize event after evolution.
0287   void finalize( Event& event);
0288 
0289   // Update dipole list after a multiparton interactions rescattering.
0290   virtual void rescatterUpdate( int iSys, Event& event);
0291 
0292   // Update dipole list after each ISR emission.
0293   virtual void update( int iSys, Event& event, bool = false);
0294 
0295   // Update dipole list after final-final splitting.
0296   void updateAfterFF( int iSysSelNow, int iSysSelRec,
0297     Event& event, int iRadBef, int iRecBef, int iRad, int iEmt, int iRec,
0298     int flavour, int colType, double pTsel);
0299 
0300   // Update dipole list after final-final splitting.
0301   void updateAfterFI( int iSysSelNow, int iSysSelRec,
0302     Event& event, int iRadBef, int iRecBef, int iRad, int iEmt, int iRec,
0303     int flavour, int colType, double pTsel, double xNew);
0304 
0305   // Select next pT in downwards evolution. Wrapper function inherited from
0306   // Pythia.
0307   virtual double pTnext( Event& event, double pTbegAll, double pTendAll,
0308     bool = false, bool = false);
0309   double newPoint( const Event& event);
0310   // Setup branching kinematics.
0311   virtual bool branch( Event& event, bool = false);
0312   bool branch_FF( Event& event, bool = false,
0313     DireSplitInfo* split = nullptr);
0314   bool branch_FI( Event& event, bool = false,
0315     DireSplitInfo* split = nullptr);
0316 
0317   pair < Vec4, Vec4 > decayWithOnshellRec( double zCS, double yCS, double phi,
0318     double m2Rec, double m2RadAft, double m2EmtAft,
0319     Vec4 pRadBef, Vec4 pRecBef);
0320   pair < Vec4, Vec4 > decayWithOffshellRec( double zCS, double yCS, double phi,
0321     double m2RadBef, double m2RadAft, double m2EmtAft,
0322     Vec4 pRadBef, Vec4 pRecBef);
0323 
0324   bool getHasWeaklyRadiated() {return false;}
0325   // Tell which system was the last processed one.
0326   int system() const {return iSysSel;};
0327 
0328   // Setup clustering kinematics.
0329   virtual Event clustered( const Event& state, int iRad, int iEmt, int iRec,
0330     string name) {
0331     pair <Event, pair<int,int> > reclus
0332       = clustered_internal(state, iRad, iEmt, iRec, name);
0333     if (reclus.first.size() > 0)
0334       reclus.first[0].mothers(reclus.second.first,reclus.second.second);
0335     return reclus.first;
0336   }
0337   pair <Event, pair<int,int> > clustered_internal( const Event& state,
0338     int iRad, int iEmt, int iRec, string name);
0339   bool cluster_FF( const Event& state, int iRad,
0340     int iEmt, int iRec, int idRadBef, Particle& radBef, Particle& recBef);
0341   bool cluster_FI( const Event& state, int iRad,
0342     int iEmt, int iRec, int idRadBef, Particle& radBef, Particle& recBef);
0343 
0344   // From Pythia version 8.215 onwards no longer virtual.
0345   double pT2Times ( const Particle& rad, const Particle& emt,
0346     const Particle& rec) {
0347     if (rec.isFinal()) return pT2_FF(rad,emt,rec);
0348     return pT2_FI(rad,emt,rec);
0349   }
0350 
0351   double pT2_FF ( const Particle& rad, const Particle& emt,
0352     const Particle& rec);
0353   double pT2_FI ( const Particle& rad, const Particle& emt,
0354     const Particle& rec);
0355 
0356   // From Pythia version 8.215 onwards no longer virtual.
0357   double zTimes ( const Particle& rad, const Particle& emt,
0358     const Particle& rec) {
0359     if (rec.isFinal()) return z_FF(rad,emt,rec);
0360     return z_FI(rad,emt,rec);
0361   }
0362 
0363   double z_FF ( const Particle& rad, const Particle& emt,
0364     const Particle& rec);
0365   double z_FI ( const Particle& rad, const Particle& emt,
0366     const Particle& rec);
0367   double z_FF_fromVec ( const Vec4& rad, const Vec4& emt, const Vec4& rec);
0368 
0369   double m2dipTimes ( const Particle& rad, const Particle& emt,
0370     const Particle& rec) {
0371     if (rec.isFinal()) return m2dip_FF(rad,emt,rec);
0372     return m2dip_FI(rad,emt,rec);
0373   }
0374 
0375   double m2dip_FF ( const Particle& rad, const Particle& emt,
0376     const Particle& rec);
0377   double m2dip_FI ( const Particle& rad, const Particle& emt,
0378     const Particle& rec);
0379 
0380   // From Pythia version 8.218 onwards.
0381   // Return the evolution variable.
0382   // Usage: getStateVariables( const Event& event,  int iRad, int iEmt,
0383   //                   int iRec, string name)
0384   // Important note:
0385   // - This map must contain an entry for the shower evolution variable,
0386   //   specified with key "t".
0387   // - This map must contain an entry for the shower evolution variable from
0388   //   which the shower would be restarted after a branching. This entry
0389   //   must have key "tRS",
0390   // - This map must contain an entry for the argument of \alpha_s used
0391   //   for the branching. This entry must have key "scaleAS".
0392   // - This map must contain an entry for the argument of the PDFs used
0393   //   for the branching. This entry must have key "scalePDF".
0394   virtual map<string, double> getStateVariables (const Event& state,
0395     int rad, int emt, int rec, string name);
0396 
0397   // From Pythia version 8.215 onwards.
0398   // Check if attempted clustering is handled by timelike shower
0399   // Usage: isTimelike( const Event& event,  int iRad, int iEmt,
0400   //                   int iRec, string name)
0401   virtual bool isTimelike(const Event& state, int iRad, int, int, string)
0402     { return state[iRad].isFinal(); }
0403 
0404   // From Pythia version 8.215 onwards.
0405   // Return a string identifier of a splitting.
0406   // Usage: getSplittingName( const Event& event, int iRad, int iEmt, int iRec)
0407   virtual vector<string> getSplittingName( const Event& state, int iRad,
0408     int iEmt, int)
0409     { return splittingsPtr->getSplittingName(state,iRad,iEmt); }
0410 
0411   // From Pythia version 8.215 onwards.
0412   // Return the splitting probability.
0413   // Usage: getSplittingProb( const Event& event, int iRad, int iEmt, int iRec)
0414   virtual double getSplittingProb( const Event& state, int iRad,
0415     int iEmt, int iRec, string name);
0416 
0417   virtual bool allowedSplitting( const Event& state, int iRad, int iEmt);
0418 
0419   virtual vector<int> getRecoilers( const Event& state, int iRad, int iEmt,
0420     string name);
0421 
0422   virtual double getCoupling( double mu2Ren, string name) {
0423     if (splits.find(name) != splits.end())
0424       return splits[name]->coupling(-1.,mu2Ren, 0., 1.);
0425     return 1.;
0426   }
0427 
0428   bool isSymmetric( string name, const Particle* rad, const Particle* emt) {
0429     if (splits.find(name) != splits.end())
0430       return splits[name]->isSymmetric(rad,emt);
0431     return false;
0432   }
0433 
0434   // Auxiliary function to return the position of a particle.
0435   // Should go int Event class eventually!
0436   int FindParticle( const Particle& particle, const Event& event,
0437     bool checkStatus = true );
0438 
0439   // Print dipole list; for debug mainly.
0440   virtual void list() const;
0441 
0442   Event makeHardEvent( int iSys, const Event& state, bool isProcess = false );
0443 
0444   // Check that particle has sensible momentum.
0445   bool validMomentum( const Vec4& p, int id, int status);
0446 
0447   // Check colour/flavour correctness of state.
0448   bool validEvent( const Event& state, bool isProcess = false,
0449     int iSysCheck = -1 );
0450 
0451   // Check that mother-daughter-relations are correctly set.
0452   bool validMotherDaughter( const Event& state );
0453 
0454   // Find index of colour partner for input colour.
0455   int FindCol(int col, vector<int> iExclude, const Event& event, int type,
0456     int iSys = -1);
0457 
0458   // Pointers to the two incoming beams.
0459   BeamParticle*  getBeamA () { return beamAPtr; }
0460   BeamParticle*  getBeamB () { return beamBPtr; }
0461 
0462   // Function to calculate the correct alphaS/2*Pi value, including
0463   // renormalisation scale variations + threshold matching.
0464   double alphasNow( double pT2, double renormMultFacNow = 1., int iSys = 0 );
0465 
0466   // Function to calculate the correct alphaEM/2*Pi value.
0467   double alphaemNow( double pT2, double renormMultFacNow = 1., int iSys = 0 );
0468 
0469   bool isInit() { return isInitSave; }
0470 
0471   // Function to calculate the absolute phase-sace boundary for emissions.
0472   double m2Max (int iDip, const Event& state) {
0473     if ( state[dipEnd[iDip].iRecoiler].isFinal()
0474       && state[dipEnd[iDip].iRadiator].isFinal())
0475       return dipEnd[iDip].m2Dip;
0476     int iSys = dipEnd[iDip].system;
0477     int inA = partonSystemsPtr->getInA(iSys);
0478     int inB = partonSystemsPtr->getInB(iSys);
0479     double x = 1.;
0480     int iRad(dipEnd[iDip].iRadiator), iRec(dipEnd[iDip].iRecoiler);
0481     if (hasPDF(state[iRad].id()) && iRad == inA)
0482       x *= state[inA].pPos() / state[0].m();
0483     if (hasPDF(state[iRad].id()) && iRad == inB)
0484       x *= state[inB].pNeg() / state[0].m();
0485     if (hasPDF(state[iRec].id()) && iRec == inA)
0486       x *= state[inA].pPos() / state[0].m();
0487     if (hasPDF(state[iRec].id()) && iRec == inB)
0488       x *= state[inB].pNeg() / state[0].m();
0489     return dipEnd[iDip].m2Dip/x;
0490   }
0491 
0492   bool dryrun;
0493 
0494 private:
0495 
0496   friend class DireSplitting;
0497   friend class DireSpace;
0498 
0499   // Number of times the same error message is repeated, unless overridden.
0500   static const int TIMESTOPRINT;
0501 
0502   // Allow conversion from mb to pb.
0503   static const double CONVERTMB2PB;
0504 
0505   // Colour factors.
0506   //static const double CA, CF, TR, NC, LEPTONZMAX;
0507   static const double LEPTONZMAX;
0508   double CA, CF, TR, NC;
0509 
0510   // Store common beam quantities.
0511   int    idASave, idBSave;
0512 
0513 protected:
0514 
0515   // Store properties to be returned by methods.
0516   int    iSysSel;
0517   double pTmaxFudge, pTLastBranch;
0518 
0519 private:
0520 
0521   // Constants: could only be changed in the code itself.
0522   static const int MAXLOOPTINYPDF;
0523   static const double MCMIN, MBMIN, SIMPLIFYROOT, XMARGIN, XMARGINCOMB,
0524          TINYPDF, LARGEM2, THRESHM2, LAMBDA3MARGIN, TINYMASS, TINYOVERESTIMATE,
0525          PT2MIN_PDF_IN_OVERESTIMATE, PT2_INCREASE_OVERESTIMATE,
0526          KERNEL_HEADROOM;
0527 
0528   // Initialization data, normally only set once.
0529   bool   isInitSave, doQCDshower, doQEDshowerByQ, doQEDshowerByL,
0530          doMEcorrections, doMEafterFirst, doPhiPolAsym,
0531          doInterleave, allowBeamRecoil, dampenBeamRecoil, recoilToColoured,
0532          useFixedFacScale, allowRescatter, canVetoEmission, hasUserHooks,
0533          doSecondHard, alphaSuseCMW, printBanner, doTrialNow;
0534   int    pTmaxMatch, pTdampMatch, alphaSorder, alphaSnfmax, alphaEMorder,
0535          nGluonToQuark, nGammaToQuark, nGammaToLepton, nFinalMax,
0536          nFinalMaxMECs,kernelOrder, kernelOrderMPI, nMPI, asScheme;
0537   double pTdampFudge, mc, mb, m2c, m2b, renormMultFac, factorMultFac,
0538          fixedFacScale2, alphaSvalue, alphaS2pi, Lambda3flav, Lambda4flav,
0539          Lambda5flav, Lambda3flav2, Lambda4flav2, Lambda5flav2,
0540          pTcolCutMin, pTcolCut,
0541          pT2colCut, m2colCut, mTolErr, mZ, gammaZ, thetaW, mW, gammaW,
0542          pTmaxFudgeMPI, sumCharge2L, sumCharge2Q, sumCharge2Tot,
0543          pT2minVariations, pT2minEnhance, pT2minMECs, Q2minMECs, pT2recombine,
0544          m2cPhys, m2bPhys;
0545   double alphaS2piOverestimate;
0546   bool usePDFalphas, usePDFmasses, useSummedPDF, usePDF, useSystems,
0547        useMassiveBeams, suppressLargeMECs;
0548 
0549   double pTchgQCut, pT2chgQCut, pTchgLCut, pT2chgLCut;
0550 
0551   unordered_map<int,double> pT2cutSave;
0552   double pT2cut(int id) {
0553     if (pT2cutSave.find(id) != pT2cutSave.end()) return pT2cutSave[id];
0554     // Else return maximal value.
0555     double ret = 0.;
0556     for ( unordered_map<int,double>::iterator it = pT2cutSave.begin();
0557       it != pT2cutSave.end(); ++it ) ret = max(ret, it->second);
0558     return ret;
0559   }
0560   double pT2cutMax(DireTimesEnd* dip) {
0561     double ret = 0.;
0562     for (int i=0; i < int(dip->allowedEmissions.size()); ++i)
0563       ret = max( ret, pT2cut(dip->allowedEmissions[i]));
0564     return ret;
0565   }
0566   double pT2cutMin(DireTimesEnd* dip) {
0567     double ret = 1e15;
0568     for (int i=0; i < int(dip->allowedEmissions.size()); ++i)
0569       ret = min( ret, pT2cut(dip->allowedEmissions[i]));
0570     return ret;
0571   }
0572 
0573   bool doDecaysAsShower;
0574 
0575   // alphaStrong and alphaEM calculations.
0576   AlphaStrong alphaS;
0577   AlphaEM     alphaEM;
0578 
0579   // Some current values.
0580   bool   dopTlimit1, dopTlimit2, dopTdamp;
0581   double pT2damp, kRad, kEmt, pdfScale2;
0582 
0583   // All dipole ends and a pointer to the selected hardest dipole end.
0584   vector<DireTimesEnd> dipEnd;
0585   DireTimesEnd* dipSel;
0586   DireSplitInfo splitInfoSel;
0587   DireSplitting* splittingSel;
0588   int iDipSel;
0589   unordered_map<string,double> kernelSel, kernelNow;
0590   double auxSel, overSel, boostSel, auxNow, overNow, boostNow;
0591 
0592   double tinypdf( double x) {
0593     double xref = 0.01;
0594     return TINYPDF*log(1-x)/log(1-xref);
0595   }
0596 
0597   // Function to check if id is part of the incoming hadron state.
0598   bool hasPDF (int id) {
0599     if ( !usePDF )                          return false;
0600     if ( particleDataPtr->colType(id) != 0) return true;
0601     if ( particleDataPtr->isLepton(id)
0602       && settingsPtr->flag("PDF:lepton"))   return true;
0603     return false;
0604   }
0605 
0606   // Wrapper around PDF calls.
0607   double getXPDF( int id, double x, double t, int iSys = 0,
0608     BeamParticle* beam = nullptr, bool finalRec = true, double z = 0.,
0609     double m2dip = 0.) {
0610     // Return one if no PDF should be used.
0611     if (!hasPDF(id)) return 1.0;
0612     // Else get PDF from beam particle.
0613     BeamParticle* b = beam;
0614     if (b == nullptr) {
0615       if (beamAPtr != nullptr || beamBPtr != nullptr) {
0616         b = (beamAPtr != nullptr && particleDataPtr->isHadron(beamAPtr->id()))
0617             ? beamAPtr
0618           : (beamBPtr != nullptr && particleDataPtr->isHadron(beamBPtr->id()))
0619             ? beamBPtr : nullptr;
0620       }
0621       if (b == nullptr && beamAPtr != 0) b = beamAPtr;
0622       if (b == nullptr && beamBPtr != 0) b = beamBPtr;
0623     }
0624 
0625     double scale2 = t;
0626     if (asScheme == 2 && z != 0. && finalRec) {
0627       double zcs = z;
0628       double xcs = m2dip * zcs * (1.-zcs) / (t + m2dip * zcs * (1.-zcs));
0629       scale2 = (1-zcs)*(1-xcs)/xcs/zcs*m2dip;
0630     }
0631 
0632     double ret =  (useSummedPDF) ? b->xf(id, x, scale2)
0633                                  : b->xfISR(iSys,id, x, scale2);
0634     // Done.
0635     return ret;
0636   }
0637 
0638   // Evolve a QCD dipole end near heavy quark threshold region.
0639   // Setup a dipole end, either QCD, QED/photon, weak or Hidden Valley one.
0640   void setupQCDdip( int iSys, int i, int colTag,  int colSign, Event& event,
0641     bool isOctetOnium = false, bool limitPTmaxIn = true);
0642   void getGenDip( int iSys, int i, int iRad, const Event& event,
0643     bool limitPTmaxIn, vector<DireTimesEnd>& dipEnds );
0644   void getQCDdip( int iRad, int colTag, int colSign, const Event& event,
0645     vector<DireTimesEnd>& dipEnds );
0646   void setupDecayDip( int iSys, int iRad, const Event& event,
0647     vector<DireTimesEnd>& dipEnds);
0648 
0649   // Function to set up and append a new dipole.
0650   bool appendDipole( const Event& state, int iRad, int iRec, double pTmax,
0651     int colType, int chgType, int gamType, int weakType, int isrType, int iSys,
0652     int MEtype, int iMEpartner, int weakPol, bool isOctetOnium,
0653     vector<DireTimesEnd>& dipEnds);
0654 
0655   vector<int> sharedColor(const Particle& rad, const Particle& rec);
0656 
0657   // Function to set up and append a new dipole.
0658   void updateDipoles(const Event& state, int iSys = -1);
0659   void checkDipoles(const Event& state);
0660   void saveSiblings(const Event& state, int iSys = -1);
0661   bool updateAllowedEmissions( const Event& state, DireTimesEnd* dip);
0662   bool appendAllowedEmissions( const Event& state, DireTimesEnd* dip);
0663 
0664   // Evolve a QCD dipole end.
0665   void pT2nextQCD( double pT2begDip, double pT2sel, DireTimesEnd& dip,
0666     Event& event, double pT2endForce = -1., double pT2freeze = 0.,
0667     bool forceBranching = false);
0668   bool pT2nextQCD_FF( double pT2begDip, double pT2sel, DireTimesEnd& dip,
0669     const Event& event, double pT2endForce = -1., double pT2freeze = 0.,
0670     bool forceBranching = false);
0671   bool pT2nextQCD_FI( double pT2begDip, double pT2sel, DireTimesEnd& dip,
0672     const Event& event, double pT2endForce = -1., double pT2freeze = 0.,
0673     bool forceBranching = false);
0674 
0675   double tNextQCD( DireTimesEnd*, double overestimateInt,
0676     double tOld, double tMin, double tFreeze=0., int algoType = 0);
0677   bool zCollNextQCD( DireTimesEnd* dip, double zMin, double zMax,
0678     double tMin = 0., double tMax = 0.);
0679   bool virtNextQCD( DireTimesEnd* dip, double tMin, double tMax,
0680     double zMin =-1., double zMax =-1.);
0681 
0682   // Function to determine how often the integrated overestimate should be
0683   // recalculated.
0684   double evalpdfstep(int idRad, double pT2, double m2cp = -1.,
0685     double m2bp = -1.) {
0686     double ret = 0.2;
0687     if (m2cp < 0.) m2cp = particleDataPtr->m0(4);
0688     if (m2bp < 0.) m2bp = particleDataPtr->m0(5);
0689     // More steps close to the thresholds.
0690     if ( abs(idRad) == 4 && pT2 < 1.2*m2cp && pT2 > m2cp) ret = 1.0;
0691     if ( abs(idRad) == 5 && pT2 < 1.2*m2bp && pT2 > m2bp) ret = 1.0;
0692     return ret;
0693   }
0694 
0695   DireSplittingLibrary* splittingsPtr;
0696 
0697   // Number of proposed splittings in hard scattering systems.
0698   unordered_map<int,int> nProposedPT;
0699 
0700   // Return headroom factors for integrated/differential overestimates.
0701   double overheadFactors(DireTimesEnd*, const Event&, string, double,
0702     double, double);
0703   double enhanceOverestimateFurther( string, int, double );
0704   double overheadFactorsMEC(const Event&, DireSplitInfo*, string);
0705 
0706   // Function to fill map of integrated overestimates.
0707   void getNewOverestimates( DireTimesEnd*, const Event&, double, double,
0708     double, double, multimap<double,string>&);
0709 
0710   // Function to sum all integrated overestimates.
0711   void addNewOverestimates( multimap<double,string>, double&);
0712 
0713   // Function to attach the correct alphaS weights to the kernels.
0714   void alphasReweight(double t, double talpha, int iSys, bool forceFixedAs,
0715     double& weight, double& fullWeight, double& overWeight,
0716     double renormMultFacNow);
0717 
0718   // Function to evaluate the accept-probability, including picking of z.
0719   void getNewSplitting( const Event&, DireTimesEnd*, double, double, double,
0720     double, double, int, string, bool, int&, int&, double&, double&,
0721     unordered_map<string,double>&, double&);
0722 
0723   pair<bool, pair<double,double> > getMEC ( const Event& state,
0724     DireSplitInfo* splitInfo);
0725   bool applyMEC ( const Event& state, DireSplitInfo* splitInfo,
0726     vector<Event> auxEvent = vector<Event>() );
0727 
0728   // Get particle masses.
0729   double getMass(int id, int strategy, double mass = 0.) {
0730     BeamParticle* beam = nullptr;
0731     if (beamAPtr != nullptr || beamBPtr != nullptr) {
0732       beam = (beamAPtr != nullptr && particleDataPtr->isHadron(beamAPtr->id()))
0733            ? beamAPtr
0734            : (beamBPtr != nullptr && particleDataPtr->isHadron(beamBPtr->id()))
0735               ? beamBPtr : nullptr;
0736     }
0737     bool usePDFmass = usePDFmasses
0738       && (toLower(settingsPtr->word("PDF:pSet")).find("lhapdf")
0739          != string::npos);
0740     double mRet = 0.;
0741     // Parton masses.
0742     if ( particleDataPtr->colType(id) != 0) {
0743       if (strategy == 1) mRet = particleDataPtr->m0(id);
0744       if (strategy == 2 &&  usePDFmass && beam != nullptr)
0745         mRet = beam->mQuarkPDF(id);
0746       if (strategy == 2 && (!usePDFmass || beam == nullptr))
0747         mRet = particleDataPtr->m0(id);
0748       if (strategy == 3) mRet = mass;
0749       if (mRet < TINYMASS) mRet = 0.;
0750     // Masses of other particles.
0751     } else {
0752       mRet = particleDataPtr->m0(id);
0753       if (strategy == 3) mRet = mass;
0754       if (mRet < TINYMASS) mRet = 0.;
0755     }
0756     return pow2(max(0.,mRet));
0757   }
0758 
0759   // Check if variables are in allowed phase space.
0760   bool inAllowedPhasespace(int kinType, double z, double pT2, double m2dip,
0761     double q2, double xOld, int splitType = 0, double m2RadBef = 0.,
0762     double m2r = 0., double m2s = 0., double m2e = 0.,
0763     vector<double> aux = vector<double>());
0764 
0765   // Auxiliary function to get number of flavours.
0766   double getNF(double pT2);
0767 
0768   // Auxiliary functions to get beta function coefficients.
0769   double beta0 (double NF)
0770     { return 11./6.*CA - 2./3.*NF*TR; }
0771   double beta1 (double NF)
0772     { return 17./6.*pow2(CA) - (5./3.*CA+CF)*NF*TR; }
0773   double beta2 (double NF)
0774     { return 2857./432.*pow(CA,3)
0775     + (-1415./216.*pow2(CA) - 205./72.*CA*CF + pow2(CF)/4.) *TR*NF
0776     + ( 79.*CA + 66.*CF)/108.*pow2(TR*NF); }
0777 
0778   // Identifier of the splitting
0779   string splittingNowName, splittingSelName;
0780 
0781   // Weighted shower book-keeping.
0782   unordered_map<string, map<double,double> > acceptProbability;
0783   unordered_map<string, multimap<double,double> > rejectProbability;
0784 
0785 public:
0786 
0787   DireWeightContainer* weights;
0788   DireInfo* direInfoPtr;
0789   ProcessLevel processLevel;
0790   unordered_map<string, DireSplitting* > splits;
0791   vector<int> bornColors;
0792 
0793 private:
0794 
0795   bool doVariations;
0796 
0797   // List of splitting kernels.
0798   //map<string, DireSplitting* > splits;
0799   unordered_map<string, double > overhead;
0800   void scaleOverheadFactor(string name, double scale) {
0801     overhead[name] *= scale;
0802     return;
0803   }
0804   void resetOverheadFactors() {
0805     for ( unordered_map<string,double>::iterator it = overhead.begin();
0806       it != overhead.end(); ++it )
0807       it->second = 1.0;
0808     return;
0809   }
0810 
0811   double octetOniumColFac;
0812   bool useLocalRecoilNow;
0813 
0814   // Map to store some settings, to be passes to splitting kernels.
0815   unordered_map<string,bool> bool_settings;
0816 
0817 };
0818 
0819 //==========================================================================
0820 
0821 } // end namespace Pythia8
0822 
0823 #endif // Pythia8_DireTimes_H
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!