Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4QGSMFragmentation.cc
Go to the documentation of this file.
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 100828 2016-11-02 15:25:59Z 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 {
55 }
56 
58 {
59 }
60 
61 //----------------------------------------------------------------------------------------------------------
62 
64 {
65  #ifdef debug_QGSMfragmentation
66  G4cout<<G4endl<<"QGSM StringFragm: String Mass " <<theString.Get4Momentum().mag()<<" Pz "
67  <<theString.Get4Momentum().pz()
68  <<"------------------------------------"<<G4endl;
69  G4cout<<"String ends Direct "<<theString.GetLeftParton()->GetPDGcode()<<" "
70  <<theString.GetRightParton()->GetPDGcode()<<" "
71  <<theString.GetDirection()<< G4endl;
72  G4cout<<"Left mom "<<theString.GetLeftParton()->Get4Momentum()<<G4endl;
73  G4cout<<"Right mom "<<theString.GetRightParton()->Get4Momentum()<<G4endl;
74  G4cout<<"Check for Fragmentation "<<G4endl;
75  #endif
76 
77  // Can no longer modify Parameters for Fragmentation.
78  PastInitPhase=true;
79 
80  // check if string has enough mass to fragment...
81 
82  G4KineticTrackVector * LeftVector=LightFragmentationTest(&theString);
83 
84  #ifdef debug_QGSMfragmentation
85  if ( LeftVector != 0 ) G4cout<<"Non fragmentable - the string is converted to one hadron "<<G4endl;
86  #endif
87 
88  if ( LeftVector != 0 ) return LeftVector;
89 
90  #ifdef debug_QGSMfragmentation
91  G4cout<<"The string will be fragmented. "<<G4endl;
92  #endif
93 
94  LeftVector = new G4KineticTrackVector;
96 
97  // this should work but its only a semi deep copy. %GF G4ExcitedString theStringInCMS(theString);
98  G4ExcitedString *theStringInCMS=CPExcited(theString);
99  G4LorentzRotation toCms=theStringInCMS->TransformToAlignedCms();
100 
101  G4bool success=false, inner_sucess=true;
102  G4int attempt=0;
103  while ( !success && attempt++ < StringLoopInterrupt ) /* Loop checking, 07.08.2015, A.Ribon */
104  {
105  #ifdef debug_QGSMfragmentation
106  G4cout<<"Loop_toFrag "<<theStringInCMS->GetLeftParton()->GetPDGcode()<<" "
107  <<theStringInCMS->GetRightParton()->GetPDGcode()<<" "
108  <<theStringInCMS->GetDirection()<< G4endl;
109  #endif
110 
111  G4FragmentingString *currentString=new G4FragmentingString(*theStringInCMS);
112 
113  std::for_each(LeftVector->begin(), LeftVector->end(), DeleteKineticTrack());
114  LeftVector->clear();
115  std::for_each(RightVector->begin(), RightVector->end(), DeleteKineticTrack());
116  RightVector->clear();
117 
118  inner_sucess=true; // set false on failure..
119  const G4int maxNumberOfLoops = 1000;
120  G4int loopCounter = -1;
121  while (! StopFragmenting(currentString) && ++loopCounter < maxNumberOfLoops ) /* Loop checking, 07.08.2015, A.Ribon */
122  { // Split current string into hadron + new string
123 
124  #ifdef debug_QGSMfragmentation
125  G4cout<<"The string can fragment. "<<G4endl;;
126  #endif
127 
128  G4FragmentingString *newString=0; // used as output from SplitUp...
129  G4KineticTrack * Hadron=Splitup(currentString,newString);
130 
131  if ( Hadron != 0 ) // && IsFragmentable(newString))
132  {
133  #ifdef debug_QGSMfragmentation
134  G4cout<<"Hadron prod at fragm. "<<Hadron->GetDefinition()->GetParticleName()<<G4endl;
135  #endif
136 
137  if ( currentString->GetDecayDirection() > 0 ) LeftVector->push_back(Hadron);
138  else RightVector->push_back(Hadron);
139 
140  delete currentString;
141  currentString=newString;
142 
143  } else {
144 
145  #ifdef debug_QGSMfragmentation
146  G4cout<<"abandon ... start from the beginning ---------------"<<G4endl;
147  #endif
148 
149  // abandon ... start from the beginning
150  if (newString) delete newString;
151  inner_sucess=false;
152  break;
153  }
154  }
155  if ( loopCounter >= maxNumberOfLoops ) inner_sucess=false;
156 
157  // Split current string into 2 final Hadrons
158  #ifdef debug_QGSMfragmentation
159  G4cout<<"Split remaining string into 2 final hadrons."<<G4endl;
160  #endif
161 
162  if ( inner_sucess && SplitLast(currentString,LeftVector, RightVector) )
163  {
164  success=true;
165  }
166  delete currentString;
167  } // End of while loop
168 
169  delete theStringInCMS;
170 
171  if ( ! success )
172  {
173  std::for_each(LeftVector->begin(), LeftVector->end(), DeleteKineticTrack());
174  LeftVector->clear();
175  std::for_each(RightVector->begin(), RightVector->end(), DeleteKineticTrack());
176  delete RightVector;
177  return LeftVector;
178  }
179 
180  // Join Left- and RightVector into LeftVector in correct order.
181  while(!RightVector->empty()) /* Loop checking, 07.08.2015, A.Ribon */
182  {
183  LeftVector->push_back(RightVector->back());
184  RightVector->erase(RightVector->end()-1);
185  }
186  delete RightVector;
187 
188  CalculateHadronTimePosition(theString.Get4Momentum().mag(), LeftVector);
189 
190  G4LorentzRotation toObserverFrame(toCms.inverse());
191 
192  for(size_t C1 = 0; C1 < LeftVector->size(); C1++)
193  {
194  G4KineticTrack* Hadron = LeftVector->operator[](C1);
195  G4LorentzVector Momentum = Hadron->Get4Momentum();
196  Momentum = toObserverFrame*Momentum;
197  Hadron->Set4Momentum(Momentum);
198  G4LorentzVector Coordinate(Hadron->GetPosition(), Hadron->GetFormationTime());
199  Momentum = toObserverFrame*Coordinate;
200  Hadron->SetFormationTime(Momentum.e());
201  G4ThreeVector aPosition(Momentum.vect());
202  Hadron->SetPosition(theString.GetPosition()+aPosition);
203  }
204  return LeftVector;
205 }
206 
207 //----------------------------------------------------------------------------------------------------------
208 
209 G4double G4QGSMFragmentation::GetLightConeZ(G4double zmin, G4double zmax, G4int PartonEncoding,
211 {
212  #ifdef debug_QGSMfragmentation
213  G4cout<<"GetLightConeZ zmin zmax Parton pHadron "<<zmin<<" "<<zmax<<" "<<PartonEncoding<<" "<<pHadron->GetParticleName()<<G4endl;
214  #endif
215 
216  G4double z;
217  G4double d1, d2, yf;
218  G4double invD1(0.),invD2(0.), r1(0.),r2(0.),r12(0.);
219 
220  G4int absCode = std::abs( PartonEncoding );
221  G4int absHadronCode=std::abs(pHadron->GetPDGEncoding());
222 
223  G4int q1, q2, q3;
224  q1 = absHadronCode/1000; q2 = (absHadronCode % 1000)/100; q3 = (absHadronCode % 100)/10;
225 
226  G4bool StrangeHadron = (q1 == 3) || (q2 == 3) || (q3 == 3);
227 
228  if (absCode < 10)
229  { // A quark fragmentation ----------------------------
230  if(absCode == 1 || absCode == 2)
231  {
232  if(absHadronCode < 1000)
233  { // Meson produced
234  if( !StrangeHadron ) {d1=2.0; d2 = -arho + alft;}
235  else {d1=1.0; d2 = -aphi + alft;}
236  } else
237  { // Baryon produced
238  if( !StrangeHadron ) {d1=0.0; d2 = arho - 2.0*an + alft;}
239  else {d1=0.0; d2 = 2.0*arho - 2.0*an - aphi + alft;}
240  }
241  }
242  else if(absCode == 3)
243  {
244  if(absHadronCode < 1000){d1=1.0 - aphi; d2 = -arho + alft;} // Meson produced s->K + u/d
245  else {d1=1.0 - aphi; d2 = arho - 2.0*an + alft;} // Baryon produced
246 
247  } else throw G4HadronicException(__FILE__, __LINE__, "Unknown PDGencoding in G4QGSMFragmentation::G4LightConeZ");
248 
249  #ifdef debug_QGSMfragmentation
250  G4cout<<"d1 d2 "<<d1<<" "<<d2<<G4endl;
251  #endif
252 
253  d1+=1.0; d2+=1.0;
254 
255  invD1=1./d1; invD2=1./d2;
256 
257  const G4int maxNumberOfLoops = 10000;
258  G4int loopCounter = 0;
259  do
260  {
261  r1=G4Pow::GetInstance()->powA(G4UniformRand(),invD1);
262  r2=G4Pow::GetInstance()->powA(G4UniformRand(),invD2);
263  r12=r1+r2;
264  z=r1/r12;
265  } while( ( (r12 > 1.0) || !((zmin <= z)&&(z <= zmax))) && ++loopCounter < maxNumberOfLoops ); /* Loop checking, 07.08.2015, A.Ribon */
266  if ( loopCounter >= maxNumberOfLoops ) {
267  z = 0.5*(zmin + zmax); // Just a value between zmin and zmax, no physics considerations at all!
268  }
269 
270  return z;
271  }
272  else
273  { // A di-quark fragmentation -------------------------
274  if(absCode == 1103 || absCode == 2101 ||
275  absCode == 2203 || absCode == 2103)
276  {
277  if(absHadronCode < 1000) // Meson production
278  {
279  if( !StrangeHadron ) {d1=1.0; d2= arho - 2.0*an + alft;}
280  else {d1=1.0; d2 = 2.*arho - 2.0*an - aphi + alft;}
281  } else { // Baryon production
282  if( !StrangeHadron ) {d1=2.0*(arho - an); d2= -arho + alft;}
283  else {d1=2.0*(arho - an); d2 =-aphi + alft;}
284  }
285 
286  #ifdef debug_QGSMfragmentation
287  G4cout<<"d1 d2 "<<d1<<" "<<d2<<G4endl;
288  #endif
289 
290  d1+=1.0; d2+=1.0;
291  invD1=1./d1; invD2=1./d2;
292 
293  const G4int maxNumberOfLoops = 10000;
294  G4int loopCounter = 0;
295  do
296  {
297  r1=G4Pow::GetInstance()->powA(G4UniformRand(),invD1);
298  r2=G4Pow::GetInstance()->powA(G4UniformRand(),invD2);
299  r12=r1+r2;
300  z=r1/r12;
301  } while( ( (r12 > 1.0) || !((zmin <= z)&&(z <= zmax))) && ++loopCounter < maxNumberOfLoops ); /* Loop checking, 07.08.2015, A.Ribon */
302  if ( loopCounter >= maxNumberOfLoops ) {
303  z = 0.5*(zmin + zmax); // Just a value between zmin and zmax, no physics considerations at all!
304  }
305 
306  return z;
307  }
308  else if(absCode == 3101 || absCode == 3103 || // For strange d-quarks
309  absCode == 3201 || absCode == 3203)
310  {
311  d2 = (alft - (2.*ala - arho));
312  }
313  else
314  {
315  d2 = (alft - (2.*aksi - arho));
316  }
317 
318  const G4int maxNumberOfLoops = 1000;
319  G4int loopCounter = 0;
320  do
321  {
322  z = zmin + G4UniformRand() * (zmax - zmin);
323  d1 = (1. - z);
324  yf = G4Pow::GetInstance()->powA(d1, d2);
325  }
326  while( (G4UniformRand() > yf) && ++loopCounter < maxNumberOfLoops ); /* Loop checking, 07.08.2015, A.Ribon */
327  }
328 
329  return z;
330 }
331 
332 //-----------------------------------------------------------------------------------------
333 
334 G4LorentzVector * G4QGSMFragmentation::SplitEandP(G4ParticleDefinition * pHadron,
335  G4FragmentingString * string,
337 {
338  G4double HadronMass = pHadron->GetPDGMass();
339 
340  //G4double MinimalStringMass= FragmentationMass(NewString,&G4HadronBuilder::BuildHighSpin);
341  G4double MinimalStringMass= FragmentationMass(NewString,&G4HadronBuilder::Build);
342  //FragmentationMass(NewString,&G4HadronBuilder::BuildLowSpin); // Uzhi 03.06.2015
343  // Uzhi 03.06.2015 It would be well to sample randomly HighSpin
344 
345  #ifdef debug_QGSMfragmentation
346  G4cout<<"G4QGSMFragmentation::SplitEandP "<<pHadron->GetParticleName()<<G4endl;
347  G4cout<<"String 4 mom, String M "<<string->Get4Momentum()<<" "<<string->Mass()<<G4endl;
348  G4cout<<"HadM MinimalStringMassLeft StringM hM+sM "<<HadronMass<<" "<<MinimalStringMass<<" "
349  <<string->Mass()<<" "<<HadronMass+MinimalStringMass<<G4endl;
350  #endif
351 
352  if(HadronMass + MinimalStringMass > string->Mass())
353  {
354  #ifdef debug_QGSMfragmentation
355  G4cout<<"Mass of the string is not sufficient to produce the hadron!"<<G4endl;
356  #endif
357  return 0;
358  } // have to start all over!
359 
360  // calculate and assign hadron transverse momentum component HadronPx andHadronPy
361  G4double StringMT2 = string->MassT2();
362  G4double StringMT = std::sqrt(StringMT2);
363 
364  G4LorentzVector String4Momentum = string->Get4Momentum();
365  String4Momentum.setPz(0.);
366  G4ThreeVector StringPt = String4Momentum.vect();
367 
368  G4ThreeVector HadronPt , RemSysPt;
369  G4double HadronMassT2, ResidualMassT2;
370 
371  //... sample Pt of the hadron
372  G4int attempt=0;
373  do
374  {
375  attempt++; if(attempt > StringLoopInterrupt) return 0;
376 
377  HadronPt =SampleQuarkPt() + string->DecayPt();
378  HadronPt.setZ(0);
379  RemSysPt = StringPt - HadronPt;
380 
381  HadronMassT2 = sqr(HadronMass) + HadronPt.mag2();
382  ResidualMassT2=sqr(MinimalStringMass) + RemSysPt.mag2();
383 
384  } while(std::sqrt(HadronMassT2) + std::sqrt(ResidualMassT2) > StringMT); /* Loop checking, 07.08.2015, A.Ribon */
385 
386  //... sample z to define hadron longitudinal momentum and energy
387  //... but first check the available phase space
388 
389  G4double Pz2 = (sqr(StringMT2 - HadronMassT2 - ResidualMassT2) -
390  4*HadronMassT2 * ResidualMassT2)/4./StringMT2;
391 
392  if(Pz2 < 0 ) {return 0;} // have to start all over!
393 
394  //... then compute allowed z region z_min <= z <= z_max
395 
396  G4double Pz = std::sqrt(Pz2);
397  G4double zMin = (std::sqrt(HadronMassT2+Pz2) - Pz)/std::sqrt(StringMT2);
398  G4double zMax = (std::sqrt(HadronMassT2+Pz2) + Pz)/std::sqrt(StringMT2);
399 
400  /*
401  G4double DecayQuarkMass2 = sqr(string->GetDecayParton()->GetPDGMass());
402  G4double HadronMass2T = sqr(HadronMass) + HadronPt.mag2();
403  if (DecayQuarkMass2 + HadronMass2T >= SmoothParam*(string->Mass2()) )
404  return 0; // have to start all over!
405  //... then compute allowed z region z_min <= z <= z_max
406  //G4double zMin = HadronMass2T/(string->Mass2());
407  //G4double zMax = 1. - DecayQuarkMass2/(string->Mass2());
408  */
409 
410  if (zMin >= zMax) return 0; // have to start all over!
411 
412  G4double z = GetLightConeZ(zMin, zMax, string->GetDecayParton()->GetPDGEncoding(),
413  pHadron, HadronPt.x(), HadronPt.y());
414 
415  //... now compute hadron longitudinal momentum and energy
416  // longitudinal hadron momentum component HadronPz
417 
418  HadronPt.setZ(0.5* string->GetDecayDirection() * (z * string->LightConeDecay() -
419  HadronMassT2/(z * string->LightConeDecay())));
420  G4double HadronE = 0.5* (z * string->LightConeDecay() +
421  HadronMassT2/(z * string->LightConeDecay()));
422 
423  G4LorentzVector * a4Momentum= new G4LorentzVector(HadronPt,HadronE);
424 
425  #ifdef debug_QGSMfragmentation
426  G4cout<<"string->GetDecayDirection() string->LightConeDecay() "
427  <<string->GetDecayDirection()<<" "<<string->LightConeDecay()<<G4endl;
428  G4cout<<"HadronPt,HadronE "<<HadronPt<<" "<<HadronE<<G4endl;
429  //G4cout<<"String4Momentum "<<String4Momentum<<G4endl;
430  //G4int Uzhi; G4cin>>Uzhi;
431  G4cout<<"Out of QGSM SplitEandP "<<G4endl;
432  #endif
433 
434  return a4Momentum;
435 }
436 
437 //-----------------------------------------------------------------------------------------
438 
439 G4bool G4QGSMFragmentation::SplitLast(G4FragmentingString * string,
440  G4KineticTrackVector * LeftVector,
441  G4KineticTrackVector * RightVector)
442 {
443  //... perform last cluster decay
444 
445  G4ThreeVector ClusterVel =string->Get4Momentum().boostVector();
446  G4double ResidualMass =string->Mass();
447 
448  #ifdef debug_QGSMfragmentation
449  G4cout<<"Split last-----------------------------------------"<<G4endl;
450  G4cout<<"StrMass "<<ResidualMass<<" q's "
451  <<string->GetLeftParton()->GetParticleName()<<" "
452  <<string->GetRightParton()->GetParticleName()<<G4endl;
453  #endif
454 
455  G4double ClusterMassCut = ClusterMass; // Taken from G4VLongitudinalStringDecay
456  G4int cClusterInterrupt = 0;
457  G4ParticleDefinition * LeftHadron, * RightHadron;
458  const G4int maxNumberOfLoops = 1000;
459  G4int loopCounter = 0;
460  do
461  {
462  if (cClusterInterrupt++ >= ClusterLoopInterrupt)
463  {
464  return false;
465  }
466 
467  G4ParticleDefinition * quark = NULL;
468  string->SetLeftPartonStable(); // to query quark contents..
469 
470  if (string->DecayIsQuark() && string->StableIsQuark() )
471  {
472  //... there are quarks on cluster ends
473  LeftHadron= QuarkSplitup(string->GetLeftParton(), quark);
474  } else {
475  //... there is a Diquark on cluster ends
476  G4int IsParticle;
477  if ( string->StableIsQuark() ) {
478  IsParticle=(string->GetLeftParton()->GetPDGEncoding()>0) ? -1 : +1;
479  } else {
480  IsParticle=(string->GetLeftParton()->GetPDGEncoding()>0) ? +1 : -1;
481  }
482 
483  pDefPair QuarkPair = CreatePartonPair(IsParticle,false); // no diquarks wanted
484  quark = QuarkPair.second;
485  LeftHadron=hadronizer->Build(QuarkPair.first, string->GetLeftParton());
486  }
487 
488  RightHadron = hadronizer->Build(string->GetRightParton(), quark);
489 
490  //... repeat procedure, if mass of cluster is too low to produce hadrons
491  //... ClusterMassCut = 0.15*GeV model parameter
492  if ( quark->GetParticleSubType()== "quark" ) {ClusterMassCut = 0.;}
493  else {ClusterMassCut = ClusterMass;}
494 
495  } while ( (ResidualMass <= LeftHadron->GetPDGMass() + RightHadron->GetPDGMass() + ClusterMassCut)
496  && ++loopCounter < maxNumberOfLoops ); /* Loop checking, 07.08.2015, A.Ribon */
497 
498  if ( loopCounter >= maxNumberOfLoops ) {
499  return false;
500  }
501 
502  //... compute hadron momenta and energies
503  G4LorentzVector LeftMom, RightMom;
505  Sample4Momentum(&LeftMom , LeftHadron->GetPDGMass() ,
506  &RightMom, RightHadron->GetPDGMass(), ResidualMass);
507  LeftMom.boost(ClusterVel);
508  RightMom.boost(ClusterVel);
509 
510  #ifdef debug_QGSMfragmentation
511  G4cout<<LeftHadron->GetParticleName()<<" "<<RightHadron->GetParticleName()<<G4endl;
512  G4cout<<"Left Hadrom P M "<<LeftMom<<" "<<LeftMom.mag()<<G4endl;
513  G4cout<<"Right Hadrom P M "<<RightMom<<" "<<RightMom.mag()<<G4endl;
514  #endif
515 
516  LeftVector->push_back(new G4KineticTrack(LeftHadron, 0, Pos, LeftMom));
517  RightVector->push_back(new G4KineticTrack(RightHadron, 0, Pos, RightMom));
518 
519  return true;
520 }
521 
522 //----------------------------------------------------------------------------------------------------------
523 
524 G4bool G4QGSMFragmentation::IsFragmentable(const G4FragmentingString * const string)
525 {
526  return sqr(FragmentationMass(string)+MassCut) < string->Mass2();
527 }
528 
529 //----------------------------------------------------------------------------------------------------------
530 
531 G4bool G4QGSMFragmentation::StopFragmenting(const G4FragmentingString * const string)
532 {
534  string->Get4Momentum().mag2();
535 }
536 
537 //----------------------------------------------------------------------------------------------------------
538 
539 void G4QGSMFragmentation::Sample4Momentum(G4LorentzVector* Mom , G4double Mass ,
540  G4LorentzVector* AntiMom, G4double AntiMass, G4double InitialMass)
541 {
542  G4double r_val = sqr(InitialMass*InitialMass - Mass*Mass - AntiMass*AntiMass) - sqr(2.*Mass*AntiMass);
543  G4double Pabs = (r_val > 0.)? std::sqrt(r_val)/(2.*InitialMass) : 0;
544 
545  //... sample unit vector
546  G4double pz = 1. - 2.*G4UniformRand();
547  G4double st = std::sqrt(1. - pz * pz)*Pabs;
548  G4double phi = 2.*pi*G4UniformRand();
549  G4double px = st*std::cos(phi);
550  G4double py = st*std::sin(phi);
551  pz *= Pabs;
552 
553  Mom->setPx(px); Mom->setPy(py); Mom->setPz(pz);
554  Mom->setE(std::sqrt(Pabs*Pabs + Mass*Mass));
555 
556  AntiMom->setPx(-px); AntiMom->setPy(-py); AntiMom->setPz(-pz);
557  AntiMom->setE (std::sqrt(Pabs*Pabs + AntiMass*AntiMass));
558 }
559 
560 //*********************************************************************************************
561 // Uzhi June 2014 Insert from G4ExcitedStringDecay.cc
562 //-----------------------------------------------------------------------------
563 
564 G4ParticleDefinition *G4QGSMFragmentation::DiQuarkSplitup(G4ParticleDefinition* decay,
565  G4ParticleDefinition *&created)
566 {
567  //... can Diquark break or not?
569  //... Diquark break
570 
571  G4int stableQuarkEncoding = decay->GetPDGEncoding()/1000;
572  G4int decayQuarkEncoding = (decay->GetPDGEncoding()/100)%10;
573 
574  if (G4UniformRand() < 0.5)
575  {
576  G4int Swap = stableQuarkEncoding;
577  stableQuarkEncoding = decayQuarkEncoding;
578  decayQuarkEncoding = Swap;
579  }
580 
581  G4int IsParticle=(decayQuarkEncoding>0) ? -1 : +1; // if we have a quark, we need antiquark)
582 
583  G4double StrSup=GetStrangeSuppress(); // for changing s-sbar production
584  StrangeSuppress=0.41; // was 0.47
585  pDefPair QuarkPair = CreatePartonPair(IsParticle,false); // no diquarks wanted
586  StrangeSuppress=StrSup;
587 
588  //... Build new Diquark
589  G4int QuarkEncoding=QuarkPair.second->GetPDGEncoding();
590  G4int i10 = std::max(std::abs(QuarkEncoding), std::abs(stableQuarkEncoding));
591  G4int i20 = std::min(std::abs(QuarkEncoding), std::abs(stableQuarkEncoding));
592  G4int spin = (i10 != i20 && G4UniformRand() <= 0.5)? 1 : 3;
593  G4int NewDecayEncoding = -1*IsParticle*(i10 * 1000 + i20 * 100 + spin);
594  created = FindParticle(NewDecayEncoding);
595  G4ParticleDefinition * decayQuark=FindParticle(decayQuarkEncoding);
596  G4ParticleDefinition * had=hadronizer->Build(QuarkPair.first, decayQuark);
597 
598  return had;
599  //return hadronizer->Build(QuarkPair.first, decayQuark);
600 
601  } else {
602  //... Diquark does not break
603 
604  G4int IsParticle=(decay->GetPDGEncoding()>0) ? +1 : -1; // if we have a diquark, we need quark)
605  G4double StrSup=GetStrangeSuppress(); // for changing s-sbar production
606  StrangeSuppress=0.41; //0.41; 0.47
607  pDefPair QuarkPair = CreatePartonPair(IsParticle,false); // no diquarks wanted
608  StrangeSuppress=StrSup;
609 
610  created = QuarkPair.second;
611 
612  G4ParticleDefinition * had=hadronizer->Build(QuarkPair.first, decay);
613  return had;
614  //return G4ParticleDefinition * had=hadronizer->Build(QuarkPair.first, decay);
615  }
616 }
617 
G4ParticleDefinition * BuildHighSpin(G4ParticleDefinition *black, G4ParticleDefinition *white)
static G4Pow * GetInstance()
Definition: G4Pow.cc:55
G4ParticleDefinition * GetRightParton(void) const
G4double powA(G4double A, G4double y) const
Definition: G4Pow.hh:259
G4ParticleDefinition * QuarkSplitup(G4ParticleDefinition *decay, G4ParticleDefinition *&created)
G4int GetPDGcode() const
Definition: G4Parton.hh:127
const double C1
const G4ThreeVector & GetPosition() const
G4ExcitedString * CPExcited(const G4ExcitedString &string)
ush Pos
Definition: deflate.h:89
static const G4double d2
void SetFormationTime(G4double aFormationTime)
void SetStrangenessSuppression(G4double aValue)
const G4LorentzVector & Get4Momentum() const
Definition: G4Parton.hh:143
G4ParticleDefinition * Build(G4ParticleDefinition *black, G4ParticleDefinition *white)
G4Parton * GetLeftParton(void) const
G4double Mass() const
G4ThreeVector SampleQuarkPt(G4double ptMax=-1.)
virtual G4KineticTrackVector * FragmentString(const G4ExcitedString &theString)
const G4String & GetParticleSubType() const
G4double FragmentationMass(const G4FragmentingString *const string, Pcreate build=0, pDefPair *pdefs=0)
int G4int
Definition: G4Types.hh:78
const G4String & GetParticleName() const
void setZ(double)
G4ParticleDefinition * GetDecayParton() const
pDefPair CreatePartonPair(G4int NeedParticle, G4bool AllowDiquarks=true)
#define NewString(str)
const G4ThreeVector & GetPosition() const
G4ParticleDefinition * GetLeftParton(void) const
void SetDiquarkBreakProbability(G4double aValue)
Hep3Vector vect() const
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
ParticleList decay(Cluster *const c)
Carries out a cluster decay.
G4KineticTrackVector * LightFragmentationTest(const G4ExcitedString *const theString)
G4LorentzVector Get4Momentum() const
G4double GetFormationTime() const
std::pair< G4ParticleDefinition *, G4ParticleDefinition * > pDefPair
double mag() const
bool G4bool
Definition: G4Types.hh:79
G4LorentzRotation TransformToAlignedCms()
HepLorentzVector & boost(double, double, double)
G4ParticleDefinition * FindParticle(G4int Encoding)
void SetPosition(const G4ThreeVector aPosition)
G4int GetDecayDirection() const
void Set4Momentum(const G4LorentzVector &a4Momentum)
static const G4double d1
G4Parton * GetRightParton(void) const
G4double GetPDGMass() const
T max(const T t1, const T t2)
brief Return the largest of the two arguments
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
double mag2() const
tuple z
Definition: test.py:28
double pz() const
G4int GetDirection(void) const
#define G4endl
Definition: G4ios.hh:61
static constexpr double pi
Definition: G4SIunits.hh:75
HepLorentzRotation inverse() const
T sqr(const T &x)
Definition: templates.hh:145
double G4double
Definition: G4Types.hh:76
const G4LorentzVector & Get4Momentum() const
void CalculateHadronTimePosition(G4double theInitialStringMass, G4KineticTrackVector *)
const G4ParticleDefinition * GetDefinition() const
G4KineticTrack * Splitup(G4FragmentingString *string, G4FragmentingString *&newString)
void SetDiquarkSuppression(G4double aValue)
CLHEP::HepLorentzVector G4LorentzVector