Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
G4QGSMFragmentation.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: G4QGSMFragmentation.cc 107968 2017-12-14 13:16:35Z gcosmo $
28 //
29 // -----------------------------------------------------------------------------
30 // GEANT 4 class implementation file
31 //
32 // History: first implementation, Maxim Komogorov, 10-Jul-1998
33 // -----------------------------------------------------------------------------
34 #include "G4QGSMFragmentation.hh"
35 #include "G4PhysicalConstants.hh"
36 #include "Randomize.hh"
37 #include "G4ios.hh"
38 #include "G4FragmentingString.hh"
39 #include "G4DiQuarks.hh"
40 #include "G4Quarks.hh"
41 
42 #include "G4Pow.hh"
43 
44 //#define debug_QGSMfragmentation
45 
46 // Class G4QGSMFragmentation
47 //****************************************************************************************
48 
50 arho(0.5), aphi(0.), an(-0.5), ala(-0.75), aksi(-1.), alft(0.5)
51  {
52  SetStrangenessSuppression((1.0 - 0.168)/2.);//It was (0.43); Uzhi Dec. 2017
53  SetDiquarkSuppression(0.299); // 0.087 std 0.07 Uzhi Dec. 2017
54  SetDiquarkBreakProbability(0.4); // 0.05 std 0.1 Uzhi Dec. 2017
55 
56  //... pspin_meson is probability to create pseudo-scalar meson
57  pspin_meson = 0.30; SetVectorMesonProbability(pspin_meson); // Uzhi Dec. 2017
58 
59  //... pspin_barion is probability to create 1/2 barion
61 
62  //... vectorMesonMix[] is quark mixing parameters for vector mesons (Variable spin = 3)
63  vectorMesonMix[0] = 0.; // 0.5 Uzhi Dec. 2017
64  vectorMesonMix[1] = 0.375;
65  vectorMesonMix[2] = 0.0;
66  vectorMesonMix[3] = 0.375; // 0.5 Uzhi Dec. 2017
67  vectorMesonMix[4] = 1.0;
68  vectorMesonMix[5] = 1.0;
69 
70  SetVectorMesonMixings(vectorMesonMix); // Uzhi Dec. 2017
71  }
72 
74  {
75  }
76 
77 //----------------------------------------------------------------------------------------------------------
78 
80 {
81 #ifdef debug_QGSMfragmentation
82  G4cout<<G4endl<<"QGSM StringFragm: String Mass "
83  <<theString.Get4Momentum().mag()<<" Pz "
84  <<theString.Get4Momentum().pz()
85  <<"------------------------------------"<<G4endl;
86  G4cout<<"String ends Direct "<<theString.GetLeftParton()->GetPDGcode()<<" "
87  <<theString.GetRightParton()->GetPDGcode()<<" "
88  <<theString.GetDirection()<< G4endl;
89  G4cout<<"Left mom "<<theString.GetLeftParton()->Get4Momentum()<<G4endl;
90  G4cout<<"Right mom "<<theString.GetRightParton()->Get4Momentum()<<G4endl;
91  G4cout<<"Check for Fragmentation "<<G4endl;
92 #endif
93 
94 // Can no longer modify Parameters for Fragmentation.
95  PastInitPhase=true;
96 
97 // check if string has enough mass to fragment...
98 
99  G4KineticTrackVector * LeftVector=LightFragmentationTest(&theString);
100 
101 #ifdef debug_QGSMfragmentation
102  if ( LeftVector != 0 ) G4cout<<"Non fragmentable - the string is converted to one hadron "<<G4endl;
103 #endif
104 
105  if ( LeftVector != 0 ) return LeftVector;
106 
107 #ifdef debug_QGSMfragmentation
108  G4cout<<"The string will be fragmented. "<<G4endl;
109 #endif
110 
111  LeftVector = new G4KineticTrackVector;
112  G4KineticTrackVector * RightVector=new G4KineticTrackVector;
113 
114 // this should work but its only a semi deep copy. %GF G4ExcitedString theStringInCMS(theString);
115  G4ExcitedString *theStringInCMS=CopyExcited(theString);
116  G4LorentzRotation toCms=theStringInCMS->TransformToAlignedCms();
117 
118  G4bool success=false, inner_sucess=true;
119  G4int attempt=0;
120  while ( !success && attempt++ < StringLoopInterrupt ) /* Loop checking, 07.08.2015, A.Ribon */
121  {
122 #ifdef debug_QGSMfragmentation
123  G4cout<<"Loop_toFrag "<<theStringInCMS->GetLeftParton()->GetPDGcode()<<" "
124  <<theStringInCMS->GetRightParton()->GetPDGcode()<<" "
125  <<theStringInCMS->GetDirection()<< G4endl;
126 #endif
127 
128  G4FragmentingString *currentString=new G4FragmentingString(*theStringInCMS);
129 
130  std::for_each(LeftVector->begin(), LeftVector->end(), DeleteKineticTrack());
131  LeftVector->clear();
132  std::for_each(RightVector->begin(), RightVector->end(), DeleteKineticTrack());
133  RightVector->clear();
134 
135  inner_sucess=true; // set false on failure..
136  const G4int maxNumberOfLoops = 1000;
137  G4int loopCounter = -1;
138  while (! StopFragmenting(currentString) && ++loopCounter < maxNumberOfLoops ) /* Loop checking, 07.08.2015, A.Ribon */
139  { // Split current string into hadron + new string
140 
141 #ifdef debug_QGSMfragmentation
142  G4cout<<"The string can fragment. "<<G4endl;;
143 #endif
144  G4FragmentingString *newString=0; // used as output from SplitUp...
145  G4KineticTrack * Hadron=Splitup(currentString,newString);
146 
147  if ( Hadron != 0 )
148  {
149 #ifdef debug_QGSMfragmentation
150  G4cout<<"Hadron prod at fragm. "<<Hadron->GetDefinition()->GetParticleName()<<G4endl;
151 #endif
152  if ( currentString->GetDecayDirection() > 0 )
153  LeftVector->push_back(Hadron);
154  else
155  RightVector->push_back(Hadron);
156 
157  delete currentString;
158  currentString=newString;
159 
160  } else {
161 
162 #ifdef debug_QGSMfragmentation
163  G4cout<<"abandon ... start from the beginning ---------------"<<G4endl;
164 #endif
165 
166  // abandon ... start from the beginning
167  if (newString) delete newString;
168  inner_sucess=false;
169  break;
170  }
171  }
172  if ( loopCounter >= maxNumberOfLoops ) {
173  inner_sucess=false;
174  }
175 
176  // Split current string into 2 final Hadrons
177 #ifdef debug_QGSMfragmentation
178  G4cout<<"Split remaining string into 2 final hadrons."<<G4endl;
179 #endif
180 
181  if ( inner_sucess &&
182  SplitLast(currentString,LeftVector, RightVector) )
183  {
184  success=true;
185  }
186  delete currentString;
187  }
188 
189  delete theStringInCMS;
190 
191  if ( ! success )
192  {
193  std::for_each(LeftVector->begin(), LeftVector->end(), DeleteKineticTrack());
194  LeftVector->clear();
195  std::for_each(RightVector->begin(), RightVector->end(), DeleteKineticTrack());
196  delete RightVector;
197  return LeftVector;
198  }
199 
200  // Join Left- and RightVector into LeftVector in correct order.
201  while(!RightVector->empty()) /* Loop checking, 07.08.2015, A.Ribon */
202  {
203  LeftVector->push_back(RightVector->back());
204  RightVector->erase(RightVector->end()-1);
205  }
206  delete RightVector;
207 
208  CalculateHadronTimePosition(theString.Get4Momentum().mag(), LeftVector);
209 
210  G4LorentzRotation toObserverFrame(toCms.inverse());
211 
212  for(size_t C1 = 0; C1 < LeftVector->size(); C1++)
213  {
214  G4KineticTrack* Hadron = LeftVector->operator[](C1);
215  G4LorentzVector Momentum = Hadron->Get4Momentum();
216  Momentum = toObserverFrame*Momentum;
217  Hadron->Set4Momentum(Momentum);
218  G4LorentzVector Coordinate(Hadron->GetPosition(), Hadron->GetFormationTime());
219  Momentum = toObserverFrame*Coordinate;
220  Hadron->SetFormationTime(Momentum.e());
221  G4ThreeVector aPosition(Momentum.vect());
222  Hadron->SetPosition(theString.GetPosition()+aPosition);
223  }
224  return LeftVector;
225 
226 
227 
228 }
229 
230 //----------------------------------------------------------------------------------------------------------
231 
234 {
235 #ifdef debug_QGSMfragmentation
236  G4cout<<"GetLightConeZ zmin zmax Parton pHadron "<<zmin<<" "<<zmax<<" "<<PartonEncoding<<" "<<pHadron->GetParticleName()<<G4endl;
237 #endif
238  G4double z;
239  G4double d1, d2, yf;
240  G4double invD1(0.),invD2(0.), r1(0.),r2(0.),r12(0.);
241 
242  G4int absCode = std::abs( PartonEncoding );
243  G4int absHadronCode=std::abs(pHadron->GetPDGEncoding());
244 
245  G4int q1, q2, q3;
246  q1 = absHadronCode/1000; q2 = (absHadronCode % 1000)/100; q3 = (absHadronCode % 100)/10;
247 
248  G4bool StrangeHadron = (q1 == 3) || (q2 == 3) || (q3 == 3);
249 
250  if (absCode < 10)
251  { // A quark fragmentation ----------------------------
252  if(absCode == 1 || absCode == 2)
253  {
254  if(absHadronCode < 1000)
255  { // Meson produced
256  if( !StrangeHadron ) {d1=2.0; d2 = -arho + alft;}
257  else {d1=1.0; d2 = -aphi + alft;}
258  } else
259  { // Baryon produced
260  if( !StrangeHadron ) {d1=0.0; d2 = arho - 2.0*an + alft;}
261  else {d1=0.0; d2 = 2.0*arho - 2.0*an - aphi + alft;}
262  }
263  }
264  else if(absCode == 3)
265  {
266  if(absHadronCode < 1000){d1=1.0 - aphi; d2 = -arho + alft;} // Meson produced s->K + u/d
267  else {d1=1.0 - aphi; d2 = arho - 2.0*an + alft;} // Baryon produced
268 
269  } else throw G4HadronicException(__FILE__, __LINE__, "Unknown PDGencoding in G4QGSMFragmentation::G4LightConeZ");
270 
271 #ifdef debug_QGSMfragmentation
272  G4cout<<"d1 d2 "<<d1<<" "<<d2<<G4endl;
273 #endif
274 
275  d1+=1.0; d2+=1.0;
276 
277  invD1=1./d1; invD2=1./d2;
278 
279  const G4int maxNumberOfLoops = 10000;
280  G4int loopCounter = 0;
281  do
282  {
283  r1=G4Pow::GetInstance()->powA(G4UniformRand(),invD1);
284  r2=G4Pow::GetInstance()->powA(G4UniformRand(),invD2);
285  r12=r1+r2;
286  z=r1/r12;
287  } while( ( (r12 > 1.0) || !((zmin <= z)&&(z <= zmax))) && ++loopCounter < maxNumberOfLoops ); /* Loop checking, 07.08.2015, A.Ribon */
288  if ( loopCounter >= maxNumberOfLoops ) {
289  z = 0.5*(zmin + zmax); // Just a value between zmin and zmax, no physics considerations at all!
290  }
291 
292  return z;
293  }
294  else
295  { // A di-quark fragmentation -------------------------
296  if(absCode == 1103 || absCode == 2101 ||
297  absCode == 2203 || absCode == 2103)
298  {
299  if(absHadronCode < 1000) // Meson production
300  {
301  if( !StrangeHadron ) {d1=1.0; d2= arho - 2.0*an + alft;}
302  else {d1=1.0; d2 = 2.*arho - 2.0*an - aphi + alft;}
303  } else // Baryon production
304  {
305  if( !StrangeHadron ) {d1=2.0*(arho - an); d2= -arho + alft;}
306  else {d1=2.0*(arho - an); d2 =-aphi + alft;}
307  }
308 
309 #ifdef debug_QGSMfragmentation
310  G4cout<<"d1 d2 "<<d1<<" "<<d2<<G4endl;
311 #endif
312 
313  d1+=1.0; d2+=1.0;
314  invD1=1./d1; invD2=1./d2;
315 
316  const G4int maxNumberOfLoops = 10000;
317  G4int loopCounter = 0;
318  do
319  {
320  r1=G4Pow::GetInstance()->powA(G4UniformRand(),invD1);
321  r2=G4Pow::GetInstance()->powA(G4UniformRand(),invD2);
322  r12=r1+r2;
323  z=r1/r12;
324  } while( ( (r12 > 1.0) || !((zmin <= z)&&(z <= zmax))) && ++loopCounter < maxNumberOfLoops ); /* Loop checking, 07.08.2015, A.Ribon */
325  if ( loopCounter >= maxNumberOfLoops ) {
326  z = 0.5*(zmin + zmax); // Just a value between zmin and zmax, no physics considerations at all!
327  }
328 
329  return z;
330  }
331  else if(absCode == 3101 || absCode == 3103 || // For strange d-quarks
332  absCode == 3201 || absCode == 3203)
333  {
334  d2 = (alft - (2.*ala - arho));
335 
336  }
337  else
338  {
339  d2 = (alft - (2.*aksi - arho));
340  }
341 
342  const G4int maxNumberOfLoops = 1000;
343  G4int loopCounter = 0;
344  do
345  {
346  z = zmin + G4UniformRand() * (zmax - zmin);
347  d1 = (1. - z);
348  yf = G4Pow::GetInstance()->powA(d1, d2);
349  }
350  while( (G4UniformRand() > yf) && ++loopCounter < maxNumberOfLoops ); /* Loop checking, 07.08.2015, A.Ribon */
351  }
352 
353  return z;
354 }
355 //-----------------------------------------------------------------------------------------
356 
358  G4FragmentingString * string,
360 {
361  G4double HadronMass = pHadron->GetPDGMass();
362 
363 // G4double MinimalStringMass= FragmentationMass(NewString,&G4HadronBuilder::BuildHighSpin);
364  G4double MinimalStringMass=
366 // FragmentationMass(NewString,&G4HadronBuilder::BuildLowSpin); // Uzhi 03.06.2015
367 
368 #ifdef debug_QGSMfragmentation
369  G4cout<<"G4QGSMFragmentation::SplitEandP "<<pHadron->GetParticleName()<<G4endl;
370  G4cout<<"String 4 mom, String M "<<string->Get4Momentum()<<" "<<string->Mass()<<G4endl;
371  G4cout<<"HadM MinimalStringMassLeft StringM hM+sM "<<HadronMass<<" "<<MinimalStringMass<<" "
372  <<string->Mass()<<" "<<HadronMass+MinimalStringMass<<G4endl;
373 #endif
374 
375  if(HadronMass + MinimalStringMass > string->Mass())
376  {
377 #ifdef debug_QGSMfragmentation
378  G4cout<<"Mass of the string is not sufficient to produce the hadron!"<<G4endl;
379 #endif
380  return 0;
381  }// have to start all over!
382 
383  // calculate and assign hadron transverse momentum component HadronPx andHadronPy
384  G4double StringMT2 = string->MassT2();
385  G4double StringMT = std::sqrt(StringMT2);
386 
387  G4LorentzVector String4Momentum = string->Get4Momentum();
388  String4Momentum.setPz(0.);
389  G4ThreeVector StringPt = String4Momentum.vect();
390 
391  G4ThreeVector HadronPt , RemSysPt;
392  G4double HadronMassT2, ResidualMassT2;
393 
394  //... sample Pt of the hadron
395  G4int attempt=0;
396  do
397  {
398  attempt++; if(attempt > StringLoopInterrupt) return 0;
399 
400  HadronPt =SampleQuarkPt() + string->DecayPt();
401  HadronPt.setZ(0);
402  RemSysPt = StringPt - HadronPt;
403 
404  HadronMassT2 = sqr(HadronMass) + HadronPt.mag2();
405  ResidualMassT2=sqr(MinimalStringMass) + RemSysPt.mag2();
406 
407  } while(std::sqrt(HadronMassT2) + std::sqrt(ResidualMassT2) > StringMT); /* Loop checking, 07.08.2015, A.Ribon */
408 
409  //... sample z to define hadron longitudinal momentum and energy
410  //... but first check the available phase space
411 
412  G4double Pz2 = (sqr(StringMT2 - HadronMassT2 - ResidualMassT2) -
413  4*HadronMassT2 * ResidualMassT2)/4./StringMT2;
414 
415  if(Pz2 < 0 ) {return 0;} // have to start all over!
416 
417  //... then compute allowed z region z_min <= z <= z_max
418 
419  G4double Pz = std::sqrt(Pz2);
420  G4double zMin = (std::sqrt(HadronMassT2+Pz2) - Pz)/std::sqrt(StringMT2);
421  G4double zMax = (std::sqrt(HadronMassT2+Pz2) + Pz)/std::sqrt(StringMT2);
422 
423  if (zMin >= zMax) return 0; // have to start all over!
424 
425  G4double z = GetLightConeZ(zMin, zMax,
426  string->GetDecayParton()->GetPDGEncoding(), pHadron,
427  HadronPt.x(), HadronPt.y());
428 
429  //... now compute hadron longitudinal momentum and energy
430  // longitudinal hadron momentum component HadronPz
431 
432  HadronPt.setZ(0.5* string->GetDecayDirection() *
433  (z * string->LightConeDecay() -
434  HadronMassT2/(z * string->LightConeDecay())));
435  G4double HadronE = 0.5* (z * string->LightConeDecay() +
436  HadronMassT2/(z * string->LightConeDecay()));
437 
438  G4LorentzVector * a4Momentum= new G4LorentzVector(HadronPt,HadronE);
439 
440 #ifdef debug_QGSMfragmentation
441  G4cout<<"string->GetDecayDirection() string->LightConeDecay() "
442  <<string->GetDecayDirection()<<" "<<string->LightConeDecay()<<G4endl;
443  G4cout<<"HadronPt,HadronE "<<HadronPt<<" "<<HadronE<<G4endl;
444  G4cout<<"Out of QGSM SplitEandP "<<G4endl;
445 #endif
446 
447  return a4Momentum;
448 }
449 
450 
451 //-----------------------------------------------------------------------------------------
452 
454  G4KineticTrackVector * LeftVector,
455  G4KineticTrackVector * RightVector)
456 {
457  //... perform last cluster decay
458 
459  G4ThreeVector ClusterVel =string->Get4Momentum().boostVector();
460  G4double ResidualMass =string->Mass();
461 
462 #ifdef debug_QGSMfragmentation
463  G4cout<<"Split last-----------------------------------------"<<G4endl;
464  G4cout<<"StrMass "<<ResidualMass<<" q's "
465  <<string->GetLeftParton()->GetParticleName()<<" "
466  <<string->GetRightParton()->GetParticleName()<<G4endl;
467 #endif
468 
469  G4double ClusterMassCut = ClusterMass; // Taken from G4VLongitudinalStringDecay
470  G4int cClusterInterrupt = 0;
471  G4ParticleDefinition * LeftHadron, * RightHadron;
472  const G4int maxNumberOfLoops = 1000;
473  G4int loopCounter = 0;
474  do
475  {
476  if (cClusterInterrupt++ >= ClusterLoopInterrupt)
477  {
478  return false;
479  }
480 
481  G4ParticleDefinition * quark = NULL;
482  string->SetLeftPartonStable(); // to query quark contents..
483 
484  if (string->DecayIsQuark() && string->StableIsQuark() )
485  {
486  //... there are quarks on cluster ends
487  LeftHadron= QuarkSplitup(string->GetLeftParton(), quark);
488  } else {
489  //... there is a Diquark on cluster ends
490  G4int IsParticle;
491  if ( string->StableIsQuark() ) {
492  IsParticle=(string->GetLeftParton()->GetPDGEncoding()>0) ? -1 : +1;
493  } else {
494  IsParticle=(string->GetLeftParton()->GetPDGEncoding()>0) ? +1 : -1;
495  }
496 
497 //G4double ProbSaS = 1.0 - 2.0 * GetStrangeSuppress();
498 //G4double ActualProb = ProbSaS * 1.4;
499 //SetStrangenessSuppression((1.0-ActualProb)/2.0);
500 
501  pDefPair QuarkPair = CreatePartonPair(IsParticle,false); // no diquarks wanted
502 //SetStrangenessSuppression((1.0-ProbSaS)/2.0);
503  quark = QuarkPair.second;
504  LeftHadron=hadronizer->Build(QuarkPair.first, string->GetLeftParton());
505  }
506 
507  RightHadron = hadronizer->Build(string->GetRightParton(), quark);
508 
509  //... repeat procedure, if mass of cluster is too low to produce hadrons
510  //... ClusterMassCut = 0.15*GeV model parameter
511  if ( quark->GetParticleSubType()== "quark" ) {ClusterMassCut = 0.;}
512  else {ClusterMassCut = ClusterMass;}
513  }
514  while ( (ResidualMass <= LeftHadron->GetPDGMass() + RightHadron->GetPDGMass() + ClusterMassCut)
515  && ++loopCounter < maxNumberOfLoops ); /* Loop checking, 07.08.2015, A.Ribon */
516 
517  if ( loopCounter >= maxNumberOfLoops ) {
518  return false;
519  }
520 
521  //... compute hadron momenta and energies
522  G4LorentzVector LeftMom, RightMom;
524  Sample4Momentum(&LeftMom , LeftHadron->GetPDGMass() ,
525  &RightMom, RightHadron->GetPDGMass(), ResidualMass);
526  LeftMom.boost(ClusterVel);
527  RightMom.boost(ClusterVel);
528 
529 #ifdef debug_QGSMfragmentation
530  G4cout<<LeftHadron->GetParticleName()<<" "<<RightHadron->GetParticleName()<<G4endl;
531  G4cout<<"Left Hadrom P M "<<LeftMom<<" "<<LeftMom.mag()<<G4endl;
532  G4cout<<"Right Hadrom P M "<<RightMom<<" "<<RightMom.mag()<<G4endl;
533 #endif
534 
535  LeftVector->push_back(new G4KineticTrack(LeftHadron, 0, Pos, LeftMom));
536  RightVector->push_back(new G4KineticTrack(RightHadron, 0, Pos, RightMom));
537 
538  return true;
539 
540 }
541 
542 //----------------------------------------------------------------------------------------------------------
543 
545 {
546  return sqr(FragmentationMass(string)+MassCut) <
547  string->Mass2();
548 }
549 
550 //----------------------------------------------------------------------------------------------------------
551 
553 {
554  return
556  string->Get4Momentum().mag2();
557 }
558 
559 //----------------------------------------------------------------------------------------------------------
560 
561 //----------------------------------------------------------------------------------------------------------
562 
564  G4LorentzVector* AntiMom, G4double AntiMass, G4double InitialMass)
565  {
566  G4double r_val = sqr(InitialMass*InitialMass - Mass*Mass - AntiMass*AntiMass) - sqr(2.*Mass*AntiMass);
567  G4double Pabs = (r_val > 0.)? std::sqrt(r_val)/(2.*InitialMass) : 0;
568 
569  //... sample unit vector
570  G4double pz = 1. - 2.*G4UniformRand();
571  G4double st = std::sqrt(1. - pz * pz)*Pabs;
572  G4double phi = 2.*pi*G4UniformRand();
573  G4double px = st*std::cos(phi);
574  G4double py = st*std::sin(phi);
575  pz *= Pabs;
576 
577  Mom->setPx(px); Mom->setPy(py); Mom->setPz(pz);
578  Mom->setE(std::sqrt(Pabs*Pabs + Mass*Mass));
579 
580  AntiMom->setPx(-px); AntiMom->setPy(-py); AntiMom->setPz(-pz);
581  AntiMom->setE (std::sqrt(Pabs*Pabs + AntiMass*AntiMass));
582  }
583 
584 //-----------------------------------------------------------------------------
585 
588  G4ParticleDefinition *&created)
589 {
590  //... can Diquark break or not?
592  //... Diquark break
593 
594  G4int stableQuarkEncoding = decay->GetPDGEncoding()/1000;
595  G4int decayQuarkEncoding = (decay->GetPDGEncoding()/100)%10;
596 
597  if (G4UniformRand() < 0.5)
598  {
599  G4int Swap = stableQuarkEncoding;
600  stableQuarkEncoding = decayQuarkEncoding;
601  decayQuarkEncoding = Swap;
602  }
603 
604  G4int IsParticle=(decayQuarkEncoding>0) ? -1 : +1;
605  // if we have a quark, we need antiquark)
606 
607  G4double StrSup=GetStrangeSuppress();
608  StrangeSuppress=0.41; // was 0.47 Uzhi Dec. 2017
609  pDefPair QuarkPair = CreatePartonPair(IsParticle,false); // no diquarks wanted
610  StrangeSuppress=StrSup;
611 
612  //... Build new Diquark
613  G4int QuarkEncoding=QuarkPair.second->GetPDGEncoding();
614  G4int i10 = std::max(std::abs(QuarkEncoding), std::abs(stableQuarkEncoding));
615  G4int i20 = std::min(std::abs(QuarkEncoding), std::abs(stableQuarkEncoding));
616  G4int spin = (i10 != i20 && G4UniformRand() <= 0.5)? 1 : 3;
617  G4int NewDecayEncoding = -1*IsParticle*(i10 * 1000 + i20 * 100 + spin);
618  created = FindParticle(NewDecayEncoding);
619  G4ParticleDefinition * decayQuark=FindParticle(decayQuarkEncoding);
620  G4ParticleDefinition * had=hadronizer->Build(QuarkPair.first, decayQuark);
621 
622  return had;
623 // return hadronizer->Build(QuarkPair.first, decayQuark);
624 
625  } else {
626  //... Diquark does not break
627 
628  G4int IsParticle=(decay->GetPDGEncoding()>0) ? +1 : -1;
629  // if we have a diquark, we need quark)
630  G4double StrSup=GetStrangeSuppress(); // for changing s-sbar production, Uzhi Oct. 2014
631  StrangeSuppress=0.41; //0.41; 0.47
632  pDefPair QuarkPair = CreatePartonPair(IsParticle,false); // no diquarks wanted
633  StrangeSuppress=StrSup;
634 
635  created = QuarkPair.second;
636 
637  G4ParticleDefinition * had=hadronizer->Build(QuarkPair.first, decay);
638  return had;
639 // return G4ParticleDefinition * had=hadronizer->Build(QuarkPair.first, decay);
640  }
641 }
642 
643 //-----------------------------------------------------------------------------
644 
646  G4FragmentingString *string,
647  G4FragmentingString *&newString)
648 {
649 #ifdef debug_QGSMfragmentation
650  G4cout<<G4endl;
651  G4cout<<"Start SplitUP (G4VLongitudinalStringDecay) ========================="<<G4endl;
652  G4cout<<"String partons: " <<string->GetLeftParton()->GetPDGEncoding()<<" "
653  <<string->GetRightParton()->GetPDGEncoding()<<" "
654  <<"Direction " <<string->GetDecayDirection()<<G4endl;
655 #endif
656 
657  //... random choice of string end to use for creating the hadron (decay)
658  G4int SideOfDecay = (G4UniformRand() < 0.5)? 1: -1;
659  if (SideOfDecay < 0)
660  {
661  string->SetLeftPartonStable();
662  } else
663  {
664  string->SetRightPartonStable();
665  }
666 
667  G4ParticleDefinition *newStringEnd;
668  G4ParticleDefinition * HadronDefinition;
669  if (string->DecayIsQuark())
670  {
671 // Uzhi Dec. 2017, Start
672 G4double ProbDqADq = GetDiquarkSuppress();
673 
674 G4int NumberOfpossibleBaryons = 2;
675 
676 if(string->GetLeftParton()->GetParticleSubType() != "quark") NumberOfpossibleBaryons++;
677 if(string->GetRightParton()->GetParticleSubType() != "quark") NumberOfpossibleBaryons++;
678 
679 G4double ActualProb = ProbDqADq ;
680  ActualProb *= (1.0-G4Exp(2.0*(1.0 - string->Mass()/(NumberOfpossibleBaryons*1400.0))));
681 
682 SetDiquarkSuppression(ActualProb);
683 // Uzhi Dec. 2017, End
684  HadronDefinition= QuarkSplitup(string->GetDecayParton(), newStringEnd);
685 
686 SetDiquarkSuppression(ProbDqADq); // Uzhi Dec. 2017
687  } else {
688  HadronDefinition= DiQuarkSplitup(string->GetDecayParton(), newStringEnd);
689  }
690 
691 #ifdef debug_QGSMfragmentation
692  G4cout<<"The parton "<<string->GetDecayParton()->GetPDGEncoding()<<" "
693  <<" produces hadron "<<HadronDefinition->GetParticleName()
694  <<" and is transformed to "<<newStringEnd->GetPDGEncoding()<<G4endl;
695  G4cout<<"The side of the string decay Left/Right (1/-1) "<<SideOfDecay<<G4endl;
696 #endif
697 // create new String from old, ie. keep Left and Right order, but replace decay
698 
699  newString=new G4FragmentingString(*string,newStringEnd); // To store possible
700  // quark containt of new string
701 
702 #ifdef debug_QGSMfragmentation
703  G4cout<<"An attempt to determine its energy (SplitEandP)"<<G4endl;
704 #endif
705  G4LorentzVector* HadronMomentum=SplitEandP(HadronDefinition, string, newString);
706 
707  delete newString; newString=0;
708 
709  G4KineticTrack * Hadron =0;
710  if ( HadronMomentum != 0 ) {
711 
712 #ifdef debug_QGSMfragmentation
713  G4cout<<"The attempt was successful"<<G4endl;
714 #endif
716  Hadron = new G4KineticTrack(HadronDefinition, 0,Pos, *HadronMomentum);
717 
718  newString=new G4FragmentingString(*string,newStringEnd,
719  HadronMomentum);
720 
721  delete HadronMomentum;
722  }
723  else
724  {
725 
726 #ifdef debug_QGSMfragmentation
727  G4cout<<"The attempt was not successful !!!"<<G4endl;
728 #endif
729  }
730 
731 #ifdef debug_VStringDecay
732  G4cout<<"End SplitUP (G4VLongitudinalStringDecay) ====================="<<G4endl;
733 #endif
734 
735  return Hadron;
736 }
737 
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
ParticleList decay(Cluster *const c)
Carries out a cluster decay.
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4double Mass() const
G4Parton * GetRightParton(void) const
G4ParticleDefinition * Build(G4ParticleDefinition *black, G4ParticleDefinition *white)
virtual G4KineticTrackVector * FragmentString(const G4ExcitedString &theString)
const G4ThreeVector & GetPosition() const
const G4LorentzVector & Get4Momentum() const
void SetDiquarkSuppression(G4double aValue)
G4ParticleDefinition * GetRightParton(void) const
const double C1
const G4String & GetParticleSubType() const
#define G4endl
Definition: G4ios.hh:61
ush Pos
Definition: deflate.h:92
virtual G4LorentzVector * SplitEandP(G4ParticleDefinition *pHadron, G4FragmentingString *string, G4FragmentingString *newString)
Double_t z
const G4String & GetParticleName() const
const G4ParticleDefinition * GetDefinition() const
void SetVectorMesonMixings(std::vector< G4double > aVector)
G4double FragmentationMass(const G4FragmentingString *const string, Pcreate build=0, pDefPair *pdefs=0)
G4double GetPDGMass() const
void SetSpinThreeHalfBarionProbability(G4double aValue)
void SetDiquarkBreakProbability(G4double aValue)
std::vector< G4double > vectorMesonMix
#define NewString(str)
G4int GetDirection(void) const
double G4double
Definition: G4Types.hh:76
bool G4bool
Definition: G4Types.hh:79
void SetVectorMesonProbability(G4double aValue)
static G4Pow * GetInstance()
Definition: G4Pow.cc:57
G4ParticleDefinition * FindParticle(G4int Encoding)
G4int GetDecayDirection() const
G4ThreeVector SampleQuarkPt(G4double ptMax=-1.)
static const G4double d2
void CalculateHadronTimePosition(G4double theInitialStringMass, G4KineticTrackVector *)
G4double powA(G4double A, G4double y) const
Definition: G4Pow.hh:242
void setZ(double)
G4KineticTrackVector * LightFragmentationTest(const G4ExcitedString *const theString)
static const G4double d1
#define G4UniformRand()
Definition: Randomize.hh:53
void SetStrangenessSuppression(G4double aValue)
G4LorentzVector Get4Momentum() const
HepLorentzRotation inverse() const
G4Parton * GetLeftParton(void) const
double mag2() const
virtual G4bool SplitLast(G4FragmentingString *string, G4KineticTrackVector *LeftVector, G4KineticTrackVector *RightVector)
G4ParticleDefinition * QuarkSplitup(G4ParticleDefinition *decay, G4ParticleDefinition *&created)
virtual void Sample4Momentum(G4LorentzVector *Mom, G4double Mass, G4LorentzVector *AntiMom, G4double AntiMass, G4double InitialMass)
G4ParticleDefinition * BuildHighSpin(G4ParticleDefinition *black, G4ParticleDefinition *white)
G4int GetPDGcode() const
Definition: G4Parton.hh:127
double mag() const
G4double GetFormationTime() const
int G4int
Definition: G4Types.hh:78
G4ParticleDefinition * GetLeftParton(void) const
void SetPosition(const G4ThreeVector aPosition)
virtual G4bool IsFragmentable(const G4FragmentingString *const string)
pDefPair CreatePartonPair(G4int NeedParticle, G4bool AllowDiquarks=true)
double pz() const
virtual G4KineticTrack * Splitup(G4FragmentingString *string, G4FragmentingString *&newString)
G4ParticleDefinition * GetDecayParton() const
G4GLOB_DLL std::ostream G4cout
double x() const
virtual G4ParticleDefinition * DiQuarkSplitup(G4ParticleDefinition *decay, G4ParticleDefinition *&created)
virtual G4double GetLightConeZ(G4double zmin, G4double zmax, G4int PartonEncoding, G4ParticleDefinition *pHadron, G4double Px, G4double Py)
Hep3Vector vect() const
T sqr(const T &x)
Definition: templates.hh:145
G4ExcitedString * CopyExcited(const G4ExcitedString &string)
static constexpr double pi
Definition: G4SIunits.hh:75
const G4ThreeVector & GetPosition() const
CLHEP::HepLorentzVector G4LorentzVector
std::pair< G4ParticleDefinition *, G4ParticleDefinition * > pDefPair
double y() const
const G4LorentzVector & Get4Momentum() const
Definition: G4Parton.hh:143
void Set4Momentum(const G4LorentzVector &a4Momentum)
void SetFormationTime(G4double aFormationTime)
G4LorentzRotation TransformToAlignedCms()
virtual G4bool StopFragmenting(const G4FragmentingString *const string)
HepLorentzVector & boost(double, double, double)
T min(const T t1, const T t2)
brief Return the smallest of the two arguments