ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4VLongitudinalStringDecay.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4VLongitudinalStringDecay.cc
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 //
28 // -----------------------------------------------------------------------------
29 // GEANT 4 class implementation file
30 //
31 // History: first implementation, Maxim Komogorov, 1-Jul-1998
32 // redesign Gunter Folger, August/September 2001
33 // -----------------------------------------------------------------------------
34 #include "G4VLongitudinalStringDecay.hh"
35 #include "G4PhysicalConstants.hh"
36 #include "G4SystemOfUnits.hh"
37 #include "G4ios.hh"
38 #include "Randomize.hh"
39 #include "G4FragmentingString.hh"
40 
41 #include "G4ParticleDefinition.hh"
42 #include "G4ParticleTypes.hh"
43 #include "G4ParticleChange.hh"
44 #include "G4VShortLivedParticle.hh"
45 #include "G4ShortLivedConstructor.hh"
46 #include "G4ParticleTable.hh"
47 #include "G4PhaseSpaceDecayChannel.hh"
48 #include "G4VDecayChannel.hh"
49 #include "G4DecayTable.hh"
50 
51 #include "G4DiQuarks.hh"
52 #include "G4Quarks.hh"
53 #include "G4Gluons.hh"
54 
55 #include "G4Exp.hh"
56 #include "G4Log.hh"
57 
58 //------------------------debug switches
59 //#define debug_VStringDecay
60 
61 //********************************************************************************
62 // Constructors
63 
64 G4VLongitudinalStringDecay::G4VLongitudinalStringDecay() : ProbCCbar(0.0), ProbBBbar(0.0)
65 {
66  MassCut = 210.0*MeV; // Mpi + Delta
67 
68  StringLoopInterrupt = 1000;
69  ClusterLoopInterrupt = 500;
70 
71  // Changable Parameters below.
72  SigmaQT = 0.5 * GeV;
73 
74  StrangeSuppress = 0.44; // =0.27/2.27 suppresion of strange quark pait prodution, ie. u:d:s=1:1:0.27
75  DiquarkSuppress = 0.07; // Probability of qq-qqbar pair production
76  DiquarkBreakProb = 0.1; // Probability of (qq)->h+(qq)'
77 
78  //... pspin_meson is probability to create pseudo-scalar meson
79  pspin_meson = 0.5;
80 
81  //... pspin_barion is probability to create 1/2 barion
82  pspin_barion = 0.5;
83 
84  //... vectorMesonMix[] is quark mixing parameters for vector mesons (Variable spin = 3)
85  vectorMesonMix.resize(6);
86  vectorMesonMix[0] = 0.0;
87  vectorMesonMix[1] = 0.375;
88  vectorMesonMix[2] = 0.0;
89  vectorMesonMix[3] = 0.375;
90  vectorMesonMix[4] = 1.0;
91  vectorMesonMix[5] = 1.0;
92 
93  //... scalarMesonMix[] is quark mixing parameters for scalar mesons (Variable spin=1)
94  scalarMesonMix.resize(6);
95  scalarMesonMix[0] = 0.5;
96  scalarMesonMix[1] = 0.25;
97  scalarMesonMix[2] = 0.5;
98  scalarMesonMix[3] = 0.25;
99  scalarMesonMix[4] = 1.0;
100  scalarMesonMix[5] = 0.5;
101 
102  // For the time being, set to 0.0 the probabilities for c-cbar and b-bbar creation.
103  SetProbCCbar(0.0); //SetProbCCbar(0.43e-11); // Probability of CCbar pair creation
104  //Pythia8 and Pythia 6.4 Comp. Phys. Commun. 191 (2015) 159; arXiv:1410.3012
105  SetProbEta_c(0.1); // Mixing of Eta_c and J/Psi
106  SetProbBBbar(0.0); // Probability of BBbar pair creation,
107  SetProbEta_b(0.0); // Mixing of Eta_b and Ipsilon_b
108 
109  // Parameters may be changed until the first fragmentation starts
110  PastInitPhase=false;
111  hadronizer = new G4HadronBuilder(pspin_meson,pspin_barion,
112  scalarMesonMix,vectorMesonMix,
113  ProbEta_c, ProbEta_b);
114 
115  MaxMass=-350.0*GeV; // If there will be a particle with mass larger than Higgs the value must be changed.
116 
117  SetMinMasses(); // Re-calculation of minimal mass of strings and weights of particles in 2-part. decays
118 
119  Kappa = 1.0 * GeV/fermi;
120 }
121 
122 
123 G4VLongitudinalStringDecay::~G4VLongitudinalStringDecay()
124 {
125  delete hadronizer;
126 }
127 
128 //=============================================================================
129 
130 // Operators
131 
132 //-----------------------------------------------------------------------------
133 
134 G4bool G4VLongitudinalStringDecay::operator==(const G4VLongitudinalStringDecay &) const
135 {
136  throw G4HadronicException(__FILE__, __LINE__, "G4VLongitudinalStringDecay::operator== forbidden");
137  return false;
138 }
139 
140 //-------------------------------------------------------------------------------------
141 
142 G4bool G4VLongitudinalStringDecay::operator!=(const G4VLongitudinalStringDecay &) const
143 {
144  throw G4HadronicException(__FILE__, __LINE__, "G4VLongitudinalStringDecay::operator!= forbidden");
145  return true;
146 }
147 
148 //***********************************************************************************
149 
150 // For changing Mass Cut used for selection of very small mass strings
151 void G4VLongitudinalStringDecay::SetMassCut(G4double aValue){ MassCut=aValue; }
152 G4double G4VLongitudinalStringDecay::GetMassCut() { return MassCut; }
153 
154 //-----------------------------------------------------------------------------
155 
156 // For handling a string with very low mass
157 
158 G4KineticTrackVector* G4VLongitudinalStringDecay::ProduceOneHadron(const G4ExcitedString * const string)
159 {
160  G4KineticTrackVector* result = nullptr;
161  pDefPair hadrons( nullptr, nullptr );
162  G4FragmentingString aString( *string );
163 
164  #ifdef debug_VStringDecay
165  G4cout<<"G4VLongitudinalStringDecay::ProduceOneHadron: PossibleHmass StrMass "
166  <<aString.Mass()<<" MassCut "<<MassCut<<G4endl;
167  #endif
168 
169  SetMinimalStringMass( &aString );
170  PossibleHadronMass( &aString, 0, &hadrons );
171  result = new G4KineticTrackVector;
172  if ( hadrons.first != nullptr ) {
173  if ( hadrons.second == nullptr ) {
174  // Substitute string by light hadron, Note that Energy is not conserved here!
175 
176  #ifdef debug_VStringDecay
177  G4cout << "VlongSD Warning replacing string by single hadron (G4VLongitudinalStringDecay)" <<G4endl;
178  G4cout << hadrons.first->GetParticleName()<<G4endl
179  << "string .. " << string->Get4Momentum() << " "
180  << string->Get4Momentum().m() << G4endl;
181  #endif
182 
183  G4ThreeVector Mom3 = string->Get4Momentum().vect();
184  G4LorentzVector Mom( Mom3, std::sqrt( Mom3.mag2() + sqr( hadrons.first->GetPDGMass() ) ) );
185  result->push_back( new G4KineticTrack( hadrons.first, 0, string->GetPosition(), Mom ) );
186  } else {
187  //... string was qq--qqbar type: Build two stable hadrons,
188 
189  #ifdef debug_VStringDecay
190  G4cout << "VlongSD Warning replacing qq-qqbar string by TWO hadrons (G4VLongitudinalStringDecay)"
191  << hadrons.first->GetParticleName() << " / "
192  << hadrons.second->GetParticleName()
193  << "string .. " << string->Get4Momentum() << " "
194  << string->Get4Momentum().m() << G4endl;
195  #endif
196 
197  G4LorentzVector Mom1, Mom2;
198  Sample4Momentum( &Mom1, hadrons.first->GetPDGMass(),
199  &Mom2, hadrons.second->GetPDGMass(),
200  string->Get4Momentum().mag() );
201 
202  result->push_back( new G4KineticTrack( hadrons.first, 0, string->GetPosition(), Mom1 ) );
203  result->push_back( new G4KineticTrack( hadrons.second, 0, string->GetPosition(), Mom2 ) );
204 
205  G4ThreeVector Velocity = string->Get4Momentum().boostVector();
206  result->Boost(Velocity);
207  }
208  }
209  return result;
210 }
211 
212 //----------------------------------------------------------------------------------------
213 
214 G4double G4VLongitudinalStringDecay::PossibleHadronMass( const G4FragmentingString * const string,
215  Pcreate build, pDefPair * pdefs )
216 {
217  G4double mass;
218 
219  if ( build==0 ) build=&G4HadronBuilder::BuildLowSpin;
220 
221  G4ParticleDefinition *Hadron1, *Hadron2=0;
222 
223  if (!string->IsAFourQuarkString() )
224  {
225  // spin 0 meson or spin 1/2 barion will be built
226 
227  Hadron1 = (hadronizer->*build)(string->GetLeftParton(), string->GetRightParton());
228  #ifdef debug_VStringDecay
229  G4cout<<"VlongSF Quarks at the string ends "<<string->GetLeftParton()->GetParticleName()
230  <<" "<<string->GetRightParton()->GetParticleName()<<G4endl;
231  if ( Hadron1 != NULL) {
232  G4cout<<"(G4VLongitudinalStringDecay) Hadron "<<Hadron1->GetParticleName()
233  <<" "<<Hadron1->GetPDGMass()<<G4endl;
234  }
235  #endif
236  if ( Hadron1 != NULL) { mass = (Hadron1)->GetPDGMass();}
237  else { mass = MaxMass;}
238  } else
239  {
240  //... string is qq--qqbar: Build two stable hadrons,
241  //... with extra uubar or ddbar quark pair
242 
243  #ifdef debug_VStringDecay
244  G4cout<<"VlongSF string is qq--qqbar: Build two stable hadrons"<<G4endl;
245  #endif
246 
247  G4int iflc = (G4UniformRand() < 0.5)? 1 : 2;
248  if (string->GetLeftParton()->GetPDGEncoding() < 0) iflc = -iflc;
249 
250  //... theSpin = 4; spin 3/2 baryons will be built
251  Hadron1 = (hadronizer->*build)(string->GetLeftParton(), FindParticle(iflc));
252  Hadron2 = (hadronizer->*build)(string->GetRightParton(), FindParticle(-iflc));
253 
254  if ( (Hadron1 != NULL) && (Hadron2 != NULL) ) { mass = (Hadron1)->GetPDGMass() + (Hadron2)->GetPDGMass();}
255  else { mass = MaxMass;}
256  }
257 
258  #ifdef debug_VStringDecay
259  G4cout<<"VlongSF *Hadrons 1 and 2, proposed mass "<<Hadron1<<" "<<Hadron2<<" "<<mass<<G4endl;
260  #endif
261 
262  if ( pdefs != 0 )
263  { // need to return hadrons as well....
264  pdefs->first = Hadron1;
265  pdefs->second = Hadron2;
266  }
267 
268  return mass;
269 }
270 
271 //----------------------------------------------------------------------------
272 
273 G4ParticleDefinition* G4VLongitudinalStringDecay::FindParticle(G4int Encoding)
274 {
275  /*
276  G4cout<<Encoding<<" G4VLongitudinalStringDecay::FindParticle Check di-quarks *******************"<<G4endl;
277  for (G4int i=4; i<6;i++){
278  for (G4int j=1;j<6;j++){
279  G4cout<<i<<" "<<j<<" ";
280  G4int Code = 1000 * i + 100 * j +1;
281  G4ParticleDefinition* ptr1 = G4ParticleTable::GetParticleTable()->FindParticle(Code);
282  Code +=2;
283  G4ParticleDefinition* ptr2 = G4ParticleTable::GetParticleTable()->FindParticle(Code);
284  G4cout<<"Code "<<Code - 2<<" ptr "<<ptr1<<" :: Code "<<Code<<" ptr "<<ptr2<<G4endl;
285  }
286  G4cout<<G4endl;
287  }
288  */
289 
290  G4ParticleDefinition* ptr = G4ParticleTable::GetParticleTable()->FindParticle(Encoding);
291 
292  if (ptr == NULL)
293  {
294  for (size_t i=0; i < NewParticles.size(); i++)
295  {
296  if ( Encoding == NewParticles[i]->GetPDGEncoding() ) { ptr = NewParticles[i]; return ptr;}
297  }
298  }
299 
300  return ptr;
301 }
302 
303 //*********************************************************************************
304 // For decision on continue or stop string fragmentation
305 // virtual G4bool StopFragmenting(const G4FragmentingString * const string)=0;
306 // virtual G4bool IsItFragmentable(const G4FragmentingString * const string)=0;
307 //
308 // If a string can not fragment, make last break into 2 hadrons
309 // virtual G4bool SplitLast(G4FragmentingString * string,
310 // G4KineticTrackVector * LeftVector,
311 // G4KineticTrackVector * RightVector)=0;
312 //-----------------------------------------------------------------------------
313 //
314 // If a string can fragment, do the following
315 //
316 // For transver of a string to its CMS frame
317 //-----------------------------------------------------------------------------
318 
319 G4ExcitedString *G4VLongitudinalStringDecay::CopyExcited(const G4ExcitedString & in)
320 {
321  G4Parton *Left=new G4Parton(*in.GetLeftParton());
322  G4Parton *Right=new G4Parton(*in.GetRightParton());
323  return new G4ExcitedString(Left,Right,in.GetDirection());
324 }
325 
326 //-----------------------------------------------------------------------------
327 
328 G4ParticleDefinition * G4VLongitudinalStringDecay::QuarkSplitup( G4ParticleDefinition* decay,
329  G4ParticleDefinition *&created )
330 {
331  #ifdef debug_VStringDecay
332  G4cout<<"VlongSF QuarkSplitup: quark ID "<<decay->GetPDGEncoding()<<G4endl;
333  #endif
334 
335  G4int IsParticle=(decay->GetPDGEncoding()>0) ? -1 : +1; // if we have a quark, we need antiquark (or diquark)
336 
337  pDefPair QuarkPair = CreatePartonPair(IsParticle);
338  created = QuarkPair.second;
339 
340  DecayQuark = decay->GetPDGEncoding();
341  NewQuark = created->GetPDGEncoding();
342 
343  #ifdef debug_VStringDecay
344  G4cout<<"VlongSF QuarkSplitup: "<<decay->GetPDGEncoding()<<" -> "<<QuarkPair.second->GetPDGEncoding()<<G4endl;
345  G4cout<<"hadronizer->Build(QuarkPair.first, decay)"<<G4endl;
346  #endif
347  return hadronizer->Build(QuarkPair.first, decay);
348 }
349 
350 //-----------------------------------------------------------------------------
351 
352 G4VLongitudinalStringDecay::pDefPair G4VLongitudinalStringDecay::
353 CreatePartonPair(G4int NeedParticle,G4bool AllowDiquarks)
354 {
355  // NeedParticle = +1 for Particle, -1 for Antiparticle
356  if ( AllowDiquarks && G4UniformRand() < DiquarkSuppress )
357  {
358  // Create a Diquark - AntiDiquark pair , first in pair is anti to IsParticle
359  #ifdef debug_VStringDecay
360  G4cout<<"VlongSF Create a Diquark - AntiDiquark pair"<<G4endl;
361  #endif
362  G4int q1(0), q2(0), spin(0), PDGcode(0);
363 
364  q1 = SampleQuarkFlavor();
365  q2 = SampleQuarkFlavor();
366 
367  spin = (q1 != q2 && G4UniformRand() <= 0.5)? 1 : 3;
368  // convention: quark with higher PDG number is first
369  PDGcode = (std::max(q1,q2) * 1000 + std::min(q1,q2) * 100 + spin) * NeedParticle;
370 
371  return pDefPair (FindParticle(-PDGcode),FindParticle(PDGcode));
372 
373  } else {
374  // Create a Quark - AntiQuark pair, first in pair IsParticle
375  #ifdef debug_VStringDecay
376  G4cout<<"VlongSF Create a Quark - AntiQuark pair"<<G4endl;
377  #endif
378  G4int PDGcode=SampleQuarkFlavor()*NeedParticle;
379  return pDefPair (FindParticle(PDGcode),FindParticle(-PDGcode));
380  }
381 }
382 
383 //-----------------------------------------------------------------------------
384 
385 G4int G4VLongitudinalStringDecay::SampleQuarkFlavor(void)
386 {
387  G4int quark(1);
388  G4double ksi = G4UniformRand();
389  if ( ksi < ProbCB ) {
390  if ( ksi < ProbCCbar ) {quark = 4;} // c quark
391  else {quark = 5;} // b quark
392  } else {
393  quark = 1 + (int)(G4UniformRand()/StrangeSuppress);
394  }
395  #ifdef debug_VStringDecay
396  G4cout<<"VlongSF SampleQuarkFlavor "<<quark<<" (ProbCB ProbCCbar ProbBBbar "<<ProbCB
397  <<" "<<ProbCCbar<<" "<<ProbBBbar<<" )"<<G4endl;
398  #endif
399  return quark;
400 }
401 
402 //-----------------------------------------------------------------------------
403 
404 G4ThreeVector G4VLongitudinalStringDecay::SampleQuarkPt(G4double ptMax)
405 {
406  G4double Pt;
407  if ( ptMax < 0 ) {
408  // sample full gaussian
409  Pt = -G4Log(G4UniformRand());
410  } else {
411  // sample in limited range
412  Pt = -G4Log(G4RandFlat::shoot(G4Exp(-sqr(ptMax)/sqr(SigmaQT)), 1.));
413  }
414  Pt = SigmaQT * std::sqrt(Pt);
415  G4double phi = 2.*pi*G4UniformRand();
416  return G4ThreeVector(Pt * std::cos(phi),Pt * std::sin(phi),0);
417 }
418 
419 //******************************************************************************
420 
421 void G4VLongitudinalStringDecay::CalculateHadronTimePosition(G4double theInitialStringMass,
422  G4KineticTrackVector* Hadrons)
423 {
424  // `yo-yo` formation time
425  // const G4double kappa = 1.0 * GeV/fermi/4.;
426  G4double kappa = GetStringTensionParameter();
427  for (size_t c1 = 0; c1 < Hadrons->size(); c1++)
428  {
429  G4double SumPz = 0;
430  G4double SumE = 0;
431  for (size_t c2 = 0; c2 < c1; c2++)
432  {
433  SumPz += Hadrons->operator[](c2)->Get4Momentum().pz();
434  SumE += Hadrons->operator[](c2)->Get4Momentum().e();
435  }
436  G4double HadronE = Hadrons->operator[](c1)->Get4Momentum().e();
437  G4double HadronPz = Hadrons->operator[](c1)->Get4Momentum().pz();
438  Hadrons->operator[](c1)->SetFormationTime(
439  (theInitialStringMass - 2.*SumPz + HadronE - HadronPz ) / (2.*kappa) / c_light );
440  G4ThreeVector aPosition( 0, 0,
441  (theInitialStringMass - 2.*SumE - HadronE + HadronPz) / (2.*kappa) );
442  Hadrons->operator[](c1)->SetPosition(aPosition);
443  }
444 }
445 
446 //-----------------------------------------------------------------------------
447 
448 void G4VLongitudinalStringDecay::SetSigmaTransverseMomentum(G4double aValue)
449 {
450  if ( PastInitPhase ) {
451  throw G4HadronicException(__FILE__, __LINE__,
452  "G4VLongitudinalStringDecay::SetSigmaTransverseMomentum after FragmentString() not allowed");
453  } else {
454  SigmaQT = aValue;
455  }
456 }
457 
458 //----------------------------------------------------------------------------------------------------------
459 
460 void G4VLongitudinalStringDecay::SetStrangenessSuppression(G4double aValue)
461 {
462  StrangeSuppress = aValue;
463 }
464 
465 //----------------------------------------------------------------------------------------------------------
466 
467 void G4VLongitudinalStringDecay::SetDiquarkSuppression(G4double aValue)
468 {
469  DiquarkSuppress = aValue;
470 }
471 
472 //----------------------------------------------------------------------------------------
473 
474 void G4VLongitudinalStringDecay::SetDiquarkBreakProbability(G4double aValue)
475 {
476  if ( PastInitPhase ) {
477  throw G4HadronicException(__FILE__, __LINE__,
478  "G4VLongitudinalStringDecay::SetDiquarkBreakProbability after FragmentString() not allowed");
479  } else {
480  DiquarkBreakProb = aValue;
481  }
482 }
483 
484 //----------------------------------------------------------------------------------------------------------
485 
486 void G4VLongitudinalStringDecay::SetVectorMesonProbability(G4double aValue)
487 {
488  if ( PastInitPhase ) {
489  throw G4HadronicException(__FILE__, __LINE__,
490  "G4VLongitudinalStringDecay::SetVectorMesonProbability after FragmentString() not allowed");
491  } else {
492  pspin_meson = aValue;
493  delete hadronizer;
494  hadronizer = new G4HadronBuilder( pspin_meson, pspin_barion, scalarMesonMix, vectorMesonMix, ProbEta_c, ProbEta_b );
495  }
496 }
497 
498 //----------------------------------------------------------------------------------------------------------
499 
500 void G4VLongitudinalStringDecay::SetSpinThreeHalfBarionProbability(G4double aValue)
501 {
502  if ( PastInitPhase ) {
503  throw G4HadronicException(__FILE__, __LINE__,
504  "G4VLongitudinalStringDecay::SetSpinThreeHalfBarionProbability after FragmentString() not allowed");
505  } else {
506  pspin_barion = aValue;
507  delete hadronizer;
508  hadronizer = new G4HadronBuilder( pspin_meson, pspin_barion, scalarMesonMix, vectorMesonMix, ProbEta_c, ProbEta_b );
509  }
510 }
511 
512 //----------------------------------------------------------------------------------------------------------
513 
514 void G4VLongitudinalStringDecay::SetScalarMesonMixings(std::vector<G4double> aVector)
515 {
516  if ( PastInitPhase ) {
517  throw G4HadronicException(__FILE__, __LINE__,
518  "G4VLongitudinalStringDecay::SetScalarMesonMixings after FragmentString() not allowed");
519  } else {
520  if ( aVector.size() < 6 )
521  throw G4HadronicException(__FILE__, __LINE__,
522  "G4VLongitudinalStringDecay::SetScalarMesonMixings( argument Vector too small");
523  scalarMesonMix[0] = aVector[0];
524  scalarMesonMix[1] = aVector[1];
525  scalarMesonMix[2] = aVector[2];
526  scalarMesonMix[3] = aVector[3];
527  scalarMesonMix[4] = aVector[4];
528  scalarMesonMix[5] = aVector[5];
529  delete hadronizer;
530  hadronizer = new G4HadronBuilder( pspin_meson, pspin_barion, scalarMesonMix, vectorMesonMix, ProbEta_c, ProbEta_b );
531  }
532 }
533 
534 //----------------------------------------------------------------------------------------------------------
535 
536 void G4VLongitudinalStringDecay::SetVectorMesonMixings(std::vector<G4double> aVector)
537 {
538  if ( PastInitPhase ) {
539  throw G4HadronicException(__FILE__, __LINE__,
540  "G4VLongitudinalStringDecay::SetVectorMesonMixings after FragmentString() not allowed");
541  } else {
542  if ( aVector.size() < 6 )
543  throw G4HadronicException(__FILE__, __LINE__,
544  "G4VLongitudinalStringDecay::SetVectorMesonMixings( argument Vector too small");
545  vectorMesonMix[0] = aVector[0];
546  vectorMesonMix[1] = aVector[1];
547  vectorMesonMix[2] = aVector[2];
548  vectorMesonMix[3] = aVector[3];
549  vectorMesonMix[4] = aVector[4];
550  vectorMesonMix[5] = aVector[5];
551  delete hadronizer;
552  hadronizer = new G4HadronBuilder( pspin_meson, pspin_barion, scalarMesonMix, vectorMesonMix, ProbEta_c, ProbEta_b );
553  }
554 }
555 
556 //-------------------------------------------------------------------------------------------
557 
558 void G4VLongitudinalStringDecay::SetProbCCbar(G4double aValue)
559 {
560  ProbCCbar = aValue;
561  ProbCB = ProbCCbar + ProbBBbar;
562 }
563 
564 //-------------------------------------------------------------------------------------------
565 
566 void G4VLongitudinalStringDecay::SetProbEta_c(G4double aValue)
567 {
568  ProbEta_c = aValue;
569 }
570 
571 //-------------------------------------------------------------------------------------------
572 
573 void G4VLongitudinalStringDecay::SetProbBBbar(G4double aValue)
574 {
575  ProbBBbar = aValue;
576  ProbCB = ProbCCbar + ProbBBbar;
577 }
578 
579 //-------------------------------------------------------------------------------------------
580 
581 void G4VLongitudinalStringDecay::SetProbEta_b(G4double aValue)
582 {
583  ProbEta_b = aValue;
584 }
585 
586 //-------------------------------------------------------------------------------------------
587 
588 void G4VLongitudinalStringDecay::SetStringTensionParameter(G4double aValue)
589 {
590  Kappa = aValue * GeV/fermi;
591 }
592 
593 //-----------------------------------------------------------------------
594 
595 void G4VLongitudinalStringDecay::SetMinMasses()
596 {
597  // ------ For estimation of a minimal string mass ---------------
598  Mass_of_light_quark =140.*MeV;
599  Mass_of_s_quark =500.*MeV;
600  Mass_of_c_quark = 0.*MeV; // ???
601  Mass_of_b_quark = 0.*MeV; // ???
602  Mass_of_string_junction=720.*MeV;
603 
604  // ---------------- Determination of minimal mass of q-qbar strings -------------------
605  G4ParticleDefinition * hadron1; G4int Code1;
606  G4ParticleDefinition * hadron2; G4int Code2;
607  for (G4int i=1; i < 6; i++) {
608  Code1 = 100*i + 10*1 + 1;
609  hadron1 = FindParticle(Code1);
610 
611  if (hadron1 != NULL) {
612  for (G4int j=1; j < 6; j++) {
613  Code2 = 100*j + 10*1 + 1;
614  hadron2 = FindParticle(Code2);
615  if (hadron2 != NULL) {
616  minMassQQbarStr[i-1][j-1] = hadron1->GetPDGMass() + hadron2->GetPDGMass() + 70.0 * MeV;
617  }
618  }
619  }
620  }
621 
622  minMassQQbarStr[1][1] = minMassQQbarStr[0][0]; // u-ubar = 0.5 Pi0 + 0.24 Eta + 0.25 Eta'
623 
624  // ---------------- Determination of minimal mass of qq-q strings -------------------
625  G4ParticleDefinition * hadron3;
626  G4int kfla, kflb;
627  // MaxMass = -350.0*GeV; // If there will be a particle with mass larger than Higgs the value must be changed.
628 
629  for (G4int i=1; i < 6; i++) { //i=1
630  Code1 = 100*i + 10*1 + 1;
631  hadron1 = FindParticle(Code1);
632  for (G4int j=1; j < 6; j++) {
633  for (G4int k=1; k < 6; k++) {
634  kfla = std::max(j,k);
635  kflb = std::min(j,k);
636 
637  // Add d-quark
638  Code2 = 1000*kfla + 100*kflb + 10*1 + 2;
639  if ( (j == 1) && (k==1)) Code2 = 1000*2 + 100*1 + 10*1 + 2; // In the case - add u-quark.
640 
641  hadron2 = G4ParticleTable::GetParticleTable()->FindParticle(Code2);
642  hadron3 = G4ParticleTable::GetParticleTable()->FindParticle(Code2 + 2);
643 
644  if ((hadron2 == NULL) && (hadron3 == NULL)) {minMassQDiQStr[i-1][j-1][k-1] = MaxMass; continue;};
645 
646  if ((hadron2 != NULL) && (hadron3 != NULL)) {
647  if (hadron2->GetPDGMass() > hadron3->GetPDGMass() ) { hadron2 = hadron3; }
648  };
649 
650  if ((hadron2 != NULL) && (hadron3 == NULL)) {};
651 
652  if ((hadron2 == NULL) && (hadron3 != NULL)) {hadron2 = hadron3;};
653 
654  minMassQDiQStr[i-1][j-1][k-1] = hadron1->GetPDGMass() + hadron2->GetPDGMass() + 70.0 * MeV;
655  }
656  }
657  }
658 
659  // ------ An estimated minimal string mass ----------------------
660  MinimalStringMass = 0.;
661  MinimalStringMass2 = 0.;
662  // q charges d u s c b
663  Qcharge[0] = -1; Qcharge[1] = 2; Qcharge[2] = -1; Qcharge[3] = 2; Qcharge[4] = -1;
664 
665  // For treating of small string decays
666  for (G4int i=0; i<5; i++)
667  { for (G4int j=0; j<5; j++)
668  { for (G4int k=0; k<7; k++)
669  {
670  Meson[i][j][k]=0; MesonWeight[i][j][k]=0.;
671  }
672  }
673  }
674  //--------------------------
675  for (G4int i=0; i<5; i++)
676  { for (G4int j=0; j<5; j++)
677  {
678  Meson[i][j][0] = 100 * (std::max(i,j)+1) + 10 * (std::min(i,j)+1) + 1; // Scalar meson
679  MesonWeight[i][j][0] = ( pspin_meson);
680  Meson[i][j][1] = 100 * (std::max(i,j)+1) + 10 * (std::min(i,j)+1) + 3; // Vector meson
681  MesonWeight[i][j][1] = (1.-pspin_meson);
682  }
683  }
684 
685  //qqs indexes
686  //dd1 -> scalarMesonMix[0] * 111 + (1-scalarMesonMix[0]-scalarMesonMix[1]) * 221 + scalarMesonMix[1] * 331 (000)
687  //dd1 -> Pi0 Eta Eta'
688 
689  Meson[0][0][0] = 111; MesonWeight[0][0][0] = ( pspin_meson) * ( scalarMesonMix[0] ); // Pi0
690  Meson[0][0][2] = 221; MesonWeight[0][0][3] = ( pspin_meson) * (1-scalarMesonMix[0]-scalarMesonMix[1]); // Eta
691  Meson[0][0][3] = 331; MesonWeight[0][0][4] = ( pspin_meson) * ( scalarMesonMix[1]); // Eta'
692 
693  //dd3 -> vectorMesonMix[0] * 113 + (1-vectorMesonMix[0]-vectorMesonMix[1]) * 223 + vectorMesonMix[1] * 333 (001)
694  //dd3 -> rho_0 omega fi
695 
696  Meson[0][0][1] = 113; MesonWeight[0][0][1] = (1.-pspin_meson) * ( vectorMesonMix[0] ); // Rho
697  Meson[0][0][4] = 223; MesonWeight[0][0][4] = (1.-pspin_meson) * (1-vectorMesonMix[0]-vectorMesonMix[1]); // omega
698  Meson[0][0][5] = 333; MesonWeight[0][0][5] = (1.-pspin_meson) * ( vectorMesonMix[1]); // fi
699 
700  //uu1 -> scalarMesonMix[0] * 111 + (1-scalarMesonMix[0]-scalarMesonMix[1]) * 221 + scalarMesonMix[1] * 331 (110)
701  //uu1 -> Pi0 Eta Eta'
702 
703  Meson[1][1][0] = 111; MesonWeight[1][1][0] = ( pspin_meson) * ( scalarMesonMix[0] ); // Pi0
704  Meson[1][1][2] = 221; MesonWeight[1][1][2] = ( pspin_meson) * (1-scalarMesonMix[0]-scalarMesonMix[1]); // Eta
705  Meson[1][1][3] = 331; MesonWeight[1][1][3] = ( pspin_meson) * ( scalarMesonMix[1]); // Eta'
706 
707  //uu3 -> vectorMesonMix[0] * 113 + (1-vectorMesonMix[0]-vectorMesonMix[1]) * 223 + vectorMesonMix[1] * 333 (111)
708  //uu3 -> rho_0 omega fi
709 
710  Meson[1][1][1] = 113; MesonWeight[1][1][1] = (1.-pspin_meson) * ( vectorMesonMix[0] ); // Rho
711  Meson[1][1][4] = 223; MesonWeight[1][1][4] = (1.-pspin_meson) * (1-vectorMesonMix[0]-vectorMesonMix[1]); // omega
712  Meson[1][1][5] = 333; MesonWeight[1][1][5] = (1.-pspin_meson) * ( vectorMesonMix[1]); // fi
713 
714  //ss1 -> (1-scalarMesonMix[5]) * 221 + scalarMesonMix[5] * 331 (220)
715  //ss1 -> Eta Eta'
716 
717  Meson[2][2][0] = 221; MesonWeight[2][2][0] = ( pspin_meson) * (1-scalarMesonMix[5] ); // Eta
718  Meson[2][2][2] = 331; MesonWeight[2][2][2] = ( pspin_meson) * ( scalarMesonMix[5]); // Eta'
719 
720  //ss3 -> (1-vectorMesonMix[5]) * 223 + vectorMesonMix[5] * 333 (221)
721  //ss3 -> omega fi
722 
723  Meson[2][2][1] = 223; MesonWeight[2][2][1] = (1.-pspin_meson) * (1-vectorMesonMix[5] ); // omega
724  Meson[2][2][3] = 333; MesonWeight[2][2][3] = (1.-pspin_meson) * ( vectorMesonMix[5]); // fi
725 
726  //cc1 -> ProbEta_c /(1-pspin_meson) 441 (330) Probability of Eta_c
727  //cc3 -> (1-ProbEta_c)/( pspin_meson) 443 (331) Probability of J/Psi
728 
729  //bb1 -> ProbEta_b /pspin_meson 551 (440) Probability of Eta_b
730  //bb3 -> (1-ProbEta_b)/pspin_meson 553 (441) Probability of ipsilon
731 
732  if ( pspin_meson != 0. ) {
733  Meson[3][3][0] *= ( ProbEta_c)/( pspin_meson); // Eta_c
734  Meson[3][3][1] *= (1.0-ProbEta_c)/(1.-pspin_meson); // J/Psi
735 
736  Meson[4][4][0] *= ( ProbEta_b)/( pspin_meson); // Eta_b
737  Meson[4][4][1] *= (1.0-ProbEta_b)/(1.-pspin_meson); // ipsilon
738  }
739 
740  //--------------------------
741 
742  for (G4int i=0; i<5; i++)
743  { for (G4int j=0; j<5; j++)
744  { for (G4int k=0; k<5; k++)
745  { for (G4int l=0; l<4; l++)
746  { Baryon[i][j][k][l]=0; BaryonWeight[i][j][k][l]=0.;}
747  }
748  }
749  }
750 
751  kfla =0; kflb =0;
752  G4int kflc(0), kfld(0), kfle(0), kflf(0);
753  for (G4int i=0; i<5; i++)
754  { for (G4int j=0; j<5; j++)
755  { for (G4int k=0; k<5; k++)
756  {
757  kfla = i+1; kflb = j+1; kflc = k+1;
758  kfld = std::max(kfla,kflb);
759  kfld = std::max(kfld,kflc);
760 
761  kflf = std::min(kfla,kflb);
762  kflf = std::min(kflf,kflc);
763 
764  kfle = kfla + kflb + kflc - kfld - kflf;
765 
766  Baryon[i][j][k][0] = 1000 * kfld + 100 * kfle + 10 * kflf + 2; // spin=1/2
767  BaryonWeight[i][j][k][0] = ( pspin_barion);
768  Baryon[i][j][k][1] = 1000 * kfld + 100 * kfle + 10 * kflf + 4; // spin=3/2
769  BaryonWeight[i][j][k][1] = (1.-pspin_barion);
770  }
771  }
772  }
773 
774  // Delta- ddd - only 1114
775  Baryon[0][0][0][0] = 1114; BaryonWeight[0][0][0][0] = 1.0;
776  Baryon[0][0][0][1] = 0; BaryonWeight[0][0][0][1] = 0.0;
777 
778  // Delta++ uuu - only 2224
779  Baryon[1][1][1][0] = 2224; BaryonWeight[1][1][1][0] = 1.0;
780  Baryon[1][1][1][1] = 0; BaryonWeight[1][1][1][1] = 0.0;
781 
782  // Omega- sss - only 3334
783  Baryon[2][2][2][0] = 3334; BaryonWeight[2][2][2][0] = 1.0;
784  Baryon[2][2][2][1] = 0; BaryonWeight[2][2][2][1] = 0.0;
785 
786  // Omega_cc++ ccc - only 4444
787  Baryon[3][3][3][0] = 4444; BaryonWeight[3][3][3][0] = 1.0;
788  Baryon[3][3][3][1] = 0; BaryonWeight[3][3][3][1] = 0.0;
789 
790  // Omega_bb- bbb - only 5454
791  Baryon[4][4][4][0] = 5554; BaryonWeight[4][4][4][0] = 1.0;
792  Baryon[4][4][4][1] = 0; BaryonWeight[4][4][4][1] = 0.0;
793 
794  // Lambda/Sigma0 sud - 3122/3212
795  Baryon[0][1][2][0] = 3122; BaryonWeight[0][1][2][0] *= 0.5; // Lambda
796  Baryon[0][2][1][0] = 3122; BaryonWeight[0][2][1][0] *= 0.5;
797  Baryon[1][0][2][0] = 3122; BaryonWeight[1][0][2][0] *= 0.5;
798  Baryon[1][2][0][0] = 3122; BaryonWeight[1][2][0][0] *= 0.5;
799  Baryon[2][0][1][0] = 3122; BaryonWeight[2][0][1][0] *= 0.5;
800  Baryon[2][1][0][0] = 3122; BaryonWeight[2][1][0][0] *= 0.5;
801 
802  Baryon[0][1][2][2] = 3212; BaryonWeight[0][1][2][2] = 0.5 * pspin_barion; // Sigma0
803  Baryon[0][2][1][2] = 3212; BaryonWeight[0][2][1][2] = 0.5 * pspin_barion;
804  Baryon[1][0][2][2] = 3212; BaryonWeight[1][0][2][2] = 0.5 * pspin_barion;
805  Baryon[1][2][0][2] = 3212; BaryonWeight[1][2][0][2] = 0.5 * pspin_barion;
806  Baryon[2][0][1][2] = 3212; BaryonWeight[2][0][1][2] = 0.5 * pspin_barion;
807  Baryon[2][1][0][2] = 3212; BaryonWeight[2][1][0][2] = 0.5 * pspin_barion;
808 
809  // Lambda_c+/Sigma_c+ cud - 4122/4212
810  Baryon[0][1][3][0] = 4122; BaryonWeight[0][1][3][0] *= 0.5; // Lambda_c+
811  Baryon[0][3][1][0] = 4122; BaryonWeight[0][3][1][0] *= 0.5;
812  Baryon[1][0][3][0] = 4122; BaryonWeight[1][0][3][0] *= 0.5;
813  Baryon[1][3][0][0] = 4122; BaryonWeight[1][3][0][0] *= 0.5;
814  Baryon[3][0][1][0] = 4122; BaryonWeight[3][0][1][0] *= 0.5;
815  Baryon[3][1][0][0] = 4122; BaryonWeight[3][1][0][0] *= 0.5;
816 
817  Baryon[0][1][3][2] = 4212; BaryonWeight[0][1][3][2] = 0.5 * pspin_barion; // SigmaC+
818  Baryon[0][3][1][2] = 4212; BaryonWeight[0][3][1][2] = 0.5 * pspin_barion;
819  Baryon[1][0][3][2] = 4212; BaryonWeight[1][0][3][2] = 0.5 * pspin_barion;
820  Baryon[1][3][0][2] = 4212; BaryonWeight[1][3][0][2] = 0.5 * pspin_barion;
821  Baryon[3][0][1][2] = 4212; BaryonWeight[3][0][1][2] = 0.5 * pspin_barion;
822  Baryon[3][1][0][2] = 4212; BaryonWeight[3][1][0][2] = 0.5 * pspin_barion;
823 
824  // Xi_c+/Xi_c+' cus - 4232/4322
825  Baryon[1][2][3][0] = 4232; BaryonWeight[1][2][3][0] *= 0.5; // Xi_c+
826  Baryon[1][3][2][0] = 4232; BaryonWeight[1][3][2][0] *= 0.5;
827  Baryon[2][1][3][0] = 4232; BaryonWeight[2][1][3][0] *= 0.5;
828  Baryon[2][3][1][0] = 4232; BaryonWeight[2][3][1][0] *= 0.5;
829  Baryon[3][1][2][0] = 4232; BaryonWeight[3][1][2][0] *= 0.5;
830  Baryon[3][2][1][0] = 4232; BaryonWeight[3][2][1][0] *= 0.5;
831 
832  Baryon[1][2][3][2] = 4322; BaryonWeight[1][2][3][2] = 0.5 * pspin_barion; // Xi_c+'
833  Baryon[1][3][2][2] = 4322; BaryonWeight[1][3][2][2] = 0.5 * pspin_barion;
834  Baryon[2][1][3][2] = 4322; BaryonWeight[2][1][3][2] = 0.5 * pspin_barion;
835  Baryon[2][3][1][2] = 4322; BaryonWeight[2][3][1][2] = 0.5 * pspin_barion;
836  Baryon[3][1][2][2] = 4322; BaryonWeight[3][1][2][2] = 0.5 * pspin_barion;
837  Baryon[3][2][1][2] = 4322; BaryonWeight[3][2][1][2] = 0.5 * pspin_barion;
838 
839  // Xi_c0/Xi_c0' cus - 4232/4322
840  Baryon[0][2][3][0] = 4132; BaryonWeight[0][2][3][0] *= 0.5; // Xi_c0
841  Baryon[0][3][2][0] = 4132; BaryonWeight[0][3][2][0] *= 0.5;
842  Baryon[2][0][3][0] = 4132; BaryonWeight[2][0][3][0] *= 0.5;
843  Baryon[2][3][0][0] = 4132; BaryonWeight[2][3][0][0] *= 0.5;
844  Baryon[3][0][2][0] = 4132; BaryonWeight[3][0][2][0] *= 0.5;
845  Baryon[3][2][0][0] = 4132; BaryonWeight[3][2][0][0] *= 0.5;
846 
847  Baryon[0][2][3][2] = 4312; BaryonWeight[0][2][3][2] = 0.5 * pspin_barion; // Xi_c0'
848  Baryon[0][3][2][2] = 4312; BaryonWeight[0][3][2][2] = 0.5 * pspin_barion;
849  Baryon[2][0][3][2] = 4312; BaryonWeight[2][0][3][2] = 0.5 * pspin_barion;
850  Baryon[2][3][0][2] = 4312; BaryonWeight[2][3][0][2] = 0.5 * pspin_barion;
851  Baryon[3][0][2][2] = 4312; BaryonWeight[3][0][2][2] = 0.5 * pspin_barion;
852  Baryon[3][2][0][2] = 4312; BaryonWeight[3][2][0][2] = 0.5 * pspin_barion;
853 
854  // Lambda_b0/Sigma_b0 bud - 5122/5212
855  Baryon[0][1][4][0] = 5122; BaryonWeight[0][1][4][0] *= 0.5; // Lambda_b0
856  Baryon[0][4][1][0] = 5122; BaryonWeight[0][4][1][0] *= 0.5;
857  Baryon[1][0][4][0] = 5122; BaryonWeight[1][0][4][0] *= 0.5;
858  Baryon[1][4][0][0] = 5122; BaryonWeight[1][4][0][0] *= 0.5;
859  Baryon[4][0][1][0] = 5122; BaryonWeight[4][0][1][0] *= 0.5;
860  Baryon[4][1][0][0] = 5122; BaryonWeight[4][1][0][0] *= 0.5;
861 
862  Baryon[0][1][4][2] = 5212; BaryonWeight[0][1][4][2] = 0.5 * pspin_barion; // Sigma_b0
863  Baryon[0][4][1][2] = 5212; BaryonWeight[0][4][1][2] = 0.5 * pspin_barion;
864  Baryon[1][0][4][2] = 5212; BaryonWeight[1][0][4][2] = 0.5 * pspin_barion;
865  Baryon[1][4][0][2] = 5212; BaryonWeight[1][4][0][2] = 0.5 * pspin_barion;
866  Baryon[4][0][1][2] = 5212; BaryonWeight[4][0][1][2] = 0.5 * pspin_barion;
867  Baryon[4][1][0][2] = 5212; BaryonWeight[4][1][0][2] = 0.5 * pspin_barion;
868 
869  // Xi_b-/Xi_b-' bus - 5232/5322
870  Baryon[1][2][4][0] = 5232; BaryonWeight[1][2][4][0] *= 0.5; // Xi_b-
871  Baryon[1][4][2][0] = 5232; BaryonWeight[1][4][2][0] *= 0.5;
872  Baryon[2][1][4][0] = 5232; BaryonWeight[2][1][4][0] *= 0.5;
873  Baryon[2][4][1][0] = 5232; BaryonWeight[2][4][1][0] *= 0.5;
874  Baryon[4][1][2][0] = 5232; BaryonWeight[4][1][2][0] *= 0.5;
875  Baryon[4][2][1][0] = 5232; BaryonWeight[4][2][1][0] *= 0.5;
876 
877  Baryon[1][2][4][2] = 5322; BaryonWeight[1][2][4][2] = 0.5 * pspin_barion; // Xi_b-'
878  Baryon[1][4][2][2] = 5322; BaryonWeight[1][4][2][2] = 0.5 * pspin_barion;
879  Baryon[2][1][4][2] = 5322; BaryonWeight[2][1][4][2] = 0.5 * pspin_barion;
880  Baryon[2][4][1][2] = 5322; BaryonWeight[2][4][1][2] = 0.5 * pspin_barion;
881  Baryon[4][1][2][2] = 5322; BaryonWeight[4][1][2][2] = 0.5 * pspin_barion;
882  Baryon[4][2][1][2] = 5322; BaryonWeight[4][2][1][2] = 0.5 * pspin_barion;
883 
884  // Xi_b0/Xi_b0' bus - 5232/5322
885  Baryon[0][2][4][0] = 5132; BaryonWeight[0][2][4][0] *= 0.5; // Xi_b0
886  Baryon[0][4][2][0] = 5132; BaryonWeight[0][4][2][0] *= 0.5;
887  Baryon[2][0][4][0] = 5132; BaryonWeight[2][0][4][0] *= 0.5;
888  Baryon[2][4][0][0] = 5132; BaryonWeight[2][4][0][0] *= 0.5;
889  Baryon[4][0][2][0] = 5132; BaryonWeight[4][0][2][0] *= 0.5;
890  Baryon[4][2][0][0] = 5132; BaryonWeight[4][2][0][0] *= 0.5;
891 
892  Baryon[0][2][4][2] = 5312; BaryonWeight[0][2][4][2] = 0.5 * pspin_barion; // Xi_b0'
893  Baryon[0][4][2][2] = 5312; BaryonWeight[0][4][2][2] = 0.5 * pspin_barion;
894  Baryon[2][0][4][2] = 5312; BaryonWeight[2][0][4][2] = 0.5 * pspin_barion;
895  Baryon[2][4][0][2] = 5312; BaryonWeight[2][4][0][2] = 0.5 * pspin_barion;
896  Baryon[4][0][2][2] = 5312; BaryonWeight[4][0][2][2] = 0.5 * pspin_barion;
897  Baryon[4][2][0][2] = 5312; BaryonWeight[4][2][0][2] = 0.5 * pspin_barion;
898 
899  for (G4int i=0; i<5; i++)
900  { for (G4int j=0; j<5; j++)
901  { for (G4int k=0; k<5; k++)
902  { for (G4int l=0; l<4; l++)
903  {
904  G4ParticleDefinition * TestHadron=
905  G4ParticleTable::GetParticleTable()->FindParticle(Baryon[i][j][k][l]);
906  /*
907  G4cout<<i<<" "<<j<<" "<<k<<" "<<l<<" "<<Baryon[i][j][k][l]<<" "<<TestHadron<<" "<<BaryonWeight[i][j][k][l];
908  if (TestHadron != NULL) G4cout<<" "<<TestHadron->GetParticleName();
909  if ((TestHadron == NULL)&&(Baryon[i][j][k][l] != 0)) G4cout<<" *****";
910  if ((TestHadron == NULL)&&(Baryon[i][j][k][l] == 0)) G4cout<<" ---------------";
911  G4cout<<G4endl;
912  */
913  if ((TestHadron == NULL)&&(Baryon[i][j][k][l] != 0)) Baryon[i][j][k][l] = 0;
914  }
915  }
916  }
917  }
918 
919  // --------- Probabilities of q-qbar pair productions for kink or gluons.
920  G4double ProbUUbar = 0.33;
921  Prob_QQbar[0]=ProbUUbar; // Probability of ddbar production
922  Prob_QQbar[1]=ProbUUbar; // Probability of uubar production
923  Prob_QQbar[2]=1.0-2.*ProbUUbar; // Probability of ssbar production
924  Prob_QQbar[3]=0.0; // Probability of ccbar production
925  Prob_QQbar[4]=0.0; // Probability of bbbar production
926 
927  for ( G4int i=0 ; i<350 ; i++ ) { // Must be checked
928  FS_LeftHadron[i] = 0;
929  FS_RightHadron[i] = 0;
930  FS_Weight[i] = 0.0;
931  }
932 
933  NumberOf_FS = 0;
934 }
935 
936 // --------------------------------------------------------------
937 
938 void G4VLongitudinalStringDecay::SetMinimalStringMass(const G4FragmentingString * const string)
939 {
940  //MaxMass = -350.0*GeV;
941  G4double EstimatedMass=0.;
942 
943  G4int Qleft =std::abs(string->GetLeftParton()->GetPDGEncoding());
944  G4int Qright=std::abs(string->GetRightParton()->GetPDGEncoding());
945 
946  if ((Qleft < 6) && (Qright < 6)) { // Q-Qbar string
947  EstimatedMass=minMassQQbarStr[Qleft-1][Qright-1];
948  MinimalStringMass=EstimatedMass;
949  SetMinimalStringMass2(EstimatedMass);
950  return;
951  }
952 
953  if ((Qleft < 6) && (Qright > 1000)) { // Q - DiQ string
954  G4int q1=Qright/1000;
955  G4int q2=(Qright/100)%10;
956  EstimatedMass=minMassQDiQStr[Qleft-1][q1-1][q2-1];
957  MinimalStringMass=EstimatedMass; // It can be negative!
958  SetMinimalStringMass2(EstimatedMass);
959  return;
960  }
961 
962  if ((Qleft > 1000) && (Qright < 6)) { // DiQ - Q string 6 6 6
963  G4int q1=Qleft/1000;
964  G4int q2=(Qleft/100)%10;
965  EstimatedMass=minMassQDiQStr[Qright-1][q1-1][q2-1];
966  MinimalStringMass=EstimatedMass; // It can be negative!
967  SetMinimalStringMass2(EstimatedMass);
968  return;
969  }
970 
971  // DiQuark - Anti DiQuark string -----------------
972 
973  G4double StringM=string->Get4Momentum().mag();
974 
975  #ifdef debug_LUNDfragmentation
976  // G4cout<<"MinStringMass// Input String mass "<<string->Get4Momentum().mag()<<" Qleft "<<Qleft<<G4endl;
977  #endif
978 
979  G4int q1= Qleft/1000 ;
980  G4int q2=(Qleft/100)%10 ;
981 
982  G4int q3= Qright/1000 ;
983  G4int q4=(Qright/100)%10;
984 
985  // -------------- 2 baryon production or 2 mesons production --------
986 
987  G4double EstimatedMass1 = minMassQDiQStr[q1-1][q2-1][0];
988  G4double EstimatedMass2 = minMassQDiQStr[q3-1][q4-1][0];
989  // Mass is negative if there is no corresponding particle.
990 
991  if ( (EstimatedMass1 > 0.) && (EstimatedMass2 > 0.)) {
992  EstimatedMass = EstimatedMass1 + EstimatedMass2;
993  if ( StringM > EstimatedMass ) { // 2 baryon production is possible.
994  MinimalStringMass=EstimatedMass1 + EstimatedMass2;
995  SetMinimalStringMass2(EstimatedMass);
996  return;
997  }
998  }
999 
1000  if ( (EstimatedMass1 < 0.) && (EstimatedMass2 > 0.)) {
1001  EstimatedMass = MaxMass;
1002  MinimalStringMass=EstimatedMass;
1003  SetMinimalStringMass2(EstimatedMass);
1004  return;
1005  }
1006 
1007  if ( (EstimatedMass1 > 0.) && (EstimatedMass2 < 0.)) {
1008  EstimatedMass = EstimatedMass1;
1009  MinimalStringMass=EstimatedMass;
1010  SetMinimalStringMass2(EstimatedMass);
1011  return;
1012  }
1013 
1014  // if ( EstimatedMass >= StringM ) {
1015  // ------------- Re-orangement ---------------
1016  EstimatedMass=std::min(minMassQQbarStr[q1-1][q3-1] + minMassQQbarStr[q2-1][q4-1],
1017  minMassQQbarStr[q1-1][q4-1] + minMassQQbarStr[q2-1][q3-1]);
1018 
1019  // In principle, re-orangement and 2 baryon production can compite.
1020  // More physics consideration is needed.
1021 
1022  MinimalStringMass=EstimatedMass;
1023  SetMinimalStringMass2(EstimatedMass);
1024 
1025  return;
1026 }
1027 
1028 //--------------------------------------------------------------------------------------
1029 
1030 void G4VLongitudinalStringDecay::SetMinimalStringMass2(const G4double aValue)
1031 {
1032  MinimalStringMass2=aValue * aValue;
1033 }
1034