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