]>
Commit | Line | Data |
---|---|---|
63ba5337 | 1 | // UserHooks.cc is a part of the PYTHIA event generator. |
2 | // Copyright (C) 2012 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 | // Note: compilation crashes if PhaseSpace.h is moved to UserHooks.h. | |
9 | #include "PhaseSpace.h" | |
10 | #include "UserHooks.h" | |
11 | ||
12 | namespace Pythia8 { | |
13 | ||
14 | //========================================================================== | |
15 | ||
16 | // The UserHooks class. | |
17 | ||
18 | //-------------------------------------------------------------------------- | |
19 | ||
20 | // multiplySigmaBy allows the user to introduce a multiplicative factor | |
21 | // that modifies the cross section of a hard process. Since it is called | |
22 | // from before the event record is generated in full, the normal analysis | |
23 | // does not work. The code here provides a rather extensive summary of | |
24 | // which methods actually do work. It is a convenient starting point for | |
25 | // writing your own derived routine. | |
26 | ||
27 | double UserHooks::multiplySigmaBy( const SigmaProcess* sigmaProcessPtr, | |
28 | const PhaseSpace* phaseSpacePtr, bool inEvent) { | |
29 | ||
30 | // Process code, necessary when some to be treated differently. | |
31 | //int code = sigmaProcessPtr->code(); | |
32 | ||
33 | // Final multiplicity, i.e. whether 2 -> 1 or 2 -> 2. | |
34 | //int nFinal = sigmaProcessPtr->nFinal(); | |
35 | ||
36 | // Incoming x1 and x2 to the hard collision, and factorization scale. | |
37 | //double x1 = phaseSpacePtr->x1(); | |
38 | //double x2 = phaseSpacePtr->x2(); | |
39 | //double Q2Fac = sigmaProcessPtr->Q2Fac(); | |
40 | ||
41 | // Renormalization scale and assumed alpha_strong and alpha_EM. | |
42 | //double Q2Ren = sigmaProcessPtr->Q2Ren(); | |
43 | //double alphaS = sigmaProcessPtr->alphaSRen(); | |
44 | //double alphaEM = sigmaProcessPtr->alphaEMRen(); | |
45 | ||
46 | // Subprocess mass-square. | |
47 | //double sHat = phaseSpacePtr->sHat(); | |
48 | ||
49 | // Now methods only relevant for 2 -> 2. | |
50 | //if (nFinal == 2) { | |
51 | ||
52 | // Mandelstam variables and hard-process pT. | |
53 | //double tHat = phaseSpacePtr->tHat(); | |
54 | //double uHat = phaseSpacePtr->uHat(); | |
55 | //double pTHat = phaseSpacePtr->pTHat(); | |
56 | ||
57 | // Masses of the final-state particles. (Here 0 for light quarks.) | |
58 | //double m3 = sigmaProcessPtr->m(3); | |
59 | //double m4 = sigmaProcessPtr->m(4); | |
60 | //} | |
61 | ||
62 | // Dummy statement to avoid compiler warnings. | |
63 | return ((inEvent && sigmaProcessPtr->code() == 0 | |
64 | && phaseSpacePtr->sHat() < 0.) ? 0. : 1.); | |
65 | ||
66 | } | |
67 | ||
68 | //-------------------------------------------------------------------------- | |
69 | ||
70 | // biasSelectionBy allows the user to introduce a multiplicative factor | |
71 | // that modifies the cross section of a hard process. The event is assigned | |
72 | // a wegith that is the inverse of the selection bias, such that the | |
73 | // cross section is unchanged. Since it is called from before the | |
74 | // event record is generated in full, the normal analysis does not work. | |
75 | // The code here provides a rather extensive summary of which methods | |
76 | // actually do work. It is a convenient starting point for writing | |
77 | // your own derived routine. | |
78 | ||
79 | double UserHooks::biasSelectionBy( const SigmaProcess* sigmaProcessPtr, | |
80 | const PhaseSpace* phaseSpacePtr, bool inEvent) { | |
81 | ||
82 | // Process code, necessary when some to be treated differently. | |
83 | //int code = sigmaProcessPtr->code(); | |
84 | ||
85 | // Final multiplicity, i.e. whether 2 -> 1 or 2 -> 2. | |
86 | //int nFinal = sigmaProcessPtr->nFinal(); | |
87 | ||
88 | // Incoming x1 and x2 to the hard collision, and factorization scale. | |
89 | //double x1 = phaseSpacePtr->x1(); | |
90 | //double x2 = phaseSpacePtr->x2(); | |
91 | //double Q2Fac = sigmaProcessPtr->Q2Fac(); | |
92 | ||
93 | // Renormalization scale and assumed alpha_strong and alpha_EM. | |
94 | //double Q2Ren = sigmaProcessPtr->Q2Ren(); | |
95 | //double alphaS = sigmaProcessPtr->alphaSRen(); | |
96 | //double alphaEM = sigmaProcessPtr->alphaEMRen(); | |
97 | ||
98 | // Subprocess mass-square. | |
99 | //double sHat = phaseSpacePtr->sHat(); | |
100 | ||
101 | // Now methods only relevant for 2 -> 2. | |
102 | //if (nFinal == 2) { | |
103 | ||
104 | // Mandelstam variables and hard-process pT. | |
105 | //double tHat = phaseSpacePtr->tHat(); | |
106 | //double uHat = phaseSpacePtr->uHat(); | |
107 | //double pTHat = phaseSpacePtr->pTHat(); | |
108 | ||
109 | // Masses of the final-state particles. (Here 0 for light quarks.) | |
110 | //double m3 = sigmaProcessPtr->m(3); | |
111 | //double m4 = sigmaProcessPtr->m(4); | |
112 | //} | |
113 | ||
114 | // Insert here your calculation of the selection bias. | |
115 | // Here illustrated by a weighting up of events at high pT. | |
116 | //selBias = pow4(phaseSpacePtr->pTHat()); | |
117 | ||
118 | // Return the selBias weight. | |
119 | // Warning: if you use another variable than selBias | |
120 | // the compensating weight will not be set correctly. | |
121 | //return selBias; | |
122 | ||
123 | // Dummy statement to avoid compiler warnings. | |
124 | return ((inEvent && sigmaProcessPtr->code() == 0 | |
125 | && phaseSpacePtr->sHat() < 0.) ? 0. : 1.); | |
126 | } | |
127 | ||
128 | //-------------------------------------------------------------------------- | |
129 | ||
130 | // omitResonanceDecays omits resonance decay chains from process record. | |
131 | ||
132 | void UserHooks::omitResonanceDecays(const Event& process, bool finalOnly) { | |
133 | ||
134 | // Reset work event to be empty | |
135 | workEvent.clear(); | |
136 | ||
137 | // Loop through all partons. Beam particles should be copied. | |
138 | for (int i = 0; i < process.size(); ++i) { | |
139 | bool doCopy = false; | |
140 | bool isFinal = false; | |
141 | if (i < 3) doCopy = true; | |
142 | ||
143 | // Daughters of beams should normally be copied. | |
144 | else { | |
145 | int iMother = process[i].mother1(); | |
146 | if (iMother == 1 || iMother == 2) doCopy = true; | |
147 | ||
148 | // Granddaughters of beams should normally be copied and are final. | |
149 | else if (iMother > 2) { | |
150 | int iGrandMother = process[iMother].mother1(); | |
151 | if (iGrandMother == 1 || iGrandMother == 2) { | |
152 | doCopy = true; | |
153 | isFinal = true; | |
154 | } | |
155 | } | |
156 | } | |
157 | ||
158 | // Optionally non-final are not copied. | |
159 | if (finalOnly && !isFinal) doCopy = false; | |
160 | ||
161 | // Do copying and modify status/daughters of final. | |
162 | if (doCopy) { | |
163 | int iNew = workEvent.append( process[i]); | |
164 | if (isFinal) { | |
165 | workEvent[iNew].statusPos(); | |
166 | workEvent[iNew].daughters( 0, 0); | |
167 | // When final only : no mothers; position in full event as daughters. | |
168 | if (finalOnly) { | |
169 | workEvent[iNew].mothers( 0, 0); | |
170 | workEvent[iNew].daughters( i, i); | |
171 | } | |
172 | } | |
173 | } | |
174 | } | |
175 | ||
176 | } | |
177 | ||
178 | //-------------------------------------------------------------------------- | |
179 | ||
180 | // subEvent extracts currently resolved partons in the hard process. | |
181 | ||
182 | void UserHooks::subEvent(const Event& event, bool isHardest) { | |
183 | ||
184 | // Reset work event to be empty. | |
185 | workEvent.clear(); | |
186 | ||
187 | // At the PartonLevel final partons are bookkept by subsystem. | |
188 | if (partonSystemsPtr->sizeSys() > 0) { | |
189 | ||
190 | // Find which subsystem to study. | |
191 | int iSys = 0; | |
192 | if (!isHardest) iSys = partonSystemsPtr->sizeSys() - 1; | |
193 | ||
194 | // Loop through all the final partons of the given subsystem. | |
195 | for (int i = 0; i < partonSystemsPtr->sizeOut(iSys); ++i) { | |
196 | int iOld = partonSystemsPtr->getOut( iSys, i); | |
197 | ||
198 | // Copy partons to work event. | |
199 | int iNew = workEvent.append( event[iOld]); | |
200 | ||
201 | // No mothers. Position in full event as daughters. | |
202 | workEvent[iNew].mothers( 0, 0); | |
203 | workEvent[iNew].daughters( iOld, iOld); | |
204 | } | |
205 | ||
206 | // At the ProcessLevel no subsystems have been defined. | |
207 | } else { | |
208 | ||
209 | // Loop through all partons, and copy all final ones. | |
210 | for (int iOld = 0; iOld < event.size(); ++iOld) | |
211 | if (event[iOld].isFinal()) { | |
212 | int iNew = workEvent.append( event[iOld]); | |
213 | ||
214 | // No mothers. Position in full event as daughters. | |
215 | workEvent[iNew].mothers( 0, 0); | |
216 | workEvent[iNew].daughters( iOld, iOld); | |
217 | } | |
218 | } | |
219 | ||
220 | } | |
221 | ||
222 | //========================================================================== | |
223 | ||
224 | // The SuppressSmallPT class, derived from UserHooks. | |
225 | ||
226 | //-------------------------------------------------------------------------- | |
227 | ||
228 | // Modify event weight at the trial level, before selection. | |
229 | ||
230 | double SuppressSmallPT::multiplySigmaBy( const SigmaProcess* sigmaProcessPtr, | |
231 | const PhaseSpace* phaseSpacePtr, bool ) { | |
232 | ||
233 | // Need to initialize first time this method is called. | |
234 | if (!isInit) { | |
235 | ||
236 | // Calculate pT0 as for multiparton interactions. | |
237 | // Fudge factor allows offset relative to MPI framework. | |
238 | double eCM = phaseSpacePtr->ecm(); | |
239 | double pT0Ref = settingsPtr->parm("MultipartonInteractions:pT0Ref"); | |
240 | double ecmRef = settingsPtr->parm("MultipartonInteractions:ecmRef"); | |
241 | double ecmPow = settingsPtr->parm("MultipartonInteractions:ecmPow"); | |
242 | double pT0 = pT0timesMPI * pT0Ref * pow(eCM / ecmRef, ecmPow); | |
243 | pT20 = pT0 * pT0; | |
244 | ||
245 | // Initialize alpha_strong object as for multiparton interactions, | |
246 | // alternatively as for hard processes. | |
247 | double alphaSvalue; | |
248 | int alphaSorder; | |
249 | if (useSameAlphaSasMPI) { | |
250 | alphaSvalue = settingsPtr->parm("MultipartonInteractions:alphaSvalue"); | |
251 | alphaSorder = settingsPtr->mode("MultipartonInteractions:alphaSorder"); | |
252 | } else { | |
253 | alphaSvalue = settingsPtr->parm("SigmaProcess:alphaSvalue"); | |
254 | alphaSorder = settingsPtr->mode("SigmaProcess:alphaSorder"); | |
255 | } | |
256 | alphaS.init( alphaSvalue, alphaSorder); | |
257 | ||
258 | // Initialization finished. | |
259 | isInit = true; | |
260 | } | |
261 | ||
262 | // Only modify 2 -> 2 processes. | |
263 | int nFinal = sigmaProcessPtr->nFinal(); | |
264 | if (nFinal != 2) return 1.; | |
265 | ||
266 | // pT scale of process. Weight pT^4 / (pT^2 + pT0^2)^2 | |
267 | double pTHat = phaseSpacePtr->pTHat(); | |
268 | double pT2 = pTHat * pTHat; | |
269 | double wt = pow2( pT2 / (pT20 + pT2) ); | |
270 | ||
271 | if (numberAlphaS > 0) { | |
272 | // Renormalization scale and assumed alpha_strong. | |
273 | double Q2RenOld = sigmaProcessPtr->Q2Ren(); | |
274 | double alphaSOld = sigmaProcessPtr->alphaSRen(); | |
275 | ||
276 | // Reweight to new alpha_strong at new scale. | |
277 | double Q2RenNew = pT20 + Q2RenOld; | |
278 | double alphaSNew = alphaS.alphaS(Q2RenNew); | |
279 | wt *= pow( alphaSNew / alphaSOld, numberAlphaS); | |
280 | } | |
281 | ||
282 | // End weight calculation. | |
283 | return wt; | |
284 | ||
285 | } | |
286 | ||
287 | ||
288 | //========================================================================== | |
289 | ||
290 | } // end namespace Pythia8 |