StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
UserHooks.cc
1 // UserHooks.cc is a part of the PYTHIA event generator.
2 // Copyright (C) 2020 Torbjorn Sjostrand.
3 // PYTHIA is licenced under the GNU GPL v2 or later, 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  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 = parm("MultipartonInteractions:pT0Ref");
240  double ecmRef = parm("MultipartonInteractions:ecmRef");
241  double ecmPow = 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  int alphaSnfmax = mode("StandardModel:alphaSnfmax");
250  if (useSameAlphaSasMPI) {
251  alphaSvalue = parm("MultipartonInteractions:alphaSvalue");
252  alphaSorder = mode("MultipartonInteractions:alphaSorder");
253  } else {
254  alphaSvalue = parm("SigmaProcess:alphaSvalue");
255  alphaSorder = mode("SigmaProcess:alphaSorder");
256  }
257  alphaS.init( alphaSvalue, alphaSorder, alphaSnfmax, false);
258 
259  // Initialization finished.
260  isInit = true;
261  }
262 
263  // Only modify 2 -> 2 processes.
264  int nFinal = sigmaProcessPtr->nFinal();
265  if (nFinal != 2) return 1.;
266 
267  // pT scale of process. Weight pT^4 / (pT^2 + pT0^2)^2
268  double pTHat = phaseSpacePtr->pTHat();
269  double pT2 = pTHat * pTHat;
270  double wt = pow2( pT2 / (pT20 + pT2) );
271 
272  if (numberAlphaS > 0) {
273  // Renormalization scale and assumed alpha_strong.
274  double Q2RenOld = sigmaProcessPtr->Q2Ren();
275  double alphaSOld = sigmaProcessPtr->alphaSRen();
276 
277  // Reweight to new alpha_strong at new scale.
278  double Q2RenNew = pT20 + Q2RenOld;
279  double alphaSNew = alphaS.alphaS(Q2RenNew);
280  wt *= pow( alphaSNew / alphaSOld, numberAlphaS);
281  }
282 
283  // End weight calculation.
284  return wt;
285 
286 }
287 
288 
289 //==========================================================================
290 
291 } // end namespace Pythia8
Definition: AgUStep.h:26