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