Line data Source code
1 : // TimeShower.h is a part of the PYTHIA event generator.
2 : // Copyright (C) 2015 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 : // Header file for the timelike final-state showers.
7 : // TimeDipoleEnd: data on a radiating dipole end.
8 : // TimeShower: handles the showering description.
9 :
10 : #ifndef Pythia8_TimeShower_H
11 : #define Pythia8_TimeShower_H
12 :
13 : #include "Pythia8/Basics.h"
14 : #include "Pythia8/BeamParticle.h"
15 : #include "Pythia8/Event.h"
16 : #include "Pythia8/Info.h"
17 : #include "Pythia8/ParticleData.h"
18 : #include "Pythia8/PartonSystems.h"
19 : #include "Pythia8/PythiaStdlib.h"
20 : #include "Pythia8/Settings.h"
21 : #include "Pythia8/StandardModel.h"
22 : #include "Pythia8/UserHooks.h"
23 : #include "Pythia8/MergingHooks.h"
24 : #include "Pythia8/WeakShowerMEs.h"
25 :
26 : namespace Pythia8 {
27 :
28 : //==========================================================================
29 :
30 : // Data on radiating dipole ends; only used inside TimeShower class.
31 :
32 : class TimeDipoleEnd {
33 :
34 : public:
35 :
36 : // Constructors.
37 0 : TimeDipoleEnd() : iRadiator(-1), iRecoiler(-1), pTmax(0.), colType(0),
38 0 : chgType(0), gamType(0), weakType(0), isrType(0), system(0), systemRec(0),
39 0 : MEtype(0), iMEpartner(-1), weakPol(0), isOctetOnium(false),
40 0 : isHiddenValley(false), colvType(0), MEmix(0.), MEorder(true),
41 0 : MEsplit(true), MEgluinoRec(false), isFlexible(false) { }
42 : TimeDipoleEnd(int iRadiatorIn, int iRecoilerIn, double pTmaxIn = 0.,
43 : int colIn = 0, int chgIn = 0, int gamIn = 0, int weakTypeIn = 0,
44 : int isrIn = 0, int systemIn = 0, int MEtypeIn = 0, int iMEpartnerIn = -1,
45 : int weakPolIn = 0, bool isOctetOniumIn = false,
46 : bool isHiddenValleyIn = false, int colvTypeIn = 0, double MEmixIn = 0.,
47 : bool MEorderIn = true, bool MEsplitIn = true, bool MEgluinoRecIn = false,
48 : bool isFlexibleIn = false) :
49 0 : iRadiator(iRadiatorIn), iRecoiler(iRecoilerIn), pTmax(pTmaxIn),
50 0 : colType(colIn), chgType(chgIn), gamType(gamIn), weakType(weakTypeIn),
51 0 : isrType(isrIn), system(systemIn), systemRec(systemIn), MEtype(MEtypeIn),
52 0 : iMEpartner(iMEpartnerIn), weakPol(weakPolIn), isOctetOnium(isOctetOniumIn),
53 0 : isHiddenValley(isHiddenValleyIn), colvType(colvTypeIn), MEmix(MEmixIn),
54 0 : MEorder (MEorderIn), MEsplit(MEsplitIn), MEgluinoRec(MEgluinoRecIn),
55 0 : isFlexible(isFlexibleIn) { }
56 :
57 : // Basic properties related to dipole and matrix element corrections.
58 : int iRadiator, iRecoiler;
59 : double pTmax;
60 : int colType, chgType, gamType, weakType, isrType, system, systemRec,
61 : MEtype, iMEpartner, weakPol;
62 : bool isOctetOnium, isHiddenValley;
63 : int colvType;
64 : double MEmix;
65 : bool MEorder, MEsplit, MEgluinoRec, isFlexible;
66 :
67 : // Properties specific to current trial emission.
68 : int flavour, iAunt;
69 : double mRad, m2Rad, mRec, m2Rec, mDip, m2Dip, m2DipCorr,
70 : pT2, m2, z, mFlavour, asymPol, flexFactor;
71 :
72 : } ;
73 :
74 : //==========================================================================
75 :
76 : // The TimeShower class does timelike showers.
77 :
78 : class TimeShower {
79 :
80 : public:
81 :
82 : // Constructor.
83 0 : TimeShower() {beamOffset = 0;}
84 :
85 : // Destructor.
86 0 : virtual ~TimeShower() {}
87 :
88 : // Initialize various pointers.
89 : // (Separated from rest of init since not virtual.)
90 : void initPtr(Info* infoPtrIn, Settings* settingsPtrIn,
91 : ParticleData* particleDataPtrIn, Rndm* rndmPtrIn,
92 : CoupSM* coupSMPtrIn, PartonSystems* partonSystemsPtrIn,
93 : UserHooks* userHooksPtrIn, MergingHooks* mergingHooksPtrIn = 0) {
94 0 : infoPtr = infoPtrIn; settingsPtr = settingsPtrIn;
95 0 : particleDataPtr = particleDataPtrIn; rndmPtr = rndmPtrIn;
96 0 : coupSMPtr = coupSMPtrIn; partonSystemsPtr = partonSystemsPtrIn;
97 0 : userHooksPtr = userHooksPtrIn; mergingHooksPtr = mergingHooksPtrIn;}
98 :
99 : // Initialize alphaStrong and related pTmin parameters.
100 : virtual void init( BeamParticle* beamAPtrIn = 0,
101 : BeamParticle* beamBPtrIn = 0);
102 :
103 : // New beams possible for handling of hard diffraction. (Not virtual.)
104 : void reassignBeamPtrs( BeamParticle* beamAPtrIn, BeamParticle* beamBPtrIn,
105 0 : int beamOffsetIn = 0) {beamAPtr = beamAPtrIn; beamBPtr = beamBPtrIn;
106 0 : beamOffset = beamOffsetIn;}
107 :
108 : // Find whether to limit maximum scale of emissions, and whether to dampen.
109 : virtual bool limitPTmax( Event& event, double Q2Fac = 0.,
110 : double Q2Ren = 0.);
111 :
112 : // Potential enhancement factor of pTmax scale for hardest emission.
113 0 : virtual double enhancePTmax() {return pTmaxFudge;}
114 :
115 : // Top-level routine to do a full time-like shower in resonance decay.
116 : virtual int shower( int iBeg, int iEnd, Event& event, double pTmax,
117 : int nBranchMax = 0);
118 :
119 : // Top-level routine for QED radiation in hadronic decay to two leptons.
120 : virtual int showerQED( int i1, int i2, Event& event, double pTmax);
121 :
122 : // Provide the pT scale of the last branching in the above shower.
123 : double pTLastInShower() {return pTLastBranch;}
124 :
125 : // Global recoil: reset counters and store locations of outgoing partons.
126 : virtual void prepareGlobal( Event& event);
127 :
128 : // Prepare system for evolution after each new interaction; identify ME.
129 : virtual void prepare( int iSys, Event& event, bool limitPTmaxIn = true);
130 :
131 : // Update dipole list after a multiparton interactions rescattering.
132 : virtual void rescatterUpdate( int iSys, Event& event);
133 :
134 : // Update dipole list after each ISR emission.
135 : virtual void update( int iSys, Event& event, bool hasWeakRad = false);
136 :
137 : // Select next pT in downwards evolution.
138 : virtual double pTnext( Event& event, double pTbegAll, double pTendAll,
139 : bool isFirstTrial = false);
140 :
141 : // ME corrections and kinematics that may give failure.
142 : virtual bool branch( Event& event, bool isInterleaved = false);
143 :
144 : // Tell which system was the last processed one.
145 0 : int system() const {return iSysSel;};
146 :
147 : // Tell whether FSR has done a weak emission.
148 0 : bool getHasWeaklyRadiated() {return hasWeaklyRadiated;}
149 :
150 : // Print dipole list; for debug mainly.
151 : virtual void list( ostream& os = cout) const;
152 :
153 : protected:
154 :
155 : // Pointer to various information on the generation.
156 : Info* infoPtr;
157 :
158 : // Pointer to the settings database.
159 : Settings* settingsPtr;
160 :
161 : // Pointer to the particle data table.
162 : ParticleData* particleDataPtr;
163 :
164 : // Pointer to the random number generator.
165 : Rndm* rndmPtr;
166 :
167 : // Pointer to Standard Model couplings.
168 : CoupSM* coupSMPtr;
169 :
170 : // Pointers to the two incoming beams. Offset their location in event.
171 : BeamParticle* beamAPtr;
172 : BeamParticle* beamBPtr;
173 : int beamOffset;
174 :
175 : // Pointer to information on subcollision parton locations.
176 : PartonSystems* partonSystemsPtr;
177 :
178 : // Pointer to userHooks object for user interaction with program.
179 : UserHooks* userHooksPtr;
180 :
181 : // Weak matrix elements used for corrections both of ISR and FSR.
182 : WeakShowerMEs weakShowerMEs;
183 :
184 : // Store properties to be returned by methods.
185 : int iSysSel;
186 : double pTmaxFudge, pTLastBranch;
187 :
188 : private:
189 :
190 : // Constants: could only be changed in the code itself.
191 : static const double MCMIN, MBMIN, SIMPLIFYROOT, XMARGIN, XMARGINCOMB,
192 : TINYPDF, LARGEM2, THRESHM2, LAMBDA3MARGIN, WEAKPSWEIGHT, WG2QEXTRA;
193 : // Rescatter: try to fix up recoil between systems
194 : static const bool FIXRESCATTER, VETONEGENERGY;
195 : static const double MAXVIRTUALITYFRACTION, MAXNEGENERGYFRACTION;
196 :
197 : // Initialization data, normally only set once.
198 : bool doQCDshower, doQEDshowerByQ, doQEDshowerByL, doQEDshowerByGamma,
199 : doWeakShower, doMEcorrections, doMEafterFirst, doPhiPolAsym,
200 : doInterleave, allowBeamRecoil, dampenBeamRecoil, recoilToColoured,
201 : useFixedFacScale, allowRescatter, canVetoEmission, doHVshower,
202 : brokenHVsym, globalRecoil, useLocalRecoilNow, doSecondHard,
203 : singleWeakEmission, alphaSuseCMW, vetoWeakJets, allowMPIdipole;
204 : int pTmaxMatch, pTdampMatch, alphaSorder, alphaSnfmax, nGluonToQuark,
205 : weightGluonToQuark, alphaEMorder, nGammaToQuark, nGammaToLepton,
206 : nCHV, idHV, nMaxGlobalRecoil, weakMode;
207 : double pTdampFudge, mc, mb, m2c, m2b, renormMultFac, factorMultFac,
208 : fixedFacScale2, alphaSvalue, alphaS2pi, Lambda3flav, Lambda4flav,
209 : Lambda5flav, Lambda3flav2, Lambda4flav2, Lambda5flav2,
210 : scaleGluonToQuark, extraGluonToQuark, pTcolCutMin, pTcolCut,
211 : pT2colCut, pTchgQCut, pT2chgQCut, pTchgLCut, pT2chgLCut,
212 : pTweakCut, pT2weakCut, mMaxGamma, m2MaxGamma, octetOniumFraction,
213 : octetOniumColFac, mZ, gammaZ, thetaWRat, mW, gammaW, CFHV,
214 : alphaHVfix, pThvCut, pT2hvCut, mHV, pTmaxFudgeMPI,
215 : weakEnhancement, vetoWeakDeltaR2;
216 :
217 : // alphaStrong and alphaEM calculations.
218 : AlphaStrong alphaS;
219 : AlphaEM alphaEM;
220 :
221 : // Some current values.
222 : bool dopTlimit1, dopTlimit2, dopTdamp, hasWeaklyRadiated;
223 : double pT2damp, kRad, kEmt, pdfScale2;
224 :
225 : // All dipole ends and a pointer to the selected hardest dipole end.
226 : vector<TimeDipoleEnd> dipEnd;
227 : TimeDipoleEnd* dipSel;
228 : int iDipSel;
229 :
230 : // Setup a dipole end, either QCD, QED/photon, weak or Hidden Valley one.
231 : void setupQCDdip( int iSys, int i, int colTag, int colSign, Event& event,
232 : bool isOctetOnium = false, bool limitPTmaxIn = true);
233 : void setupQEDdip( int iSys, int i, int chgType, int gamType, Event& event,
234 : bool limitPTmaxIn = true);
235 : void setupWeakdip( int iSys, int i,int weakType, Event& event,
236 : bool limitPTmaxIn = true);
237 : void setupHVdip( int iSys, int i, Event& event, bool limitPTmaxIn = true);
238 :
239 : // Evolve a QCD dipole end.
240 : void pT2nextQCD( double pT2begDip, double pT2sel, TimeDipoleEnd& dip,
241 : Event& event);
242 :
243 : // Evolve a QED dipole end, either charged or photon.
244 : void pT2nextQED( double pT2begDip, double pT2sel, TimeDipoleEnd& dip,
245 : Event& event);
246 :
247 : // Evolve a weak dipole end.
248 : void pT2nextWeak( double pT2begDip, double pT2sel, TimeDipoleEnd& dip,
249 : Event& event);
250 :
251 : // Evolve a Hidden Valley dipole end.
252 : void pT2nextHV( double pT2begDip, double pT2sel, TimeDipoleEnd& dip,
253 : Event& );
254 :
255 : // Find kind of QCD ME correction.
256 : void findMEtype( Event& event, TimeDipoleEnd& dip);
257 :
258 : // Find type of particle; used by findMEtype.
259 : int findMEparticle( int id, bool isHiddenColour = false);
260 :
261 : // Find mixture of V and A in gamma/Z: energy- and flavour-dependent.
262 : double gammaZmix( Event& event, int iRes, int iDau1, int iDau2);
263 :
264 : // Set up to calculate QCD ME correction with calcMEcorr.
265 : double findMEcorr(TimeDipoleEnd* dip, Particle& rad, Particle& partner,
266 : Particle& emt, bool cutEdge = true);
267 :
268 : // Set up to calculate weak ME corrections for t-channel processes.
269 : double findMEcorrWeak(TimeDipoleEnd* dip, Vec4 rad, Vec4 rec,
270 : Vec4 emt, Vec4 p3, Vec4 p4, Vec4 radBef, Vec4 recBef);
271 :
272 : // Calculate value of QCD ME correction.
273 : double calcMEcorr( int kind, int combiIn, double mixIn, double x1,
274 : double x2, double r1, double r2, double r3 = 0., bool cutEdge = true);
275 :
276 : // Find coefficient of azimuthal asymmetry from gluon polarization.
277 : void findAsymPol( Event& event, TimeDipoleEnd* dip);
278 :
279 : // Rescatter: propagate dipole recoil to internal lines connecting systems.
280 : bool rescatterPropagateRecoil( Event& event, Vec4& pNew);
281 :
282 : // Properties stored for (some) global recoil schemes.
283 : // Vectors of event indices defining the hard process.
284 : vector<int> hardPartons;
285 : // Number of partons in current hard event, number of partons in Born-type
286 : // hard event (to distinguish between S and H), maximally allowed number of
287 : // global recoil branchings.
288 : int nHard, nFinalBorn, nMaxGlobalBranch;
289 : // Number of proposed splittings in hard scattering systems.
290 : map<int,int> nProposed;
291 : // Number of splittings with global recoil (currently only 1).
292 : int nGlobal, globalRecoilMode;
293 : // Switch to constrain recoiling system.
294 : bool limitMUQ;
295 :
296 : // Pointer to MergingHooks object for NLO merging.
297 : MergingHooks* mergingHooksPtr;
298 :
299 : };
300 :
301 : //==========================================================================
302 :
303 : } // end namespace Pythia8
304 :
305 : #endif // Pythia8_TimeShower_H
|
