9 #include "Pythia8/ResonanceDecays.h"
24 const int ResonanceDecays::NTRYCHANNEL = 10;
27 const int ResonanceDecays::NTRYMASSES = 10000;
30 const double ResonanceDecays::MSAFETY = 0.1;
33 const double ResonanceDecays::WIDTHCUT = 5.;
36 const double ResonanceDecays::TINY = 1e-10;
39 const double ResonanceDecays::TINYBWRANGE = 1e-8;
43 const double ResonanceDecays::WTCORRECTION[11] = { 1., 1., 1.,
44 2., 5., 15., 60., 250., 1250., 7000., 50000. };
48 bool ResonanceDecays::next(
Event& process,
int iDecNow) {
54 Particle& decayer = process[iDec];
55 if (decayer.isFinal() && decayer.canDecay() && decayer.mayDecay()
56 && decayer.isResonance() ) {
63 idProd.push_back( id0 );
64 mProd.push_back( m0 );
67 int idIn = process[decayer.mother1()].id();
70 if (!decayer.particleDataEntry().preparePick(id0, m0, idIn)) {
72 osWarn <<
"for id = " << id0;
73 infoPtr->errorMsg(
"Error in ResonanceDecays::next:"
74 " no open decay channel", osWarn.str());
79 bool foundChannel =
false;
80 for (
int iTryChannel = 0; iTryChannel < NTRYCHANNEL; ++iTryChannel) {
83 DecayChannel& channel = decayer.particleDataEntry().pickChannel();
84 mult = channel.multiplicity();
89 for (
int i = 1; i <= mult; ++i) {
90 idNow = channel.product(i - 1);
91 if (id0 < 0 && particleDataPtr->hasAnti(idNow)) idNow = -idNow;
92 idProd.push_back( idNow);
97 if (!pickMasses())
continue;
104 ostringstream osWarn;
105 osWarn <<
"for id = " << id0;
106 infoPtr->errorMsg(
"Error in ResonanceDecays::next:"
107 " failed to find workable decay channel", osWarn.str());
112 if (!pickColours(iDec, process))
return false;
116 pProd.push_back( decayer.p() );
117 if (!pickKinematics())
return false;
120 int iFirst = process.size();
121 for (
int i = 1; i <= mult; ++i) {
122 process.append( idProd[i], 23, iDec, 0, 0, 0, cols[i], acols[i],
123 pProd[i], mProd[i], m0);
125 int iLast = process.size() - 1;
128 if (process[iDec].hasVertex() || process[iDec].tau() > 0.) {
129 Vec4 vDec = process[iDec].vDec();
130 for (
int i = iFirst; i <= iLast; ++i) process[i].vProd( vDec );
134 for (
int i = iFirst; i <= iLast; ++i)
135 process[i].tau( process[i].tau0() * rndmPtr->exp() );
139 decayer.daughters(iFirst, iLast);
143 }
while (iDecNow == 0 && ++iDec < process.size());
154 bool ResonanceDecays::pickMasses() {
158 vector<double> m0BW, mMinBW, mMaxBW, widthBW;
159 double mMother = mProd[0];
160 double m2Mother = mMother * mMother;
161 useBW.push_back(
false );
162 m0BW.push_back( mMother );
163 mMinBW.push_back( mMother );
164 mMaxBW.push_back( mMother );
165 widthBW.push_back( 0. );
169 double m0Now, mMinNow, mMaxNow, widthNow;
170 for (
int i = 1; i <= mult; ++i) {
171 useBWNow = particleDataPtr->useBreitWigner( idProd[i] );
172 m0Now = particleDataPtr->m0( idProd[i] );
173 mMinNow = particleDataPtr->m0Min( idProd[i] );
174 mMaxNow = particleDataPtr->m0Max( idProd[i] );
175 if (useBWNow && mMaxNow < mMinNow) mMaxNow = mMother;
176 widthNow = particleDataPtr->mWidth( idProd[i] );
177 useBW.push_back( useBWNow );
178 m0BW.push_back( m0Now );
179 mMinBW.push_back( mMinNow );
180 mMaxBW.push_back( mMaxNow );
181 widthBW.push_back( widthNow );
182 mProd.push_back( m0Now );
189 for (
int i = 1; i <= mult; ++i) {
191 mSum += max( m0BW[i], mMinBW[i]);
192 mSumMin += mMinBW[i];
196 if (mSumMin + MSAFETY > mMother)
return false;
199 if (mSum + 0.5 * MSAFETY < mMother && nBW == 0)
return true;
202 if (mSum + MSAFETY < mMother) {
203 double wtMax = 2. * sqrtpos(1. - mSum*mSum / m2Mother);
205 for (
int iTryMasses = 0; iTryMasses <= NTRYMASSES; ++ iTryMasses) {
206 if (iTryMasses == NTRYMASSES)
return false;
208 for (
int i = 1; i <= mult; ++i) {
209 if (useBW[i]) mProd[i] = particleDataPtr->mSel( idProd[i] );
212 wt = (mSum + 0.5 * MSAFETY < mMother)
213 ? sqrtpos(1. - mSum*mSum / m2Mother) : 0.;
214 if (wt > rndmPtr->flat() * wtMax)
break;
224 for (
int i = 1; i <= mult; ++i) {
225 if (useBW[i]) iBW.push_back(i);
226 else mSum0 += mProd[i];
228 for (
int i = 1; i < nBW; ++i) {
229 for (
int j = i - 1; j >= 0; --j) {
230 if (widthBW[iBW[j+1]] > widthBW[iBW[j]]) swap (iBW[j+1], iBW[j]);
237 int iMin = (nBW == 1) ? 0 : 2;
238 for (
int i = nBW - 1; i >= iMin; --i) {
242 double mMax = mMother - mSum0 - MSAFETY;
243 if (nBW != 1)
for (
int j = 0; j < i; ++j) mMax -= mMinBW[iBW[j]];
244 mMax = min( mMaxBW[iBWi], mMax );
245 double mMin = min( mMinBW[iBWi], mMax - MSAFETY);
246 if (mMin < 0.)
return false;
249 double m2Nom = pow2( m0BW[iBWi] );
250 double m2Max = mMax * mMax;
251 double m2Min = mMin * mMin;
252 double mmWid = m0BW[iBWi] * widthBW[iBWi];
253 double atanMin = atan( (m2Min - m2Nom) / mmWid );
254 double atanMax = atan( (m2Max - m2Nom) / mmWid );
255 double atanDif = atanMax - atanMin;
258 if (atanDif < TINYBWRANGE)
return false;
261 double mr1 = mSum0*mSum0 / m2Mother;
262 double mr2 = m2Min / m2Mother;
263 double wtMax = sqrtpos( pow2(1. - mr1 - mr2) - 4. * mr1 * mr2 );
265 for (
int iTryMasses = 0; iTryMasses <= NTRYMASSES; ++ iTryMasses) {
266 if (iTryMasses == NTRYMASSES)
return false;
267 m2Now = m2Nom + mmWid * tan(atanMin + rndmPtr->flat() * atanDif);
268 mr2 = m2Now / m2Mother;
269 wt = sqrtpos( pow2(1. - mr1 - mr2) - 4. * mr1 * mr2 );
270 if (wt > rndmPtr->flat() * wtMax)
break;
274 mProd[iBWi] = sqrt(m2Now);
275 mSum0 += mProd[iBWi];
277 if (nBW == 1)
return true;
283 int idMother = abs(idProd[0]);
284 int idDau1 = abs(idProd[iBW1]);
285 int idDau2 = abs(idProd[iBW2]);
289 if ( (idMother == 25 || idMother == 35) && idDau1 < 19
290 && idDau2 == idDau1 ) psMode = 3;
291 if ( (idMother == 25 || idMother == 35 )
292 && (idDau1 == 23 || idDau1 == 24) && idDau2 == idDau1 ) psMode = 5;
294 && (idDau1 == 23 || idDau1 == 24) && idDau2 == idDau1 ) psMode = 6;
297 double mRem = mMother - mSum0 - MSAFETY;
298 double mMax1 = min( mMaxBW[iBW1], mRem - mMinBW[iBW2] );
299 double mMin1 = min( mMinBW[iBW1], mMax1 - MSAFETY);
300 double mMax2 = min( mMaxBW[iBW2], mRem - mMinBW[iBW1] );
301 double mMin2 = min( mMinBW[iBW2], mMax2 - MSAFETY);
304 if (mMin1 + mMin2 > mRem)
return false;
305 double mMid = 0.5 * mRem;
306 bool hasMid1 = (mMin1 < mMid);
307 bool hasMid2 = (mMin2 < mMid);
308 if (!hasMid1 && !hasMid2)
return false;
311 double m2Nom1 = pow2( m0BW[iBW1] );
312 double m2Max1 = mMax1 * mMax1;
313 double m2Min1 = mMin1 * mMin1;
314 double m2Mid1 = min( mMid * mMid, m2Max1);
315 double mmWid1 = m0BW[iBW1] * widthBW[iBW1];
316 double atanMin1 = atan( (m2Min1 - m2Nom1) / mmWid1 );
317 double atanMax1 = atan( (m2Max1 - m2Nom1) / mmWid1 );
318 double atanMid1 = (hasMid1) ? atan( (m2Mid1 - m2Nom1) / mmWid1 ) : 0.;
319 double m2Nom2 = pow2( m0BW[iBW2] );
320 double m2Max2 = mMax2 * mMax2;
321 double m2Min2 = mMin2 * mMin2;
322 double m2Mid2 = min( mMid * mMid, m2Max2);
323 double mmWid2 = m0BW[iBW2] * widthBW[iBW2];
324 double atanMin2 = atan( (m2Min2 - m2Nom2) / mmWid2 );
325 double atanMax2 = atan( (m2Max2 - m2Nom2) / mmWid2 );
326 double atanMid2 = (hasMid2) ? atan( (m2Mid2 - m2Nom2) / mmWid2 ) : 0.;
329 double probLow1 = (hasMid1) ? 1. : 0.;
330 if (hasMid1 && hasMid2) {
331 double intLow1 = (atanMid1 - atanMin1) * (atanMax2 - atanMin2);
332 double intLow2 = (atanMax1 - atanMin1) * (atanMid2 - atanMin2);
333 probLow1 = intLow1 / (intLow1 + intLow2);
337 double m2Rem = mRem * mRem;
338 double mr1 = m2Min1 / m2Rem;
339 double mr2 = m2Min2 / m2Rem;
340 double psMax = sqrtpos( pow2(1. - mr1 - mr2) - 4. * mr1 * mr2 );
342 if (psMode == 1) wtMax = psMax;
343 else if (psMode == 2) wtMax = psMax * psMax;
344 else if (psMode == 3) wtMax = pow3(psMax);
345 else if (psMode == 5) wtMax = psMax
346 * (pow2(1. - mr1 - mr2) + 8. * mr1 * mr2);
347 else if (psMode == 6) wtMax = pow3(psMax);
350 double atanDif1, atanDif2, m2Now1, m2Now2, mNow1, mNow2, ps, wt;
351 for (
int iTryMasses = 0; iTryMasses <= NTRYMASSES; ++ iTryMasses) {
352 if (iTryMasses == NTRYMASSES)
return false;
355 bool pickLow1 =
false;
356 if (rndmPtr->flat() < probLow1) {
357 atanDif1 = atanMid1 - atanMin1;
358 atanDif2 = atanMax2 - atanMin2;
361 atanDif1 = atanMax1 - atanMin1;
362 atanDif2 = atanMid2 - atanMin2;
364 m2Now1 = m2Nom1 + mmWid1 * tan(atanMin1 + rndmPtr->flat() * atanDif1);
365 m2Now2 = m2Nom2 + mmWid2 * tan(atanMin2 + rndmPtr->flat() * atanDif2);
366 mNow1 = sqrt(m2Now1);
367 mNow2 = sqrt(m2Now2);
370 bool rejectRegion = (pickLow1) ? (mNow1 > mNow2) : (mNow2 > mNow1);
371 if (rejectRegion)
continue;
374 mr1 = m2Now1 / m2Rem;
375 mr2 = m2Now2 / m2Rem;
377 if (mNow1 + mNow2 + MSAFETY < mMother) {
378 ps = sqrtpos( pow2(1. - mr1 - mr2) - 4. * mr1 * mr2 );
380 if (psMode == 1) wt = ps;
381 else if (psMode == 2) wt = ps * ps;
382 else if (psMode == 3) wt = pow3(ps);
383 else if (psMode == 5) wt = ps
384 * (pow2(1. - mr1 - mr2) + 8. * mr1 * mr2);
385 else if (psMode == 6) wt = pow3(ps)*mr1*mr2;
387 if (wt > rndmPtr->flat() * wtMax)
break;
401 bool ResonanceDecays::pickColours(
int iDec,
Event& process) {
406 vector<int> iTriplet, iAtriplet, iOctet, iDipCol, iDipAcol;
409 int col0 = process[iDec].col();
410 int acol0 = process[iDec].acol();
411 int colType0 = process[iDec].colType();
412 cols.push_back( col0);
413 acols.push_back(acol0);
417 for (
int i = 1; i <= mult; ++i) {
422 colTypeNow = particleDataPtr->colType( idProd[i] );
423 if (colTypeNow == 0);
424 else if (colTypeNow == 1) iTriplet.push_back(i);
425 else if (colTypeNow == -1) iAtriplet.push_back(i);
426 else if (colTypeNow == 2) iOctet.push_back(i);
428 else if (colTypeNow == 3) {
429 iTriplet.push_back(i);
430 iTriplet.push_back(i);
431 }
else if (colTypeNow == -3) {
432 iAtriplet.push_back(i);
433 iAtriplet.push_back(i);
435 infoPtr->errorMsg(
"Error in ResonanceDecays::pickColours:"
436 " unknown colour type encountered");
442 int nCol = iTriplet.size();
443 if (colType0 == 1 || colType0 == 2) nCol -= 1;
444 else if (colType0 == 3) nCol -= 2;
445 int nAcol = iAtriplet.size();
446 if (colType0 == -1 || colType0 == 2) nAcol -= 1;
447 else if (colType0 == -3) nAcol -= 2;
453 if (nCol - nAcol == 3) {
454 int kindJun = (colType0 == 0 || colType0 == 1 || colType0 == 3) ? 1 : 3;
458 colJun[0] = (kindJun == 1) ? process.nextColTag() : acol0;
459 colJun[1] = process.nextColTag();
460 colJun[2] = process.nextColTag();
461 process.appendJunction( kindJun, colJun[0], colJun[1], colJun[2]);
464 for (
int leg = 0; leg < 3; ++leg) {
465 if (leg == 0 && kindJun != 1) acol0 = 0;
469 int pickT = (iTriplet.size() == 1) ? 0
470 :
int( TINY + rndmPtr->flat() * (iTriplet.size() - TINY) );
471 int iPickT = iTriplet[pickT];
472 cols[iPickT] = colJun[leg];
475 iTriplet[pickT] = iTriplet.back();
477 iDipCol.push_back(iPickT);
478 iDipAcol.push_back(0);
490 if (nAcol - nCol == 3) {
491 int kindJun = (colType0 == 0 || colType0 == -1 || colType0 == -3) ? 2 : 4;
495 acolJun[0] = (kindJun == 2) ? process.nextColTag() : col0;
496 acolJun[1] = process.nextColTag();
497 acolJun[2] = process.nextColTag();
498 process.appendJunction( kindJun, acolJun[0], acolJun[1], acolJun[2]);
501 for (
int leg = 0; leg < 3; ++leg) {
502 if (leg == 0 && kindJun != 2) col0 = 0;
506 int pickA = (iAtriplet.size() == 1) ? 0
507 :
int( TINY + rndmPtr->flat() * (iAtriplet.size() - TINY) );
508 int iPickA = iAtriplet[pickA];
509 acols[iPickA] = acolJun[leg];
512 iAtriplet[pickA] = iAtriplet.back();
513 iAtriplet.pop_back();
514 iDipCol.push_back(0);
515 iDipAcol.push_back(iPickA);
525 infoPtr->errorMsg(
"Error in ResonanceDecays::pickColours:"
526 " inconsistent colour tags");
531 if (col0 > 0 && iTriplet.size() > 0) {
532 int pickT = (iTriplet.size() == 1) ? 0
533 :
int( TINY + rndmPtr->flat() * (iTriplet.size() - TINY) );
534 int iPickT = iTriplet[pickT];
539 iTriplet[pickT] = iTriplet.back();
541 iDipCol.push_back(iPickT);
542 iDipAcol.push_back(0);
546 if (acol0 > 0 && iAtriplet.size() > 0) {
547 int pickA = (iAtriplet.size() == 1) ? 0
548 :
int( TINY + rndmPtr->flat() * (iAtriplet.size() - TINY) );
549 int iPickA = iAtriplet[pickA];
550 acols[iPickA] = acol0;
554 iAtriplet[pickA] = iAtriplet.back();
555 iAtriplet.pop_back();
556 iDipCol.push_back(0);
557 iDipAcol.push_back(iPickA);
561 if (acol0 < 0 && iTriplet.size() > 0) {
562 int pickT = (iTriplet.size() == 1) ? 0
563 :
int( TINY + rndmPtr->flat() * (iTriplet.size() - TINY) );
564 int iPickT = iTriplet[pickT];
565 cols[iPickT] = -acol0;
569 iTriplet[pickT] = iTriplet.back();
571 iDipCol.push_back(iPickT);
572 iDipAcol.push_back(0);
576 if (col0 < 0 && iAtriplet.size() > 0) {
577 int pickA = (iAtriplet.size() == 1) ? 0
578 :
int( TINY + rndmPtr->flat() * (iAtriplet.size() - TINY) );
579 int iPickA = iAtriplet[pickA];
580 acols[iPickA] = -col0;
584 iAtriplet[pickA] = iAtriplet.back();
585 iAtriplet.pop_back();
586 iDipCol.push_back(0);
587 iDipAcol.push_back(iPickA);
591 if ( (iTriplet.size() != iAtriplet.size())
592 || (col0 != 0 && acol0 == 0) || (col0 == 0 && acol0 != 0) ) {
593 infoPtr->errorMsg(
"Error in ResonanceDecays::pickColours:"
594 " inconsistent colour tags");
599 for (
int pickT = 0; pickT < int(iTriplet.size()); ++pickT) {
600 int iPickT = iTriplet[pickT];
601 int pickA = (iAtriplet.size() == 1) ? 0
602 :
int( TINY + rndmPtr->flat() * (iAtriplet.size() - TINY) );
603 int iPickA = iAtriplet[pickA];
606 cols[iPickT] = process.nextColTag();
607 acols[iPickA] = cols[iPickT];
610 iAtriplet[pickA] = iAtriplet.back();
611 iAtriplet.pop_back();
612 iDipCol.push_back(iPickT);
613 iDipAcol.push_back(iPickA);
617 if (col0 == 0 && acol0 == 0 && iOctet.size() == 0)
return true;
621 iDipCol.push_back(0);
622 iDipAcol.push_back(0);
626 for (
int i = 0; i < int(iOctet.size()); ++i) {
627 int iOct = iOctet[i];
630 if (iDipCol.size() == 0) {
631 cols[iOct] = process.nextColTag();
632 acols[iOct] = cols[iOct] ;
633 iDipCol.push_back(iOct);
634 iDipAcol.push_back(iOct);
639 int pickDip = (iDipCol.size() == 1) ? 0
640 :
int( TINY + rndmPtr->flat() * (iDipCol.size() - TINY) );
643 if (iDipCol[pickDip] == 0 && iDipAcol[pickDip] == 0) {
646 iDipAcol[pickDip] = iOct;
647 iDipCol.push_back(iOct);
648 iDipAcol.push_back(0);
651 }
else if (iDipAcol[pickDip] == 0) {
652 int iPickCol = iDipCol[pickDip];
653 cols[iOct] = cols[iPickCol];
654 acols[iOct] = process.nextColTag();
655 cols[iPickCol] = acols[iOct];
656 iDipCol[pickDip] = iOct;
657 iDipCol.push_back(iPickCol);
658 iDipAcol.push_back(iOct);
663 int iPickAcol = iDipAcol[pickDip];
664 acols[iOct] = acols[iPickAcol];
665 cols[iOct] = process.nextColTag();
666 acols[iPickAcol] = cols[iOct];
667 iDipAcol[pickDip] = iOct;
668 iDipCol.push_back(iOct);
669 iDipAcol.push_back(iPickAcol);
675 if (iDipCol.size() < 2) {
676 infoPtr->errorMsg(
"Error in ResonanceDecays::pickColours:"
677 " inconsistent colour tags");
691 bool ResonanceDecays::pickKinematics() {
698 double m1 = mProd[1];
699 double m2 = mProd[2];
702 double e1 = 0.5 * (m0*m0 + m1*m1 - m2*m2) / m0;
703 double e2 = 0.5 * (m0*m0 + m2*m2 - m1*m1) / m0;
704 double pAbs = 0.5 * sqrtpos( (m0 - m1 - m2) * (m0 + m1 + m2)
705 * (m0 + m1 - m2) * (m0 - m1 + m2) ) / m0;
708 double cosTheta = 2. * rndmPtr->flat() - 1.;
709 double sinTheta = sqrt(1. - cosTheta*cosTheta);
710 double phi = 2. * M_PI * rndmPtr->flat();
711 double pX = pAbs * sinTheta * cos(phi);
712 double pY = pAbs * sinTheta * sin(phi);
713 double pZ = pAbs * cosTheta;
716 pProd.push_back( Vec4( pX, pY, pZ, e1) );
717 pProd.push_back( Vec4( -pX, -pY, -pZ, e2) );
718 pProd[1].bst( pProd[0] );
719 pProd[2].bst( pProd[0] );
730 double m1 = mProd[1];
731 double m2 = mProd[2];
732 double m3 = mProd[3];
733 double mDiff = m0 - (m1 + m2 + m3);
736 double m23Min = m2 + m3;
737 double m23Max = m0 - m1;
738 double p1Max = 0.5 * sqrtpos( (m0 - m1 - m23Min) * (m0 + m1 + m23Min)
739 * (m0 + m1 - m23Min) * (m0 - m1 + m23Min) ) / m0;
740 double p23Max = 0.5 * sqrtpos( (m23Max - m2 - m3) * (m23Max + m2 + m3)
741 * (m23Max + m2 - m3) * (m23Max - m2 + m3) ) / m23Max;
742 double wtPSmax = 0.5 * p1Max * p23Max;
745 double wtPS, m23, p1Abs, p23Abs;
747 m23 = m23Min + rndmPtr->flat() * mDiff;
750 p1Abs = 0.5 * sqrtpos( (m0 - m1 - m23) * (m0 + m1 + m23)
751 * (m0 + m1 - m23) * (m0 - m1 + m23) ) / m0;
752 p23Abs = 0.5 * sqrtpos( (m23 - m2 - m3) * (m23 + m2 + m3)
753 * (m23 + m2 - m3) * (m23 - m2 + m3) ) / m23;
754 wtPS = p1Abs * p23Abs;
757 }
while ( wtPS < rndmPtr->flat() * wtPSmax );
760 double cosTheta = 2. * rndmPtr->flat() - 1.;
761 double sinTheta = sqrt(1. - cosTheta*cosTheta);
762 double phi = 2. * M_PI * rndmPtr->flat();
763 double pX = p23Abs * sinTheta * cos(phi);
764 double pY = p23Abs * sinTheta * sin(phi);
765 double pZ = p23Abs * cosTheta;
766 double e2 = sqrt( m2*m2 + p23Abs*p23Abs);
767 double e3 = sqrt( m3*m3 + p23Abs*p23Abs);
768 Vec4 p2( pX, pY, pZ, e2);
769 Vec4 p3( -pX, -pY, -pZ, e3);
772 cosTheta = 2. * rndmPtr->flat() - 1.;
773 sinTheta = sqrt(1. - cosTheta*cosTheta);
774 phi = 2. * M_PI * rndmPtr->flat();
775 pX = p1Abs * sinTheta * cos(phi);
776 pY = p1Abs * sinTheta * sin(phi);
777 pZ = p1Abs * cosTheta;
778 double e1 = sqrt( m1*m1 + p1Abs*p1Abs);
779 double e23 = sqrt( m23*m23 + p1Abs*p1Abs);
780 pProd.push_back( Vec4( pX, pY, pZ, e1) );
783 Vec4 p23( -pX, -pY, -pZ, e23);
786 pProd.push_back( p2 );
787 pProd.push_back( p3 );
788 pProd[1].bst( pProd[0] );
789 pProd[2].bst( pProd[0] );
790 pProd[3].bst( pProd[0] );
800 double mSum = mProd[1];
801 for (
int i = 2; i <= mult; ++i) mSum += mProd[i];
802 double mDiff = m0 - mSum;
806 for (
int i = 0; i <= mult; ++i) mInv.push_back( mProd[i]);
809 double wtPSmax = 1. / WTCORRECTION[mult];
810 double mMax = mDiff + mProd[mult];
812 for (
int i = mult - 1; i > 0; --i) {
815 double mNow = mProd[i];
816 wtPSmax *= 0.5 * sqrtpos( (mMax - mMin - mNow) * (mMax + mMin + mNow)
817 * (mMax + mMin - mNow) * (mMax - mMin + mNow) ) / mMax;
821 vector<double> rndmOrd;
828 rndmOrd.push_back(1.);
829 for (
int i = 1; i < mult - 1; ++i) {
830 double rndm = rndmPtr->flat();
831 rndmOrd.push_back(rndm);
832 for (
int j = i - 1; j > 0; --j) {
833 if (rndm > rndmOrd[j]) swap( rndmOrd[j], rndmOrd[j+1] );
837 rndmOrd.push_back(0.);
840 for (
int i = mult - 1; i > 0; --i) {
841 mInv[i] = mInv[i+1] + mProd[i] + (rndmOrd[i-1] - rndmOrd[i]) * mDiff;
842 wtPS *= 0.5 * sqrtpos( (mInv[i] - mInv[i+1] - mProd[i])
843 * (mInv[i] + mInv[i+1] + mProd[i]) * (mInv[i] + mInv[i+1] - mProd[i])
844 * (mInv[i] - mInv[i+1] + mProd[i]) ) / mInv[i];
848 }
while ( wtPS < rndmPtr->flat() * wtPSmax );
852 pInv.resize(mult + 1);
853 for (
int i = 1; i < mult; ++i) {
854 double pAbs = 0.5 * sqrtpos( (mInv[i] - mInv[i+1] - mProd[i])
855 * (mInv[i] + mInv[i+1] + mProd[i]) * (mInv[i] + mInv[i+1] - mProd[i])
856 * (mInv[i] - mInv[i+1] + mProd[i]) ) / mInv[i];
859 double cosTheta = 2. * rndmPtr->flat() - 1.;
860 double sinTheta = sqrt(1. - cosTheta*cosTheta);
861 double phi = 2. * M_PI * rndmPtr->flat();
862 double pX = pAbs * sinTheta * cos(phi);
863 double pY = pAbs * sinTheta * sin(phi);
864 double pZ = pAbs * cosTheta;
867 double eHad = sqrt( mProd[i]*mProd[i] + pAbs*pAbs);
868 double eInv = sqrt( mInv[i+1]*mInv[i+1] + pAbs*pAbs);
869 pProd.push_back( Vec4( pX, pY, pZ, eHad) );
870 pInv[i+1].p( -pX, -pY, -pZ, eInv);
872 pProd.push_back( pInv[mult] );
876 for (
int iFrame = mult - 1; iFrame > 0; --iFrame)
877 for (
int i = iFrame; i <= mult; ++i) pProd[i].bst(pInv[iFrame]);