1498{
1500
1501 const G4double aHugeValue = 1.0e+10;
1502
1503 #ifdef debugQGSParticipants
1506 #endif
1507
1509
1512
1513 G4double VqM_pr(0.), VaqM_pr(0.), VqM_tr(350.), VqqM_tr(700);
1515
1516 #ifdef debugQGSParticipants
1518 <<
"Target nucleon momenta at start"<<
G4endl;
1519 #endif
1520
1521 std::vector<G4VSplitableHadron*>::iterator i;
1523
1525 {
1526 Psum += (*i)->Get4Momentum();
1527 #ifdef debugQGSParticipants
1528 G4cout<<
"Nusleus nucleon # and its 4Mom. "<<NuclNo<<
" "<<(*i)->Get4Momentum()<<
G4endl;
1529 #endif
1530 NuclNo++;
1531 }
1532
1534
1536
1540 Projectile4Momentum.transform( toCms );
1541
1542
1543 #ifdef debugQGSParticipants
1545 G4cout<<
"Projectile 4 Mom "<<Projectile4Momentum<<
G4endl;
1546 #endif
1547
1548 NuclNo=0;
1551 {
1553 (*i)->Set4Momentum( tmp );
1554 #ifdef debugQGSParticipants
1555 G4cout<<
"Target nucleon # and 4Mom "<<
" "<<NuclNo<<
" "<<(*i)->Get4Momentum()<<
G4endl;
1556 #endif
1557 Target4Momentum += tmp;
1558 NuclNo++;
1559 }
1560
1563
1564 #ifdef debugQGSParticipants
1566 G4cout<<
"Target nucleons mom: px, py, z_1, m_i"<<
G4endl;
1567 #endif
1568
1569
1570 G4double PminusTarget = Target4Momentum.minus();
1571 NuclNo=0;
1572
1574 {
1576
1577
1578
1579
1580
1581
1582
1583
1586 if ( Mass2 < 0.0 ) {
1589 <<
"Â 4-momentum " << Psum <<
G4endl;
1590 ed <<
"LorentzVector tmp " << tmp <<
" with Mass2 " << Mass2 <<
G4endl;
1591 G4Exception(
"G4QGSParticipants::DeterminePartonMomenta(): 4-momentum with negative mass!",
1593 } else {
1594 Mass = std::sqrt( Mass2 );
1595 }
1596
1598 (*i)->Set4Momentum(tmp);
1599 #ifdef debugQGSParticipants
1600 G4cout<<
"Target nucleons # and mom: "<<NuclNo<<
" "<<(*i)->Get4Momentum()<<
G4endl;
1601 #endif
1602 NuclNo++;
1603 }
1604
1605
1606
1611
1613 G4double ProjSumMt(0.), ProjSumMt2perX(0.);
1614 G4double TargSumMt(0.), TargSumMt2perX(0.);
1615
1616
1618
1626
1628
1630 const G4int maxNumberOfAttempts = 1000;
1631 do
1632 {
1633 attempt++; if( attempt == 100*(attempt/100) ) {SigPt/=2.;}
1634
1635 ProjSumMt=0.; ProjSumMt2perX=0.;
1636 TargSumMt=0.; TargSumMt2perX=0.;
1637
1638 Success = true;
1640 #ifdef debugQGSParticipants
1641 G4cout<<
"attempt ------------------------ "<<attempt<<
G4endl;
1643 #endif
1644
1648 G4int NumberOfUnsampledSeaQuarks = 2*nSeaPair;
1649
1651 for (
G4int aSeaPair = 0; aSeaPair < nSeaPair; aSeaPair++)
1652 {
1654 #ifdef debugQGSParticipants
1655 G4cout<<
"Sea quarks: "<<aSeaPair<<
" "<<aParton->GetDefinition()->GetParticleName();
1656 #endif
1657 aPtVector = GaussianPt(SigPt, aHugeValue);
1658 tmp.
setPx(aPtVector.x()); tmp.
setPy(aPtVector.y());
1659 SumPx += aPtVector.x(); SumPy += aPtVector.y();
1660 Mt = std::sqrt(aPtVector.mag2()+
sqr(Qmass));
1661 ProjSumMt += Mt;
1662
1663
1664 tmp.
setPz(
SampleX(Xmin, NumberOfUnsampledSeaQuarks, 2*nSeaPair, aBeta)*(1.0-SumZ));
1666
1667 NumberOfUnsampledSeaQuarks--;
1668 ProjSumMt2perX +=
sqr(Mt)/tmp.
pz();
1670 aParton->Set4Momentum(tmp);
1671
1673 #ifdef debugQGSParticipants
1674 G4cout<<
" "<<aParton->GetDefinition()->GetParticleName()<<
G4endl;
1676 #endif
1677 aPtVector = GaussianPt(SigPt, aHugeValue);
1678 tmp.
setPx(aPtVector.x()); tmp.
setPy(aPtVector.y());
1679 SumPx += aPtVector.x(); SumPy += aPtVector.y();
1680 Mt = std::sqrt(aPtVector.mag2()+
sqr(Qmass));
1681 ProjSumMt += Mt;
1682
1683
1684 tmp.
setPz(
SampleX(Xmin, NumberOfUnsampledSeaQuarks, 2*nSeaPair, aBeta)*(1.0-SumZ));
1686
1687 NumberOfUnsampledSeaQuarks--;
1688 ProjSumMt2perX +=
sqr(Mt)/tmp.
pz();
1690 aParton->Set4Momentum(tmp);
1691 #ifdef debugQGSParticipants
1693 #endif
1694 }
1695
1696
1698 #ifdef debugQGSParticipants
1699 G4cout<<
"Val quark of Pr"<<
" "<<aParton->GetDefinition()->GetParticleName();
1700 #endif
1701 aPtVector = GaussianPt(SigPt, aHugeValue);
1702 tmp.
setPx(aPtVector.x()); tmp.
setPy(aPtVector.y());
1703 SumPx += aPtVector.x(); SumPy += aPtVector.y();
1704 Mt = std::sqrt(aPtVector.mag2()+
sqr(VqM_pr));
1705 ProjSumMt += Mt;
1706
1707
1708 tmp.
setPz(
SampleX(Xmin, NumberOfUnsampledSeaQuarks, 2*nSeaPair, aBeta)*(1.0-SumZ));
1710
1711 ProjSumMt2perX +=
sqr(Mt)/tmp.
pz();
1713 aParton->Set4Momentum(tmp);
1714
1715
1717 #ifdef debugQGSParticipants
1718 G4cout<<
" "<<aParton->GetDefinition()->GetParticleName()<<
G4endl;
1720 #endif
1722
1723 Mt = std::sqrt(
sqr(SumPx) +
sqr(SumPy) +
sqr(VaqM_pr) );
1724 ProjSumMt += Mt;
1726
1727 ProjSumMt2perX +=
sqr(Mt)/tmp.
pz();
1729 aParton->Set4Momentum(tmp);
1730 #ifdef debugQGSParticipants
1731 G4cout<<
" "<<tmp<<
" "<<SumZ+(1.-SumZ)<<
" (z-fraction)"<<
G4endl;
1732 #endif
1733
1734
1735
1736
1737
1738 NuclNo=0;
1740 {
1741 nSeaPair = (*i)->GetSoftCollisionCount()-1;
1742 #ifdef debugQGSParticipants
1744 <<
"Target nucleon 4Mom "<<(*i)->Get4Momentum()<<
G4endl;
1745 #endif
1746
1747 SumPx = (*i)->Get4Momentum().px() * (-1.);
1748 SumPy = (*i)->Get4Momentum().py() * (-1.);
1749 SumZ = 0.;
1750
1752 NumberOfUnsampledSeaQuarks = 2*nSeaPair;
1753
1754 Qmass=0;
1755 for (
G4int aSeaPair = 0; aSeaPair < nSeaPair; aSeaPair++)
1756 {
1757 aParton = (*i)->GetNextParton();
1758 #ifdef debugQGSParticipants
1759 G4cout<<
"Sea quarks: "<<aSeaPair<<
" "<<aParton->GetDefinition()->GetParticleName();
1760 #endif
1761 aPtVector = GaussianPt(SigPt, aHugeValue);
1762 tmp.
setPx(aPtVector.x()); tmp.
setPy(aPtVector.y());
1763 SumPx += aPtVector.x(); SumPy += aPtVector.y();
1764 Mt=std::sqrt(aPtVector.mag2()+
sqr(Qmass));
1765 TargSumMt += Mt;
1766
1767
1768 tmp.
setPz(
SampleX(Xmin, NumberOfUnsampledSeaQuarks, 2*nSeaPair, aBeta)*(1.0-SumZ));
1770 tmp.
setPz((*i)->Get4Momentum().pz()*tmp.
pz());
1772 NumberOfUnsampledSeaQuarks--;
1773 TargSumMt2perX +=
sqr(Mt)/tmp.
pz();
1775 aParton->Set4Momentum(tmp);
1776
1777 aParton = (*i)->GetNextAntiParton();
1778 #ifdef debugQGSParticipants
1779 G4cout<<
" "<<aParton->GetDefinition()->GetParticleName()<<
G4endl;
1781 #endif
1782 aPtVector = GaussianPt(SigPt, aHugeValue);
1783 tmp.
setPx(aPtVector.x()); tmp.
setPy(aPtVector.y());
1784 SumPx += aPtVector.x(); SumPy += aPtVector.y();
1785 Mt=std::sqrt(aPtVector.mag2()+
sqr(Qmass));
1786 TargSumMt += Mt;
1787
1788
1789 tmp.
setPz(
SampleX(Xmin, NumberOfUnsampledSeaQuarks, 2*nSeaPair, aBeta)*(1.0-SumZ));
1791 tmp.
setPz((*i)->Get4Momentum().pz()*tmp.
pz());
1793 NumberOfUnsampledSeaQuarks--;
1794 TargSumMt2perX +=
sqr(Mt)/tmp.
pz();
1796 aParton->Set4Momentum(tmp);
1797 #ifdef debugQGSParticipants
1799 #endif
1800 }
1801
1802
1803 aParton = (*i)->GetNextParton();
1804 #ifdef debugQGSParticipants
1805 G4cout<<
"Val quark of Tr"<<
" "<<aParton->GetDefinition()->GetParticleName();
1806 #endif
1807 aPtVector = GaussianPt(SigPt, aHugeValue);
1808 tmp.
setPx(aPtVector.x()); tmp.
setPy(aPtVector.y());
1809 SumPx += aPtVector.x(); SumPy += aPtVector.y();
1810 Mt=std::sqrt(aPtVector.mag2()+
sqr(VqM_tr));
1811 TargSumMt += Mt;
1812
1813
1814 tmp.
setPz(
SampleX(Xmin, NumberOfUnsampledSeaQuarks, 2*nSeaPair, aBeta)*(1.0-SumZ));
1816 tmp.
setPz((*i)->Get4Momentum().pz()*tmp.
pz());
1818 TargSumMt2perX +=
sqr(Mt)/tmp.
pz();
1820 aParton->Set4Momentum(tmp);
1821
1822
1823 aParton = (*i)->GetNextAntiParton();
1824 #ifdef debugQGSParticipants
1825 G4cout<<
" "<<aParton->GetDefinition()->GetParticleName()<<
G4endl;
1826 G4cout<<
" "<<tmp<<
" "<<SumZw<<
" (sum z-fracs) "<<SumZ<<
" (total z-sum) "<<
G4endl;
1827 #endif
1829
1830 Mt=std::sqrt(
sqr(SumPx) +
sqr(SumPy) +
sqr(VqqM_tr) );
1831 TargSumMt += Mt;
1832
1833 tmp.
setPz((*i)->Get4Momentum().pz()*(1.0 - SumZ));
1835 TargSumMt2perX +=
sqr(Mt)/tmp.
pz();
1837 aParton->Set4Momentum(tmp);
1838 #ifdef debugQGSParticipants
1839 G4cout<<
" "<<tmp<<
" "<<SumZw<<
" "<<1.0<<
" "<<(*i)->Get4Momentum().
pz()<<
G4endl;
1840 #endif
1841
1842 }
1843
1844 if( ProjSumMt + TargSumMt > SqrtS ) {
1845 Success = false; continue;}
1846 if( std::sqrt(ProjSumMt2perX) + std::sqrt(TargSumMt2perX) > SqrtS ) {
1847 Success = false; continue;}
1848
1849 } while( (!Success) &&
1850 attempt < maxNumberOfAttempts );
1851
1852 if ( attempt >= maxNumberOfAttempts ) {
1853 return false;
1854 }
1855
1856
1857
1859 - 2.0*
S*ProjSumMt2perX - 2.0*
S*TargSumMt2perX - 2.0*ProjSumMt2perX*TargSumMt2perX;
1860
1861 G4double targetWminus=(
S - ProjSumMt2perX + TargSumMt2perX + std::sqrt( DecayMomentum2 ))/2.0/SqrtS;
1862 G4double projectileWplus = SqrtS - TargSumMt2perX/targetWminus;
1863
1866
1868 #ifdef debugQGSParticipants
1869 G4cout<<
"Backward transformation ===================="<<
G4endl;
1871 #endif
1872
1873 for (
G4int aSeaPair = 0; aSeaPair < nSeaPair; aSeaPair++)
1874 {
1876 #ifdef debugQGSParticipants
1877 G4cout<<
"Sea quarks: "<<aSeaPair<<
" "<<aParton->GetDefinition()->GetParticleName();
1878 #endif
1879 Tmp =aParton->Get4Momentum(); z=Tmp.z();
1880
1881 Tmp.setPz(projectileWplus*z/2.0 - Tmp.e()/(2.0*z*projectileWplus));
1882 Tmp.setE( projectileWplus*z/2.0 + Tmp.e()/(2.0*z*projectileWplus));
1883 Tmp.transform( toLab );
1884
1885 aParton->Set4Momentum(Tmp);
1886
1888 #ifdef debugQGSParticipants
1889 G4cout<<
" "<<aParton->GetDefinition()->GetParticleName()<<
G4endl;
1891 #endif
1892 Tmp =aParton->Get4Momentum(); z=Tmp.z();
1893 Tmp.setPz(projectileWplus*z/2.0 - Tmp.e()/(2.0*z*projectileWplus));
1894 Tmp.setE( projectileWplus*z/2.0 + Tmp.e()/(2.0*z*projectileWplus));
1895 Tmp.transform( toLab );
1896
1897 aParton->Set4Momentum(Tmp);
1898 #ifdef debugQGSParticipants
1900 #endif
1901 }
1902
1903
1905 #ifdef debugQGSParticipants
1906 G4cout<<
"Val quark of Pr"<<
" "<<aParton->GetDefinition()->GetParticleName();
1907 #endif
1908 Tmp =aParton->Get4Momentum(); z=Tmp.z();
1909 Tmp.setPz(projectileWplus*z/2.0 - Tmp.e()/(2.0*z*projectileWplus));
1910 Tmp.setE( projectileWplus*z/2.0 + Tmp.e()/(2.0*z*projectileWplus));
1911 Tmp.transform( toLab );
1912
1913 aParton->Set4Momentum(Tmp);
1914
1915
1917 #ifdef debugQGSParticipants
1918 G4cout<<
" "<<aParton->GetDefinition()->GetParticleName()<<
G4endl;
1920 #endif
1921 Tmp =aParton->Get4Momentum(); z=Tmp.z();
1922 Tmp.setPz(projectileWplus*z/2.0 - Tmp.e()/(2.0*z*projectileWplus));
1923 Tmp.setE( projectileWplus*z/2.0 + Tmp.e()/(2.0*z*projectileWplus));
1924 Tmp.transform( toLab );
1925
1926 aParton->Set4Momentum(Tmp);
1927
1928 #ifdef debugQGSParticipants
1930 #endif
1931
1932
1933
1934
1935 NuclNo=0;
1937 {
1938 nSeaPair = (*i)->GetSoftCollisionCount()-1;
1939 #ifdef debugQGSParticipants
1940 G4cout<<
"nSeaPair of target and N# "<<nSeaPair<<
" "<<NuclNo<<
G4endl;
1941 #endif
1942 NuclNo++;
1943 for (
G4int aSeaPair = 0; aSeaPair < nSeaPair; aSeaPair++)
1944 {
1945 aParton = (*i)->GetNextParton();
1946 #ifdef debugQGSParticipants
1947 G4cout<<
"Sea quarks: "<<aSeaPair<<
" "<<aParton->GetDefinition()->GetParticleName();
1948 #endif
1949 Tmp =aParton->Get4Momentum(); z=Tmp.z();
1950 Tmp.setPz(-targetWminus*z/2.0 + Tmp.e()/(2.0*z*targetWminus));
1951 Tmp.setE( targetWminus*z/2.0 + Tmp.e()/(2.0*z*targetWminus));
1952 Tmp.transform( toLab );
1953
1954 aParton->Set4Momentum(Tmp);
1955
1956 aParton = (*i)->GetNextAntiParton();
1957 #ifdef debugQGSParticipants
1958 G4cout<<
" "<<aParton->GetDefinition()->GetParticleName()<<
G4endl;
1960 #endif
1961 Tmp =aParton->Get4Momentum(); z=Tmp.z();
1962 Tmp.setPz(-targetWminus*z/2.0 + Tmp.e()/(2.0*z*targetWminus));
1963 Tmp.setE( targetWminus*z/2.0 + Tmp.e()/(2.0*z*targetWminus));
1964 Tmp.transform( toLab );
1965
1966 aParton->Set4Momentum(Tmp);
1967 #ifdef debugQGSParticipants
1969 #endif
1970 }
1971
1972
1973
1974 aParton = (*i)->GetNextParton();
1975 #ifdef debugQGSParticipants
1976 G4cout<<
"Val quark of Tr"<<
" "<<aParton->GetDefinition()->GetParticleName();
1977 #endif
1978 Tmp =aParton->Get4Momentum(); z=Tmp.z();
1979 Tmp.setPz(-targetWminus*z/2.0 + Tmp.e()/(2.0*z*targetWminus));
1980 Tmp.setE( targetWminus*z/2.0 + Tmp.e()/(2.0*z*targetWminus));
1981 Tmp.transform( toLab );
1982
1983 aParton->Set4Momentum(Tmp);
1984
1985
1986 aParton = (*i)->GetNextAntiParton();
1987 #ifdef debugQGSParticipants
1988 G4cout<<
" "<<aParton->GetDefinition()->GetParticleName()<<
G4endl;
1990 #endif
1991 Tmp =aParton->Get4Momentum(); z=Tmp.z();
1992 Tmp.setPz(-targetWminus*z/2.0 + Tmp.e()/(2.0*z*targetWminus));
1993 Tmp.setE( targetWminus*z/2.0 + Tmp.e()/(2.0*z*targetWminus));
1994 Tmp.transform( toLab );
1995
1996 aParton->Set4Momentum(Tmp);
1997 #ifdef debugQGSParticipants
1999 #endif
2000 NuclNo++;
2001 }
2002
2003 return true;
2004}
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
std::ostringstream G4ExceptionDescription
Hep3Vector boostVector() const
HepLorentzVector & rotateZ(double)
HepLorentzVector & rotateY(double)
HepLorentzVector & transform(const HepRotation &)
static G4Gamma * GammaDefinition()
static G4KaonMinus * KaonMinusDefinition()
static G4KaonPlus * KaonPlusDefinition()
static G4PionMinus * PionMinusDefinition()
static G4PionPlus * PionPlusDefinition()
static G4PionZero * PionZeroDefinition()
virtual G4Parton * GetNextParton()
virtual G4Parton * GetNextAntiParton()
G4double SampleX(G4double anXmin, G4int nSea, G4int theTotalSea, G4double aBeta)
const G4ParticleDefinition * GetDefinition() const
const G4LorentzVector & Get4Momentum() const
G4int GetSoftCollisionCount()