| 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 // $Id: G4VLongitudinalStringDecay.cc,v 1.7 2007/01/24 10:00:11 gunter Exp $
>> 28 // GEANT4 tag $Name: geant4-08-02-patch-01 $
27 // 29 //
28 // ------------------------------------------- 30 // -----------------------------------------------------------------------------
29 // GEANT 4 class implementation file 31 // GEANT 4 class implementation file
30 // 32 //
31 // History: first implementation, Maxim K 33 // History: first implementation, Maxim Komogorov, 1-Jul-1998
32 // redesign Gunter Folger, Augu 34 // redesign Gunter Folger, August/September 2001
33 // ------------------------------------------- 35 // -----------------------------------------------------------------------------
34 #include "G4VLongitudinalStringDecay.hh" <<
35 #include "G4PhysicalConstants.hh" <<
36 #include "G4SystemOfUnits.hh" <<
37 #include "G4ios.hh" 36 #include "G4ios.hh"
38 #include "Randomize.hh" 37 #include "Randomize.hh"
>> 38 #include "G4VLongitudinalStringDecay.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 "G4ShortLivedTable.hh"
47 #include "G4PhaseSpaceDecayChannel.hh" 48 #include "G4PhaseSpaceDecayChannel.hh"
48 #include "G4VDecayChannel.hh" 49 #include "G4VDecayChannel.hh"
49 #include "G4DecayTable.hh" 50 #include "G4DecayTable.hh"
50 51
51 #include "G4DiQuarks.hh" 52 #include "G4DiQuarks.hh"
52 #include "G4Quarks.hh" 53 #include "G4Quarks.hh"
53 #include "G4Gluons.hh" 54 #include "G4Gluons.hh"
54 55
55 #include "G4Exp.hh" <<
56 #include "G4Log.hh" <<
57 <<
58 #include "G4HadronicException.hh" <<
59 <<
60 //------------------------debug switches 56 //------------------------debug switches
61 //#define debug_VStringDecay << 57 //#define DEBUG_LightFragmentationTest 1
62 //#define debug_heavyHadrons << 58
63 59
64 //******************************************** << 60 //********************************************************************************
65 // Constructors 61 // Constructors
66 62
67 G4VLongitudinalStringDecay::G4VLongitudinalStr << 63 G4VLongitudinalStringDecay::G4VLongitudinalStringDecay()
68 : G4HadronicInteraction(name), ProbCCbar(0.0 <<
69 { 64 {
70 MassCut = 210.0*MeV; // Mpi + Delta << 65 MassCut = 0.35*GeV;
>> 66 ClusterMass = 0.15*GeV;
71 67
72 StringLoopInterrupt = 1000; << 68 SmoothParam = 0.9;
73 ClusterLoopInterrupt = 500; << 69 StringLoopInterrupt = 1000;
>> 70 ClusterLoopInterrupt = 500;
74 71
75 // Changable Parameters below. << 72 // Changable Parameters below.
76 SigmaQT = 0.5 * GeV; <<
77 73
78 StrangeSuppress = 0.44; // =0.27/2.27 s << 74 SigmaQT = 0.5 * GeV;
79 DiquarkSuppress = 0.07; // Probability <<
80 DiquarkBreakProb = 0.1; // Probability <<
81 75
82 //... pspin_meson is probability to create << 76 StrangeSuppress = 0.44; // 27 % strange quarks produced, ie. u:d:s=1:1:0.27
83 pspin_meson.resize(3); << 77 DiquarkSuppress = 0.1;
84 pspin_meson[0] = 0.5; // u or d + anti-u o << 78 DiquarkBreakProb = 0.1;
85 pspin_meson[1] = 0.4; // one of the quark <<
86 pspin_meson[2] = 0.3; // both of the quark <<
87 79
88 //... pspin_barion is probability to create << 80 //... pspin_meson is probability to create vector meson
>> 81 pspin_meson = 0.5;
>> 82
>> 83 //... pspin_barion is probability to create 3/2 barion
89 pspin_barion = 0.5; 84 pspin_barion = 0.5;
90 85
91 //... vectorMesonMix[] is quark mixing para 86 //... vectorMesonMix[] is quark mixing parameters for vector mesons (Variable spin = 3)
92 vectorMesonMix.resize(6); 87 vectorMesonMix.resize(6);
93 vectorMesonMix[0] = 0.0; << 88 vectorMesonMix[0] = 0.5;
94 vectorMesonMix[1] = 0.5; << 89 vectorMesonMix[1] = 0.0;
95 vectorMesonMix[2] = 0.0; << 90 vectorMesonMix[2] = 0.5;
96 vectorMesonMix[3] = 0.5; << 91 vectorMesonMix[3] = 0.0;
97 vectorMesonMix[4] = 1.0; 92 vectorMesonMix[4] = 1.0;
98 vectorMesonMix[5] = 1.0; << 93 vectorMesonMix[5] = 1.0;
99 94
100 //... scalarMesonMix[] is quark mixing para 95 //... scalarMesonMix[] is quark mixing parameters for scalar mesons (Variable spin=1)
101 scalarMesonMix.resize(6); 96 scalarMesonMix.resize(6);
102 scalarMesonMix[0] = 0.5; << 97 scalarMesonMix[0] = 0.5;
103 scalarMesonMix[1] = 0.25; << 98 scalarMesonMix[1] = 0.25;
104 scalarMesonMix[2] = 0.5; << 99 scalarMesonMix[2] = 0.5;
105 scalarMesonMix[3] = 0.25; << 100 scalarMesonMix[3] = 0.25;
106 scalarMesonMix[4] = 1.0; << 101 scalarMesonMix[4] = 1.0;
107 scalarMesonMix[5] = 0.5; << 102 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 103
114 // Parameters may be changed until the firs << 104 // Parameters may be changed until the first fragmentation starts
115 PastInitPhase=false; 105 PastInitPhase=false;
116 hadronizer = new G4HadronBuilder( pspin_mes << 106 hadronizer = new G4HadronBuilder(pspin_meson,pspin_barion,
117 ProbEta_c << 107 scalarMesonMix,vectorMesonMix);
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 } 108 }
>> 109
126 110
127 G4VLongitudinalStringDecay::~G4VLongitudinalSt 111 G4VLongitudinalStringDecay::~G4VLongitudinalStringDecay()
128 { << 112 {
129 delete hadronizer; 113 delete hadronizer;
130 } << 114 }
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 115
211 G4ParticleDefinition* Hadron1 = nullpt << 116 //=============================================================================================-------------
212 G4ParticleDefinition* Hadron2 = nullptr; <<
213 117
214 if (!string->IsAFourQuarkString() ) << 118 // Operators
215 { <<
216 // spin 0 meson or spin 1/2 barion <<
217 119
218 Hadron1 = (hadronizer->*build)(stri << 120 //const & G4VLongitudinalStringDecay::operator=(const G4VLongitudinalStringDecay &)
219 #ifdef debug_VStringDecay << 121 // {
220 G4cout<<"VlongSD PossibleHadronMass"<<G4e << 122 // }
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 123
264 mass = (Hadron1)->GetPDGMass() + (Hadron2 << 124 //----------------------------------------------------------------------------------------------------------
265 } <<
266 125
267 #ifdef debug_VStringDecay << 126 int G4VLongitudinalStringDecay::operator==(const G4VLongitudinalStringDecay &) const
268 G4cout<<"VlongSD *Hadrons 1 and 2, pro << 127 {
269 #endif << 128 throw G4HadronicException(__FILE__, __LINE__, "G4VLongitudinalStringDecay::operator== forbidden");
270 << 129 return false;
271 if ( pdefs != 0 ) << 130 }
272 { // need to return hadrons as well.... <<
273 pdefs->first = Hadron1; <<
274 pdefs->second = Hadron2; <<
275 } <<
276 <<
277 return mass; <<
278 } <<
279 131
280 //-------------------------------------------- << 132 //----------------------------------------------------------------------------------------------------------
281 133
282 G4ParticleDefinition* G4VLongitudinalStringDec << 134 int G4VLongitudinalStringDecay::operator!=(const G4VLongitudinalStringDecay &) const
283 { << 135 {
284 /* << 136 throw G4HadronicException(__FILE__, __LINE__, "G4VLongitudinalStringDecay::operator!= forbidden");
285 G4cout<<Encoding<<" G4VLongitudinalStringDec << 137 return true;
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 } 138 }
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 139
312 //******************************************** << 140 //==========================================================================================================
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 141
328 G4ExcitedString *G4VLongitudinalStringDecay::C << 142 G4int G4VLongitudinalStringDecay::SampleQuarkFlavor(void)
329 { << 143 {
330 G4Parton *Left=new G4Parton(*in.GetLeftParto << 144 return (1 + (int)(G4UniformRand()/StrangeSuppress));
331 G4Parton *Right=new G4Parton(*in.GetRightPar << 145 }
332 return new G4ExcitedString(Left,Right,in.Get <<
333 } <<
334 146
335 //-------------------------------------------- << 147 //----------------------------------------------------------------------------------------------------------
336 148
337 G4ParticleDefinition * G4VLongitudinalStringDe << 149 G4VLongitudinalStringDecay::pDefPair G4VLongitudinalStringDecay::CreatePartonPair(G4int NeedParticle,G4bool AllowDiquarks)
338 <<
339 { 150 {
340 #ifdef debug_VStringDecay << 151 // NeedParticle = +1 for Particle, -1 for Antiparticle
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 152
362 G4VLongitudinalStringDecay::pDefPair G4VLongit <<
363 CreatePartonPair(G4int NeedParticle,G4bool All <<
364 { <<
365 // NeedParticle = +1 for Particle, -1 for <<
366 if ( AllowDiquarks && G4UniformRand() < Di 153 if ( AllowDiquarks && G4UniformRand() < DiquarkSuppress )
367 { 154 {
368 // Create a Diquark - AntiDiquark pair , 155 // Create a Diquark - AntiDiquark pair , first in pair is anti to IsParticle
369 #ifdef debug_VStringDecay << 156 G4int q1 = SampleQuarkFlavor();
370 G4cout<<"VlongSD Create a Diquark - Anti << 157 G4int q2 = SampleQuarkFlavor();
371 #endif << 158 G4int spin = (q1 != q2 && G4UniformRand() <= 0.5)? 1 : 3;
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 159 // convention: quark with higher PDG number is first
379 PDGcode = (std::max(q1,q2) * 1000 + std: << 160 G4int PDGcode = (std::max(q1,q2) * 1000 + std::min(q1,q2) * 100 + spin) * NeedParticle;
380 <<
381 return pDefPair (FindParticle(-PDGcode), 161 return pDefPair (FindParticle(-PDGcode),FindParticle(PDGcode));
>> 162
382 163
383 } else { 164 } else {
384 // Create a Quark - AntiQuark pair, firs 165 // Create a Quark - AntiQuark pair, first in pair IsParticle
385 #ifdef debug_VStringDecay <<
386 G4cout<<"VlongSD Create a Quark - AntiQu <<
387 #endif <<
388 G4int PDGcode=SampleQuarkFlavor()*NeedPa 166 G4int PDGcode=SampleQuarkFlavor()*NeedParticle;
389 return pDefPair (FindParticle(PDGcode),F 167 return pDefPair (FindParticle(PDGcode),FindParticle(-PDGcode));
390 } 168 }
391 } <<
392 169
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 } 170 }
415 171
416 //-------------------------------------------- << 172 //----------------------------------------------------------------------------------------------------------
417 173
418 G4ThreeVector G4VLongitudinalStringDecay::Samp << 174 // G4ThreeVector G4VLongitudinalStringDecay::SampleQuarkPt()
419 { << 175 // {
420 G4double Pt; << 176 // G4double width_param= 2.0 * GeV*GeV;
421 if ( ptMax < 0 ) { << 177 // G4double R = G4UniformRand();
422 // sample full gaussian << 178 // G4double Pt = std::sqrt(width_param*R/(1-R));
423 Pt = -G4Log(G4UniformRand()); << 179 // G4double phi = 2.*pi*G4UniformRand();
424 } else { << 180 // return G4ThreeVector(Pt * std::cos(phi),Pt * std::sin(phi),0);
425 // sample in limited range << 181 // }
426 G4double q = ptMax/SigmaQT; << 182 G4ThreeVector G4VLongitudinalStringDecay::SampleQuarkPt()
427 G4double ymin = (q > 20.) ? 0.0 : G4Exp( << 183 {
428 Pt = -G4Log(G4RandFlat::shoot(ymin, 1.)) << 184 G4double Pt = -std::log(G4UniformRand());
429 } <<
430 Pt = SigmaQT * std::sqrt(Pt); 185 Pt = SigmaQT * std::sqrt(Pt);
431 G4double phi = 2.*pi*G4UniformRand(); 186 G4double phi = 2.*pi*G4UniformRand();
432 return G4ThreeVector(Pt * std::cos(phi),Pt 187 return G4ThreeVector(Pt * std::cos(phi),Pt * std::sin(phi),0);
433 } << 188 }
434 189
435 //******************************************** << 190 //----------------------------------------------------------------------------------------------------------
436 191
437 void G4VLongitudinalStringDecay::CalculateHadr << 192 void G4VLongitudinalStringDecay::CalculateHadronTimePosition(G4double theInitialStringMass, G4KineticTrackVector* Hadrons)
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 { 193 {
>> 194 // `yo-yo` formation time
>> 195 const G4double kappa = 1.0 * GeV/fermi;
>> 196 for(size_t c1 = 0; c1 < Hadrons->size(); c1++)
>> 197 {
445 G4double SumPz = 0; 198 G4double SumPz = 0;
446 G4double SumE = 0; 199 G4double SumE = 0;
447 for (size_t c2 = 0; c2 < c1; c2++) << 200 for(size_t c2 = 0; c2 < c1; c2++)
448 { << 201 {
449 SumPz += Hadrons->operator[](c2)->Get 202 SumPz += Hadrons->operator[](c2)->Get4Momentum().pz();
450 SumE += Hadrons->operator[](c2)->Get 203 SumE += Hadrons->operator[](c2)->Get4Momentum().e();
451 } << 204 }
452 G4double HadronE = Hadrons->operator[]( 205 G4double HadronE = Hadrons->operator[](c1)->Get4Momentum().e();
453 G4double HadronPz = Hadrons->operator[]( 206 G4double HadronPz = Hadrons->operator[](c1)->Get4Momentum().pz();
454 Hadrons->operator[](c1)->SetFormationTim << 207 Hadrons->operator[](c1)->SetFormationTime((theInitialStringMass - 2.*SumPz + HadronE - HadronPz)/(2.*kappa));
455 (theInitialStringMass - 2.*SumPz + Had << 208 G4ThreeVector aPosition(0, 0, (theInitialStringMass - 2.*SumE - HadronE + HadronPz)/(2.*kappa));
456 G4ThreeVector aPosition( 0, 0, <<
457 (theInitialStringMass - 2.*SumE - Had <<
458 Hadrons->operator[](c1)->SetPosition(aPo 209 Hadrons->operator[](c1)->SetPosition(aPosition);
>> 210 }
459 } 211 }
460 } <<
461 212
462 //-------------------------------------------- << 213 //----------------------------------------------------------------------------------------------------------
463 214
464 void G4VLongitudinalStringDecay::SetSigmaTrans << 215 /*
465 { << 216 void G4VLongitudinalStringDecay::CalculateHadronTimePosition(G4double theInitialStringMass, G4KineticTrackVector* Hadrons)
466 if ( PastInitPhase ) { << 217 {
467 throw G4HadronicException(__FILE__, __LIN << 218 // 'constituent' formation time
468 "G4VLongitudinalStringDecay::SetSigmaTr << 219 const G4double kappa = 1.0 * GeV/fermi;
469 } else { << 220 for(G4int c1 = 0; c1 < Hadrons->length(); c1++)
470 SigmaQT = aValue; << 221 {
>> 222 G4double SumPz = 0;
>> 223 G4double SumE = 0;
>> 224 for(G4int c2 = 0; c2 <= c1; c2++)
>> 225 {
>> 226 SumPz += Hadrons->at(c2)->Get4Momentum().pz();
>> 227 SumE += Hadrons->at(c2)->Get4Momentum().e();
>> 228 }
>> 229 Hadrons->at(c1)->SetFormationTime((theInitialStringMass - 2.*SumPz)/(2.*kappa));
>> 230 G4ThreeVector aPosition(0, 0, (theInitialStringMass - 2.*SumE)/(2.*kappa));
>> 231 Hadrons->at(c1)->SetPosition(aPosition);
>> 232 }
>> 233 c1 = Hadrons->length()-1;
>> 234 Hadrons->at(c1)->SetFormationTime(Hadrons->at(c1-1)->GetFormationTime());
>> 235 Hadrons->at(c1)->SetPosition(Hadrons->at(c1-1)->GetPosition());
471 } 236 }
472 } << 237 */
473 238
474 //-------------------------------------------- 239 //----------------------------------------------------------------------------------------------------------
475 240
476 void G4VLongitudinalStringDecay::SetStrangenes << 241 G4ParticleDefinition *
477 { << 242 G4VLongitudinalStringDecay::QuarkSplitup(G4ParticleDefinition*
478 StrangeSuppress = aValue; << 243 decay, G4ParticleDefinition *&created)
>> 244 {
>> 245 G4int IsParticle=(decay->GetPDGEncoding()>0) ? -1 : +1; // if we have a quark, we need antiquark (or diquark)
>> 246 pDefPair QuarkPair = CreatePartonPair(IsParticle);
>> 247 created = QuarkPair.second;
>> 248 return hadronizer->Build(QuarkPair.first, decay);
>> 249
479 } 250 }
480 251
481 //-------------------------------------------- 252 //----------------------------------------------------------------------------------------------------------
482 253
483 void G4VLongitudinalStringDecay::SetDiquarkSup << 254 G4ParticleDefinition *G4VLongitudinalStringDecay::DiQuarkSplitup(
484 { << 255 G4ParticleDefinition* decay,
485 DiquarkSuppress = aValue; << 256 G4ParticleDefinition *&created)
>> 257 {
>> 258 //... can Diquark break or not?
>> 259 if (G4UniformRand() < DiquarkBreakProb ){
>> 260 //... Diquark break
>> 261
>> 262 G4int stableQuarkEncoding = decay->GetPDGEncoding()/1000;
>> 263 G4int decayQuarkEncoding = (decay->GetPDGEncoding()/100)%10;
>> 264 if (G4UniformRand() < 0.5)
>> 265 {
>> 266 G4int Swap = stableQuarkEncoding;
>> 267 stableQuarkEncoding = decayQuarkEncoding;
>> 268 decayQuarkEncoding = Swap;
>> 269 }
>> 270
>> 271 G4int IsParticle=(decayQuarkEncoding>0) ? -1 : +1;
>> 272 // if we have a quark, we need antiquark)
>> 273 pDefPair QuarkPair = CreatePartonPair(IsParticle,false); // no diquarks wanted
>> 274 //... Build new Diquark
>> 275 G4int QuarkEncoding=QuarkPair.second->GetPDGEncoding();
>> 276 G4int i10 = std::max(std::abs(QuarkEncoding), std::abs(stableQuarkEncoding));
>> 277 G4int i20 = std::min(std::abs(QuarkEncoding), std::abs(stableQuarkEncoding));
>> 278 G4int spin = (i10 != i20 && G4UniformRand() <= 0.5)? 1 : 3;
>> 279 G4int NewDecayEncoding = -1*IsParticle*(i10 * 1000 + i20 * 100 + spin);
>> 280 created = FindParticle(NewDecayEncoding);
>> 281 G4ParticleDefinition * decayQuark=FindParticle(decayQuarkEncoding);
>> 282 G4ParticleDefinition * had=hadronizer->Build(QuarkPair.first, decayQuark);
>> 283 return had;
>> 284 // return hadronizer->Build(QuarkPair.first, decayQuark);
>> 285
>> 286 } else {
>> 287 //... Diquark does not break
>> 288
>> 289 G4int IsParticle=(decay->GetPDGEncoding()>0) ? +1 : -1;
>> 290 // if we have a diquark, we need quark)
>> 291 pDefPair QuarkPair = CreatePartonPair(IsParticle,false); // no diquarks wanted
>> 292 created = QuarkPair.second;
>> 293
>> 294 G4ParticleDefinition * had=hadronizer->Build(QuarkPair.first, decay);
>> 295 return had;
>> 296 // return G4ParticleDefinition * had=hadronizer->Build(QuarkPair.first, decay);
>> 297 }
486 } 298 }
487 299
488 //-------------------------------------------- << 300 //-----------------------------------------------------------------------------------------
489 301
490 void G4VLongitudinalStringDecay::SetDiquarkBre << 302 G4KineticTrack * G4VLongitudinalStringDecay::Splitup(
>> 303 G4FragmentingString *string,
>> 304 G4FragmentingString *&newString)
491 { 305 {
492 if ( PastInitPhase ) { <<
493 throw G4HadronicException(__FILE__, __LINE <<
494 "G4VLongitudinalStringDecay::SetDiquarkB <<
495 } else { <<
496 DiquarkBreakProb = aValue; <<
497 } <<
498 } <<
499 306
500 //-------------------------------------------- << 307 //... random choice of string end to use for creating the hadron (decay)
>> 308 SideOfDecay = (G4UniformRand() < 0.5)? 1: -1;
>> 309 if (SideOfDecay < 0)
>> 310 {
>> 311 string->SetLeftPartonStable();
>> 312 } else
>> 313 {
>> 314 string->SetRightPartonStable();
>> 315 }
501 316
502 void G4VLongitudinalStringDecay::SetSpinThreeH << 317 G4ParticleDefinition *newStringEnd;
503 { << 318 G4ParticleDefinition * HadronDefinition;
504 if ( PastInitPhase ) { << 319 if (string->DecayIsQuark())
505 throw G4HadronicException(__FILE__, __LINE << 320 {
506 "G4VLongitudinalStringDecay::SetSpinThre << 321 HadronDefinition= QuarkSplitup(string->GetDecayParton(), newStringEnd);
507 } else { << 322 } else {
508 pspin_barion = aValue; << 323 HadronDefinition= DiQuarkSplitup(string->GetDecayParton(), newStringEnd);
509 delete hadronizer; << 324 }
510 hadronizer = new G4HadronBuilder( pspin_me << 325 // create new String from old, ie. keep Left and Right order, but replace decay
511 ProbEta_ << 326 G4LorentzVector* HadronMomentum=SplitEandP(HadronDefinition, string);
512 } << 327
>> 328 G4KineticTrack * Hadron =0;
>> 329 if ( HadronMomentum != 0 ) {
>> 330
>> 331 G4ThreeVector Pos;
>> 332 Hadron = new G4KineticTrack(HadronDefinition, 0,Pos, *HadronMomentum);
>> 333
>> 334 newString=new G4FragmentingString(*string,newStringEnd,
>> 335 HadronMomentum);
>> 336
>> 337 delete HadronMomentum;
>> 338 }
>> 339 return Hadron;
>> 340 }
>> 341
>> 342 //-----------------------------------------------------------------------------------------
>> 343
>> 344 G4LorentzVector * G4VLongitudinalStringDecay::SplitEandP(G4ParticleDefinition * pHadron,
>> 345 G4FragmentingString * string)
>> 346 {
>> 347 G4double HadronMass = pHadron->GetPDGMass();
>> 348
>> 349 // calculate and assign hadron transverse momentum component HadronPx andHadronPy
>> 350 G4ThreeVector thePt;
>> 351 thePt=SampleQuarkPt();
>> 352 G4ThreeVector HadronPt = thePt +string->DecayPt();
>> 353 HadronPt.setZ(0);
>> 354 //... sample z to define hadron longitudinal momentum and energy
>> 355 //... but first check the available phase space
>> 356 G4double DecayQuarkMass2 = sqr(string->GetDecayParton()->GetPDGMass());
>> 357 G4double HadronMass2T = sqr(HadronMass) + HadronPt.mag2();
>> 358 if (DecayQuarkMass2 + HadronMass2T >= SmoothParam*(string->Mass2()) )
>> 359 return 0; // have to start all over!
>> 360
>> 361 //... then compute allowed z region z_min <= z <= z_max
>> 362
>> 363 G4double zMin = HadronMass2T/(string->Mass2());
>> 364 G4double zMax = 1. - DecayQuarkMass2/(string->Mass2());
>> 365 if (zMin >= zMax) return 0; // have to start all over!
>> 366
>> 367 G4double z = GetLightConeZ(zMin, zMax,
>> 368 string->GetDecayParton()->GetPDGEncoding(), pHadron,
>> 369 HadronPt.x(), HadronPt.y());
>> 370
>> 371 //... now compute hadron longitudinal momentum and energy
>> 372 // longitudinal hadron momentum component HadronPz
>> 373
>> 374 HadronPt.setZ(0.5* string->GetDecayDirection() *
>> 375 (z * string->LightConeDecay() -
>> 376 HadronMass2T/(z * string->LightConeDecay())));
>> 377 G4double HadronE = 0.5* (z * string->LightConeDecay() +
>> 378 HadronMass2T/(z * string->LightConeDecay()));
>> 379
>> 380 G4LorentzVector * a4Momentum= new G4LorentzVector(HadronPt,HadronE);
>> 381
>> 382 return a4Momentum;
513 } 383 }
514 384
515 //-------------------------------------------- <<
516 385
517 void G4VLongitudinalStringDecay::SetScalarMeso << 386 //-----------------------------------------------------------------------------------------
518 { << 387
519 if ( PastInitPhase ) { << 388 G4bool G4VLongitudinalStringDecay::SplitLast(G4FragmentingString * string,
520 throw G4HadronicException(__FILE__, __LINE << 389 G4KineticTrackVector * LeftVector,
521 "G4VLongitudinalStringDecay::SetScalarMe << 390 G4KineticTrackVector * RightVector)
522 } else { << 391 {
523 if ( aVector.size() < 6 ) << 392 //... perform last cluster decay
524 throw G4HadronicException(__FILE__, __LI << 393 G4ThreeVector ClusterVel =string->Get4Momentum().boostVector();
525 "G4VLongitudinalStringDecay::SetScalar << 394 G4double ResidualMass =string->Mass();
526 scalarMesonMix[0] = aVector[0]; << 395 G4double ClusterMassCut = ClusterMass;
527 scalarMesonMix[1] = aVector[1]; << 396 G4int cClusterInterrupt = 0;
528 scalarMesonMix[2] = aVector[2]; << 397 G4ParticleDefinition * LeftHadron, * RightHadron;
529 scalarMesonMix[3] = aVector[3]; << 398 do
530 scalarMesonMix[4] = aVector[4]; << 399 {
531 scalarMesonMix[5] = aVector[5]; << 400 if (cClusterInterrupt++ >= ClusterLoopInterrupt)
532 delete hadronizer; << 401 {
533 hadronizer = new G4HadronBuilder( pspin_me << 402 return false;
534 ProbEta_ << 403 }
535 } << 404 G4ParticleDefinition * quark = NULL;
536 } << 405 string->SetLeftPartonStable(); // to query quark contents..
>> 406 if (string->DecayIsQuark() && string->StableIsQuark() )
>> 407 {
>> 408 //... there are quarks on cluster ends
>> 409 LeftHadron= QuarkSplitup(string->GetLeftParton(), quark);
>> 410 } else {
>> 411 //... there is a Diquark on cluster ends
>> 412 G4int IsParticle;
>> 413 if ( string->StableIsQuark() ) {
>> 414 IsParticle=(string->GetLeftParton()->GetPDGEncoding()>0) ? -1 : +1;
>> 415 } else {
>> 416 IsParticle=(string->GetLeftParton()->GetPDGEncoding()>0) ? +1 : -1;
>> 417 }
>> 418 pDefPair QuarkPair = CreatePartonPair(IsParticle,false); // no diquarks wanted
>> 419 quark = QuarkPair.second;
>> 420 LeftHadron=hadronizer->Build(QuarkPair.first, string->GetLeftParton());
>> 421 }
>> 422 RightHadron = hadronizer->Build(string->GetRightParton(), quark);
537 423
538 //-------------------------------------------- << 424 //... repeat procedure, if mass of cluster is too low to produce hadrons
>> 425 //... ClusterMassCut = 0.15*GeV model parameter
>> 426 if ( quark->GetParticleSubType()== "quark" ) {ClusterMassCut = 0.;}
>> 427 else {ClusterMassCut = ClusterMass;}
>> 428 }
>> 429 while (ResidualMass <= LeftHadron->GetPDGMass() + RightHadron->GetPDGMass() + ClusterMassCut);
>> 430
>> 431 //... compute hadron momenta and energies
>> 432 G4LorentzVector LeftMom, RightMom;
>> 433 G4ThreeVector Pos;
>> 434 Sample4Momentum(&LeftMom, LeftHadron->GetPDGMass(), &RightMom, RightHadron->GetPDGMass(), ResidualMass);
>> 435 LeftMom.boost(ClusterVel);
>> 436 RightMom.boost(ClusterVel);
>> 437 LeftVector->push_back(new G4KineticTrack(LeftHadron, 0, Pos, LeftMom));
>> 438 RightVector->push_back(new G4KineticTrack(RightHadron, 0, Pos, RightMom));
>> 439
>> 440 return true;
539 441
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 } 442 }
560 443
561 //-------------------------------------------- << 444 //----------------------------------------------------------------------------------------------------------
562 445
563 void G4VLongitudinalStringDecay::SetProbCCbar( << 446 G4KineticTrackVector* G4VLongitudinalStringDecay::FragmentString(const G4ExcitedString& theString)
564 { 447 {
565 ProbCCbar = aValue; << 448 // Can no longer modify Parameters for Fragmentation.
566 ProbCB = ProbCCbar + ProbBBbar; << 449 PastInitPhase=true;
567 } << 450
>> 451 // check if string has enough mass to fragment...
>> 452 G4KineticTrackVector * LeftVector=LightFragmentationTest(&theString);
>> 453 if ( LeftVector != 0 ) return LeftVector;
>> 454
>> 455 LeftVector = new G4KineticTrackVector;
>> 456 G4KineticTrackVector * RightVector=new G4KineticTrackVector;
568 457
569 //-------------------------------------------- << 458 // this should work but its only a semi deep copy. %GF G4ExcitedString theStringInCMS(theString);
>> 459 G4ExcitedString *theStringInCMS=CPExcited(theString);
>> 460 G4LorentzRotation toCms=theStringInCMS->TransformToAlignedCms();
>> 461
>> 462 G4bool success=false, inner_sucess=true;
>> 463 G4int attempt=0;
>> 464 while ( !success && attempt++ < StringLoopInterrupt )
>> 465 {
>> 466 G4FragmentingString *currentString=new G4FragmentingString(*theStringInCMS);
>> 467
>> 468 std::for_each(LeftVector->begin(), LeftVector->end(), DeleteKineticTrack());
>> 469 LeftVector->clear();
>> 470 std::for_each(RightVector->begin(), RightVector->end(), DeleteKineticTrack());
>> 471 RightVector->clear();
>> 472
>> 473 inner_sucess=true; // set false on failure..
>> 474 while (! StopFragmenting(currentString) )
>> 475 { // Split current string into hadron + new string
>> 476 G4FragmentingString *newString=0; // used as output from SplitUp...
>> 477 G4KineticTrack * Hadron=Splitup(currentString,newString);
>> 478 if ( Hadron != 0 && IsFragmentable(newString))
>> 479 {
>> 480 if ( currentString->GetDecayDirection() > 0 )
>> 481 LeftVector->push_back(Hadron);
>> 482 else
>> 483 RightVector->push_back(Hadron);
>> 484 delete currentString;
>> 485 currentString=newString;
>> 486 } else {
>> 487 // abandon ... start from the beginning
>> 488 if (newString) delete newString;
>> 489 if (Hadron) delete Hadron;
>> 490 inner_sucess=false;
>> 491 break;
>> 492 }
>> 493 }
>> 494 // Split current string into 2 final Hadrons
>> 495 if ( inner_sucess &&
>> 496 SplitLast(currentString,LeftVector, RightVector) )
>> 497 {
>> 498 success=true;
>> 499 }
>> 500 delete currentString;
>> 501 }
>> 502
>> 503 delete theStringInCMS;
>> 504
>> 505 if ( ! success )
>> 506 {
>> 507 std::for_each(LeftVector->begin(), LeftVector->end(), DeleteKineticTrack());
>> 508 LeftVector->clear();
>> 509 std::for_each(RightVector->begin(), RightVector->end(), DeleteKineticTrack());
>> 510 delete RightVector;
>> 511 return LeftVector;
>> 512 }
>> 513
>> 514 // Join Left- and RightVector into LeftVector in correct order.
>> 515 while(!RightVector->empty())
>> 516 {
>> 517 LeftVector->push_back(RightVector->back());
>> 518 RightVector->erase(RightVector->end()-1);
>> 519 }
>> 520 delete RightVector;
570 521
571 void G4VLongitudinalStringDecay::SetProbEta_c( << 522 CalculateHadronTimePosition(theString.Get4Momentum().mag(), LeftVector);
572 { <<
573 ProbEta_c = aValue; <<
574 } <<
575 523
576 //-------------------------------------------- << 524 G4LorentzRotation toObserverFrame(toCms.inverse());
577 525
578 void G4VLongitudinalStringDecay::SetProbBBbar( << 526 for(size_t C1 = 0; C1 < LeftVector->size(); C1++)
579 { << 527 {
580 ProbBBbar = aValue; << 528 G4KineticTrack* Hadron = LeftVector->operator[](C1);
581 ProbCB = ProbCCbar + ProbBBbar; << 529 G4LorentzVector Momentum = Hadron->Get4Momentum();
582 } << 530 Momentum = toObserverFrame*Momentum;
>> 531 Hadron->Set4Momentum(Momentum);
>> 532 G4LorentzVector Coordinate(Hadron->GetPosition(), Hadron->GetFormationTime());
>> 533 Momentum = toObserverFrame*Coordinate;
>> 534 Hadron->SetFormationTime(Momentum.e());
>> 535 G4ThreeVector aPosition(Momentum.vect());
>> 536 Hadron->SetPosition(theString.GetPosition()+aPosition);
>> 537 }
>> 538 return LeftVector;
>> 539
583 540
584 //-------------------------------------------- <<
585 541
586 void G4VLongitudinalStringDecay::SetProbEta_b( <<
587 { <<
588 ProbEta_b = aValue; <<
589 } 542 }
590 543
591 //-------------------------------------------- << 544 //----------------------------------------------------------------------------------------------------------
592 545
593 void G4VLongitudinalStringDecay::SetStringTens << 546 G4ExcitedString *G4VLongitudinalStringDecay::CPExcited(const G4ExcitedString & in)
594 { 547 {
595 Kappa = aValue * GeV/fermi; << 548 G4Parton *Left=new G4Parton(*in.GetLeftParton());
596 } << 549 G4Parton *Right=new G4Parton(*in.GetRightParton());
597 << 550 return new G4ExcitedString(Left,Right,in.GetDirection());
598 //-------------------------------------------- << 551 }
599 552
600 void G4VLongitudinalStringDecay::SetMinMasses( << 553 G4double G4VLongitudinalStringDecay::FragmentationMass(
>> 554 const G4FragmentingString *
>> 555 const string,
>> 556 Pcreate build,
>> 557 pDefPair * pdefs)
601 { 558 {
602 // ------ For estimation of a minimal stri << 559
603 Mass_of_light_quark =140.*MeV; << 560 G4double mass;
604 Mass_of_s_quark =500.*MeV; << 561 static G4bool NeedInit(true);
605 Mass_of_c_quark =1600.*MeV; << 562 static std::vector<double> nomix;
606 Mass_of_b_quark =4500.*MeV; << 563 static G4HadronBuilder * minMassHadronizer;
607 Mass_of_string_junction=720.*MeV; << 564 if ( NeedInit )
608 << 565 {
609 // ---------------- Determination of minim << 566 NeedInit = false;
610 G4ParticleDefinition * hadron1; G4int C << 567 nomix.resize(6);
611 G4ParticleDefinition * hadron2; G4int C << 568 for ( G4int i=0; i<6 ; i++ ) nomix[i]=0;
612 for (G4int i=1; i < 6; i++) { << 569 // minMassHadronizer=new G4HadronBuilder(pspin_meson,pspin_barion,nomix,nomix);
613 Code1 = 100*i + 10*1 + 1; << 570 minMassHadronizer=hadronizer;
614 hadron1 = FindParticle(Code1); << 571 }
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 572
651 if ((hadron2 != nullptr) && (h << 573 if ( build==0 ) build=&G4HadronBuilder::BuildLowSpin;
652 if (hadron2->GetPDGMass() > <<
653 }; <<
654 574
655 if ((hadron2 != nullptr) && (h << 575 G4ParticleDefinition *Hadron1, *Hadron2=0;
656 576
657 if ((hadron2 == nullptr) && (h << 577 if (!string->FourQuarkString() )
>> 578 {
>> 579 // spin 0 meson or spin 1/2 barion will be built
658 580
659 minMassQDiQStr[i-1][j-1][k-1] << 581 Hadron1 = (minMassHadronizer->*build)(string->GetLeftParton(),
660 } << 582 string->GetRightParton());
>> 583 mass= (Hadron1)->GetPDGMass();
>> 584 } else
>> 585 {
>> 586 //... string is qq--qqbar: Build two stable hadrons,
>> 587 //... with extra uubar or ddbar quark pair
>> 588 G4int iflc = (G4UniformRand() < 0.5)? 1 : 2;
>> 589 if (string->GetLeftParton()->GetPDGEncoding() < 0) iflc = -iflc;
>> 590
>> 591 //... theSpin = 4; spin 3/2 baryons will be built
>> 592 Hadron1 = (minMassHadronizer->*build)(string->GetLeftParton(),FindParticle(iflc));
>> 593 Hadron2 =(minMassHadronizer->*build)(string->GetRightParton(),FindParticle(-iflc));
>> 594 mass = (Hadron1)->GetPDGMass() + (Hadron2)->GetPDGMass();
661 } 595 }
662 } << 596
663 << 597 if ( pdefs != 0 )
664 // ------ An estimated minimal string mass << 598 { // need to return hadrons as well....
665 MinimalStringMass = 0.; << 599 pdefs->first = Hadron1;
666 MinimalStringMass2 = 0.; << 600 pdefs->second = Hadron2;
667 // q charges d u << 601 }
668 Qcharge[0] = -1; Qcharge[1] = 2; Qcharge[2 << 602
669 << 603 return mass;
670 // For treating of small string decays << 604 }
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 605
759 kfla =0; kflb =0; << 606 //----------------------------------------------------------------------------------------------------------
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 607
782 // Delta- ddd - only 1114 << 608 G4bool G4VLongitudinalStringDecay::IsFragmentable(const G4FragmentingString * const string)
783 Baryon[0][0][0][0] = 1114; BaryonWeight << 609 {
784 Baryon[0][0][0][1] = 0; BaryonWeight << 610 return sqr(FragmentationMass(string)+MassCut) <
785 << 611 string->Mass2();
786 // Delta++ uuu - only 2224 << 612 }
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 613
927 // --------- Probabilities of q-qbar pair << 614 //----------------------------------------------------------------------------------------------------------
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 615
941 NumberOf_FS = 0; << 616 G4bool G4VLongitudinalStringDecay::StopFragmenting(const G4FragmentingString * const string)
>> 617 {
>> 618 return
>> 619 sqr(FragmentationMass(string,&G4HadronBuilder::BuildHighSpin)+MassCut) >
>> 620 string->Get4Momentum().mag2();
942 } 621 }
943 622
944 // ------------------------------------------- << 623 //----------------------------------------------------------------------------------------------------------
945 624
946 void G4VLongitudinalStringDecay::SetMinimalStr << 625 G4KineticTrackVector* G4VLongitudinalStringDecay::LightFragmentationTest(const
>> 626 G4ExcitedString * const string)
947 { 627 {
948 //MaxMass = -350.0*GeV; << 628 // Check string decay threshold
949 G4double EstimatedMass=MaxMass; << 629
>> 630 G4KineticTrackVector * result=0; // return 0 when string exceeds the mass cut
>> 631
>> 632 pDefPair hadrons((G4ParticleDefinition *)0,(G4ParticleDefinition *)0);
>> 633 G4FragmentingString aString(*string);
>> 634 if ( sqr(FragmentationMass(&aString,0,&hadrons)+MassCut) < aString.Mass2()) {
>> 635 return 0;
>> 636 }
>> 637
>> 638 result=new G4KineticTrackVector;
>> 639
>> 640 if ( hadrons.second ==0 )
>> 641 {
>> 642 // Substitute string by light hadron, Note that Energy is not conserved here!
>> 643
>> 644 #ifdef DEBUG_LightFragmentationTest
>> 645 G4cout << "VlongSF Warning replacing string by single hadron "
>> 646 << hadrons.first->GetParticleName()
>> 647 << "string .. " << string->Get4Momentum() << " "
>> 648 << string->Get4Momentum().m() << G4endl;
>> 649 #endif
>> 650
>> 651 G4ThreeVector Mom3 = string->Get4Momentum().vect();
>> 652 G4LorentzVector Mom(Mom3,
>> 653 std::sqrt(Mom3.mag2() + sqr(hadrons.first->GetPDGMass())));
>> 654 result->push_back(new G4KineticTrack(hadrons.first, 0, string->GetPosition(), Mom));
>> 655 } else
>> 656 {
>> 657 //... string was qq--qqbar type: Build two stable hadrons,
>> 658
>> 659 #ifdef DEBUG_LightFragmentationTest
>> 660 G4cout << "VlongSF Warning replacing qq-qqbar string by TWO hadrons "
>> 661 << hadrons.first->GetParticleName() << " / "
>> 662 << hadrons.second->GetParticleName()
>> 663 << "string .. " << string->Get4Momentum() << " "
>> 664 << string->Get4Momentum().m() << G4endl;
>> 665 #endif
>> 666
>> 667 G4LorentzVector Mom1, Mom2;
>> 668 Sample4Momentum(&Mom1, hadrons.first->GetPDGMass(),
>> 669 &Mom2,hadrons.second->GetPDGMass(),
>> 670 string->Get4Momentum().mag());
>> 671 result->push_back(new G4KineticTrack(hadrons.first, 0, string->GetPosition(), Mom1));
>> 672 result->push_back(new G4KineticTrack(hadrons.second, 0, string->GetPosition(), Mom2));
>> 673 G4ThreeVector Velocity = string->Get4Momentum().boostVector();
>> 674 result->Boost(Velocity);
>> 675 }
950 676
951 G4ParticleDefinition* LeftParton = st << 677 return result;
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 678
968 G4int Qleft =std::abs(string->GetLeftP << 679 }
969 G4int Qright=std::abs(string->GetRight <<
970 680
971 if ((Qleft < 6) && (Qright < 6)) { / << 681 //----------------------------------------------------------------------------------------------------------
972 EstimatedMass=minMassQQbarStr[Qleft- <<
973 MinimalStringMass=EstimatedMass; <<
974 SetMinimalStringMass2(EstimatedMass) <<
975 return; <<
976 } <<
977 682
978 if ((Qleft < 6) && (Qright > 1000)) { << 683 G4ParticleDefinition* G4VLongitudinalStringDecay::FindParticle(G4int Encoding)
979 G4int q1=Qright/1000; << 684 {
980 G4int q2=(Qright/100)%10; << 685 G4ParticleDefinition* ptr = G4ParticleTable::GetParticleTable()->FindParticle(Encoding);
981 EstimatedMass=minMassQDiQStr[Qleft-1 << 686 if (ptr == NULL)
982 MinimalStringMass=EstimatedMass; << 687 {
983 SetMinimalStringMass2(EstimatedMass) << 688 G4cout << "Particle with encoding "<<Encoding<<" does not exist!!!"<<G4endl;
984 return; << 689 throw G4HadronicException(__FILE__, __LINE__, "Check your particle table");
985 } << 690 }
>> 691 return ptr;
>> 692 }
986 693
987 if ((Qleft > 1000) && (Qright < 6)) { << 694 //----------------------------------------------------------------------------------------------------------
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 695
996 // DiQuark - Anti DiQuark string ----- << 696 void G4VLongitudinalStringDecay::Sample4Momentum(G4LorentzVector* Mom, G4double Mass, G4LorentzVector* AntiMom, G4double AntiMass, G4double InitialMass)
>> 697 {
>> 698 G4double r_val = sqr(InitialMass*InitialMass - Mass*Mass - AntiMass*AntiMass) - sqr(2.*Mass*AntiMass);
>> 699 G4double Pabs = (r_val > 0.)? std::sqrt(r_val)/(2.*InitialMass) : 0;
997 700
998 G4double StringM=string->Get4Momentum().mag( << 701 //... sample unit vector
>> 702 G4double pz = 1. - 2.*G4UniformRand();
>> 703 G4double st = std::sqrt(1. - pz * pz)*Pabs;
>> 704 G4double phi = 2.*pi*G4UniformRand();
>> 705 G4double px = st*std::cos(phi);
>> 706 G4double py = st*std::sin(phi);
>> 707 pz *= Pabs;
>> 708
>> 709 Mom->setPx(px); Mom->setPy(py); Mom->setPz(pz);
>> 710 Mom->setE(std::sqrt(Pabs*Pabs + Mass*Mass));
999 711
1000 #ifdef debug_LUNDfragmentation << 712 AntiMom->setPx(-px); AntiMom->setPy(-py); AntiMom->setPz(-pz);
1001 // G4cout<<"MinStringMass// Input Str << 713 AntiMom->setE (std::sqrt(Pabs*Pabs + AntiMass*AntiMass));
1002 #endif << 714 }
1003 715
1004 G4int q1= Qleft/1000 ; << 716 //----------------------------------------------------------------------------------------------------------
1005 G4int q2=(Qleft/100)%10 ; <<
1006 717
1007 G4int q3= Qright/1000 ; << 718 void G4VLongitudinalStringDecay::SetSigmaTransverseMomentum(G4double aValue)
1008 G4int q4=(Qright/100)%10; << 719 {
>> 720 if ( PastInitPhase ) {
>> 721 throw G4HadronicException(__FILE__, __LINE__, "4VLongitudinalStringDecay::SetSigmaTransverseMomentum after FragmentString() not allowed");
>> 722 } else {
>> 723 SigmaQT = aValue;
>> 724 }
>> 725 }
1009 726
1010 // -------------- 2 baryon production << 727 //----------------------------------------------------------------------------------------------------------
1011 728
1012 G4double EstimatedMass1 = minMassQDiQ << 729 void G4VLongitudinalStringDecay::SetStrangenessSuppression(G4double aValue)
1013 G4double EstimatedMass2 = minMassQDiQ << 730 {
1014 // Mass is negative if there is no co << 731 if ( PastInitPhase ) {
>> 732 throw G4HadronicException(__FILE__, __LINE__, "4VLongitudinalStringDecay::SetStrangenessSuppression after FragmentString() not allowed");
>> 733 } else {
>> 734 StrangeSuppress = aValue;
>> 735 }
>> 736 }
1015 737
1016 if ( (EstimatedMass1 > 0.) && (Estima << 738 //----------------------------------------------------------------------------------------------------------
1017 EstimatedMass = EstimatedMass1 + E <<
1018 if ( StringM > EstimatedMass ) { <<
1019 MinimalStringMass=EstimatedMass <<
1020 SetMinimalStringMass2(Estimated <<
1021 return; <<
1022 } <<
1023 } <<
1024 739
1025 if ( (EstimatedMass1 < 0.) && (Estima << 740 void G4VLongitudinalStringDecay::SetDiquarkSuppression(G4double aValue)
1026 EstimatedMass = MaxMass; << 741 {
1027 MinimalStringMass=EstimatedMass; << 742 if ( PastInitPhase ) {
1028 SetMinimalStringMass2(EstimatedMas << 743 throw G4HadronicException(__FILE__, __LINE__, "4VLongitudinalStringDecay::SetDiquarkSuppression after FragmentString() not allowed");
1029 return; << 744 } else {
1030 } << 745 DiquarkSuppress = aValue;
>> 746 }
>> 747 }
1031 748
1032 if ( (EstimatedMass1 > 0.) && (Estima << 749 void G4VLongitudinalStringDecay::SetDiquarkBreakProbability(G4double aValue)
1033 EstimatedMass = EstimatedMass1; << 750 {
1034 MinimalStringMass=EstimatedMass; << 751 if ( PastInitPhase ) {
1035 SetMinimalStringMass2(EstimatedMas << 752 throw G4HadronicException(__FILE__, __LINE__, "4VLongitudinalStringDecay::SetDiquarkBreakProbability after FragmentString() not allowed");
1036 return; << 753 } else {
1037 } << 754 DiquarkBreakProb = aValue;
>> 755 }
>> 756 }
1038 757
1039 // if ( EstimatedMass >= StringM << 758 //----------------------------------------------------------------------------------------------------------
1040 // ------------- Re-orangement ------ <<
1041 EstimatedMass=std::min(minMassQQbarSt <<
1042 minMassQQbarSt <<
1043 759
1044 // In principle, re-arrangement and 2 << 760 void G4VLongitudinalStringDecay::SetVectorMesonProbability(G4double aValue)
1045 // More physics consideration is need << 761 {
>> 762 if ( PastInitPhase ) {
>> 763 throw G4HadronicException(__FILE__, __LINE__, "G4VLongitudinalStringDecay::SetVectorMesonProbability after FragmentString() not allowed");
>> 764 } else {
>> 765 pspin_meson = aValue;
>> 766 delete hadronizer;
>> 767 hadronizer = new G4HadronBuilder(pspin_meson,pspin_barion,
>> 768 scalarMesonMix,vectorMesonMix);
>> 769 }
>> 770 }
1046 771
1047 MinimalStringMass=EstimatedMass; << 772 //----------------------------------------------------------------------------------------------------------
1048 SetMinimalStringMass2(EstimatedMass); <<
1049 773
1050 return; << 774 void G4VLongitudinalStringDecay::SetSpinThreeHalfBarionProbability(G4double aValue)
>> 775 {
>> 776 if ( PastInitPhase ) {
>> 777 throw G4HadronicException(__FILE__, __LINE__, "G4VLongitudinalStringDecay::SetSpinThreeHalfBarionProbability after FragmentString() not allowed");
>> 778 } else {
>> 779 pspin_barion = aValue;
>> 780 delete hadronizer;
>> 781 hadronizer = new G4HadronBuilder(pspin_meson,pspin_barion,
>> 782 scalarMesonMix,vectorMesonMix);
>> 783 }
1051 } 784 }
1052 785
1053 //------------------------------------------- << 786 //----------------------------------------------------------------------------------------------------------
1054 787
1055 void G4VLongitudinalStringDecay::SetMinimalSt << 788 void G4VLongitudinalStringDecay::SetScalarMesonMixings(std::vector<G4double> aVector)
1056 { 789 {
1057 MinimalStringMass2=aValue * aValue; << 790 if ( PastInitPhase ) {
>> 791 throw G4HadronicException(__FILE__, __LINE__, "G4VLongitudinalStringDecay::SetScalarMesonMixings after FragmentString() not allowed");
>> 792 } else {
>> 793 if ( aVector.size() < 6 )
>> 794 throw G4HadronicException(__FILE__, __LINE__, "G4VLongitudinalStringDecay::SetScalarMesonMixings( argument Vector too small");
>> 795 scalarMesonMix[0] = aVector[0];
>> 796 scalarMesonMix[1] = aVector[1];
>> 797 scalarMesonMix[2] = aVector[2];
>> 798 scalarMesonMix[3] = aVector[3];
>> 799 scalarMesonMix[4] = aVector[4];
>> 800 scalarMesonMix[5] = aVector[5];
>> 801 delete hadronizer;
>> 802 hadronizer = new G4HadronBuilder(pspin_meson,pspin_barion,
>> 803 scalarMesonMix,vectorMesonMix);
>> 804 }
1058 } 805 }
1059 806
>> 807 //----------------------------------------------------------------------------------------------------------
>> 808
>> 809 void G4VLongitudinalStringDecay::SetVectorMesonMixings(std::vector<G4double> aVector)
>> 810 {
>> 811 if ( PastInitPhase ) {
>> 812 throw G4HadronicException(__FILE__, __LINE__, "G4VLongitudinalStringDecay::SetVectorMesonMixings after FragmentString() not allowed");
>> 813 } else {
>> 814 if ( aVector.size() < 6 )
>> 815 throw G4HadronicException(__FILE__, __LINE__, "G4VLongitudinalStringDecay::SetVectorMesonMixings( argument Vector too small");
>> 816 vectorMesonMix[0] = aVector[0];
>> 817 vectorMesonMix[1] = aVector[1];
>> 818 vectorMesonMix[2] = aVector[2];
>> 819 vectorMesonMix[3] = aVector[3];
>> 820 vectorMesonMix[4] = aVector[4];
>> 821 vectorMesonMix[5] = aVector[5];
>> 822 delete hadronizer;
>> 823 hadronizer = new G4HadronBuilder(pspin_meson,pspin_barion,
>> 824 scalarMesonMix,vectorMesonMix);
>> 825
>> 826 }
>> 827 }
>> 828 //*******************************************************************************************************
1060 829