| 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 "G4QGSMFragmentation.hh" 33 #include "G4QGSMFragmentation.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 "G4ios.hh" 37 #include "G4ios.hh"
38 #include "G4FragmentingString.hh" 38 #include "G4FragmentingString.hh"
39 #include "G4DiQuarks.hh" 39 #include "G4DiQuarks.hh"
40 #include "G4Quarks.hh" 40 #include "G4Quarks.hh"
41 #include "G4HadronicParameters.hh" 41 #include "G4HadronicParameters.hh"
42 #include "G4Pow.hh" 42 #include "G4Pow.hh"
43 43
44 //#define debug_QGSMfragmentation 44 //#define debug_QGSMfragmentation
45 45
46 // Class G4QGSMFragmentation 46 // Class G4QGSMFragmentation
47 //******************************************** 47 //****************************************************************************************
48 48
49 G4QGSMFragmentation::G4QGSMFragmentation() 49 G4QGSMFragmentation::G4QGSMFragmentation()
50 { 50 {
51 SigmaQT = 0.45 * GeV; << 51 SigmaQT = 0.45 * GeV; // Uzhi June 2020
52 52
53 MassCut = 0.35*GeV; 53 MassCut = 0.35*GeV;
54 54
55 SetStrangenessSuppression((1.0 - 0.16)/2.) << 55 SetStrangenessSuppression((1.0 - 0.12)/2.); // Uzhi June 2020 0.16 -> 0.12
56 56
57 // Check if charmed and bottom hadrons are 57 // Check if charmed and bottom hadrons are enabled: if this is the case, then
58 // set the non-zero probabilities for c-cb 58 // set the non-zero probabilities for c-cbar and b-bbar creation from the vacuum,
59 // else set them to 0.0. If these probabil 59 // else set them to 0.0. If these probabilities are/aren't zero then charmed or bottom
60 // hadrons can't/can be created during the 60 // hadrons can't/can be created during the string fragmentation of ordinary
61 // (i.e. not heavy) projectile hadron nucl 61 // (i.e. not heavy) projectile hadron nuclear reactions.
62 if ( G4HadronicParameters::Instance()->Ena 62 if ( G4HadronicParameters::Instance()->EnableBCParticles() ) {
63 SetProbCCbar(0.0002); // According to O << 63 SetProbCCbar(0.005); // According to O.I. Piskunova Yad. Fiz. 56 (1993) 1094
64 SetProbBBbar(5.0e-5); // According to O 64 SetProbBBbar(5.0e-5); // According to O.I. Piskunova Yad. Fiz. 56 (1993) 1094
65 } else { 65 } else {
66 SetProbCCbar(0.0); 66 SetProbCCbar(0.0);
67 SetProbBBbar(0.0); 67 SetProbBBbar(0.0);
68 } 68 }
69 69
70 SetDiquarkSuppression(0.32); << 70 SetDiquarkSuppression(0.195); // Uzhi June 2020 0.32 -> 0.195
71 SetDiquarkBreakProbability(0.7); << 71 SetDiquarkBreakProbability(0.0); // Uzhi June 2020 0.7 -> 0.0
72 72
73 SetMinMasses(); 73 SetMinMasses();
74 74
75 arho = 0.5; // alpha_rho0 75 arho = 0.5; // alpha_rho0
76 aphi = 0.0; // alpha_fi 76 aphi = 0.0; // alpha_fi
77 aJPs =-2.2; // alpha_J/Psi 77 aJPs =-2.2; // alpha_J/Psi
78 aUps =-8.0; // alpha_Y ??? O. Pisk 78 aUps =-8.0; // alpha_Y ??? O. Piskunova Yad. Phys. 56 (1993) 1094.
79 79
80 aksi =-1.0; 80 aksi =-1.0;
81 alft = 0.5; // 2 * alpha'_R *<Pt^2> 81 alft = 0.5; // 2 * alpha'_R *<Pt^2>
82 82
83 an = -0.5 ; 83 an = -0.5 ;
84 ala = -0.75; // an - arho/2 + aphi/2 84 ala = -0.75; // an - arho/2 + aphi/2
85 alaC = an - arho/2.0 + aJPs/2.0; 85 alaC = an - arho/2.0 + aJPs/2.0;
86 alaB = an - arho/2.0 + aUps/2.0; 86 alaB = an - arho/2.0 + aUps/2.0;
87 aXi = 0.0; // ?? 87 aXi = 0.0; // ??
88 aXiC = 0.0; // ?? 88 aXiC = 0.0; // ??
89 aXiB = 0.0; // ?? 89 aXiB = 0.0; // ??
90 aXiCC = 0.0; // ?? 90 aXiCC = 0.0; // ??
91 aXiCB = 0.0; // ?? 91 aXiCB = 0.0; // ??
92 aXiBB = 0.0; // ?? 92 aXiBB = 0.0; // ??
93 93
94 SetFFq2q(); 94 SetFFq2q();
95 SetFFq2qq(); 95 SetFFq2qq();
96 SetFFqq2q(); 96 SetFFqq2q();
97 SetFFqq2qq(); 97 SetFFqq2qq();
98 // d u s c b 98 // d u s c b
99 G4int Index[5][5] = { { 0, 1, 2, 3, 4 } 99 G4int Index[5][5] = { { 0, 1, 2, 3, 4 }, // d
100 { 1, 5, 6, 7, 8 } 100 { 1, 5, 6, 7, 8 }, // u
101 { 2, 6, 9, 10, 11 } 101 { 2, 6, 9, 10, 11 }, // s
102 { 3, 7, 10, 12, 13 } 102 { 3, 7, 10, 12, 13 }, // c
103 { 4, 8, 11, 13, 14 } 103 { 4, 8, 11, 13, 14 } }; // b
104 for (G4int i = 0; i < 5; i++ ) { 104 for (G4int i = 0; i < 5; i++ ) {
105 for ( G4int j = 0; j < 5; j++ ) { 105 for ( G4int j = 0; j < 5; j++ ) {
106 IndexDiQ[i][j] = Index[i][j]; 106 IndexDiQ[i][j] = Index[i][j];
107 } 107 }
108 }; 108 };
109 } 109 }
110 110
111 G4QGSMFragmentation::~G4QGSMFragmentation() 111 G4QGSMFragmentation::~G4QGSMFragmentation()
112 {} 112 {}
113 113
114 //-------------------------------------------- 114 //----------------------------------------------------------------------------------------------------------
115 115
116 G4KineticTrackVector* G4QGSMFragmentation::Fra 116 G4KineticTrackVector* G4QGSMFragmentation::FragmentString(const G4ExcitedString& theString)
117 { 117 {
118 118
119 G4FragmentingString aString(theString); 119 G4FragmentingString aString(theString);
120 SetMinimalStringMass(&aString); 120 SetMinimalStringMass(&aString);
121 121
122 #ifdef debug_QGSMfragmentation 122 #ifdef debug_QGSMfragmentation
123 G4cout<<G4endl<<"QGSM StringFragm: String Ma 123 G4cout<<G4endl<<"QGSM StringFragm: String Mass "
124 <<theString.Get4M 124 <<theString.Get4Momentum().mag()<<" Pz "
125 <<theString.Get4M 125 <<theString.Get4Momentum().pz()
126 <<"-------------- 126 <<"------------------------------------"<<G4endl;
127 G4cout<<"String ends Direct "<<theString.Get 127 G4cout<<"String ends Direct "<<theString.GetLeftParton()->GetPDGcode()<<" "
128 <<theString.Get 128 <<theString.GetRightParton()->GetPDGcode()<<" "
129 <<theString.Get 129 <<theString.GetDirection()<< G4endl;
130 G4cout<<"Left mom "<<theString.GetLeftParto 130 G4cout<<"Left mom "<<theString.GetLeftParton()->Get4Momentum()<<G4endl;
131 G4cout<<"Right mom "<<theString.GetRightPart 131 G4cout<<"Right mom "<<theString.GetRightParton()->Get4Momentum()<<G4endl;
132 G4cout<<"Check for Fragmentation "<<G4endl; 132 G4cout<<"Check for Fragmentation "<<G4endl;
133 #endif 133 #endif
134 134
135 // Can no longer modify Parameters for Fragm 135 // Can no longer modify Parameters for Fragmentation.
136 PastInitPhase=true; 136 PastInitPhase=true;
137 137
138 // Check if string has enough mass to fragme 138 // Check if string has enough mass to fragment...
139 G4KineticTrackVector * LeftVector=NULL; 139 G4KineticTrackVector * LeftVector=NULL;
140 140
141 if ( !IsItFragmentable(&aString) ) { 141 if ( !IsItFragmentable(&aString) ) {
142 LeftVector=ProduceOneHadron(&theString); 142 LeftVector=ProduceOneHadron(&theString);
143 143
144 #ifdef debug_QGSMfragmentation 144 #ifdef debug_QGSMfragmentation
145 if ( LeftVector != 0 ) G4cout<<"Non fragm 145 if ( LeftVector != 0 ) G4cout<<"Non fragmentable - the string is converted to one hadron "<<G4endl;
146 #endif 146 #endif
147 147
148 if ( LeftVector == nullptr ) LeftVector = 148 if ( LeftVector == nullptr ) LeftVector = new G4KineticTrackVector;
149 return LeftVector; 149 return LeftVector;
150 } 150 }
151 151
152 #ifdef debug_QGSMfragmentation 152 #ifdef debug_QGSMfragmentation
153 G4cout<<"The string will be fragmented. "<<G 153 G4cout<<"The string will be fragmented. "<<G4endl;
154 #endif 154 #endif
155 155
156 LeftVector = new G4KineticTrackVector; 156 LeftVector = new G4KineticTrackVector;
157 G4KineticTrackVector * RightVector=new G4Kin 157 G4KineticTrackVector * RightVector=new G4KineticTrackVector;
158 158
159 G4ExcitedString *theStringInCMS=CopyExcited( 159 G4ExcitedString *theStringInCMS=CopyExcited(theString);
160 G4LorentzRotation toCms=theStringInCMS->Tran 160 G4LorentzRotation toCms=theStringInCMS->TransformToAlignedCms();
161 161
162 G4bool success=false, inner_sucess=true; 162 G4bool success=false, inner_sucess=true;
163 G4int attempt=0; 163 G4int attempt=0;
164 while ( !success && attempt++ < StringLoopIn 164 while ( !success && attempt++ < StringLoopInterrupt ) /* Loop checking, 07.08.2015, A.Ribon */
165 { 165 {
166 #ifdef debug_QGSMfragmentation 166 #ifdef debug_QGSMfragmentation
167 G4cout<<"Loop_toFrag "<<theStr 167 G4cout<<"Loop_toFrag "<<theStringInCMS->GetLeftParton()->GetPDGcode()<<" "
168 <<theStr 168 <<theStringInCMS->GetRightParton()->GetPDGcode()<<" "
169 <<theStr 169 <<theStringInCMS->GetDirection()<< G4endl;
170 #endif 170 #endif
171 171
172 G4FragmentingString *currentString=new G4F 172 G4FragmentingString *currentString=new G4FragmentingString(*theStringInCMS);
173 173
174 std::for_each(LeftVector->begin(), LeftVec 174 std::for_each(LeftVector->begin(), LeftVector->end(), DeleteKineticTrack());
175 LeftVector->clear(); 175 LeftVector->clear();
176 std::for_each(RightVector->begin(), RightV 176 std::for_each(RightVector->begin(), RightVector->end(), DeleteKineticTrack());
177 RightVector->clear(); 177 RightVector->clear();
178 178
179 inner_sucess=true; // set false on failur 179 inner_sucess=true; // set false on failure..
180 const G4int maxNumberOfLoops = 180 const G4int maxNumberOfLoops = 1000;
181 G4int loopCounter = -1; 181 G4int loopCounter = -1;
182 while (! StopFragmenting(currentString) && 182 while (! StopFragmenting(currentString) && ++loopCounter < maxNumberOfLoops ) /* Loop checking, 07.08.2015, A.Ribon */
183 { // Split current string into hadron + n 183 { // Split current string into hadron + new string
184 184
185 #ifdef debug_QGSMfragm 185 #ifdef debug_QGSMfragmentation
186 G4cout<<"The string ca 186 G4cout<<"The string can fragment. "<<G4endl;;
187 #endif 187 #endif
188 G4FragmentingString *newString=0; // us 188 G4FragmentingString *newString=0; // used as output from SplitUp...
189 G4KineticTrack * Hadron=Splitup(currentS 189 G4KineticTrack * Hadron=Splitup(currentString,newString);
190 190
191 if ( Hadron != 0 ) 191 if ( Hadron != 0 )
192 { 192 {
193 #ifdef debug_QGSMfr 193 #ifdef debug_QGSMfragmentation
194 G4cout<<"Hadron pro 194 G4cout<<"Hadron prod at fragm. "<<Hadron->GetDefinition()->GetParticleName()<<G4endl;
195 #endif 195 #endif
196 // To close the pro 196 // To close the production of hadrons at fragmentation stage
197 if ( currentString->GetDecayDirection 197 if ( currentString->GetDecayDirection() > 0 )
198 LeftVector->push_back(Hadron); 198 LeftVector->push_back(Hadron);
199 else 199 else
200 RightVector->push_back(Hadron); 200 RightVector->push_back(Hadron);
201 201
202 delete currentString; 202 delete currentString;
203 currentString=newString; 203 currentString=newString;
204 204
205 } else { 205 } else {
206 206
207 #ifdef debug_QGSMfr 207 #ifdef debug_QGSMfragmentation
208 G4cout<<"abandon .. 208 G4cout<<"abandon ... start from the beginning ---------------"<<G4endl;
209 #endif 209 #endif
210 210
211 // Abandon ... start from the beginni 211 // Abandon ... start from the beginning
212 if (newString) delete newString; 212 if (newString) delete newString;
213 inner_sucess=false; 213 inner_sucess=false;
214 break; 214 break;
215 } 215 }
216 } 216 }
217 if ( loopCounter >= maxNumberO 217 if ( loopCounter >= maxNumberOfLoops ) {
218 inner_sucess=false; 218 inner_sucess=false;
219 } 219 }
220 220
221 // Split current string into 2 final Hadro 221 // Split current string into 2 final Hadrons
222 #ifdef debug_QGSMfragmentation 222 #ifdef debug_QGSMfragmentation
223 if( inner_sucess ) { << 223 if( inner_sucess ) { // Uzhi June 2020
224 G4cout<<"Split remaining str 224 G4cout<<"Split remaining string into 2 final hadrons."<<G4endl;
225 } else { 225 } else {
226 G4cout<<" New attempt to fragment string 226 G4cout<<" New attempt to fragment string"<<G4endl;
227 } << 227 } // Uzhi June 2020
228 #endif 228 #endif
229 // To the close production of 229 // To the close production of hadrons at last string decay
230 if ( inner_sucess && 230 if ( inner_sucess &&
231 SplitLast(currentString,LeftVector, R 231 SplitLast(currentString,LeftVector, RightVector) )
232 { 232 {
233 success=true; 233 success=true;
234 } 234 }
235 delete currentString; 235 delete currentString;
236 } 236 }
237 237
238 delete theStringInCMS; 238 delete theStringInCMS;
239 239
240 if ( ! success ) 240 if ( ! success )
241 { 241 {
242 std::for_each(LeftVector->begin(), LeftVec 242 std::for_each(LeftVector->begin(), LeftVector->end(), DeleteKineticTrack());
243 LeftVector->clear(); 243 LeftVector->clear();
244 std::for_each(RightVector->begin(), RightV 244 std::for_each(RightVector->begin(), RightVector->end(), DeleteKineticTrack());
245 delete RightVector; 245 delete RightVector;
246 return LeftVector; 246 return LeftVector;
247 } 247 }
248 248
249 // Join Left- and RightVector into LeftVecto 249 // Join Left- and RightVector into LeftVector in correct order.
250 while(!RightVector->empty()) /* Loop checki 250 while(!RightVector->empty()) /* Loop checking, 07.08.2015, A.Ribon */
251 { 251 {
252 LeftVector->push_back(RightVector->back( 252 LeftVector->push_back(RightVector->back());
253 RightVector->erase(RightVector->end()-1) 253 RightVector->erase(RightVector->end()-1);
254 } 254 }
255 delete RightVector; 255 delete RightVector;
256 256
257 CalculateHadronTimePosition(theString.Get4Mo 257 CalculateHadronTimePosition(theString.Get4Momentum().mag(), LeftVector);
258 258
259 G4LorentzRotation toObserverFrame(toCms.inve 259 G4LorentzRotation toObserverFrame(toCms.inverse());
260 260
261 for (size_t C1 = 0; C1 < LeftVector->size(); 261 for (size_t C1 = 0; C1 < LeftVector->size(); C1++)
262 { 262 {
263 G4KineticTrack* Hadron = LeftVector->oper 263 G4KineticTrack* Hadron = LeftVector->operator[](C1);
264 G4LorentzVector Momentum = Hadron->Get4Mo 264 G4LorentzVector Momentum = Hadron->Get4Momentum();
265 Momentum = toObserverFrame*Momentum; 265 Momentum = toObserverFrame*Momentum;
266 Hadron->Set4Momentum(Momentum); 266 Hadron->Set4Momentum(Momentum);
267 G4LorentzVector Coordinate(Hadron->GetPos 267 G4LorentzVector Coordinate(Hadron->GetPosition(), Hadron->GetFormationTime());
268 Momentum = toObserverFrame*Coordinate; 268 Momentum = toObserverFrame*Coordinate;
269 Hadron->SetFormationTime(Momentum.e()); 269 Hadron->SetFormationTime(Momentum.e());
270 G4ThreeVector aPosition(Momentum.vect()); 270 G4ThreeVector aPosition(Momentum.vect());
271 Hadron->SetPosition(theString.GetPosition 271 Hadron->SetPosition(theString.GetPosition()+aPosition);
272 } 272 }
273 return LeftVector; 273 return LeftVector;
274 } 274 }
275 275
276 //-------------------------------------------- 276 //----------------------------------------------------------------------------------------------------------
277 277
278 G4bool G4QGSMFragmentation::IsItFragmentable(c 278 G4bool G4QGSMFragmentation::IsItFragmentable(const G4FragmentingString * string)
279 { 279 {
280 return sqr( MinimalStringMass + MassCut ) < << 280 //Uzhi June 2020 return sqr( PossibleHadronMass(string) + MassCut ) < string->Mass2();
>> 281 return sqr( MinimalStringMass + MassCut ) < string->Mass2(); // Uzhi June 2020
281 } 282 }
282 283
283 //-------------------------------------------- 284 //----------------------------------------------------------------------------------------------------------
284 285
285 G4bool G4QGSMFragmentation::StopFragmenting(co 286 G4bool G4QGSMFragmentation::StopFragmenting(const G4FragmentingString * string)
286 { 287 {
287 SetMinimalStringMass(string); 288 SetMinimalStringMass(string);
288 if ( MinimalStringMass < 0.0 ) return 289 if ( MinimalStringMass < 0.0 ) return true;
289 290
290 G4double smass = string->Mass(); 291 G4double smass = string->Mass();
291 G4double x = (string->IsAFourQuarkString()) 292 G4double x = (string->IsAFourQuarkString()) ? 0.005*(smass - MinimalStringMass)
292 : 0.66e-6*(smass - MinimalStringMass)*(sma 293 : 0.66e-6*(smass - MinimalStringMass)*(smass + MinimalStringMass);
293 294
294 G4bool res = true; 295 G4bool res = true;
295 if(x > 0.0) { 296 if(x > 0.0) {
296 res = (x < 200.) ? (G4UniformRand() 297 res = (x < 200.) ? (G4UniformRand() < G4Exp(-x)) : false;
297 } 298 }
298 return res; 299 return res;
299 } 300 }
300 301
301 //-------------------------------------------- 302 //-----------------------------------------------------------------------------
302 303
303 G4KineticTrack * G4QGSMFragmentation::Splitup( 304 G4KineticTrack * G4QGSMFragmentation::Splitup( G4FragmentingString *string,
304 G4FragmentingStri 305 G4FragmentingString *&newString )
305 { 306 {
306 #ifdef debug_QGSMfragmentation 307 #ifdef debug_QGSMfragmentation
307 G4cout<<G4endl; 308 G4cout<<G4endl;
308 G4cout<<"Start SplitUP (G4VLongitudinal 309 G4cout<<"Start SplitUP (G4VLongitudinalStringDecay) ========================="<<G4endl;
309 G4cout<<"String partons: " <<string->Ge 310 G4cout<<"String partons: " <<string->GetLeftParton()->GetPDGEncoding()<<" "
310 <<string->Ge 311 <<string->GetRightParton()->GetPDGEncoding()<<" "
311 <<"Direction " <<string->Ge 312 <<"Direction " <<string->GetDecayDirection()<<G4endl;
312 #endif 313 #endif
313 314
314 //... random choice of string end to us 315 //... random choice of string end to use for creating the hadron (decay)
315 G4int SideOfDecay = (G4UniformRand() < 316 G4int SideOfDecay = (G4UniformRand() < 0.5)? 1: -1;
316 if (SideOfDecay < 0) 317 if (SideOfDecay < 0)
317 { 318 {
318 string->SetLeftPartonStable(); 319 string->SetLeftPartonStable();
319 } else 320 } else
320 { 321 {
321 string->SetRightPartonStable(); 322 string->SetRightPartonStable();
322 } 323 }
323 324
324 G4ParticleDefinition *newStringEnd; 325 G4ParticleDefinition *newStringEnd;
325 G4ParticleDefinition * HadronDefinition 326 G4ParticleDefinition * HadronDefinition;
326 if (string->DecayIsQuark()) 327 if (string->DecayIsQuark())
327 { 328 {
328 G4double ProbDqADq = GetDiquarkSuppress(); 329 G4double ProbDqADq = GetDiquarkSuppress();
329 330
330 G4int NumberOfpossibleBaryons = 2; 331 G4int NumberOfpossibleBaryons = 2;
331 332
332 if (string->GetLeftParton()->GetParticleSu 333 if (string->GetLeftParton()->GetParticleSubType() != "quark") NumberOfpossibleBaryons++;
333 if (string->GetRightParton()->GetParticleS 334 if (string->GetRightParton()->GetParticleSubType() != "quark") NumberOfpossibleBaryons++;
334 335
335 G4double ActualProb = ProbDqADq ; 336 G4double ActualProb = ProbDqADq ;
336 ActualProb *= (1.0-G4Exp(2.0*(1.0 - 337 ActualProb *= (1.0-G4Exp(2.0*(1.0 - string->Mass()/(NumberOfpossibleBaryons*1400.0))));
337 338
338 SetDiquarkSuppression(ActualProb); 339 SetDiquarkSuppression(ActualProb);
339 340
340 HadronDefinition= QuarkSplitup(strin 341 HadronDefinition= QuarkSplitup(string->GetDecayParton(), newStringEnd);
341 342
342 SetDiquarkSuppression(ProbDqADq); 343 SetDiquarkSuppression(ProbDqADq);
343 } else { 344 } else {
344 HadronDefinition= DiQuarkSplitup(str 345 HadronDefinition= DiQuarkSplitup(string->GetDecayParton(), newStringEnd);
345 } 346 }
346 347
347 if ( HadronDefinition == NULL ) return 348 if ( HadronDefinition == NULL ) return NULL;
348 349
349 #ifdef debug_QGSMfragmentation 350 #ifdef debug_QGSMfragmentation
350 G4cout<<"The parton "<<string->GetDecay 351 G4cout<<"The parton "<<string->GetDecayParton()->GetPDGEncoding()<<" "
351 <<" produces hadron "<<HadronDefi 352 <<" produces hadron "<<HadronDefinition->GetParticleName()
352 <<" and is transformed to "<<newS 353 <<" and is transformed to "<<newStringEnd->GetPDGEncoding()<<G4endl;
353 G4cout<<"The side of the string decay L 354 G4cout<<"The side of the string decay Left/Right (1/-1) "<<SideOfDecay<<G4endl;
354 #endif 355 #endif
355 // create new String from old, ie. keep 356 // create new String from old, ie. keep Left and Right order, but replace decay
356 357
357 newString=new G4FragmentingString(*stri 358 newString=new G4FragmentingString(*string,newStringEnd); // To store possible
358 359 // quark containt of new string
359 360
360 #ifdef debug_QGSMfragmentation 361 #ifdef debug_QGSMfragmentation
361 G4cout<<"An attempt to determine its en 362 G4cout<<"An attempt to determine its energy (SplitEandP)"<<G4endl;
362 #endif 363 #endif
363 G4LorentzVector* HadronMomentum=SplitEa 364 G4LorentzVector* HadronMomentum=SplitEandP(HadronDefinition, string, newString);
364 365
365 delete newString; newString=0; 366 delete newString; newString=0;
366 367
367 G4KineticTrack * Hadron =0; 368 G4KineticTrack * Hadron =0;
368 if ( HadronMomentum != 0 ) { 369 if ( HadronMomentum != 0 ) {
369 370
370 #ifdef debug_QGSMfragmentation 371 #ifdef debug_QGSMfragmentation
371 G4cout<<"The attempt was successful 372 G4cout<<"The attempt was successful"<<G4endl;
372 #endif 373 #endif
373 G4ThreeVector Pos; 374 G4ThreeVector Pos;
374 Hadron = new G4KineticTrack(HadronDefinit 375 Hadron = new G4KineticTrack(HadronDefinition, 0,Pos, *HadronMomentum);
375 376
376 newString=new G4FragmentingString(*st 377 newString=new G4FragmentingString(*string,newStringEnd,HadronMomentum);
377 378
378 delete HadronMomentum; 379 delete HadronMomentum;
379 } 380 }
380 else 381 else
381 { 382 {
382 #ifdef debug_QGSMfragmentation 383 #ifdef debug_QGSMfragmentation
383 G4cout<<"The attempt was not successf 384 G4cout<<"The attempt was not successful !!!"<<G4endl;
384 #endif 385 #endif
385 } 386 }
386 387
387 #ifdef debug_VStringDecay 388 #ifdef debug_VStringDecay
388 G4cout<<"End SplitUP (G4VLongitudinalSt 389 G4cout<<"End SplitUP (G4VLongitudinalStringDecay) ====================="<<G4endl;
389 #endif 390 #endif
390 391
391 return Hadron; 392 return Hadron;
392 } 393 }
393 394
394 //-------------------------------------------- 395 //-----------------------------------------------------------------------------
395 396
396 G4ParticleDefinition *G4QGSMFragmentation::DiQ 397 G4ParticleDefinition *G4QGSMFragmentation::DiQuarkSplitup( G4ParticleDefinition* decay,
397 398 G4ParticleDefinition *&created )
398 { 399 {
399 //... can Diquark break or not? 400 //... can Diquark break or not?
400 if (G4UniformRand() < DiquarkBreakProb ) / 401 if (G4UniformRand() < DiquarkBreakProb ) //... Diquark break
401 { 402 {
402 G4int stableQuarkEncoding = decay->GetPD 403 G4int stableQuarkEncoding = decay->GetPDGEncoding()/1000;
403 G4int decayQuarkEncoding = (decay->GetPD 404 G4int decayQuarkEncoding = (decay->GetPDGEncoding()/100)%10;
404 405
405 if (G4UniformRand() < 0.5) 406 if (G4UniformRand() < 0.5)
406 { 407 {
407 G4int Swap = stableQuarkEncoding; 408 G4int Swap = stableQuarkEncoding;
408 stableQuarkEncoding = decayQuarkEncod 409 stableQuarkEncoding = decayQuarkEncoding;
409 decayQuarkEncoding = Swap; 410 decayQuarkEncoding = Swap;
410 } 411 }
411 412
412 G4int IsParticle=(decayQuarkEncoding>0) 413 G4int IsParticle=(decayQuarkEncoding>0) ? -1 : +1; // if we have a quark, we need antiquark
413 414
414 G4double StrSup=GetStrangeSuppress(); 415 G4double StrSup=GetStrangeSuppress();
415 SetStrangenessSuppression((1.0 - 0.07)/2 416 SetStrangenessSuppression((1.0 - 0.07)/2.); // Prob qq->K qq' 0.07
416 pDefPair QuarkPair = CreatePartonPair(Is 417 pDefPair QuarkPair = CreatePartonPair(IsParticle,false); // no diquarks wanted
417 SetStrangenessSuppression(StrSup); 418 SetStrangenessSuppression(StrSup);
418 419
419 //... Build new Diquark 420 //... Build new Diquark
420 G4int QuarkEncoding=QuarkPair.second->Ge 421 G4int QuarkEncoding=QuarkPair.second->GetPDGEncoding();
421 G4int i10 = std::max(std::abs(QuarkEnco 422 G4int i10 = std::max(std::abs(QuarkEncoding), std::abs(stableQuarkEncoding));
422 G4int i20 = std::min(std::abs(QuarkEnco 423 G4int i20 = std::min(std::abs(QuarkEncoding), std::abs(stableQuarkEncoding));
423 G4int spin = (i10 != i20 && G4UniformRan 424 G4int spin = (i10 != i20 && G4UniformRand() <= 0.5)? 1 : 3;
424 G4int NewDecayEncoding = -1*IsParticle*( 425 G4int NewDecayEncoding = -1*IsParticle*(i10 * 1000 + i20 * 100 + spin);
425 426
426 created = FindParticle(NewDecayEncoding) 427 created = FindParticle(NewDecayEncoding);
427 G4ParticleDefinition * decayQuark=FindPa 428 G4ParticleDefinition * decayQuark=FindParticle(decayQuarkEncoding);
428 G4ParticleDefinition * had=hadronizer->B 429 G4ParticleDefinition * had=hadronizer->Build(QuarkPair.first, decayQuark);
429 430
430 DecayQuark = decay->GetPDGEncoding(); << 431 DecayQuark = decay->GetPDGEncoding(); //Uzhi June 2020 decayQuarkEncoding;
431 NewQuark = NewDecayEncoding; << 432 NewQuark = NewDecayEncoding; //Uzhi June 2020 QuarkPair.first->GetPDGEncoding();
432 433
433 return had; 434 return had;
434 435
435 } else { //... Diquark does not break 436 } else { //... Diquark does not break
436 437
437 G4int IsParticle=(decay->GetPDGEncoding( 438 G4int IsParticle=(decay->GetPDGEncoding()>0) ? +1 : -1; // if we have a diquark, we need quark)
438 G4double StrSup=GetStrangeSuppress(); / 439 G4double StrSup=GetStrangeSuppress(); // for changing s-sbar production
439 SetStrangenessSuppression((1.0 - 0.07)/2 440 SetStrangenessSuppression((1.0 - 0.07)/2.);
440 pDefPair QuarkPair = CreatePartonPair(Is 441 pDefPair QuarkPair = CreatePartonPair(IsParticle,false); // no diquarks wanted
441 SetStrangenessSuppression(StrSup); 442 SetStrangenessSuppression(StrSup);
442 443
443 created = QuarkPair.second; 444 created = QuarkPair.second;
444 445
445 DecayQuark = decay->GetPDGEncoding(); 446 DecayQuark = decay->GetPDGEncoding();
446 NewQuark = created->GetPDGEncoding(); 447 NewQuark = created->GetPDGEncoding();
447 448
448 G4ParticleDefinition * had=hadronizer->B 449 G4ParticleDefinition * had=hadronizer->Build(QuarkPair.first, decay);
449 return had; 450 return had;
450 } 451 }
451 } 452 }
452 453
453 //-------------------------------------------- 454 //-----------------------------------------------------------------------------------------
454 455
455 G4LorentzVector * G4QGSMFragmentation::SplitEa 456 G4LorentzVector * G4QGSMFragmentation::SplitEandP(G4ParticleDefinition * pHadron,
456 457 G4FragmentingString * string,
457 458 G4FragmentingString * NewString)
458 { 459 {
459 G4double HadronMass = pHadron->GetPDGMa 460 G4double HadronMass = pHadron->GetPDGMass();
460 461
461 SetMinimalStringMass(NewString); 462 SetMinimalStringMass(NewString);
462 463
463 if ( MinimalStringMass < 0.0 ) return n 464 if ( MinimalStringMass < 0.0 ) return nullptr;
464 465
465 #ifdef debug_QGSMfragmentation 466 #ifdef debug_QGSMfragmentation
466 G4cout<<"G4QGSMFragmentation::SplitEand 467 G4cout<<"G4QGSMFragmentation::SplitEandP "<<pHadron->GetParticleName()<<G4endl;
467 G4cout<<"String 4 mom, String M "<<stri 468 G4cout<<"String 4 mom, String M "<<string->Get4Momentum()<<" "<<string->Mass()<<G4endl;
468 G4cout<<"HadM MinimalStringMassLeft Str 469 G4cout<<"HadM MinimalStringMassLeft StringM hM+sM "<<HadronMass<<" "<<MinimalStringMass<<" "
469 <<string->Mass()<<" "<<HadronMass 470 <<string->Mass()<<" "<<HadronMass+MinimalStringMass<<G4endl;
470 #endif 471 #endif
471 472
472 if (HadronMass + MinimalStringMass > st 473 if (HadronMass + MinimalStringMass > string->Mass())
473 { 474 {
474 #ifdef debug_QGSMfragmentation 475 #ifdef debug_QGSMfragmentation
475 G4cout<<"Mass of the string is not su 476 G4cout<<"Mass of the string is not sufficient to produce the hadron!"<<G4endl;
476 #endif 477 #endif
477 return 0; 478 return 0;
478 } // have to start all over! 479 } // have to start all over!
479 480
480 // calculate and assign hadron transver 481 // calculate and assign hadron transverse momentum component HadronPx andHadronPy
481 G4double StringMT2 = string->MassT2(); 482 G4double StringMT2 = string->MassT2();
482 G4double StringMT = std::sqrt(StringMT 483 G4double StringMT = std::sqrt(StringMT2);
483 484
484 G4LorentzVector String4Momentum = strin 485 G4LorentzVector String4Momentum = string->Get4Momentum();
485 String4Momentum.setPz(0.); 486 String4Momentum.setPz(0.);
486 G4ThreeVector StringPt = String4Momentu 487 G4ThreeVector StringPt = String4Momentum.vect();
487 488
488 G4ThreeVector HadronPt , RemSysPt; 489 G4ThreeVector HadronPt , RemSysPt;
489 G4double HadronMassT2, ResidualMas 490 G4double HadronMassT2, ResidualMassT2;
490 491
491 // Mt distribution is implemented << 492 //Uzhi June 2020 Mt distribution is implemented
492 G4double HadronMt, Pt, Pt2, phi; << 493 G4double HadronMt, Pt, Pt2, phi; // Uzhi June 2020
493 494
494 //... sample Pt of the hadron 495 //... sample Pt of the hadron
495 G4int attempt=0; 496 G4int attempt=0;
496 do 497 do
497 { 498 {
498 attempt++; if (attempt > StringLoopIn 499 attempt++; if (attempt > StringLoopInterrupt) return 0;
499 500
500 HadronMt = HadronMass - 200.0*G4Log(G << 501 HadronMt = HadronMass - 200.0*G4Log(G4UniformRand()); // Uzhi June 2020, 200.0 must be tuned
501 Pt2 = sqr(HadronMt)-sqr(HadronMass); << 502 Pt2 = sqr(HadronMt)-sqr(HadronMass); Pt=std::sqrt(Pt2); // Uzhi June 2020
502 phi = 2.*pi*G4UniformRand(); 503 phi = 2.*pi*G4UniformRand();
503 G4ThreeVector SampleQuarkPtw= G4Three 504 G4ThreeVector SampleQuarkPtw= G4ThreeVector(Pt*std::cos(phi), Pt*std::sin(phi), 0);
504 HadronPt =SampleQuarkPtw + string->D << 505 HadronPt =SampleQuarkPtw + string->DecayPt(); // Uzhi June 2020
>> 506 //Uzhi June 2020 HadronPt =SampleQuarkPt() + string->DecayPt(); // Save this for possible return
505 HadronPt.setZ(0); 507 HadronPt.setZ(0);
506 RemSysPt = StringPt - HadronPt; 508 RemSysPt = StringPt - HadronPt;
507 509
508 HadronMassT2 = sqr(HadronMass) + Hadr 510 HadronMassT2 = sqr(HadronMass) + HadronPt.mag2();
509 ResidualMassT2=sqr(MinimalStringMass) 511 ResidualMassT2=sqr(MinimalStringMass) + RemSysPt.mag2();
510 512
511 } while (std::sqrt(HadronMassT2) + std: 513 } while (std::sqrt(HadronMassT2) + std::sqrt(ResidualMassT2) > StringMT); /* Loop checking, 07.08.2015, A.Ribon */
512 514
513 //... sample z to define hadron longit 515 //... sample z to define hadron longitudinal momentum and energy
514 //... but first check the available pha 516 //... but first check the available phase space
515 517
516 G4double Pz2 = (sqr(StringMT2 - HadronM 518 G4double Pz2 = (sqr(StringMT2 - HadronMassT2 - ResidualMassT2) -
517 4*HadronMassT2 * ResidualMassT2)/4./ 519 4*HadronMassT2 * ResidualMassT2)/4./StringMT2;
518 520
519 if ( Pz2 < 0 ) {return 0;} // 521 if ( Pz2 < 0 ) {return 0;} // have to start all over!
520 522
521 //... then compute allowed z region z_ 523 //... then compute allowed z region z_min <= z <= z_max
522 524
523 G4double Pz = std::sqrt(Pz2); 525 G4double Pz = std::sqrt(Pz2);
524 G4double zMin = (std::sqrt(HadronMassT2 526 G4double zMin = (std::sqrt(HadronMassT2+Pz2) - Pz)/std::sqrt(StringMT2);
525 G4double zMax = (std::sqrt(HadronMassT2 527 G4double zMax = (std::sqrt(HadronMassT2+Pz2) + Pz)/std::sqrt(StringMT2);
526 528
527 if (zMin >= zMax) return 0; // have 529 if (zMin >= zMax) return 0; // have to start all over!
528 530
529 G4double z = GetLightConeZ(zMin, zMax, 531 G4double z = GetLightConeZ(zMin, zMax,
530 string->GetDecayParton() 532 string->GetDecayParton()->GetPDGEncoding(), pHadron,
531 HadronPt.x(), HadronPt.y 533 HadronPt.x(), HadronPt.y());
532 534
533 //... now compute hadron longitudinal m 535 //... now compute hadron longitudinal momentum and energy
534 // longitudinal hadron momentum compone 536 // longitudinal hadron momentum component HadronPz
535 537
536 HadronPt.setZ( 0.5* string->GetDecayDir 538 HadronPt.setZ( 0.5* string->GetDecayDirection() *
537 (z * string->LightConeDecay() - 539 (z * string->LightConeDecay() -
538 HadronMassT2/(z * string->LightCone 540 HadronMassT2/(z * string->LightConeDecay())) );
539 G4double HadronE = 0.5* (z * string->L 541 G4double HadronE = 0.5* (z * string->LightConeDecay() +
540 HadronMassT2/(z * string->LightConeDe 542 HadronMassT2/(z * string->LightConeDecay()) );
541 543
542 G4LorentzVector * a4Momentum= new G4Lor 544 G4LorentzVector * a4Momentum= new G4LorentzVector(HadronPt,HadronE);
543 545
544 #ifdef debug_QGSMfragmentation 546 #ifdef debug_QGSMfragmentation
545 G4cout<<"string->GetDecayDirection() st 547 G4cout<<"string->GetDecayDirection() string->LightConeDecay() "
546 <<string->GetDecayDirection()<<" 548 <<string->GetDecayDirection()<<" "<<string->LightConeDecay()<<G4endl;
547 G4cout<<"HadronPt,HadronE "<<HadronPt<< 549 G4cout<<"HadronPt,HadronE "<<HadronPt<<" "<<HadronE<<G4endl;
548 G4cout<<"Out of QGSM SplitEandP "<<G4en 550 G4cout<<"Out of QGSM SplitEandP "<<G4endl;
549 #endif 551 #endif
550 552
551 return a4Momentum; 553 return a4Momentum;
552 } 554 }
553 555
554 //-------------------------------------------- 556 //----------------------------------------------------------------------------------------------------------
555 557
556 G4double G4QGSMFragmentation::GetLightConeZ(G4 558 G4double G4QGSMFragmentation::GetLightConeZ(G4double zmin, G4double zmax, G4int /* PartonEncoding */ ,
557 G4 559 G4ParticleDefinition* /* pHadron */, G4double ptx , G4double pty)
558 { 560 {
559 G4double lambda = 2.0*(sqr(ptx)+sqr(pty))/sq << 561 G4double lambda = 2.0*(sqr(ptx)+sqr(pty))/sqr(GeV); // Uzhi June 2020
560 562
561 #ifdef debug_QGSMfragmentation 563 #ifdef debug_QGSMfragmentation
562 G4cout<<"GetLightConeZ zmin zmax Parton pHad 564 G4cout<<"GetLightConeZ zmin zmax Parton pHadron "<<zmin<<" "<<zmax<<" "/*<< PartonEncoding */
563 <<" "/*<< pHadron->GetParticleName() * 565 <<" "/*<< pHadron->GetParticleName() */ <<G4endl;
564 #endif 566 #endif
565 567
566 G4double z(0.); 568 G4double z(0.);
567 G4int DiQold(0), DiQnew(0); 569 G4int DiQold(0), DiQnew(0);
568 G4double d1(-1.0), d2(0.); 570 G4double d1(-1.0), d2(0.);
569 G4double invD1(0.),invD2(0.), r1(0.),r2(0.), 571 G4double invD1(0.),invD2(0.), r1(0.),r2(0.),r12(0.);
570 572
571 G4int absDecayQuarkCode = std::abs( DecayQua 573 G4int absDecayQuarkCode = std::abs( DecayQuark );
572 G4int absNewQuarkCode = std::abs( NewQuark 574 G4int absNewQuarkCode = std::abs( NewQuark );
573 575
574 G4int q1(0), q2(0); 576 G4int q1(0), q2(0);
575 // q1 = absDecayQuarkCode/1000; q2 = (absDe 577 // q1 = absDecayQuarkCode/1000; q2 = (absDecayQuarkCode % 1000)/100;
576 578
577 G4int qA(0), qB(0); 579 G4int qA(0), qB(0);
578 // qA = absNewQuarkCode/1000; qB = (absNe 580 // qA = absNewQuarkCode/1000; qB = (absNewQuarkCode % 1000)/100;
579 581
580 if ( (absDecayQuarkCode < 6) && (absNewQuark 582 if ( (absDecayQuarkCode < 6) && (absNewQuarkCode < 6) ) {
581 d1 = FFq2q[absDecayQuarkCode-1][absNewQuar 583 d1 = FFq2q[absDecayQuarkCode-1][absNewQuarkCode-1][0]; d2 = FFq2q[absDecayQuarkCode-1][absNewQuarkCode-1][1];
582 } 584 }
583 585
584 if ( (absDecayQuarkCode < 6) && (absNewQuark 586 if ( (absDecayQuarkCode < 6) && (absNewQuarkCode > 6) ) {
585 qA = absNewQuarkCode/1000; qB = (absNewQu 587 qA = absNewQuarkCode/1000; qB = (absNewQuarkCode % 1000)/100; DiQnew = IndexDiQ[qA-1][qB-1];
586 d1 = FFq2qq[absDecayQuarkCode-1][DiQnew][0] 588 d1 = FFq2qq[absDecayQuarkCode-1][DiQnew][0]; d2 = FFq2q[absDecayQuarkCode-1][DiQnew][1];
587 } 589 }
588 590
589 if ( (absDecayQuarkCode > 6) && (absNewQuark 591 if ( (absDecayQuarkCode > 6) && (absNewQuarkCode < 6) ) {
590 q1 = absDecayQuarkCode/1000; q2 = (absDeca 592 q1 = absDecayQuarkCode/1000; q2 = (absDecayQuarkCode % 1000)/100; DiQold = IndexDiQ[q1-1][q2-1];
591 d1 = FFqq2q[DiQold][absNewQuarkCode-1][0]; 593 d1 = FFqq2q[DiQold][absNewQuarkCode-1][0]; d2 = FFqq2q[DiQold][absNewQuarkCode-1][1];
592 } 594 }
593 595
594 if ( d1 < 0. ) { 596 if ( d1 < 0. ) {
595 q1 = absDecayQuarkCode/1000; q2 = (absDeca 597 q1 = absDecayQuarkCode/1000; q2 = (absDecayQuarkCode % 1000)/100; DiQold = IndexDiQ[q1-1][q2-1];
596 qA = absNewQuarkCode/1000; qB = (absNewQ << 598 d1 = FFqq2qq[DiQold][absNewQuarkCode-1][0]; d2 = FFqq2qq[DiQold][absNewQuarkCode-1][1];
597 d1 = FFqq2qq[DiQold][DiQnew][0]; d2 = FFqq <<
598 } 599 }
599 600
600 d2 +=lambda; << 601 d2 +=lambda; // Uzhi June 2020
601 d1+=1.0; d2+=1.0; 602 d1+=1.0; d2+=1.0;
602 603
603 invD1=1./d1; invD2=1./d2; 604 invD1=1./d1; invD2=1./d2;
604 605
605 const G4int maxNumberOfLoops = 10000; 606 const G4int maxNumberOfLoops = 10000;
606 G4int loopCounter = 0; 607 G4int loopCounter = 0;
607 do // Jong's algorithm 608 do // Jong's algorithm
608 { 609 {
609 r1=G4Pow::GetInstance()->powA(G4UniformRan 610 r1=G4Pow::GetInstance()->powA(G4UniformRand(),invD1);
610 r2=G4Pow::GetInstance()->powA(G4UniformRan 611 r2=G4Pow::GetInstance()->powA(G4UniformRand(),invD2);
611 r12=r1+r2; 612 r12=r1+r2;
612 z=r1/r12; 613 z=r1/r12;
613 } while ( ( (r12 > 1.0) || !((zmin <= z)&&(z 614 } while ( ( (r12 > 1.0) || !((zmin <= z)&&(z <= zmax))) &&
614 ++loopCounter < maxNumberOfLoops ) 615 ++loopCounter < maxNumberOfLoops ); /* Loop checking, 07.08.2015, A.Ribon */
615 616
616 if ( loopCounter >= maxNumberOfLoops ) { 617 if ( loopCounter >= maxNumberOfLoops ) {
617 z = 0.5*(zmin + zmax); // Just a value be 618 z = 0.5*(zmin + zmax); // Just a value between zmin and zmax, no physics considerations at all!
618 } 619 }
619 620
620 return z; 621 return z;
621 } 622 }
622 623
623 //-------------------------------------------- 624 //-----------------------------------------------------------------------------------------
624 625
625 G4bool G4QGSMFragmentation::SplitLast(G4Fragme 626 G4bool G4QGSMFragmentation::SplitLast(G4FragmentingString * string,
626 G4KineticTrackVector * Lef 627 G4KineticTrackVector * LeftVector,
627 G4KineticTrackVector * Right 628 G4KineticTrackVector * RightVector)
628 { 629 {
629 //... perform last cluster decay 630 //... perform last cluster decay
630 631
631 G4ThreeVector ClusterVel =string->Get4Mome 632 G4ThreeVector ClusterVel =string->Get4Momentum().boostVector();
632 G4double ResidualMass =string->Mass(); 633 G4double ResidualMass =string->Mass();
633 634
634 #ifdef debug_QGSMfragmentation 635 #ifdef debug_QGSMfragmentation
635 G4cout<<"Split last----------------------- 636 G4cout<<"Split last-----------------------------------------"<<G4endl;
636 G4cout<<"StrMass "<<ResidualMass<<" q's " 637 G4cout<<"StrMass "<<ResidualMass<<" q's "
637 <<string->GetLeftParton()->GetPartic 638 <<string->GetLeftParton()->GetParticleName()<<" "
638 <<string->GetRightParton()->GetParti 639 <<string->GetRightParton()->GetParticleName()<<G4endl;
639 #endif 640 #endif
640 641
641 G4int cClusterInterrupt = 0; 642 G4int cClusterInterrupt = 0;
642 G4ParticleDefinition *LeftHadron = nullptr 643 G4ParticleDefinition *LeftHadron = nullptr;
643 G4ParticleDefinition *RightHadron = nullpt 644 G4ParticleDefinition *RightHadron = nullptr;
644 const G4int maxNumberOfLoops = 1000; 645 const G4int maxNumberOfLoops = 1000;
645 G4int loopCounter = 0; 646 G4int loopCounter = 0;
646 647
647 G4double LeftHadronMass(0.); G4double Righ 648 G4double LeftHadronMass(0.); G4double RightHadronMass(0.);
648 do 649 do
649 { 650 {
650 if (cClusterInterrupt++ >= ClusterLoop << 651 if (cClusterInterrupt++ >= ClusterLoopInterrupt) return false; // Uzhi June 2020
651 LeftHadronMass = -MaxMass; RightHadron 652 LeftHadronMass = -MaxMass; RightHadronMass = -MaxMass;
652 653
653 G4ParticleDefinition * quark = nullptr; 654 G4ParticleDefinition * quark = nullptr;
654 string->SetLeftPartonStable(); // to query q 655 string->SetLeftPartonStable(); // to query quark contents..
655 656
656 if (string->DecayIsQuark() && string->Stable 657 if (string->DecayIsQuark() && string->StableIsQuark() )
657 { 658 {
658 //... there are quarks on cluster ends 659 //... there are quarks on cluster ends
659 660
660 G4int IsParticle=(string->GetLeftParton()- 661 G4int IsParticle=(string->GetLeftParton()->GetPDGEncoding()>0) ? -1 : +1;
661 // if we have a quark, we need 662 // if we have a quark, we need antiquark or diquark
662 663
663 pDefPair QuarkPair = CreatePartonPair(IsPa 664 pDefPair QuarkPair = CreatePartonPair(IsParticle);
664 quark = QuarkPair.second; 665 quark = QuarkPair.second;
665 666
666 LeftHadron= hadronizer->Build(QuarkPair.fi << 667 LeftHadron= hadronizer->BuildLowSpin(QuarkPair.first, string->GetLeftParton());
667 if ( LeftHadron == NULL ) continue; << 668 if ( LeftHadron == NULL ) continue; // Uzhi June 2020
668 RightHadron = hadronizer->Build(stri << 669 RightHadron = hadronizer->BuildLowSpin(string->GetRightParton(), quark); // Uzhi June 2020
669 if ( RightHadron == NULL ) continue; << 670 if ( RightHadron == NULL ) continue; // Uzhi June 2020
670 } else if( (!string->DecayIsQuark() && stri << 671 } else if( (!string->DecayIsQuark() && string->StableIsQuark() ) || // Uzhi June 2020
671 ( string->DecayIsQuark() && !stri << 672 ( string->DecayIsQuark() && !string->StableIsQuark() ) ) { // Uzhi June 2020
672 //... there is a Diquark on one of cluster 673 //... there is a Diquark on one of cluster ends
673 G4int IsParticle; 674 G4int IsParticle;
674 if ( string->StableIsQuark() ) { 675 if ( string->StableIsQuark() ) {
675 IsParticle=(string->GetLeftParton()->Get 676 IsParticle=(string->GetLeftParton()->GetPDGEncoding()>0) ? -1 : +1;
676 } else { 677 } else {
677 IsParticle=(string->GetLeftParton()->Get 678 IsParticle=(string->GetLeftParton()->GetPDGEncoding()>0) ? +1 : -1;
678 } 679 }
679 680
680 //G4double ProbSaS = 1.0 - 2.0 * G 681 //G4double ProbSaS = 1.0 - 2.0 * GetStrangeSuppress();
681 //G4double ActualProb = ProbSaS * 1. 682 //G4double ActualProb = ProbSaS * 1.4;
682 //SetStrangenessSuppression((1.0-Act 683 //SetStrangenessSuppression((1.0-ActualProb)/2.0);
683 684
684 pDefPair QuarkPair = CreatePartonPai 685 pDefPair QuarkPair = CreatePartonPair(IsParticle,false); // no diquarks wanted
685 //SetStrangenessSuppression((1.0-Pro 686 //SetStrangenessSuppression((1.0-ProbSaS)/2.0);
686 quark = QuarkPair.second; 687 quark = QuarkPair.second;
687 LeftHadron=hadronizer->Build(QuarkPa << 688 LeftHadron=hadronizer->BuildLowSpin(QuarkPair.first, string->GetLeftParton());
688 if ( LeftHadron == NULL ) continue; << 689 if ( LeftHadron == NULL ) continue; // Uzhi June 2020
689 RightHadron = hadronizer->Build(stri << 690 RightHadron = hadronizer->BuildLowSpin(string->GetRightParton(), quark); // Uzhi June 2020
690 if ( RightHadron == NULL ) continue; << 691 if ( RightHadron == NULL ) continue; // Uzhi June 2020
691 } else { // Diquark and anti-diquark << 692 } else { // Diquark and anti-diquark are on the string ends // Uzhi June 2020
>> 693 //+++++++++++++++++++++++++++++++ Inserted from FTF // Uzhi June 2020
>> 694 // Uzhi G4double StringMass = string->Mass();
692 if (cClusterInterrupt++ >= ClusterLo 695 if (cClusterInterrupt++ >= ClusterLoopInterrupt) return false;
693 G4int LeftQuark1= string->GetLeftPar 696 G4int LeftQuark1= string->GetLeftParton()->GetPDGEncoding()/1000;
694 G4int LeftQuark2=(string->GetLeftPar 697 G4int LeftQuark2=(string->GetLeftParton()->GetPDGEncoding()/100)%10;
695 G4int RightQuark1= string->GetRightP 698 G4int RightQuark1= string->GetRightParton()->GetPDGEncoding()/1000;
696 G4int RightQuark2=(string->GetRightP 699 G4int RightQuark2=(string->GetRightParton()->GetPDGEncoding()/100)%10;
697 if (G4UniformRand()<0.5) { 700 if (G4UniformRand()<0.5) {
698 LeftHadron =hadronizer->Build(Fin 701 LeftHadron =hadronizer->Build(FindParticle( LeftQuark1), FindParticle(RightQuark1));
699 RightHadron =hadronizer->Build(Fin 702 RightHadron =hadronizer->Build(FindParticle( LeftQuark2), FindParticle(RightQuark2));
700 } else { 703 } else {
701 LeftHadron =hadronizer->Build(Fin 704 LeftHadron =hadronizer->Build(FindParticle( LeftQuark1), FindParticle(RightQuark2));
702 RightHadron =hadronizer->Build(Fin 705 RightHadron =hadronizer->Build(FindParticle( LeftQuark2), FindParticle(RightQuark1));
703 } 706 }
704 if ( (LeftHadron == NULL) || (RightHadron 707 if ( (LeftHadron == NULL) || (RightHadron == NULL) ) continue;
>> 708 // End of inserting from FTF Uzhi June 2020
705 } 709 }
706 LeftHadronMass = LeftHadron->GetPDGMa 710 LeftHadronMass = LeftHadron->GetPDGMass();
707 RightHadronMass = RightHadron->GetPDGM 711 RightHadronMass = RightHadron->GetPDGMass();
708 //... repeat procedure, if mass of clu 712 //... repeat procedure, if mass of cluster is too low to produce hadrons
709 } while ( ( ResidualMass <= LeftHadronMass 713 } while ( ( ResidualMass <= LeftHadronMass + RightHadronMass )
710 && ++loopCounter < maxNumberOfLo 714 && ++loopCounter < maxNumberOfLoops ); /* Loop checking, 07.08.2015, A.Ribon */
711 715
712 if ( loopCounter >= maxNumberOfLoops ) { 716 if ( loopCounter >= maxNumberOfLoops ) {
713 return false; 717 return false;
714 } 718 }
715 719
716 //... compute hadron momenta and energies 720 //... compute hadron momenta and energies
717 G4LorentzVector LeftMom, RightMom; 721 G4LorentzVector LeftMom, RightMom;
718 G4ThreeVector Pos; 722 G4ThreeVector Pos;
719 Sample4Momentum(&LeftMom , LeftHadron->Get 723 Sample4Momentum(&LeftMom , LeftHadron->GetPDGMass() ,
720 &RightMom, RightHadron->Ge 724 &RightMom, RightHadron->GetPDGMass(), ResidualMass);
721 LeftMom.boost(ClusterVel); 725 LeftMom.boost(ClusterVel);
722 RightMom.boost(ClusterVel); 726 RightMom.boost(ClusterVel);
723 727
724 #ifdef debug_QGSMfragmentation 728 #ifdef debug_QGSMfragmentation
725 G4cout<<LeftHadron->GetParticleName()<<" " 729 G4cout<<LeftHadron->GetParticleName()<<" "<<RightHadron->GetParticleName()<<G4endl;
726 G4cout<<"Left Hadrom P M "<<LeftMom<<" "< 730 G4cout<<"Left Hadrom P M "<<LeftMom<<" "<<LeftMom.mag()<<G4endl;
727 G4cout<<"Right Hadrom P M "<<RightMom<<" " 731 G4cout<<"Right Hadrom P M "<<RightMom<<" "<<RightMom.mag()<<G4endl;
728 #endif 732 #endif
729 733
730 LeftVector->push_back(new G4KineticTrack(L 734 LeftVector->push_back(new G4KineticTrack(LeftHadron, 0, Pos, LeftMom));
731 RightVector->push_back(new G4KineticTrack( 735 RightVector->push_back(new G4KineticTrack(RightHadron, 0, Pos, RightMom));
732 736
733 return true; 737 return true;
734 } 738 }
735 739
736 //-------------------------------------------- 740 //----------------------------------------------------------------------------------------------------------
737 741
738 void G4QGSMFragmentation::Sample4Momentum(G4Lo 742 void G4QGSMFragmentation::Sample4Momentum(G4LorentzVector* Mom , G4double Mass ,
739 G4Lo 743 G4LorentzVector* AntiMom, G4double AntiMass, G4double InitialMass)
740 { 744 {
741 #ifdef debug_QGSMfragmentation << 745 #ifdef debug_QGSMfragmentation // Uzhi June 2020
742 G4cout<<"Sample4Momentum Last------------- 746 G4cout<<"Sample4Momentum Last-----------------------------------------"<<G4endl;
743 G4cout<<" StrMass "<<InitialMass<<" Mass1 747 G4cout<<" StrMass "<<InitialMass<<" Mass1 "<<Mass<<" Mass2 "<<AntiMass<<G4endl;
744 G4cout<<" SumMass "<<Mass+AntiMass<<G4end 748 G4cout<<" SumMass "<<Mass+AntiMass<<G4endl;
745 #endif 749 #endif
746 750
747 G4double r_val = sqr(InitialMass*InitialMa 751 G4double r_val = sqr(InitialMass*InitialMass - Mass*Mass - AntiMass*AntiMass) - sqr(2.*Mass*AntiMass);
748 G4double Pabs = (r_val > 0.)? std::sqrt(r_ 752 G4double Pabs = (r_val > 0.)? std::sqrt(r_val)/(2.*InitialMass) : 0;
749 753
750 //... sample unit vector 754 //... sample unit vector
751 G4double pz = 1. - 2.*G4UniformRand(); 755 G4double pz = 1. - 2.*G4UniformRand();
752 G4double st = std::sqrt(1. - pz * pz)* 756 G4double st = std::sqrt(1. - pz * pz)*Pabs;
753 G4double phi = 2.*pi*G4UniformRand(); 757 G4double phi = 2.*pi*G4UniformRand();
754 G4double px = st*std::cos(phi); 758 G4double px = st*std::cos(phi);
755 G4double py = st*std::sin(phi); 759 G4double py = st*std::sin(phi);
756 pz *= Pabs; 760 pz *= Pabs;
757 761
758 Mom->setPx(px); Mom->setPy(py); Mom->setPz 762 Mom->setPx(px); Mom->setPy(py); Mom->setPz(pz);
759 Mom->setE(std::sqrt(Pabs*Pabs + Mass*Mass) 763 Mom->setE(std::sqrt(Pabs*Pabs + Mass*Mass));
760 764
761 AntiMom->setPx(-px); AntiMom->setPy(-py); 765 AntiMom->setPx(-px); AntiMom->setPy(-py); AntiMom->setPz(-pz);
762 AntiMom->setE (std::sqrt(Pabs*Pabs + AntiM 766 AntiMom->setE (std::sqrt(Pabs*Pabs + AntiMass*AntiMass));
763 } 767 }
764 768
765 //-------------------------------------------- 769 //----------------------------------------------------------------------------------------------------------
766 770
767 void G4QGSMFragmentation::SetFFq2q() // q-> q 771 void G4QGSMFragmentation::SetFFq2q() // q-> q' + Meson (q anti q')
768 { 772 {
769 for (G4int i=0; i < 5; i++) { 773 for (G4int i=0; i < 5; i++) {
770 FFq2q[i][0][0] = 2.0 ; FFq2q[i][0][1] = -a 774 FFq2q[i][0][0] = 2.0 ; FFq2q[i][0][1] = -arho + alft; // q->d + (q dbar) Pi0, Eta, Eta', Rho0, omega
771 FFq2q[i][1][0] = 2.0 ; FFq2q[i][1][1] = -a 775 FFq2q[i][1][0] = 2.0 ; FFq2q[i][1][1] = -arho + alft; // q->u + (q ubar) Pi-, Rho-
772 FFq2q[i][2][0] = 1.0 ; FFq2q[i][2][1] = -a 776 FFq2q[i][2][0] = 1.0 ; FFq2q[i][2][1] = -aphi + alft; // q->s + (q sbar) K0, K*0
773 FFq2q[i][3][0] = 1.0 ; FFq2q[i][3][1] = -a 777 FFq2q[i][3][0] = 1.0 ; FFq2q[i][3][1] = -aJPs + alft; // q->c + (q+cbar) D-, D*-
774 FFq2q[i][4][0] = 1.0 ; FFq2q[i][4][1] = -a 778 FFq2q[i][4][0] = 1.0 ; FFq2q[i][4][1] = -aUps + alft; // q->b + (q bbar) EtaB, Upsilon
775 } 779 }
776 } 780 }
777 781
778 //-------------------------------------------- 782 //----------------------------------------------------------------------------------------------------------
779 783
780 void G4QGSMFragmentation::SetFFq2qq() // q-> 784 void G4QGSMFragmentation::SetFFq2qq() // q-> anti (q1'q2') + Baryon (q + q1 + q2)
781 { 785 {
782 for (G4int i=0; i < 5; i++) { 786 for (G4int i=0; i < 5; i++) {
783 FFq2qq[i][ 0][0] = 0.0 ; FFq2qq[i][ 0][1] 787 FFq2qq[i][ 0][0] = 0.0 ; FFq2qq[i][ 0][1] = arho - 2.0*an + alft; // q->dd bar + (q dd)
784 FFq2qq[i][ 1][0] = 0.0 ; FFq2qq[i][ 1][1] 788 FFq2qq[i][ 1][0] = 0.0 ; FFq2qq[i][ 1][1] = arho - 2.0*an + alft; // q->ud bar + (q ud)
785 FFq2qq[i][ 2][0] = 0.0 ; FFq2qq[i][ 2][1] 789 FFq2qq[i][ 2][0] = 0.0 ; FFq2qq[i][ 2][1] = arho - 2.0*ala + alft; // q->sd bar + (q sd)
786 FFq2qq[i][ 3][0] = 0.0 ; FFq2qq[i][ 3][1] 790 FFq2qq[i][ 3][0] = 0.0 ; FFq2qq[i][ 3][1] = arho - 2.0*alaC + alft; // q->cd bar + (q cd)
787 FFq2qq[i][ 4][0] = 0.0 ; FFq2qq[i][ 4][1] 791 FFq2qq[i][ 4][0] = 0.0 ; FFq2qq[i][ 4][1] = arho - 2.0*alaB + alft; // q->bd bar + (q bd)
788 FFq2qq[i][ 5][0] = 0.0 ; FFq2qq[i][ 5][1] 792 FFq2qq[i][ 5][0] = 0.0 ; FFq2qq[i][ 5][1] = arho - 2.0*an + alft; // q->uu bar + (q uu)
789 FFq2qq[i][ 6][0] = 0.0 ; FFq2qq[i][ 6][1] 793 FFq2qq[i][ 6][0] = 0.0 ; FFq2qq[i][ 6][1] = arho - 2.0*ala + alft; // q->su bar + (q su)
790 FFq2qq[i][ 7][0] = 0.0 ; FFq2qq[i][ 7][1] 794 FFq2qq[i][ 7][0] = 0.0 ; FFq2qq[i][ 7][1] = arho - 2.0*alaC + alft; // q->cu bar + (q cu)
791 FFq2qq[i][ 8][0] = 0.0 ; FFq2qq[i][ 8][1] 795 FFq2qq[i][ 8][0] = 0.0 ; FFq2qq[i][ 8][1] = arho - 2.0*alaB + alft; // q->bu bar + (q bu)
792 FFq2qq[i][ 9][0] = 0.0 ; FFq2qq[i][ 9][1] 796 FFq2qq[i][ 9][0] = 0.0 ; FFq2qq[i][ 9][1] = arho - 2.0*aXi + alft; // q->ss bar + (q ss)
793 FFq2qq[i][10][0] = 0.0 ; FFq2qq[i][10][1] 797 FFq2qq[i][10][0] = 0.0 ; FFq2qq[i][10][1] = arho - 2.0*aXiC + alft; // q->cs bar + (q cs)
794 FFq2qq[i][11][0] = 0.0 ; FFq2qq[i][11][1] 798 FFq2qq[i][11][0] = 0.0 ; FFq2qq[i][11][1] = arho - 2.0*aXiB + alft; // q->bs bar + (q bc)
795 FFq2qq[i][12][0] = 0.0 ; FFq2qq[i][12][1] 799 FFq2qq[i][12][0] = 0.0 ; FFq2qq[i][12][1] = arho - 2.0*aXiCC + alft; // q->cc bar + (q cc)
796 FFq2qq[i][13][0] = 0.0 ; FFq2qq[i][13][1] 800 FFq2qq[i][13][0] = 0.0 ; FFq2qq[i][13][1] = arho - 2.0*aXiCB + alft; // q->bc bar + (q bc)
797 FFq2qq[i][14][0] = 0.0 ; FFq2qq[i][14][1] 801 FFq2qq[i][14][0] = 0.0 ; FFq2qq[i][14][1] = arho - 2.0*aXiBB + alft; // q->bb bar + (q bb)
798 } 802 }
799 } 803 }
800 804
801 //-------------------------------------------- 805 //----------------------------------------------------------------------------------------------------------
802 806
803 void G4QGSMFragmentation::SetFFqq2q() // q1q2 807 void G4QGSMFragmentation::SetFFqq2q() // q1q2-> anti(q') + Baryon (q1 + q2 + q')
804 { 808 {
805 for (G4int i=0; i < 15; i++) { 809 for (G4int i=0; i < 15; i++) {
806 FFqq2q[i][0][0] = 2.0*(arho - an); FFqq2q[ 810 FFqq2q[i][0][0] = 2.0*(arho - an); FFqq2q[i][0][1] = -arho + alft;
807 FFqq2q[i][1][0] = 2.0*(arho - an); FFqq2q[ 811 FFqq2q[i][1][0] = 2.0*(arho - an); FFqq2q[i][1][1] = -arho + alft;
808 FFqq2q[i][2][0] = 2.0*(arho - an); FFqq2q[ 812 FFqq2q[i][2][0] = 2.0*(arho - an); FFqq2q[i][2][1] = -aphi + alft;
809 FFqq2q[i][3][0] = 2.0*(arho - an); FFqq2q[ 813 FFqq2q[i][3][0] = 2.0*(arho - an); FFqq2q[i][3][1] = -aJPs + alft;
810 FFqq2q[i][4][0] = 2.0*(arho - an); FFqq2q[ 814 FFqq2q[i][4][0] = 2.0*(arho - an); FFqq2q[i][4][1] = -aUps + alft;
811 } 815 }
812 } 816 }
813 817
814 //-------------------------------------------- 818 //----------------------------------------------------------------------------------------------------------
815 819
816 void G4QGSMFragmentation::SetFFqq2qq() // q1( 820 void G4QGSMFragmentation::SetFFqq2qq() // q1(q2)-> q'(q2) + Meson(q1 anti q')
817 { 821 {
818 for (G4int i=0; i < 15; i++) { 822 for (G4int i=0; i < 15; i++) {
819 FFqq2qq[i][0][0] = 0. ; FFqq2qq[i][0][1] 823 FFqq2qq[i][0][0] = 0. ; FFqq2qq[i][0][1] = 2.0*arho - 2.0*an -arho + alft; // dd -> dd + Pi0 (d d bar)
820 FFqq2qq[i][1][0] = 0. ; FFqq2qq[i][1][1] 824 FFqq2qq[i][1][0] = 0. ; FFqq2qq[i][1][1] = 2.0*arho - 2.0*an -arho + alft; // dd -> ud + Pi- (d u bar)
821 FFqq2qq[i][2][0] = 0. ; FFqq2qq[i][2][1] 825 FFqq2qq[i][2][0] = 0. ; FFqq2qq[i][2][1] = 2.0*arho - 2.0*an -aphi + alft; // dd -> sd + K0 (d s bar)
822 FFqq2qq[i][3][0] = 0. ; FFqq2qq[i][3][1] 826 FFqq2qq[i][3][0] = 0. ; FFqq2qq[i][3][1] = 2.0*arho - 2.0*an -aJPs + alft; // dd -> cd + D- (d c bar)
823 FFqq2qq[i][4][0] = 0. ; FFqq2qq[i][4][1] 827 FFqq2qq[i][4][0] = 0. ; FFqq2qq[i][4][1] = 2.0*arho - 2.0*an -aUps + alft; // dd -> bd + B0 (d b bar)
824 } 828 }
825 } 829 }
826 830
827 831