72 for(
G4int i=0; i<3; i++)
73 {
for(
G4int j=0; j<3; j++)
74 {
for(
G4int k=0; k<6; k++)
166 for(
G4int i=0; i<3; i++)
167 {
for(
G4int j=0; j<3; j++)
168 {
for(
G4int k=0; k<3; k++)
169 {
for(
G4int l=0; l<4; l++)
398 for (
G4int i=0 ; i<35 ; i++ ) {
416 G4int Number_of_quarks=0;
417 G4int Number_of_squarks=0;
419 G4double StringM=
string->Get4Momentum().mag();
423 #ifdef debug_LUNDfragmentation
434 G4int q2=(Qleft/100)%10;
447 #ifdef debug_LUNDfragmentation
454 G4int q1=Qright/1000;
458 G4int q2=(Qright/100)%10;
471 #ifdef debug_LUNDfragmentation
476 if(Number_of_quarks==2){EstimatedMass += 70.*
MeV;}
478 if(Number_of_quarks==3)
480 if(Number_of_squarks==0) {EstimatedMass += 740.*
MeV;}
481 if(Number_of_squarks==1) {EstimatedMass += 740.*
MeV;}
482 if(Number_of_squarks==2) {EstimatedMass += 400.*
MeV;}
483 if(Number_of_squarks==3) {EstimatedMass += 382.*
MeV;}
485 if(Number_of_quarks==4)
487 if((StringM > 1880.) && ( EstimatedMass < 2100)) {EstimatedMass = 2020.;}
489 else if((StringM > 2232.) && ( EstimatedMass < 2730)){EstimatedMass = 2570.;}
490 else if((StringM > 5130.) && ( EstimatedMass < 3450)){EstimatedMass = 5130.;}
493 if(EstimatedMass <= 1600.*
MeV){EstimatedMass-=200.*
MeV;}
494 else {EstimatedMass+=100.*
MeV;}
498 #ifdef debug_LUNDfragmentation
522 #ifdef debug_LUNDfragmentation
525 <<
"------------------------------------"<<
G4endl;
539 #ifdef debug_LUNDfragmentation
540 G4cout<<
"Non fragmentable - the string is converted to one hadron "<<
G4endl;
548 LeftVector->operator[](0)->SetPosition(theString.
GetPosition());
550 if(LeftVector->size() > 1)
554 LeftVector->operator[](1)->SetPosition(theString.
GetPosition());
559 #ifdef debug_LUNDfragmentation
569 #ifdef debug_LUNDfragmentation
576 #ifdef debug_LUNDfragmentation
579 G4cout<<G4endl<<
"LUND StringFragm: String Mass " <<theStringInCMS->
Get4Momentum().mag()<<
" Pz Lab "
581 <<
"------------------------------------"<<
G4endl;
589 delete theStringInCMS;
601 while(!RightVector->empty())
603 LeftVector->push_back(RightVector->back());
604 RightVector->erase(RightVector->end()-1);
615 for(
size_t C1 = 0;
C1 < LeftVector->size();
C1++)
620 Momentum = toObserverFrame*Momentum;
624 Momentum = toObserverFrame*Coordinate;
627 Hadron->
SetPosition(PositionOftheStringCreation+aPosition);
664 #ifdef debug_LUNDfragmentation
665 G4cout<<
"Split last-----------------------------------------"<<
G4endl;
669 G4ThreeVector ClusterVel=
string->Get4Momentum().boostVector();
677 #ifdef debug_LUNDfragmentation
678 G4cout<<
"StrMass "<<StringMass<<
" q "
679 <<
string->GetLeftParton()->GetParticleName()<<
" "
680 <<
string->GetRightParton()->GetParticleName()<<
G4endl;
683 string->SetLeftPartonStable();
713 if (string->
DecayIsQuark() &&
string->StableIsQuark() ) {
715 #ifdef debug_LUNDfragmentation
722 #ifdef debug_LUNDfragmentation
734 #ifdef debug_LUNDfragmentation
735 G4cout<<
"Selected LeftHad RightHad "<<sampledState<<
" "
746 LeftMom.boost(ClusterVel);
747 RightMom.boost(ClusterVel);
749 LeftVector->push_back(
new G4KineticTrack(LeftHadron, 0, Pos, LeftMom));
750 RightVector->push_back(
new G4KineticTrack(RightHadron, 0, Pos, RightMom));
766 #ifdef debug_LUNDfragmentation
767 G4cout<<
"Sampling of momenta of 2 last produced hadrons ----------------"<<
G4endl;
768 G4cout<<
"Masses "<<InitialMass<<
" "<<Mass<<
" "<<AntiMass<<
G4endl;
771 if((Mass > 930. || AntiMass > 930.) && (
G4UniformRand() < ProbIsotropy)) {
773 G4double r_val =
sqr(InitialMass*InitialMass - Mass*Mass - AntiMass*AntiMass) -
774 sqr(2.*Mass*AntiMass);
775 G4double Pabs = (r_val > 0.)? std::sqrt(r_val)/(2.*InitialMass) : 0;
780 G4double st = std::sqrt(1. - pz * pz)*Pabs;
786 Mom->setPx(px); Mom->setPy(py); Mom->setPz(pz);
787 Mom->setE(std::sqrt(Pabs*Pabs + Mass*Mass));
789 AntiMom->setPx(-px); AntiMom->setPy(-py); AntiMom->setPz(-pz);
790 AntiMom->setE (std::sqrt(Pabs*Pabs + AntiMass*AntiMass));
793 const G4int maxNumberOfLoops = 1000;
794 G4int loopCounter = 0;
799 G4double termD = InitialMass*InitialMass -Mass*Mass - AntiMass*AntiMass;
801 G4double termN = 2*termD + 4*Mass*Mass + 4*AntiMass*AntiMass;
802 G4double pt2max=(termD*termD - termab )/ termN ;
807 MassMt2 = Mass * Mass + Pt2;
808 AntiMassMt2= AntiMass * AntiMass + Pt2;
810 AvailablePz2=
sqr(InitialMass*InitialMass - MassMt2 - AntiMassMt2) -
811 4.*MassMt2*AntiMassMt2;
813 }
while( (AvailablePz2 < 0.) &&
814 ++loopCounter < maxNumberOfLoops );
816 if ( loopCounter >= maxNumberOfLoops ) {
820 AvailablePz2 /=(4.*InitialMass*InitialMass);
821 AvailablePz = std::sqrt(AvailablePz2);
826 Mom->setPx(Px); Mom->setPy(Py); Mom->setPz(AvailablePz);
827 Mom->setE(std::sqrt(MassMt2+AvailablePz2));
829 AntiMom->setPx(-Px); AntiMom->setPy(-Py); AntiMom->setPz(-AvailablePz);
830 AntiMom->setE (std::sqrt(AntiMassMt2+AvailablePz2));
840 G4double StringMT2=
string->MassT2();
841 G4double StringMT =std::sqrt(StringMT2);
847 #ifdef debug_LUNDfragmentation
849 G4cout<<
"String 4 mom, String M and Mt "<<String4Momentum<<
" "<<String4Momentum.mag()<<
" "<<std::sqrt(StringMT2)<<
G4endl;
857 #ifdef debug_LUNDfragmentation
858 G4cout<<
"Mass of the string is not sufficient to produce the hadron!"<<
G4endl;
864 String4Momentum.setPz(0.);
869 G4double HadronMassT2, ResidualMassT2;
879 RemSysPt = StringPt - HadronPt;
881 HadronMassT2 =
sqr(HadronMass) + HadronPt.mag2();
884 }
while(std::sqrt(HadronMassT2) + std::sqrt(ResidualMassT2) > StringMT);
889 G4double Pz2 = (
sqr(StringMT2 - HadronMassT2 - ResidualMassT2) -
890 4*HadronMassT2 * ResidualMassT2)/4./StringMT2;
892 if(Pz2 < 0 ) {
return 0;}
897 G4double zMin = (std::sqrt(HadronMassT2+Pz2) - Pz)/std::sqrt(StringMT2);
899 G4double zMax = (std::sqrt(HadronMassT2+Pz2) + Pz)/std::sqrt(StringMT2);
901 if (zMin >= zMax)
return 0;
904 pHadron, HadronPt.x(), HadronPt.y());
910 (z *
string->LightConeDecay() - HadronMassT2/(z *
string->LightConeDecay())));
911 G4double HadronE = 0.5* (z *
string->LightConeDecay() +
912 HadronMassT2/(z *
string->LightConeDecay()));
916 #ifdef debug_LUNDfragmentation
917 G4cout<<
"string->LightConeDecay() "<<
string->LightConeDecay()<<
G4endl;
918 G4cout<<
"HadronPt,HadronE "<<HadronPt<<
" "<<HadronE<<
G4endl;
935 G4double Mt2 = Px*Px + Py*Py + Mass*Mass;
938 G4double zOfMaxyf(0.), maxYf(1.), z(0.), yf(1.);
939 if(std::abs(PDGEncodingOfDecayParton) < 1000)
945 zOfMaxyf=alund*Mt2/(alund*Mt2 + 1.);
946 maxYf=(1-zOfMaxyf)/zOfMaxyf *
G4Exp(-alund*Mt2/zOfMaxyf);
948 const G4int maxNumberOfLoops = 1000;
949 G4int loopCounter = 0;
953 yf = (1-z)/z *
G4Exp(-alund*Mt2/z);
955 }
while ( (
G4UniformRand()*maxYf > yf) && ++loopCounter < maxNumberOfLoops );
956 if ( loopCounter >= maxNumberOfLoops ) {
957 z = 0.5*(zmin + zmax);
962 if(std::abs(PDGEncodingOfDecayParton) > 1000)
1002 #ifdef debug_LUNDfragmentation
1008 G4bool final_success=
false;
1009 G4bool inner_success=
true;
1018 LeftVector->clear();
1020 RightVector->clear();
1024 const G4int maxNumberOfLoops = 1000;
1025 G4int loopCounter = -1;
1026 while ( (!
StopFragmenting(currentString)) && ++loopCounter < maxNumberOfLoops )
1028 #ifdef debug_LUNDfragmentation
1036 #ifdef debug_LUNDfragmentation
1042 LeftVector->push_back(Hadron);
1044 RightVector->push_back(Hadron);
1046 delete currentString;
1047 currentString=newString;
1050 if ( loopCounter >= maxNumberOfLoops ) {
1051 inner_success=
false;
1055 #ifdef debug_LUNDfragmentation
1056 G4cout<<
"Split remaining string into 2 final hadrons."<<
G4endl;
1059 if ( inner_success &&
SplitLast(currentString, LeftVector, RightVector) )
1065 delete currentString;
1067 return final_success;
1076 G4double StringMass =
string->Mass();
1077 G4int cClusterInterrupt = 0;
1083 G4int LeftQuark1=
string->GetLeftParton()->GetPDGEncoding()/1000;
1084 G4int LeftQuark2=(
string->GetLeftParton()->GetPDGEncoding()/100)%10;
1086 G4int RightQuark1=
string->GetRightParton()->GetPDGEncoding()/1000;
1087 G4int RightQuark2=(
string->GetRightParton()->GetPDGEncoding()/100)%10;
1100 while ((StringMass <= LeftHadron->GetPDGMass() + RightHadron->
GetPDGMass()));
1113 G4double StringMass =
string->Mass();
1119 Anti_Di_Quark =
string->GetLeftParton();
1120 Di_Quark=
string->GetRightParton();
1122 Anti_Di_Quark =
string->GetRightParton();
1123 Di_Quark=
string->GetLeftParton();
1127 G4int AbsIDAnti_di_quark =std::abs(IDAnti_di_quark);
1129 G4int AbsIDdi_quark =std::abs(IDdi_quark);
1131 G4int ADi_q1=AbsIDAnti_di_quark/1000;
1132 G4int ADi_q2=(AbsIDAnti_di_quark-ADi_q1*1000)/100;
1134 G4int Di_q1=AbsIDdi_quark/1000;
1135 G4int Di_q2=(AbsIDdi_quark-Di_q1*1000)/100;
1138 for(
G4int ProdQ=1; ProdQ < 4; ProdQ++)
1141 const G4int maxNumberOfLoops = 1000;
1142 G4int loopCounter = 0;
1146 -
Baryon[ADi_q1-1][ADi_q2-1][ProdQ-1][StateADiQ]);
1152 const G4int maxNumberOfInternalLoops = 1000;
1153 G4int internalLoopCounter = 0;
1157 +
Baryon[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]);
1160 if(StringMass > LeftHadronMass + RightHadronMass)
1165 G4Exception(
"G4LundStringFragmentation::Diquark_AntiDiquark_aboveThreshold_lastSplitting ",
1185 }
while( (
Baryon[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]!=0) &&
1186 ++internalLoopCounter < maxNumberOfInternalLoops );
1187 if ( internalLoopCounter >= maxNumberOfInternalLoops ) {
1192 }
while( (
Baryon[ADi_q1-1][ADi_q2-1][ProdQ-1][StateADiQ]!=0) &&
1193 ++loopCounter < maxNumberOfLoops );
1194 if ( loopCounter >= maxNumberOfLoops ) {
1208 G4double StringMass =
string->Mass();
1215 Quark =
string->GetLeftParton();
1216 Di_Quark=
string->GetRightParton();
1218 Quark =
string->GetRightParton();
1219 Di_Quark=
string->GetLeftParton();
1223 G4int AbsIDquark =std::abs(IDquark);
1225 G4int AbsIDdi_quark=std::abs(IDdi_quark);
1226 G4int Di_q1=AbsIDdi_quark/1000;
1227 G4int Di_q2=(AbsIDdi_quark-Di_q1*1000)/100;
1230 if(IDdi_quark < 0) SignDiQ=-1;
1233 for(
G4int ProdQ=1; ProdQ < 4; ProdQ++)
1238 if(IDquark == 2) SignQ= 1;
1239 if((IDquark == 1) && (ProdQ == 3)) SignQ= 1;
1240 if((IDquark == 3) && (ProdQ == 1)) SignQ=-1;
1243 if(IDquark == -2) SignQ=-1;
1244 if((IDquark ==-1) && (ProdQ == 3)) SignQ=-1;
1245 if((IDquark ==-3) && (ProdQ == 1)) SignQ= 1;
1248 if(AbsIDquark == ProdQ) SignQ= 1;
1254 const G4int maxNumberOfLoops = 1000;
1255 G4int loopCounter = 0;
1259 SignQ*
Meson[AbsIDquark-1][ProdQ-1][StateQ]);
1263 const G4int maxNumberOfInternalLoops = 1000;
1264 G4int internalLoopCounter = 0;
1268 SignDiQ*
Baryon[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]);
1271 if(StringMass > LeftHadronMass + RightHadronMass)
1276 G4Exception(
"G4LundStringFragmentation::Quark_Diquark_lastSplitting ",
1294 }
while( (
Baryon[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]!=0) &&
1295 ++internalLoopCounter < maxNumberOfInternalLoops );
1296 if ( internalLoopCounter >= maxNumberOfInternalLoops ) {
1301 }
while( (
Meson[AbsIDquark-1][ProdQ-1][StateQ]!=0) &&
1302 ++loopCounter < maxNumberOfLoops );
1303 if ( loopCounter >= maxNumberOfLoops ) {
1317 G4double StringMass =
string->Mass();
1324 Quark =
string->GetLeftParton();
1325 Anti_Quark=
string->GetRightParton();
1327 Quark =
string->GetRightParton();
1328 Anti_Quark=
string->GetLeftParton();
1332 G4int AbsIDquark =std::abs(IDquark);
1334 G4int AbsIDanti_quark=std::abs(IDanti_quark);
1337 for(
G4int ProdQ=1; ProdQ < 4; ProdQ++)
1340 if(IDquark == 2) SignQ= 1;
1341 if((IDquark == 1) && (ProdQ == 3)) SignQ= 1;
1342 if((IDquark == 3) && (ProdQ == 1)) SignQ=-1;
1343 if(IDquark == ProdQ) SignQ= 1;
1346 if(IDanti_quark == -2) SignAQ=-1;
1347 if((IDanti_quark ==-1) && (ProdQ == 3)) SignAQ=-1;
1348 if((IDanti_quark ==-3) && (ProdQ == 1)) SignAQ= 1;
1349 if(AbsIDanti_quark == ProdQ) SignAQ= 1;
1352 const G4int maxNumberOfLoops = 1000;
1353 G4int loopCounter = 0;
1357 SignQ*
Meson[AbsIDquark-1][ProdQ-1][StateQ]);
1361 const G4int maxNumberOfInternalLoops = 1000;
1362 G4int internalLoopCounter = 0;
1366 SignAQ*
Meson[AbsIDanti_quark-1][ProdQ-1][StateAQ]);
1369 if(StringMass > LeftHadronMass + RightHadronMass)
1374 G4Exception(
"G4LundStringFragmentation::Quark_AntiQuark_lastSplitting ",
1397 }
while( (
Meson[AbsIDanti_quark-1][ProdQ-1][StateAQ]!=0) &&
1398 ++internalLoopCounter < maxNumberOfInternalLoops );
1399 if ( internalLoopCounter >= maxNumberOfInternalLoops ) {
1404 }
while( (
Meson[AbsIDquark-1][ProdQ-1][StateQ]!=0) &&
1405 ++loopCounter < maxNumberOfLoops );
1406 if ( loopCounter >= maxNumberOfLoops ) {
1429 G4int indexPosition = 0;
1435 if(Sum >= ksi)
break;
1437 return indexPosition;
1455 G4int Swap = stableQuarkEncoding;
1456 stableQuarkEncoding = decayQuarkEncoding;
1457 decayQuarkEncoding = Swap;
1460 G4int IsParticle=(decayQuarkEncoding>0) ? -1 : +1;
1470 G4int QuarkEncoding=QuarkPair.second->GetPDGEncoding();
1471 G4int i10 =
std::max(std::abs(QuarkEncoding), std::abs(stableQuarkEncoding));
1472 G4int i20 =
std::min(std::abs(QuarkEncoding), std::abs(stableQuarkEncoding));
1474 G4int NewDecayEncoding = -1*IsParticle*(i10 * 1000 + i20 * 100 + spin);
1491 created = QuarkPair.second;
virtual G4KineticTrackVector * FragmentString(const G4ExcitedString &theString)
G4double Mass_of_heavy_quark
G4ParticleDefinition * GetRightParton(void) const
G4ParticleDefinition * FindParticle(G4int PDGEncoding)
G4bool Quark_AntiQuark_lastSplitting(G4FragmentingString *&string, G4ParticleDefinition *&LeftHadron, G4ParticleDefinition *&RightHadron)
std::ostringstream G4ExceptionDescription
CLHEP::Hep3Vector G4ThreeVector
virtual G4LorentzVector * SplitEandP(G4ParticleDefinition *pHadron, G4FragmentingString *string, G4FragmentingString *newString)
CLHEP::HepLorentzRotation G4LorentzRotation
G4double Mass_of_string_junction
const G4ThreeVector & GetPosition() const
G4ExcitedString * CPExcited(const G4ExcitedString &string)
virtual G4bool IsFragmentable(const G4FragmentingString *const string)
void SetFormationTime(G4double aFormationTime)
void SetStrangenessSuppression(G4double aValue)
const G4LorentzVector & Get4Momentum() const
G4double MesonWeight[3][3][6]
G4ParticleDefinition * Build(G4ParticleDefinition *black, G4ParticleDefinition *white)
G4bool Quark_Diquark_lastSplitting(G4FragmentingString *&string, G4ParticleDefinition *&LeftHadron, G4ParticleDefinition *&RightHadron)
G4int StringLoopInterrupt
G4Parton * GetLeftParton(void) const
G4int ClusterLoopInterrupt
G4int GetPDGEncoding() const
G4ThreeVector SampleQuarkPt(G4double ptMax=-1.)
std::vector< G4double > vectorMesonMix
const G4String & GetParticleSubType() const
virtual void SetMassCut(G4double aValue)
const G4String & GetParticleName() const
G4ParticleDefinition * GetDecayParton() const
virtual G4double GetLightConeZ(G4double zmin, G4double zmax, G4int PartonEncoding, G4ParticleDefinition *pHadron, G4double Px, G4double Py)
pDefPair CreatePartonPair(G4int NeedParticle, G4bool AllowDiquarks=true)
void SetStringTensionParameter(G4double aValue)
const G4ThreeVector & GetPosition() const
virtual void Sample4Momentum(G4LorentzVector *Mom, G4double Mass, G4LorentzVector *AntiMom, G4double AntiMass, G4double InitialMass)
G4ParticleDefinition * GetLeftParton(void) const
void SetDiquarkBreakProbability(G4double aValue)
G4double GetStrangeSuppress()
void SetMinimalStringMass2(const G4double aValue)
G4GLOB_DLL std::ostream G4cout
ParticleList decay(Cluster *const c)
Carries out a cluster decay.
G4KineticTrackVector * LightFragmentationTest(const G4ExcitedString *const theString)
G4LorentzVector Get4Momentum() const
G4bool Loop_toFragmentString(G4ExcitedString *&theStringInCMS, G4KineticTrackVector *&LeftVector, G4KineticTrackVector *&RightVector)
G4double GetFormationTime() const
std::pair< G4ParticleDefinition *, G4ParticleDefinition * > pDefPair
void SetMinimalStringMass(const G4FragmentingString *const string)
G4LorentzRotation TransformToAlignedCms()
G4ParticleDefinition * FindParticle(G4int Encoding)
void SetPosition(const G4ThreeVector aPosition)
virtual G4bool StopFragmenting(const G4FragmentingString *const string)
G4int GetDecayDirection() const
void Set4Momentum(const G4LorentzVector &a4Momentum)
G4LundStringFragmentation()
G4double lambda(G4double s, G4double m1_Sqr, G4double m2_Sqr)
G4bool FourQuarkString(void) const
std::vector< G4double > scalarMesonMix
virtual G4ParticleDefinition * DiQuarkSplitup(G4ParticleDefinition *decay, G4ParticleDefinition *&created)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
G4Parton * GetRightParton(void) const
virtual ~G4LundStringFragmentation()
G4HadronBuilder * hadronizer
G4double GetPDGMass() const
G4bool Diquark_AntiDiquark_aboveThreshold_lastSplitting(G4FragmentingString *&string, G4ParticleDefinition *&LeftHadron, G4ParticleDefinition *&RightHadron)
static G4ParticleTable * GetParticleTable()
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4double MinimalStringMass
G4double Mass_of_light_quark
G4double GetTimeOfCreation() const
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
static constexpr double GeV
G4double DiquarkBreakProb
G4int GetDirection(void) const
G4ParticleDefinition * FS_LeftHadron[35]
static constexpr double MeV
G4bool Diquark_AntiDiquark_belowThreshold_lastSplitting(G4FragmentingString *&string, G4ParticleDefinition *&LeftHadron, G4ParticleDefinition *&RightHadron)
static constexpr double pi
virtual G4bool SplitLast(G4FragmentingString *string, G4KineticTrackVector *LeftVector, G4KineticTrackVector *RightVector)
static constexpr double fermi
const G4LorentzVector & Get4Momentum() const
void CalculateHadronTimePosition(G4double theInitialStringMass, G4KineticTrackVector *)
const G4ParticleDefinition * GetDefinition() const
G4double MinimalStringMass2
G4double BaryonWeight[3][3][3][4]
G4KineticTrack * Splitup(G4FragmentingString *string, G4FragmentingString *&newString)
G4ParticleDefinition * FS_RightHadron[35]
void SetDiquarkSuppression(G4double aValue)
CLHEP::HepLorentzVector G4LorentzVector