Line data Source code
1 : // UserHooks.cc 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 : // 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 "Pythia8/PhaseSpace.h"
10 : #include "Pythia8/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 0 : return ((inEvent && sigmaProcessPtr->code() == 0
64 0 : && 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 0 : return ((inEvent && sigmaProcessPtr->code() == 0
125 0 : && 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 0 : workEvent.clear();
136 :
137 : // Loop through all partons. Beam particles should be copied.
138 0 : for (int i = 0; i < process.size(); ++i) {
139 : bool doCopy = false;
140 : bool isFinal = false;
141 0 : if (i < 3) doCopy = true;
142 :
143 : // Daughters of beams should normally be copied.
144 : else {
145 0 : int iMother = process[i].mother1();
146 0 : if (iMother == 1 || iMother == 2) doCopy = true;
147 :
148 : // Granddaughters of beams should normally be copied and are final.
149 0 : else if (iMother > 2) {
150 0 : int iGrandMother = process[iMother].mother1();
151 0 : if (iGrandMother == 1 || iGrandMother == 2) {
152 : doCopy = true;
153 : isFinal = true;
154 0 : }
155 0 : }
156 : }
157 :
158 : // Optionally non-final are not copied.
159 0 : if (finalOnly && !isFinal) doCopy = false;
160 :
161 : // Do copying and modify status/daughters of final.
162 0 : if (doCopy) {
163 0 : int iNew = workEvent.append( process[i]);
164 0 : if (isFinal) {
165 0 : workEvent[iNew].statusPos();
166 0 : workEvent[iNew].daughters( 0, 0);
167 : // When final only : no mothers; position in full event as daughters.
168 0 : if (finalOnly) {
169 0 : workEvent[iNew].mothers( 0, 0);
170 0 : workEvent[iNew].daughters( i, i);
171 0 : }
172 : }
173 0 : }
174 : }
175 :
176 0 : }
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 0 : workEvent.clear();
186 :
187 : // At the PartonLevel final partons are bookkept by subsystem.
188 0 : if (partonSystemsPtr->sizeSys() > 0) {
189 :
190 : // Find which subsystem to study.
191 : int iSys = 0;
192 0 : if (!isHardest) iSys = partonSystemsPtr->sizeSys() - 1;
193 :
194 : // Loop through all the final partons of the given subsystem.
195 0 : for (int i = 0; i < partonSystemsPtr->sizeOut(iSys); ++i) {
196 0 : int iOld = partonSystemsPtr->getOut( iSys, i);
197 :
198 : // Copy partons to work event.
199 0 : int iNew = workEvent.append( event[iOld]);
200 :
201 : // No mothers. Position in full event as daughters.
202 0 : workEvent[iNew].mothers( 0, 0);
203 0 : workEvent[iNew].daughters( iOld, iOld);
204 : }
205 :
206 : // At the ProcessLevel no subsystems have been defined.
207 0 : } else {
208 :
209 : // Loop through all partons, and copy all final ones.
210 0 : for (int iOld = 0; iOld < event.size(); ++iOld)
211 0 : if (event[iOld].isFinal()) {
212 0 : int iNew = workEvent.append( event[iOld]);
213 :
214 : // No mothers. Position in full event as daughters.
215 0 : workEvent[iNew].mothers( 0, 0);
216 0 : workEvent[iNew].daughters( iOld, iOld);
217 0 : }
218 : }
219 :
220 0 : }
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 0 : if (!isInit) {
235 :
236 : // Calculate pT0 as for multiparton interactions.
237 : // Fudge factor allows offset relative to MPI framework.
238 0 : double eCM = phaseSpacePtr->ecm();
239 0 : double pT0Ref = settingsPtr->parm("MultipartonInteractions:pT0Ref");
240 0 : double ecmRef = settingsPtr->parm("MultipartonInteractions:ecmRef");
241 0 : double ecmPow = settingsPtr->parm("MultipartonInteractions:ecmPow");
242 0 : double pT0 = pT0timesMPI * pT0Ref * pow(eCM / ecmRef, ecmPow);
243 0 : 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 0 : int alphaSnfmax = settingsPtr->mode("StandardModel:alphaSnfmax");
250 0 : if (useSameAlphaSasMPI) {
251 0 : alphaSvalue = settingsPtr->parm("MultipartonInteractions:alphaSvalue");
252 0 : alphaSorder = settingsPtr->mode("MultipartonInteractions:alphaSorder");
253 0 : } else {
254 0 : alphaSvalue = settingsPtr->parm("SigmaProcess:alphaSvalue");
255 0 : alphaSorder = settingsPtr->mode("SigmaProcess:alphaSorder");
256 : }
257 0 : alphaS.init( alphaSvalue, alphaSorder, alphaSnfmax, false);
258 :
259 : // Initialization finished.
260 0 : isInit = true;
261 0 : }
262 :
263 : // Only modify 2 -> 2 processes.
264 0 : int nFinal = sigmaProcessPtr->nFinal();
265 0 : if (nFinal != 2) return 1.;
266 :
267 : // pT scale of process. Weight pT^4 / (pT^2 + pT0^2)^2
268 0 : double pTHat = phaseSpacePtr->pTHat();
269 0 : double pT2 = pTHat * pTHat;
270 0 : double wt = pow2( pT2 / (pT20 + pT2) );
271 :
272 0 : if (numberAlphaS > 0) {
273 : // Renormalization scale and assumed alpha_strong.
274 0 : double Q2RenOld = sigmaProcessPtr->Q2Ren();
275 0 : double alphaSOld = sigmaProcessPtr->alphaSRen();
276 :
277 : // Reweight to new alpha_strong at new scale.
278 0 : double Q2RenNew = pT20 + Q2RenOld;
279 0 : double alphaSNew = alphaS.alphaS(Q2RenNew);
280 0 : wt *= pow( alphaSNew / alphaSOld, numberAlphaS);
281 0 : }
282 :
283 : // End weight calculation.
284 : return wt;
285 :
286 0 : }
287 :
288 :
289 : //==========================================================================
290 :
291 : } // end namespace Pythia8
|
