Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
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 // $Id: G4LundStringFragmentation.cc 108072 2017-12-19 15:31:21Z gcosmo $
28 // GEANT4 tag $Name: $ 1.8
29 //
30 // -----------------------------------------------------------------------------
31 // GEANT 4 class implementation file
32 //
33 // History: first implementation, Maxim Komogorov, 10-Jul-1998
34 // -----------------------------------------------------------------------------
35 #include "G4LundStringFragmentation.hh"
36 #include "G4PhysicalConstants.hh"
37 #include "G4SystemOfUnits.hh"
38 #include "Randomize.hh"
39 #include "G4FragmentingString.hh"
40 #include "G4DiQuarks.hh"
41 #include "G4Quarks.hh"
42 
43 #include "G4Exp.hh"
44 #include "G4Pow.hh"
45 
46 //#define debug_LUNDfragmentation
47 
48 // Class G4LundStringFragmentation
49 //*************************************************************************************
50 
51 G4LundStringFragmentation::G4LundStringFragmentation()
52 {
53  SetMassCut(160.*MeV); // For LightFragmentationTest it is required
54  // that no one pi-meson can be produced.
55 
56  SigmaQT = 0.435 * GeV;
57 
58  SetStringTensionParameter(1.*GeV/fermi);
59  SetDiquarkBreakProbability(0.5); //(0.05);
60  SetStrangenessSuppression((1.0 - 0.10)/2.0);// 0.446 -> 0.424 Uzhi Nov. 2017 0.12
61 
62  SetDiquarkSuppression(0.150); // 0.067 -> 0.15 Uzhi Nov. 2017
63 
64 
65 // For treating of small string decays
66  SetMinMasses();
67 
68  //... pspin_meson is probability to create pseudo-scalar meson
69  pspin_meson = 0.55; SetVectorMesonProbability(pspin_meson); // Uzhi Dec. 2017 0.5
70 
71  //... pspin_barion is probability to create 1/2 barion
72  pspin_barion = 0.5; SetSpinThreeHalfBarionProbability(pspin_barion); // Uzhi Dec. 2017
73 
74  //... vectorMesonMix[] is quark mixing parameters for vector mesons (Variable spin = 3)
75  vectorMesonMix[0] = 0.; // 0.5 Uzhi Dec. 2017
76  vectorMesonMix[1] = 0.45;
77  vectorMesonMix[2] = 0.0;
78  vectorMesonMix[3] = 0.45; // 0.5 Uzhi Dec. 2017
79  vectorMesonMix[4] = 1.0;
80  vectorMesonMix[5] = 1.0;
81 
82  SetVectorMesonMixings(vectorMesonMix); // Uzhi Dec. 2017
83 
84 }
85 
86 //--------------------------------------------------------------------------------------
87 G4KineticTrackVector* G4LundStringFragmentation::FragmentString(const G4ExcitedString& theString)
88 {
89  // Can no longer modify Parameters for Fragmentation.
90 
91  PastInitPhase=true;
92 
93  G4FragmentingString aString(theString);
94  SetMinimalStringMass(&aString);
95 
96 #ifdef debug_LUNDfragmentation
97  G4cout<<G4endl<<"LUND StringFragmentation ------------------------------------"<<G4endl;
98  G4cout<<G4endl<<"LUND StringFragm: String Mass "
99  <<theString.Get4Momentum().mag()<<G4endl
100  <<" 4Mom "<<theString.Get4Momentum()
101  <<"------------------------------------"<<G4endl;
102  G4cout<<"String ends Direct "<<theString.GetLeftParton()->GetPDGcode()<<" "
103  <<theString.GetRightParton()->GetPDGcode()<<" "
104  <<theString.GetDirection()<< G4endl;
105  G4cout<<"Left mom "<<theString.GetLeftParton()->Get4Momentum()<<G4endl;
106  G4cout<<"Right mom "<<theString.GetRightParton()->Get4Momentum()<<G4endl;
107  G4cout<<"Check for Fragmentation "<<G4endl;
108 #endif
109 
110  G4KineticTrackVector * LeftVector(0);
111 
112  if(!aString.FourQuarkString() && !IsFragmentable(&aString))
113  {
114 #ifdef debug_LUNDfragmentation
115  G4cout<<"Non fragmentable - the string is converted to one hadron "<<G4endl;
116 #endif
117 // SetMassCut(160.*MeV); // For LightFragmentationTest it is required
118 // // that no one pi-meson can be produced.
119 
120  G4double Mcut = GetMassCut();
121  SetMassCut(10000.*MeV);
122  LeftVector=LightFragmentationTest(&theString);
123  SetMassCut(Mcut);
124 
125  LeftVector->operator[](0)->SetFormationTime(theString.GetTimeOfCreation());
126  LeftVector->operator[](0)->SetPosition(theString.GetPosition());
127 
128  if(LeftVector->size() > 1)
129  {
130  // 2 hadrons created from qq-qqbar are stored
131  LeftVector->operator[](1)->SetFormationTime(theString.GetTimeOfCreation());
132  LeftVector->operator[](1)->SetPosition(theString.GetPosition());
133  }
134  return LeftVector;
135  } // end of if(!IsFragmentable(&aString))
136 
137 #ifdef debug_LUNDfragmentation
138  G4cout<<"The string will be fragmented. "<<G4endl;
139 #endif
140 
141  // The string can fragment. At least two particles can be produced.
142  LeftVector =new G4KineticTrackVector;
143  G4KineticTrackVector * RightVector=new G4KineticTrackVector;
144 
145  G4bool success = Loop_toFragmentString(theString, LeftVector, RightVector);
146 
147  if ( ! success )
148  {
149  std::for_each(LeftVector->begin(), LeftVector->end(), DeleteKineticTrack());
150  LeftVector->clear();
151  std::for_each(RightVector->begin(), RightVector->end(), DeleteKineticTrack());
152  delete RightVector;
153  return LeftVector;
154  }
155 
156  // Join Left- and RightVector into LeftVector in correct order.
157  while(!RightVector->empty())
158  {
159  LeftVector->push_back(RightVector->back());
160  RightVector->erase(RightVector->end()-1);
161  }
162  delete RightVector;
163 
164  return LeftVector;
165 }
166 
167 //----------------------------------------------------------------------------------
168 G4bool G4LundStringFragmentation::IsFragmentable(const G4FragmentingString * const string)
169 {
170  SetMinimalStringMass(string);
171 // G4cout<<"MinM StrM "<<MinimalStringMass<<" "<< string->Get4Momentum().mag()<<G4endl;
172  return MinimalStringMass < string->Get4Momentum().mag();
173 }
174 
175 //----------------------------------------------------------------------------------------
176 G4bool G4LundStringFragmentation::Loop_toFragmentString(
177  const G4ExcitedString &theString,
178  G4KineticTrackVector * & LeftVector,
179  G4KineticTrackVector * & RightVector )
180 {
181 #ifdef debug_LUNDfragmentation
182  G4cout<<"Loop_toFrag "<<theString.GetLeftParton()->GetPDGcode()<<" "
183  <<theString.GetLeftParton()->Get4Momentum()<<G4endl
184  <<" "<<theString.GetRightParton()->GetPDGcode()<<" "
185  <<theString.GetRightParton()->Get4Momentum()<<G4endl
186  <<"Direction "<<theString.GetDirection()<< G4endl;
187 #endif
188 
189  G4bool final_success=false;
190  G4bool inner_success=true;
191 
192  G4int attempt=0;
193 
194  while ( ! final_success && attempt++ < StringLoopInterrupt )
195  { // If the string fragmentation does not be happend,
196  // repeat the fragmentation.
197 
198  G4FragmentingString *currentString=new G4FragmentingString(theString);
199  G4LorentzRotation toCms, toObserverFrame;
200 
201  //G4cout<<"Main loop start whilecounter "<<attempt<<G4endl;
202 
203  // Cleaning up the previously produced hadrons
204  std::for_each(LeftVector->begin(), LeftVector->end(), DeleteKineticTrack());
205  LeftVector->clear();
206  std::for_each(RightVector->begin(), RightVector->end(), DeleteKineticTrack());
207  RightVector->clear();
208 
209  // Main fragmentation loop until the string will not be able to fragment
210  inner_success=true; // set false on failure.
211  const G4int maxNumberOfLoops = 1000;
212  G4int loopCounter = -1;
213 
214  while ( (! StopFragmenting(currentString)) && ++loopCounter < maxNumberOfLoops )
215  { // Split current string into hadron + new string
216 #ifdef debug_LUNDfragmentation
217  G4cout<<"The string can fragment. "<<G4endl;;
218 //G4cout<<"1 "<<currentString->GetDecayDirection()<<G4endl;
219 #endif
220  G4FragmentingString *newString=0; // used as output from SplitUp.
221 
222  toCms=currentString->TransformToAlignedCms();
223  toObserverFrame= toCms.inverse();
224 #ifdef debug_LUNDfragmentation
225  G4cout<<"CMS Left mom "<<currentString->GetPleft()<<G4endl;
226  G4cout<<"CMS Right mom "<<currentString->GetPright()<<G4endl;
227  G4cout<<"CMS String M "<<currentString->GetPstring()<<G4endl;
228 #endif
229 
230  G4KineticTrack * Hadron=Splitup(currentString,newString);
231 
232  if ( Hadron != 0 ) // Store the hadron
233  {
234 #ifdef debug_LUNDfragmentation
235  G4cout<<"Hadron prod at fragm. "<<Hadron->GetDefinition()->GetParticleName()<<G4endl;
236 //G4cout<<"2 "<<currentString->GetDecayDirection()<<G4endl;
237 #endif
238 
239  Hadron->Set4Momentum(toObserverFrame*Hadron->Get4Momentum());
240 
241  G4double TimeOftheStringCreation=theString.GetTimeOfCreation();
242  G4ThreeVector PositionOftheStringCreation(theString.GetPosition());
243 
244  G4LorentzVector Coordinate(Hadron->GetPosition(), Hadron->GetFormationTime());
245  G4LorentzVector Momentum = toObserverFrame*Coordinate;
246  Hadron->SetFormationTime(TimeOftheStringCreation + Momentum.e() - fermi/c_light);
247  G4ThreeVector aPosition(Momentum.vect());
248  Hadron->SetPosition(PositionOftheStringCreation+aPosition);
249 
250 
251  if ( currentString->GetDecayDirection() > 0 )
252  {
253  LeftVector->push_back(Hadron);
254  } else
255  {
256  RightVector->push_back(Hadron);
257  }
258  delete currentString;
259  currentString=newString;
260  } else {
261  if ( newString ) delete newString; //AR-09Mar2017 Fixed memory leak
262  }
263 
264  currentString->LorentzRotate(toObserverFrame);
265  }; // while ( (! StopFragmenting(currentString)) ...
266 
267  if ( loopCounter >= maxNumberOfLoops ) {
268  inner_success=false;
269  }
270 
271  // Split remaining string into 2 final hadrons.
272 #ifdef debug_LUNDfragmentation
273  if(inner_success) G4cout<<"Split remaining string into 2 final hadrons."<<G4endl;
274 #endif
275  if ( inner_success && SplitLast(currentString, LeftVector, RightVector) )
276  {
277  final_success = true;
278  }
279 
280  delete currentString;
281  } // End of the loop where we try to fragment the string.
282  return final_success;
283 }
284 
285 //----------------------------------------------------------------------------------------
286 G4bool G4LundStringFragmentation::StopFragmenting(const G4FragmentingString * const string)
287 {
288  SetMinimalStringMass(string);
289  if (string->FourQuarkString())
290  {
291  return G4UniformRand() < G4Exp(-0.0005*(string->Mass() - MinimalStringMass));
292  } else {
293 
294  G4bool Result = G4UniformRand() <
295  G4Exp(-0.66e-6*(string->Mass()*string->Mass() - MinimalStringMass*MinimalStringMass));
296 #ifdef debug_LUNDfragmentation
297  G4cout<<"StopFragmenting MinimalStringMass string->Mass() "<<MinimalStringMass<<" "<<string->Mass()<<G4endl;
298  G4cout<<"StopFragmenting - Yes/No "<<Result<<G4endl;
299 #endif
300  return Result;
301  }
302 }
303 
304 //-----------------------------------------------------------------------------
305 G4KineticTrack * G4LundStringFragmentation::Splitup(
306  G4FragmentingString *string,
307  G4FragmentingString *&newString)
308 {
309 #ifdef debug_LUNDfragmentation
310  G4cout<<G4endl;
311  G4cout<<"Start SplitUP (G4VLongitudinalStringDecay) ========================="<<G4endl;
312  G4cout<<"String partons: " <<string->GetLeftParton()->GetPDGEncoding()<<" "
313  <<string->GetRightParton()->GetPDGEncoding()<<" "
314  <<"Direction " <<string->GetDecayDirection()<<G4endl;
315 #endif
316 
317  //... random choice of string end to use for creating the hadron (decay)
318  G4int SideOfDecay = (G4UniformRand() < 0.5)? 1: -1;
319  if (SideOfDecay < 0)
320  {
321  string->SetLeftPartonStable();
322  } else
323  {
324  string->SetRightPartonStable();
325  }
326 
327  G4ParticleDefinition *newStringEnd;
328  G4ParticleDefinition * HadronDefinition;
329 
330  G4double StringMass=string->Mass();
331  G4double HadronMass=0.;
332  G4double MinMass =0.;
333  G4int Iter(0), MaxIter(1);//000); // Uzhi Jan. 2017
334 
335 // Uzhi Nov. 2017 Start
336 G4double ProbDqADq = GetDiquarkSuppress();
337 G4double ProbSaS = 1.0 - 2.0 * GetStrangeSuppress();
338 
339 #ifdef debug_LUNDfragmentation
340 G4cout<<"StrMass DiquarkSuppression "<<StringMass<<" "<<GetDiquarkSuppress()<<G4endl;
341 #endif
342 
343 G4int NumberOfpossibleBaryons = 2;
344 
345 if(string->GetLeftParton()->GetParticleSubType() != "quark") NumberOfpossibleBaryons++;
346 if(string->GetRightParton()->GetParticleSubType() != "quark") NumberOfpossibleBaryons++;
347 
348 G4double ActualProb = ProbDqADq ;
349  ActualProb *= (1.0-sqr(NumberOfpossibleBaryons*1400.0/StringMass));
350 
351 SetDiquarkSuppression(ActualProb);
352 
353 G4double Mth = 1250.0; // 2 Mk + Mpi 1130.0 * 1.5 1500
354 if ( NumberOfpossibleBaryons == 3 ){Mth = 2000.0;} // Mlambda/Msigma + Mk + Mpi 2520 3400
355 else if( NumberOfpossibleBaryons == 4 ){Mth = 2550.0;} // 2 Mlambda/Msigma + Mk + Mpi 4620
356 else {}
357 
358 // ActualProb = ProbSaS * (1.0 - sqr(sqr(Mth/StringMass))); if( ActualProb < 0.0 ) ActualProb = 0.0;
359  ActualProb = ProbSaS * (1.0 - G4Pow::GetInstance()->powA( Mth/StringMass , 4.0 ));
360 // ActualProb = ProbSaS * (1.0 - G4Exp( -(StringMass -Mth)*500. ));
361 if( ActualProb < 0.0 ) ActualProb = 0.0;
362  SetStrangenessSuppression((1.0-ActualProb)/2.0);
363 
364 #ifdef debug_LUNDfragmentation
365 G4cout<<"StrMass DiquarkSuppression Pow "<<StringMass<<" "<<GetDiquarkSuppress()<<G4endl;
366 #endif
367 // Uzhi Nov. 2017 End
368 
369  do {
370 
371  if (string->DecayIsQuark())
372  {
373  HadronDefinition= QuarkSplitup(string->GetDecayParton(), newStringEnd);
374  } else {
375  HadronDefinition= DiQuarkSplitup(string->GetDecayParton(), newStringEnd);
376  }
377 
378 #ifdef debug_LUNDfragmentation
379 G4cout<<"Iter "<<Iter<<G4endl;
380  G4cout<<"The parton "<<string->GetDecayParton()->GetPDGEncoding()<<" "
381  <<" produces hadron "<<HadronDefinition->GetParticleName()
382  <<" and is transformed to "<<newStringEnd->GetPDGEncoding()<<G4endl;
383  G4cout<<"The side of the string decay Left/Right (1/-1) "<<SideOfDecay<<G4endl;
384 #endif
385 // create new String from old, ie. keep Left and Right order, but replace decay
386 
387  if ( newString ) delete newString; //AR-09Mar2017 Fixed memory leak
388 
389  newString=new G4FragmentingString(*string,newStringEnd); // To store possible
390  // quark containt of new string
391  SetMinimalStringMass(newString); MinMass = MinimalStringMass;
392  HadronMass = HadronDefinition->GetPDGMass();
393 
394  if(Iter >= MaxIter) {SetDiquarkSuppression(ProbDqADq);
395  SetStrangenessSuppression((1.0-ProbSaS)/2.0);
396  G4KineticTrack * Hadron =0;
397  return Hadron;} // Uzhi Nov. 2017
398  Iter++;
399 
400  } while(HadronMass + MinMass > StringMass);
401 
402 SetDiquarkSuppression(ProbDqADq); // Uzhi Nov. 2017
403 SetStrangenessSuppression((1.0-ProbSaS)/2.0); // Uzhi Nov. 2017
404 
405 #ifdef debug_LUNDfragmentation
406  G4cout<<"An attempt to determine its energy (SplitEandP)"<<G4endl;
407 #endif
408  G4LorentzVector* HadronMomentum=SplitEandP(HadronDefinition, string, newString);
409 
410  delete newString; newString=0;
411 
412  G4KineticTrack * Hadron =0;
413  if ( HadronMomentum != 0 ) {
414 
415 #ifdef debug_LUNDfragmentation
416  G4cout<<"The attempt was successful"<<G4endl;
417 #endif
418  G4ThreeVector Pos;
419  Hadron = new G4KineticTrack(HadronDefinition, 0,Pos, *HadronMomentum);
420 
421  if ( newString ) delete newString; //AR-09Mar2017 Fixed memory leak
422 
423  newString=new G4FragmentingString(*string,newStringEnd,
424  HadronMomentum);
425  delete HadronMomentum;
426  }
427  else
428  {
429 
430 #ifdef debug_LUNDfragmentation
431  G4cout<<"The attempt was not successful !!!"<<G4endl;
432 #endif
433  }
434 
435 #ifdef debug_LUNDfragmentation
436  G4cout<<"End SplitUP (G4VLongitudinalStringDecay) ====================="<<G4endl;
437 #endif
438 
439  return Hadron;
440 }
441 
442 //-----------------------------------------------------------------------------
443 G4ParticleDefinition * G4LundStringFragmentation::DiQuarkSplitup(
444  G4ParticleDefinition* decay,
445  G4ParticleDefinition *&created)
446 {
447 G4double StrSup=GetStrangeSuppress(); // Uzhi Dec. 2017
448  //... can Diquark break or not?
449  if (G4UniformRand() < DiquarkBreakProb ){
450 
451  //... Diquark break
452  G4int stableQuarkEncoding = decay->GetPDGEncoding()/1000;
453  G4int decayQuarkEncoding = (decay->GetPDGEncoding()/100)%10;
454  if (G4UniformRand() < 0.5)
455  {
456  G4int Swap = stableQuarkEncoding;
457  stableQuarkEncoding = decayQuarkEncoding;
458  decayQuarkEncoding = Swap;
459  }
460 
461  G4int IsParticle=(decayQuarkEncoding>0) ? -1 : +1;
462  // if we have a quark, we need antiquark)
463 
464 StrangeSuppress=(1.0 - 0.16)/2.0; // Uzhi Dec. 2017
465  pDefPair QuarkPair = CreatePartonPair(IsParticle,false); // no diquarks wanted
466 StrangeSuppress=StrSup; // Uzhi Dec. 2017
467 
468  //... Build new Diquark
469  G4int QuarkEncoding=QuarkPair.second->GetPDGEncoding();
470  G4int i10 = std::max(std::abs(QuarkEncoding), std::abs(stableQuarkEncoding));
471  G4int i20 = std::min(std::abs(QuarkEncoding), std::abs(stableQuarkEncoding));
472  G4int spin = (i10 != i20 && G4UniformRand() <= 0.5)? 1 : 3;
473  G4int NewDecayEncoding = -1*IsParticle*(i10 * 1000 + i20 * 100 + spin);
474  created = FindParticle(NewDecayEncoding);
475  G4ParticleDefinition * decayQuark=FindParticle(decayQuarkEncoding);
476  G4ParticleDefinition * had=hadronizer->Build(QuarkPair.first, decayQuark);
477  return had;
478 // return hadronizer->Build(QuarkPair.first, decayQuark);
479 
480  } else {
481  //... Diquark does not break
482 
483  G4int IsParticle=(decay->GetPDGEncoding()>0) ? +1 : -1;
484  // if we have a diquark, we need quark)
485 
486  pDefPair QuarkPair = CreatePartonPair(IsParticle,false); // no diquarks wanted
487 
488  created = QuarkPair.second;
489 
490  G4ParticleDefinition * had=hadronizer->Build(QuarkPair.first, decay);
491  return had;
492  }
493 }
494 
495 //-----------------------------------------------------------------------------
496 G4LorentzVector * G4LundStringFragmentation::SplitEandP(G4ParticleDefinition * pHadron,
497  G4FragmentingString * string,
498  G4FragmentingString * newString)
499 {
500  G4LorentzVector String4Momentum=string->Get4Momentum();
501  G4double StringMT2=string->MassT2();
502  G4double StringMT =std::sqrt(StringMT2);
503 
504  G4double HadronMass = pHadron->GetPDGMass();
505 
506  SetMinimalStringMass(newString);
507 
508 #ifdef debug_LUNDfragmentation
509  G4cout<<G4endl<<"Start LUND SplitEandP "<<G4endl;
510  G4cout<<"String 4 mom, String M and Mt "<<String4Momentum<<" "<<String4Momentum.mag()<<" "<<std::sqrt(StringMT2)<<G4endl;
511  G4cout<<"Hadron "<<pHadron->GetParticleName()<<G4endl;
512  G4cout<<"HadM MinimalStringMassLeft StringM hM+sM "<<HadronMass<<" "<<MinimalStringMass<<" "
513  <<String4Momentum.mag()<<" "<<HadronMass+MinimalStringMass<<G4endl;
514 #endif
515 
516  if(HadronMass + MinimalStringMass > string->Mass())
517  {
518 #ifdef debug_LUNDfragmentation
519  G4cout<<"Mass of the string is not sufficient to produce the hadron!"<<G4endl;
520 #endif
521  return 0;
522  }// have to start all over!
523 
524  String4Momentum.setPz(0.);
525  G4ThreeVector StringPt=String4Momentum.vect();
526 
527  // calculate and assign hadron transverse momentum component HadronPx and HadronPy
528  G4ThreeVector HadronPt , RemSysPt;
529  G4double HadronMassT2, ResidualMassT2;
530  G4double HadronMt, Pt, Pt2, phi; // Uzhi Nov. 2016 Introduction of Mt distribution
531  //... sample Pt of the hadron
532  G4int attempt=0;
533  do
534  {
535  attempt++; if(attempt > StringLoopInterrupt) {return 0;}
536  // Uzhi Nov. 2016 Introduction of Mt distribution
537  HadronMt = HadronMass - 300.0*G4Log(G4UniformRand()); //200 -> 300 // Uzhi Dec. 2017
538  Pt2 = sqr(HadronMt)-sqr(HadronMass); Pt=std::sqrt(Pt2);
539  phi = 2.*pi*G4UniformRand();
540  G4ThreeVector SampleQuarkPtw= G4ThreeVector(Pt * std::cos(phi),Pt * std::sin(phi),0);
541  // Uzhi Nov. 2016 Introduction of Mt distribution
542  HadronPt =SampleQuarkPtw + string->DecayPt(); // SampleQuarkPt() // Uzhi Nov. 2016
543  HadronPt.setZ(0);
544  RemSysPt = StringPt - HadronPt;
545 
546  HadronMassT2 = sqr(HadronMass) + HadronPt.mag2();
547  ResidualMassT2=sqr(MinimalStringMass) + RemSysPt.mag2();
548 
549  } while(std::sqrt(HadronMassT2) + std::sqrt(ResidualMassT2) > StringMT);
550 
551  //... sample z to define hadron longitudinal momentum and energy
552  //... but first check the available phase space
553 
554  G4double Pz2 = (sqr(StringMT2 - HadronMassT2 - ResidualMassT2) -
555  4*HadronMassT2 * ResidualMassT2)/4./StringMT2;
556 
557  if(Pz2 < 0 ) {return 0;} // have to start all over!
558 
559  //... then compute allowed z region z_min <= z <= z_max
560 
561  G4double Pz = std::sqrt(Pz2);
562  G4double zMin = (std::sqrt(HadronMassT2+Pz2) - Pz)/std::sqrt(StringMT2);
563 // G4double zMin = (std::sqrt(HadronMassT2+Pz2) - 0.)/std::sqrt(StringMT2); // For testing purposes
564  G4double zMax = (std::sqrt(HadronMassT2+Pz2) + Pz)/std::sqrt(StringMT2);
565 
566  if (zMin >= zMax) return 0; // have to start all over!
567 
568  G4double z = GetLightConeZ(zMin, zMax,
569  string->GetDecayParton()->GetPDGEncoding(), pHadron,
570  HadronPt.x(), HadronPt.y());
571 
572  //... now compute hadron longitudinal momentum and energy
573  // longitudinal hadron momentum component HadronPz
574 
575  HadronPt.setZ(0.5* string->GetDecayDirection() *
576  (z * string->LightConeDecay() -
577  HadronMassT2/(z * string->LightConeDecay())));
578  G4double HadronE = 0.5* (z * string->LightConeDecay() +
579  HadronMassT2/(z * string->LightConeDecay()));
580 
581  G4LorentzVector * a4Momentum= new G4LorentzVector(HadronPt,HadronE);
582 
583 #ifdef debug_LUNDfragmentation
584  G4cout<<G4endl<<" string->GetDecayDirection() "<<string->GetDecayDirection()<<G4endl<<G4endl;
585  G4cout<<"string->LightConeDecay() "<<string->LightConeDecay()<<G4endl;
586  G4cout<<"HadronPt,HadronE "<<HadronPt<<" "<<HadronE<<G4endl;
587  G4cout<<"String4Momentum "<<String4Momentum<<G4endl;
588  G4cout<<"Out of LUND SplitEandP "<<G4endl<<G4endl;
589 #endif
590 
591  return a4Momentum;
592 }
593 
594 //-----------------------------------------------------------------------------------------
595 G4double G4LundStringFragmentation::GetLightConeZ(G4double zmin, G4double zmax,
596  G4int PDGEncodingOfDecayParton,
597  G4ParticleDefinition* pHadron,
598  G4double Px, G4double Py)
599 {
600  G4double Mass = pHadron->GetPDGMass();
601  G4int HadronEncoding=std::abs(pHadron->GetPDGEncoding());
602 
603  G4double Mt2 = Px*Px + Py*Py + Mass*Mass;
604 
605  G4double Alund, Blund;
606  G4double zOfMaxyf(0.), maxYf(1.), z(0.), yf(1.);
607 
608 // if(! ( ((std::abs(PDGEncodingOfDecayParton) > 1000) && (HadronEncoding > 1000)) || (HadronEncoding == 113) ) ) // Uzhi Dec. 2017
609  if(!((std::abs(PDGEncodingOfDecayParton) > 1000) && (HadronEncoding > 1000)) )
610  { // ---------------- Quark fragmentation and qq-> meson ----------------------
611  Alund=1.;
612  Blund=0.7/GeV/GeV;
613 
614  G4double BMt2 = Blund*Mt2;
615  if(Alund == 1.0) {
616  zOfMaxyf=BMt2/(Blund*Mt2 + 1.);}
617  else {
618  zOfMaxyf = ((1.0+BMt2) - std::sqrt(sqr(1.0-BMt2) + 4.0*BMt2*Alund))/2.0/(1.-Alund);
619  }
620 
621  if(zOfMaxyf < zmin) {zOfMaxyf=zmin;}
622  if(zOfMaxyf > zmax) {zOfMaxyf=zmax;}
623  maxYf=(1-zOfMaxyf)/zOfMaxyf * G4Exp(-Blund*Mt2/zOfMaxyf);
624 
625  const G4int maxNumberOfLoops = 1000;
626  G4int loopCounter = 0;
627  do
628  {
629  z = zmin + G4UniformRand()*(zmax-zmin);
630 // yf = (1-z)/z * G4Exp(-Blund*Mt2/z);
631  yf = G4Pow::GetInstance()->powA(1.0-z,Alund)/z*G4Exp(-BMt2/z);
632  }
633  while ( (G4UniformRand()*maxYf > yf) && ++loopCounter < maxNumberOfLoops );
634  if ( loopCounter >= maxNumberOfLoops ) {
635  z = 0.5*(zmin + zmax); // Just a value between zmin and zmax, no physics considerations at all!
636  }
637  return z;
638  }
639 
640 // if((std::abs(PDGEncodingOfDecayParton) > 1000) || (HadronEncoding == 113)) // Uzhi dec. 2017
641  if(std::abs(PDGEncodingOfDecayParton) > 1000)
642  {
643  G4double an = 2.5;
644 //Uzhi if(pHadron->GetPDGIsospin() > 0.5) an=0.5;//0.75;
645  an +=(sqr(Px)+sqr(Py))/sqr(GeV)-0.5;
646  z=zmin + (zmax-zmin)*G4Pow::GetInstance()->powA(G4UniformRand(),1./an);
647  }
648 
649  return z;
650 
651 }
652 //
653 
654 // ####################################################
655 
656 //----------------------------------------------------------------------------------------------------------
657 G4bool G4LundStringFragmentation::SplitLast(G4FragmentingString * string,
658  G4KineticTrackVector * LeftVector,
659  G4KineticTrackVector * RightVector)
660 {
661  //... perform last cluster decay
662 
663 #ifdef debug_LUNDfragmentation
664  G4cout<<G4endl<<"Split last-----------------------------------------"<<G4endl;
665  G4cout<<"Left "<<string->GetLeftParton()->GetPDGEncoding()<<" "<<string->GetPleft()<<G4endl;
666  G4cout<<"Right "<<string->GetRightParton()->GetPDGEncoding()<<" "<<string->GetPright()<<G4endl;
667  G4cout<<"String4mom "<<string->GetPstring()<<" "<<string->GetPstring().mag()<<G4endl;
668 #endif
669 
670  G4LorentzVector Str4Mom=string->Get4Momentum();
671 
672  G4LorentzRotation toCms=string->TransformToAlignedCms();
673  G4LorentzRotation toObserverFrame= toCms.inverse();
674 
675  G4double StringMass=string->Mass();
676 
677  G4ParticleDefinition * LeftHadron(0), * RightHadron(0);
678 
679  NumberOf_FS=0;
680  for(G4int i=0; i<35; i++) {FS_Weight[i]=0.;}
681 
682  G4int sampledState = 0;
683 
684 #ifdef debug_LUNDfragmentation
685  G4cout<<"StrMass "<<StringMass<<" q's "
686  <<string->GetLeftParton()->GetParticleName()<<" "
687  <<string->GetRightParton()->GetParticleName()<<G4endl;
688 #endif
689 
690  string->SetLeftPartonStable(); // to query quark contents..
691 
692  if (string->FourQuarkString() )
693  {
694  // The string is qq-qqbar type. Diquarks are on the string ends
695  //G4cout<<"The string is qq-qqbar type. Diquarks are on the string ends"<<G4endl;
696  if(StringMass-MinimalStringMass < 0.)
697  {
698  if (! Diquark_AntiDiquark_belowThreshold_lastSplitting(string, LeftHadron, RightHadron) )
699  {
700  return false;
701  }
702  } else
703  {
704  Diquark_AntiDiquark_aboveThreshold_lastSplitting(string, LeftHadron, RightHadron);
705 
706  if(NumberOf_FS == 0) return false;
707 
708  sampledState = SampleState();
709  if(string->GetLeftParton()->GetPDGEncoding() < 0)
710  {
711  LeftHadron =FS_LeftHadron[sampledState];
712  RightHadron=FS_RightHadron[sampledState];
713  } else
714  {
715  LeftHadron =FS_RightHadron[sampledState];
716  RightHadron=FS_LeftHadron[sampledState];
717  }
718  }
719  } else
720  {
721  if (string->DecayIsQuark() && string->StableIsQuark() )
722  { //... there are quarks on cluster ends
723 #ifdef debug_LUNDfragmentation
724  G4cout<<"Q Q string LastSplit"<<G4endl;
725 #endif
726  Quark_AntiQuark_lastSplitting(string, LeftHadron, RightHadron);
727 
728  if(NumberOf_FS == 0) return false;
729  sampledState = SampleState();
730  if(string->GetLeftParton()->GetPDGEncoding() < 0)
731  {
732  LeftHadron =FS_RightHadron[sampledState];
733  RightHadron=FS_LeftHadron[sampledState];
734  } else
735  {
736  LeftHadron =FS_LeftHadron[sampledState];
737  RightHadron=FS_RightHadron[sampledState];
738  }
739  } else
740  { //... there is a Diquark on one of the cluster ends
741 #ifdef debug_LUNDfragmentation
742  G4cout<<"DiQ Q string Last Split"<<G4endl;
743 #endif
744 
745  Quark_Diquark_lastSplitting(string, LeftHadron, RightHadron);
746 
747  if(NumberOf_FS == 0) return false;
748  sampledState = SampleState();
749 
750  if(string->GetLeftParton()->GetParticleSubType() == "quark")
751  {
752  LeftHadron =FS_LeftHadron[sampledState];
753  RightHadron=FS_RightHadron[sampledState];
754  } else
755  {
756  LeftHadron =FS_RightHadron[sampledState];
757  RightHadron=FS_LeftHadron[sampledState];
758  }
759  }
760 
761  } // End of if(!string->FourQuarkString())
762 
763 #ifdef debug_LUNDfragmentation
764  G4cout<<"Sampled hadrons: "<<LeftHadron->GetParticleName()<<" "<<RightHadron->GetParticleName()<<G4endl;
765 #endif
766 
767  G4LorentzVector P_left =string->GetPleft(), P_right = string->GetPright();
768 
769  G4LorentzVector LeftMom, RightMom;
770  G4ThreeVector Pos;
771 
772  Sample4Momentum(&LeftMom, LeftHadron->GetPDGMass(),
773  &RightMom, RightHadron->GetPDGMass(),
774  StringMass);
775 
776 // Uzhi 24 June 2017: Sample4Momentum ascribes LeftMom.pz() along positive Z axis for baryons in many cases.
777 // It must be negative in case when the system is moving against Z axis.
778 
779  if(!(string->DecayIsQuark() && string->StableIsQuark() ))
780  { // Only for qq - q, q - qq, and qq - qqbar -------------------
781 
782  if(std::abs(string->GetLeftParton()->GetPDGEncoding()) > 1000)
783  {
784  if(P_left.z() <= 0.) {G4LorentzVector tmp = LeftMom; LeftMom=RightMom; RightMom=tmp;}
785  }
786  else
787  {
788  if(P_right.z() >= 0.) {G4LorentzVector tmp = LeftMom; LeftMom=RightMom; RightMom=tmp;}
789  }
790  }
791 
792  LeftMom *=toObserverFrame;
793  RightMom*=toObserverFrame;
794 
795  LeftVector->push_back(new G4KineticTrack(LeftHadron, 0, Pos, LeftMom));
796  RightVector->push_back(new G4KineticTrack(RightHadron, 0, Pos, RightMom));
797 
798  string->LorentzRotate(toObserverFrame);
799  return true;
800 
801 }
802 
803 //----------------------------------------------------------------------------------------
804 G4bool G4LundStringFragmentation::
805 Diquark_AntiDiquark_belowThreshold_lastSplitting(G4FragmentingString * & string,
806  G4ParticleDefinition * & LeftHadron,
807  G4ParticleDefinition * & RightHadron)
808 {
809  G4double StringMass = string->Mass();
810 
811  G4int cClusterInterrupt = 0;
812  do
813  {
814  if (cClusterInterrupt++ >= ClusterLoopInterrupt)
815  {
816  return false;
817  }
818 
819  G4int LeftQuark1= string->GetLeftParton()->GetPDGEncoding()/1000;
820  G4int LeftQuark2=(string->GetLeftParton()->GetPDGEncoding()/100)%10;
821 
822  G4int RightQuark1= string->GetRightParton()->GetPDGEncoding()/1000;
823  G4int RightQuark2=(string->GetRightParton()->GetPDGEncoding()/100)%10;
824 
825  if(G4UniformRand()<0.5)
826  {
827  LeftHadron =hadronizer->Build(FindParticle( LeftQuark1),
828  FindParticle(RightQuark1));
829  RightHadron=hadronizer->Build(FindParticle( LeftQuark2),
830  FindParticle(RightQuark2));
831  } else
832  {
833  LeftHadron =hadronizer->Build(FindParticle( LeftQuark1),
834  FindParticle(RightQuark2));
835  RightHadron=hadronizer->Build(FindParticle( LeftQuark2),
836  FindParticle(RightQuark1));
837  }
838 
839  //... repeat procedure, if mass of cluster is too low to produce hadrons
840  //... ClusterMassCut = 0.15*GeV model parameter
841  }
842  while ((StringMass <= LeftHadron->GetPDGMass() + RightHadron->GetPDGMass()));
843 
844  return true;
845 }
846 
847 //----------------------------------------------------------------------------------------
848 G4bool G4LundStringFragmentation::
849 Diquark_AntiDiquark_aboveThreshold_lastSplitting(G4FragmentingString * & string,
850  G4ParticleDefinition * & LeftHadron,
851  G4ParticleDefinition * & RightHadron)
852 {
853  // StringMass-MinimalStringMass > 0. Creation of 2 baryons is possible ----
854 
855  G4double StringMass = string->Mass();
856  G4double StringMassSqr= sqr(StringMass);
857  G4ParticleDefinition * Di_Quark;
858  G4ParticleDefinition * Anti_Di_Quark;
859 
860  if(string->GetLeftParton()->GetPDGEncoding() < 0)
861  {
862  Anti_Di_Quark =string->GetLeftParton();
863  Di_Quark=string->GetRightParton();
864  } else
865  {
866  Anti_Di_Quark =string->GetRightParton();
867  Di_Quark=string->GetLeftParton();
868  }
869 
870  G4int IDAnti_di_quark =Anti_Di_Quark->GetPDGEncoding();
871  G4int AbsIDAnti_di_quark =std::abs(IDAnti_di_quark);
872  G4int IDdi_quark =Di_Quark->GetPDGEncoding();
873  G4int AbsIDdi_quark =std::abs(IDdi_quark);
874 
875  G4int ADi_q1=AbsIDAnti_di_quark/1000;
876  G4int ADi_q2=(AbsIDAnti_di_quark-ADi_q1*1000)/100;
877 
878  G4int Di_q1=AbsIDdi_quark/1000;
879  G4int Di_q2=(AbsIDdi_quark-Di_q1*1000)/100;
880 
881  NumberOf_FS=0;
882  for(G4int ProdQ=1; ProdQ < 4; ProdQ++)
883  {
884  G4int StateADiQ=0;
885  const G4int maxNumberOfLoops = 1000;
886  G4int loopCounter = 0;
887  do // while(Meson[AbsIDquark-1][ProdQ-1][StateQ]<>0);
888  {
889  LeftHadron=G4ParticleTable::GetParticleTable()->FindParticle(
890  -Baryon[ADi_q1-1][ADi_q2-1][ProdQ-1][StateADiQ]);
891  G4double LeftHadronMass=LeftHadron->GetPDGMass();
892 
893  G4int StateDiQ=0;
894  const G4int maxNumberOfInternalLoops = 1000;
895  G4int internalLoopCounter = 0;
896  do // while(Baryon[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]<>0);
897  {
898  RightHadron=G4ParticleTable::GetParticleTable()->FindParticle(
899  +Baryon[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]);
900 
901  G4double RightHadronMass=RightHadron->GetPDGMass();
902 
903  if(StringMass > LeftHadronMass + RightHadronMass)
904  {
905  if ( NumberOf_FS > 34 ) {
906  G4ExceptionDescription ed;
907  ed << " NumberOf_FS exceeds its limit: NumberOf_FS=" << NumberOf_FS << G4endl;
908  G4Exception( "G4LundStringFragmentation::Diquark_AntiDiquark_aboveThreshold_lastSplitting ",
909  "HAD_LUND_001", JustWarning, ed );
910  NumberOf_FS = 34;
911  }
912 
913  G4double FS_Psqr=lambda(StringMassSqr,sqr(LeftHadronMass),
914  sqr(RightHadronMass));
915  //FS_Psqr=1.;
916  FS_Weight[NumberOf_FS]=std::sqrt(FS_Psqr)*FS_Psqr*
917  BaryonWeight[ADi_q1-1][ADi_q2-1][ProdQ-1][StateADiQ]*
918  BaryonWeight[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]*
919  Prob_QQbar[ProdQ-1];
920 
921  FS_LeftHadron[NumberOf_FS] = LeftHadron;
922  FS_RightHadron[NumberOf_FS]= RightHadron;
923  NumberOf_FS++;
924  } // End of if(StringMass > LeftHadronMass + RightHadronMass)
925 
926  StateDiQ++;
927 
928  } while( (Baryon[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]!=0) &&
929  ++internalLoopCounter < maxNumberOfInternalLoops );
930  if ( internalLoopCounter >= maxNumberOfInternalLoops ) {
931  return false;
932  }
933 
934  StateADiQ++;
935  } while( (Baryon[ADi_q1-1][ADi_q2-1][ProdQ-1][StateADiQ]!=0) &&
936  ++loopCounter < maxNumberOfLoops );
937  if ( loopCounter >= maxNumberOfLoops ) {
938  return false;
939  }
940  } // End of for(G4int ProdQ=1; ProdQ < 4; ProdQ++)
941 
942  return true;
943 }
944 
945 //----------------------------------------------------------------------------------------
946 G4bool G4LundStringFragmentation::
947 Quark_Diquark_lastSplitting(G4FragmentingString * & string,
948  G4ParticleDefinition * & LeftHadron,
949  G4ParticleDefinition * & RightHadron)
950 {
951  G4double StringMass = string->Mass();
952  G4double StringMassSqr= sqr(StringMass);
953 
954  G4ParticleDefinition * Di_Quark;
955  G4ParticleDefinition * Quark;
956 
957  if(string->GetLeftParton()->GetParticleSubType()== "quark")
958  {
959  Quark =string->GetLeftParton();
960  Di_Quark=string->GetRightParton();
961  } else
962  {
963  Quark =string->GetRightParton();
964  Di_Quark=string->GetLeftParton();
965  }
966 
967  G4int IDquark =Quark->GetPDGEncoding();
968  G4int AbsIDquark =std::abs(IDquark);
969  G4int IDdi_quark =Di_Quark->GetPDGEncoding();
970  G4int AbsIDdi_quark=std::abs(IDdi_quark);
971  G4int Di_q1=AbsIDdi_quark/1000;
972  G4int Di_q2=(AbsIDdi_quark-Di_q1*1000)/100;
973 
974  G4int SignDiQ= 1;
975  if(IDdi_quark < 0) SignDiQ=-1;
976 
977  NumberOf_FS=0;
978  for(G4int ProdQ=1; ProdQ < 4; ProdQ++)
979  {
980  G4int SignQ;
981  if(IDquark > 0)
982  { SignQ=-1;
983  if(IDquark == 2) SignQ= 1;
984  if((IDquark == 1) && (ProdQ == 3)) SignQ= 1; // K0
985  if((IDquark == 3) && (ProdQ == 1)) SignQ=-1; // K0bar
986  } else
987  {
988  SignQ= 1;
989  if(IDquark == -2) SignQ=-1;
990  if((IDquark ==-1) && (ProdQ == 3)) SignQ=-1; // K0bar
991  if((IDquark ==-3) && (ProdQ == 1)) SignQ= 1; // K0
992  }
993 
994  if(AbsIDquark == ProdQ) SignQ= 1;
995 
996  G4int StateQ=0;
997  const G4int maxNumberOfLoops = 1000;
998  G4int loopCounter = 0;
999  do // while(Meson[AbsIDquark-1][ProdQ-1][StateQ]<>0);
1000  {
1001  LeftHadron=G4ParticleTable::GetParticleTable()->FindParticle(SignQ*
1002  Meson[AbsIDquark-1][ProdQ-1][StateQ]);
1003  G4double LeftHadronMass=LeftHadron->GetPDGMass();
1004 
1005  G4int StateDiQ=0;
1006  const G4int maxNumberOfInternalLoops = 1000;
1007  G4int internalLoopCounter = 0;
1008  do // while(Baryon[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]<>0);
1009  {
1010  RightHadron=G4ParticleTable::GetParticleTable()->FindParticle(SignDiQ*
1011  Baryon[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]);
1012  G4double RightHadronMass=RightHadron->GetPDGMass();
1013 
1014  if(StringMass > LeftHadronMass + RightHadronMass)
1015  {
1016  if ( NumberOf_FS > 34 ) {
1017  G4ExceptionDescription ed;
1018  ed << " NumberOf_FS exceeds its limit: NumberOf_FS=" << NumberOf_FS << G4endl;
1019  G4Exception( "G4LundStringFragmentation::Quark_Diquark_lastSplitting ",
1020  "HAD_LUND_002", JustWarning, ed );
1021  NumberOf_FS = 34;
1022  }
1023 
1024  G4double FS_Psqr=lambda(StringMassSqr,sqr(LeftHadronMass),
1025  sqr(RightHadronMass));
1026  FS_Weight[NumberOf_FS]=std::sqrt(FS_Psqr)*
1027  MesonWeight[AbsIDquark-1][ProdQ-1][StateQ]*
1028  BaryonWeight[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]*
1029  Prob_QQbar[ProdQ-1];
1030 
1031  FS_LeftHadron[NumberOf_FS] = LeftHadron;
1032  FS_RightHadron[NumberOf_FS]= RightHadron;
1033 
1034  NumberOf_FS++;
1035  } // End of if(StringMass > LeftHadronMass + RightHadronMass)
1036 
1037  StateDiQ++;
1038 
1039  } while( (Baryon[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]!=0) &&
1040  ++internalLoopCounter < maxNumberOfInternalLoops );
1041  if ( internalLoopCounter >= maxNumberOfInternalLoops ) {
1042  return false;
1043  }
1044 
1045  StateQ++;
1046  } while( (Meson[AbsIDquark-1][ProdQ-1][StateQ]!=0) &&
1047  ++loopCounter < maxNumberOfLoops ); /* Loop checking, 07.08.2015, A.Ribon */
1048 
1049  if ( loopCounter >= maxNumberOfLoops ) {
1050  return false;
1051  }
1052  } // End of for(G4int ProdQ=1; ProdQ < 4; ProdQ++)
1053 
1054  return true;
1055 }
1056 
1057 //----------------------------------------------------------------------------------------
1058 G4bool G4LundStringFragmentation::
1059 Quark_AntiQuark_lastSplitting(G4FragmentingString * & string,
1060  G4ParticleDefinition * & LeftHadron,
1061  G4ParticleDefinition * & RightHadron)
1062 {
1063  G4double StringMass = string->Mass();
1064  G4double StringMassSqr= sqr(StringMass);
1065 
1066  G4ParticleDefinition * Quark;
1067  G4ParticleDefinition * Anti_Quark;
1068 
1069  if(string->GetLeftParton()->GetPDGEncoding()>0)
1070  {
1071  Quark =string->GetLeftParton();
1072  Anti_Quark=string->GetRightParton();
1073  } else
1074  {
1075  Quark =string->GetRightParton();
1076  Anti_Quark=string->GetLeftParton();
1077  }
1078 
1079  G4int IDquark =Quark->GetPDGEncoding();
1080  G4int AbsIDquark =std::abs(IDquark);
1081  G4int IDanti_quark =Anti_Quark->GetPDGEncoding();
1082  G4int AbsIDanti_quark=std::abs(IDanti_quark);
1083 
1084  NumberOf_FS=0;
1085  for(G4int ProdQ=1; ProdQ < 4; ProdQ++)
1086  {
1087  G4int SignQ=-1;
1088  if(IDquark == 2) SignQ= 1;
1089  if((IDquark == 1) && (ProdQ == 3)) SignQ= 1; // K0
1090  if((IDquark == 3) && (ProdQ == 1)) SignQ=-1; // K0bar
1091  if(IDquark == ProdQ) SignQ= 1;
1092 
1093  G4int SignAQ= 1;
1094  if(IDanti_quark == -2) SignAQ=-1;
1095  if((IDanti_quark ==-1) && (ProdQ == 3)) SignAQ=-1; // K0bar
1096  if((IDanti_quark ==-3) && (ProdQ == 1)) SignAQ= 1; // K0
1097  if(AbsIDanti_quark == ProdQ) SignAQ= 1;
1098 
1099  G4int StateQ=0;
1100  const G4int maxNumberOfLoops = 1000;
1101  G4int loopCounter = 0;
1102  do // while(Meson[AbsIDquark-1][ProdQ-1][StateQ]<>0);
1103  {
1104  LeftHadron=G4ParticleTable::GetParticleTable()->FindParticle(SignQ*
1105  Meson[AbsIDquark-1][ProdQ-1][StateQ]);
1106  G4double LeftHadronMass=LeftHadron->GetPDGMass();
1107 
1108  G4int StateAQ=0;
1109  const G4int maxNumberOfInternalLoops = 1000;
1110  G4int internalLoopCounter = 0;
1111  do // while(Meson[AbsIDanti_quark-1][ProdQ-1][StateAQ]<>0);
1112  {
1113  RightHadron=G4ParticleTable::GetParticleTable()->FindParticle(SignAQ*
1114  Meson[AbsIDanti_quark-1][ProdQ-1][StateAQ]);
1115  G4double RightHadronMass=RightHadron->GetPDGMass();
1116 
1117  if(StringMass > LeftHadronMass + RightHadronMass)
1118  {
1119  if ( NumberOf_FS > 34 ) {
1120  G4ExceptionDescription ed;
1121  ed << " NumberOf_FS exceeds its limit: NumberOf_FS=" << NumberOf_FS << G4endl;
1122  G4Exception( "G4LundStringFragmentation::Quark_AntiQuark_lastSplitting ",
1123  "HAD_LUND_003", JustWarning, ed );
1124  NumberOf_FS = 34;
1125  }
1126 
1127  G4double FS_Psqr=lambda(StringMassSqr,sqr(LeftHadronMass),
1128  sqr(RightHadronMass));
1129  //FS_Psqr=1.;
1130  FS_Weight[NumberOf_FS]=std::sqrt(FS_Psqr)*
1131  MesonWeight[AbsIDquark-1][ProdQ-1][StateQ]*
1132  MesonWeight[AbsIDanti_quark-1][ProdQ-1][StateAQ]*
1133  Prob_QQbar[ProdQ-1];
1134 
1135  if(string->GetLeftParton()->GetPDGEncoding()>0)
1136  {
1137  FS_LeftHadron[NumberOf_FS] = RightHadron;
1138  FS_RightHadron[NumberOf_FS]= LeftHadron;
1139  } else
1140  {
1141  FS_LeftHadron[NumberOf_FS] = LeftHadron;
1142  FS_RightHadron[NumberOf_FS]= RightHadron;
1143  }
1144 
1145  NumberOf_FS++;
1146 
1147  } // End of if(StringMass > LeftHadronMass + RightHadronMass)
1148 
1149  StateAQ++;
1150  } while( (Meson[AbsIDanti_quark-1][ProdQ-1][StateAQ]!=0) &&
1151  ++internalLoopCounter < maxNumberOfInternalLoops );
1152  if ( internalLoopCounter >= maxNumberOfInternalLoops ) {
1153  return false;
1154  }
1155 
1156  StateQ++;
1157  } while( (Meson[AbsIDquark-1][ProdQ-1][StateQ]!=0) &&
1158  ++loopCounter < maxNumberOfLoops );
1159  if ( loopCounter >= maxNumberOfLoops ) {
1160  return false;
1161  }
1162  } // End of for(G4int ProdQ=1; ProdQ < 4; ProdQ++)
1163 
1164  return true;
1165 }
1166 
1167 //----------------------------------------------------------------------------------------------------------
1168 G4int G4LundStringFragmentation::SampleState(void)
1169 {
1170  if ( NumberOf_FS > 34 ) {
1171  G4ExceptionDescription ed;
1172  ed << " NumberOf_FS exceeds its limit: NumberOf_FS=" << NumberOf_FS << G4endl;
1173  G4Exception( "G4LundStringFragmentation::SampleState ", "HAD_LUND_004", JustWarning, ed );
1174  NumberOf_FS = 34;
1175  }
1176 
1177  G4double SumWeights=0.;
1178 
1179  for(G4int i=0; i<NumberOf_FS; i++) {SumWeights+=FS_Weight[i];}// G4cout<<i<<" "<<FS_Weight[i]<<G4endl;}
1180 
1181  G4double ksi=G4UniformRand();
1182  G4double Sum=0.;
1183  G4int indexPosition = 0;
1184 
1185  for(G4int i=0; i<NumberOf_FS; i++)
1186  {
1187  Sum+=(FS_Weight[i]/SumWeights);
1188  indexPosition=i;
1189  if(Sum >= ksi) break;
1190  }
1191  return indexPosition;
1192 }
1193 
1194 //----------------------------------------------------------------------------------------------------------
1195 void G4LundStringFragmentation::Sample4Momentum(G4LorentzVector* Mom, G4double Mass,
1196  G4LorentzVector* AntiMom, G4double AntiMass,
1197  G4double InitialMass)
1198 {
1199  // ------ Sampling of momenta of 2 last produced hadrons --------------------
1200  G4ThreeVector Pt;
1201  G4double MassMt, AntiMassMt;
1202  G4double AvailablePz, AvailablePz2;
1203 // G4double ProbIsotropy = sqr(sqr(938.0/InitialMass)); // Uzhi Oct. 2017
1204 #ifdef debug_LUNDfragmentation
1205  G4cout<<"Sampling of momenta of 2 last produced hadrons ----------------"<<G4endl;
1206  G4cout<<"Init Mass "<<InitialMass<<" FirstM "<<Mass<<" SecondM "<<AntiMass<<" ProbIsotropy "<<G4endl;
1207 #endif
1208 
1209  G4double r_val = sqr(InitialMass*InitialMass - Mass*Mass - AntiMass*AntiMass) -
1210  sqr(2.*Mass*AntiMass);
1211  G4double Pabs = (r_val > 0.)? std::sqrt(r_val)/(2.*InitialMass) : 0;
1212 
1213 // if((Mass > 930. || AntiMass > 930.) && (G4UniformRand() < ProbIsotropy)) // Uzhi Oct. 2017
1214 // { //If there is a baryon
1215 // // ----------------- Isotropic decay ------------------------------------
1216 // //... sample unit vector
1217 // G4double pz =1. - 2.*G4UniformRand();
1218 // G4double st = std::sqrt(1. - pz * pz)*Pabs;
1219 // G4double phi = 2.*pi*G4UniformRand();
1220 // G4double px = st*std::cos(phi);
1221 // G4double py = st*std::sin(phi);
1222 // pz *= Pabs;
1223 //
1224 // Mom->setPx(px); Mom->setPy(py); Mom->setPz(pz);
1225 // Mom->setE(std::sqrt(Pabs*Pabs + Mass*Mass));
1226 //
1227 // AntiMom->setPx(-px); AntiMom->setPy(-py); AntiMom->setPz(-pz);
1228 // AntiMom->setE (std::sqrt(Pabs*Pabs + AntiMass*AntiMass));
1229 // }
1230 // else
1231 // {
1232 
1233  const G4int maxNumberOfLoops = 1000;
1234 G4double SigmaQTw=SigmaQT; // Uzhi Oct. 2017
1235 if(Mass > 930. || AntiMass > 930.) SigmaQT *=(1.0-0.55*sqr((Mass+AntiMass)/InitialMass)); // Uzhi Oct. 2017
1236 
1237  G4int loopCounter = 0;
1238  do
1239  {
1240  // GF 22-May-09, limit sampled pt to allowed range
1241  Pt=SampleQuarkPt(Pabs); Pt.setZ(0); G4double Pt2=Pt.mag2();
1242 
1243  MassMt = std::sqrt( Mass * Mass + Pt2);
1244  AntiMassMt= std::sqrt(AntiMass * AntiMass + Pt2);
1245  }
1246  while( (InitialMass < MassMt + AntiMassMt) && ++loopCounter < maxNumberOfLoops );
1247 
1248 if(Mass > 930. || AntiMass > 930.) SigmaQT=SigmaQTw;; // Uzhi Oct. 2017
1249 
1250  if ( loopCounter >= maxNumberOfLoops ) {
1251  AvailablePz2 = 0.0;
1252  }
1253 
1254  AvailablePz2= sqr(InitialMass*InitialMass - sqr(MassMt) - sqr(AntiMassMt)) -
1255  4.*sqr(MassMt*AntiMassMt);
1256 
1257  AvailablePz2 /=(4.*InitialMass*InitialMass);
1258  AvailablePz = std::sqrt(AvailablePz2);
1259 
1260  G4double Px=Pt.getX();
1261  G4double Py=Pt.getY();
1262 
1263  Mom->setPx(Px); Mom->setPy(Py); Mom->setPz(AvailablePz);
1264  Mom->setE(std::sqrt(sqr(MassMt)+AvailablePz2));
1265 
1266  AntiMom->setPx(-Px); AntiMom->setPy(-Py); AntiMom->setPz(-AvailablePz);
1267  AntiMom->setE (std::sqrt(sqr(AntiMassMt)+AvailablePz2));
1268 // } // Uzhi Oct. 2017
1269 #ifdef debug_LUNDfragmentation
1270  G4cout<<"Fmass Mom "<<Mom->getX()<<" "<<Mom->getY()<<" "<<Mom->getZ()<<" "<<Mom->getT()<<G4endl;
1271  G4cout<<"Smass Mom "<<AntiMom->getX()<<" "<<AntiMom->getY()<<" "<<AntiMom->getZ()<<" "<<AntiMom->getT()<<G4endl;
1272 #endif
1273 }
1274 
1275 //------------------------------------------------------------------------
1276 G4double G4LundStringFragmentation::lambda(G4double S, G4double m1_Sqr, G4double m2_Sqr)
1277 {
1278  G4double lam = sqr(S - m1_Sqr - m2_Sqr) - 4.*m1_Sqr*m2_Sqr;
1279  return lam;
1280 }
1281 
1282 //-----------------------------------------------------------------------
1283 void G4LundStringFragmentation::SetMinMasses()
1284 {
1285 // ------ For estimation of a minimal string mass ---------------
1286  Mass_of_light_quark =140.*MeV;
1287  Mass_of_heavy_quark =500.*MeV;
1288  Mass_of_string_junction=720.*MeV;
1289 
1290  G4double minMQQbarStr[3][3] ={ {350.0, 350.0, 710.0}, //DDbar, DUbar, DSbar in MeV
1291  {350.0, 350.0, 710.0}, //UDbar, UUbar, USbar in Mev
1292  {710.0, 710.0,1070.0 }};//SDbar, SUbar, SSbar in MeV
1293  for(G4int i=0; i<3; i++){ for(G4int j=0; j<3; j++){minMassQQbarStr[i][j]=minMQQbarStr[i][j];};};
1294 
1295 
1296  G4double minMQDiQStr[3][3][3] = {{{1160., 1160., 1340.}, {1160., 1160., 1340.}, {1340., 1340., 1540.},}, //d-dd, d-du, d-ds, d-ud, d-uu, d-us, d-sd, d-su, d-ss
1297  {{1160., 1160., 1340.}, {1160., 1160., 1340.}, {1340., 1340., 1540.},}, //u-dd, u-du, u-ds, u-ud, u-uu, u-us, u-sd, u-su, u-ss
1298  {{1520., 1520., 1690.}, {1520., 1520., 1690.}, {1690., 1690., 1890. }}};//s-dd, s-du, s-ds, s-ud, s-uu, s-us, s-sd, s-su, s-ss
1299  for(G4int i=0; i<3; i++){ for(G4int j=0; j<3; j++){ for(G4int k=0; k<3; k++){minMassQDiQStr[i][j][k]=minMQDiQStr[i][j][k];};};};
1300 
1301 // ------ An estimated minimal string mass ----------------------
1302  MinimalStringMass = 0.;
1303  MinimalStringMass2 = 0.;
1304 
1305 // For treating of small string decays
1306  for(G4int i=0; i<3; i++)
1307  { for(G4int j=0; j<3; j++)
1308  { for(G4int k=0; k<6; k++)
1309  {
1310  Meson[i][j][k]=0; MesonWeight[i][j][k]=0.;
1311  }
1312  }
1313  }
1314 //--------------------------
1315  Meson[0][0][0]=111; // dbar-d Pi0
1316  MesonWeight[0][0][0]=(1.-pspin_meson)*(1.-scalarMesonMix[0]);
1317 
1318  Meson[0][0][1]=221; // dbar-d Eta
1319  MesonWeight[0][0][1]=(1.-pspin_meson)*(scalarMesonMix[0]-scalarMesonMix[1]);
1320 
1321  Meson[0][0][2]=331; // dbar-d EtaPrime
1322  MesonWeight[0][0][2]=(1.-pspin_meson)*(scalarMesonMix[1]);
1323 
1324  Meson[0][0][3]=113; // dbar-d Rho0
1325  MesonWeight[0][0][3]=pspin_meson*(1.-vectorMesonMix[0]);
1326 
1327  Meson[0][0][4]=223; // dbar-d Omega
1328  MesonWeight[0][0][4]=pspin_meson*(vectorMesonMix[0]);
1329 //--------------------------
1330 
1331  Meson[0][1][0]=211; // dbar-u Pi+
1332  MesonWeight[0][1][0]=(1.-pspin_meson);
1333 
1334  Meson[0][1][1]=213; // dbar-u Rho+
1335  MesonWeight[0][1][1]=pspin_meson;
1336 //--------------------------
1337 
1338  Meson[0][2][0]=311; // dbar-s K0bar
1339  MesonWeight[0][2][0]=(1.-pspin_meson);
1340 
1341  Meson[0][2][1]=313; // dbar-s K*0bar
1342  MesonWeight[0][2][1]=pspin_meson;
1343 //--------------------------
1344 //--------------------------
1345  Meson[1][0][0]=211; // ubar-d Pi-
1346  MesonWeight[1][0][0]=(1.-pspin_meson);
1347 
1348  Meson[1][0][1]=213; // ubar-d Rho-
1349  MesonWeight[1][0][1]=pspin_meson;
1350 //--------------------------
1351 
1352  Meson[1][1][0]=111; // ubar-u Pi0
1353  MesonWeight[1][1][0]=(1.-pspin_meson)*(1.-scalarMesonMix[0]);
1354 
1355  Meson[1][1][1]=221; // ubar-u Eta
1356  MesonWeight[1][1][1]=(1.-pspin_meson)*(scalarMesonMix[0]-scalarMesonMix[1]);
1357 
1358  Meson[1][1][2]=331; // ubar-u EtaPrime
1359  MesonWeight[1][1][2]=(1.-pspin_meson)*(scalarMesonMix[1]);
1360 
1361  Meson[1][1][3]=113; // ubar-u Rho0
1362  MesonWeight[1][1][3]=pspin_meson*(1.-vectorMesonMix[0]);
1363 
1364  Meson[1][1][4]=223; // ubar-u Omega
1365  MesonWeight[1][1][4]=pspin_meson*(vectorMesonMix[0]);
1366 //--------------------------
1367 
1368  Meson[1][2][0]=321; // ubar-s K-
1369  MesonWeight[1][2][0]=(1.-pspin_meson);
1370 
1371  Meson[1][2][1]=323; // ubar-s K*-bar -
1372  MesonWeight[1][2][1]=pspin_meson;
1373 //--------------------------
1374 //--------------------------
1375 
1376  Meson[2][0][0]=311; // sbar-d K0
1377  MesonWeight[2][0][0]=(1.-pspin_meson);
1378 
1379  Meson[2][0][1]=313; // sbar-d K*0
1380  MesonWeight[2][0][1]=pspin_meson;
1381 //--------------------------
1382 
1383  Meson[2][1][0]=321; // sbar-u K+
1384  MesonWeight[2][1][0]=(1.-pspin_meson);
1385 
1386  Meson[2][1][1]=323; // sbar-u K*+
1387  MesonWeight[2][1][1]=pspin_meson;
1388 //--------------------------
1389 
1390  Meson[2][2][0]=221; // sbar-s Eta
1391  MesonWeight[2][2][0]=(1.-pspin_meson)*(1.-scalarMesonMix[5]);
1392 
1393  Meson[2][2][1]=331; // sbar-s EtaPrime
1394  MesonWeight[2][2][1]=(1.-pspin_meson)*(1.-scalarMesonMix[5]);
1395 
1396  Meson[2][2][3]=333; // sbar-s EtaPrime
1397  MesonWeight[2][2][3]=pspin_meson*(vectorMesonMix[5]);
1398 //--------------------------
1399 
1400  for(G4int i=0; i<3; i++)
1401  { for(G4int j=0; j<3; j++)
1402  { for(G4int k=0; k<3; k++)
1403  { for(G4int l=0; l<4; l++)
1404  { Baryon[i][j][k][l]=0; BaryonWeight[i][j][k][l]=0.;}
1405  }
1406  }
1407  }
1408 
1409  G4double pspin_barion_in=pspin_barion;
1410  //pspin_barion=0.75;
1411 //---------------------------------------
1412  Baryon[0][0][0][0]=1114; // Delta-
1413  BaryonWeight[0][0][0][0]=1.;
1414 
1415 //---------------------------------------
1416  Baryon[0][0][1][0]=2112; // neutron
1417  BaryonWeight[0][0][1][0]=1.-pspin_barion;
1418 
1419  Baryon[0][0][1][1]=2114; // Delta0
1420  BaryonWeight[0][0][1][1]=pspin_barion;
1421 
1422 //---------------------------------------
1423  Baryon[0][0][2][0]=3112; // Sigma-
1424  BaryonWeight[0][0][2][0]=1.-pspin_barion;
1425 
1426  Baryon[0][0][2][1]=3114; // Sigma*-
1427  BaryonWeight[0][0][2][1]=pspin_barion;
1428 
1429 //---------------------------------------
1430  Baryon[0][1][0][0]=2112; // neutron
1431  BaryonWeight[0][1][0][0]=1.-pspin_barion;
1432 
1433  Baryon[0][1][0][1]=2114; // Delta0
1434  BaryonWeight[0][1][0][1]=pspin_barion;
1435 
1436 //---------------------------------------
1437  Baryon[0][1][1][0]=2212; // proton
1438  BaryonWeight[0][1][1][0]=1.-pspin_barion;
1439 
1440  Baryon[0][1][1][1]=2214; // Delta+
1441  BaryonWeight[0][1][1][1]=pspin_barion;
1442 
1443 //---------------------------------------
1444  Baryon[0][1][2][0]=3122; // Lambda
1445  BaryonWeight[0][1][2][0]=(1.-pspin_barion)*0.5;
1446 
1447  Baryon[0][1][2][1]=3212; // Sigma0
1448  BaryonWeight[0][1][2][1]=(1.-pspin_barion)*0.5;
1449 
1450  Baryon[0][1][2][2]=3214; // Sigma*0
1451  BaryonWeight[0][1][2][2]=pspin_barion;
1452 
1453 //---------------------------------------
1454  Baryon[0][2][0][0]=3112; // Sigma-
1455  BaryonWeight[0][2][0][0]=1.-pspin_barion;
1456 
1457  Baryon[0][2][0][1]=3114; // Sigma*-
1458  BaryonWeight[0][2][0][1]=pspin_barion;
1459 
1460 //---------------------------------------
1461  Baryon[0][2][1][0]=3122; // Lambda
1462  BaryonWeight[0][2][1][0]=(1.-pspin_barion)*0.5;
1463 
1464  Baryon[0][2][1][1]=3212; // Sigma0
1465  BaryonWeight[0][2][1][1]=(1.-pspin_barion)*0.5;
1466 
1467  Baryon[0][2][1][2]=3214; // Sigma*0
1468  BaryonWeight[0][2][1][2]=pspin_barion;
1469 
1470 //---------------------------------------
1471  Baryon[0][2][2][0]=3312; // Theta-
1472  BaryonWeight[0][2][2][0]=1.-pspin_barion;
1473 
1474  Baryon[0][2][2][1]=3314; // Theta*-
1475  BaryonWeight[0][2][2][1]=pspin_barion;
1476 
1477 //---------------------------------------
1478 //---------------------------------------
1479  Baryon[1][0][0][0]=2112; // neutron
1480  BaryonWeight[1][0][0][0]=1.-pspin_barion;
1481 
1482  Baryon[1][0][0][1]=2114; // Delta0
1483  BaryonWeight[1][0][0][1]=pspin_barion;
1484 
1485 //---------------------------------------
1486  Baryon[1][0][1][0]=2212; // proton
1487  BaryonWeight[1][0][1][0]=1.-pspin_barion;
1488 
1489  Baryon[1][0][1][1]=2214; // Delta+
1490  BaryonWeight[1][0][1][1]=pspin_barion;
1491 
1492 //---------------------------------------
1493  Baryon[1][0][2][0]=3122; // Lambda
1494  BaryonWeight[1][0][2][0]=(1.-pspin_barion)*0.5;
1495 
1496  Baryon[1][0][2][1]=3212; // Sigma0
1497  BaryonWeight[1][0][2][1]=(1.-pspin_barion)*0.5;
1498 
1499  Baryon[1][0][2][2]=3214; // Sigma*0
1500  BaryonWeight[1][0][2][2]=pspin_barion;
1501 
1502 //---------------------------------------
1503  Baryon[1][1][0][0]=2212; // proton
1504  BaryonWeight[1][1][0][0]=1.-pspin_barion;
1505 
1506  Baryon[1][1][0][1]=2214; // Delta+
1507  BaryonWeight[1][1][0][1]=pspin_barion;
1508 
1509 //---------------------------------------
1510  Baryon[1][1][1][0]=2224; // Delta++
1511  BaryonWeight[1][1][1][0]=1.;
1512 
1513 //---------------------------------------
1514  Baryon[1][1][2][0]=3222; // Sigma+
1515  BaryonWeight[1][1][2][0]=1.-pspin_barion;
1516 
1517  Baryon[1][1][2][1]=3224; // Sigma*+
1518  BaryonWeight[1][1][2][1]=pspin_barion;
1519 
1520 //---------------------------------------
1521  Baryon[1][2][0][0]=3122; // Lambda
1522  BaryonWeight[1][2][0][0]=(1.-pspin_barion)*0.5;
1523 
1524  Baryon[1][2][0][1]=3212; // Sigma0
1525  BaryonWeight[1][2][0][1]=(1.-pspin_barion)*0.5;
1526 
1527  Baryon[1][2][0][2]=3214; // Sigma*0
1528  BaryonWeight[1][2][0][2]=pspin_barion;
1529 
1530 //---------------------------------------
1531  Baryon[1][2][1][0]=3222; // Sigma+
1532  BaryonWeight[1][2][1][0]=1.-pspin_barion;
1533 
1534  Baryon[1][2][1][1]=3224; // Sigma*+
1535  BaryonWeight[1][2][1][1]=pspin_barion;
1536 
1537 //---------------------------------------
1538  Baryon[1][2][2][0]=3322; // Theta0
1539  BaryonWeight[1][2][2][0]=1.-pspin_barion;
1540 
1541  Baryon[1][2][2][1]=3324; // Theta*0
1542  BaryonWeight[1][2][2][1]=pspin_barion;
1543 
1544 //---------------------------------------
1545 //---------------------------------------
1546  Baryon[2][0][0][0]=3112; // Sigma-
1547  BaryonWeight[2][0][0][0]=1.-pspin_barion;
1548 
1549  Baryon[2][0][0][1]=3114; // Sigma*-
1550  BaryonWeight[2][0][0][1]=pspin_barion;
1551 
1552 //---------------------------------------
1553  Baryon[2][0][1][0]=3122; // Lambda
1554  BaryonWeight[2][0][1][0]=(1.-pspin_barion)*0.5;
1555 
1556  Baryon[2][0][1][1]=3212; // Sigma0
1557  BaryonWeight[2][0][1][1]=(1.-pspin_barion)*0.5;
1558 
1559  Baryon[2][0][1][2]=3214; // Sigma*0
1560  BaryonWeight[2][0][1][2]=pspin_barion;
1561 
1562 //---------------------------------------
1563  Baryon[2][0][2][0]=3312; // Sigma-
1564  BaryonWeight[2][0][2][0]=1.-pspin_barion;
1565 
1566  Baryon[2][0][2][1]=3314; // Sigma*-
1567  BaryonWeight[2][0][2][1]=pspin_barion;
1568 
1569 //---------------------------------------
1570  Baryon[2][1][0][0]=3122; // Lambda
1571  BaryonWeight[2][1][0][0]=(1.-pspin_barion)*0.5;
1572 
1573  Baryon[2][1][0][1]=3212; // Sigma0
1574  BaryonWeight[2][1][0][1]=(1.-pspin_barion)*0.5;
1575 
1576  Baryon[2][1][0][2]=3214; // Sigma*0
1577  BaryonWeight[2][1][0][2]=pspin_barion;
1578 
1579 //---------------------------------------
1580  Baryon[2][1][1][0]=3222; // Sigma+
1581  BaryonWeight[2][1][1][0]=1.-pspin_barion;
1582 
1583  Baryon[2][1][1][1]=3224; // Sigma*+
1584  BaryonWeight[2][1][1][1]=pspin_barion;
1585 
1586 //---------------------------------------
1587  Baryon[2][1][2][0]=3322; // Theta0
1588  BaryonWeight[2][1][2][0]=1.-pspin_barion;
1589 
1590  Baryon[2][1][2][1]=3324; // Theta*0
1591  BaryonWeight[2][1][2][1]=pspin_barion;
1592 
1593 //---------------------------------------
1594  Baryon[2][2][0][0]=3312; // Theta-
1595  BaryonWeight[2][2][0][0]=1.-pspin_barion;
1596 
1597  Baryon[2][2][0][1]=3314; // Theta*-
1598  BaryonWeight[2][2][0][1]=pspin_barion;
1599 
1600 //---------------------------------------
1601  Baryon[2][2][1][0]=3322; // Theta0
1602  BaryonWeight[2][2][1][0]=1.-pspin_barion;
1603 
1604  Baryon[2][2][1][1]=3324; // Theta*0
1605  BaryonWeight[2][2][1][1]=pspin_barion;
1606 
1607 //---------------------------------------
1608  Baryon[2][2][2][0]=3334; // Omega
1609  BaryonWeight[2][2][2][0]=1.;
1610 
1611 //---------------------------------------
1612  pspin_barion=pspin_barion_in;
1613  /*
1614  for(G4int i=0; i<3; i++)
1615  { for(G4int j=0; j<3; j++)
1616  { for(G4int k=0; k<3; k++)
1617  { for(G4int l=0; l<4; l++)
1618  { G4cout<<i<<" "<<j<<" "<<k<<" "<<l<<" "<<Baryon[i][j][k][l]<<G4endl;}
1619  }
1620  }
1621  }
1622  G4int Uzhi;
1623  G4cin>>Uzhi;
1624  */
1625 
1626  G4double ProbUUbar = 0.33; // 0.40 -> 0.33 Uzhi Nov. 2017
1627  Prob_QQbar[0]=ProbUUbar; // Probability of ddbar production
1628  Prob_QQbar[1]=ProbUUbar; // Probability of uubar production
1629  Prob_QQbar[2]=1.0-2.*ProbUUbar; // Probability of ssbar production
1630 
1631  for ( G4int i=0 ; i<35 ; i++ ) {
1632  FS_LeftHadron[i] = 0;
1633  FS_RightHadron[i] = 0;
1634  FS_Weight[i] = 0.0;
1635  }
1636 
1637  NumberOf_FS = 0;
1638 
1639 }
1640 
1641 // --------------------------------------------------------------
1642 G4LundStringFragmentation::~G4LundStringFragmentation()
1643 {}
1644 
1645 //--------------------------------------------------------------------------------------
1646 void G4LundStringFragmentation::SetMinimalStringMass(const G4FragmentingString * const string)
1647 {
1648  G4double EstimatedMass=0.;
1649 
1650  G4int Qleft =std::abs(string->GetLeftParton()->GetPDGEncoding());
1651  G4int Qright=std::abs(string->GetRightParton()->GetPDGEncoding());
1652 
1653  if((Qleft < 4) && (Qright < 4)) { // Q-Qbar string
1654  EstimatedMass=minMassQQbarStr[Qleft-1][Qright-1];
1655  MinimalStringMass=EstimatedMass;
1656  SetMinimalStringMass2(EstimatedMass);
1657  return;
1658  }
1659 
1660  if((Qleft < 4) && (Qright > 1000)) { // Q - DiQ string
1661  G4int q1=Qright/1000;
1662  G4int q2=(Qright/100)%10;
1663  EstimatedMass=minMassQDiQStr[Qleft-1][q1-1][q2-1];
1664  MinimalStringMass=EstimatedMass;
1665  SetMinimalStringMass2(EstimatedMass);
1666  return;
1667  }
1668 
1669  if((Qleft > 1000) && (Qright < 4)) { // DiQ - Q string
1670  G4int q1=Qleft/1000;
1671  G4int q2=(Qleft/100)%10;
1672  EstimatedMass=minMassQDiQStr[Qright-1][q1-1][q2-1];
1673  MinimalStringMass=EstimatedMass;
1674  SetMinimalStringMass2(EstimatedMass);
1675  return;
1676  }
1677 
1678 // DiQuark - Anti DiQuark string -----------------
1679  G4int Number_of_quarks=0;
1680  G4int Number_of_squarks=0;
1681 
1682  G4double StringM=string->Get4Momentum().mag();
1683 
1684 #ifdef debug_LUNDfragmentation
1685 // G4cout<<"MinStringMass// Input String mass "<<string->Get4Momentum().mag()<<" Qleft "<<Qleft<<G4endl;
1686 #endif
1687 
1688  if( Qleft > 1000)
1689  {
1690  Number_of_quarks+=2;
1691  G4int q1=Qleft/1000;
1692  if( q1 < 3) {EstimatedMass +=Mass_of_light_quark;}
1693  if( q1 > 2) {EstimatedMass +=Mass_of_heavy_quark; Number_of_squarks++;}
1694 
1695  G4int q2=(Qleft/100)%10;
1696  if( q2 < 3) {EstimatedMass +=Mass_of_light_quark;}
1697  if( q2 > 2) {EstimatedMass +=Mass_of_heavy_quark; Number_of_squarks++;}
1698 // EstimatedMass +=Mass_of_string_junction;
1699  }
1700 
1701 #ifdef debug_LUNDfragmentation
1702 // G4cout<<"Min mass with Qleft "<<Qleft<<" "<<EstimatedMass<<G4endl;
1703 #endif
1704 // G4int Qright=std::abs(string->GetRightParton()->GetPDGEncoding());
1705 
1706  if( Qright > 1000)
1707  {
1708  Number_of_quarks+=2;
1709  G4int q1=Qright/1000;
1710  if( q1 < 3) {EstimatedMass +=Mass_of_light_quark;}
1711  if( q1 > 2) {EstimatedMass +=Mass_of_heavy_quark; Number_of_squarks++;}
1712 
1713  G4int q2=(Qright/100)%10;
1714  if( q2 < 3) {EstimatedMass +=Mass_of_light_quark;}
1715  if( q2 > 2) {EstimatedMass +=Mass_of_heavy_quark; Number_of_squarks++;}
1716  //EstimatedMass +=Mass_of_string_junction;
1717  }
1718 
1719 #ifdef debug_LUNDfragmentation
1720 // G4cout<<"Min mass with Qleft and Qright "<<Qright<<" "<<EstimatedMass<<G4endl;
1721 // G4cout<<"Number_of_quarks "<<Number_of_quarks<<" Number_of_squarks "<<Number_of_squarks<<G4endl;
1722 #endif
1723 
1724  if(Number_of_quarks==4)
1725  {
1726  if(StringM > 1880.) { // 2*Mn = 1880
1727  if(Number_of_squarks==0) {EstimatedMass += 1320.*MeV;}//560+1320=1880=2*Mn
1728  else if(Number_of_squarks==1) {EstimatedMass += 1150.*MeV;}//920+1150=2070=M(Lam+N)
1729  else if(Number_of_squarks==2) {EstimatedMass += 960.*MeV;}//1280+960=2240= 2*M Lam
1730  else if(Number_of_squarks==3) {EstimatedMass += 800.*MeV;}//1640+800=2440=Mxi+Mlam
1731  else if(Number_of_squarks==4) {EstimatedMass += 640.*MeV;}//2000+640=2640=2*Mxi
1732  else {}
1733  }
1734 /*
1735  if((StringM > 1880.) && ( EstimatedMass < 2100)) {EstimatedMass = 2020.;}//1880.;}
1736  // if((StringM > 1880.) && ( EstimatedMass < 2100)) {EstimatedMass = 2051.;}
1737  else if((StringM > 2232.) && ( EstimatedMass < 2730)){EstimatedMass = 2570.;}
1738  else if((StringM > 5130.) && ( EstimatedMass < 3450)){EstimatedMass = 5130.;}
1739 */
1740  else
1741  {
1742  if(Number_of_squarks < 3) {EstimatedMass -= 200.*MeV;}
1743  else if(Number_of_squarks==3) {EstimatedMass -= 50.*MeV;}
1744  else if(Number_of_squarks==4) {EstimatedMass -= 40.*MeV;}
1745  else {}
1746  }
1747  }
1748 
1749 #ifdef debug_LUNDfragmentation
1750 // G4cout<<"EstimatedMass -------------------- "<<EstimatedMass <<G4endl;
1751 #endif
1752 
1753  MinimalStringMass=EstimatedMass;
1754  SetMinimalStringMass2(EstimatedMass);
1755 }
1756 
1757 //--------------------------------------------------------------------------------------
1758 void G4LundStringFragmentation::SetMinimalStringMass2(const G4double aValue)
1759 {
1760  MinimalStringMass2=aValue * aValue;
1761 }
G4ParticleTable::FindParticle
G4ParticleDefinition * FindParticle(G4int PDGEncoding)
Definition: G4ParticleTable.cc:577
G4FragmentingString::GetPleft
G4LorentzVector GetPleft()
Definition: G4FragmentingString.cc:326
G4LundStringFragmentation::Quark_Diquark_lastSplitting
G4bool Quark_Diquark_lastSplitting(G4FragmentingString *&string, G4ParticleDefinition *&LeftHadron, G4ParticleDefinition *&RightHadron)
Definition: G4LundStringFragmentation.cc:947
G4Exp
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
G4FragmentingString::GetPright
G4LorentzVector GetPright()
Definition: G4FragmentingString.cc:329
G4INCL::ClusterDecay::decay
ParticleList decay(Cluster *const c)
Carries out a cluster decay.
Definition: G4INCLClusterDecay.cc:553
G4INCL::Math::max
T max(const T t1, const T t2)
brief Return the largest of the two arguments
Definition: G4INCLGlobals.hh:112
G4FragmentingString::Mass
G4double Mass() const
Definition: G4FragmentingString.cc:312
G4ExcitedString::GetRightParton
G4Parton * GetRightParton(void) const
Definition: G4ExcitedString.cc:168
G4FragmentingString::LorentzRotate
void LorentzRotate(const G4LorentzRotation &rotation)
Definition: G4FragmentingString.hh:162
G4LundStringFragmentation::FS_LeftHadron
G4ParticleDefinition * FS_LeftHadron[35]
Definition: G4LundStringFragmentation.hh:140
G4HadronBuilder::Build
G4ParticleDefinition * Build(G4ParticleDefinition *black, G4ParticleDefinition *white)
Definition: G4HadronBuilder.cc:53
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:42
G4ExcitedString::GetPosition
const G4ThreeVector & GetPosition() const
Definition: G4ExcitedString.hh:149
G4ExceptionDescription
std::ostringstream G4ExceptionDescription
Definition: G4Exception.hh:45
G4KineticTrack::Get4Momentum
const G4LorentzVector & Get4Momentum() const
Definition: G4KineticTrack.hh:240
G4LundStringFragmentation::DiQuarkSplitup
virtual G4ParticleDefinition * DiQuarkSplitup(G4ParticleDefinition *decay, G4ParticleDefinition *&created)
Definition: G4LundStringFragmentation.cc:443
G4ParticleTable::GetParticleTable
static G4ParticleTable * GetParticleTable()
Definition: G4ParticleTable.cc:105
MeV
static constexpr double MeV
Definition: G4SIunits.hh:214
G4VLongitudinalStringDecay::SetDiquarkSuppression
void SetDiquarkSuppression(G4double aValue)
Definition: G4VLongitudinalStringDecay.cc:410
G4FragmentingString::GetRightParton
G4ParticleDefinition * GetRightParton(void) const
Definition: G4FragmentingString.hh:155
G4VLongitudinalStringDecay::SetMassCut
virtual void SetMassCut(G4double aValue)
Definition: G4VLongitudinalStringDecay.cc:143
G4ParticleDefinition::GetParticleSubType
const G4String & GetParticleSubType() const
Definition: G4ParticleDefinition.hh:143
G4endl
#define G4endl
Definition: G4ios.hh:61
Pos
ush Pos
Definition: deflate.h:92
G4VLongitudinalStringDecay::SetStringTensionParameter
void SetStringTensionParameter(G4double aValue)
Definition: G4VLongitudinalStringDecay.cc:502
z
Double_t z
Definition: advanced/microbeam/plot.C:277
G4LundStringFragmentation::FragmentString
virtual G4KineticTrackVector * FragmentString(const G4ExcitedString &theString)
Definition: G4LundStringFragmentation.cc:87
G4ParticleDefinition::GetParticleName
const G4String & GetParticleName() const
Definition: G4ParticleDefinition.hh:121
G4KineticTrack::GetDefinition
const G4ParticleDefinition * GetDefinition() const
Definition: G4KineticTrack.hh:210
CLHEP::HepLorentzVector::getX
double getX() const
G4VLongitudinalStringDecay::SetVectorMesonMixings
void SetVectorMesonMixings(std::vector< G4double > aVector)
Definition: G4VLongitudinalStringDecay.cc:481
tmp
Float_t tmp
Definition: advanced/microbeam/plot.C:36
G4ParticleDefinition::GetPDGMass
G4double GetPDGMass() const
Definition: G4ParticleDefinition.hh:123
G4VLongitudinalStringDecay::SetSpinThreeHalfBarionProbability
void SetSpinThreeHalfBarionProbability(G4double aValue)
Definition: G4VLongitudinalStringDecay.cc:446
G4VLongitudinalStringDecay::SetDiquarkBreakProbability
void SetDiquarkBreakProbability(G4double aValue)
Definition: G4VLongitudinalStringDecay.cc:421
G4Log
G4double G4Log(G4double x)
Definition: G4Log.hh:230
S
double S(double temp)
Definition: G4DNAElectronHoleRecombination.cc:77
G4VLongitudinalStringDecay::vectorMesonMix
std::vector< G4double > vectorMesonMix
Definition: G4VLongitudinalStringDecay.hh:190
G4LundStringFragmentation::lambda
G4double lambda(G4double s, G4double m1_Sqr, G4double m2_Sqr)
Definition: G4LundStringFragmentation.cc:1276
G4LundStringFragmentation::FS_RightHadron
G4ParticleDefinition * FS_RightHadron[35]
Definition: G4LundStringFragmentation.hh:140
G4ExcitedString::GetDirection
G4int GetDirection(void) const
Definition: G4ExcitedString.cc:154
G4double
double G4double
Definition: G4Types.hh:76
G4bool
bool G4bool
Definition: G4Types.hh:79
fermi
static constexpr double fermi
Definition: G4SIunits.hh:103
G4VLongitudinalStringDecay::SetVectorMesonProbability
void SetVectorMesonProbability(G4double aValue)
Definition: G4VLongitudinalStringDecay.cc:432
G4Pow::GetInstance
static G4Pow * GetInstance()
Definition: G4Pow.cc:57
G4VLongitudinalStringDecay::FindParticle
G4ParticleDefinition * FindParticle(G4int Encoding)
Definition: G4VLongitudinalStringDecay.cc:262
G4LundStringFragmentation::Splitup
virtual G4KineticTrack * Splitup(G4FragmentingString *string, G4FragmentingString *&newString)
Definition: G4LundStringFragmentation.cc:305
G4FragmentingString::GetDecayDirection
G4int GetDecayDirection() const
Definition: G4FragmentingString.cc:231
G4VLongitudinalStringDecay::SampleQuarkPt
G4ThreeVector SampleQuarkPt(G4double ptMax=-1.)
Definition: G4VLongitudinalStringDecay.cc:342
G4Pow::powA
G4double powA(G4double A, G4double y) const
Definition: G4Pow.hh:242
CLHEP::Hep3Vector::setZ
void setZ(double)
G4VLongitudinalStringDecay::LightFragmentationTest
G4KineticTrackVector * LightFragmentationTest(const G4ExcitedString *const theString)
Definition: G4VLongitudinalStringDecay.cc:149
G4ExcitedString::GetTimeOfCreation
G4double GetTimeOfCreation() const
Definition: G4ExcitedString.hh:137
G4FragmentingString::TransformToAlignedCms
G4LorentzRotation TransformToAlignedCms()
Definition: G4FragmentingString.hh:190
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:53
G4VLongitudinalStringDecay::SetStrangenessSuppression
void SetStrangenessSuppression(G4double aValue)
Definition: G4VLongitudinalStringDecay.cc:399
CLHEP::Hep3Vector::getX
double getX() const
G4ExcitedString::Get4Momentum
G4LorentzVector Get4Momentum() const
Definition: G4ExcitedString.hh:161
CLHEP::HepLorentzRotation::inverse
HepLorentzRotation inverse() const
G4ExcitedString::GetLeftParton
G4Parton * GetLeftParton(void) const
Definition: G4ExcitedString.cc:161
CLHEP::Hep3Vector::mag2
double mag2() const
CLHEP::Hep3Vector::getY
double getY() const
G4LundStringFragmentation::SplitEandP
virtual G4LorentzVector * SplitEandP(G4ParticleDefinition *pHadron, G4FragmentingString *string, G4FragmentingString *newString)
Definition: G4LundStringFragmentation.cc:496
G4LundStringFragmentation::SetMinimalStringMass
void SetMinimalStringMass(const G4FragmentingString *const string)
Definition: G4LundStringFragmentation.cc:1646
G4LundStringFragmentation::Diquark_AntiDiquark_belowThreshold_lastSplitting
G4bool Diquark_AntiDiquark_belowThreshold_lastSplitting(G4FragmentingString *&string, G4ParticleDefinition *&LeftHadron, G4ParticleDefinition *&RightHadron)
Definition: G4LundStringFragmentation.cc:805
G4VLongitudinalStringDecay::QuarkSplitup
G4ParticleDefinition * QuarkSplitup(G4ParticleDefinition *decay, G4ParticleDefinition *&created)
Definition: G4VLongitudinalStringDecay.cc:298
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.hh:65
G4LundStringFragmentation::Sample4Momentum
virtual void Sample4Momentum(G4LorentzVector *Mom, G4double Mass, G4LorentzVector *AntiMom, G4double AntiMass, G4double InitialMass)
Definition: G4LundStringFragmentation.cc:1195
G4Parton::GetPDGcode
G4int GetPDGcode() const
Definition: G4Parton.hh:127
CLHEP::HepLorentzVector::mag
double mag() const
G4KineticTrack::GetFormationTime
G4double GetFormationTime() const
Definition: G4KineticTrack.hh:220
G4int
int G4int
Definition: G4Types.hh:78
CLHEP::c_light
static constexpr double c_light
Definition: PhysicalConstants.h:55
G4FragmentingString::GetLeftParton
G4ParticleDefinition * GetLeftParton(void) const
Definition: G4FragmentingString.hh:149
G4KineticTrack::SetPosition
void SetPosition(const G4ThreeVector aPosition)
Definition: G4KineticTrack.hh:235
G4FragmentingString::FourQuarkString
G4bool FourQuarkString(void) const
Definition: G4FragmentingString.cc:241
G4VLongitudinalStringDecay::CreatePartonPair
pDefPair CreatePartonPair(G4int NeedParticle, G4bool AllowDiquarks=true)
Definition: G4VLongitudinalStringDecay.cc:319
CLHEP::HepLorentzVector::getY
double getY() const
G4FragmentingString::GetDecayParton
G4ParticleDefinition * GetDecayParton() const
Definition: G4FragmentingString.hh:143
G4cout
G4GLOB_DLL std::ostream G4cout
CLHEP::Hep3Vector::x
double x() const
G4LundStringFragmentation::GetLightConeZ
virtual G4double GetLightConeZ(G4double zmin, G4double zmax, G4int PartonEncoding, G4ParticleDefinition *pHadron, G4double Px, G4double Py)
Definition: G4LundStringFragmentation.cc:595
G4LundStringFragmentation::Loop_toFragmentString
G4bool Loop_toFragmentString(const G4ExcitedString &theStringInCMS, G4KineticTrackVector *&LeftVector, G4KineticTrackVector *&RightVector)
Definition: G4LundStringFragmentation.cc:176
CLHEP::HepLorentzVector::vect
Hep3Vector vect() const
G4LundStringFragmentation::SplitLast
virtual G4bool SplitLast(G4FragmentingString *string, G4KineticTrackVector *LeftVector, G4KineticTrackVector *RightVector)
Definition: G4LundStringFragmentation.cc:657
sqr
T sqr(const T &x)
Definition: templates.hh:145
pi
static constexpr double pi
Definition: G4SIunits.hh:75
G4KineticTrack::GetPosition
const G4ThreeVector & GetPosition() const
Definition: G4KineticTrack.hh:230
G4LorentzVector
CLHEP::HepLorentzVector G4LorentzVector
Definition: G4LorentzVector.hh:35
G4VLongitudinalStringDecay::pDefPair
std::pair< G4ParticleDefinition *, G4ParticleDefinition * > pDefPair
Definition: G4VLongitudinalStringDecay.hh:72
CLHEP::Hep3Vector::y
double y() const
G4LundStringFragmentation::IsFragmentable
virtual G4bool IsFragmentable(const G4FragmentingString *const string)
Definition: G4LundStringFragmentation.cc:168
CLHEP::HepLorentzVector::getT
double getT() const
G4LundStringFragmentation::Diquark_AntiDiquark_aboveThreshold_lastSplitting
G4bool Diquark_AntiDiquark_aboveThreshold_lastSplitting(G4FragmentingString *&string, G4ParticleDefinition *&LeftHadron, G4ParticleDefinition *&RightHadron)
Definition: G4LundStringFragmentation.cc:849
G4VLongitudinalStringDecay::scalarMesonMix
std::vector< G4double > scalarMesonMix
Definition: G4VLongitudinalStringDecay.hh:191
G4Parton::Get4Momentum
const G4LorentzVector & Get4Momentum() const
Definition: G4Parton.hh:143
G4KineticTrack::Set4Momentum
void Set4Momentum(const G4LorentzVector &a4Momentum)
Definition: G4KineticTrack.hh:250
G4LundStringFragmentation::Quark_AntiQuark_lastSplitting
G4bool Quark_AntiQuark_lastSplitting(G4FragmentingString *&string, G4ParticleDefinition *&LeftHadron, G4ParticleDefinition *&RightHadron)
Definition: G4LundStringFragmentation.cc:1059
G4KineticTrack::SetFormationTime
void SetFormationTime(G4double aFormationTime)
Definition: G4KineticTrack.hh:225
GeV
static constexpr double GeV
Definition: G4SIunits.hh:217
CLHEP::HepLorentzVector::getZ
double getZ() const
G4FragmentingString::GetPstring
G4LorentzVector GetPstring()
Definition: G4FragmentingString.cc:323
G4INCL::Math::min
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
Definition: G4INCLGlobals.hh:117
G4LundStringFragmentation::StopFragmenting
virtual G4bool StopFragmenting(const G4FragmentingString *const string)
Definition: G4LundStringFragmentation.cc:286
G4LundStringFragmentation::SetMinimalStringMass2
void SetMinimalStringMass2(const G4double aValue)
Definition: G4LundStringFragmentation.cc:1758
다음에 의해 생성됨 :   doxygen 1.8.5