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