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