]>
Commit | Line | Data |
---|---|---|
c6b60c38 | 1 | // HiddenValleyFragmentation.cc is a part of the PYTHIA event generator. |
2 | // Copyright (C) 2013 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 |