50 arho(0.5), aphi(0.),
an(-0.5), ala(-0.75), aksi(-1.), alft(0.5)
81 #ifdef debug_QGSMfragmentation
85 <<
"------------------------------------"<<
G4endl;
101 #ifdef debug_QGSMfragmentation
102 if ( LeftVector != 0 )
G4cout<<
"Non fragmentable - the string is converted to one hadron "<<
G4endl;
105 if ( LeftVector != 0 )
return LeftVector;
107 #ifdef debug_QGSMfragmentation
118 G4bool success=
false, inner_sucess=
true;
122 #ifdef debug_QGSMfragmentation
133 RightVector->clear();
136 const G4int maxNumberOfLoops = 1000;
137 G4int loopCounter = -1;
138 while (!
StopFragmenting(currentString) && ++loopCounter < maxNumberOfLoops )
141 #ifdef debug_QGSMfragmentation
149 #ifdef debug_QGSMfragmentation
153 LeftVector->push_back(Hadron);
155 RightVector->push_back(Hadron);
157 delete currentString;
158 currentString=newString;
162 #ifdef debug_QGSMfragmentation
163 G4cout<<
"abandon ... start from the beginning ---------------"<<
G4endl;
167 if (newString)
delete newString;
172 if ( loopCounter >= maxNumberOfLoops ) {
177 #ifdef debug_QGSMfragmentation
178 G4cout<<
"Split remaining string into 2 final hadrons."<<
G4endl;
182 SplitLast(currentString,LeftVector, RightVector) )
186 delete currentString;
189 delete theStringInCMS;
201 while(!RightVector->empty())
203 LeftVector->push_back(RightVector->back());
204 RightVector->erase(RightVector->end()-1);
212 for(
size_t C1 = 0;
C1 < LeftVector->size();
C1++)
216 Momentum = toObserverFrame*Momentum;
219 Momentum = toObserverFrame*Coordinate;
235 #ifdef debug_QGSMfragmentation
240 G4double invD1(0.),invD2(0.), r1(0.),r2(0.),r12(0.);
242 G4int absCode = std::abs( PartonEncoding );
246 q1 = absHadronCode/1000; q2 = (absHadronCode % 1000)/100; q3 = (absHadronCode % 100)/10;
248 G4bool StrangeHadron = (q1 == 3) || (q2 == 3) || (q3 == 3);
252 if(absCode == 1 || absCode == 2)
254 if(absHadronCode < 1000)
256 if( !StrangeHadron ) {d1=2.0; d2 = -
arho +
alft;}
260 if( !StrangeHadron ) {d1=0.0; d2 =
arho - 2.0*
an +
alft;}
264 else if(absCode == 3)
266 if(absHadronCode < 1000){d1=1.0 -
aphi; d2 = -
arho +
alft;}
269 }
else throw G4HadronicException(__FILE__, __LINE__,
"Unknown PDGencoding in G4QGSMFragmentation::G4LightConeZ");
271 #ifdef debug_QGSMfragmentation
277 invD1=1./
d1; invD2=1./
d2;
279 const G4int maxNumberOfLoops = 10000;
280 G4int loopCounter = 0;
287 }
while( ( (r12 > 1.0) || !((zmin <= z)&&(z <= zmax))) && ++loopCounter < maxNumberOfLoops );
288 if ( loopCounter >= maxNumberOfLoops ) {
289 z = 0.5*(zmin + zmax);
296 if(absCode == 1103 || absCode == 2101 ||
297 absCode == 2203 || absCode == 2103)
299 if(absHadronCode < 1000)
301 if( !StrangeHadron ) {d1=1.0; d2=
arho - 2.0*
an +
alft;}
309 #ifdef debug_QGSMfragmentation
314 invD1=1./
d1; invD2=1./
d2;
316 const G4int maxNumberOfLoops = 10000;
317 G4int loopCounter = 0;
324 }
while( ( (r12 > 1.0) || !((zmin <= z)&&(z <= zmax))) && ++loopCounter < maxNumberOfLoops );
325 if ( loopCounter >= maxNumberOfLoops ) {
326 z = 0.5*(zmin + zmax);
331 else if(absCode == 3101 || absCode == 3103 ||
332 absCode == 3201 || absCode == 3203)
342 const G4int maxNumberOfLoops = 1000;
343 G4int loopCounter = 0;
350 while( (
G4UniformRand() > yf) && ++loopCounter < maxNumberOfLoops );
368 #ifdef debug_QGSMfragmentation
370 G4cout<<
"String 4 mom, String M "<<
string->Get4Momentum()<<
" "<<
string->Mass()<<
G4endl;
371 G4cout<<
"HadM MinimalStringMassLeft StringM hM+sM "<<HadronMass<<
" "<<MinimalStringMass<<
" "
372 <<
string->Mass()<<
" "<<HadronMass+MinimalStringMass<<
G4endl;
375 if(HadronMass + MinimalStringMass > string->
Mass())
377 #ifdef debug_QGSMfragmentation
378 G4cout<<
"Mass of the string is not sufficient to produce the hadron!"<<
G4endl;
384 G4double StringMT2 =
string->MassT2();
385 G4double StringMT = std::sqrt(StringMT2);
388 String4Momentum.
setPz(0.);
392 G4double HadronMassT2, ResidualMassT2;
402 RemSysPt = StringPt - HadronPt;
404 HadronMassT2 =
sqr(HadronMass) + HadronPt.mag2();
405 ResidualMassT2=
sqr(MinimalStringMass) + RemSysPt.
mag2();
407 }
while(std::sqrt(HadronMassT2) + std::sqrt(ResidualMassT2) > StringMT);
412 G4double Pz2 = (
sqr(StringMT2 - HadronMassT2 - ResidualMassT2) -
413 4*HadronMassT2 * ResidualMassT2)/4./StringMT2;
415 if(Pz2 < 0 ) {
return 0;}
420 G4double zMin = (std::sqrt(HadronMassT2+Pz2) - Pz)/std::sqrt(StringMT2);
421 G4double zMax = (std::sqrt(HadronMassT2+Pz2) + Pz)/std::sqrt(StringMT2);
423 if (zMin >= zMax)
return 0;
427 HadronPt.
x(), HadronPt.
y());
433 (z *
string->LightConeDecay() -
434 HadronMassT2/(z *
string->LightConeDecay())));
435 G4double HadronE = 0.5* (z *
string->LightConeDecay() +
436 HadronMassT2/(z *
string->LightConeDecay()));
440 #ifdef debug_QGSMfragmentation
441 G4cout<<
"string->GetDecayDirection() string->LightConeDecay() "
442 <<
string->GetDecayDirection()<<
" "<<
string->LightConeDecay()<<
G4endl;
443 G4cout<<
"HadronPt,HadronE "<<HadronPt<<
" "<<HadronE<<
G4endl;
459 G4ThreeVector ClusterVel =
string->Get4Momentum().boostVector();
460 G4double ResidualMass =
string->Mass();
462 #ifdef debug_QGSMfragmentation
463 G4cout<<
"Split last-----------------------------------------"<<
G4endl;
464 G4cout<<
"StrMass "<<ResidualMass<<
" q's "
465 <<
string->GetLeftParton()->GetParticleName()<<
" "
466 <<
string->GetRightParton()->GetParticleName()<<
G4endl;
470 G4int cClusterInterrupt = 0;
472 const G4int maxNumberOfLoops = 1000;
473 G4int loopCounter = 0;
482 string->SetLeftPartonStable();
492 IsParticle=(
string->GetLeftParton()->GetPDGEncoding()>0) ? -1 : +1;
494 IsParticle=(
string->GetLeftParton()->GetPDGEncoding()>0) ? +1 : -1;
503 quark = QuarkPair.second;
514 while ( (ResidualMass <= LeftHadron->GetPDGMass() + RightHadron->
GetPDGMass() + ClusterMassCut)
515 && ++loopCounter < maxNumberOfLoops );
517 if ( loopCounter >= maxNumberOfLoops ) {
525 &RightMom, RightHadron->
GetPDGMass(), ResidualMass);
526 LeftMom.
boost(ClusterVel);
527 RightMom.
boost(ClusterVel);
529 #ifdef debug_QGSMfragmentation
535 LeftVector->push_back(
new G4KineticTrack(LeftHadron, 0, Pos, LeftMom));
536 RightVector->push_back(
new G4KineticTrack(RightHadron, 0, Pos, RightMom));
556 string->Get4Momentum().mag2();
566 G4double r_val =
sqr(InitialMass*InitialMass - Mass*Mass - AntiMass*AntiMass) -
sqr(2.*Mass*AntiMass);
567 G4double Pabs = (r_val > 0.)? std::sqrt(r_val)/(2.*InitialMass) : 0;
571 G4double st = std::sqrt(1. - pz * pz)*Pabs;
578 Mom->
setE(std::sqrt(Pabs*Pabs + Mass*Mass));
581 AntiMom->
setE (std::sqrt(Pabs*Pabs + AntiMass*AntiMass));
599 G4int Swap = stableQuarkEncoding;
600 stableQuarkEncoding = decayQuarkEncoding;
601 decayQuarkEncoding = Swap;
604 G4int IsParticle=(decayQuarkEncoding>0) ? -1 : +1;
613 G4int QuarkEncoding=QuarkPair.second->GetPDGEncoding();
614 G4int i10 =
std::max(std::abs(QuarkEncoding), std::abs(stableQuarkEncoding));
615 G4int i20 =
std::min(std::abs(QuarkEncoding), std::abs(stableQuarkEncoding));
617 G4int NewDecayEncoding = -1*IsParticle*(i10 * 1000 + i20 * 100 + spin);
635 created = QuarkPair.second;
649 #ifdef debug_QGSMfragmentation
651 G4cout<<
"Start SplitUP (G4VLongitudinalStringDecay) ========================="<<
G4endl;
652 G4cout<<
"String partons: " <<
string->GetLeftParton()->GetPDGEncoding()<<
" "
653 <<
string->GetRightParton()->GetPDGEncoding()<<
" "
654 <<
"Direction " <<
string->GetDecayDirection()<<
G4endl;
661 string->SetLeftPartonStable();
664 string->SetRightPartonStable();
674 G4int NumberOfpossibleBaryons = 2;
680 ActualProb *= (1.0-
G4Exp(2.0*(1.0 - string->
Mass()/(NumberOfpossibleBaryons*1400.0))));
691 #ifdef debug_QGSMfragmentation
692 G4cout<<
"The parton "<<
string->GetDecayParton()->GetPDGEncoding()<<
" "
695 G4cout<<
"The side of the string decay Left/Right (1/-1) "<<SideOfDecay<<
G4endl;
702 #ifdef debug_QGSMfragmentation
703 G4cout<<
"An attempt to determine its energy (SplitEandP)"<<
G4endl;
707 delete newString; newString=0;
710 if ( HadronMomentum != 0 ) {
712 #ifdef debug_QGSMfragmentation
716 Hadron =
new G4KineticTrack(HadronDefinition, 0,Pos, *HadronMomentum);
721 delete HadronMomentum;
726 #ifdef debug_QGSMfragmentation
731 #ifdef debug_VStringDecay
732 G4cout<<
"End SplitUP (G4VLongitudinalStringDecay) ====================="<<
G4endl;
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
ParticleList decay(Cluster *const c)
Carries out a cluster decay.
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4Parton * GetRightParton(void) const
G4ParticleDefinition * Build(G4ParticleDefinition *black, G4ParticleDefinition *white)
virtual G4KineticTrackVector * FragmentString(const G4ExcitedString &theString)
const G4ThreeVector & GetPosition() const
const G4LorentzVector & Get4Momentum() const
void SetDiquarkSuppression(G4double aValue)
G4ParticleDefinition * GetRightParton(void) const
const G4String & GetParticleSubType() const
virtual G4LorentzVector * SplitEandP(G4ParticleDefinition *pHadron, G4FragmentingString *string, G4FragmentingString *newString)
G4int StringLoopInterrupt
const G4String & GetParticleName() const
const G4ParticleDefinition * GetDefinition() const
void SetVectorMesonMixings(std::vector< G4double > aVector)
G4double GetDiquarkSuppress()
G4double FragmentationMass(const G4FragmentingString *const string, Pcreate build=0, pDefPair *pdefs=0)
G4double GetPDGMass() const
void SetSpinThreeHalfBarionProbability(G4double aValue)
void SetDiquarkBreakProbability(G4double aValue)
std::vector< G4double > vectorMesonMix
G4int GetDirection(void) const
void SetVectorMesonProbability(G4double aValue)
static G4Pow * GetInstance()
G4ParticleDefinition * FindParticle(G4int Encoding)
G4int GetDecayDirection() const
G4ThreeVector SampleQuarkPt(G4double ptMax=-1.)
void CalculateHadronTimePosition(G4double theInitialStringMass, G4KineticTrackVector *)
G4double powA(G4double A, G4double y) const
G4KineticTrackVector * LightFragmentationTest(const G4ExcitedString *const theString)
void SetStrangenessSuppression(G4double aValue)
G4LorentzVector Get4Momentum() const
HepLorentzRotation inverse() const
G4Parton * GetLeftParton(void) const
virtual G4bool SplitLast(G4FragmentingString *string, G4KineticTrackVector *LeftVector, G4KineticTrackVector *RightVector)
G4ParticleDefinition * QuarkSplitup(G4ParticleDefinition *decay, G4ParticleDefinition *&created)
virtual void Sample4Momentum(G4LorentzVector *Mom, G4double Mass, G4LorentzVector *AntiMom, G4double AntiMass, G4double InitialMass)
G4ParticleDefinition * BuildHighSpin(G4ParticleDefinition *black, G4ParticleDefinition *white)
G4double GetFormationTime() const
G4ParticleDefinition * GetLeftParton(void) const
G4double GetStrangeSuppress()
void SetPosition(const G4ThreeVector aPosition)
virtual G4bool IsFragmentable(const G4FragmentingString *const string)
pDefPair CreatePartonPair(G4int NeedParticle, G4bool AllowDiquarks=true)
virtual G4KineticTrack * Splitup(G4FragmentingString *string, G4FragmentingString *&newString)
G4ParticleDefinition * GetDecayParton() const
G4int GetPDGEncoding() const
G4GLOB_DLL std::ostream G4cout
virtual G4ParticleDefinition * DiQuarkSplitup(G4ParticleDefinition *decay, G4ParticleDefinition *&created)
G4int ClusterLoopInterrupt
virtual G4double GetLightConeZ(G4double zmin, G4double zmax, G4int PartonEncoding, G4ParticleDefinition *pHadron, G4double Px, G4double Py)
G4ExcitedString * CopyExcited(const G4ExcitedString &string)
static constexpr double pi
const G4ThreeVector & GetPosition() const
CLHEP::HepLorentzVector G4LorentzVector
G4HadronBuilder * hadronizer
std::pair< G4ParticleDefinition *, G4ParticleDefinition * > pDefPair
const G4LorentzVector & Get4Momentum() const
void Set4Momentum(const G4LorentzVector &a4Momentum)
void SetFormationTime(G4double aFormationTime)
G4LorentzRotation TransformToAlignedCms()
virtual G4bool StopFragmenting(const G4FragmentingString *const string)
G4double DiquarkBreakProb
HepLorentzVector & boost(double, double, double)
T min(const T t1, const T t2)
brief Return the smallest of the two arguments