81 #ifdef debug_LUNDfragmentation
82 G4cout<<
G4endl<<
"LUND StringFragmentation ------------------------------------"<<
G4endl;
83 G4cout<<G4endl<<
"LUND StringFragm: String Mass "
86 <<
"------------------------------------"<<
G4endl;
99 #ifdef debug_LUNDfragmentation
100 G4cout<<
"Non fragmentable - the string is converted to one hadron "<<
G4endl;
111 LeftVector->operator[](0)->SetPosition(theString.
GetPosition());
113 if (LeftVector->size() > 1)
117 LeftVector->operator[](1)->SetPosition(theString.
GetPosition());
122 #ifdef debug_LUNDfragmentation
142 while (!RightVector->empty())
144 LeftVector->push_back(RightVector->back());
145 RightVector->erase(RightVector->end()-1);
170 #ifdef debug_LUNDfragmentation
178 G4bool final_success=
false;
179 G4bool inner_success=
true;
196 RightVector->clear();
200 const G4int maxNumberOfLoops = 1000;
201 G4int loopCounter = -1;
203 while ( (!
StopFragmenting(currentString)) && ++loopCounter < maxNumberOfLoops )
205 #ifdef debug_LUNDfragmentation
212 toObserverFrame= toCms.
inverse();
213 #ifdef debug_LUNDfragmentation
223 #ifdef debug_LUNDfragmentation
237 Hadron->
SetPosition(PositionOftheStringCreation+aPosition);
242 LeftVector->push_back(Hadron);
245 RightVector->push_back(Hadron);
247 delete currentString;
248 currentString=newString;
250 if ( newString )
delete newString;
256 if ( loopCounter >= maxNumberOfLoops ) {
261 #ifdef debug_LUNDfragmentation
262 if (inner_success)
G4cout<<
"Split remaining string into 2 final hadrons."<<
G4endl;
264 if ( inner_success &&
SplitLast(currentString, LeftVector, RightVector) )
266 final_success =
true;
269 delete currentString;
271 return final_success;
293 #ifdef debug_LUNDfragmentation
295 <<
" "<<
string->Mass()<<
G4endl;
307 #ifdef debug_LUNDfragmentation
309 G4cout<<
"Start SplitUP ========================="<<
G4endl;
310 G4cout<<
"String partons: " <<
string->GetLeftParton()->GetPDGEncoding()<<
" "
311 <<
string->GetRightParton()->GetPDGEncoding()<<
" "
312 <<
"Direction " <<
string->GetDecayDirection()<<
G4endl;
319 string->SetLeftPartonStable();
322 string->SetRightPartonStable();
333 #ifdef debug_LUNDfragmentation
337 G4int NumberOfpossibleBaryons = 2;
343 ActualProb *= (1.0-
sqr(NumberOfpossibleBaryons*1400.0/StringMass));
348 if ( NumberOfpossibleBaryons == 3 ){Mth = 2520.0;}
349 else if ( NumberOfpossibleBaryons == 4 ){Mth = 2380.0;}
353 if ( ActualProb < 0.0 ) ActualProb = 0.0;
356 #ifdef debug_LUNDfragmentation
370 if ( HadronDefinition == NULL ) {
G4KineticTrack * Hadron =0;
return Hadron; }
372 #ifdef debug_LUNDfragmentation
373 G4cout<<
"The parton "<<
string->GetDecayParton()->GetPDGEncoding()<<
" "
376 G4cout<<
"The side of the string decay Left/Right (1/-1) "<<SideOfDecay<<
G4endl;
380 if ( newString )
delete newString;
384 #ifdef debug_LUNDfragmentation
385 G4cout<<
"An attempt to determine its energy (SplitEandP)"<<
G4endl;
389 delete newString; newString=0;
392 if ( HadronMomentum != 0 ) {
394 #ifdef debug_LUNDfragmentation
398 Hadron =
new G4KineticTrack(HadronDefinition, 0,Pos, *HadronMomentum);
400 if ( newString )
delete newString;
404 delete HadronMomentum;
408 #ifdef debug_LUNDfragmentation
413 #ifdef debug_LUNDfragmentation
414 G4cout<<
"End SplitUP (G4VLongitudinalStringDecay) ====================="<<
G4endl;
426 G4double ProbQQbar = (1.0 - 2.0*StrSup);
436 G4int Swap = stableQuarkEncoding;
437 stableQuarkEncoding = decayQuarkEncoding;
438 decayQuarkEncoding = Swap;
441 G4int IsParticle=(decayQuarkEncoding>0) ? -1 : +1;
446 G4int QuarkEncoding=QuarkPair.second->GetPDGEncoding();
450 G4int NewDecayEncoding = -1*IsParticle*(i10 * 1000 + i20 * 100 + spin);
466 created = QuarkPair.second;
482 G4double StringMT2=
string->MassT2();
483 G4double StringMT =std::sqrt(StringMT2);
490 #ifdef debug_LUNDfragmentation
492 G4cout<<
"String 4 mom, String M and Mt "<<String4Momentum<<
" "<<String4Momentum.
mag()
493 <<
" "<<std::sqrt(StringMT2)<<
G4endl;
501 #ifdef debug_LUNDfragmentation
502 G4cout<<
"Mass of the string is not sufficient to produce the hadron!"<<
G4endl;
507 String4Momentum.
setPz(0.);
512 G4double HadronMassT2, ResidualMassT2;
521 Pt2 =
sqr(HadronMt)-
sqr(HadronMass); Pt=std::sqrt(Pt2);
525 HadronPt =SampleQuarkPtw +
string->DecayPt();
527 RemSysPt = StringPt - HadronPt;
529 HadronMassT2 =
sqr(HadronMass) + HadronPt.mag2();
532 }
while (std::sqrt(HadronMassT2) + std::sqrt(ResidualMassT2) > StringMT);
537 G4double Pz2 = (
sqr(StringMT2 - HadronMassT2 - ResidualMassT2) -
538 4*HadronMassT2 * ResidualMassT2)/4./StringMT2;
540 if (Pz2 < 0 ) {
return 0;}
545 G4double zMin = (std::sqrt(HadronMassT2+Pz2) - Pz)/std::sqrt(StringMT2);
547 G4double zMax = (std::sqrt(HadronMassT2+Pz2) + Pz)/std::sqrt(StringMT2);
549 if (zMin >= zMax)
return 0;
553 HadronPt.
x(), HadronPt.
y());
559 (z *
string->LightConeDecay() -
560 HadronMassT2/(z *
string->LightConeDecay())));
561 G4double HadronE = 0.5* (z *
string->LightConeDecay() +
562 HadronMassT2/(z *
string->LightConeDecay()));
566 #ifdef debug_LUNDfragmentation
567 G4cout<<G4endl<<
" string->GetDecayDirection() "<<
string->GetDecayDirection()<<G4endl<<
G4endl;
568 G4cout<<
"string->LightConeDecay() "<<
string->LightConeDecay()<<
G4endl;
569 G4cout<<
"HadronPt,HadronE "<<HadronPt<<
" "<<HadronE<<
G4endl;
580 G4int PDGEncodingOfDecayParton,
587 G4double Mt2 = Px*Px + Py*Py + Mass*Mass;
590 G4double zOfMaxyf(0.), maxYf(1.),
z(0.), yf(1.);
592 if (!((
std::abs(PDGEncodingOfDecayParton) > 1000) && (HadronEncoding > 1000)) )
599 zOfMaxyf=BMt2/(Blund*Mt2 + 1.);}
601 zOfMaxyf = ((1.0+BMt2) - std::sqrt(
sqr(1.0-BMt2) + 4.0*BMt2*Alund))/2.0/(1.-Alund);
604 if (zOfMaxyf < zmin) {zOfMaxyf=zmin;}
605 if (zOfMaxyf > zmax) {zOfMaxyf=zmax;}
606 maxYf=(1-zOfMaxyf)/zOfMaxyf *
G4Exp(-Blund*Mt2/zOfMaxyf);
608 const G4int maxNumberOfLoops = 1000;
609 G4int loopCounter = 0;
616 while ( (
G4UniformRand()*maxYf > yf) && ++loopCounter < maxNumberOfLoops );
617 if ( loopCounter >= maxNumberOfLoops ) {
618 z = 0.5*(zmin + zmax);
623 if (
std::abs(PDGEncodingOfDecayParton) > 1000)
642 #ifdef debug_LUNDfragmentation
645 G4cout<<
"Left "<<
string->GetLeftParton()->GetPDGEncoding()<<
" "<<
string->GetPleft()<<
G4endl;
646 G4cout<<
"Right "<<
string->GetRightParton()->GetPDGEncoding()<<
" "<<
string->GetPright()<<
G4endl;
647 G4cout<<
"String4mom "<<
string->GetPstring()<<
" "<<
string->GetPstring().mag()<<
G4endl;
662 G4int sampledState = 0;
664 #ifdef debug_LUNDfragmentation
665 G4cout<<
"StrMass "<<StringMass<<
" q's "
666 <<
string->GetLeftParton()->GetParticleName()<<
" "
667 <<
string->GetRightParton()->GetParticleName()<<
G4endl;
670 string->SetLeftPartonStable();
702 #ifdef debug_LUNDfragmentation
720 #ifdef debug_LUNDfragmentation
742 #ifdef debug_LUNDfragmentation
746 G4LorentzVector P_left =
string->GetPleft(), P_right =
string->GetPright();
758 if (!(string->
DecayIsQuark() &&
string->StableIsQuark() ))
771 LeftMom *=toObserverFrame;
772 RightMom*=toObserverFrame;
774 LeftVector->push_back(
new G4KineticTrack(LeftHadron, 0, Pos, LeftMom));
775 RightVector->push_back(
new G4KineticTrack(RightHadron, 0, Pos, RightMom));
777 string->LorentzRotate(toObserverFrame);
788 G4double StringMass =
string->Mass();
791 G4int cClusterInterrupt = 0;
794 G4int LeftQuark1=
string->GetLeftParton()->GetPDGEncoding()/1000;
795 G4int LeftQuark2=(
string->GetLeftParton()->GetPDGEncoding()/100)%10;
797 G4int RightQuark1=
string->GetRightParton()->GetPDGEncoding()/1000;
798 G4int RightQuark2=(
string->GetRightParton()->GetPDGEncoding()/100)%10;
804 RightHadron= (LeftHadron ==
nullptr) ?
nullptr :
811 RightHadron=(LeftHadron ==
nullptr) ?
nullptr :
816 isOK = (LeftHadron !=
nullptr) && (RightHadron !=
nullptr);
836 G4double StringMass =
string->Mass();
843 Anti_Di_Quark =
string->GetLeftParton();
844 Di_Quark=
string->GetRightParton();
847 Anti_Di_Quark =
string->GetRightParton();
848 Di_Quark=
string->GetLeftParton();
856 G4int ADi_q1=AbsIDAnti_di_quark/1000;
857 G4int ADi_q2=(AbsIDAnti_di_quark-ADi_q1*1000)/100;
859 G4int Di_q1=AbsIDdi_quark/1000;
860 G4int Di_q2=(AbsIDdi_quark-Di_q1*1000)/100;
863 for (
G4int ProdQ=1; ProdQ < 6; ProdQ++)
866 const G4int maxNumberOfLoops = 1000;
867 G4int loopCounter = 0;
871 -
Baryon[ADi_q1-1][ADi_q2-1][ProdQ-1][StateADiQ]);
873 if (LeftHadron == NULL)
continue;
877 const G4int maxNumberOfInternalLoops = 1000;
878 G4int internalLoopCounter = 0;
882 +
Baryon[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]);
884 if (RightHadron == NULL)
continue;
887 if (StringMass > LeftHadronMass + RightHadronMass)
892 G4Exception(
"G4LundStringFragmentation::Diquark_AntiDiquark_aboveThreshold_lastSplitting ",
898 sqr(RightHadronMass));
912 }
while( (
Baryon[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]!=0) &&
913 ++internalLoopCounter < maxNumberOfInternalLoops );
914 if ( internalLoopCounter >= maxNumberOfInternalLoops ) {
919 }
while( (
Baryon[ADi_q1-1][ADi_q2-1][ProdQ-1][StateADiQ]!=0) &&
920 ++loopCounter < maxNumberOfLoops );
921 if ( loopCounter >= maxNumberOfLoops ) {
935 G4double StringMass =
string->Mass();
943 Quark =
string->GetLeftParton();
944 Di_Quark=
string->GetRightParton();
947 Quark =
string->GetRightParton();
948 Di_Quark=
string->GetLeftParton();
955 G4int Di_q1=AbsIDdi_quark/1000;
956 G4int Di_q2=(AbsIDdi_quark-Di_q1*1000)/100;
959 if (IDdi_quark < 0) SignDiQ=-1;
962 for (
G4int ProdQ=1; ProdQ < 4; ProdQ++)
967 if (IDquark == 2) SignQ= 1;
968 if ((IDquark == 1) && (ProdQ == 3)) SignQ= 1;
969 if ((IDquark == 3) && (ProdQ == 1)) SignQ=-1;
973 if (IDquark == -2) SignQ=-1;
974 if ((IDquark ==-1) && (ProdQ == 3)) SignQ=-1;
975 if ((IDquark ==-3) && (ProdQ == 1)) SignQ= 1;
978 if (AbsIDquark == ProdQ) SignQ= 1;
981 const G4int maxNumberOfLoops = 1000;
982 G4int loopCounter = 0;
986 Meson[AbsIDquark-1][ProdQ-1][StateQ]);
987 if (LeftHadron == NULL)
continue;
991 const G4int maxNumberOfInternalLoops = 1000;
992 G4int internalLoopCounter = 0;
996 Baryon[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]);
997 if (RightHadron == NULL)
continue;
1000 if (StringMass > LeftHadronMass + RightHadronMass)
1005 G4Exception(
"G4LundStringFragmentation::Quark_Diquark_lastSplitting ",
1011 sqr(RightHadronMass));
1025 }
while( (
Baryon[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]!=0) &&
1026 ++internalLoopCounter < maxNumberOfInternalLoops );
1027 if ( internalLoopCounter >= maxNumberOfInternalLoops ) {
1032 }
while( (
Meson[AbsIDquark-1][ProdQ-1][StateQ]!=0) &&
1033 ++loopCounter < maxNumberOfLoops );
1035 if ( loopCounter >= maxNumberOfLoops ) {
1049 G4double StringMass =
string->Mass();
1057 Quark =
string->GetLeftParton();
1058 Anti_Quark=
string->GetRightParton();
1061 Quark =
string->GetRightParton();
1062 Anti_Quark=
string->GetLeftParton();
1073 G4int LeftHadronCharge(0), RightHadronCharge(0);
1078 for (
G4int ProdQ=1; ProdQ < 4; ProdQ++)
1080 LeftHadronCharge = QuarkCharge -
Qcharge[ProdQ-1];
1081 G4int SignQ = LeftHadronCharge/3;
if (SignQ == 0) SignQ = 1;
1083 if ((IDquark == 1) && (ProdQ == 3)) SignQ= 1;
1084 if ((IDquark == 3) && (ProdQ == 1)) SignQ=-1;
1086 RightHadronCharge = AntiQuarkCharge + Qcharge[ProdQ-1];
1087 G4int SignAQ = RightHadronCharge/3;
if (SignAQ == 0) SignAQ = 1;
1089 if ((IDanti_quark ==-1) && (ProdQ == 3)) SignAQ=-1;
1090 if ((IDanti_quark ==-3) && (ProdQ == 1)) SignAQ= 1;
1095 const G4int maxNumberOfLoops = 1000;
1096 G4int loopCounter = 0;
1102 Meson[AbsIDquark-1][ProdQ-1][StateQ]);
1104 if (LeftHadron == NULL) { StateQ++;
continue; }
1109 const G4int maxNumberOfInternalLoops = 1000;
1110 G4int internalLoopCounter = 0;
1116 Meson[AbsIDanti_quark-1][ProdQ-1][StateAQ]);
1118 if(RightHadron == NULL) { StateAQ++;
continue; }
1122 if (StringMass > LeftHadronMass + RightHadronMass)
1127 G4Exception(
"G4LundStringFragmentation::Quark_AntiQuark_lastSplitting ",
1133 sqr(RightHadronMass));
1157 }
while ( (
Meson[AbsIDanti_quark-1][ProdQ-1][StateAQ]!=0) &&
1158 ++internalLoopCounter < maxNumberOfInternalLoops );
1159 if ( internalLoopCounter >= maxNumberOfInternalLoops ) {
1167 }
while ( (
Meson[AbsIDquark-1][ProdQ-1][StateQ]!=0) &&
1168 ++loopCounter < maxNumberOfLoops );
1169 if ( loopCounter >= maxNumberOfLoops ) {
1194 G4int indexPosition = 0;
1200 if (Sum >= ksi)
break;
1202 return indexPosition;
1214 G4double AvailablePz, AvailablePz2;
1215 #ifdef debug_LUNDfragmentation
1216 G4cout<<
"Sampling of momenta of 2 last produced hadrons ----------------"<<
G4endl;
1217 G4cout<<
"Init Mass "<<InitialMass<<
" FirstM "<<Mass<<
" SecondM "<<AntiMass<<
" ProbIsotropy "<<
G4endl;
1220 G4double r_val =
sqr(InitialMass*InitialMass - Mass*Mass - AntiMass*AntiMass) -
1221 sqr(2.*Mass*AntiMass);
1222 G4double Pabs = (r_val > 0.)? std::sqrt(r_val)/(2.*InitialMass) : 0;
1224 const G4int maxNumberOfLoops = 1000;
1226 if (Mass > 930. || AntiMass > 930.)
SigmaQT *=(1.0-0.55*
sqr((Mass+AntiMass)/InitialMass));
1228 G4int loopCounter = 0;
1232 MassMt = std::sqrt( Mass * Mass + Pt2);
1233 AntiMassMt= std::sqrt(AntiMass * AntiMass + Pt2);
1235 while ( (InitialMass < MassMt + AntiMassMt) && ++loopCounter < maxNumberOfLoops );
1237 if (Mass > 930. || AntiMass > 930.)
SigmaQT=SigmaQTw;
1239 if ( loopCounter >= maxNumberOfLoops ) {
1243 AvailablePz2=
sqr(InitialMass*InitialMass -
sqr(MassMt) -
sqr(AntiMassMt)) -
1244 4.*
sqr(MassMt*AntiMassMt);
1246 AvailablePz2 /=(4.*InitialMass*InitialMass);
1247 AvailablePz = std::sqrt(AvailablePz2);
1253 Mom->
setE(std::sqrt(
sqr(MassMt)+AvailablePz2));
1255 AntiMom->
setPx(-Px); AntiMom->
setPy(-Py); AntiMom->
setPz(-AvailablePz);
1256 AntiMom->
setE (std::sqrt(
sqr(AntiMassMt)+AvailablePz2));
1258 #ifdef debug_LUNDfragmentation