1 | // HiddenValleyFragmentation.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 |
---|
7 | // HiddenValleyFragmentation class and its helper classes. |
---|
8 | |
---|
9 | #include "HiddenValleyFragmentation.h" |
---|
10 | |
---|
11 | namespace Pythia8 { |
---|
12 | |
---|
13 | //========================================================================== |
---|
14 | |
---|
15 | // The HVStringFlav class is used to select HV-quark and HV-hadron flavours. |
---|
16 | |
---|
17 | //-------------------------------------------------------------------------- |
---|
18 | |
---|
19 | // Initialize data members of the flavour generation. |
---|
20 | |
---|
21 | void HVStringFlav::init(Settings& settings, Rndm* rndmPtrIn) { |
---|
22 | |
---|
23 | // Save pointer. |
---|
24 | rndmPtr = rndmPtrIn; |
---|
25 | |
---|
26 | // Read in data from Settings. |
---|
27 | nFlav = settings.mode("HiddenValley:nFlav"); |
---|
28 | probVector = settings.parm("HiddenValley:probVector"); |
---|
29 | |
---|
30 | } |
---|
31 | |
---|
32 | //-------------------------------------------------------------------------- |
---|
33 | |
---|
34 | // Pick a new HV-flavour given an incoming one. |
---|
35 | |
---|
36 | FlavContainer HVStringFlav::pick(FlavContainer& flavOld) { |
---|
37 | |
---|
38 | // Initial values for new flavour. |
---|
39 | FlavContainer flavNew; |
---|
40 | flavNew.rank = flavOld.rank + 1; |
---|
41 | |
---|
42 | // Pick new HV-flavour at random; keep track of sign. |
---|
43 | flavNew.id = 4900100 + min( 1 + int(nFlav * rndmPtr->flat()), nFlav); |
---|
44 | if (flavOld.id > 0) flavNew.id = -flavNew.id; |
---|
45 | |
---|
46 | // Done. |
---|
47 | return flavNew; |
---|
48 | |
---|
49 | } |
---|
50 | |
---|
51 | //-------------------------------------------------------------------------- |
---|
52 | |
---|
53 | // Combine two HV-flavours to produce an HV-hadron. |
---|
54 | // This is simplified procedure, assuming only two HV mesons defined. |
---|
55 | |
---|
56 | int HVStringFlav::combine(FlavContainer& flav1, FlavContainer& flav2) { |
---|
57 | |
---|
58 | // Positive and negative flavour. Note that with kinetic mixing |
---|
59 | // the Fv are really intended to represent qv, so remap. |
---|
60 | int idMeson = 0; |
---|
61 | int idPos = max( flav1.id, flav2.id) - 4900000; |
---|
62 | int idNeg = -min( flav1.id, flav2.id) - 4900000; |
---|
63 | if (idPos < 20) idPos = 101; |
---|
64 | if (idNeg < 20) idNeg = 101; |
---|
65 | |
---|
66 | // Pick HV-meson code, spin either 0 or 1. |
---|
67 | if (idNeg == idPos) idMeson = 4900111; |
---|
68 | else if (idPos > idNeg) idMeson = 4900211; |
---|
69 | else idMeson = -4900211; |
---|
70 | if (rndmPtr->flat() < probVector) idMeson += ((idMeson > 0) ? 2 : -2); |
---|
71 | |
---|
72 | // Done. |
---|
73 | return idMeson; |
---|
74 | |
---|
75 | } |
---|
76 | |
---|
77 | //========================================================================== |
---|
78 | |
---|
79 | // The HVStringPT class is used to select pT in HV fragmentation. |
---|
80 | |
---|
81 | //-------------------------------------------------------------------------- |
---|
82 | |
---|
83 | // Initialize data members of the string pT selection. |
---|
84 | |
---|
85 | void HVStringPT::init(Settings& settings, ParticleData& particleData, |
---|
86 | Rndm* rndmPtrIn) { |
---|
87 | |
---|
88 | // Save pointer. |
---|
89 | rndmPtr = rndmPtrIn; |
---|
90 | |
---|
91 | // Parameter of the pT width. No enhancement, since this is finetuning. |
---|
92 | double sigmamqv = settings.parm("HiddenValley:sigmamqv"); |
---|
93 | double sigma = sigmamqv * particleData.m0( 4900101); |
---|
94 | sigmaQ = sigma / sqrt(2.); |
---|
95 | enhancedFraction = 0.; |
---|
96 | enhancedWidth = 0.; |
---|
97 | |
---|
98 | // Parameter for pT suppression in MiniStringFragmentation. |
---|
99 | sigma2Had = 2. * pow2( max( SIGMAMIN, sigma) ); |
---|
100 | |
---|
101 | } |
---|
102 | |
---|
103 | //========================================================================== |
---|
104 | |
---|
105 | // The HVStringZ class is used to select z in HV fragmentation. |
---|
106 | |
---|
107 | //-------------------------------------------------------------------------- |
---|
108 | |
---|
109 | // Initialize data members of the string z selection. |
---|
110 | |
---|
111 | void HVStringZ::init(Settings& settings, ParticleData& particleData, |
---|
112 | Rndm* rndmPtrIn) { |
---|
113 | |
---|
114 | // Save pointer. |
---|
115 | rndmPtr = rndmPtrIn; |
---|
116 | |
---|
117 | // Paramaters of Lund/Bowler symmetric fragmentation function. |
---|
118 | aLund = settings.parm("HiddenValley:aLund"); |
---|
119 | bmqv2 = settings.parm("HiddenValley:bmqv2"); |
---|
120 | rFactqv = settings.parm("HiddenValley:rFactqv"); |
---|
121 | |
---|
122 | // Use qv mass to set scale of bEff = b * m^2; |
---|
123 | mqv2 = pow2( particleData.m0( 4900101) ); |
---|
124 | bLund = bmqv2 / mqv2; |
---|
125 | |
---|
126 | // Mass of qv meson used to set stop scale for fragmentation iteration. |
---|
127 | mhvMeson = particleData.m0( 4900111); |
---|
128 | |
---|
129 | } |
---|
130 | |
---|
131 | //-------------------------------------------------------------------------- |
---|
132 | |
---|
133 | // Generate the fraction z that the next hadron will take using Lund/Bowler. |
---|
134 | |
---|
135 | double HVStringZ::zFrag( int , int , double mT2) { |
---|
136 | |
---|
137 | // Shape parameters of Lund symmetric fragmentation function. |
---|
138 | double bShape = bLund * mT2; |
---|
139 | double cShape = 1. + rFactqv * bmqv2; |
---|
140 | return zLund( aLund, bShape, cShape); |
---|
141 | |
---|
142 | } |
---|
143 | |
---|
144 | //========================================================================== |
---|
145 | |
---|
146 | // The HiddenValleyFragmentation class. |
---|
147 | |
---|
148 | //-------------------------------------------------------------------------- |
---|
149 | |
---|
150 | // Initialize and save pointers. |
---|
151 | |
---|
152 | bool HiddenValleyFragmentation::init(Info* infoPtrIn, Settings& settings, |
---|
153 | ParticleData* particleDataPtrIn, Rndm* rndmPtrIn) { |
---|
154 | |
---|
155 | // Save pointers. |
---|
156 | infoPtr = infoPtrIn; |
---|
157 | particleDataPtr = particleDataPtrIn; |
---|
158 | rndmPtr = rndmPtrIn; |
---|
159 | |
---|
160 | // Check whether Hidden Valley fragmentation switched on, and SU(N). |
---|
161 | doHVfrag = settings.flag("HiddenValley:fragment"); |
---|
162 | if (settings.mode("HiddenValley:Ngauge") < 2) doHVfrag = false; |
---|
163 | if (!doHVfrag) return false; |
---|
164 | |
---|
165 | // Several copies of qv may be needed. Taken to have same mass. |
---|
166 | nFlav = settings.mode("HiddenValley:nFlav"); |
---|
167 | if (nFlav > 1) { |
---|
168 | int spinType = particleDataPtr->spinType(4900101); |
---|
169 | double m0 = particleDataPtr->m0(4900101); |
---|
170 | for (int iFlav = 2; iFlav <= nFlav; ++iFlav) |
---|
171 | particleDataPtr->addParticle( 4900100 + iFlav, "qv", "qvbar", |
---|
172 | spinType, 0, 0, m0); |
---|
173 | } |
---|
174 | |
---|
175 | // Hidden Valley meson mass used to choose hadronization mode. |
---|
176 | mhvMeson = particleDataPtr->m0(4900111); |
---|
177 | |
---|
178 | // Initialize the hvEvent instance of an event record. |
---|
179 | hvEvent.init( "(Hidden Valley fragmentation)", particleDataPtr); |
---|
180 | |
---|
181 | // Create HVStringFlav instance for HV-flavour selection. |
---|
182 | hvFlavSelPtr = new HVStringFlav(); |
---|
183 | hvFlavSelPtr->init( settings, rndmPtr); |
---|
184 | |
---|
185 | // Create HVStringPT instance for pT selection in HV fragmentation. |
---|
186 | hvPTSelPtr = new HVStringPT(); |
---|
187 | hvPTSelPtr->init( settings, *particleDataPtr, rndmPtr); |
---|
188 | |
---|
189 | // Create HVStringZ instance for z selection in HV fragmentation. |
---|
190 | hvZSelPtr = new HVStringZ(); |
---|
191 | hvZSelPtr->init( settings, *particleDataPtr, rndmPtr); |
---|
192 | |
---|
193 | // Initialize auxiliary administrative class. |
---|
194 | hvColConfig.init(infoPtr, settings, hvFlavSelPtr); |
---|
195 | |
---|
196 | // Initialize HV-string and HV-ministring fragmentation. |
---|
197 | hvStringFrag.init(infoPtr, settings, particleDataPtr, rndmPtr, |
---|
198 | hvFlavSelPtr, hvPTSelPtr, hvZSelPtr); |
---|
199 | hvMinistringFrag.init(infoPtr, settings, particleDataPtr, rndmPtr, |
---|
200 | hvFlavSelPtr, hvPTSelPtr, hvZSelPtr); |
---|
201 | |
---|
202 | // Done. |
---|
203 | return true; |
---|
204 | |
---|
205 | } |
---|
206 | |
---|
207 | //-------------------------------------------------------------------------- |
---|
208 | |
---|
209 | // Perform the fragmentation. |
---|
210 | |
---|
211 | bool HiddenValleyFragmentation::fragment(Event& event) { |
---|
212 | |
---|
213 | // Reset containers for next event. |
---|
214 | hvEvent.reset(); |
---|
215 | hvColConfig.clear(); |
---|
216 | ihvParton.resize(0); |
---|
217 | |
---|
218 | // Extract HV-particles from event to hvEvent. Assign HV-colours. |
---|
219 | // Done if no HV-particles found. |
---|
220 | if (!extractHVevent(event)) return true; |
---|
221 | |
---|
222 | // Store found string system. Analyze its properties. |
---|
223 | if (!hvColConfig.insert(ihvParton, hvEvent)) return false; |
---|
224 | |
---|
225 | // Collect sequentially all partons in the HV subsystem. |
---|
226 | // Copy also if already in order, or else history tracing may fail. |
---|
227 | hvColConfig.collect(0, hvEvent, false); |
---|
228 | |
---|
229 | // Mass used to decide how to fragment system. |
---|
230 | mSys = hvColConfig[0].mass; |
---|
231 | |
---|
232 | // HV-string fragmentation when enough mass to produce >= 3 HV-mesons. |
---|
233 | if (mSys > 3.5 * mhvMeson) { |
---|
234 | if (!hvStringFrag.fragment( 0, hvColConfig, hvEvent)) return false; |
---|
235 | |
---|
236 | // HV-ministring fragmentation when enough mass to produce 2 HV-mesons. |
---|
237 | } else if (mSys > 2.1 * mhvMeson) { |
---|
238 | if (!hvMinistringFrag.fragment( 0, hvColConfig, hvEvent, true)) |
---|
239 | return false; |
---|
240 | |
---|
241 | // If only enough mass for one HV-meson assume HV-glueballs emitted. |
---|
242 | } else if (!collapseToMeson()) return false; |
---|
243 | |
---|
244 | // Insert HV particles from hvEvent to event. |
---|
245 | insertHVevent(event); |
---|
246 | |
---|
247 | // Done. |
---|
248 | return true; |
---|
249 | |
---|
250 | } |
---|
251 | |
---|
252 | //-------------------------------------------------------------------------- |
---|
253 | |
---|
254 | // Extract HV-particles from event to hvEvent. Assign HV-colours. |
---|
255 | |
---|
256 | bool HiddenValleyFragmentation::extractHVevent(Event& event) { |
---|
257 | |
---|
258 | // Copy Hidden-Valley particles to special event record. |
---|
259 | for (int i = 0; i < event.size(); ++i) { |
---|
260 | int idAbs = event[i].idAbs(); |
---|
261 | bool isHV = (idAbs > 4900000 && idAbs < 4900007) |
---|
262 | || (idAbs > 4900010 && idAbs < 4900017) |
---|
263 | || idAbs == 4900021 || idAbs == 4900101; |
---|
264 | if (isHV) { |
---|
265 | int iHV = hvEvent.append( event[i]); |
---|
266 | // Convert HV-gluons into normal ones so as to use normal machinery. |
---|
267 | if (event[i].id() == 4900021) hvEvent[iHV].id(21); |
---|
268 | // Second mother points back to position in complete event; |
---|
269 | // otherwise construct the HV history inside hvEvent. |
---|
270 | hvEvent[iHV].mothers( 0, i); |
---|
271 | hvEvent[iHV].daughters( 0, 0); |
---|
272 | int iMother = event[i].mother1(); |
---|
273 | for (int iHVM = 1; iHVM < hvEvent.size(); ++iHVM) |
---|
274 | if (hvEvent[iHVM].mother2() == iMother) { |
---|
275 | hvEvent[iHV].mother1( iHVM); |
---|
276 | if (hvEvent[iHVM].daughter1() == 0) hvEvent[iHVM].daughter1(iHV); |
---|
277 | else hvEvent[iHVM].daughter2(iHV); |
---|
278 | } |
---|
279 | } |
---|
280 | } |
---|
281 | |
---|
282 | // Done if no HV particles found. |
---|
283 | hvOldSize = hvEvent.size(); |
---|
284 | if (hvOldSize == 1) return false; |
---|
285 | |
---|
286 | // Initial colour - anticolour parton pair. |
---|
287 | int colBeg = hvEvent.nextColTag(); |
---|
288 | for (int iHV = 1; iHV < hvOldSize; ++iHV) |
---|
289 | if (hvEvent[iHV].mother1() == 0) { |
---|
290 | if (hvEvent[iHV].id() > 0) hvEvent[iHV].col( colBeg); |
---|
291 | else hvEvent[iHV].acol( colBeg); |
---|
292 | } |
---|
293 | |
---|
294 | // Then trace colours down to daughters; new colour if two daughters. |
---|
295 | for (int iHV = 1; iHV < hvOldSize; ++iHV) { |
---|
296 | int dau1 = hvEvent[iHV].daughter1(); |
---|
297 | int dau2 = hvEvent[iHV].daughter2(); |
---|
298 | if (dau1 > 0 && dau2 == 0) |
---|
299 | hvEvent[dau1].cols( hvEvent[iHV].col(), hvEvent[iHV].acol()); |
---|
300 | else if (dau2 > 0) { |
---|
301 | int colHV = hvEvent[iHV].col(); |
---|
302 | int acolHV = hvEvent[iHV].acol(); |
---|
303 | int colNew = hvEvent.nextColTag(); |
---|
304 | if (acolHV == 0) { |
---|
305 | hvEvent[dau1].cols( colNew, 0); |
---|
306 | hvEvent[dau2].cols( colHV, colNew); |
---|
307 | } else if (colHV == 0) { |
---|
308 | hvEvent[dau1].cols( 0, colNew); |
---|
309 | hvEvent[dau2].cols( colNew, acolHV); |
---|
310 | // Temporary: should seek recoiling dipole end!?? |
---|
311 | } else if (rndmPtr->flat() > 0.5) { |
---|
312 | hvEvent[dau1].cols( colHV, colNew); |
---|
313 | hvEvent[dau2].cols( colNew, acolHV); |
---|
314 | } else { |
---|
315 | hvEvent[dau1].cols( colNew, acolHV); |
---|
316 | hvEvent[dau2].cols( colHV, colNew); |
---|
317 | } |
---|
318 | } |
---|
319 | } |
---|
320 | |
---|
321 | // Pick up the colour end. |
---|
322 | int colNow = 0; |
---|
323 | for (int iHV = 1; iHV < hvOldSize; ++iHV) |
---|
324 | if (hvEvent[iHV].isFinal() && hvEvent[iHV].acol() == 0) { |
---|
325 | ihvParton.push_back( iHV); |
---|
326 | colNow = hvEvent[iHV].col(); |
---|
327 | } |
---|
328 | |
---|
329 | // Trace colour by colour until reached anticolour end. |
---|
330 | while (colNow > 0) { |
---|
331 | for (int iHV = 1; iHV < hvOldSize; ++iHV) |
---|
332 | if (hvEvent[iHV].isFinal() && hvEvent[iHV].acol() == colNow) { |
---|
333 | ihvParton.push_back( iHV); |
---|
334 | colNow = hvEvent[iHV].col(); |
---|
335 | break; |
---|
336 | } |
---|
337 | } |
---|
338 | |
---|
339 | // Done. |
---|
340 | return true; |
---|
341 | |
---|
342 | } |
---|
343 | |
---|
344 | //-------------------------------------------------------------------------- |
---|
345 | |
---|
346 | // Collapse of light system to one HV-meson, by the emission of HV-glueballs. |
---|
347 | |
---|
348 | bool HiddenValleyFragmentation::collapseToMeson() { |
---|
349 | |
---|
350 | // If too low mass then cannot do anything. Should not happen. |
---|
351 | if (mSys < 1.001 * mhvMeson) { |
---|
352 | infoPtr->errorMsg("Error in HiddenValleyFragmentation::collapseToMeson:" |
---|
353 | " too low mass to do anything"); |
---|
354 | return false; |
---|
355 | } |
---|
356 | |
---|
357 | // Choose mass of collective HV-glueball states flat between limits. |
---|
358 | double mhvGlue = (0.001 + 0.998 * rndmPtr->flat()) * (mSys - mhvMeson); |
---|
359 | |
---|
360 | // Find momentum in rest frame, with isotropic "decay" angles. |
---|
361 | double pAbs = 0.5 * sqrtpos( pow2(mSys*mSys - mhvMeson*mhvMeson |
---|
362 | - mhvGlue*mhvGlue) - pow2(2. * mhvMeson * mhvGlue) ) / mSys; |
---|
363 | double pz = (2 * rndmPtr->flat() - 1.) * pAbs; |
---|
364 | double pT = sqrtpos( pAbs*pAbs - pz*pz); |
---|
365 | double phi = 2. * M_PI * rndmPtr->flat(); |
---|
366 | double px = pT * cos(phi); |
---|
367 | double py = pT * sin(phi); |
---|
368 | |
---|
369 | // Construct four-vectors and boost them to event frame. |
---|
370 | Vec4 phvMeson( px, py, pz, sqrt(mhvMeson*mhvMeson + pAbs*pAbs) ); |
---|
371 | Vec4 phvGlue( -px, -py, -pz, sqrt(mhvGlue*mhvGlue + pAbs*pAbs) ); |
---|
372 | phvMeson.bst( hvColConfig[0].pSum ); |
---|
373 | phvGlue.bst( hvColConfig[0].pSum ); |
---|
374 | |
---|
375 | // Add produced particles to the event record. |
---|
376 | vector<int> iParton = hvColConfig[0].iParton; |
---|
377 | int iFirst = hvEvent.append( 4900111, 82, iParton.front(), |
---|
378 | iParton.back(), 0, 0, 0, 0, phvMeson, mhvMeson); |
---|
379 | int iLast = hvEvent.append( 4900991, 82, iParton.front(), |
---|
380 | iParton.back(), 0, 0, 0, 0, phvGlue, mhvGlue); |
---|
381 | |
---|
382 | // Mark original partons as hadronized and set their daughter range. |
---|
383 | for (int i = 0; i < int(iParton.size()); ++i) { |
---|
384 | hvEvent[ iParton[i] ].statusNeg(); |
---|
385 | hvEvent[ iParton[i] ].daughters(iFirst, iLast); |
---|
386 | } |
---|
387 | |
---|
388 | // Done. |
---|
389 | return true; |
---|
390 | |
---|
391 | } |
---|
392 | |
---|
393 | //-------------------------------------------------------------------------- |
---|
394 | |
---|
395 | // Insert HV-particles from hvEvent to event. |
---|
396 | |
---|
397 | bool HiddenValleyFragmentation::insertHVevent(Event& event) { |
---|
398 | |
---|
399 | // Offset for mother/daughter indices. |
---|
400 | hvNewSize = hvEvent.size(); |
---|
401 | int nOffset = event.size() - hvOldSize; |
---|
402 | |
---|
403 | // Copy back HV-particles. |
---|
404 | int iNew, iMot1, iMot2, iDau1, iDau2; |
---|
405 | for (int iHV = hvOldSize; iHV < hvNewSize; ++iHV) { |
---|
406 | iNew = event.append( hvEvent[iHV]); |
---|
407 | |
---|
408 | // Restore HV-gluon codes. Do not keep HV-colours, to avoid confusion. |
---|
409 | if (hvEvent[iHV].id() == 21) event[iNew].id(4900021); |
---|
410 | event[iNew].cols( 0, 0); |
---|
411 | |
---|
412 | // Begin history construction. |
---|
413 | iMot1 = hvEvent[iHV].mother1(); |
---|
414 | iMot2 = hvEvent[iHV].mother2(); |
---|
415 | iDau1 = hvEvent[iHV].daughter1(); |
---|
416 | iDau2 = hvEvent[iHV].daughter2(); |
---|
417 | // Special mother for partons copied from event, else simple offset. |
---|
418 | // Also set daughters of mothers in original record. |
---|
419 | if (iMot1 > 0 && iMot1 < hvOldSize) { |
---|
420 | iMot1 = hvEvent[iMot1].mother2(); |
---|
421 | event[iMot1].statusNeg(); |
---|
422 | event[iMot1].daughter1(iNew); |
---|
423 | } else if (iMot1 > 0) iMot1 += nOffset; |
---|
424 | if (iMot2 > 0 && iMot2 < hvOldSize) { |
---|
425 | iMot2 = hvEvent[iMot2].mother2(); |
---|
426 | event[iMot2].statusNeg(); |
---|
427 | if (event[iMot2].daughter1() == 0) event[iMot2].daughter1(iNew); |
---|
428 | else event[iMot2].daughter2(iNew); |
---|
429 | } else if (iMot2 > 0) iMot2 += nOffset; |
---|
430 | if (iDau1 > 0) iDau1 += nOffset; |
---|
431 | if (iDau2 > 0) iDau2 += nOffset; |
---|
432 | event[iNew].mothers( iMot1, iMot2); |
---|
433 | event[iNew].daughters( iDau1, iDau2); |
---|
434 | } |
---|
435 | |
---|
436 | // Done. |
---|
437 | return true; |
---|
438 | |
---|
439 | } |
---|
440 | |
---|
441 | //========================================================================== |
---|
442 | |
---|
443 | } // end namespace Pythia8 |
---|