Update to pythi8.170
[u/mrichter/AliRoot.git] / PYTHIA8 / pythia8170 / src / UserHooks.cxx
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