| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 //
27 //
28 // -------------------------------------------
29 // GEANT 4 class implementation file
30 //
31 // History: first implementation, Maxim K
32 // redesign Gunter Folger, Augu
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
67 G4VLongitudinalStringDecay::G4VLongitudinalStr
68 : G4HadronicInteraction(name), ProbCCbar(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 s
79 DiquarkSuppress = 0.07; // Probability
80 DiquarkBreakProb = 0.1; // Probability
81
82 //... pspin_meson is probability to create
83 pspin_meson.resize(3);
84 pspin_meson[0] = 0.5; // u or d + anti-u o
85 pspin_meson[1] = 0.4; // one of the quark
86 pspin_meson[2] = 0.3; // both of the quark
87
88 //... pspin_barion is probability to create
89 pspin_barion = 0.5;
90
91 //... vectorMesonMix[] is quark mixing para
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 para
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
110 SetProbEta_c(0.1); // Mixing of Eta_c and
111 SetProbBBbar(0.0); // Probability of BBbar
112 SetProbEta_b(0.0); // Mixing of Eta_b and
113
114 // Parameters may be changed until the firs
115 PastInitPhase=false;
116 hadronizer = new G4HadronBuilder( pspin_mes
117 ProbEta_c
118
119 MaxMass=-350.0*GeV; // If there will be a
120
121 SetMinMasses(); // Re-calculation of minim
122
123 Kappa = 1.0 * GeV/fermi;
124 DecayQuark = NewQuark = 0;
125 }
126
127 G4VLongitudinalStringDecay::~G4VLongitudinalSt
128 {
129 delete hadronizer;
130 }
131
132 G4HadFinalState*
133 G4VLongitudinalStringDecay::ApplyYourself(cons
134 {
135 return nullptr;
136 }
137
138 //============================================
139
140 // For changing Mass Cut used for selection of
141 void G4VLongitudinalStringDecay::SetMassCut(G4
142 G4double G4VLongitudinalStringDecay::GetMassCu
143
144 //--------------------------------------------
145
146 // For handling a string with very low mass
147
148 G4KineticTrackVector* G4VLongitudinalStringDec
149 {
150 G4KineticTrackVector* result = nullptr
151 pDefPair hadrons( nullptr, nullptr );
152 G4FragmentingString aString( *string )
153
154 #ifdef debug_VStringDecay
155 G4cout<<"G4VLongitudinalStringDecay::P
156 <<aString.Mass()<<" MassCut "<<M
157 #endif
158
159 SetMinimalStringMass( &aString );
160 PossibleHadronMass( &aString, 0, &hadr
161 result = new G4KineticTrackVector;
162 if ( hadrons.first != nullptr ) {
163 if ( hadrons.second == nullptr ) {
164 // Substitute string by light h
165
166 #ifdef debug_VStringDecay
167 G4cout << "VlongSD Warning repl
168 G4cout << hadrons.first->GetPar
169 << "string .. " << strin
170 << string->Get4Momentum(
171 #endif
172
173 G4ThreeVector Mom3 = string->
174 G4LorentzVector Mom( Mom3, std:
175 result->push_back( new G4Kineti
176 } else {
177 //... string was qq--qqbar type
178
179 #ifdef debug_VStringDecay
180 G4cout << "VlongSD Warning repl
181 << hadrons.first->GetPar
182 << hadrons.second->GetPa
183 << "string .. " << strin
184 << string->Get4Momentum(
185 #endif
186
187 G4LorentzVector Mom1, Mom2;
188 Sample4Momentum( &Mom1, hadrons
189 &Mom2, hadrons
190 string->Get4Mo
191
192 result->push_back( new G4Kineti
193 result->push_back( new G4Kineti
194
195 G4ThreeVector Velocity = string
196 result->Boost(Velocity);
197 }
198 }
199 return result;
200 }
201
202 //--------------------------------------------
203
204 G4double G4VLongitudinalStringDecay::PossibleH
205
206 {
207 G4double mass = 0.0;
208
209 if ( build==0 ) build=&G4HadronBuilder::Buil
210
211 G4ParticleDefinition* Hadron1 = nullpt
212 G4ParticleDefinition* Hadron2 = nullptr;
213
214 if (!string->IsAFourQuarkString() )
215 {
216 // spin 0 meson or spin 1/2 barion
217
218 Hadron1 = (hadronizer->*build)(stri
219 #ifdef debug_VStringDecay
220 G4cout<<"VlongSD PossibleHadronMass"<<G4e
221 G4cout<<"VlongSD Quarks at the stri
222 <<" "<<string->GetRightParton
223 if ( Hadron1 != nullptr) {
224 G4cout<<"(G4VLongitudinalStringDe
225 <<" "<<Hadron1->GetPDGMass(
226 }
227 #endif
228 if ( Hadron1 != nullptr) { mass = (
229 else { mass = MaxM
230 } else
231 {
232 //... string is qq--qqbar: Build tw
233
234 #ifdef debug_VStringDecay
235 G4cout<<"VlongSD PossibleHadronMass
236 G4cout<<"VlongSD string is qq--qqba
237 #endif
238
239 G4double StringMass = string->Mas
240 G4int cClusterInterrupt = 0;
241 do
242 {
243 if (cClusterInterrupt++ >= Cluste
244
245 G4int LeftQuark1= string->GetLeft
246 G4int LeftQuark2=(string->GetLeft
247
248 G4int RightQuark1= string->GetRig
249 G4int RightQuark2=(string->GetRig
250
251 if (G4UniformRand()<0.5) {
252 Hadron1 =hadronizer->Build(Find
253 Hadron2 =hadronizer->Build(Find
254 } else {
255 Hadron1 =hadronizer->Build(Find
256 Hadron2 =hadronizer->Build(Find
257 }
258 //... repeat procedure, if mass o
259 //... ClusterMassCut = 0.15*GeV m
260 }
261 while ( Hadron1 == nullptr || Hadro
262 ( StringMass <= Hadron1->Ge
263
264 mass = (Hadron1)->GetPDGMass() + (Hadron2
265 }
266
267 #ifdef debug_VStringDecay
268 G4cout<<"VlongSD *Hadrons 1 and 2, pro
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
282 G4ParticleDefinition* G4VLongitudinalStringDec
283 {
284 /*
285 G4cout<<Encoding<<" G4VLongitudinalStringDec
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 = G4ParticleT
291 Code +=2;
292 G4ParticleDefinition* ptr2 = G4ParticleT
293 G4cout<<"Code "<<Code - 2<<" ptr "<<ptr1
294 }
295 G4cout<<G4endl;
296 }
297 */
298
299 G4ParticleDefinition* ptr = G4ParticleTable:
300
301 if (ptr == nullptr)
302 {
303 for (size_t i=0; i < NewParticles.size();
304 {
305 if ( Encoding == NewParticles[i]->GetPD
306 }
307 }
308
309 return ptr;
310 }
311
312 //********************************************
313 // For decision on continue or stop string f
314 // virtual G4bool StopFragmenting(const G4Fr
315 // virtual G4bool IsItFragmentable(const G4F
316 //
317 // If a string can not fragment, make last b
318 // virtual G4bool SplitLast(G4FragmentingStr
319 // G4KineticTrackVe
320 // G4KineticTrackVe
321 //--------------------------------------------
322 //
323 // If a string can fragment, do the followin
324 //
325 // For transver of a string to its CMS frame
326 //--------------------------------------------
327
328 G4ExcitedString *G4VLongitudinalStringDecay::C
329 {
330 G4Parton *Left=new G4Parton(*in.GetLeftParto
331 G4Parton *Right=new G4Parton(*in.GetRightPar
332 return new G4ExcitedString(Left,Right,in.Get
333 }
334
335 //--------------------------------------------
336
337 G4ParticleDefinition * G4VLongitudinalStringDe
338
339 {
340 #ifdef debug_VStringDecay
341 G4cout<<"VlongSD QuarkSplitup: quark ID "<<
342 #endif
343
344 G4int IsParticle=(decay->GetPDGEncoding()>0
345
346 pDefPair QuarkPair = CreatePartonPair(IsPar
347 created = QuarkPair.second;
348
349 DecayQuark = decay->GetPDGEncoding();
350 NewQuark = created->GetPDGEncoding();
351
352 #ifdef debug_VStringDecay
353 G4cout<<"VlongSD QuarkSplitup: "<<decay->Ge
354 G4cout<<"hadronizer->Build(QuarkPair.first,
355 #endif
356
357 return hadronizer->Build(QuarkPair.first, d
358 }
359
360 //--------------------------------------------
361
362 G4VLongitudinalStringDecay::pDefPair G4VLongit
363 CreatePartonPair(G4int NeedParticle,G4bool All
364 {
365 // NeedParticle = +1 for Particle, -1 for
366 if ( AllowDiquarks && G4UniformRand() < Di
367 {
368 // Create a Diquark - AntiDiquark pair ,
369 #ifdef debug_VStringDecay
370 G4cout<<"VlongSD Create a Diquark - Anti
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
378 // conv
379 PDGcode = (std::max(q1,q2) * 1000 + std:
380
381 return pDefPair (FindParticle(-PDGcode),
382
383 } else {
384 // Create a Quark - AntiQuark pair, firs
385 #ifdef debug_VStringDecay
386 G4cout<<"VlongSD Create a Quark - AntiQu
387 #endif
388 G4int PDGcode=SampleQuarkFlavor()*NeedPa
389 return pDefPair (FindParticle(PDGcode),F
390 }
391 }
392
393 //--------------------------------------------
394
395 G4int G4VLongitudinalStringDecay::SampleQuarkF
396 {
397 G4int quark(1);
398 G4double ksi = G4UniformRand();
399 if ( ksi < ProbCB ) {
400 if ( ksi < ProbCCbar ) {quark = 4;} //
401 else {quark = 5;} //
402 #ifdef debug_heavyHadrons
403 G4cout << "G4VLongitudinalStringDecay::S
404 << quark << G4endl;
405 #endif
406 } else {
407 quark = 1 + (int)(G4UniformRand()/Strange
408 }
409 #ifdef debug_VStringDecay
410 G4cout<<"VlongSD SampleQuarkFlavor "<<quark
411 <<" "<<ProbCCbar<<" "<<ProbBBbar<<" )
412 #endif
413 return quark;
414 }
415
416 //--------------------------------------------
417
418 G4ThreeVector G4VLongitudinalStringDecay::Samp
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(
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
433 }
434
435 //********************************************
436
437 void G4VLongitudinalStringDecay::CalculateHadr
438
439 {
440 // `yo-yo` formation time
441 // const G4double kappa = 1.0 * GeV/fermi
442 G4double kappa = GetStringTensionParameter(
443 for (size_t c1 = 0; c1 < Hadrons->size(); c
444 {
445 G4double SumPz = 0;
446 G4double SumE = 0;
447 for (size_t c2 = 0; c2 < c1; c2++)
448 {
449 SumPz += Hadrons->operator[](c2)->Get
450 SumE += Hadrons->operator[](c2)->Get
451 }
452 G4double HadronE = Hadrons->operator[](
453 G4double HadronPz = Hadrons->operator[](
454 Hadrons->operator[](c1)->SetFormationTim
455 (theInitialStringMass - 2.*SumPz + Had
456 G4ThreeVector aPosition( 0, 0,
457 (theInitialStringMass - 2.*SumE - Had
458 Hadrons->operator[](c1)->SetPosition(aPo
459 }
460 }
461
462 //--------------------------------------------
463
464 void G4VLongitudinalStringDecay::SetSigmaTrans
465 {
466 if ( PastInitPhase ) {
467 throw G4HadronicException(__FILE__, __LIN
468 "G4VLongitudinalStringDecay::SetSigmaTr
469 } else {
470 SigmaQT = aValue;
471 }
472 }
473
474 //--------------------------------------------
475
476 void G4VLongitudinalStringDecay::SetStrangenes
477 {
478 StrangeSuppress = aValue;
479 }
480
481 //--------------------------------------------
482
483 void G4VLongitudinalStringDecay::SetDiquarkSup
484 {
485 DiquarkSuppress = aValue;
486 }
487
488 //--------------------------------------------
489
490 void G4VLongitudinalStringDecay::SetDiquarkBre
491 {
492 if ( PastInitPhase ) {
493 throw G4HadronicException(__FILE__, __LINE
494 "G4VLongitudinalStringDecay::SetDiquarkB
495 } else {
496 DiquarkBreakProb = aValue;
497 }
498 }
499
500 //--------------------------------------------
501
502 void G4VLongitudinalStringDecay::SetSpinThreeH
503 {
504 if ( PastInitPhase ) {
505 throw G4HadronicException(__FILE__, __LINE
506 "G4VLongitudinalStringDecay::SetSpinThre
507 } else {
508 pspin_barion = aValue;
509 delete hadronizer;
510 hadronizer = new G4HadronBuilder( pspin_me
511 ProbEta_
512 }
513 }
514
515 //--------------------------------------------
516
517 void G4VLongitudinalStringDecay::SetScalarMeso
518 {
519 if ( PastInitPhase ) {
520 throw G4HadronicException(__FILE__, __LINE
521 "G4VLongitudinalStringDecay::SetScalarMe
522 } else {
523 if ( aVector.size() < 6 )
524 throw G4HadronicException(__FILE__, __LI
525 "G4VLongitudinalStringDecay::SetScalar
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_me
534 ProbEta_
535 }
536 }
537
538 //--------------------------------------------
539
540 void G4VLongitudinalStringDecay::SetVectorMeso
541 {
542 if ( PastInitPhase ) {
543 throw G4HadronicException(__FILE__, __LINE
544 "G4VLongitudinalStringDecay::SetVectorMe
545 } else {
546 if ( aVector.size() < 6 )
547 throw G4HadronicException(__FILE__, __LI
548 "G4VLongitudinalStringDecay::SetVector
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_me
557 ProbEta_
558 }
559 }
560
561 //--------------------------------------------
562
563 void G4VLongitudinalStringDecay::SetProbCCbar(
564 {
565 ProbCCbar = aValue;
566 ProbCB = ProbCCbar + ProbBBbar;
567 }
568
569 //--------------------------------------------
570
571 void G4VLongitudinalStringDecay::SetProbEta_c(
572 {
573 ProbEta_c = aValue;
574 }
575
576 //--------------------------------------------
577
578 void G4VLongitudinalStringDecay::SetProbBBbar(
579 {
580 ProbBBbar = aValue;
581 ProbCB = ProbCCbar + ProbBBbar;
582 }
583
584 //--------------------------------------------
585
586 void G4VLongitudinalStringDecay::SetProbEta_b(
587 {
588 ProbEta_b = aValue;
589 }
590
591 //--------------------------------------------
592
593 void G4VLongitudinalStringDecay::SetStringTens
594 {
595 Kappa = aValue * GeV/fermi;
596 }
597
598 //--------------------------------------------
599
600 void G4VLongitudinalStringDecay::SetMinMasses(
601 {
602 // ------ For estimation of a minimal stri
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 minim
610 G4ParticleDefinition * hadron1; G4int C
611 G4ParticleDefinition * hadron2; G4int C
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] = h
622 }
623 }
624 }
625 }
626
627 minMassQQbarStr[1][1] = minMassQQbarStr[0]
628
629 // ---------------- Determination of minim
630 G4ParticleDefinition * hadron3;
631 G4int kfla, kflb;
632 // MaxMass = -350.0*GeV; // If there wi
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 +
644 if ( (j == 1) && (k==1)) Code2 = 1000*2 +
645
646 hadron2 = G4ParticleTable::Get
647 hadron3 = G4ParticleTable::Get
648
649 if ((hadron2 == nullptr) && (h
650
651 if ((hadron2 != nullptr) && (h
652 if (hadron2->GetPDGMass() >
653 };
654
655 if ((hadron2 != nullptr) && (h
656
657 if ((hadron2 == nullptr) && (h
658
659 minMassQDiQStr[i-1][j-1][k-1]
660 }
661 }
662 }
663
664 // ------ An estimated minimal string mass
665 MinimalStringMass = 0.;
666 MinimalStringMass2 = 0.;
667 // q charges d u
668 Qcharge[0] = -1; Qcharge[1] = 2; Qcharge[2
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
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::ma
690 MesonWeight[i][j][0] = ( pspin_meso
691 Meson[i][j][1] = 100 * (std::ma
692 MesonWeight[i][j][1] = (1.-pspin_meso
693 }
694 }
695
696 //qqs
697 //dd1 -> scalarMesonMix[0] * 111 + (1-scal
698 //dd1 -> Pi0
699
700 Meson[0][0][0] = 111; MesonWeight[0][0][0]
701 Meson[0][0][2] = 221; MesonWeight[0][0][3]
702 Meson[0][0][3] = 331; MesonWeight[0][0][4]
703
704 //dd3 -> (1-vectorMesonMix[1] * 113 + vect
705 //dd3 -> rho_0
706
707 Meson[0][0][1] = 113; MesonWeight[0][0][1]
708 Meson[0][0][4] = 223; MesonWeight[0][0][4]
709
710 //uu1 -> scalarMesonMix[0] * 111 + (1-scal
711 //uu1 -> Pi0
712
713 Meson[1][1][0] = 111; MesonWeight[1][1][0]
714 Meson[1][1][2] = 221; MesonWeight[1][1][2]
715 Meson[1][1][3] = 331; MesonWeight[1][1][3]
716
717 //uu3 -> (1-vectorMesonMix[1]) * 113 + vec
718 //uu3 -> rho_0
719
720 Meson[1][1][1] = 113; MesonWeight[1][1][1]
721 Meson[1][1][4] = 223; MesonWeight[1][1][4]
722
723 //ss1 ->
724 //ss1 ->
725
726 Meson[2][2][0] = 221; MesonWeight[2][2][0]
727 Meson[2][2][2] = 331; MesonWeight[2][2][2]
728
729 //ss3 ->
730 //ss3 ->
731
732 Meson[2][2][1] = 333; MesonWeight[2][2][1]
733
734 //cc1 -> ProbEta_c /(1-pspin_meson)
735 //cc3 -> (1-ProbEta_c)/( pspin_meson)
736
737 //bb1 -> ProbEta_b /pspin_meson 551
738 //bb3 -> (1-ProbEta_b)/pspin_meson 553
739
740 if ( pspin_meson[2] != 0. ) {
741 Meson[3][3][0] *= ( ProbEta_c)/( p
742 Meson[3][3][1] *= (1.0-ProbEta_c)/(1.-p
743
744 Meson[4][4][0] *= ( ProbEta_b)/( p
745 Meson[4][4][1] *= (1.0-ProbEta_b)/(1.-p
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; BaryonWei
755 }
756 }
757 }
758
759 kfla =0; kflb =0;
760 G4int kflc(0), kfld(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 -
773
774 Baryon[i][j][k][0] = 1000 * k
775 BaryonWeight[i][j][k][0] = ( pspi
776 Baryon[i][j][k][1] = 1000 * k
777 BaryonWeight[i][j][k][1] = (1.-pspi
778 }
779 }
780 }
781
782 // Delta- ddd - only 1114
783 Baryon[0][0][0][0] = 1114; BaryonWeight
784 Baryon[0][0][0][1] = 0; BaryonWeight
785
786 // Delta++ uuu - only 2224
787 Baryon[1][1][1][0] = 2224; BaryonWeight
788 Baryon[1][1][1][1] = 0; BaryonWeight
789
790 // Omega- sss - only 3334
791 Baryon[2][2][2][0] = 3334; BaryonWeight
792 Baryon[2][2][2][1] = 0; BaryonWeight
793
794 // Omega_cc++ ccc - only 4444
795 Baryon[3][3][3][0] = 4444; BaryonWeight
796 Baryon[3][3][3][1] = 0; BaryonWeight
797
798 // Omega_bb- bbb - only 5554
799 Baryon[4][4][4][0] = 5554; BaryonWeight
800 Baryon[4][4][4][1] = 0; BaryonWeight
801
802 // Lambda/Sigma0 sud - 3122/3212
803 Baryon[0][1][2][0] = 3122; BaryonWeight
804 Baryon[0][2][1][0] = 3122; BaryonWeight
805 Baryon[1][0][2][0] = 3122; BaryonWeight
806 Baryon[1][2][0][0] = 3122; BaryonWeight
807 Baryon[2][0][1][0] = 3122; BaryonWeight
808 Baryon[2][1][0][0] = 3122; BaryonWeight
809
810 Baryon[0][1][2][2] = 3212; BaryonWeight
811 Baryon[0][2][1][2] = 3212; BaryonWeight
812 Baryon[1][0][2][2] = 3212; BaryonWeight
813 Baryon[1][2][0][2] = 3212; BaryonWeight
814 Baryon[2][0][1][2] = 3212; BaryonWeight
815 Baryon[2][1][0][2] = 3212; BaryonWeight
816
817 // Lambda_c+/Sigma_c+ cud - 4122/4212
818 Baryon[0][1][3][0] = 4122; BaryonWeight
819 Baryon[0][3][1][0] = 4122; BaryonWeight
820 Baryon[1][0][3][0] = 4122; BaryonWeight
821 Baryon[1][3][0][0] = 4122; BaryonWeight
822 Baryon[3][0][1][0] = 4122; BaryonWeight
823 Baryon[3][1][0][0] = 4122; BaryonWeight
824
825 Baryon[0][1][3][2] = 4212; BaryonWeight
826 Baryon[0][3][1][2] = 4212; BaryonWeight
827 Baryon[1][0][3][2] = 4212; BaryonWeight
828 Baryon[1][3][0][2] = 4212; BaryonWeight
829 Baryon[3][0][1][2] = 4212; BaryonWeight
830 Baryon[3][1][0][2] = 4212; BaryonWeight
831
832 // Xi_c+/Xi_c+' cus - 4232/4322
833 Baryon[1][2][3][0] = 4232; BaryonWeight
834 Baryon[1][3][2][0] = 4232; BaryonWeight
835 Baryon[2][1][3][0] = 4232; BaryonWeight
836 Baryon[2][3][1][0] = 4232; BaryonWeight
837 Baryon[3][1][2][0] = 4232; BaryonWeight
838 Baryon[3][2][1][0] = 4232; BaryonWeight
839
840 Baryon[1][2][3][2] = 4322; BaryonWeight
841 Baryon[1][3][2][2] = 4322; BaryonWeight
842 Baryon[2][1][3][2] = 4322; BaryonWeight
843 Baryon[2][3][1][2] = 4322; BaryonWeight
844 Baryon[3][1][2][2] = 4322; BaryonWeight
845 Baryon[3][2][1][2] = 4322; BaryonWeight
846
847 // Xi_c0/Xi_c0' cus - 4132/4312
848 Baryon[0][2][3][0] = 4132; BaryonWeight
849 Baryon[0][3][2][0] = 4132; BaryonWeight
850 Baryon[2][0][3][0] = 4132; BaryonWeight
851 Baryon[2][3][0][0] = 4132; BaryonWeight
852 Baryon[3][0][2][0] = 4132; BaryonWeight
853 Baryon[3][2][0][0] = 4132; BaryonWeight
854
855 Baryon[0][2][3][2] = 4312; BaryonWeight
856 Baryon[0][3][2][2] = 4312; BaryonWeight
857 Baryon[2][0][3][2] = 4312; BaryonWeight
858 Baryon[2][3][0][2] = 4312; BaryonWeight
859 Baryon[3][0][2][2] = 4312; BaryonWeight
860 Baryon[3][2][0][2] = 4312; BaryonWeight
861
862 // Lambda_b0/Sigma_b0 bud - 5122/5212
863 Baryon[0][1][4][0] = 5122; BaryonWeight
864 Baryon[0][4][1][0] = 5122; BaryonWeight
865 Baryon[1][0][4][0] = 5122; BaryonWeight
866 Baryon[1][4][0][0] = 5122; BaryonWeight
867 Baryon[4][0][1][0] = 5122; BaryonWeight
868 Baryon[4][1][0][0] = 5122; BaryonWeight
869
870 Baryon[0][1][4][2] = 5212; BaryonWeight
871 Baryon[0][4][1][2] = 5212; BaryonWeight
872 Baryon[1][0][4][2] = 5212; BaryonWeight
873 Baryon[1][4][0][2] = 5212; BaryonWeight
874 Baryon[4][0][1][2] = 5212; BaryonWeight
875 Baryon[4][1][0][2] = 5212; BaryonWeight
876
877 // Xi_b0/Xi_b0' bus - 5232/5322
878 Baryon[1][2][4][0] = 5232; BaryonWeight
879 Baryon[1][4][2][0] = 5232; BaryonWeight
880 Baryon[2][1][4][0] = 5232; BaryonWeight
881 Baryon[2][4][1][0] = 5232; BaryonWeight
882 Baryon[4][1][2][0] = 5232; BaryonWeight
883 Baryon[4][2][1][0] = 5232; BaryonWeight
884
885 Baryon[1][2][4][2] = 5322; BaryonWeight
886 Baryon[1][4][2][2] = 5322; BaryonWeight
887 Baryon[2][1][4][2] = 5322; BaryonWeight
888 Baryon[2][4][1][2] = 5322; BaryonWeight
889 Baryon[4][1][2][2] = 5322; BaryonWeight
890 Baryon[4][2][1][2] = 5322; BaryonWeight
891
892 // Xi_b-/Xi_b-' bus - 5132/5312
893 Baryon[0][2][4][0] = 5132; BaryonWeight
894 Baryon[0][4][2][0] = 5132; BaryonWeight
895 Baryon[2][0][4][0] = 5132; BaryonWeight
896 Baryon[2][4][0][0] = 5132; BaryonWeight
897 Baryon[4][0][2][0] = 5132; BaryonWeight
898 Baryon[4][2][0][0] = 5132; BaryonWeight
899
900 Baryon[0][2][4][2] = 5312; BaryonWeight
901 Baryon[0][4][2][2] = 5312; BaryonWeight
902 Baryon[2][0][4][2] = 5312; BaryonWeight
903 Baryon[2][4][0][2] = 5312; BaryonWeight
904 Baryon[4][0][2][2] = 5312; BaryonWeight
905 Baryon[4][2][0][2] = 5312; BaryonWeight
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 * Te
913 G4ParticleTable::GetPar
914 /*
915 G4cout<<i<<" "<<j<<" "<<k
916 if (TestHadron != nullptr
917 if ((TestHadron == nullpt
918 if ((TestHadron == nullpt
919 G4cout<<G4endl;
920 */
921 if ((TestHadron == nullpt
922 }
923 }
924 }
925 }
926
927 // --------- Probabilities of q-qbar pair
928 G4double ProbUUbar = 0.33;
929 Prob_QQbar[0]=ProbUUbar; // Probab
930 Prob_QQbar[1]=ProbUUbar; // Probab
931 Prob_QQbar[2]=1.0-2.*ProbUUbar; // Probab
932 Prob_QQbar[3]=0.0; // Probab
933 Prob_QQbar[4]=0.0; // Probab
934
935 for ( G4int i=0 ; i<350 ; i++ ) { // Must
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
946 void G4VLongitudinalStringDecay::SetMinimalStr
947 {
948 //MaxMass = -350.0*GeV;
949 G4double EstimatedMass=MaxMass;
950
951 G4ParticleDefinition* LeftParton = st
952 G4ParticleDefinition* RightParton = st
953 if( LeftParton->GetParticleSubType() =
954 if( LeftParton->GetPDGEncoding() * R
955 // Not allowed combination of the
956 throw G4HadronicException(__FILE__
957 "G4VLongitudinalStringDecay::Set
958 }
959 }
960 if( LeftParton->GetParticleSubType() !
961 if( LeftParton->GetPDGEncoding() * R
962 // Not allowed combination of the
963 throw G4HadronicException(__FILE__
964 "G4VLongitudinalStringDecay::Set
965 }
966 }
967
968 G4int Qleft =std::abs(string->GetLeftP
969 G4int Qright=std::abs(string->GetRight
970
971 if ((Qleft < 6) && (Qright < 6)) { /
972 EstimatedMass=minMassQQbarStr[Qleft-
973 MinimalStringMass=EstimatedMass;
974 SetMinimalStringMass2(EstimatedMass)
975 return;
976 }
977
978 if ((Qleft < 6) && (Qright > 1000)) {
979 G4int q1=Qright/1000;
980 G4int q2=(Qright/100)%10;
981 EstimatedMass=minMassQDiQStr[Qleft-1
982 MinimalStringMass=EstimatedMass;
983 SetMinimalStringMass2(EstimatedMass)
984 return;
985 }
986
987 if ((Qleft > 1000) && (Qright < 6)) {
988 G4int q1=Qleft/1000;
989 G4int q2=(Qleft/100)%10;
990 EstimatedMass=minMassQDiQStr[Qright-
991 MinimalStringMass=EstimatedMass;
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 Str
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
1011
1012 G4double EstimatedMass1 = minMassQDiQ
1013 G4double EstimatedMass2 = minMassQDiQ
1014 // Mass is negative if there is no co
1015
1016 if ( (EstimatedMass1 > 0.) && (Estima
1017 EstimatedMass = EstimatedMass1 + E
1018 if ( StringM > EstimatedMass ) {
1019 MinimalStringMass=EstimatedMass
1020 SetMinimalStringMass2(Estimated
1021 return;
1022 }
1023 }
1024
1025 if ( (EstimatedMass1 < 0.) && (Estima
1026 EstimatedMass = MaxMass;
1027 MinimalStringMass=EstimatedMass;
1028 SetMinimalStringMass2(EstimatedMas
1029 return;
1030 }
1031
1032 if ( (EstimatedMass1 > 0.) && (Estima
1033 EstimatedMass = EstimatedMass1;
1034 MinimalStringMass=EstimatedMass;
1035 SetMinimalStringMass2(EstimatedMas
1036 return;
1037 }
1038
1039 // if ( EstimatedMass >= StringM
1040 // ------------- Re-orangement ------
1041 EstimatedMass=std::min(minMassQQbarSt
1042 minMassQQbarSt
1043
1044 // In principle, re-arrangement and 2
1045 // More physics consideration is need
1046
1047 MinimalStringMass=EstimatedMass;
1048 SetMinimalStringMass2(EstimatedMass);
1049
1050 return;
1051 }
1052
1053 //-------------------------------------------
1054
1055 void G4VLongitudinalStringDecay::SetMinimalSt
1056 {
1057 MinimalStringMass2=aValue * aValue;
1058 }
1059
1060