Line data Source code
1 : // HiddenValleyFragmentation.cc is a part of the PYTHIA event generator.
2 : // Copyright (C) 2015 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 "Pythia8/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 0 : rndmPtr = rndmPtrIn;
25 :
26 : // Read in data from Settings.
27 0 : nFlav = settings.mode("HiddenValley:nFlav");
28 0 : probVector = settings.parm("HiddenValley:probVector");
29 :
30 0 : }
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 0 : FlavContainer flavNew;
40 0 : flavNew.rank = flavOld.rank + 1;
41 :
42 : // Pick new HV-flavour at random; keep track of sign.
43 0 : flavNew.id = 4900100 + min( 1 + int(nFlav * rndmPtr->flat()), nFlav);
44 0 : if (flavOld.id > 0) flavNew.id = -flavNew.id;
45 :
46 : // Done.
47 0 : 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 0 : int idPos = max( flav1.id, flav2.id) - 4900000;
62 0 : int idNeg = -min( flav1.id, flav2.id) - 4900000;
63 0 : if (idPos < 20) idPos = 101;
64 0 : if (idNeg < 20) idNeg = 101;
65 :
66 : // Pick HV-meson code, spin either 0 or 1.
67 0 : if (idNeg == idPos) idMeson = 4900111;
68 0 : else if (idPos > idNeg) idMeson = 4900211;
69 : else idMeson = -4900211;
70 0 : if (rndmPtr->flat() < probVector) idMeson += ((idMeson > 0) ? 2 : -2);
71 :
72 : // Done.
73 0 : 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 0 : rndmPtr = rndmPtrIn;
90 :
91 : // Parameter of the pT width. No enhancement, since this is finetuning.
92 0 : double sigmamqv = settings.parm("HiddenValley:sigmamqv");
93 0 : double sigma = sigmamqv * particleData.m0( 4900101);
94 0 : sigmaQ = sigma / sqrt(2.);
95 0 : enhancedFraction = 0.;
96 0 : enhancedWidth = 0.;
97 :
98 : // Parameter for pT suppression in MiniStringFragmentation.
99 0 : sigma2Had = 2. * pow2( max( SIGMAMIN, sigma) );
100 :
101 0 : }
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 0 : rndmPtr = rndmPtrIn;
116 :
117 : // Paramaters of Lund/Bowler symmetric fragmentation function.
118 0 : aLund = settings.parm("HiddenValley:aLund");
119 0 : bmqv2 = settings.parm("HiddenValley:bmqv2");
120 0 : rFactqv = settings.parm("HiddenValley:rFactqv");
121 :
122 : // Use qv mass to set scale of bEff = b * m^2;
123 0 : mqv2 = pow2( particleData.m0( 4900101) );
124 0 : bLund = bmqv2 / mqv2;
125 :
126 : // Mass of qv meson used to set stop scale for fragmentation iteration.
127 0 : mhvMeson = particleData.m0( 4900111);
128 :
129 0 : }
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 0 : double bShape = bLund * mT2;
139 0 : double cShape = 1. + rFactqv * bmqv2;
140 0 : 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 0 : infoPtr = infoPtrIn;
157 0 : particleDataPtr = particleDataPtrIn;
158 0 : rndmPtr = rndmPtrIn;
159 :
160 : // Check whether Hidden Valley fragmentation switched on, and SU(N).
161 0 : doHVfrag = settings.flag("HiddenValley:fragment");
162 0 : if (settings.mode("HiddenValley:Ngauge") < 2) doHVfrag = false;
163 0 : if (!doHVfrag) return false;
164 :
165 : // Several copies of qv may be needed. Taken to have same mass.
166 0 : nFlav = settings.mode("HiddenValley:nFlav");
167 0 : if (nFlav > 1) {
168 0 : int spinType = particleDataPtr->spinType(4900101);
169 0 : double m0 = particleDataPtr->m0(4900101);
170 0 : for (int iFlav = 2; iFlav <= nFlav; ++iFlav)
171 0 : particleDataPtr->addParticle( 4900100 + iFlav, "qv", "qvbar",
172 : spinType, 0, 0, m0);
173 0 : }
174 :
175 : // Hidden Valley meson mass used to choose hadronization mode.
176 0 : mhvMeson = particleDataPtr->m0(4900111);
177 :
178 : // Initialize the hvEvent instance of an event record.
179 0 : hvEvent.init( "(Hidden Valley fragmentation)", particleDataPtr);
180 :
181 : // Create HVStringFlav instance for HV-flavour selection.
182 0 : hvFlavSelPtr = new HVStringFlav();
183 0 : hvFlavSelPtr->init( settings, rndmPtr);
184 :
185 : // Create HVStringPT instance for pT selection in HV fragmentation.
186 0 : hvPTSelPtr = new HVStringPT();
187 0 : hvPTSelPtr->init( settings, *particleDataPtr, rndmPtr);
188 :
189 : // Create HVStringZ instance for z selection in HV fragmentation.
190 0 : hvZSelPtr = new HVStringZ();
191 0 : hvZSelPtr->init( settings, *particleDataPtr, rndmPtr);
192 :
193 : // Initialize auxiliary administrative class.
194 0 : hvColConfig.init(infoPtr, settings, hvFlavSelPtr);
195 :
196 : // Initialize HV-string and HV-ministring fragmentation.
197 0 : hvStringFrag.init(infoPtr, settings, particleDataPtr, rndmPtr,
198 0 : hvFlavSelPtr, hvPTSelPtr, hvZSelPtr);
199 0 : hvMinistringFrag.init(infoPtr, settings, particleDataPtr, rndmPtr,
200 0 : hvFlavSelPtr, hvPTSelPtr, hvZSelPtr);
201 :
202 : // Done.
203 0 : return true;
204 :
205 0 : }
206 :
207 : //--------------------------------------------------------------------------
208 :
209 : // Perform the fragmentation.
210 :
211 : bool HiddenValleyFragmentation::fragment(Event& event) {
212 :
213 : // Reset containers for next event.
214 0 : hvEvent.reset();
215 0 : hvColConfig.clear();
216 0 : ihvParton.resize(0);
217 :
218 : // Extract HV-particles from event to hvEvent. Assign HV-colours.
219 : // Done if no HV-particles found.
220 0 : if (!extractHVevent(event)) return true;
221 :
222 : // Store found string system. Analyze its properties.
223 0 : 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 0 : hvColConfig.collect(0, hvEvent, false);
228 :
229 : // Mass used to decide how to fragment system.
230 0 : mSys = hvColConfig[0].mass;
231 :
232 : // HV-string fragmentation when enough mass to produce >= 3 HV-mesons.
233 0 : if (mSys > 3.5 * mhvMeson) {
234 0 : if (!hvStringFrag.fragment( 0, hvColConfig, hvEvent)) return false;
235 :
236 : // HV-ministring fragmentation when enough mass to produce 2 HV-mesons.
237 0 : } else if (mSys > 2.1 * mhvMeson) {
238 0 : if (!hvMinistringFrag.fragment( 0, hvColConfig, hvEvent, true))
239 0 : return false;
240 :
241 : // If only enough mass for one HV-meson assume HV-glueballs emitted.
242 0 : } else if (!collapseToMeson()) return false;
243 :
244 : // Insert HV particles from hvEvent to event.
245 0 : insertHVevent(event);
246 :
247 : // Done.
248 0 : return true;
249 :
250 0 : }
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 0 : for (int i = 0; i < event.size(); ++i) {
260 0 : int idAbs = event[i].idAbs();
261 0 : bool isHV = (idAbs > 4900000 && idAbs < 4900007)
262 0 : || (idAbs > 4900010 && idAbs < 4900017)
263 0 : || idAbs == 4900021 || idAbs == 4900101;
264 0 : if (isHV) {
265 0 : int iHV = hvEvent.append( event[i]);
266 : // Convert HV-gluons into normal ones so as to use normal machinery.
267 0 : 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 0 : hvEvent[iHV].mothers( 0, i);
271 0 : hvEvent[iHV].daughters( 0, 0);
272 0 : int iMother = event[i].mother1();
273 0 : for (int iHVM = 1; iHVM < hvEvent.size(); ++iHVM)
274 0 : if (hvEvent[iHVM].mother2() == iMother) {
275 0 : hvEvent[iHV].mother1( iHVM);
276 0 : if (hvEvent[iHVM].daughter1() == 0) hvEvent[iHVM].daughter1(iHV);
277 0 : else hvEvent[iHVM].daughter2(iHV);
278 : }
279 0 : }
280 : }
281 :
282 : // Done if no HV particles found.
283 0 : hvOldSize = hvEvent.size();
284 0 : if (hvOldSize == 1) return false;
285 :
286 : // Initial colour - anticolour parton pair.
287 0 : int colBeg = hvEvent.nextColTag();
288 0 : for (int iHV = 1; iHV < hvOldSize; ++iHV)
289 0 : if (hvEvent[iHV].mother1() == 0) {
290 0 : if (hvEvent[iHV].id() > 0) hvEvent[iHV].col( colBeg);
291 0 : else hvEvent[iHV].acol( colBeg);
292 : }
293 :
294 : // Then trace colours down to daughters; new colour if two daughters.
295 0 : for (int iHV = 1; iHV < hvOldSize; ++iHV) {
296 0 : int dau1 = hvEvent[iHV].daughter1();
297 0 : int dau2 = hvEvent[iHV].daughter2();
298 0 : if (dau1 > 0 && dau2 == 0)
299 0 : hvEvent[dau1].cols( hvEvent[iHV].col(), hvEvent[iHV].acol());
300 0 : else if (dau2 > 0) {
301 0 : int colHV = hvEvent[iHV].col();
302 0 : int acolHV = hvEvent[iHV].acol();
303 0 : int colNew = hvEvent.nextColTag();
304 0 : if (acolHV == 0) {
305 0 : hvEvent[dau1].cols( colNew, 0);
306 0 : hvEvent[dau2].cols( colHV, colNew);
307 0 : } else if (colHV == 0) {
308 0 : hvEvent[dau1].cols( 0, colNew);
309 0 : hvEvent[dau2].cols( colNew, acolHV);
310 : // Temporary: should seek recoiling dipole end!??
311 0 : } else if (rndmPtr->flat() > 0.5) {
312 0 : hvEvent[dau1].cols( colHV, colNew);
313 0 : hvEvent[dau2].cols( colNew, acolHV);
314 0 : } else {
315 0 : hvEvent[dau1].cols( colNew, acolHV);
316 0 : hvEvent[dau2].cols( colHV, colNew);
317 : }
318 0 : }
319 : }
320 :
321 : // Pick up the colour end.
322 : int colNow = 0;
323 0 : for (int iHV = 1; iHV < hvOldSize; ++iHV)
324 0 : if (hvEvent[iHV].isFinal() && hvEvent[iHV].acol() == 0) {
325 0 : ihvParton.push_back( iHV);
326 0 : colNow = hvEvent[iHV].col();
327 0 : }
328 :
329 : // Trace colour by colour until reached anticolour end.
330 0 : while (colNow > 0) {
331 0 : for (int iHV = 1; iHV < hvOldSize; ++iHV)
332 0 : if (hvEvent[iHV].isFinal() && hvEvent[iHV].acol() == colNow) {
333 0 : ihvParton.push_back( iHV);
334 0 : colNow = hvEvent[iHV].col();
335 0 : break;
336 : }
337 : }
338 :
339 : // Done.
340 : return true;
341 :
342 0 : }
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 0 : if (mSys < 1.001 * mhvMeson) {
352 0 : infoPtr->errorMsg("Error in HiddenValleyFragmentation::collapseToMeson:"
353 : " too low mass to do anything");
354 0 : return false;
355 : }
356 :
357 : // Choose mass of collective HV-glueball states flat between limits.
358 0 : double mhvGlue = (0.001 + 0.998 * rndmPtr->flat()) * (mSys - mhvMeson);
359 :
360 : // Find momentum in rest frame, with isotropic "decay" angles.
361 0 : double pAbs = 0.5 * sqrtpos( pow2(mSys*mSys - mhvMeson*mhvMeson
362 0 : - mhvGlue*mhvGlue) - pow2(2. * mhvMeson * mhvGlue) ) / mSys;
363 0 : double pz = (2 * rndmPtr->flat() - 1.) * pAbs;
364 0 : double pT = sqrtpos( pAbs*pAbs - pz*pz);
365 0 : double phi = 2. * M_PI * rndmPtr->flat();
366 0 : double px = pT * cos(phi);
367 0 : double py = pT * sin(phi);
368 :
369 : // Construct four-vectors and boost them to event frame.
370 0 : Vec4 phvMeson( px, py, pz, sqrt(mhvMeson*mhvMeson + pAbs*pAbs) );
371 0 : Vec4 phvGlue( -px, -py, -pz, sqrt(mhvGlue*mhvGlue + pAbs*pAbs) );
372 0 : phvMeson.bst( hvColConfig[0].pSum );
373 0 : phvGlue.bst( hvColConfig[0].pSum );
374 :
375 : // Add produced particles to the event record.
376 0 : vector<int> iParton = hvColConfig[0].iParton;
377 0 : int iFirst = hvEvent.append( 4900111, 82, iParton.front(),
378 0 : iParton.back(), 0, 0, 0, 0, phvMeson, mhvMeson);
379 0 : int iLast = hvEvent.append( 4900991, 82, iParton.front(),
380 0 : iParton.back(), 0, 0, 0, 0, phvGlue, mhvGlue);
381 :
382 : // Mark original partons as hadronized and set their daughter range.
383 0 : for (int i = 0; i < int(iParton.size()); ++i) {
384 0 : hvEvent[ iParton[i] ].statusNeg();
385 0 : hvEvent[ iParton[i] ].daughters(iFirst, iLast);
386 : }
387 :
388 : // Done.
389 : return true;
390 :
391 0 : }
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 0 : hvNewSize = hvEvent.size();
401 0 : int nOffset = event.size() - hvOldSize;
402 :
403 : // Copy back HV-particles.
404 : int iNew, iMot1, iMot2, iDau1, iDau2;
405 0 : for (int iHV = hvOldSize; iHV < hvNewSize; ++iHV) {
406 0 : iNew = event.append( hvEvent[iHV]);
407 :
408 : // Restore HV-gluon codes. Do not keep HV-colours, to avoid confusion.
409 0 : if (hvEvent[iHV].id() == 21) event[iNew].id(4900021);
410 0 : event[iNew].cols( 0, 0);
411 :
412 : // Begin history construction.
413 0 : iMot1 = hvEvent[iHV].mother1();
414 0 : iMot2 = hvEvent[iHV].mother2();
415 0 : iDau1 = hvEvent[iHV].daughter1();
416 0 : 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 0 : if (iMot1 > 0 && iMot1 < hvOldSize) {
420 0 : iMot1 = hvEvent[iMot1].mother2();
421 0 : event[iMot1].statusNeg();
422 0 : event[iMot1].daughter1(iNew);
423 0 : } else if (iMot1 > 0) iMot1 += nOffset;
424 0 : if (iMot2 > 0 && iMot2 < hvOldSize) {
425 0 : iMot2 = hvEvent[iMot2].mother2();
426 0 : event[iMot2].statusNeg();
427 0 : if (event[iMot2].daughter1() == 0) event[iMot2].daughter1(iNew);
428 0 : else event[iMot2].daughter2(iNew);
429 0 : } else if (iMot2 > 0) iMot2 += nOffset;
430 0 : if (iDau1 > 0) iDau1 += nOffset;
431 0 : if (iDau2 > 0) iDau2 += nOffset;
432 0 : event[iNew].mothers( iMot1, iMot2);
433 0 : event[iNew].daughters( iDau1, iDau2);
434 : }
435 :
436 : // Done.
437 0 : return true;
438 :
439 : }
440 :
441 : //==========================================================================
442 :
443 : } // end namespace Pythia8
|
