]>
Commit | Line | Data |
---|---|---|
9419eeef | 1 | // UserHooks.cc is a part of the PYTHIA event generator. |
2 | // Copyright (C) 2010 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 | // omitResonanceDecays omits resonance decay chains from process record. | |
69 | ||
70 | void UserHooks::omitResonanceDecays(const Event& process) { | |
71 | ||
72 | // Reset work event to be empty | |
73 | workEvent.clear(); | |
74 | ||
75 | // Loop through all partons. Beam particles should be copied. | |
76 | for (int i = 0; i < process.size(); ++i) { | |
77 | bool doCopy = false; | |
78 | bool isFinal = false; | |
79 | if (i < 3) doCopy = true; | |
80 | ||
81 | // Daughters of beams should be copied. | |
82 | else { | |
83 | int iMother = process[i].mother1(); | |
84 | if (iMother == 1 || iMother == 2) doCopy = true; | |
85 | ||
86 | // Granddaughters of beams should be copied and are final. | |
87 | else if (iMother > 2) { | |
88 | int iGrandMother = process[iMother].mother1(); | |
89 | if (iGrandMother == 1 || iGrandMother == 2) { | |
90 | doCopy = true; | |
91 | isFinal = true; | |
92 | } | |
93 | } | |
94 | } | |
95 | ||
96 | // Do copying and modify status/daughters of final. | |
97 | if (doCopy) { | |
98 | int iNew = workEvent.append( process[i]); | |
99 | if (isFinal) { | |
100 | workEvent[iNew].statusPos(); | |
101 | workEvent[iNew].daughters( 0, 0); | |
102 | } | |
103 | } | |
104 | } | |
105 | ||
106 | } | |
107 | ||
108 | //-------------------------------------------------------------------------- | |
109 | ||
110 | // subEvent extracts currently resolved partons in the hard process. | |
111 | ||
112 | void UserHooks::subEvent(const Event& event, bool isHardest) { | |
113 | ||
114 | // Reset work event to be empty. | |
115 | workEvent.clear(); | |
116 | ||
117 | // Find which subsystem to study. | |
118 | int iSys = 0; | |
119 | if (!isHardest) iSys = partonSystemsPtr->sizeSys() - 1; | |
120 | ||
121 | // Loop through all the final partons of the given subsystem. | |
122 | for (int i = 0; i < partonSystemsPtr->sizeOut(iSys); ++i) { | |
123 | int iOld = partonSystemsPtr->getOut( iSys, i); | |
124 | ||
125 | // Copy partons to work event. | |
126 | int iNew = workEvent.append( event[iOld]); | |
127 | ||
128 | // No mothers. Position in full event as daughters. | |
129 | workEvent[iNew].mothers( 0, 0); | |
130 | workEvent[iNew].daughters( iOld, iOld); | |
131 | } | |
132 | ||
133 | } | |
134 | ||
135 | //========================================================================== | |
136 | ||
137 | // The SuppressSmallPT class, derived from UserHooks. | |
138 | ||
139 | //-------------------------------------------------------------------------- | |
140 | ||
141 | // Modify event weight at the trial level, before selection. | |
142 | ||
143 | double SuppressSmallPT::multiplySigmaBy( const SigmaProcess* sigmaProcessPtr, | |
144 | const PhaseSpace* phaseSpacePtr, bool ) { | |
145 | ||
146 | // Need to initialize first time this method is called. | |
147 | if (!isInit) { | |
148 | ||
149 | // Calculate pT0 as for multiple interactions. | |
150 | // Fudge factor allows offset relative to MI framework. | |
151 | double eCM = phaseSpacePtr->ecm(); | |
152 | double pT0Ref = settingsPtr->parm("MultipleInteractions:pT0Ref"); | |
153 | double ecmRef = settingsPtr->parm("MultipleInteractions:ecmRef"); | |
154 | double ecmPow = settingsPtr->parm("MultipleInteractions:ecmPow"); | |
155 | double pT0 = pT0timesMI * pT0Ref * pow(eCM / ecmRef, ecmPow); | |
156 | pT20 = pT0 * pT0; | |
157 | ||
158 | // Initialize alpha_strong object as for multiple interactions, | |
159 | // alternatively as for hard processes. | |
160 | double alphaSvalue; | |
161 | int alphaSorder; | |
162 | if (useSameAlphaSasMI) { | |
163 | alphaSvalue = settingsPtr->parm("MultipleInteractions:alphaSvalue"); | |
164 | alphaSorder = settingsPtr->mode("MultipleInteractions:alphaSorder"); | |
165 | } else { | |
166 | alphaSvalue = settingsPtr->parm("SigmaProcess:alphaSvalue"); | |
167 | alphaSorder = settingsPtr->mode("SigmaProcess:alphaSorder"); | |
168 | } | |
169 | alphaS.init( alphaSvalue, alphaSorder); | |
170 | ||
171 | // Initialization finished. | |
172 | isInit = true; | |
173 | } | |
174 | ||
175 | // Only modify 2 -> 2 processes. | |
176 | int nFinal = sigmaProcessPtr->nFinal(); | |
177 | if (nFinal != 2) return 1.; | |
178 | ||
179 | // pT scale of process. Weight pT^4 / (pT^2 + pT0^2)^2 | |
180 | double pTHat = phaseSpacePtr->pTHat(); | |
181 | double pT2 = pTHat * pTHat; | |
182 | double wt = pow2( pT2 / (pT20 + pT2) ); | |
183 | ||
184 | if (numberAlphaS > 0) { | |
185 | // Renormalization scale and assumed alpha_strong. | |
186 | double Q2RenOld = sigmaProcessPtr->Q2Ren(); | |
187 | double alphaSOld = sigmaProcessPtr->alphaSRen(); | |
188 | ||
189 | // Reweight to new alpha_strong at new scale. | |
190 | double Q2RenNew = pT20 + Q2RenOld; | |
191 | double alphaSNew = alphaS.alphaS(Q2RenNew); | |
192 | wt *= pow( alphaSNew / alphaSOld, numberAlphaS); | |
193 | } | |
194 | ||
195 | // End weight calculation. | |
196 | return wt; | |
197 | ||
198 | } | |
199 | ||
200 | ||
201 | //========================================================================== | |
202 | ||
203 | } // end namespace Pythia8 |