9 #include "Pythia8/HiddenValleyFragmentation.h"
21 void HVStringFlav::init(Settings& settings, ParticleData* particleDataPtrIn,
22 Rndm* rndmPtrIn, Info* infoPtrIn) {
25 particleDataPtr = particleDataPtrIn;
30 nFlav = settings.mode(
"HiddenValley:nFlav");
31 probVector = settings.parm(
"HiddenValley:probVector");
35 mT2suppression =
false;
43 FlavContainer HVStringFlav::pick(FlavContainer& flavOld,
double,
double) {
46 FlavContainer flavNew;
47 flavNew.rank = flavOld.rank + 1;
50 flavNew.id = 4900100 + min( 1 +
int(nFlav * rndmPtr->flat()), nFlav);
51 if (flavOld.id > 0) flavNew.id = -flavNew.id;
63 int HVStringFlav::combine(FlavContainer& flav1, FlavContainer& flav2) {
68 int idPos = max( flav1.id, flav2.id) - 4900000;
69 int idNeg = -min( flav1.id, flav2.id) - 4900000;
70 if (idPos < 20) idPos = 101;
71 if (idNeg < 20) idNeg = 101;
74 if (idNeg == idPos) idMeson = 4900111;
75 else if (idPos > idNeg) idMeson = 4900211;
76 else idMeson = -4900211;
77 if (rndmPtr->flat() < probVector) idMeson += ((idMeson > 0) ? 2 : -2);
92 void HVStringPT::init(Settings& settings, ParticleData* particleDataPtrIn,
93 Rndm* rndmPtrIn, Info* infoPtrIn) {
96 particleDataPtr = particleDataPtrIn;
101 double sigmamqv = settings.parm(
"HiddenValley:sigmamqv");
102 double sigma = sigmamqv * particleDataPtrIn->m0( 4900101);
103 sigmaQ = sigma / sqrt(2.);
104 enhancedFraction = 0.;
108 sigma2Had = 2. * pow2( max( SIGMAMIN, sigma) );
109 thermalModel =
false;
111 closePacking =
false;
123 void HVStringZ::init(Settings& settings, ParticleData& particleData,
124 Rndm* rndmPtrIn, Info* infoPtrIn) {
131 aLund = settings.parm(
"HiddenValley:aLund");
132 bmqv2 = settings.parm(
"HiddenValley:bmqv2");
133 rFactqv = settings.parm(
"HiddenValley:rFactqv");
136 mqv2 = pow2( particleData.m0( 4900101) );
137 bLund = bmqv2 / mqv2;
140 mhvMeson = particleData.m0( 4900111);
148 double HVStringZ::zFrag(
int ,
int ,
double mT2) {
151 double bShape = bLund * mT2;
152 double cShape = 1. + rFactqv * bmqv2;
153 return zLund( aLund, bShape, cShape);
165 bool HiddenValleyFragmentation::init(Info* infoPtrIn, Settings& settings,
166 ParticleData* particleDataPtrIn, Rndm* rndmPtrIn) {
170 particleDataPtr = particleDataPtrIn;
174 doHVfrag = settings.flag(
"HiddenValley:fragment");
175 if (settings.mode(
"HiddenValley:Ngauge") < 2) doHVfrag =
false;
176 if (!doHVfrag)
return false;
179 nFlav = settings.mode(
"HiddenValley:nFlav");
181 int spinType = particleDataPtr->spinType(4900101);
182 double m0 = particleDataPtr->m0(4900101);
183 for (
int iFlav = 2; iFlav <= nFlav; ++iFlav)
184 particleDataPtr->addParticle( 4900100 + iFlav,
"qv",
"qvbar",
189 mhvMeson = particleDataPtr->m0(4900111);
192 hvEvent.init(
"(Hidden Valley fragmentation)", particleDataPtr);
195 hvFlavSelPtr =
new HVStringFlav();
196 hvFlavSelPtr->init( settings, particleDataPtr, rndmPtr, infoPtr);
199 hvPTSelPtr =
new HVStringPT();
200 hvPTSelPtr->init( settings, particleDataPtr, rndmPtr, infoPtr);
203 hvZSelPtr =
new HVStringZ();
204 hvZSelPtr->init( settings, *particleDataPtr, rndmPtr, infoPtr);
207 hvColConfig.init(infoPtr, settings, hvFlavSelPtr);
210 hvStringFrag.init(infoPtr, settings, particleDataPtr, rndmPtr,
211 hvFlavSelPtr, hvPTSelPtr, hvZSelPtr);
212 hvMinistringFrag.init(infoPtr, settings, particleDataPtr, rndmPtr,
213 hvFlavSelPtr, hvPTSelPtr, hvZSelPtr);
224 bool HiddenValleyFragmentation::fragment(
Event& event) {
233 if (!extractHVevent(event))
return true;
236 if (!hvColConfig.insert(ihvParton, hvEvent))
return false;
240 hvColConfig.collect(0, hvEvent,
false);
243 mSys = hvColConfig[0].mass;
246 if (mSys > 3.5 * mhvMeson) {
247 if (!hvStringFrag.fragment( 0, hvColConfig, hvEvent))
return false;
250 }
else if (mSys > 2.1 * mhvMeson) {
251 if (!hvMinistringFrag.fragment( 0, hvColConfig, hvEvent,
true))
255 }
else if (!collapseToMeson())
return false;
258 insertHVevent(event);
269 bool HiddenValleyFragmentation::extractHVevent(
Event& event) {
272 for (
int i = 0; i <
event.size(); ++i) {
273 int idAbs =
event[i].idAbs();
274 bool isHV = (idAbs > 4900000 && idAbs < 4900007)
275 || (idAbs > 4900010 && idAbs < 4900017)
277 || (idAbs > 4900100 && idAbs < 4900109);
279 int iHV = hvEvent.append( event[i]);
281 if (event[i].
id() == 4900021) hvEvent[iHV].id(21);
284 hvEvent[iHV].mothers( 0, i);
285 hvEvent[iHV].daughters( 0, 0);
286 int iMother =
event[i].mother1();
287 for (
int iHVM = 1; iHVM < hvEvent.size(); ++iHVM)
288 if (hvEvent[iHVM].mother2() == iMother) {
289 hvEvent[iHV].mother1( iHVM);
290 if (hvEvent[iHVM].daughter1() == 0) hvEvent[iHVM].daughter1(iHV);
291 else hvEvent[iHVM].daughter2(iHV);
297 hvOldSize = hvEvent.size();
298 if (hvOldSize == 1)
return false;
301 int colBeg = hvEvent.nextColTag();
302 for (
int iHV = 1; iHV < hvOldSize; ++iHV)
303 if (hvEvent[iHV].mother1() == 0) {
304 if (hvEvent[iHV].
id() > 0) hvEvent[iHV].col( colBeg);
305 else hvEvent[iHV].acol( colBeg);
309 for (
int iHV = 1; iHV < hvOldSize; ++iHV) {
310 int dau1 = hvEvent[iHV].daughter1();
311 int dau2 = hvEvent[iHV].daughter2();
312 if (dau1 > 0 && dau2 == 0)
313 hvEvent[dau1].cols( hvEvent[iHV].col(), hvEvent[iHV].acol());
315 int colHV = hvEvent[iHV].col();
316 int acolHV = hvEvent[iHV].acol();
317 int colNew = hvEvent.nextColTag();
319 hvEvent[dau1].cols( colNew, 0);
320 hvEvent[dau2].cols( colHV, colNew);
321 }
else if (colHV == 0) {
322 hvEvent[dau1].cols( 0, colNew);
323 hvEvent[dau2].cols( colNew, acolHV);
325 }
else if (rndmPtr->flat() > 0.5) {
326 hvEvent[dau1].cols( colHV, colNew);
327 hvEvent[dau2].cols( colNew, acolHV);
329 hvEvent[dau1].cols( colNew, acolHV);
330 hvEvent[dau2].cols( colHV, colNew);
337 for (
int iHV = 1; iHV < hvOldSize; ++iHV)
338 if (hvEvent[iHV].isFinal() && hvEvent[iHV].acol() == 0) {
339 ihvParton.push_back( iHV);
340 colNow = hvEvent[iHV].col();
345 for (
int iHV = 1; iHV < hvOldSize; ++iHV)
346 if (hvEvent[iHV].isFinal() && hvEvent[iHV].acol() == colNow) {
347 ihvParton.push_back( iHV);
348 colNow = hvEvent[iHV].col();
362 bool HiddenValleyFragmentation::collapseToMeson() {
365 if (mSys < 1.001 * mhvMeson) {
366 infoPtr->errorMsg(
"Error in HiddenValleyFragmentation::collapseToMeson:"
367 " too low mass to do anything");
372 double mhvGlue = (0.001 + 0.998 * rndmPtr->flat()) * (mSys - mhvMeson);
375 double pAbs = 0.5 * sqrtpos( pow2(mSys*mSys - mhvMeson*mhvMeson
376 - mhvGlue*mhvGlue) - pow2(2. * mhvMeson * mhvGlue) ) / mSys;
377 double pz = (2 * rndmPtr->flat() - 1.) * pAbs;
378 double pT = sqrtpos( pAbs*pAbs - pz*pz);
379 double phi = 2. * M_PI * rndmPtr->flat();
380 double px = pT * cos(phi);
381 double py = pT * sin(phi);
384 Vec4 phvMeson( px, py, pz, sqrt(mhvMeson*mhvMeson + pAbs*pAbs) );
385 Vec4 phvGlue( -px, -py, -pz, sqrt(mhvGlue*mhvGlue + pAbs*pAbs) );
386 phvMeson.bst( hvColConfig[0].pSum );
387 phvGlue.bst( hvColConfig[0].pSum );
390 vector<int> iParton = hvColConfig[0].iParton;
391 int iFirst = hvEvent.append( 4900111, 82, iParton.front(),
392 iParton.back(), 0, 0, 0, 0, phvMeson, mhvMeson);
393 int iLast = hvEvent.append( 4900991, 82, iParton.front(),
394 iParton.back(), 0, 0, 0, 0, phvGlue, mhvGlue);
397 for (
int i = 0; i < int(iParton.size()); ++i) {
398 hvEvent[ iParton[i] ].statusNeg();
399 hvEvent[ iParton[i] ].daughters(iFirst, iLast);
411 bool HiddenValleyFragmentation::insertHVevent(
Event& event) {
414 hvNewSize = hvEvent.size();
415 int nOffset =
event.size() - hvOldSize;
418 int iNew, iMot1, iMot2, iDau1, iDau2;
419 for (
int iHV = hvOldSize; iHV < hvNewSize; ++iHV) {
420 iNew =
event.append( hvEvent[iHV]);
423 if (hvEvent[iHV].
id() == 21)
event[iNew].id(4900021);
424 event[iNew].cols( 0, 0);
427 iMot1 = hvEvent[iHV].mother1();
428 iMot2 = hvEvent[iHV].mother2();
429 iDau1 = hvEvent[iHV].daughter1();
430 iDau2 = hvEvent[iHV].daughter2();
433 if (iMot1 > 0 && iMot1 < hvOldSize) {
434 iMot1 = hvEvent[iMot1].mother2();
435 event[iMot1].statusNeg();
436 event[iMot1].daughter1(iNew);
437 }
else if (iMot1 > 0) iMot1 += nOffset;
438 if (iMot2 > 0 && iMot2 < hvOldSize) {
439 iMot2 = hvEvent[iMot2].mother2();
440 event[iMot2].statusNeg();
441 if (event[iMot2].daughter1() == 0)
event[iMot2].daughter1(iNew);
442 else event[iMot2].daughter2(iNew);
443 }
else if (iMot2 > 0) iMot2 += nOffset;
444 if (iDau1 > 0) iDau1 += nOffset;
445 if (iDau2 > 0) iDau2 += nOffset;
446 event[iNew].mothers( iMot1, iMot2);
447 event[iNew].daughters( iDau1, iDau2);