ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4LundStringFragmentation.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4LundStringFragmentation.cc
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 //
28 // -----------------------------------------------------------------------------
29 // GEANT 4 class implementation file
30 //
31 // History: first implementation, Maxim Komogorov, 10-Jul-1998
32 // -----------------------------------------------------------------------------
34 #include "G4PhysicalConstants.hh"
35 #include "G4SystemOfUnits.hh"
36 #include "Randomize.hh"
37 #include "G4FragmentingString.hh"
38 #include "G4DiQuarks.hh"
39 #include "G4Quarks.hh"
40 
41 #include "G4Exp.hh"
42 #include "G4Pow.hh"
43 
44 //#define debug_LUNDfragmentation
45 
46 // Class G4LundStringFragmentation
47 //*************************************************************************************
48 
50 {
51  SetMassCut(210.*MeV); // Mpi + Delta
52  // For ProduceOneHadron it is required
53  // that no one pi-meson can be produced.
54 
55  SigmaQT = 0.435 * GeV;
56 
59  SetStrangenessSuppression((1.0 - 0.13)/2.0);
60 
62 
63  // For the time being, set to 0.0 the probabilities for c-cbar and b-bbar creation.
64  SetProbCCbar(0.0); //(0.005); // According to O.I. Piskunova Yad. Fiz. 56 (1993) 1094
65  SetProbBBbar(0.0); //(5.0e-5); // According to O.I. Piskunova Yad. Fiz. 56 (1993) 1094
66 
67  SetMinMasses(); // For treating of small string decays
68 }
69 
70 //--------------------------------------------------------------------------------------
71 
73 {
74  // Can no longer modify Parameters for Fragmentation.
75 
76  PastInitPhase=true;
77 
78  G4FragmentingString aString(theString);
79  SetMinimalStringMass(&aString);
80 
81  #ifdef debug_LUNDfragmentation
82  G4cout<<G4endl<<"LUND StringFragmentation ------------------------------------"<<G4endl;
83  G4cout<<G4endl<<"LUND StringFragm: String Mass "
84  <<theString.Get4Momentum().mag()<<G4endl
85  <<"4Mom "<<theString.Get4Momentum()<<G4endl
86  <<"------------------------------------"<<G4endl;
87  G4cout<<"String ends Direct "<<theString.GetLeftParton()->GetPDGcode()<<" "
88  <<theString.GetRightParton()->GetPDGcode()<<" "
89  <<theString.GetDirection()<< G4endl;
90  G4cout<<"Left mom "<<theString.GetLeftParton()->Get4Momentum()<<G4endl;
91  G4cout<<"Right mom "<<theString.GetRightParton()->Get4Momentum()<<G4endl<<G4endl;
92  G4cout<<"Check for Fragmentation "<<G4endl;
93  #endif
94 
95  G4KineticTrackVector * LeftVector(0);
96 
97  if (!aString.IsAFourQuarkString() && !IsItFragmentable(&aString))
98  {
99  #ifdef debug_LUNDfragmentation
100  G4cout<<"Non fragmentable - the string is converted to one hadron "<<G4endl;
101  #endif
102  // SetMassCut(210.*MeV); // For ProduceOneHadron it is required
103  // that no one pi-meson can be produced.
104 
105  G4double Mcut = GetMassCut();
106  SetMassCut(10000.*MeV);
107  LeftVector=ProduceOneHadron(&theString);
108  SetMassCut(Mcut);
109 
110  LeftVector->operator[](0)->SetFormationTime(theString.GetTimeOfCreation());
111  LeftVector->operator[](0)->SetPosition(theString.GetPosition());
112 
113  if (LeftVector->size() > 1)
114  {
115  // 2 hadrons created from qq-qqbar are stored
116  LeftVector->operator[](1)->SetFormationTime(theString.GetTimeOfCreation());
117  LeftVector->operator[](1)->SetPosition(theString.GetPosition());
118  }
119  return LeftVector;
120  }
121 
122  #ifdef debug_LUNDfragmentation
123  G4cout<<"The string will be fragmented. "<<G4endl;
124  #endif
125 
126  // The string can fragment. At least two particles can be produced.
127  LeftVector =new G4KineticTrackVector;
128  G4KineticTrackVector * RightVector=new G4KineticTrackVector;
129 
130  G4bool success = Loop_toFragmentString(theString, LeftVector, RightVector);
131 
132  if ( ! success )
133  {
134  std::for_each(LeftVector->begin(), LeftVector->end(), DeleteKineticTrack());
135  LeftVector->clear();
136  std::for_each(RightVector->begin(), RightVector->end(), DeleteKineticTrack());
137  delete RightVector;
138  return LeftVector;
139  }
140 
141  // Join Left- and RightVector into LeftVector in correct order.
142  while (!RightVector->empty())
143  {
144  LeftVector->push_back(RightVector->back());
145  RightVector->erase(RightVector->end()-1);
146  }
147  delete RightVector;
148 
149  return LeftVector;
150 }
151 
152 //----------------------------------------------------------------------------------
153 
155 {
156  SetMinimalStringMass(string);
157  //G4cout<<"MinM StrM "<<MinimalStringMass<<" "<< string->Get4Momentum().mag()<<G4endl;
158 
159  return std::abs(MinimalStringMass) < string->Get4Momentum().mag();
160 
161  //MinimalStringMass is negative and large for a string with unknown particles in a final 2-particle decay.
162 }
163 
164 //----------------------------------------------------------------------------------------
165 
167  G4KineticTrackVector * & LeftVector,
168  G4KineticTrackVector * & RightVector )
169 {
170  #ifdef debug_LUNDfragmentation
171  G4cout<<"Loop_toFrag "<<theString.GetLeftParton()->GetPDGcode()<<" "
172  <<theString.GetLeftParton()->Get4Momentum()<<G4endl
173  <<" "<<theString.GetRightParton()->GetPDGcode()<<" "
174  <<theString.GetRightParton()->Get4Momentum()<<G4endl
175  <<"Direction "<<theString.GetDirection()<< G4endl;
176  #endif
177 
178  G4bool final_success=false;
179  G4bool inner_success=true;
180 
181  G4int attempt=0;
182 
183  while ( ! final_success && attempt++ < StringLoopInterrupt )
184  { // If the string fragmentation does not be happend,
185  // repeat the fragmentation.
186 
187  G4FragmentingString *currentString=new G4FragmentingString(theString);
188  G4LorentzRotation toCms, toObserverFrame;
189 
190  //G4cout<<"Main loop start whilecounter "<<attempt<<G4endl;
191 
192  // Cleaning up the previously produced hadrons
193  std::for_each(LeftVector->begin(), LeftVector->end(), DeleteKineticTrack());
194  LeftVector->clear();
195  std::for_each(RightVector->begin(), RightVector->end(), DeleteKineticTrack());
196  RightVector->clear();
197 
198  // Main fragmentation loop until the string will not be able to fragment
199  inner_success=true; // set false on failure.
200  const G4int maxNumberOfLoops = 1000;
201  G4int loopCounter = -1;
202 
203  while ( (! StopFragmenting(currentString)) && ++loopCounter < maxNumberOfLoops )
204  { // Split current string into hadron + new string
205  #ifdef debug_LUNDfragmentation
206  G4cout<<"The string will fragment. "<<G4endl;;
207  //G4cout<<"1 "<<currentString->GetDecayDirection()<<G4endl;
208  #endif
209  G4FragmentingString *newString=0; // used as output from SplitUp.
210 
211  toCms=currentString->TransformToAlignedCms();
212  toObserverFrame= toCms.inverse();
213  #ifdef debug_LUNDfragmentation
214  //G4cout<<"CMS Left mom "<<currentString->GetPleft()<<G4endl;
215  //G4cout<<"CMS Right mom "<<currentString->GetPright()<<G4endl;
216  //G4cout<<"CMS String M "<<currentString->GetPstring()<<G4endl;
217  #endif
218 
219  G4KineticTrack * Hadron=Splitup(currentString,newString);
220 
221  if ( Hadron != 0 ) // Store the hadron
222  {
223  #ifdef debug_LUNDfragmentation
224  G4cout<<"Hadron prod at fragm. "<<Hadron->GetDefinition()->GetParticleName()<<G4endl;
225  //G4cout<<"2 "<<currentString->GetDecayDirection()<<G4endl;
226  #endif
227 
228  Hadron->Set4Momentum(toObserverFrame*Hadron->Get4Momentum());
229 
230  G4double TimeOftheStringCreation=theString.GetTimeOfCreation();
231  G4ThreeVector PositionOftheStringCreation(theString.GetPosition());
232 
233  G4LorentzVector Coordinate(Hadron->GetPosition(), Hadron->GetFormationTime());
234  G4LorentzVector Momentum = toObserverFrame*Coordinate;
235  Hadron->SetFormationTime(TimeOftheStringCreation + Momentum.e() - fermi/c_light);
236  G4ThreeVector aPosition(Momentum.vect());
237  Hadron->SetPosition(PositionOftheStringCreation+aPosition);
238 
239  // Open to protect hadron production at fragmentation
240  if ( currentString->GetDecayDirection() > 0 )
241  {
242  LeftVector->push_back(Hadron);
243  } else
244  {
245  RightVector->push_back(Hadron);
246  }
247  delete currentString;
248  currentString=newString;
249  } else {
250  if ( newString ) delete newString;
251  }
252 
253  currentString->LorentzRotate(toObserverFrame);
254  };
255 
256  if ( loopCounter >= maxNumberOfLoops ) {
257  inner_success=false;
258  }
259 
260  // Split remaining string into 2 final hadrons.
261  #ifdef debug_LUNDfragmentation
262  if (inner_success) G4cout<<"Split remaining string into 2 final hadrons."<<G4endl;
263  #endif
264  if ( inner_success && SplitLast(currentString, LeftVector, RightVector) ) // Close to protect Last Str. Decay
265  {
266  final_success = true;
267  }
268 
269  delete currentString;
270  } // End of the loop where we try to fragment the string.
271  return final_success;
272 }
273 
274 //----------------------------------------------------------------------------------------
275 
277 {
278  SetMinimalStringMass(string);
279 
280  if ( MinimalStringMass < 0.) return true;
281 
282  if (string->IsAFourQuarkString())
283  {
284  return G4UniformRand() < G4Exp(-0.0005*(string->Mass() - MinimalStringMass));
285  } else {
286 
287  if (MinimalStringMass < 0.0 ) return false; // For a string with di-quark having c or b quarks and s, c, b quarks
288 
289  G4bool Result = G4UniformRand() <
290  G4Exp(-0.66e-6*(string->Mass()*string->Mass() - MinimalStringMass*MinimalStringMass));
291  // G4bool Result = string->Mass() < MinimalStringMass + 150.*MeV*G4UniformRand(); // a'la LUND
292 
293  #ifdef debug_LUNDfragmentation
294  G4cout<<"StopFragmenting MinimalStringMass string->Mass() "<<MinimalStringMass
295  <<" "<<string->Mass()<<G4endl;
296  G4cout<<"StopFragmenting - Yes/No "<<Result<<G4endl;
297  #endif
298  return Result;
299  }
300 }
301 
302 //-----------------------------------------------------------------------------
303 
305  G4FragmentingString *&newString)
306 {
307  #ifdef debug_LUNDfragmentation
308  G4cout<<G4endl;
309  G4cout<<"Start SplitUP ========================="<<G4endl;
310  G4cout<<"String partons: " <<string->GetLeftParton()->GetPDGEncoding()<<" "
311  <<string->GetRightParton()->GetPDGEncoding()<<" "
312  <<"Direction " <<string->GetDecayDirection()<<G4endl;
313  #endif
314 
315  //... random choice of string end to use for creating the hadron (decay)
316  G4int SideOfDecay = (G4UniformRand() < 0.5)? 1: -1;
317  if (SideOfDecay < 0)
318  {
319  string->SetLeftPartonStable();
320  } else
321  {
322  string->SetRightPartonStable();
323  }
324 
325  G4ParticleDefinition *newStringEnd;
326  G4ParticleDefinition * HadronDefinition;
327 
328  G4double StringMass=string->Mass();
329 
330  G4double ProbDqADq = GetDiquarkSuppress();
331  G4double ProbSaS = 1.0 - 2.0 * GetStrangeSuppress();
332 
333  #ifdef debug_LUNDfragmentation
334  G4cout<<"StrMass DiquarkSuppression "<<StringMass<<" "<<GetDiquarkSuppress()<<G4endl;
335  #endif
336 
337  G4int NumberOfpossibleBaryons = 2;
338 
339  if (string->GetLeftParton()->GetParticleSubType() != "quark") NumberOfpossibleBaryons++;
340  if (string->GetRightParton()->GetParticleSubType() != "quark") NumberOfpossibleBaryons++;
341 
342  G4double ActualProb = ProbDqADq ;
343  ActualProb *= (1.0-sqr(NumberOfpossibleBaryons*1400.0/StringMass));
344 
345  SetDiquarkSuppression(ActualProb);
346 
347  G4double Mth = 1250.0; // 2 Mk + Mpi
348  if ( NumberOfpossibleBaryons == 3 ){Mth = 2520.0;} // Mlambda/Msigma + Mk + Mpi
349  else if ( NumberOfpossibleBaryons == 4 ){Mth = 2380.0;} // 2 Mlambda/Msigma + Mk + Mpi
350  else {}
351 
352  ActualProb = ProbSaS * (1.0 - G4Pow::GetInstance()->powA( Mth/StringMass , 4.0 ));
353  if ( ActualProb < 0.0 ) ActualProb = 0.0;
354  SetStrangenessSuppression((1.0-ActualProb)/2.0);
355 
356  #ifdef debug_LUNDfragmentation
357  G4cout<<"StrMass DiquarkSuppression corrected "<<StringMass<<" "<<GetDiquarkSuppress()<<G4endl;
358  #endif
359 
360  if (string->DecayIsQuark())
361  {
362  HadronDefinition= QuarkSplitup(string->GetDecayParton(), newStringEnd);
363  } else {
364  HadronDefinition= DiQuarkSplitup(string->GetDecayParton(), newStringEnd);
365  }
366 
367  SetDiquarkSuppression(ProbDqADq);
368  SetStrangenessSuppression((1.0-ProbSaS)/2.0);
369 
370  if ( HadronDefinition == NULL ) { G4KineticTrack * Hadron =0; return Hadron; }
371 
372  #ifdef debug_LUNDfragmentation
373  G4cout<<"The parton "<<string->GetDecayParton()->GetPDGEncoding()<<" "
374  <<" produces hadron "<<HadronDefinition->GetParticleName()
375  <<" and is transformed to "<<newStringEnd->GetPDGEncoding()<<G4endl;
376  G4cout<<"The side of the string decay Left/Right (1/-1) "<<SideOfDecay<<G4endl;
377  #endif
378  // create new String from old, ie. keep Left and Right order, but replace decay
379 
380  if ( newString ) delete newString;
381 
382  newString=new G4FragmentingString(*string,newStringEnd); // To store possible quark containt of new string
383 
384  #ifdef debug_LUNDfragmentation
385  G4cout<<"An attempt to determine its energy (SplitEandP)"<<G4endl;
386  #endif
387  G4LorentzVector* HadronMomentum=SplitEandP(HadronDefinition, string, newString);
388 
389  delete newString; newString=0;
390 
391  G4KineticTrack * Hadron =0;
392  if ( HadronMomentum != 0 ) {
393 
394  #ifdef debug_LUNDfragmentation
395  G4cout<<"The attempt was successful"<<G4endl;
396  #endif
398  Hadron = new G4KineticTrack(HadronDefinition, 0,Pos, *HadronMomentum);
399 
400  if ( newString ) delete newString;
401 
402  newString=new G4FragmentingString(*string,newStringEnd,
403  HadronMomentum);
404  delete HadronMomentum;
405  }
406  else
407  {
408  #ifdef debug_LUNDfragmentation
409  G4cout<<"The attempt was not successful !!!"<<G4endl;
410  #endif
411  }
412 
413  #ifdef debug_LUNDfragmentation
414  G4cout<<"End SplitUP (G4VLongitudinalStringDecay) ====================="<<G4endl;
415  #endif
416 
417  return Hadron;
418 }
419 
420 //-----------------------------------------------------------------------------
421 
423  G4ParticleDefinition *&created)
424 {
425  G4double StrSup=GetStrangeSuppress();
426  G4double ProbQQbar = (1.0 - 2.0*StrSup);
427 
428  //... can Diquark break or not?
430 
431  //... Diquark break
432  G4int stableQuarkEncoding = decay->GetPDGEncoding()/1000;
433  G4int decayQuarkEncoding = (decay->GetPDGEncoding()/100)%10;
434  if (G4UniformRand() < 0.5)
435  {
436  G4int Swap = stableQuarkEncoding;
437  stableQuarkEncoding = decayQuarkEncoding;
438  decayQuarkEncoding = Swap;
439  }
440 
441  G4int IsParticle=(decayQuarkEncoding>0) ? -1 : +1; // if we have a quark, we need antiquark)
442 
443  pDefPair QuarkPair = CreatePartonPair(IsParticle,false); // no diquarks wanted
444 
445  //... Build new Diquark
446  G4int QuarkEncoding=QuarkPair.second->GetPDGEncoding();
447  G4int i10 = std::max(std::abs(QuarkEncoding), std::abs(stableQuarkEncoding));
448  G4int i20 = std::min(std::abs(QuarkEncoding), std::abs(stableQuarkEncoding));
449  G4int spin = (i10 != i20 && G4UniformRand() <= 0.5)? 1 : 3;
450  G4int NewDecayEncoding = -1*IsParticle*(i10 * 1000 + i20 * 100 + spin);
451  created = FindParticle(NewDecayEncoding);
452  G4ParticleDefinition * decayQuark=FindParticle(decayQuarkEncoding);
453  G4ParticleDefinition * had=hadronizer->Build(QuarkPair.first, decayQuark);
454  StrangeSuppress=StrSup;
455 
456  return had;
457 
458  } else {
459  //... Diquark does not break
460 
461  G4int IsParticle=(decay->GetPDGEncoding()>0) ? +1 : -1; // if we have a diquark, we need quark
462 
463  StrangeSuppress=(1.0 - ProbQQbar * 0.9)/2.0;
464  pDefPair QuarkPair = CreatePartonPair(IsParticle,false); // no diquarks wanted
465 
466  created = QuarkPair.second;
467 
468  G4ParticleDefinition * had=hadronizer->Build(QuarkPair.first, decay);
469  StrangeSuppress=StrSup;
470 
471  return had;
472  }
473 }
474 
475 //-----------------------------------------------------------------------------
476 
478  G4FragmentingString * string,
479  G4FragmentingString * newString)
480 {
481  G4LorentzVector String4Momentum=string->Get4Momentum();
482  G4double StringMT2=string->MassT2();
483  G4double StringMT =std::sqrt(StringMT2);
484 
485  G4double HadronMass = pHadron->GetPDGMass();
486  SetMinimalStringMass(newString);
487 
488  if ( MinimalStringMass < 0.0 ) return nullptr;
489 
490  #ifdef debug_LUNDfragmentation
491  G4cout<<G4endl<<"Start LUND SplitEandP "<<G4endl;
492  G4cout<<"String 4 mom, String M and Mt "<<String4Momentum<<" "<<String4Momentum.mag()
493  <<" "<<std::sqrt(StringMT2)<<G4endl;
494  G4cout<<"Hadron "<<pHadron->GetParticleName()<<G4endl;
495  G4cout<<"HadM MinimalStringMassLeft StringM hM+sM "<<HadronMass<<" "<<MinimalStringMass<<" "
496  <<String4Momentum.mag()<<" "<<HadronMass+MinimalStringMass<<G4endl;
497  #endif
498 
499  if ((HadronMass + MinimalStringMass > string->Mass()) || MinimalStringMass < 0.)
500  {
501  #ifdef debug_LUNDfragmentation
502  G4cout<<"Mass of the string is not sufficient to produce the hadron!"<<G4endl;
503  #endif
504  return 0;
505  } // have to start all over!
506 
507  String4Momentum.setPz(0.);
508  G4ThreeVector StringPt=String4Momentum.vect();
509 
510  // calculate and assign hadron transverse momentum component HadronPx and HadronPy
511  G4ThreeVector HadronPt , RemSysPt;
512  G4double HadronMassT2, ResidualMassT2;
513  G4double HadronMt, Pt, Pt2, phi;
514  //... sample Pt of the hadron
515  G4int attempt=0;
516  do
517  {
518  attempt++; if (attempt > StringLoopInterrupt) {return 0;}
519 
520  HadronMt = HadronMass - 300.0*G4Log(G4UniformRand());
521  Pt2 = sqr(HadronMt)-sqr(HadronMass); Pt=std::sqrt(Pt2);
522  phi = 2.*pi*G4UniformRand();
523  G4ThreeVector SampleQuarkPtw= G4ThreeVector(Pt * std::cos(phi),Pt * std::sin(phi),0);
524 
525  HadronPt =SampleQuarkPtw + string->DecayPt();
526  HadronPt.setZ(0);
527  RemSysPt = StringPt - HadronPt;
528 
529  HadronMassT2 = sqr(HadronMass) + HadronPt.mag2();
530  ResidualMassT2=sqr(MinimalStringMass) + RemSysPt.mag2();
531 
532  } while (std::sqrt(HadronMassT2) + std::sqrt(ResidualMassT2) > StringMT);
533 
534  //... sample z to define hadron longitudinal momentum and energy
535  //... but first check the available phase space
536 
537  G4double Pz2 = (sqr(StringMT2 - HadronMassT2 - ResidualMassT2) -
538  4*HadronMassT2 * ResidualMassT2)/4./StringMT2;
539 
540  if (Pz2 < 0 ) {return 0;} // have to start all over!
541 
542  //... then compute allowed z region z_min <= z <= z_max
543 
544  G4double Pz = std::sqrt(Pz2);
545  G4double zMin = (std::sqrt(HadronMassT2+Pz2) - Pz)/std::sqrt(StringMT2);
546  // G4double zMin = (std::sqrt(HadronMassT2+Pz2) - 0.)/std::sqrt(StringMT2); // For testing purposes
547  G4double zMax = (std::sqrt(HadronMassT2+Pz2) + Pz)/std::sqrt(StringMT2);
548 
549  if (zMin >= zMax) return 0; // have to start all over!
550 
551  G4double z = GetLightConeZ(zMin, zMax,
552  string->GetDecayParton()->GetPDGEncoding(), pHadron,
553  HadronPt.x(), HadronPt.y());
554 
555  //... now compute hadron longitudinal momentum and energy
556  // longitudinal hadron momentum component HadronPz
557 
558  HadronPt.setZ(0.5* string->GetDecayDirection() *
559  (z * string->LightConeDecay() -
560  HadronMassT2/(z * string->LightConeDecay())));
561  G4double HadronE = 0.5* (z * string->LightConeDecay() +
562  HadronMassT2/(z * string->LightConeDecay()));
563 
564  G4LorentzVector * a4Momentum= new G4LorentzVector(HadronPt,HadronE);
565 
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;
570  G4cout<<"String4Momentum "<<String4Momentum<<G4endl;
571  G4cout<<"Out of LUND SplitEandP "<<G4endl<<G4endl;
572  #endif
573 
574  return a4Momentum;
575 }
576 
577 //-----------------------------------------------------------------------------------------
578 
580  G4int PDGEncodingOfDecayParton,
581  G4ParticleDefinition* pHadron,
582  G4double Px, G4double Py)
583 {
584  G4double Mass = pHadron->GetPDGMass();
585  G4int HadronEncoding=std::abs(pHadron->GetPDGEncoding());
586 
587  G4double Mt2 = Px*Px + Py*Py + Mass*Mass;
588 
589  G4double Alund, Blund;
590  G4double zOfMaxyf(0.), maxYf(1.), z(0.), yf(1.);
591 
592  if (!((std::abs(PDGEncodingOfDecayParton) > 1000) && (HadronEncoding > 1000)) )
593  { // ---------------- Quark fragmentation and qq-> meson ----------------------
594  Alund=1.;
595  Blund=0.7/GeV/GeV;
596 
597  G4double BMt2 = Blund*Mt2;
598  if (Alund == 1.0) {
599  zOfMaxyf=BMt2/(Blund*Mt2 + 1.);}
600  else {
601  zOfMaxyf = ((1.0+BMt2) - std::sqrt(sqr(1.0-BMt2) + 4.0*BMt2*Alund))/2.0/(1.-Alund);
602  }
603 
604  if (zOfMaxyf < zmin) {zOfMaxyf=zmin;}
605  if (zOfMaxyf > zmax) {zOfMaxyf=zmax;}
606  maxYf=(1-zOfMaxyf)/zOfMaxyf * G4Exp(-Blund*Mt2/zOfMaxyf);
607 
608  const G4int maxNumberOfLoops = 1000;
609  G4int loopCounter = 0;
610  do
611  {
612  z = zmin + G4UniformRand()*(zmax-zmin);
613  //yf = (1-z)/z * G4Exp(-Blund*Mt2/z);
614  yf = G4Pow::GetInstance()->powA(1.0-z,Alund)/z*G4Exp(-BMt2/z);
615  }
616  while ( (G4UniformRand()*maxYf > yf) && ++loopCounter < maxNumberOfLoops );
617  if ( loopCounter >= maxNumberOfLoops ) {
618  z = 0.5*(zmin + zmax); // Just a value between zmin and zmax, no physics considerations at all!
619  }
620  return z;
621  }
622 
623  if (std::abs(PDGEncodingOfDecayParton) > 1000)
624  {
625  G4double an = 2.5;
626  an +=(sqr(Px)+sqr(Py))/sqr(GeV)-0.5;
627  z=zmin + (zmax-zmin)*G4Pow::GetInstance()->powA(G4UniformRand(),1./an);
628  }
629 
630  return z;
631 }
632 
633 //----------------------------------------------------------------------------------------------------------
634 
636  G4KineticTrackVector * LeftVector,
637  G4KineticTrackVector * RightVector)
638 {
639  //... perform last cluster decay
640  SetMinimalStringMass( string);
641  if ( MinimalStringMass < 0.) return false;
642  #ifdef debug_LUNDfragmentation
643  G4cout<<G4endl<<"Split last-----------------------------------------"<<G4endl;
644  G4cout<<"MinimalStringMass "<<MinimalStringMass<<G4endl;
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;
648  #endif
649 
650  G4LorentzVector Str4Mom=string->Get4Momentum();
651 
652  G4LorentzRotation toCms=string->TransformToAlignedCms();
653  G4LorentzRotation toObserverFrame= toCms.inverse();
654 
655  G4double StringMass=string->Mass();
656 
657  G4ParticleDefinition * LeftHadron(0), * RightHadron(0);
658 
659  NumberOf_FS=0;
660  for (G4int i=0; i<350; i++) {FS_Weight[i]=0.;}
661 
662  G4int sampledState = 0;
663 
664  #ifdef debug_LUNDfragmentation
665  G4cout<<"StrMass "<<StringMass<<" q's "
666  <<string->GetLeftParton()->GetParticleName()<<" "
667  <<string->GetRightParton()->GetParticleName()<<G4endl;
668  #endif
669 
670  string->SetLeftPartonStable(); // to query quark contents..
671 
672  if (string->IsAFourQuarkString() )
673  {
674  // The string is qq-qqbar type. Diquarks are on the string ends
675  if (StringMass-MinimalStringMass < 0.)
676  {
677  if (! Diquark_AntiDiquark_belowThreshold_lastSplitting(string, LeftHadron, RightHadron) )
678  {
679  return false;
680  }
681  } else
682  {
683  Diquark_AntiDiquark_aboveThreshold_lastSplitting(string, LeftHadron, RightHadron);
684 
685  if (NumberOf_FS == 0) return false;
686 
687  sampledState = SampleState();
688  if (string->GetLeftParton()->GetPDGEncoding() < 0)
689  {
690  LeftHadron =FS_LeftHadron[sampledState];
691  RightHadron=FS_RightHadron[sampledState];
692  } else
693  {
694  LeftHadron =FS_RightHadron[sampledState];
695  RightHadron=FS_LeftHadron[sampledState];
696  }
697  }
698  } else
699  {
700  if (string->DecayIsQuark() && string->StableIsQuark() )
701  { //... there are quarks on cluster ends
702  #ifdef debug_LUNDfragmentation
703  G4cout<<"Q Q string LastSplit"<<G4endl;
704  #endif
705  Quark_AntiQuark_lastSplitting(string, LeftHadron, RightHadron);
706 
707  if (NumberOf_FS == 0) return false;
708  sampledState = SampleState();
709  if (string->GetLeftParton()->GetPDGEncoding() < 0)
710  {
711  LeftHadron =FS_RightHadron[sampledState];
712  RightHadron=FS_LeftHadron[sampledState];
713  } else
714  {
715  LeftHadron =FS_LeftHadron[sampledState];
716  RightHadron=FS_RightHadron[sampledState];
717  }
718  } else
719  { //... there is a Diquark on one of the cluster ends
720  #ifdef debug_LUNDfragmentation
721  G4cout<<"DiQ Q string Last Split"<<G4endl;
722  #endif
723 
724  Quark_Diquark_lastSplitting(string, LeftHadron, RightHadron);
725 
726  if (NumberOf_FS == 0) return false;
727  sampledState = SampleState();
728 
729  if (string->GetLeftParton()->GetParticleSubType() == "quark")
730  {
731  LeftHadron =FS_LeftHadron[sampledState];
732  RightHadron=FS_RightHadron[sampledState];
733  } else
734  {
735  LeftHadron =FS_RightHadron[sampledState];
736  RightHadron=FS_LeftHadron[sampledState];
737  }
738  }
739 
740  }
741 
742  #ifdef debug_LUNDfragmentation
743  G4cout<<"Sampled hadrons: "<<LeftHadron->GetParticleName()<<" "<<RightHadron->GetParticleName()<<G4endl;
744  #endif
745 
746  G4LorentzVector P_left =string->GetPleft(), P_right = string->GetPright();
747 
748  G4LorentzVector LeftMom, RightMom;
750 
751  Sample4Momentum(&LeftMom, LeftHadron->GetPDGMass(),
752  &RightMom, RightHadron->GetPDGMass(),
753  StringMass);
754 
755  // Sample4Momentum ascribes LeftMom.pz() along positive Z axis for baryons in many cases.
756  // It must be negative in case when the system is moving against Z axis.
757 
758  if (!(string->DecayIsQuark() && string->StableIsQuark() ))
759  { // Only for qq - q, q - qq, and qq - qqbar -------------------
760 
761  if (std::abs(string->GetLeftParton()->GetPDGEncoding()) > 1000)
762  {
763  if (P_left.z() <= 0.) {G4LorentzVector tmp = LeftMom; LeftMom=RightMom; RightMom=tmp;}
764  }
765  else
766  {
767  if (P_right.z() >= 0.) {G4LorentzVector tmp = LeftMom; LeftMom=RightMom; RightMom=tmp;}
768  }
769  }
770 
771  LeftMom *=toObserverFrame;
772  RightMom*=toObserverFrame;
773 
774  LeftVector->push_back(new G4KineticTrack(LeftHadron, 0, Pos, LeftMom));
775  RightVector->push_back(new G4KineticTrack(RightHadron, 0, Pos, RightMom));
776 
777  string->LorentzRotate(toObserverFrame);
778  return true;
779 }
780 
781 //----------------------------------------------------------------------------------------
782 
785  G4ParticleDefinition * & LeftHadron,
786  G4ParticleDefinition * & RightHadron)
787 {
788  G4double StringMass = string->Mass();
789 
790  G4bool isOK = false;
791  G4int cClusterInterrupt = 0;
792  do
793  {
794  G4int LeftQuark1= string->GetLeftParton()->GetPDGEncoding()/1000;
795  G4int LeftQuark2=(string->GetLeftParton()->GetPDGEncoding()/100)%10;
796 
797  G4int RightQuark1= string->GetRightParton()->GetPDGEncoding()/1000;
798  G4int RightQuark2=(string->GetRightParton()->GetPDGEncoding()/100)%10;
799 
800  if (G4UniformRand()<0.5)
801  {
802  LeftHadron =hadronizer->Build(FindParticle( LeftQuark1),
803  FindParticle(RightQuark1));
804  RightHadron= (LeftHadron == nullptr) ? nullptr :
805  hadronizer->Build(FindParticle( LeftQuark2),
806  FindParticle(RightQuark2));
807  } else
808  {
809  LeftHadron =hadronizer->Build(FindParticle( LeftQuark1),
810  FindParticle(RightQuark2));
811  RightHadron=(LeftHadron == nullptr) ? nullptr :
812  hadronizer->Build(FindParticle( LeftQuark2),
813  FindParticle(RightQuark1));
814  }
815 
816  isOK = (LeftHadron != nullptr) && (RightHadron != nullptr);
817  if(isOK) { isOK = (StringMass > LeftHadron->GetPDGMass() + RightHadron->GetPDGMass()); }
818  ++cClusterInterrupt;
819  //... repeat procedure, if mass of cluster is too low to produce hadrons
820  //... ClusterMassCut = 0.15*GeV model parameter
821  }
822  while (isOK == false || cClusterInterrupt < ClusterLoopInterrupt);
823  /* Loop checking, 07.08.2015, A.Ribon */
824  return isOK;
825 }
826 
827 //----------------------------------------------------------------------------------------
828 
831  G4ParticleDefinition * & LeftHadron,
832  G4ParticleDefinition * & RightHadron)
833 {
834  // StringMass-MinimalStringMass > 0. Creation of 2 baryons is possible ----
835 
836  G4double StringMass = string->Mass();
837  G4double StringMassSqr= sqr(StringMass);
838  G4ParticleDefinition * Di_Quark;
839  G4ParticleDefinition * Anti_Di_Quark;
840 
841  if (string->GetLeftParton()->GetPDGEncoding() < 0)
842  {
843  Anti_Di_Quark =string->GetLeftParton();
844  Di_Quark=string->GetRightParton();
845  } else
846  {
847  Anti_Di_Quark =string->GetRightParton();
848  Di_Quark=string->GetLeftParton();
849  }
850 
851  G4int IDAnti_di_quark =Anti_Di_Quark->GetPDGEncoding();
852  G4int AbsIDAnti_di_quark =std::abs(IDAnti_di_quark);
853  G4int IDdi_quark =Di_Quark->GetPDGEncoding();
854  G4int AbsIDdi_quark =std::abs(IDdi_quark);
855 
856  G4int ADi_q1=AbsIDAnti_di_quark/1000;
857  G4int ADi_q2=(AbsIDAnti_di_quark-ADi_q1*1000)/100;
858 
859  G4int Di_q1=AbsIDdi_quark/1000;
860  G4int Di_q2=(AbsIDdi_quark-Di_q1*1000)/100;
861 
862  NumberOf_FS=0;
863  for (G4int ProdQ=1; ProdQ < 6; ProdQ++)
864  {
865  G4int StateADiQ=0;
866  const G4int maxNumberOfLoops = 1000;
867  G4int loopCounter = 0;
868  do // while(Meson[AbsIDquark-1][ProdQ-1][StateQ]<>0);
869  {
871  -Baryon[ADi_q1-1][ADi_q2-1][ProdQ-1][StateADiQ]);
872 
873  if (LeftHadron == NULL) continue;
874  G4double LeftHadronMass=LeftHadron->GetPDGMass();
875 
876  G4int StateDiQ=0;
877  const G4int maxNumberOfInternalLoops = 1000;
878  G4int internalLoopCounter = 0;
879  do // while(Baryon[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]<>0);
880  {
882  +Baryon[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]);
883 
884  if (RightHadron == NULL) continue;
885  G4double RightHadronMass=RightHadron->GetPDGMass();
886 
887  if (StringMass > LeftHadronMass + RightHadronMass)
888  {
889  if ( NumberOf_FS > 349 ) {
891  ed << " NumberOf_FS exceeds its limit: NumberOf_FS=" << NumberOf_FS << G4endl;
892  G4Exception( "G4LundStringFragmentation::Diquark_AntiDiquark_aboveThreshold_lastSplitting ",
893  "HAD_LUND_001", JustWarning, ed );
894  NumberOf_FS = 349;
895  }
896 
897  G4double FS_Psqr=lambda(StringMassSqr,sqr(LeftHadronMass),
898  sqr(RightHadronMass));
899  //FS_Psqr=1.;
900  FS_Weight[NumberOf_FS]=std::sqrt(FS_Psqr)*FS_Psqr*
901  BaryonWeight[ADi_q1-1][ADi_q2-1][ProdQ-1][StateADiQ]*
902  BaryonWeight[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]*
903  Prob_QQbar[ProdQ-1];
904 
905  FS_LeftHadron[NumberOf_FS] = LeftHadron;
906  FS_RightHadron[NumberOf_FS]= RightHadron;
907  NumberOf_FS++;
908  } // End of if (StringMass > LeftHadronMass + RightHadronMass)
909 
910  StateDiQ++;
911 
912  } while( (Baryon[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]!=0) &&
913  ++internalLoopCounter < maxNumberOfInternalLoops );
914  if ( internalLoopCounter >= maxNumberOfInternalLoops ) {
915  return false;
916  }
917 
918  StateADiQ++;
919  } while( (Baryon[ADi_q1-1][ADi_q2-1][ProdQ-1][StateADiQ]!=0) &&
920  ++loopCounter < maxNumberOfLoops );
921  if ( loopCounter >= maxNumberOfLoops ) {
922  return false;
923  }
924  } // End of for (G4int ProdQ=1; ProdQ < 4; ProdQ++)
925 
926  return true;
927 }
928 
929 //----------------------------------------------------------------------------------------
930 
932  G4ParticleDefinition * & LeftHadron,
933  G4ParticleDefinition * & RightHadron)
934 {
935  G4double StringMass = string->Mass();
936  G4double StringMassSqr= sqr(StringMass);
937 
938  G4ParticleDefinition * Di_Quark;
939  G4ParticleDefinition * Quark;
940 
941  if (string->GetLeftParton()->GetParticleSubType()== "quark")
942  {
943  Quark =string->GetLeftParton();
944  Di_Quark=string->GetRightParton();
945  } else
946  {
947  Quark =string->GetRightParton();
948  Di_Quark=string->GetLeftParton();
949  }
950 
951  G4int IDquark =Quark->GetPDGEncoding();
952  G4int AbsIDquark =std::abs(IDquark);
953  G4int IDdi_quark =Di_Quark->GetPDGEncoding();
954  G4int AbsIDdi_quark=std::abs(IDdi_quark);
955  G4int Di_q1=AbsIDdi_quark/1000;
956  G4int Di_q2=(AbsIDdi_quark-Di_q1*1000)/100;
957 
958  G4int SignDiQ= 1;
959  if (IDdi_quark < 0) SignDiQ=-1;
960 
961  NumberOf_FS=0;
962  for (G4int ProdQ=1; ProdQ < 4; ProdQ++)
963  {
964  G4int SignQ;
965  if (IDquark > 0)
966  { SignQ=-1;
967  if (IDquark == 2) SignQ= 1;
968  if ((IDquark == 1) && (ProdQ == 3)) SignQ= 1; // K0
969  if ((IDquark == 3) && (ProdQ == 1)) SignQ=-1; // K0bar
970  } else
971  {
972  SignQ= 1;
973  if (IDquark == -2) SignQ=-1;
974  if ((IDquark ==-1) && (ProdQ == 3)) SignQ=-1; // K0bar
975  if ((IDquark ==-3) && (ProdQ == 1)) SignQ= 1; // K0
976  }
977 
978  if (AbsIDquark == ProdQ) SignQ= 1;
979 
980  G4int StateQ=0;
981  const G4int maxNumberOfLoops = 1000;
982  G4int loopCounter = 0;
983  do // while(Meson[AbsIDquark-1][ProdQ-1][StateQ]<>0);
984  {
986  Meson[AbsIDquark-1][ProdQ-1][StateQ]);
987  if (LeftHadron == NULL) continue;
988  G4double LeftHadronMass=LeftHadron->GetPDGMass();
989 
990  G4int StateDiQ=0;
991  const G4int maxNumberOfInternalLoops = 1000;
992  G4int internalLoopCounter = 0;
993  do // while(Baryon[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]<>0);
994  {
995  RightHadron=G4ParticleTable::GetParticleTable()->FindParticle(SignDiQ*
996  Baryon[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]);
997  if (RightHadron == NULL) continue;
998  G4double RightHadronMass=RightHadron->GetPDGMass();
999 
1000  if (StringMass > LeftHadronMass + RightHadronMass)
1001  {
1002  if ( NumberOf_FS > 349 ) {
1004  ed << " NumberOf_FS exceeds its limit: NumberOf_FS=" << NumberOf_FS << G4endl;
1005  G4Exception( "G4LundStringFragmentation::Quark_Diquark_lastSplitting ",
1006  "HAD_LUND_002", JustWarning, ed );
1007  NumberOf_FS = 349;
1008  }
1009 
1010  G4double FS_Psqr=lambda(StringMassSqr,sqr(LeftHadronMass),
1011  sqr(RightHadronMass));
1012  FS_Weight[NumberOf_FS]=std::sqrt(FS_Psqr)*
1013  MesonWeight[AbsIDquark-1][ProdQ-1][StateQ]*
1014  BaryonWeight[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]*
1015  Prob_QQbar[ProdQ-1];
1016 
1017  FS_LeftHadron[NumberOf_FS] = LeftHadron;
1018  FS_RightHadron[NumberOf_FS]= RightHadron;
1019 
1020  NumberOf_FS++;
1021  } // End of if (StringMass > LeftHadronMass + RightHadronMass)
1022 
1023  StateDiQ++;
1024 
1025  } while( (Baryon[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]!=0) &&
1026  ++internalLoopCounter < maxNumberOfInternalLoops );
1027  if ( internalLoopCounter >= maxNumberOfInternalLoops ) {
1028  return false;
1029  }
1030 
1031  StateQ++;
1032  } while( (Meson[AbsIDquark-1][ProdQ-1][StateQ]!=0) &&
1033  ++loopCounter < maxNumberOfLoops ); /* Loop checking, 07.08.2015, A.Ribon */
1034 
1035  if ( loopCounter >= maxNumberOfLoops ) {
1036  return false;
1037  }
1038  }
1039 
1040  return true;
1041 }
1042 
1043 //----------------------------------------------------------------------------------------
1044 
1046  G4ParticleDefinition * & LeftHadron,
1047  G4ParticleDefinition * & RightHadron)
1048 {
1049  G4double StringMass = string->Mass();
1050  G4double StringMassSqr= sqr(StringMass);
1051 
1052  G4ParticleDefinition * Quark;
1053  G4ParticleDefinition * Anti_Quark;
1054 
1055  if (string->GetLeftParton()->GetPDGEncoding()>0)
1056  {
1057  Quark =string->GetLeftParton();
1058  Anti_Quark=string->GetRightParton();
1059  } else
1060  {
1061  Quark =string->GetRightParton();
1062  Anti_Quark=string->GetLeftParton();
1063  }
1064 
1065  G4int IDquark =Quark->GetPDGEncoding();
1066  G4int AbsIDquark =std::abs(IDquark);
1067  G4int QuarkCharge =Qcharge[IDquark-1];
1068 
1069  G4int IDanti_quark =Anti_Quark->GetPDGEncoding();
1070  G4int AbsIDanti_quark =std::abs(IDanti_quark);
1071  G4int AntiQuarkCharge =-Qcharge[AbsIDanti_quark-1];
1072 
1073  G4int LeftHadronCharge(0), RightHadronCharge(0);
1074 
1075  //G4cout<<"Q Qbar "<<IDquark<<" "<<IDanti_quark<<G4endl;
1076 
1077  NumberOf_FS=0;
1078  for (G4int ProdQ=1; ProdQ < 4; ProdQ++)
1079  {
1080  LeftHadronCharge = QuarkCharge - Qcharge[ProdQ-1];
1081  G4int SignQ = LeftHadronCharge/3; if (SignQ == 0) SignQ = 1;
1082 
1083  if ((IDquark == 1) && (ProdQ == 3)) SignQ= 1; // K0
1084  if ((IDquark == 3) && (ProdQ == 1)) SignQ=-1; // K0bar
1085 
1086  RightHadronCharge = AntiQuarkCharge + Qcharge[ProdQ-1];
1087  G4int SignAQ = RightHadronCharge/3; if (SignAQ == 0) SignAQ = 1;
1088 
1089  if ((IDanti_quark ==-1) && (ProdQ == 3)) SignAQ=-1; // K0bar
1090  if ((IDanti_quark ==-3) && (ProdQ == 1)) SignAQ= 1; // K0
1091 
1092  //G4cout<<"ProQ signs "<<ProdQ<<" "<<SignQ<<" "<<SignAQ<<G4endl;
1093 
1094  G4int StateQ=0;
1095  const G4int maxNumberOfLoops = 1000;
1096  G4int loopCounter = 0;
1097  do
1098  {
1099  //G4cout<<"[AbsIDquark-1][ProdQ-1][StateQ "<<AbsIDquark-1<<" "
1100  //<<ProdQ-1<<" "<<StateQ<<" "<<SignQ*Meson[AbsIDquark-1][ProdQ-1][StateQ]<<G4endl;
1101  LeftHadron=G4ParticleTable::GetParticleTable()->FindParticle(SignQ*
1102  Meson[AbsIDquark-1][ProdQ-1][StateQ]);
1103  //G4cout<<"LeftHadron "<<LeftHadron<<G4endl;
1104  if (LeftHadron == NULL) { StateQ++; continue; }
1105  //G4cout<<"LeftHadron "<<LeftHadron->GetParticleName()<<G4endl;
1106  G4double LeftHadronMass=LeftHadron->GetPDGMass();
1107 
1108  G4int StateAQ=0;
1109  const G4int maxNumberOfInternalLoops = 1000;
1110  G4int internalLoopCounter = 0;
1111  do
1112  {
1113  //G4cout<<" [AbsIDanti_quark-1][ProdQ-1][StateAQ] "<<AbsIDanti_quark-1<<" "
1114  // <<ProdQ-1<<" "<<StateAQ<<" "<<SignAQ*Meson[AbsIDanti_quark-1][ProdQ-1][StateAQ]<<G4endl;
1115  RightHadron=G4ParticleTable::GetParticleTable()->FindParticle(SignAQ*
1116  Meson[AbsIDanti_quark-1][ProdQ-1][StateAQ]);
1117  //G4cout<<"RightHadron "<<RightHadron<<G4endl;
1118  if(RightHadron == NULL) { StateAQ++; continue; }
1119  //G4cout<<"RightHadron "<<RightHadron->GetParticleName()<<G4endl;
1120  G4double RightHadronMass=RightHadron->GetPDGMass();
1121 
1122  if (StringMass > LeftHadronMass + RightHadronMass)
1123  {
1124  if ( NumberOf_FS > 349 ) {
1126  ed << " NumberOf_FS exceeds its limit: NumberOf_FS=" << NumberOf_FS << G4endl;
1127  G4Exception( "G4LundStringFragmentation::Quark_AntiQuark_lastSplitting ",
1128  "HAD_LUND_003", JustWarning, ed );
1129  NumberOf_FS = 349;
1130  }
1131 
1132  G4double FS_Psqr=lambda(StringMassSqr,sqr(LeftHadronMass),
1133  sqr(RightHadronMass));
1134  //FS_Psqr=1.;
1135  FS_Weight[NumberOf_FS]=std::sqrt(FS_Psqr)*
1136  MesonWeight[AbsIDquark-1][ProdQ-1][StateQ]*
1137  MesonWeight[AbsIDanti_quark-1][ProdQ-1][StateAQ]*
1138  Prob_QQbar[ProdQ-1];
1139 
1140  if (string->GetLeftParton()->GetPDGEncoding()>0)
1141  {
1142  FS_LeftHadron[NumberOf_FS] = RightHadron;
1143  FS_RightHadron[NumberOf_FS]= LeftHadron;
1144  } else
1145  {
1146  FS_LeftHadron[NumberOf_FS] = LeftHadron;
1147  FS_RightHadron[NumberOf_FS]= RightHadron;
1148  }
1149 
1150  NumberOf_FS++;
1151 
1152  }
1153 
1154  StateAQ++;
1155  //G4cout<<" StateAQ Meson[AbsIDanti_quark-1][ProdQ-1][StateAQ] "<<StateAQ<<" "
1156  // <<Meson[AbsIDanti_quark-1][ProdQ-1][StateAQ]<<" "<<internalLoopCounter<<G4endl;
1157  } while ( (Meson[AbsIDanti_quark-1][ProdQ-1][StateAQ]!=0) &&
1158  ++internalLoopCounter < maxNumberOfInternalLoops );
1159  if ( internalLoopCounter >= maxNumberOfInternalLoops ) {
1160  return false;
1161  }
1162 
1163  StateQ++;
1164  //G4cout<<"StateQ Meson[AbsIDquark-1][ProdQ-1][StateQ] "<<StateQ<<" "
1165  // <<Meson[AbsIDquark-1][ProdQ-1][StateQ]<<" "<<loopCounter<<G4endl;
1166 
1167  } while ( (Meson[AbsIDquark-1][ProdQ-1][StateQ]!=0) &&
1168  ++loopCounter < maxNumberOfLoops );
1169  if ( loopCounter >= maxNumberOfLoops ) {
1170  return false;
1171  }
1172  } // End of for (G4int ProdQ=1; ProdQ < 4; ProdQ++)
1173 
1174  return true;
1175 }
1176 
1177 //----------------------------------------------------------------------------------------------------------
1178 
1180 {
1181  if ( NumberOf_FS > 349 ) {
1183  ed << " NumberOf_FS exceeds its limit: NumberOf_FS=" << NumberOf_FS << G4endl;
1184  G4Exception( "G4LundStringFragmentation::SampleState ", "HAD_LUND_004", JustWarning, ed );
1185  NumberOf_FS = 349;
1186  }
1187 
1188  G4double SumWeights=0.;
1189 
1190  for (G4int i=0; i<NumberOf_FS; i++) {SumWeights+=FS_Weight[i];}
1191 
1192  G4double ksi=G4UniformRand();
1193  G4double Sum=0.;
1194  G4int indexPosition = 0;
1195 
1196  for (G4int i=0; i<NumberOf_FS; i++)
1197  {
1198  Sum+=(FS_Weight[i]/SumWeights);
1199  indexPosition=i;
1200  if (Sum >= ksi) break;
1201  }
1202  return indexPosition;
1203 }
1204 
1205 //----------------------------------------------------------------------------------------------------------
1206 
1208  G4LorentzVector* AntiMom, G4double AntiMass,
1209  G4double InitialMass)
1210 {
1211  // ------ Sampling of momenta of 2 last produced hadrons --------------------
1212  G4ThreeVector Pt;
1213  G4double MassMt, AntiMassMt;
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;
1218  #endif
1219 
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;
1223 
1224  const G4int maxNumberOfLoops = 1000;
1225  G4double SigmaQTw=SigmaQT;
1226  if (Mass > 930. || AntiMass > 930.) SigmaQT *=(1.0-0.55*sqr((Mass+AntiMass)/InitialMass));
1227 
1228  G4int loopCounter = 0;
1229  do
1230  {
1231  Pt=SampleQuarkPt(Pabs); Pt.setZ(0); G4double Pt2=Pt.mag2();
1232  MassMt = std::sqrt( Mass * Mass + Pt2);
1233  AntiMassMt= std::sqrt(AntiMass * AntiMass + Pt2);
1234  }
1235  while ( (InitialMass < MassMt + AntiMassMt) && ++loopCounter < maxNumberOfLoops );
1236 
1237  if (Mass > 930. || AntiMass > 930.) SigmaQT=SigmaQTw;
1238 
1239  if ( loopCounter >= maxNumberOfLoops ) {
1240  AvailablePz2 = 0.0;
1241  }
1242 
1243  AvailablePz2= sqr(InitialMass*InitialMass - sqr(MassMt) - sqr(AntiMassMt)) -
1244  4.*sqr(MassMt*AntiMassMt);
1245 
1246  AvailablePz2 /=(4.*InitialMass*InitialMass);
1247  AvailablePz = std::sqrt(AvailablePz2);
1248 
1249  G4double Px=Pt.getX();
1250  G4double Py=Pt.getY();
1251 
1252  Mom->setPx(Px); Mom->setPy(Py); Mom->setPz(AvailablePz);
1253  Mom->setE(std::sqrt(sqr(MassMt)+AvailablePz2));
1254 
1255  AntiMom->setPx(-Px); AntiMom->setPy(-Py); AntiMom->setPz(-AvailablePz);
1256  AntiMom->setE (std::sqrt(sqr(AntiMassMt)+AvailablePz2));
1257 
1258  #ifdef debug_LUNDfragmentation
1259  G4cout<<"Fmass Mom "<<Mom->getX()<<" "<<Mom->getY()<<" "<<Mom->getZ()<<" "<<Mom->getT()<<G4endl;
1260  G4cout<<"Smass Mom "<<AntiMom->getX()<<" "<<AntiMom->getY()<<" "<<AntiMom->getZ()
1261  <<" "<<AntiMom->getT()<<G4endl;
1262  #endif
1263 }
1264 
1265 //------------------------------------------------------------------------
1266 
1268 {
1269  G4double lam = sqr(S - m1_Sqr - m2_Sqr) - 4.*m1_Sqr*m2_Sqr;
1270  return lam;
1271 }
1272 
1273 // --------------------------------------------------------------
1275 {}
1276