| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // 26 //
27 // 27 //
28 // ------------------------------------------- 28 // -----------------------------------------------------------------------------
29 // GEANT 4 class implementation file 29 // GEANT 4 class implementation file
30 // 30 //
31 // History: first implementation, Maxim K 31 // History: first implementation, Maxim Komogorov, 10-Jul-1998
32 // ------------------------------------------- 32 // -----------------------------------------------------------------------------
33 #include "G4LundStringFragmentation.hh" 33 #include "G4LundStringFragmentation.hh"
34 #include "G4PhysicalConstants.hh" 34 #include "G4PhysicalConstants.hh"
35 #include "G4SystemOfUnits.hh" 35 #include "G4SystemOfUnits.hh"
36 #include "Randomize.hh" 36 #include "Randomize.hh"
37 #include "G4FragmentingString.hh" 37 #include "G4FragmentingString.hh"
38 #include "G4DiQuarks.hh" 38 #include "G4DiQuarks.hh"
39 #include "G4Quarks.hh" 39 #include "G4Quarks.hh"
40 #include "G4HadronicParameters.hh" 40 #include "G4HadronicParameters.hh"
41 #include "G4Exp.hh" 41 #include "G4Exp.hh"
42 #include "G4Pow.hh" 42 #include "G4Pow.hh"
43 43
44 //#define debug_LUNDfragmentation 44 //#define debug_LUNDfragmentation
45 45
46 // Class G4LundStringFragmentation 46 // Class G4LundStringFragmentation
47 //******************************************** 47 //*************************************************************************************
48 48
49 G4LundStringFragmentation::G4LundStringFragmen 49 G4LundStringFragmentation::G4LundStringFragmentation()
50 : G4VLongitudinalStringDecay("LundStringFrag 50 : G4VLongitudinalStringDecay("LundStringFragmentation")
51 { 51 {
52 SetMassCut(210.*MeV); // Mpi + Delta 52 SetMassCut(210.*MeV); // Mpi + Delta
53 // For ProduceOneH 53 // For ProduceOneHadron it is required
54 // that no one pi- 54 // that no one pi-meson can be produced.
55 SigmaQT = 0.435 * GeV; 55 SigmaQT = 0.435 * GeV;
56 Tmt = 190.0 * MeV; << 56 Tmt = 190.0 * MeV;
57 <<
58 SetStringTensionParameter(1.*GeV/fermi); 57 SetStringTensionParameter(1.*GeV/fermi);
59 SetDiquarkBreakProbability(0.3); << 58 SetDiquarkBreakProbability(0.5);
60 <<
61 SetStrangenessSuppression((1.0 - 0.12)/2.0 59 SetStrangenessSuppression((1.0 - 0.12)/2.0);
62 SetDiquarkSuppression(0.07); << 60 SetDiquarkSuppression(0.15);
63 61
64 // Check if charmed and bottom hadrons are 62 // Check if charmed and bottom hadrons are enabled: if this is the case, then
65 // set the non-zero probabilities for c-cb 63 // set the non-zero probabilities for c-cbar and b-bbar creation from the vacuum,
66 // else set them to 0.0. If these probabil 64 // else set them to 0.0. If these probabilities are/aren't zero then charmed or bottom
67 // hadrons can't/can be created during the 65 // hadrons can't/can be created during the string fragmentation of ordinary
68 // (i.e. not heavy) projectile hadron nucl 66 // (i.e. not heavy) projectile hadron nuclear reactions.
69 if ( G4HadronicParameters::Instance()->Ena 67 if ( G4HadronicParameters::Instance()->EnableBCParticles() ) {
70 SetProbCCbar(0.0002); // According to O << 68 SetProbCCbar(0.005); // According to O.I. Piskunova Yad. Fiz. 56 (1993) 1094
71 SetProbBBbar(5.0e-5); // According to O 69 SetProbBBbar(5.0e-5); // According to O.I. Piskunova Yad. Fiz. 56 (1993) 1094
72 } else { 70 } else {
73 SetProbCCbar(0.0); 71 SetProbCCbar(0.0);
74 SetProbBBbar(0.0); 72 SetProbBBbar(0.0);
75 } 73 }
76 74
77 SetMinMasses(); // For treating of small 75 SetMinMasses(); // For treating of small string decays
78 } 76 }
79 77
80 //-------------------------------------------- 78 //--------------------------------------------------------------------------------------
81 79
82 G4KineticTrackVector* G4LundStringFragmentatio 80 G4KineticTrackVector* G4LundStringFragmentation::FragmentString(const G4ExcitedString& theString)
83 { 81 {
84 // Can no longer modify Parameters for Fragm 82 // Can no longer modify Parameters for Fragmentation.
85 83
86 PastInitPhase=true; 84 PastInitPhase=true;
87 85
88 G4FragmentingString aString(theString); 86 G4FragmentingString aString(theString);
89 SetMinimalStringMass(&aString); 87 SetMinimalStringMass(&aString);
90 88
91 #ifdef debug_LUNDfragmentation 89 #ifdef debug_LUNDfragmentation
92 G4cout<<G4endl<<"LUND StringFragmentat 90 G4cout<<G4endl<<"LUND StringFragmentation ------------------------------------"<<G4endl;
93 G4cout<<G4endl<<"LUND StringFragm: Str 91 G4cout<<G4endl<<"LUND StringFragm: String Mass "
94 <<theString.Get4Momentum 92 <<theString.Get4Momentum().mag()<<G4endl
95 <<"4Mom "<<theString.Get 93 <<"4Mom "<<theString.Get4Momentum()<<G4endl
96 <<"--------------------- 94 <<"------------------------------------"<<G4endl;
97 G4cout<<"String ends Direct "<<theStri 95 G4cout<<"String ends Direct "<<theString.GetLeftParton()->GetPDGcode()<<" "
98 <<theStri 96 <<theString.GetRightParton()->GetPDGcode()<<" "
99 <<theStri 97 <<theString.GetDirection()<< G4endl;
100 G4cout<<"Left mom "<<theString.GetLef 98 G4cout<<"Left mom "<<theString.GetLeftParton()->Get4Momentum()<<G4endl;
101 G4cout<<"Right mom "<<theString.GetRig 99 G4cout<<"Right mom "<<theString.GetRightParton()->Get4Momentum()<<G4endl<<G4endl;
102 G4cout<<"Check for Fragmentation "<<G4 100 G4cout<<"Check for Fragmentation "<<G4endl;
103 #endif 101 #endif
104 102
105 G4KineticTrackVector * LeftVector(0); 103 G4KineticTrackVector * LeftVector(0);
106 104
107 if (!aString.IsAFourQuarkString() && !IsItFr 105 if (!aString.IsAFourQuarkString() && !IsItFragmentable(&aString))
108 { 106 {
109 #ifdef debug_LUNDfragmentation 107 #ifdef debug_LUNDfragmentation
110 G4cout<<"Non fragmentable - th 108 G4cout<<"Non fragmentable - the string is converted to one hadron "<<G4endl;
111 #endif 109 #endif
112 // SetMassCut(210.*MeV); // F 110 // SetMassCut(210.*MeV); // For ProduceOneHadron it is required
113 // t 111 // that no one pi-meson can be produced.
114 112
115 G4double Mcut = GetMassCut(); 113 G4double Mcut = GetMassCut();
116 SetMassCut(10000.*MeV); 114 SetMassCut(10000.*MeV);
117 LeftVector=ProduceOneHadron(&theString); 115 LeftVector=ProduceOneHadron(&theString);
118 SetMassCut(Mcut); 116 SetMassCut(Mcut);
119 117
120 if ( LeftVector ) 118 if ( LeftVector )
121 { 119 {
122 if ( LeftVector->size() > 0) 120 if ( LeftVector->size() > 0)
123 { 121 {
124 LeftVector->operator[](0)->SetForm 122 LeftVector->operator[](0)->SetFormationTime(theString.GetTimeOfCreation());
125 LeftVector->operator[](0)->SetPosi 123 LeftVector->operator[](0)->SetPosition(theString.GetPosition());
126 } 124 }
127 if (LeftVector->size() > 1) 125 if (LeftVector->size() > 1)
128 { 126 {
129 // 2 hadrons created from qq-qqbar 127 // 2 hadrons created from qq-qqbar are stored
130 LeftVector->operator[](1)->SetFormationT 128 LeftVector->operator[](1)->SetFormationTime(theString.GetTimeOfCreation());
131 LeftVector->operator[](1)->SetPosition(t 129 LeftVector->operator[](1)->SetPosition(theString.GetPosition());
132 } 130 }
133 } 131 }
134 return LeftVector; 132 return LeftVector;
135 } 133 }
136 134
137 #ifdef debug_LUNDfragmentation 135 #ifdef debug_LUNDfragmentation
138 G4cout<<"The string will be fragmented 136 G4cout<<"The string will be fragmented. "<<G4endl;
139 #endif 137 #endif
140 138
141 // The string can fragment. At least two par 139 // The string can fragment. At least two particles can be produced.
142 LeftVector =new G4KineticTrackVec 140 LeftVector =new G4KineticTrackVector;
143 G4KineticTrackVector * RightVector=new G4Kin 141 G4KineticTrackVector * RightVector=new G4KineticTrackVector;
144 142
145 G4bool success = Loop_toFragmentString 143 G4bool success = Loop_toFragmentString(theString, LeftVector, RightVector);
146 144
147 if ( ! success ) 145 if ( ! success )
148 { 146 {
149 std::for_each(LeftVector->begin(), LeftVec 147 std::for_each(LeftVector->begin(), LeftVector->end(), DeleteKineticTrack());
150 LeftVector->clear(); 148 LeftVector->clear();
151 std::for_each(RightVector->begin(), RightV 149 std::for_each(RightVector->begin(), RightVector->end(), DeleteKineticTrack());
152 delete RightVector; 150 delete RightVector;
153 return LeftVector; 151 return LeftVector;
154 } 152 }
155 153
156 // Join Left- and RightVector into LeftVecto 154 // Join Left- and RightVector into LeftVector in correct order.
157 while (!RightVector->empty()) 155 while (!RightVector->empty())
158 { 156 {
159 LeftVector->push_back(RightVector->back()) 157 LeftVector->push_back(RightVector->back());
160 RightVector->erase(RightVector->end()-1); 158 RightVector->erase(RightVector->end()-1);
161 } 159 }
162 delete RightVector; 160 delete RightVector;
163 161
164 return LeftVector; 162 return LeftVector;
165 } 163 }
166 164
167 //-------------------------------------------- 165 //----------------------------------------------------------------------------------
168 166
169 G4bool G4LundStringFragmentation::IsItFragment 167 G4bool G4LundStringFragmentation::IsItFragmentable(const G4FragmentingString * const string)
170 { 168 {
171 SetMinimalStringMass(string); 169 SetMinimalStringMass(string);
172 //G4cout<<"MinM StrM "<<MinimalStringM 170 //G4cout<<"MinM StrM "<<MinimalStringMass<<" "<< string->Get4Momentum().mag()<<G4endl;
173 171
174 return std::abs(MinimalStringMass) < string- 172 return std::abs(MinimalStringMass) < string->Get4Momentum().mag();
175 173
176 //MinimalStringMass is negative and la 174 //MinimalStringMass is negative and large for a string with unknown particles in a final 2-particle decay.
177 } 175 }
178 176
179 //-------------------------------------------- 177 //----------------------------------------------------------------------------------------
180 178
181 G4bool G4LundStringFragmentation::Loop_toFragm 179 G4bool G4LundStringFragmentation::Loop_toFragmentString( const G4ExcitedString &theString,
182 180 G4KineticTrackVector * & LeftVector,
183 181 G4KineticTrackVector * & RightVector )
184 { 182 {
185 #ifdef debug_LUNDfragmentation 183 #ifdef debug_LUNDfragmentation
186 G4cout<<"Loop_toFrag "<<theString.GetL 184 G4cout<<"Loop_toFrag "<<theString.GetLeftParton()->GetPDGcode()<<" "
187 <<theString.GetL 185 <<theString.GetLeftParton()->Get4Momentum()<<G4endl
188 <<" "<<theString.GetR 186 <<" "<<theString.GetRightParton()->GetPDGcode()<<" "
189 <<theString.GetR 187 <<theString.GetRightParton()->Get4Momentum()<<G4endl
190 <<"Direction "<<theString.GetD 188 <<"Direction "<<theString.GetDirection()<< G4endl;
191 #endif 189 #endif
192 190
193 G4LorentzRotation toCmsI, toObserverFr 191 G4LorentzRotation toCmsI, toObserverFrameI;
194 192
195 G4bool final_success=false; 193 G4bool final_success=false;
196 G4bool inner_success=true; 194 G4bool inner_success=true;
197 195
198 G4int attempt=0; 196 G4int attempt=0;
199 197
200 while ( ! final_success && attempt++ < Strin 198 while ( ! final_success && attempt++ < StringLoopInterrupt )
201 { // If the string fragmentation does 199 { // If the string fragmentation does not be happend,
202 // repeat the fragmentation. 200 // repeat the fragmentation.
203 201
204 G4FragmentingString *currentSt 202 G4FragmentingString *currentString = new G4FragmentingString(theString);
205 toCmsI = currentString->Transf 203 toCmsI = currentString->TransformToAlignedCms();
206 toObserverFrameI = toCmsI.inve 204 toObserverFrameI = toCmsI.inverse();
207 205
208 G4LorentzRotation toCms, toObserverFrame; 206 G4LorentzRotation toCms, toObserverFrame;
209 207
210 //G4cout<<"Main loop start whilecounter "< 208 //G4cout<<"Main loop start whilecounter "<<attempt<<G4endl;
211 209
212 // Cleaning up the previously produced had 210 // Cleaning up the previously produced hadrons
213 std::for_each(LeftVector->begin(), LeftVec 211 std::for_each(LeftVector->begin(), LeftVector->end(), DeleteKineticTrack());
214 LeftVector->clear(); 212 LeftVector->clear();
215 std::for_each(RightVector->begin(), RightV 213 std::for_each(RightVector->begin(), RightVector->end(), DeleteKineticTrack());
216 RightVector->clear(); 214 RightVector->clear();
217 215
218 // Main fragmentation loop until the strin 216 // Main fragmentation loop until the string will not be able to fragment
219 inner_success=true; // set false on failu 217 inner_success=true; // set false on failure.
220 const G4int maxNumberOfLoops = 218 const G4int maxNumberOfLoops = 1000;
221 G4int loopCounter = -1; 219 G4int loopCounter = -1;
222 220
223 while ( (! StopFragmenting(currentString)) 221 while ( (! StopFragmenting(currentString)) && ++loopCounter < maxNumberOfLoops )
224 { // Split current string into hadro 222 { // Split current string into hadron + new string
225 #ifdef debug_LUNDfragm 223 #ifdef debug_LUNDfragmentation
226 G4cout<<"The string wi 224 G4cout<<"The string will fragment. "<<G4endl;;
227 //G4cout<<"1 "<<curren 225 //G4cout<<"1 "<<currentString->GetDecayDirection()<<G4endl;
228 #endif 226 #endif
229 G4FragmentingString *newString=0; // us 227 G4FragmentingString *newString=0; // used as output from SplitUp.
230 228
231 toCms=currentString->TransformToAlignedC 229 toCms=currentString->TransformToAlignedCms();
232 toObserverFrame= toCms 230 toObserverFrame= toCms.inverse();
233 231
234 #ifdef debug_LUNDfragm 232 #ifdef debug_LUNDfragmentation
235 //G4cout<<"CMS Left m 233 //G4cout<<"CMS Left mom "<<currentString->GetPleft()<<G4endl;
236 //G4cout<<"CMS Right m 234 //G4cout<<"CMS Right mom "<<currentString->GetPright()<<G4endl;
237 //G4cout<<"CMS String 235 //G4cout<<"CMS String M "<<currentString->GetPstring()<<G4endl;
238 #endif 236 #endif
239 237
240 G4KineticTrack * Hadron=Splitup(currentS 238 G4KineticTrack * Hadron=Splitup(currentString,newString);
241 239
242 if ( Hadron != 0 ) // Store the hadron 240 if ( Hadron != 0 ) // Store the hadron
243 { 241 {
244 #ifdef debug_L 242 #ifdef debug_LUNDfragmentation
245 G4cout<<"Hadro 243 G4cout<<"Hadron prod at fragm. "<<Hadron->GetDefinition()->GetParticleName()<<G4endl;
246 //G4cout<<"2 " 244 //G4cout<<"2 "<<currentString->GetDecayDirection()<<G4endl;
247 #endif 245 #endif
248 246
249 Hadron->Set4Momentum(toObserverFrame*H 247 Hadron->Set4Momentum(toObserverFrame*Hadron->Get4Momentum());
250 248
251 G4double TimeOftheStringCreation=theSt 249 G4double TimeOftheStringCreation=theString.GetTimeOfCreation();
252 G4ThreeVector PositionOftheStringCreat 250 G4ThreeVector PositionOftheStringCreation(theString.GetPosition());
253 251
254 G4LorentzVector Coordinate(Hadron->Get 252 G4LorentzVector Coordinate(Hadron->GetPosition(), Hadron->GetFormationTime());
255 G4LorentzVector Momentum = toObserverF 253 G4LorentzVector Momentum = toObserverFrame*Coordinate;
256 Hadron->SetFormationTime(TimeOftheStri 254 Hadron->SetFormationTime(TimeOftheStringCreation + Momentum.e() - fermi/c_light);
257 G4ThreeVector aPosition(Momentum.vect( 255 G4ThreeVector aPosition(Momentum.vect());
258 Hadron->SetPosition(PositionOftheStrin 256 Hadron->SetPosition(PositionOftheStringCreation+aPosition);
259 257
260 // Open to pro 258 // Open to protect hadron production at fragmentation
261 if ( currentString->GetDecayDirection( 259 if ( currentString->GetDecayDirection() > 0 )
262 { 260 {
263 LeftVector->push_back(Hadron); 261 LeftVector->push_back(Hadron);
264 } else 262 } else
265 { 263 {
266 RightVector->push_back(Hadron); 264 RightVector->push_back(Hadron);
267 } 265 }
268 delete currentString; 266 delete currentString;
269 currentString=newString; 267 currentString=newString;
270 } else { 268 } else {
271 if ( newString ) del 269 if ( newString ) delete newString;
272 } 270 }
273 271
274 currentString->Lorentz 272 currentString->LorentzRotate(toObserverFrame);
275 }; 273 };
276 274
277 if ( loopCounter >= maxNumberO 275 if ( loopCounter >= maxNumberOfLoops ) {
278 inner_success=false; 276 inner_success=false;
279 } 277 }
280 278
281 // Split remaining string into 2 final had 279 // Split remaining string into 2 final hadrons.
282 #ifdef debug_LUNDfragmentation 280 #ifdef debug_LUNDfragmentation
283 if (inner_success) G4cout<<"Sp 281 if (inner_success) G4cout<<"Split remaining string into 2 final hadrons."<<G4endl;
284 #endif 282 #endif
285 283
286 if ( inner_success && SplitLast(currentStr 284 if ( inner_success && SplitLast(currentString, LeftVector, RightVector) ) // Close to protect Last Str. Decay
287 { 285 {
288 final_success = true; 286 final_success = true;
289 } 287 }
290 288
291 delete currentString; 289 delete currentString;
292 } // End of the loop where we try to fragme 290 } // End of the loop where we try to fragment the string.
293 291
294 G4int sign = +1; 292 G4int sign = +1;
295 if ( theString.GetDirection() < 0 ) si 293 if ( theString.GetDirection() < 0 ) sign = -1;
296 for ( unsigned int hadronI = 0; hadron 294 for ( unsigned int hadronI = 0; hadronI < LeftVector->size(); ++hadronI ) {
297 G4LorentzVector Tmp = LeftVector->o 295 G4LorentzVector Tmp = LeftVector->operator[](hadronI)->Get4Momentum();
298 Tmp.setZ(sign*Tmp.getZ()); 296 Tmp.setZ(sign*Tmp.getZ());
299 Tmp *= toObserverFrameI; 297 Tmp *= toObserverFrameI;
300 LeftVector->operator[](hadronI)->Se 298 LeftVector->operator[](hadronI)->Set4Momentum(Tmp);
301 } 299 }
302 for ( unsigned int hadronI = 0; hadron 300 for ( unsigned int hadronI = 0; hadronI < RightVector->size(); ++hadronI ) {
303 G4LorentzVector Tmp = RightVector-> 301 G4LorentzVector Tmp = RightVector->operator[](hadronI)->Get4Momentum();
304 Tmp.setZ(sign*Tmp.getZ()); 302 Tmp.setZ(sign*Tmp.getZ());
305 Tmp *= toObserverFrameI; 303 Tmp *= toObserverFrameI;
306 RightVector->operator[](hadronI)->S 304 RightVector->operator[](hadronI)->Set4Momentum(Tmp);
307 } 305 }
308 306
309 return final_success; 307 return final_success;
310 } 308 }
311 309
312 //-------------------------------------------- 310 //----------------------------------------------------------------------------------------
313 311
314 G4bool G4LundStringFragmentation::StopFragment 312 G4bool G4LundStringFragmentation::StopFragmenting(const G4FragmentingString * const string)
315 { 313 {
316 SetMinimalStringMass(string); 314 SetMinimalStringMass(string);
317 315
318 if ( MinimalStringMass < 0.) return true; 316 if ( MinimalStringMass < 0.) return true;
319 317
320 if (string->IsAFourQuarkString()) 318 if (string->IsAFourQuarkString())
321 { 319 {
322 return G4UniformRand() < G4Exp(-0.0005*(st 320 return G4UniformRand() < G4Exp(-0.0005*(string->Mass() - MinimalStringMass));
323 } else { 321 } else {
324 322
325 if (MinimalStringMass < 0.0 ) 323 if (MinimalStringMass < 0.0 ) return false; // For a string with di-quark having c or b quarks and s, c, b quarks
326 324
327 G4bool Result = G4UniformRand() < 325 G4bool Result = G4UniformRand() <
328 G4Exp(-0.66e-6*(string->Mass()*string- 326 G4Exp(-0.66e-6*(string->Mass()*string->Mass() - MinimalStringMass*MinimalStringMass));
329 // G4bool Result = string->Mas 327 // G4bool Result = string->Mass() < MinimalStringMass + 150.*MeV*G4UniformRand(); // a'la LUND
330 328
331 #ifdef debug_LUNDfragmentation 329 #ifdef debug_LUNDfragmentation
332 G4cout<<"StopFragmenting Minim 330 G4cout<<"StopFragmenting MinimalStringMass string->Mass() "<<MinimalStringMass
333 <<" "<<string->Mass()<<G 331 <<" "<<string->Mass()<<G4endl;
334 G4cout<<"StopFragmenting - Yes 332 G4cout<<"StopFragmenting - Yes/No "<<Result<<G4endl;
335 #endif 333 #endif
336 return Result; 334 return Result;
337 } 335 }
338 } 336 }
339 337
340 //-------------------------------------------- 338 //-----------------------------------------------------------------------------
341 339
342 G4KineticTrack * G4LundStringFragmentation::Sp 340 G4KineticTrack * G4LundStringFragmentation::Splitup(G4FragmentingString *string,
343 G4Fragmentin 341 G4FragmentingString *&newString)
344 { 342 {
345 #ifdef debug_LUNDfragmentation 343 #ifdef debug_LUNDfragmentation
346 G4cout<<G4endl; 344 G4cout<<G4endl;
347 G4cout<<"Start SplitUP ================ 345 G4cout<<"Start SplitUP ========================="<<G4endl;
348 G4cout<<"String partons: " <<string->Ge 346 G4cout<<"String partons: " <<string->GetLeftParton()->GetPDGEncoding()<<" "
349 <<string->Ge 347 <<string->GetRightParton()->GetPDGEncoding()<<" "
350 <<"Direction " <<string->Ge 348 <<"Direction " <<string->GetDecayDirection()<<G4endl;
351 #endif 349 #endif
352 350
353 //... random choice of string end to us 351 //... random choice of string end to use for creating the hadron (decay)
354 G4int SideOfDecay = (G4UniformRand() < 352 G4int SideOfDecay = (G4UniformRand() < 0.5)? 1: -1;
355 if (SideOfDecay < 0) 353 if (SideOfDecay < 0)
356 { 354 {
357 string->SetLeftPartonStable(); 355 string->SetLeftPartonStable();
358 } else 356 } else
359 { 357 {
360 string->SetRightPartonStable(); 358 string->SetRightPartonStable();
361 } 359 }
362 360
363 G4ParticleDefinition *newStringEnd; 361 G4ParticleDefinition *newStringEnd;
364 G4ParticleDefinition * HadronDefinition 362 G4ParticleDefinition * HadronDefinition;
365 363
366 G4double StringMass=string->Mass(); 364 G4double StringMass=string->Mass();
367 365
368 G4double ProbDqADq = GetDiquarkSuppress 366 G4double ProbDqADq = GetDiquarkSuppress();
369 G4double ProbSaS = 1.0 - 2.0 * GetStr 367 G4double ProbSaS = 1.0 - 2.0 * GetStrangeSuppress();
370 368
371 #ifdef debug_LUNDfragmentation 369 #ifdef debug_LUNDfragmentation
372 G4cout<<"StrMass DiquarkSuppression 370 G4cout<<"StrMass DiquarkSuppression "<<StringMass<<" "<<GetDiquarkSuppress()<<G4endl;
373 #endif 371 #endif
374 372
375 G4int NumberOfpossibleBaryons = 2; 373 G4int NumberOfpossibleBaryons = 2;
376 374
377 if (string->GetLeftParton()->GetParticl 375 if (string->GetLeftParton()->GetParticleSubType() != "quark") NumberOfpossibleBaryons++;
378 if (string->GetRightParton()->GetPartic 376 if (string->GetRightParton()->GetParticleSubType() != "quark") NumberOfpossibleBaryons++;
379 377
380 G4double ActualProb = ProbDqADq ; 378 G4double ActualProb = ProbDqADq ;
381 ActualProb *= (1.0-G4Pow::GetInstance() << 379 ActualProb *= (1.0-sqr(NumberOfpossibleBaryons*1400.0/StringMass));
382 if(ActualProb <0.0) ActualProb = 0.; <<
383 380
384 SetDiquarkSuppression(ActualProb); 381 SetDiquarkSuppression(ActualProb);
385 382
386 G4double Mth = 1250.0; 383 G4double Mth = 1250.0; // 2 Mk + Mpi
387 if ( NumberOfpossibleBaryons == 3 ){Mth 384 if ( NumberOfpossibleBaryons == 3 ){Mth = 2520.0;} // Mlambda/Msigma + Mk + Mpi
388 else if ( NumberOfpossibleBaryons == 4 385 else if ( NumberOfpossibleBaryons == 4 ){Mth = 2380.0;} // 2 Mlambda/Msigma + Mk + Mpi
389 else {} 386 else {}
390 387
391 ActualProb = ProbSaS; << 388 ActualProb = ProbSaS * (1.0 - G4Pow::GetInstance()->powA( Mth/StringMass , 4.0 ));
392 ActualProb *= (1.0 - G4Pow::GetInstance <<
393 if ( ActualProb < 0.0 ) ActualProb = 0. 389 if ( ActualProb < 0.0 ) ActualProb = 0.0;
394 SetStrangenessSuppression((1.0-ActualPr 390 SetStrangenessSuppression((1.0-ActualProb)/2.0);
395 391
396 #ifdef debug_LUNDfragmentation 392 #ifdef debug_LUNDfragmentation
397 G4cout<<"StrMass DiquarkSuppression cor 393 G4cout<<"StrMass DiquarkSuppression corrected "<<StringMass<<" "<<GetDiquarkSuppress()<<G4endl;
398 #endif 394 #endif
399 395
400 if (string->DecayIsQuark()) 396 if (string->DecayIsQuark())
401 { 397 {
402 HadronDefinition= QuarkSplitup(strin 398 HadronDefinition= QuarkSplitup(string->GetDecayParton(), newStringEnd);
403 } else { 399 } else {
404 HadronDefinition= DiQuarkSplitup(str 400 HadronDefinition= DiQuarkSplitup(string->GetDecayParton(), newStringEnd);
405 } 401 }
406 402
407 SetDiquarkSuppression(ProbDqADq); 403 SetDiquarkSuppression(ProbDqADq);
408 SetStrangenessSuppression((1.0-ProbSaS) 404 SetStrangenessSuppression((1.0-ProbSaS)/2.0);
409 405
410 if ( HadronDefinition == NULL ) { G4Kin 406 if ( HadronDefinition == NULL ) { G4KineticTrack * Hadron =0; return Hadron; }
411 407
412 #ifdef debug_LUNDfragmentation 408 #ifdef debug_LUNDfragmentation
413 G4cout<<"The parton "<<string->GetDecay 409 G4cout<<"The parton "<<string->GetDecayParton()->GetPDGEncoding()<<" "
414 <<" produces hadron "<<HadronDefi 410 <<" produces hadron "<<HadronDefinition->GetParticleName()
415 <<" and is transformed to "<<newS 411 <<" and is transformed to "<<newStringEnd->GetPDGEncoding()<<G4endl;
416 G4cout<<"The side of the string decay L 412 G4cout<<"The side of the string decay Left/Right (1/-1) "<<SideOfDecay<<G4endl;
417 #endif 413 #endif
418 // create new String from old, ie. keep 414 // create new String from old, ie. keep Left and Right order, but replace decay
419 415
420 if ( newString ) delete newString; 416 if ( newString ) delete newString;
421 417
422 newString=new G4FragmentingString(*stri 418 newString=new G4FragmentingString(*string,newStringEnd); // To store possible quark containt of new string
423 419
424 #ifdef debug_LUNDfragmentation 420 #ifdef debug_LUNDfragmentation
425 G4cout<<"An attempt to determine its en 421 G4cout<<"An attempt to determine its energy (SplitEandP)"<<G4endl;
426 #endif 422 #endif
427 G4LorentzVector* HadronMomentum=SplitEa 423 G4LorentzVector* HadronMomentum=SplitEandP(HadronDefinition, string, newString);
428 424
429 delete newString; newString=0; 425 delete newString; newString=0;
430 426
431 G4KineticTrack * Hadron =0; 427 G4KineticTrack * Hadron =0;
432 if ( HadronMomentum != 0 ) { 428 if ( HadronMomentum != 0 ) {
433 429
434 #ifdef debug_LUNDfragmentation 430 #ifdef debug_LUNDfragmentation
435 G4cout<<"The attempt was successful 431 G4cout<<"The attempt was successful"<<G4endl;
436 #endif 432 #endif
437 G4ThreeVector Pos; 433 G4ThreeVector Pos;
438 Hadron = new G4KineticTrack(HadronDefinit 434 Hadron = new G4KineticTrack(HadronDefinition, 0,Pos, *HadronMomentum);
439 435
440 if ( newString ) delete newString; 436 if ( newString ) delete newString;
441 437
442 newString=new G4FragmentingString(*string 438 newString=new G4FragmentingString(*string,newStringEnd,
443 HadronMomentum); 439 HadronMomentum);
444 delete HadronMomentum; 440 delete HadronMomentum;
445 } 441 }
446 else 442 else
447 { 443 {
448 #ifdef debug_LUNDfragmentation 444 #ifdef debug_LUNDfragmentation
449 G4cout<<"The attempt was not succes 445 G4cout<<"The attempt was not successful !!!"<<G4endl;
450 #endif 446 #endif
451 } 447 }
452 448
453 #ifdef debug_LUNDfragmentation 449 #ifdef debug_LUNDfragmentation
454 G4cout<<"End SplitUP (G4VLongitudinalSt 450 G4cout<<"End SplitUP (G4VLongitudinalStringDecay) ====================="<<G4endl;
455 #endif 451 #endif
456 452
457 return Hadron; 453 return Hadron;
458 } 454 }
459 455
460 //-------------------------------------------- 456 //-----------------------------------------------------------------------------
461 457
462 G4ParticleDefinition * G4LundStringFragmentati 458 G4ParticleDefinition * G4LundStringFragmentation::DiQuarkSplitup(G4ParticleDefinition* decay,
463 459 G4ParticleDefinition *&created)
464 { 460 {
465 G4double StrSup=GetStrangeSuppress(); 461 G4double StrSup=GetStrangeSuppress();
466 G4double ProbQQbar = (1.0 - 2.0*StrSup)*1.2 << 462 G4double ProbQQbar = (1.0 - 2.0*StrSup);
467 463
468 //... can Diquark break or not? 464 //... can Diquark break or not?
469 if (G4UniformRand() < DiquarkBreakProb ){ 465 if (G4UniformRand() < DiquarkBreakProb ){
470 466
471 //... Diquark break 467 //... Diquark break
472 G4int stableQuarkEncoding = decay->GetPD 468 G4int stableQuarkEncoding = decay->GetPDGEncoding()/1000;
473 G4int decayQuarkEncoding = (decay->GetPD 469 G4int decayQuarkEncoding = (decay->GetPDGEncoding()/100)%10;
474 if (G4UniformRand() < 0.5) 470 if (G4UniformRand() < 0.5)
475 { 471 {
476 G4int Swap = stableQuarkEncoding; 472 G4int Swap = stableQuarkEncoding;
477 stableQuarkEncoding = decayQuarkEncod 473 stableQuarkEncoding = decayQuarkEncoding;
478 decayQuarkEncoding = Swap; 474 decayQuarkEncoding = Swap;
479 } 475 }
480 476
481 G4int IsParticle=(decayQuarkEncoding>0) 477 G4int IsParticle=(decayQuarkEncoding>0) ? -1 : +1; // if we have a quark, we need antiquark
482 478
483 SetStrangenessSuppression((1.0-ProbQQbar <<
484 pDefPair QuarkPair = CreatePartonPair(Is 479 pDefPair QuarkPair = CreatePartonPair(IsParticle,false); // no diquarks wanted
485 SetStrangenessSuppression((1.0-StrSup)/2 <<
486 480
487 //... Build new Diquark 481 //... Build new Diquark
488 G4int QuarkEncoding=QuarkPair.second->Ge 482 G4int QuarkEncoding=QuarkPair.second->GetPDGEncoding();
489 G4int i10 = std::max(std::abs(QuarkEnco 483 G4int i10 = std::max(std::abs(QuarkEncoding), std::abs(stableQuarkEncoding));
490 G4int i20 = std::min(std::abs(QuarkEnco 484 G4int i20 = std::min(std::abs(QuarkEncoding), std::abs(stableQuarkEncoding));
491 G4int spin = (i10 != i20 && G4UniformRan 485 G4int spin = (i10 != i20 && G4UniformRand() <= 0.5)? 1 : 3;
492 G4int NewDecayEncoding = -1*IsParticle*( 486 G4int NewDecayEncoding = -1*IsParticle*(i10 * 1000 + i20 * 100 + spin);
493 created = FindParticle(NewDecayEncoding) 487 created = FindParticle(NewDecayEncoding);
494 G4ParticleDefinition * decayQuark=FindPa 488 G4ParticleDefinition * decayQuark=FindParticle(decayQuarkEncoding);
495 G4ParticleDefinition * had=hadronizer->B 489 G4ParticleDefinition * had=hadronizer->Build(QuarkPair.first, decayQuark);
496 StrangeSuppress=StrSup; 490 StrangeSuppress=StrSup;
497 491
498 return had; 492 return had;
499 493
500 } else { 494 } else {
501 //... Diquark does not break 495 //... Diquark does not break
502 496
503 G4int IsParticle=(decay->GetPDGEncoding( 497 G4int IsParticle=(decay->GetPDGEncoding()>0) ? +1 : -1; // if we have a diquark, we need quark
504 498
505 StrangeSuppress=(1.0 - ProbQQbar)/2.0; << 499 StrangeSuppress=(1.0 - ProbQQbar * 0.9)/2.0;
506 pDefPair QuarkPair = CreatePartonPair(Is 500 pDefPair QuarkPair = CreatePartonPair(IsParticle,false); // no diquarks wanted
507 501
508 created = QuarkPair.second; 502 created = QuarkPair.second;
509 503
510 G4ParticleDefinition * had=hadronizer->B 504 G4ParticleDefinition * had=hadronizer->Build(QuarkPair.first, decay);
511 StrangeSuppress=StrSup; 505 StrangeSuppress=StrSup;
512 506
513 return had; 507 return had;
514 } 508 }
515 } 509 }
516 510
517 //-------------------------------------------- 511 //-----------------------------------------------------------------------------
518 512
519 G4LorentzVector * G4LundStringFragmentation::S 513 G4LorentzVector * G4LundStringFragmentation::SplitEandP(G4ParticleDefinition * pHadron,
520 514 G4FragmentingString * string,
521 515 G4FragmentingString * newString)
522 { 516 {
523 G4LorentzVector String4Momentum=string->Get4 517 G4LorentzVector String4Momentum=string->Get4Momentum();
524 G4double StringMT2=string->MassT2(); 518 G4double StringMT2=string->MassT2();
525 G4double StringMT =std::sqrt(StringMT2); 519 G4double StringMT =std::sqrt(StringMT2);
526 520
527 G4double HadronMass = pHadron->GetPDGMass(); 521 G4double HadronMass = pHadron->GetPDGMass();
528 SetMinimalStringMass(newString); 522 SetMinimalStringMass(newString);
529 523
530 if ( MinimalStringMass < 0.0 ) return n 524 if ( MinimalStringMass < 0.0 ) return nullptr;
531 525
532 #ifdef debug_LUNDfragmentation 526 #ifdef debug_LUNDfragmentation
533 G4cout<<G4endl<<"Start LUND SplitEandP 527 G4cout<<G4endl<<"Start LUND SplitEandP "<<G4endl;
534 G4cout<<"String 4 mom, String M and Mt 528 G4cout<<"String 4 mom, String M and Mt "<<String4Momentum<<" "<<String4Momentum.mag()
535 <<" "<<std::sqrt(StringMT2)<<G4e 529 <<" "<<std::sqrt(StringMT2)<<G4endl;
536 G4cout<<"Hadron "<<pHadron->GetParticl 530 G4cout<<"Hadron "<<pHadron->GetParticleName()<<G4endl;
537 G4cout<<"HadM MinimalStringMassLeft St 531 G4cout<<"HadM MinimalStringMassLeft StringM hM+sM "<<HadronMass<<" "<<MinimalStringMass<<" "
538 <<String4Momentum.mag()<<" "<<Ha 532 <<String4Momentum.mag()<<" "<<HadronMass+MinimalStringMass<<G4endl;
539 #endif 533 #endif
540 534
541 if ((HadronMass + MinimalStringMass > string 535 if ((HadronMass + MinimalStringMass > string->Mass()) || MinimalStringMass < 0.)
542 { 536 {
543 #ifdef debug_LUNDfragmentation 537 #ifdef debug_LUNDfragmentation
544 G4cout<<"Mass of the string is not s 538 G4cout<<"Mass of the string is not sufficient to produce the hadron!"<<G4endl;
545 #endif 539 #endif
546 return 0; 540 return 0;
547 } // have to start all over! 541 } // have to start all over!
548 542
549 String4Momentum.setPz(0.); 543 String4Momentum.setPz(0.);
550 G4ThreeVector StringPt=String4Momentum.vect( 544 G4ThreeVector StringPt=String4Momentum.vect();
551 StringPt.setZ(0.); 545 StringPt.setZ(0.);
552 546
553 // calculate and assign hadron transverse mo 547 // calculate and assign hadron transverse momentum component HadronPx and HadronPy
554 G4ThreeVector HadronPt , RemSysPt; 548 G4ThreeVector HadronPt , RemSysPt;
555 G4double HadronMassT2, ResidualMassT2; 549 G4double HadronMassT2, ResidualMassT2;
556 G4double HadronMt, Pt, Pt2, phi; 550 G4double HadronMt, Pt, Pt2, phi;
557 551
558 G4double TmtCur = Tmt; 552 G4double TmtCur = Tmt;
559 <<
560 if ( (string->GetDecayParton()->GetPar 553 if ( (string->GetDecayParton()->GetParticleSubType()== "quark") &&
561 (pHadron->GetBaryonNumber() != 0) 554 (pHadron->GetBaryonNumber() != 0) ) {
562 TmtCur = Tmt*0.37; // q 555 TmtCur = Tmt*0.37; // q->B
563 } else if ( (string->GetDecayParton()- 556 } else if ( (string->GetDecayParton()->GetParticleSubType()== "quark") &&
564 (pHadron->GetBaryonNumber( 557 (pHadron->GetBaryonNumber() == 0) ) {
565 //TmtCur = Tmt; 558 //TmtCur = Tmt; // q->M
566 } else if ( (string->GetDecayParton()->GetPa 559 } else if ( (string->GetDecayParton()->GetParticleSubType()== "di_quark") &&
567 (pHadron->GetBaryonNumber( 560 (pHadron->GetBaryonNumber() == 0) ) {
568 //TmtCur = Tmt*0.89; 561 //TmtCur = Tmt*0.89; // qq -> M
569 } else if ( (string->GetDecayParton()- 562 } else if ( (string->GetDecayParton()->GetParticleSubType()== "di_quark") &&
570 (pHadron->GetBaryonNumber( 563 (pHadron->GetBaryonNumber() != 0) ) {
571 TmtCur = Tmt*1.35; 564 TmtCur = Tmt*1.35; // qq -> B
572 } 565 }
573 566
574 //... sample Pt of the hadron 567 //... sample Pt of the hadron
575 G4int attempt=0; 568 G4int attempt=0;
576 do 569 do
577 { 570 {
578 attempt++; if (attempt > StringLoopI 571 attempt++; if (attempt > StringLoopInterrupt) {return 0;}
579 572
580 HadronMt = HadronMass - TmtCur*G4Log 573 HadronMt = HadronMass - TmtCur*G4Log(G4UniformRand());
581 Pt2 = sqr(HadronMt)-sqr(HadronMass); Pt=st 574 Pt2 = sqr(HadronMt)-sqr(HadronMass); Pt=std::sqrt(Pt2);
582 phi = 2.*pi*G4UniformRand(); 575 phi = 2.*pi*G4UniformRand();
583 HadronPt = G4ThreeVector( Pt*std::co 576 HadronPt = G4ThreeVector( Pt*std::cos(phi), Pt*std::sin(phi), 0. );
584 RemSysPt = StringPt - HadronPt; 577 RemSysPt = StringPt - HadronPt;
585 HadronMassT2 = sqr(HadronMass) + Had 578 HadronMassT2 = sqr(HadronMass) + HadronPt.mag2();
586 ResidualMassT2=sqr(MinimalStringMass 579 ResidualMassT2=sqr(MinimalStringMass) + RemSysPt.mag2();
587 580
588 } while (std::sqrt(HadronMassT2) + std 581 } while (std::sqrt(HadronMassT2) + std::sqrt(ResidualMassT2) > StringMT);
589 582
590 //... sample z to define hadron longitudina 583 //... sample z to define hadron longitudinal momentum and energy
591 //... but first check the available phase sp 584 //... but first check the available phase space
592 585
593 G4double Pz2 = (sqr(StringMT2 - HadronMassT2 586 G4double Pz2 = (sqr(StringMT2 - HadronMassT2 - ResidualMassT2) -
594 4*HadronMassT2 * ResidualMassT2)/4./Stri 587 4*HadronMassT2 * ResidualMassT2)/4./StringMT2;
595 588
596 if (Pz2 < 0 ) {return 0;} // have t 589 if (Pz2 < 0 ) {return 0;} // have to start all over!
597 590
598 //... then compute allowed z region z_min < 591 //... then compute allowed z region z_min <= z <= z_max
599 592
600 G4double Pz = std::sqrt(Pz2); 593 G4double Pz = std::sqrt(Pz2);
601 G4double zMin = (std::sqrt(HadronMassT2+Pz2) 594 G4double zMin = (std::sqrt(HadronMassT2+Pz2) - Pz)/std::sqrt(StringMT2);
602 // G4double zMin = (std::sqrt(HadronMa 595 // G4double zMin = (std::sqrt(HadronMassT2+Pz2) - 0.)/std::sqrt(StringMT2); // For testing purposes
603 G4double zMax = (std::sqrt(HadronMassT2+Pz2) 596 G4double zMax = (std::sqrt(HadronMassT2+Pz2) + Pz)/std::sqrt(StringMT2);
604 597
605 if (zMin >= zMax) return 0; // have to sta 598 if (zMin >= zMax) return 0; // have to start all over!
606 599
607 G4double z = GetLightConeZ(zMin, zMax, 600 G4double z = GetLightConeZ(zMin, zMax,
608 string->GetDecayParton()->Get 601 string->GetDecayParton()->GetPDGEncoding(), pHadron,
609 HadronPt.x(), HadronPt.y()); 602 HadronPt.x(), HadronPt.y());
610 603
611 //... now compute hadron longitudinal moment 604 //... now compute hadron longitudinal momentum and energy
612 // longitudinal hadron momentum component Ha 605 // longitudinal hadron momentum component HadronPz
613 606
614 HadronPt.setZ(0.5* string->GetDecayDirection 607 HadronPt.setZ(0.5* string->GetDecayDirection() *
615 (z * string->LightConeDecay() - 608 (z * string->LightConeDecay() -
616 HadronMassT2/(z * string->LightConeD 609 HadronMassT2/(z * string->LightConeDecay())));
617 G4double HadronE = 0.5* (z * string->LightC 610 G4double HadronE = 0.5* (z * string->LightConeDecay() +
618 HadronMassT2/(z * string->LightConeD 611 HadronMassT2/(z * string->LightConeDecay()));
619 612
620 G4LorentzVector * a4Momentum= new G4LorentzV 613 G4LorentzVector * a4Momentum= new G4LorentzVector(HadronPt,HadronE);
621 614
622 #ifdef debug_LUNDfragmentation 615 #ifdef debug_LUNDfragmentation
623 G4cout<<G4endl<<" string->GetDecayDire 616 G4cout<<G4endl<<" string->GetDecayDirection() "<<string->GetDecayDirection()<<G4endl<<G4endl;
624 G4cout<<"string->LightConeDecay() "<<s 617 G4cout<<"string->LightConeDecay() "<<string->LightConeDecay()<<G4endl;
625 G4cout<<"HadronPt,HadronE "<<HadronPt< 618 G4cout<<"HadronPt,HadronE "<<HadronPt<<" "<<HadronE<<G4endl;
626 G4cout<<"String4Momentum "<<String4Mom 619 G4cout<<"String4Momentum "<<String4Momentum<<G4endl;
627 G4cout<<"Out of LUND SplitEandP "<<G4e 620 G4cout<<"Out of LUND SplitEandP "<<G4endl<<G4endl;
628 #endif 621 #endif
629 622
630 return a4Momentum; 623 return a4Momentum;
631 } 624 }
632 625
633 //-------------------------------------------- 626 //-----------------------------------------------------------------------------------------
634 627
635 G4double G4LundStringFragmentation::GetLightCo 628 G4double G4LundStringFragmentation::GetLightConeZ(G4double zmin, G4double zmax,
636 G4int PD 629 G4int PDGEncodingOfDecayParton,
637 G4Partic 630 G4ParticleDefinition* pHadron,
638 G4double 631 G4double Px, G4double Py)
639 { 632 {
640 G4double Mass = pHadron->GetPDGMass(); 633 G4double Mass = pHadron->GetPDGMass();
641 G4int HadronEncoding=std::abs(pHadron- 634 G4int HadronEncoding=std::abs(pHadron->GetPDGEncoding());
642 635
643 G4double Mt2 = Px*Px + Py*Py + Mass*Mass; 636 G4double Mt2 = Px*Px + Py*Py + Mass*Mass;
644 637
645 G4double Alund, Blund; 638 G4double Alund, Blund;
646 G4double zOfMaxyf(0.), maxYf(1.), z(0.), yf( 639 G4double zOfMaxyf(0.), maxYf(1.), z(0.), yf(1.);
647 640
648 if (!((std::abs(PDGEncodingOfDecayParton) > 641 if (!((std::abs(PDGEncodingOfDecayParton) > 1000) && (HadronEncoding > 1000)) )
649 { // ---------------- Quark fragmentation 642 { // ---------------- Quark fragmentation and qq-> meson ----------------------
650 Alund=1.; 643 Alund=1.;
651 Blund=0.7/GeV/GeV; 644 Blund=0.7/GeV/GeV;
652 645
653 G4double BMt2 = Blund*Mt2; 646 G4double BMt2 = Blund*Mt2;
654 if (Alund == 1.0) { 647 if (Alund == 1.0) {
655 zOfMaxyf=BMt2/(Blund*Mt2 + 1.);} 648 zOfMaxyf=BMt2/(Blund*Mt2 + 1.);}
656 else { 649 else {
657 zOfMaxyf = ((1.0+BMt2) - std::sqrt(sq 650 zOfMaxyf = ((1.0+BMt2) - std::sqrt(sqr(1.0-BMt2) + 4.0*BMt2*Alund))/2.0/(1.-Alund);
658 } 651 }
659 652
660 if (zOfMaxyf < zmin) {zOfMaxyf=zmin;} 653 if (zOfMaxyf < zmin) {zOfMaxyf=zmin;}
661 if (zOfMaxyf > zmax) {zOfMaxyf=zmax;} 654 if (zOfMaxyf > zmax) {zOfMaxyf=zmax;}
662 maxYf=(1-zOfMaxyf)/zOfMaxyf * G4Exp(-Blun 655 maxYf=(1-zOfMaxyf)/zOfMaxyf * G4Exp(-Blund*Mt2/zOfMaxyf);
663 656
664 const G4int maxNumberOfLoops = 1000 657 const G4int maxNumberOfLoops = 1000;
665 G4int loopCounter = 0; 658 G4int loopCounter = 0;
666 do 659 do
667 { 660 {
668 z = zmin + G4UniformRand()*(zmax-zmin); 661 z = zmin + G4UniformRand()*(zmax-zmin);
669 //yf = (1-z)/z * G4Exp(-Blund* 662 //yf = (1-z)/z * G4Exp(-Blund*Mt2/z);
670 yf = G4Pow::GetInstance()->powA(1.0-z,Alun 663 yf = G4Pow::GetInstance()->powA(1.0-z,Alund)/z*G4Exp(-BMt2/z);
671 } 664 }
672 while ( (G4UniformRand()*maxYf > yf) && + 665 while ( (G4UniformRand()*maxYf > yf) && ++loopCounter < maxNumberOfLoops );
673 if ( loopCounter >= maxNumberOfLoop 666 if ( loopCounter >= maxNumberOfLoops ) {
674 z = 0.5*(zmin + zmax); // Just a 667 z = 0.5*(zmin + zmax); // Just a value between zmin and zmax, no physics considerations at all!
675 } 668 }
676 return z; 669 return z;
677 } 670 }
678 671
679 if (std::abs(PDGEncodingOfDecayParton) > 100 672 if (std::abs(PDGEncodingOfDecayParton) > 1000)
680 { 673 {
681 G4double an = 2.5; 674 G4double an = 2.5;
682 an +=(sqr(Px)+sqr(Py))/sqr(GeV 675 an +=(sqr(Px)+sqr(Py))/sqr(GeV)-0.5;
683 z=zmin + (zmax-zmin)*G4Pow::Ge 676 z=zmin + (zmax-zmin)*G4Pow::GetInstance()->powA(G4UniformRand(),1./an);
684 if( PDGEncodingOfDecayParton > <<
685 } 677 }
686 678
687 return z; 679 return z;
688 } 680 }
689 681
690 //-------------------------------------------- 682 //----------------------------------------------------------------------------------------------------------
691 683
692 G4bool G4LundStringFragmentation::SplitLast(G4 684 G4bool G4LundStringFragmentation::SplitLast(G4FragmentingString * string,
693 G4 685 G4KineticTrackVector * LeftVector,
694 G4 686 G4KineticTrackVector * RightVector)
695 { 687 {
696 //... perform last cluster decay 688 //... perform last cluster decay
697 SetMinimalStringMass( string); 689 SetMinimalStringMass( string);
698 if ( MinimalStringMass < 0.) return fa 690 if ( MinimalStringMass < 0.) return false;
699 #ifdef debug_LUNDfragmentation 691 #ifdef debug_LUNDfragmentation
700 G4cout<<G4endl<<"Split last----------- 692 G4cout<<G4endl<<"Split last-----------------------------------------"<<G4endl;
701 G4cout<<"MinimalStringMass "<<MinimalS 693 G4cout<<"MinimalStringMass "<<MinimalStringMass<<G4endl;
702 G4cout<<"Left "<<string->GetLeftParto 694 G4cout<<"Left "<<string->GetLeftParton()->GetPDGEncoding()<<" "<<string->GetPleft()<<G4endl;
703 G4cout<<"Right "<<string->GetRightPart 695 G4cout<<"Right "<<string->GetRightParton()->GetPDGEncoding()<<" "<<string->GetPright()<<G4endl;
704 G4cout<<"String4mom "<<string->GetPstr 696 G4cout<<"String4mom "<<string->GetPstring()<<" "<<string->GetPstring().mag()<<G4endl;
705 #endif 697 #endif
706 698
707 G4LorentzVector Str4Mom=string->Get4Mo 699 G4LorentzVector Str4Mom=string->Get4Momentum();
708 G4LorentzRotation toCms(-1*Str4Mom.boo 700 G4LorentzRotation toCms(-1*Str4Mom.boostVector());
709 G4LorentzVector Pleft = toCms * string 701 G4LorentzVector Pleft = toCms * string->GetPleft();
710 toCms.rotateZ(-1*Pleft.phi()); 702 toCms.rotateZ(-1*Pleft.phi());
711 toCms.rotateY(-1*Pleft.theta()); 703 toCms.rotateY(-1*Pleft.theta());
712 704
713 G4LorentzRotation toObserverFrame= toC 705 G4LorentzRotation toObserverFrame= toCms.inverse();
714 706
715 G4double StringMass=string->Mass(); 707 G4double StringMass=string->Mass();
716 708
717 G4ParticleDefinition * LeftHadron(0), * Righ 709 G4ParticleDefinition * LeftHadron(0), * RightHadron(0);
718 710
719 NumberOf_FS=0; 711 NumberOf_FS=0;
720 for (G4int i=0; i<350; i++) {FS_Weight[i]=0. 712 for (G4int i=0; i<350; i++) {FS_Weight[i]=0.;}
721 713
722 G4int sampledState = 0; 714 G4int sampledState = 0;
723 715
724 #ifdef debug_LUNDfragmentation 716 #ifdef debug_LUNDfragmentation
725 G4cout<<"StrMass "<<StringMass<<" q's 717 G4cout<<"StrMass "<<StringMass<<" q's "
726 <<string->GetLeftParton()->GetPa 718 <<string->GetLeftParton()->GetParticleName()<<" "
727 <<string->GetRightParton()->GetP 719 <<string->GetRightParton()->GetParticleName()<<G4endl;
728 #endif 720 #endif
729 721
730 string->SetLeftPartonStable(); // to query q 722 string->SetLeftPartonStable(); // to query quark contents..
731 723
732 if (string->IsAFourQuarkString() ) 724 if (string->IsAFourQuarkString() )
733 { 725 {
734 G4int IDleft =std::abs(string->GetLe <<
735 G4int IDright=std::abs(string->GetRi <<
736 <<
737 if ( (IDleft > 3000) || (IDright > 3 <<
738 if ( ! Diquark_AntiDiquark_belowTh <<
739 { <<
740 return false; <<
741 } <<
742 } else { <<
743 // The string is qq-qqbar type. Diquarks a 726 // The string is qq-qqbar type. Diquarks are on the string ends
744 if (StringMass-MinimalStringMass < 0 727 if (StringMass-MinimalStringMass < 0.)
745 { 728 {
746 if (! Diquark_AntiDiquark_belowThreshold 729 if (! Diquark_AntiDiquark_belowThreshold_lastSplitting(string, LeftHadron, RightHadron) )
747 { 730 {
748 return false; 731 return false;
749 } 732 }
750 } else 733 } else
751 { 734 {
752 Diquark_AntiDiquark_aboveThreshold_lastS 735 Diquark_AntiDiquark_aboveThreshold_lastSplitting(string, LeftHadron, RightHadron);
753 736
754 if (NumberOf_FS == 0) return false; 737 if (NumberOf_FS == 0) return false;
755 738
756 sampledState = SampleS 739 sampledState = SampleState();
757 if (string->GetLeftParton()->GetPDGEncod 740 if (string->GetLeftParton()->GetPDGEncoding() < 0)
758 { 741 {
759 LeftHadron =FS_LeftHadron[sampledState 742 LeftHadron =FS_LeftHadron[sampledState];
760 RightHadron=FS_RightHadron[sampledStat 743 RightHadron=FS_RightHadron[sampledState];
761 } else 744 } else
762 { 745 {
763 LeftHadron =FS_RightHadron[sampledStat 746 LeftHadron =FS_RightHadron[sampledState];
764 RightHadron=FS_LeftHadron[sampledState 747 RightHadron=FS_LeftHadron[sampledState];
765 } 748 }
766 } 749 }
767 } // ID > 3300 << 750 } else
768 } else { << 751 {
769 if (string->DecayIsQuark() && string->Stab 752 if (string->DecayIsQuark() && string->StableIsQuark() )
770 { //... there are quarks on cluster 753 { //... there are quarks on cluster ends
771 #ifdef debug_LUNDfragm 754 #ifdef debug_LUNDfragmentation
772 G4cout<<"Q Q string La 755 G4cout<<"Q Q string LastSplit"<<G4endl;
773 #endif 756 #endif
774 757
775 Quark_AntiQuark_lastSplitting(string, Le 758 Quark_AntiQuark_lastSplitting(string, LeftHadron, RightHadron);
776 759
777 if (NumberOf_FS == 0) return false; 760 if (NumberOf_FS == 0) return false;
778 sampledState = SampleState() 761 sampledState = SampleState();
779 if (string->GetLeftParton()->GetPDGEncod 762 if (string->GetLeftParton()->GetPDGEncoding() < 0)
780 { 763 {
781 LeftHadron =FS_RightHadron[sampledStat 764 LeftHadron =FS_RightHadron[sampledState];
782 RightHadron=FS_LeftHadron[sampledState 765 RightHadron=FS_LeftHadron[sampledState];
783 } else 766 } else
784 { 767 {
785 LeftHadron =FS_LeftHadron[sampledState 768 LeftHadron =FS_LeftHadron[sampledState];
786 RightHadron=FS_RightHadron[sampledStat 769 RightHadron=FS_RightHadron[sampledState];
787 } 770 }
788 } else 771 } else
789 { //... there is a Diquark on one of 772 { //... there is a Diquark on one of the cluster ends
790 #ifdef debug_LUNDfragm 773 #ifdef debug_LUNDfragmentation
791 G4cout<<"DiQ Q string 774 G4cout<<"DiQ Q string Last Split"<<G4endl;
792 #endif 775 #endif
793 776
794 Quark_Diquark_lastSplitting(string, Left 777 Quark_Diquark_lastSplitting(string, LeftHadron, RightHadron);
795 778
796 if (NumberOf_FS == 0) return false; 779 if (NumberOf_FS == 0) return false;
797 sampledState = SampleState() 780 sampledState = SampleState();
798 781
799 if (string->GetLeftParton()->GetParticle 782 if (string->GetLeftParton()->GetParticleSubType() == "quark")
800 { 783 {
801 LeftHadron =FS_LeftHadron[sampledState 784 LeftHadron =FS_LeftHadron[sampledState];
802 RightHadron=FS_RightHadron[sampledStat 785 RightHadron=FS_RightHadron[sampledState];
803 } else 786 } else
804 { 787 {
805 LeftHadron =FS_RightHadron[sampledStat 788 LeftHadron =FS_RightHadron[sampledState];
806 RightHadron=FS_LeftHadron[sampledState 789 RightHadron=FS_LeftHadron[sampledState];
807 } 790 }
808 } 791 }
809 792
810 } 793 }
811 794
812 #ifdef debug_LUNDfragmentation 795 #ifdef debug_LUNDfragmentation
813 G4cout<<"Sampled hadrons: "<<LeftHadro 796 G4cout<<"Sampled hadrons: "<<LeftHadron->GetParticleName()<<" "<<RightHadron->GetParticleName()<<G4endl;
814 #endif 797 #endif
815 798
816 G4LorentzVector P_left =string->GetPleft(), 799 G4LorentzVector P_left =string->GetPleft(), P_right = string->GetPright();
817 800
818 G4LorentzVector LeftMom, RightMom; 801 G4LorentzVector LeftMom, RightMom;
819 G4ThreeVector Pos; 802 G4ThreeVector Pos;
820 803
821 Sample4Momentum(&LeftMom, LeftHadron->GetPD 804 Sample4Momentum(&LeftMom, LeftHadron->GetPDGMass(),
822 &RightMom, RightHadron->GetPDGMass(), 805 &RightMom, RightHadron->GetPDGMass(),
823 StringMass); 806 StringMass);
824 807
825 // Sample4Momentum ascribes LeftMom.pz 808 // Sample4Momentum ascribes LeftMom.pz() along positive Z axis for baryons in many cases.
826 // It must be negative in case when th 809 // It must be negative in case when the system is moving against Z axis.
827 810
828 if (!(string->DecayIsQuark() && string->Stab 811 if (!(string->DecayIsQuark() && string->StableIsQuark() ))
829 { // Only for qq - q, q - qq, and qq - qqbar 812 { // Only for qq - q, q - qq, and qq - qqbar -------------------
830 813
831 if ( G4UniformRand() <= 0.5 ) 814 if ( G4UniformRand() <= 0.5 )
832 { 815 {
833 if (P_left.z() <= 0.) {G4LorentzVector t 816 if (P_left.z() <= 0.) {G4LorentzVector tmp = LeftMom; LeftMom=RightMom; RightMom=tmp;}
834 } 817 }
835 else 818 else
836 { 819 {
837 if (P_right.z() >= 0.) {G4LorentzVector 820 if (P_right.z() >= 0.) {G4LorentzVector tmp = LeftMom; LeftMom=RightMom; RightMom=tmp;}
838 } 821 }
839 } 822 }
840 823
841 LeftMom *=toObserverFrame; 824 LeftMom *=toObserverFrame;
842 RightMom*=toObserverFrame; 825 RightMom*=toObserverFrame;
843 826
844 LeftVector->push_back(new G4KineticTrack(Lef 827 LeftVector->push_back(new G4KineticTrack(LeftHadron, 0, Pos, LeftMom));
845 RightVector->push_back(new G4KineticTrack(Ri 828 RightVector->push_back(new G4KineticTrack(RightHadron, 0, Pos, RightMom));
846 829
847 string->LorentzRotate(toObserverFrame); 830 string->LorentzRotate(toObserverFrame);
848 return true; 831 return true;
849 } 832 }
850 833
851 //-------------------------------------------- 834 //----------------------------------------------------------------------------------------
852 835
853 G4bool G4LundStringFragmentation:: 836 G4bool G4LundStringFragmentation::
854 Diquark_AntiDiquark_belowThreshold_lastSplitti 837 Diquark_AntiDiquark_belowThreshold_lastSplitting(G4FragmentingString * & string,
855 838 G4ParticleDefinition * & LeftHadron,
856 839 G4ParticleDefinition * & RightHadron)
857 { 840 {
858 G4double StringMass = string->Mass(); 841 G4double StringMass = string->Mass();
859 842
860 G4int cClusterInterrupt = 0; 843 G4int cClusterInterrupt = 0;
861 G4bool isOK = false; 844 G4bool isOK = false;
862 do 845 do
863 { 846 {
864 G4int LeftQuark1= string->GetLeftParton()- 847 G4int LeftQuark1= string->GetLeftParton()->GetPDGEncoding()/1000;
865 G4int LeftQuark2=(string->GetLeftParton()- 848 G4int LeftQuark2=(string->GetLeftParton()->GetPDGEncoding()/100)%10;
866 849
867 G4int RightQuark1= string->GetRightParton( 850 G4int RightQuark1= string->GetRightParton()->GetPDGEncoding()/1000;
868 G4int RightQuark2=(string->GetRightParton( 851 G4int RightQuark2=(string->GetRightParton()->GetPDGEncoding()/100)%10;
869 852
870 if (G4UniformRand()<0.5) 853 if (G4UniformRand()<0.5)
871 { 854 {
872 LeftHadron =hadronizer->Build(FindPartic 855 LeftHadron =hadronizer->Build(FindParticle( LeftQuark1),
873 FindParticle(RightQuark1)); 856 FindParticle(RightQuark1));
874 RightHadron= (LeftHadron == nullptr) ? n 857 RightHadron= (LeftHadron == nullptr) ? nullptr :
875 858 hadronizer->Build(FindParticle( LeftQuark2),
876 FindParticle(RightQuark2)); 859 FindParticle(RightQuark2));
877 } else 860 } else
878 { 861 {
879 LeftHadron =hadronizer->Build(FindPartic 862 LeftHadron =hadronizer->Build(FindParticle( LeftQuark1),
880 FindParticle(RightQuark2)); 863 FindParticle(RightQuark2));
881 RightHadron=(LeftHadron == nullptr) ? nu 864 RightHadron=(LeftHadron == nullptr) ? nullptr :
882 hadronizer 865 hadronizer->Build(FindParticle( LeftQuark2),
883 FindParticle(RightQuark1)); 866 FindParticle(RightQuark1));
884 } 867 }
885 868
886 isOK = (LeftHadron != nullptr) && (RightHa 869 isOK = (LeftHadron != nullptr) && (RightHadron != nullptr);
887 <<
888 if(isOK) { isOK = (StringMass > LeftHadron 870 if(isOK) { isOK = (StringMass > LeftHadron->GetPDGMass() + RightHadron->GetPDGMass()); }
889 ++cClusterInterrupt; 871 ++cClusterInterrupt;
890 //... repeat procedure, if mass of cluster 872 //... repeat procedure, if mass of cluster is too low to produce hadrons
891 //... ClusterMassCut = 0.15*GeV model para 873 //... ClusterMassCut = 0.15*GeV model parameter
892 } 874 }
893 while (isOK == false && cClusterInterrupt < << 875 while (isOK == false || cClusterInterrupt < ClusterLoopInterrupt);
894 /* Loop checking, 07.08.2015, A.Ribon */ 876 /* Loop checking, 07.08.2015, A.Ribon */
895 return isOK; 877 return isOK;
896 } 878 }
897 879
898 //-------------------------------------------- 880 //----------------------------------------------------------------------------------------
899 881
900 G4bool G4LundStringFragmentation:: 882 G4bool G4LundStringFragmentation::
901 Diquark_AntiDiquark_aboveThreshold_lastSplitti 883 Diquark_AntiDiquark_aboveThreshold_lastSplitting(G4FragmentingString * & string,
902 884 G4ParticleDefinition * & LeftHadron,
903 885 G4ParticleDefinition * & RightHadron)
904 { 886 {
905 // StringMass-MinimalStringMass > 0. Creatio 887 // StringMass-MinimalStringMass > 0. Creation of 2 baryons is possible ----
906 888
907 G4double StringMass = string->Mass(); 889 G4double StringMass = string->Mass();
908 G4double StringMassSqr= sqr(StringMass); 890 G4double StringMassSqr= sqr(StringMass);
909 G4ParticleDefinition * Di_Quark; 891 G4ParticleDefinition * Di_Quark;
910 G4ParticleDefinition * Anti_Di_Quark; 892 G4ParticleDefinition * Anti_Di_Quark;
911 893
912 if (string->GetLeftParton()->GetPDGEncoding( 894 if (string->GetLeftParton()->GetPDGEncoding() < 0)
913 { 895 {
914 Anti_Di_Quark =string->GetLeftParton(); 896 Anti_Di_Quark =string->GetLeftParton();
915 Di_Quark=string->GetRightParton(); 897 Di_Quark=string->GetRightParton();
916 } else 898 } else
917 { 899 {
918 Anti_Di_Quark =string->GetRightParton(); 900 Anti_Di_Quark =string->GetRightParton();
919 Di_Quark=string->GetLeftParton(); 901 Di_Quark=string->GetLeftParton();
920 } 902 }
921 903
922 G4int IDAnti_di_quark =Anti_Di_Quark->Get 904 G4int IDAnti_di_quark =Anti_Di_Quark->GetPDGEncoding();
923 G4int AbsIDAnti_di_quark =std::abs(IDAnti_di 905 G4int AbsIDAnti_di_quark =std::abs(IDAnti_di_quark);
924 G4int IDdi_quark =Di_Quark->GetPDGEn 906 G4int IDdi_quark =Di_Quark->GetPDGEncoding();
925 G4int AbsIDdi_quark =std::abs(IDdi_quar 907 G4int AbsIDdi_quark =std::abs(IDdi_quark);
926 908
927 G4int ADi_q1=AbsIDAnti_di_quark/1000; 909 G4int ADi_q1=AbsIDAnti_di_quark/1000;
928 G4int ADi_q2=(AbsIDAnti_di_quark-ADi_q1*1000 910 G4int ADi_q2=(AbsIDAnti_di_quark-ADi_q1*1000)/100;
929 911
930 G4int Di_q1=AbsIDdi_quark/1000; 912 G4int Di_q1=AbsIDdi_quark/1000;
931 G4int Di_q2=(AbsIDdi_quark-Di_q1*1000)/100; 913 G4int Di_q2=(AbsIDdi_quark-Di_q1*1000)/100;
932 914
933 NumberOf_FS=0; 915 NumberOf_FS=0;
934 for (G4int ProdQ=1; ProdQ < 6; ProdQ++) 916 for (G4int ProdQ=1; ProdQ < 6; ProdQ++)
935 { 917 {
936 G4int StateADiQ=0; 918 G4int StateADiQ=0;
937 const G4int maxNumberOfLoops = 919 const G4int maxNumberOfLoops = 1000;
938 G4int loopCounter = 0; 920 G4int loopCounter = 0;
939 do // while(Meson[AbsIDquark-1][ProdQ-1][ 921 do // while(Meson[AbsIDquark-1][ProdQ-1][StateQ]<>0);
940 { 922 {
941 LeftHadron=G4ParticleTable::GetParticleT 923 LeftHadron=G4ParticleTable::GetParticleTable()->FindParticle(
942 -Baryon[ADi_q1-1][ADi_q2-1][Prod 924 -Baryon[ADi_q1-1][ADi_q2-1][ProdQ-1][StateADiQ]);
943 925
944 if (LeftHadron == NULL) continue; 926 if (LeftHadron == NULL) continue;
945 G4double LeftHadronMass=LeftHadron->GetP 927 G4double LeftHadronMass=LeftHadron->GetPDGMass();
946 928
947 G4int StateDiQ=0; 929 G4int StateDiQ=0;
948 const G4int maxNumberO 930 const G4int maxNumberOfInternalLoops = 1000;
949 G4int internalLoopCoun 931 G4int internalLoopCounter = 0;
950 do // while(Baryon[Di_q1-1][Di_q2-1][Pro 932 do // while(Baryon[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]<>0);
951 { 933 {
952 RightHadron=G4ParticleTable::GetPartic 934 RightHadron=G4ParticleTable::GetParticleTable()->FindParticle(
953 +Baryon[Di_q1-1][Di_q2-1][Prod 935 +Baryon[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]);
954 936
955 if (RightHadron == NULL) continue; 937 if (RightHadron == NULL) continue;
956 G4double RightHadronMass=RightHadron-> 938 G4double RightHadronMass=RightHadron->GetPDGMass();
957 939
958 if (StringMass > LeftHadronMass + Righ 940 if (StringMass > LeftHadronMass + RightHadronMass)
959 { 941 {
960 if ( N 942 if ( NumberOf_FS > 349 ) {
961 G4Ex 943 G4ExceptionDescription ed;
962 ed < 944 ed << " NumberOf_FS exceeds its limit: NumberOf_FS=" << NumberOf_FS << G4endl;
963 G4Ex 945 G4Exception( "G4LundStringFragmentation::Diquark_AntiDiquark_aboveThreshold_lastSplitting ",
964 946 "HAD_LUND_001", JustWarning, ed );
965 Numb 947 NumberOf_FS = 349;
966 } 948 }
967 949
968 G4double FS_Psqr=lambda(StringMassSq 950 G4double FS_Psqr=lambda(StringMassSqr,sqr(LeftHadronMass),
969 sqr(RightHadronMass)); 951 sqr(RightHadronMass));
970 //FS_Psqr=1.; 952 //FS_Psqr=1.;
971 FS_Weight[NumberOf_FS]=std::sqrt(FS_ 953 FS_Weight[NumberOf_FS]=std::sqrt(FS_Psqr)*FS_Psqr*
972 BaryonWeight[ADi_q1-1][AD 954 BaryonWeight[ADi_q1-1][ADi_q2-1][ProdQ-1][StateADiQ]*
973 BaryonWeight[Di_q1-1][Di_ 955 BaryonWeight[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]*
974 Prob_QQbar[ProdQ-1]; 956 Prob_QQbar[ProdQ-1];
975 957
976 FS_LeftHadron[NumberOf_FS] = LeftHad 958 FS_LeftHadron[NumberOf_FS] = LeftHadron;
977 FS_RightHadron[NumberOf_FS]= RightHa 959 FS_RightHadron[NumberOf_FS]= RightHadron;
978 NumberOf_FS++; 960 NumberOf_FS++;
979 } // End of if (StringMass > LeftHadro 961 } // End of if (StringMass > LeftHadronMass + RightHadronMass)
980 962
981 StateDiQ++; 963 StateDiQ++;
982 964
983 } while( (Baryon[Di_q1-1][Di_q2-1][ProdQ 965 } while( (Baryon[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]!=0) &&
984 ++internalLoo 966 ++internalLoopCounter < maxNumberOfInternalLoops );
985 if ( internalLoopCount 967 if ( internalLoopCounter >= maxNumberOfInternalLoops ) {
986 return false; 968 return false;
987 } 969 }
988 970
989 StateADiQ++; 971 StateADiQ++;
990 } while( (Baryon[ADi_q1-1][ADi_q2-1][ProdQ 972 } while( (Baryon[ADi_q1-1][ADi_q2-1][ProdQ-1][StateADiQ]!=0) &&
991 ++loopCounter < maxNu 973 ++loopCounter < maxNumberOfLoops );
992 if ( loopCounter >= maxNumberO 974 if ( loopCounter >= maxNumberOfLoops ) {
993 return false; 975 return false;
994 } 976 }
995 } // End of for (G4int ProdQ=1; ProdQ < 4; P 977 } // End of for (G4int ProdQ=1; ProdQ < 4; ProdQ++)
996 978
997 return true; 979 return true;
998 } 980 }
999 981
1000 //------------------------------------------- 982 //----------------------------------------------------------------------------------------
1001 983
1002 G4bool G4LundStringFragmentation::Quark_Diqua 984 G4bool G4LundStringFragmentation::Quark_Diquark_lastSplitting(G4FragmentingString * & string,
1003 985 G4ParticleDefinition * & LeftHadron,
1004 986 G4ParticleDefinition * & RightHadron)
1005 { 987 {
1006 G4double StringMass = string->Mass(); 988 G4double StringMass = string->Mass();
1007 G4double StringMassSqr= sqr(StringMass); 989 G4double StringMassSqr= sqr(StringMass);
1008 990
1009 G4ParticleDefinition * Di_Quark; 991 G4ParticleDefinition * Di_Quark;
1010 G4ParticleDefinition * Quark; 992 G4ParticleDefinition * Quark;
1011 993
1012 if (string->GetLeftParton()->GetParticleSub 994 if (string->GetLeftParton()->GetParticleSubType()== "quark")
1013 { 995 {
1014 Quark =string->GetLeftParton(); 996 Quark =string->GetLeftParton();
1015 Di_Quark=string->GetRightParton(); 997 Di_Quark=string->GetRightParton();
1016 } else 998 } else
1017 { 999 {
1018 Quark =string->GetRightParton(); 1000 Quark =string->GetRightParton();
1019 Di_Quark=string->GetLeftParton(); 1001 Di_Quark=string->GetLeftParton();
1020 } 1002 }
1021 1003
1022 G4int IDquark =Quark->GetPDGEncoding 1004 G4int IDquark =Quark->GetPDGEncoding();
1023 G4int AbsIDquark =std::abs(IDquark); 1005 G4int AbsIDquark =std::abs(IDquark);
1024 G4int IDdi_quark =Di_Quark->GetPDGEncodin 1006 G4int IDdi_quark =Di_Quark->GetPDGEncoding();
1025 G4int AbsIDdi_quark=std::abs(IDdi_quark); 1007 G4int AbsIDdi_quark=std::abs(IDdi_quark);
1026 G4int Di_q1=AbsIDdi_quark/1000; 1008 G4int Di_q1=AbsIDdi_quark/1000;
1027 G4int Di_q2=(AbsIDdi_quark-Di_q1*1000)/100; 1009 G4int Di_q2=(AbsIDdi_quark-Di_q1*1000)/100;
1028 G4int SignDiQ= 1; 1010 G4int SignDiQ= 1;
1029 if (IDdi_quark < 0) SignDiQ=-1; 1011 if (IDdi_quark < 0) SignDiQ=-1;
1030 1012
1031 NumberOf_FS=0; 1013 NumberOf_FS=0;
1032 for (G4int ProdQ=1; ProdQ < 4; ProdQ++) // 1014 for (G4int ProdQ=1; ProdQ < 4; ProdQ++) // Loop over quark-antiquark cases: u-ubar, d-dbar, s-sbar
1033 { // 1015 { // (as last splitting, do not consider c-cbar and b-bbar cases)
1034 G4int SignQ; 1016 G4int SignQ;
1035 if (IDquark > 0) 1017 if (IDquark > 0)
1036 { 1018 {
1037 SignQ=-1; 1019 SignQ=-1;
1038 if (IDquark == 2) 1020 if (IDquark == 2) SignQ= 1;
1039 if ((IDquark == 1) && (ProdQ == 3)) Sig 1021 if ((IDquark == 1) && (ProdQ == 3)) SignQ= 1; // K0
1040 if ((IDquark == 3) && (ProdQ == 1)) Sig 1022 if ((IDquark == 3) && (ProdQ == 1)) SignQ=-1; // K0bar
1041 if (IDquark == 4) 1023 if (IDquark == 4) SignQ= 1; // D+, D0, Ds+
1042 if (IDquark == 5) Sig 1024 if (IDquark == 5) SignQ=-1; // B-, anti_B0, anti_Bs0
1043 } else 1025 } else
1044 { 1026 {
1045 SignQ= 1; 1027 SignQ= 1;
1046 if (IDquark == -2) Sig 1028 if (IDquark == -2) SignQ=-1;
1047 if ((IDquark ==-1) && (ProdQ == 3)) Sig 1029 if ((IDquark ==-1) && (ProdQ == 3)) SignQ=-1; // K0bar
1048 if ((IDquark ==-3) && (ProdQ == 1)) Sig 1030 if ((IDquark ==-3) && (ProdQ == 1)) SignQ= 1; // K0
1049 if (IDquark == -4) Sig 1031 if (IDquark == -4) SignQ=-1; // D-, anti_D0, anti_Ds+
1050 if (IDquark == -5) Sig 1032 if (IDquark == -5) SignQ= 1; // B+, B0, Bs0
1051 } 1033 }
1052 1034
1053 if (AbsIDquark == ProdQ) SignQ 1035 if (AbsIDquark == ProdQ) SignQ= 1;
1054 1036
1055 G4int StateQ=0; 1037 G4int StateQ=0;
1056 const G4int maxNumberOfLoops 1038 const G4int maxNumberOfLoops = 1000;
1057 G4int loopCounter = 0; 1039 G4int loopCounter = 0;
1058 do // while(Meson[AbsIDquark-1][ProdQ-1] 1040 do // while(Meson[AbsIDquark-1][ProdQ-1][StateQ]<>0);
1059 { 1041 {
1060 LeftHadron=G4ParticleTable::GetParticle 1042 LeftHadron=G4ParticleTable::GetParticleTable()->FindParticle(SignQ*
1061 Meson[AbsIDquark-1][ProdQ-1][St 1043 Meson[AbsIDquark-1][ProdQ-1][StateQ]);
1062 if (LeftHadron == NULL) continue; 1044 if (LeftHadron == NULL) continue;
1063 G4double LeftHadronMass=LeftHadron->Get 1045 G4double LeftHadronMass=LeftHadron->GetPDGMass();
1064 1046
1065 G4int StateDiQ=0; 1047 G4int StateDiQ=0;
1066 const G4int maxNumber 1048 const G4int maxNumberOfInternalLoops = 1000;
1067 G4int internalLoopCou 1049 G4int internalLoopCounter = 0;
1068 do // while(Baryon[Di_q1-1][Di_q2-1][Pr 1050 do // while(Baryon[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]<>0);
1069 { 1051 {
1070 RightHadron=G4ParticleTable::GetParti 1052 RightHadron=G4ParticleTable::GetParticleTable()->FindParticle(SignDiQ*
1071 Baryon[Di_q1-1][Di_q2-1][Prod 1053 Baryon[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]);
1072 if (RightHadron == NULL) continue; 1054 if (RightHadron == NULL) continue;
1073 G4double RightHadronMass=RightHadron- 1055 G4double RightHadronMass=RightHadron->GetPDGMass();
1074 1056
1075 if (StringMass > LeftHadronMass + Rig 1057 if (StringMass > LeftHadronMass + RightHadronMass)
1076 { 1058 {
1077 if ( 1059 if ( NumberOf_FS > 349 ) {
1078 G4E 1060 G4ExceptionDescription ed;
1079 ed 1061 ed << " NumberOf_FS exceeds its limit: NumberOf_FS=" << NumberOf_FS << G4endl;
1080 G4E 1062 G4Exception( "G4LundStringFragmentation::Quark_Diquark_lastSplitting ",
1081 1063 "HAD_LUND_002", JustWarning, ed );
1082 Num 1064 NumberOf_FS = 349;
1083 } 1065 }
1084 1066
1085 G4double FS_Psqr=lambda(StringMassS 1067 G4double FS_Psqr=lambda(StringMassSqr,sqr(LeftHadronMass),
1086 sqr(RightHadronMass)); 1068 sqr(RightHadronMass));
1087 FS_Weight[NumberOf_FS]=std::sqrt(FS 1069 FS_Weight[NumberOf_FS]=std::sqrt(FS_Psqr)*
1088 MesonWeight[AbsIDquark-1 1070 MesonWeight[AbsIDquark-1][ProdQ-1][StateQ]*
1089 BaryonWeight[Di_q1-1][Di 1071 BaryonWeight[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]*
1090 Prob_QQbar[ProdQ-1]; 1072 Prob_QQbar[ProdQ-1];
1091 1073
1092 FS_LeftHadron[NumberOf_FS] = LeftHa 1074 FS_LeftHadron[NumberOf_FS] = LeftHadron;
1093 FS_RightHadron[NumberOf_FS]= RightH 1075 FS_RightHadron[NumberOf_FS]= RightHadron;
1094 1076
1095 NumberOf_FS++; 1077 NumberOf_FS++;
1096 } // End of if (StringMass > LeftHadr 1078 } // End of if (StringMass > LeftHadronMass + RightHadronMass)
1097 1079
1098 StateDiQ++; 1080 StateDiQ++;
1099 1081
1100 } while( (Baryon[Di_q1-1][Di_q2-1][Prod 1082 } while( (Baryon[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]!=0) &&
1101 ++internalLo 1083 ++internalLoopCounter < maxNumberOfInternalLoops );
1102 if ( internalLoopCoun 1084 if ( internalLoopCounter >= maxNumberOfInternalLoops ) {
1103 return false; 1085 return false;
1104 } 1086 }
1105 1087
1106 StateQ++; 1088 StateQ++;
1107 } while( (Meson[AbsIDquark-1][ProdQ-1][St 1089 } while( (Meson[AbsIDquark-1][ProdQ-1][StateQ]!=0) &&
1108 ++loopCounter < maxN 1090 ++loopCounter < maxNumberOfLoops ); /* Loop checking, 07.08.2015, A.Ribon */
1109 1091
1110 if ( loopCounter >= maxNumb 1092 if ( loopCounter >= maxNumberOfLoops ) {
1111 return false; 1093 return false;
1112 } 1094 }
1113 } 1095 }
1114 1096
1115 return true; 1097 return true;
1116 } 1098 }
1117 1099
1118 //------------------------------------------- 1100 //----------------------------------------------------------------------------------------
1119 1101
1120 G4bool G4LundStringFragmentation::Quark_AntiQ 1102 G4bool G4LundStringFragmentation::Quark_AntiQuark_lastSplitting(G4FragmentingString * & string,
1121 1103 G4ParticleDefinition * & LeftHadron,
1122 1104 G4ParticleDefinition * & RightHadron)
1123 { 1105 {
1124 G4double StringMass = string->Mass(); 1106 G4double StringMass = string->Mass();
1125 G4double StringMassSqr= sqr(StringMass); 1107 G4double StringMassSqr= sqr(StringMass);
1126 1108
1127 G4ParticleDefinition * Quark; 1109 G4ParticleDefinition * Quark;
1128 G4ParticleDefinition * Anti_Quark; 1110 G4ParticleDefinition * Anti_Quark;
1129 1111
1130 if (string->GetLeftParton()->GetPDGEncoding 1112 if (string->GetLeftParton()->GetPDGEncoding()>0)
1131 { 1113 {
1132 Quark =string->GetLeftParton(); 1114 Quark =string->GetLeftParton();
1133 Anti_Quark=string->GetRightParton(); 1115 Anti_Quark=string->GetRightParton();
1134 } else 1116 } else
1135 { 1117 {
1136 Quark =string->GetRightParton(); 1118 Quark =string->GetRightParton();
1137 Anti_Quark=string->GetLeftParton(); 1119 Anti_Quark=string->GetLeftParton();
1138 } 1120 }
1139 1121
1140 G4int IDquark =Quark->GetPDGEncodin 1122 G4int IDquark =Quark->GetPDGEncoding();
1141 G4int AbsIDquark =std::abs(IDquark); 1123 G4int AbsIDquark =std::abs(IDquark);
1142 G4int QuarkCharge =Qcharge[IDquar 1124 G4int QuarkCharge =Qcharge[IDquark-1];
1143 1125
1144 G4int IDanti_quark =Anti_Quark->GetPDGEn 1126 G4int IDanti_quark =Anti_Quark->GetPDGEncoding();
1145 G4int AbsIDanti_quark =std::abs(IDanti_quar 1127 G4int AbsIDanti_quark =std::abs(IDanti_quark);
1146 G4int AntiQuarkCharge =-Qcharge[AbsID 1128 G4int AntiQuarkCharge =-Qcharge[AbsIDanti_quark-1];
1147 1129
1148 G4int LeftHadronCharge(0), RightHadro 1130 G4int LeftHadronCharge(0), RightHadronCharge(0);
1149 1131
1150 //G4cout<<"Q Qbar "<<IDquark<<" "<<ID 1132 //G4cout<<"Q Qbar "<<IDquark<<" "<<IDanti_quark<<G4endl;
1151 1133
1152 NumberOf_FS=0; 1134 NumberOf_FS=0;
1153 for (G4int ProdQ=1; ProdQ < 4; ProdQ++) // 1135 for (G4int ProdQ=1; ProdQ < 4; ProdQ++) // Loop over quark-antiquark cases: u-ubar, d-dbar, s-sbar
1154 { // 1136 { // (as last splitting, do not consider c-cbar and b-bbar cases)
1155 //G4cout<<"NumberOf_FS ProdQ <<
1156 LeftHadronCharge = QuarkCharge - Qcharge[ 1137 LeftHadronCharge = QuarkCharge - Qcharge[ProdQ-1];
1157 G4int SignQ = LeftHadronCharge/3; if (Sig 1138 G4int SignQ = LeftHadronCharge/3; if (SignQ == 0) SignQ = 1;
1158 1139
1159 if ((IDquark == 1) && (ProdQ == 3)) SignQ 1140 if ((IDquark == 1) && (ProdQ == 3)) SignQ= 1; // K0 (d,sbar)
1160 if ((IDquark == 3) && (ProdQ == 1)) SignQ 1141 if ((IDquark == 3) && (ProdQ == 1)) SignQ=-1; // K0bar (s,dbar)
1161 if ((IDquark == 4) && (ProdQ 1142 if ((IDquark == 4) && (ProdQ == 2)) SignQ= 1; // D0 (c,ubar)
1162 if ((IDquark == 5) && (ProdQ 1143 if ((IDquark == 5) && (ProdQ == 1)) SignQ=-1; // anti_B0 (b,dbar)
1163 if ((IDquark == 5) && (ProdQ 1144 if ((IDquark == 5) && (ProdQ == 3)) SignQ=-1; // anti_Bs0 (b,sbar)
1164 1145
1165 RightHadronCharge = AntiQuark 1146 RightHadronCharge = AntiQuarkCharge + Qcharge[ProdQ-1];
1166 G4int SignAQ = RightHadronCharge/3; if (S 1147 G4int SignAQ = RightHadronCharge/3; if (SignAQ == 0) SignAQ = 1;
1167 1148
1168 if ((IDanti_quark ==-1) && (ProdQ == 3)) 1149 if ((IDanti_quark ==-1) && (ProdQ == 3)) SignAQ=-1; // K0bar (dbar,s)
1169 if ((IDanti_quark ==-3) && (ProdQ == 1)) 1150 if ((IDanti_quark ==-3) && (ProdQ == 1)) SignAQ= 1; // K0 (sbar,d)
1170 if ((IDanti_quark ==-4) && (ProdQ == 2)) 1151 if ((IDanti_quark ==-4) && (ProdQ == 2)) SignAQ=-1; // anti_D0 (cbar,u)
1171 if ((IDanti_quark ==-5) && (ProdQ == 1)) 1152 if ((IDanti_quark ==-5) && (ProdQ == 1)) SignAQ= 1; // B0 (bbar,d)
1172 if ((IDanti_quark ==-5) && (ProdQ == 3)) 1153 if ((IDanti_quark ==-5) && (ProdQ == 3)) SignAQ= 1; // Bs0 (bbar,s)
1173 1154
1174 //G4cout<<"ProQ signs "<<Prod 1155 //G4cout<<"ProQ signs "<<ProdQ<<" "<<SignQ<<" "<<SignAQ<<G4endl;
1175 1156
1176 G4int StateQ=0; 1157 G4int StateQ=0;
1177 const G4int maxNumberOfLoops 1158 const G4int maxNumberOfLoops = 1000;
1178 G4int loopCounter = 0; 1159 G4int loopCounter = 0;
1179 do 1160 do
1180 { 1161 {
1181 //G4cout<<"[AbsIDquar 1162 //G4cout<<"[AbsIDquark-1][ProdQ-1][StateQ "<<AbsIDquark-1<<" "
1182 //<<ProdQ-1<<" "<<Sta 1163 //<<ProdQ-1<<" "<<StateQ<<" "<<SignQ*Meson[AbsIDquark-1][ProdQ-1][StateQ]<<G4endl;
1183 LeftHadron=G4ParticleTable::GetParticle 1164 LeftHadron=G4ParticleTable::GetParticleTable()->FindParticle(SignQ*
1184 Meson[AbsIDquark-1][ProdQ- 1165 Meson[AbsIDquark-1][ProdQ-1][StateQ]);
1185 //G4cout<<"LeftHadron << 1166 //G4cout<<"LeftHadron "<<LeftHadron<<G4endl;
1186 if (LeftHadron == NULL) { StateQ++; con 1167 if (LeftHadron == NULL) { StateQ++; continue; }
1187 //G4cout<<"LeftHadron 1168 //G4cout<<"LeftHadron "<<LeftHadron->GetParticleName()<<G4endl;
1188 G4double LeftHadronMass=LeftHadron->Get 1169 G4double LeftHadronMass=LeftHadron->GetPDGMass();
1189 1170
1190 G4int StateAQ=0; 1171 G4int StateAQ=0;
1191 const G4int maxNumber 1172 const G4int maxNumberOfInternalLoops = 1000;
1192 G4int internalLoopCou 1173 G4int internalLoopCounter = 0;
1193 do 1174 do
1194 { 1175 {
1195 //G4cout<<" 1176 //G4cout<<" [AbsIDanti_quark-1][ProdQ-1][StateAQ] "<<AbsIDanti_quark-1<<" "
1196 // <<Pro 1177 // <<ProdQ-1<<" "<<StateAQ<<" "<<SignAQ*Meson[AbsIDanti_quark-1][ProdQ-1][StateAQ]<<G4endl;
1197 RightHadron=G4ParticleTable::GetParti 1178 RightHadron=G4ParticleTable::GetParticleTable()->FindParticle(SignAQ*
1198 Meson[AbsIDanti_quark-1][Prod 1179 Meson[AbsIDanti_quark-1][ProdQ-1][StateAQ]);
1199 //G4cout<<"Ri 1180 //G4cout<<"RightHadron "<<RightHadron<<G4endl;
1200 if(RightHadron == NULL) { StateAQ++; 1181 if(RightHadron == NULL) { StateAQ++; continue; }
1201 //G4cout<<"Ri 1182 //G4cout<<"RightHadron "<<RightHadron->GetParticleName()<<G4endl;
1202 G4double RightHadronMass=RightHadron- 1183 G4double RightHadronMass=RightHadron->GetPDGMass();
1203 1184
1204 if (StringMass > LeftHadronMass + Rig 1185 if (StringMass > LeftHadronMass + RightHadronMass)
1205 { 1186 {
1206 if ( 1187 if ( NumberOf_FS > 349 ) {
1207 G4E 1188 G4ExceptionDescription ed;
1208 ed 1189 ed << " NumberOf_FS exceeds its limit: NumberOf_FS=" << NumberOf_FS << G4endl;
1209 G4E 1190 G4Exception( "G4LundStringFragmentation::Quark_AntiQuark_lastSplitting ",
1210 1191 "HAD_LUND_003", JustWarning, ed );
1211 Num 1192 NumberOf_FS = 349;
1212 } 1193 }
1213 1194
1214 G4dou 1195 G4double FS_Psqr=lambda(StringMassSqr,sqr(LeftHadronMass),
1215 sqr(RightHadronMass)); 1196 sqr(RightHadronMass));
1216 //FS_Psqr=1.; 1197 //FS_Psqr=1.;
1217 FS_Weight[NumberOf_FS]=std::sqrt(FS 1198 FS_Weight[NumberOf_FS]=std::sqrt(FS_Psqr)*
1218 MesonWeight[AbsIDquark-1 1199 MesonWeight[AbsIDquark-1][ProdQ-1][StateQ]*
1219 MesonWeight[AbsIDanti_qu 1200 MesonWeight[AbsIDanti_quark-1][ProdQ-1][StateAQ]*
1220 Prob_QQbar[ProdQ-1]; 1201 Prob_QQbar[ProdQ-1];
>> 1202
1221 if (string->GetLeftParton()->GetPDG 1203 if (string->GetLeftParton()->GetPDGEncoding()>0)
1222 { 1204 {
1223 FS_LeftHadron[NumberOf_FS] = Righ 1205 FS_LeftHadron[NumberOf_FS] = RightHadron;
1224 FS_RightHadron[NumberOf_FS]= Left 1206 FS_RightHadron[NumberOf_FS]= LeftHadron;
1225 } else 1207 } else
1226 { 1208 {
1227 FS_LeftHadron[NumberOf_FS] = Left 1209 FS_LeftHadron[NumberOf_FS] = LeftHadron;
1228 FS_RightHadron[NumberOf_FS]= Righ 1210 FS_RightHadron[NumberOf_FS]= RightHadron;
1229 } 1211 }
1230 1212
1231 NumberOf_FS++; 1213 NumberOf_FS++;
1232 1214
1233 } 1215 }
1234 1216
1235 StateAQ++; 1217 StateAQ++;
1236 //G4cout<<" 1218 //G4cout<<" StateAQ Meson[AbsIDanti_quark-1][ProdQ-1][StateAQ] "<<StateAQ<<" "
1237 // <<Mes 1219 // <<Meson[AbsIDanti_quark-1][ProdQ-1][StateAQ]<<" "<<internalLoopCounter<<G4endl;
1238 } while ( (Meson[AbsIDanti_quark-1][Pro 1220 } while ( (Meson[AbsIDanti_quark-1][ProdQ-1][StateAQ]!=0) &&
1239 ++internalL 1221 ++internalLoopCounter < maxNumberOfInternalLoops );
1240 if ( internalLoopCo 1222 if ( internalLoopCounter >= maxNumberOfInternalLoops ) {
1241 return false; 1223 return false;
1242 } 1224 }
1243 1225
1244 StateQ++; 1226 StateQ++;
1245 //G4cout<<"StateQ Mes 1227 //G4cout<<"StateQ Meson[AbsIDquark-1][ProdQ-1][StateQ] "<<StateQ<<" "
1246 // <<Meson[AbsID 1228 // <<Meson[AbsIDquark-1][ProdQ-1][StateQ]<<" "<<loopCounter<<G4endl;
1247 1229
1248 } while ( (Meson[AbsIDquark-1][ProdQ-1][S 1230 } while ( (Meson[AbsIDquark-1][ProdQ-1][StateQ]!=0) &&
1249 ++loopCounter < maxN 1231 ++loopCounter < maxNumberOfLoops );
1250 if ( loopCounter >= maxNumb 1232 if ( loopCounter >= maxNumberOfLoops ) {
1251 return false; 1233 return false;
1252 } 1234 }
1253 } // End of for (G4int ProdQ=1; ProdQ < 4; 1235 } // End of for (G4int ProdQ=1; ProdQ < 4; ProdQ++)
1254 1236
1255 return true; 1237 return true;
1256 } 1238 }
1257 1239
1258 //------------------------------------------- 1240 //----------------------------------------------------------------------------------------------------------
1259 1241
1260 G4int G4LundStringFragmentation::SampleState( 1242 G4int G4LundStringFragmentation::SampleState(void)
1261 { 1243 {
1262 if ( NumberOf_FS > 349 ) { 1244 if ( NumberOf_FS > 349 ) {
1263 G4ExceptionDescription ed; 1245 G4ExceptionDescription ed;
1264 ed << " NumberOf_FS exceeds its lim 1246 ed << " NumberOf_FS exceeds its limit: NumberOf_FS=" << NumberOf_FS << G4endl;
1265 G4Exception( "G4LundStringFragmenta 1247 G4Exception( "G4LundStringFragmentation::SampleState ", "HAD_LUND_004", JustWarning, ed );
1266 NumberOf_FS = 349; 1248 NumberOf_FS = 349;
1267 } 1249 }
1268 1250
1269 G4double SumWeights=0.; 1251 G4double SumWeights=0.;
>> 1252
1270 for (G4int i=0; i<NumberOf_FS; i++) {SumWei 1253 for (G4int i=0; i<NumberOf_FS; i++) {SumWeights+=FS_Weight[i];}
1271 1254
1272 G4double ksi=G4UniformRand(); 1255 G4double ksi=G4UniformRand();
1273 G4double Sum=0.; 1256 G4double Sum=0.;
1274 G4int indexPosition = 0; 1257 G4int indexPosition = 0;
1275 1258
1276 for (G4int i=0; i<NumberOf_FS; i++) 1259 for (G4int i=0; i<NumberOf_FS; i++)
1277 { 1260 {
1278 Sum+=(FS_Weight[i]/SumWeights); 1261 Sum+=(FS_Weight[i]/SumWeights);
1279 indexPosition=i; 1262 indexPosition=i;
1280 if (Sum >= ksi) break; 1263 if (Sum >= ksi) break;
1281 } 1264 }
1282 return indexPosition; 1265 return indexPosition;
1283 } 1266 }
1284 1267
1285 //------------------------------------------- 1268 //----------------------------------------------------------------------------------------------------------
1286 1269
1287 void G4LundStringFragmentation::Sample4Moment 1270 void G4LundStringFragmentation::Sample4Momentum(G4LorentzVector* Mom, G4double Mass,
1288 1271 G4LorentzVector* AntiMom, G4double AntiMass,
1289 1272 G4double InitialMass)
1290 { 1273 {
1291 // ------ Sampling of momenta of 2 last pro 1274 // ------ Sampling of momenta of 2 last produced hadrons --------------------
1292 G4ThreeVector Pt; 1275 G4ThreeVector Pt;
1293 G4double MassMt, AntiMassMt; 1276 G4double MassMt, AntiMassMt;
1294 G4double AvailablePz, AvailablePz2; 1277 G4double AvailablePz, AvailablePz2;
1295 #ifdef debug_LUNDfragmentation 1278 #ifdef debug_LUNDfragmentation
1296 G4cout<<"Sampling of momenta of 2 las 1279 G4cout<<"Sampling of momenta of 2 last produced hadrons ----------------"<<G4endl;
1297 G4cout<<"Init Mass "<<InitialMass<<" 1280 G4cout<<"Init Mass "<<InitialMass<<" FirstM "<<Mass<<" SecondM "<<AntiMass<<" ProbIsotropy "<<G4endl;
1298 #endif 1281 #endif
1299 1282
1300 G4double r_val = sqr(InitialMass*InitialMas 1283 G4double r_val = sqr(InitialMass*InitialMass - Mass*Mass - AntiMass*AntiMass) -
1301 sqr(2.*Mass*AntiMass); 1284 sqr(2.*Mass*AntiMass);
1302 G4double Pabs = (r_val > 0.)? std::sqrt(r_v 1285 G4double Pabs = (r_val > 0.)? std::sqrt(r_val)/(2.*InitialMass) : 0;
1303 1286
1304 const G4int maxNumberOfLoops = 1000; 1287 const G4int maxNumberOfLoops = 1000;
1305 G4double SigmaQTw = SigmaQT; 1288 G4double SigmaQTw = SigmaQT;
1306 if ( Mass > 930. || AntiMass > 930. ) { 1289 if ( Mass > 930. || AntiMass > 930. ) {
1307 SigmaQT *= ( 1.0 - 0.55*sqr( (Mass+AntiMa 1290 SigmaQT *= ( 1.0 - 0.55*sqr( (Mass+AntiMass)/InitialMass ) );
1308 } 1291 }
1309 if ( Mass < 930. && AntiMass < 930. ) 1292 if ( Mass < 930. && AntiMass < 930. ) {} // q-qbar string
1310 if ( ( Mass < 930. && AntiMass > 930. 1293 if ( ( Mass < 930. && AntiMass > 930. ) ||
1311 ( Mass > 930. && AntiMass < 930. ) ) { 1294 ( Mass > 930. && AntiMass < 930. ) ) { // q-di_q string
1312 //SigmaQT = -1.; // isotropical de 1295 //SigmaQT = -1.; // isotropical decay
1313 } 1296 }
1314 if ( Mass > 930. && AntiMass > 930. ) 1297 if ( Mass > 930. && AntiMass > 930. ) { // qq-qqbar string
1315 SigmaQT *= ( 1.0 - 0.55*sqr( (Mass+ 1298 SigmaQT *= ( 1.0 - 0.55*sqr( (Mass+AntiMass)/InitialMass ) );
1316 } 1299 }
1317 1300
1318 G4int loopCounter = 0; 1301 G4int loopCounter = 0;
1319 do 1302 do
1320 { 1303 {
1321 Pt=SampleQuarkPt(Pabs); Pt.setZ(0); G4dou 1304 Pt=SampleQuarkPt(Pabs); Pt.setZ(0); G4double Pt2=Pt.mag2();
1322 MassMt = std::sqrt( Mass * Mass 1305 MassMt = std::sqrt( Mass * Mass + Pt2);
1323 AntiMassMt= std::sqrt(AntiMass * AntiMass 1306 AntiMassMt= std::sqrt(AntiMass * AntiMass + Pt2);
1324 } 1307 }
1325 while ( (InitialMass < MassMt + AntiMassMt) 1308 while ( (InitialMass < MassMt + AntiMassMt) && ++loopCounter < maxNumberOfLoops );
1326 1309
1327 SigmaQT = SigmaQTw; 1310 SigmaQT = SigmaQTw;
>> 1311
>> 1312 if ( loopCounter >= maxNumberOfLoops ) {
>> 1313 AvailablePz2 = 0.0;
>> 1314 }
1328 1315
1329 AvailablePz2= sqr(InitialMass*InitialMass - 1316 AvailablePz2= sqr(InitialMass*InitialMass - sqr(MassMt) - sqr(AntiMassMt)) -
1330 4.*sqr(MassMt*AntiMassMt); 1317 4.*sqr(MassMt*AntiMassMt);
1331 1318
1332 AvailablePz2 /=(4.*InitialMass*InitialMass) 1319 AvailablePz2 /=(4.*InitialMass*InitialMass);
1333 AvailablePz = std::sqrt(AvailablePz2); 1320 AvailablePz = std::sqrt(AvailablePz2);
1334 1321
1335 G4double Px=Pt.getX(); 1322 G4double Px=Pt.getX();
1336 G4double Py=Pt.getY(); 1323 G4double Py=Pt.getY();
1337 1324
1338 Mom->setPx(Px); Mom->setPy(Py); Mom->setPz( 1325 Mom->setPx(Px); Mom->setPy(Py); Mom->setPz(AvailablePz);
1339 Mom->setE(std::sqrt(sqr(MassMt)+AvailablePz 1326 Mom->setE(std::sqrt(sqr(MassMt)+AvailablePz2));
1340 1327
1341 AntiMom->setPx(-Px); AntiMom->setPy(-Py); A 1328 AntiMom->setPx(-Px); AntiMom->setPy(-Py); AntiMom->setPz(-AvailablePz);
1342 AntiMom->setE (std::sqrt(sqr(AntiMassMt)+Av 1329 AntiMom->setE (std::sqrt(sqr(AntiMassMt)+AvailablePz2));
1343 1330
1344 #ifdef debug_LUNDfragmentation 1331 #ifdef debug_LUNDfragmentation
1345 G4cout<<"Fmass Mom "<<Mom->getX()<<" 1332 G4cout<<"Fmass Mom "<<Mom->getX()<<" "<<Mom->getY()<<" "<<Mom->getZ()<<" "<<Mom->getT()<<G4endl;
1346 G4cout<<"Smass Mom "<<AntiMom->getX() 1333 G4cout<<"Smass Mom "<<AntiMom->getX()<<" "<<AntiMom->getY()<<" "<<AntiMom->getZ()
1347 <<" "<<AntiMom->getT()<<G4endl; 1334 <<" "<<AntiMom->getT()<<G4endl;
1348 #endif 1335 #endif
1349 } 1336 }
1350 1337
1351 //------------------------------------------- 1338 //------------------------------------------------------------------------
1352 1339
1353 G4double G4LundStringFragmentation::lambda(G4 1340 G4double G4LundStringFragmentation::lambda(G4double S, G4double m1_Sqr, G4double m2_Sqr)
1354 { 1341 {
1355 G4double lam = sqr(S - m1_Sqr - m2_Sqr) - 4 1342 G4double lam = sqr(S - m1_Sqr - m2_Sqr) - 4.*m1_Sqr*m2_Sqr;
1356 return lam; 1343 return lam;
1357 } 1344 }
1358 1345
1359 // ------------------------------------------ 1346 // --------------------------------------------------------------
1360 G4LundStringFragmentation::~G4LundStringFragm 1347 G4LundStringFragmentation::~G4LundStringFragmentation()
1361 {} 1348 {}
1362 1349
1363 1350