]>
Commit | Line | Data |
---|---|---|
5ad4eb21 | 1 | // UserHooks.cc is a part of the PYTHIA event generator. |
2 | // Copyright (C) 2008 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 | // Function definitions (not found in the header) for the UserHooks class. | |
7 | ||
8 | #include "UserHooks.h" | |
9 | ||
10 | namespace Pythia8 { | |
11 | ||
12 | //************************************************************************** | |
13 | ||
14 | // The UserHooks class. | |
15 | ||
16 | //********* | |
17 | ||
18 | // multiplySigmaBy allows the user to introduce a multiplicative factor | |
19 | // that modifies the cross section of a hard process. Since it is called | |
20 | // from before the event record is generated in full, the normal analysis | |
21 | // does not work. The code here provides a rather extensive summary of | |
22 | // which methods actually do work. It is a convenient starting point for | |
23 | // writing your own derived routine. | |
24 | ||
25 | double UserHooks::multiplySigmaBy( const SigmaProcess* sigmaProcessPtr, | |
26 | const PhaseSpace* phaseSpacePtr, bool inEvent) { | |
27 | ||
28 | // Process code, necessary when some to be treated differently. | |
29 | //int code = sigmaProcessPtr->code(); | |
30 | ||
31 | // Final multiplicity, i.e. whether 2 -> 1 or 2 -> 2. | |
32 | //int nFinal = sigmaProcessPtr->nFinal(); | |
33 | ||
34 | // Incoming x1 and x2 to the hard collision, and factorization scale. | |
35 | //double x1 = phaseSpacePtr->x1(); | |
36 | //double x2 = phaseSpacePtr->x2(); | |
37 | //double Q2Fac = sigmaProcessPtr->Q2Fac(); | |
38 | ||
39 | // Renormalization scale and assumed alpha_strong and alpha_EM. | |
40 | //double Q2Ren = sigmaProcessPtr->Q2Ren(); | |
41 | //double alphaS = sigmaProcessPtr->alphaSRen(); | |
42 | //double alphaEM = sigmaProcessPtr->alphaEMRen(); | |
43 | ||
44 | // Subprocess mass-square. | |
45 | //double sHat = phaseSpacePtr->sHat(); | |
46 | ||
47 | // Now methods only relevant for 2 -> 2. | |
48 | //if (nFinal == 2) { | |
49 | ||
50 | // Mandelstam variables and hard-process pT. | |
51 | //double tHat = phaseSpacePtr->tHat(); | |
52 | //double uHat = phaseSpacePtr->uHat(); | |
53 | //double pTHat = phaseSpacePtr->pTHat(); | |
54 | ||
55 | // Masses of the final-state particles. (Here 0 for light quarks.) | |
56 | //double m3 = sigmaProcessPtr->m(3); | |
57 | //double m4 = sigmaProcessPtr->m(4); | |
58 | //} | |
59 | ||
60 | // Dummy statement to avoid compiler warnings. | |
61 | return ((inEvent && sigmaProcessPtr->code() == 0 | |
62 | && phaseSpacePtr->sHat() < 0.) ? 0. : 1.); | |
63 | ||
64 | } | |
65 | ||
66 | //********* | |
67 | ||
68 | // subEvent extracts currently resolved partons in the hard process. | |
69 | ||
70 | void UserHooks::subEvent(const Event& event, bool isHardest) { | |
71 | ||
72 | // Reset work event to be empty. | |
73 | workEvent.clear(); | |
74 | ||
75 | // Find which subsystem to study. | |
76 | int iSys = 0; | |
77 | if (!isHardest) iSys = event.sizeSystems() - 1; | |
78 | ||
79 | // Loop through all the final partons of the given subsystem. | |
80 | for (int i = 2; i < event.sizeSystem(iSys); ++i) { | |
81 | int iOld = event.getInSystem( iSys, i); | |
82 | ||
83 | // Copy partons to work event. | |
84 | int iNew = workEvent.append( event[iOld]); | |
85 | ||
86 | // No mothers. Position in full event as daughters. | |
87 | workEvent[iNew].mothers( 0, 0); | |
88 | workEvent[iNew].daughters( iOld, iOld); | |
89 | } | |
90 | ||
91 | } | |
92 | ||
93 | //************************************************************************** | |
94 | ||
95 | // The SuppressSmallPT class, derived from UserHooks. | |
96 | ||
97 | //********* | |
98 | ||
99 | // Modify event weight at the trial level, before selection. | |
100 | ||
101 | double SuppressSmallPT::multiplySigmaBy( const SigmaProcess* sigmaProcessPtr, | |
102 | const PhaseSpace* phaseSpacePtr, bool ) { | |
103 | ||
104 | // Need to initialize first time this method is called. | |
105 | if (!isInit) { | |
106 | ||
107 | // Calculate pT0 as for multiple interactions. | |
108 | // Fudge factor allows offset relative to MI framework. | |
109 | double eCM = phaseSpacePtr->ecm(); | |
110 | double pT0Ref = Settings::parm("MultipleInteractions:pT0Ref"); | |
111 | double ecmRef = Settings::parm("MultipleInteractions:ecmRef"); | |
112 | double ecmPow = Settings::parm("MultipleInteractions:ecmPow"); | |
113 | double pT0 = pT0timesMI * pT0Ref * pow(eCM / ecmRef, ecmPow); | |
114 | pT20 = pT0 * pT0; | |
115 | ||
116 | // Initialize alpha_strong object as for multiple interactions, | |
117 | // alternatively as for hard processes. | |
118 | double alphaSvalue; | |
119 | int alphaSorder; | |
120 | if (useSameAlphaSasMI) { | |
121 | alphaSvalue = Settings::parm("MultipleInteractions:alphaSvalue"); | |
122 | alphaSorder = Settings::mode("MultipleInteractions:alphaSorder"); | |
123 | } else { | |
124 | alphaSvalue = Settings::parm("SigmaProcess:alphaSvalue"); | |
125 | alphaSorder = Settings::mode("SigmaProcess:alphaSorder"); | |
126 | } | |
127 | alphaS.init( alphaSvalue, alphaSorder); | |
128 | ||
129 | // Initialization finished. | |
130 | isInit = true; | |
131 | } | |
132 | ||
133 | // Only modify 2 -> 2 processes. | |
134 | int nFinal = sigmaProcessPtr->nFinal(); | |
135 | if (nFinal != 2) return 1.; | |
136 | ||
137 | // pT scale of process. Weight pT^4 / (pT^2 + pT0^2)^2 | |
138 | double pTHat = phaseSpacePtr->pTHat(); | |
139 | double pT2 = pTHat * pTHat; | |
140 | double wt = pow2( pT2 / (pT20 + pT2) ); | |
141 | ||
142 | if (numberAlphaS > 0) { | |
143 | // Renormalization scale and assumed alpha_strong. | |
144 | double Q2RenOld = sigmaProcessPtr->Q2Ren(); | |
145 | double alphaSOld = sigmaProcessPtr->alphaSRen(); | |
146 | ||
147 | // Reweight to new alpha_strong at new scale. | |
148 | double Q2RenNew = pT20 + Q2RenOld; | |
149 | double alphaSNew = alphaS.alphaS(Q2RenNew); | |
150 | wt *= pow( alphaSNew / alphaSOld, numberAlphaS); | |
151 | } | |
152 | ||
153 | // End weight calculation. | |
154 | return wt; | |
155 | ||
156 | } | |
157 | ||
158 | ||
159 | //************************************************************************** | |
160 | ||
161 | } // end namespace Pythia8 |