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