| 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,v 1.6 2007/04/24 14:55:23 gunter Exp $
>> 28 // GEANT4 tag $Name: geant4-09-00 $
27 // 29 //
28 // ------------------------------------------- 30 // -----------------------------------------------------------------------------
29 // GEANT 4 class implementation file 31 // GEANT 4 class implementation file
30 // 32 //
31 // History: first implementation, Maxim K 33 // History: first implementation, Maxim Komogorov, 10-Jul-1998
32 // ------------------------------------------- 34 // -----------------------------------------------------------------------------
33 #include "G4QGSMFragmentation.hh" 35 #include "G4QGSMFragmentation.hh"
34 #include "G4PhysicalConstants.hh" <<
35 #include "G4SystemOfUnits.hh" <<
36 #include "Randomize.hh" <<
37 #include "G4ios.hh" <<
38 #include "G4FragmentingString.hh" 36 #include "G4FragmentingString.hh"
39 #include "G4DiQuarks.hh" 37 #include "G4DiQuarks.hh"
40 #include "G4Quarks.hh" 38 #include "G4Quarks.hh"
41 #include "G4HadronicParameters.hh" <<
42 #include "G4Pow.hh" <<
43 39
44 //#define debug_QGSMfragmentation << 40 #include "Randomize.hh"
>> 41 #include "G4ios.hh"
45 42
46 // Class G4QGSMFragmentation 43 // Class G4QGSMFragmentation
47 //******************************************** 44 //****************************************************************************************
48 45
49 G4QGSMFragmentation::G4QGSMFragmentation() << 46 G4QGSMFragmentation::G4QGSMFragmentation() :
50 { << 47 arho(0.5), aphi(0.), an(-0.5), ala(-0.75), aksi(-1.), alft(0.5)
51 SigmaQT = 0.45 * GeV; << 48 {
52 << 49 }
53 MassCut = 0.35*GeV; <<
54 50
55 SetStrangenessSuppression((1.0 - 0.16)/2.) << 51 G4QGSMFragmentation::G4QGSMFragmentation(const G4QGSMFragmentation &) : G4VLongitudinalStringDecay(),
>> 52 arho(0.5), aphi(0.), an(-0.5), ala(-0.75), aksi(-1.), alft(0.5)
>> 53 {
>> 54 }
56 55
57 // Check if charmed and bottom hadrons are << 56 G4QGSMFragmentation::~G4QGSMFragmentation()
58 // set the non-zero probabilities for c-cb << 57 {
59 // else set them to 0.0. If these probabil << 58 }
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 59
73 SetMinMasses(); << 60 //****************************************************************************************
74 61
75 arho = 0.5; // alpha_rho0 << 62 const G4QGSMFragmentation & G4QGSMFragmentation::operator=(const G4QGSMFragmentation &)
76 aphi = 0.0; // alpha_fi << 63 {
77 aJPs =-2.2; // alpha_J/Psi << 64 throw G4HadronicException(__FILE__, __LINE__, "G4QGSMFragmentation::operator= meant to not be accessable");
78 aUps =-8.0; // alpha_Y ??? O. Pisk << 65 return *this;
79 << 66 }
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 } <<
110 67
111 G4QGSMFragmentation::~G4QGSMFragmentation() << 68 int G4QGSMFragmentation::operator==(const G4QGSMFragmentation &right) const
112 {} << 69 {
>> 70 return !memcmp(this, &right, sizeof(G4QGSMFragmentation));
>> 71 }
113 72
>> 73 int G4QGSMFragmentation::operator!=(const G4QGSMFragmentation &right) const
>> 74 {
>> 75 return memcmp(this, &right, sizeof(G4QGSMFragmentation));
>> 76 }
>> 77
>> 78 //****************************************************************************************
114 //-------------------------------------------- 79 //----------------------------------------------------------------------------------------------------------
115 80
116 G4KineticTrackVector* G4QGSMFragmentation::Fra 81 G4KineticTrackVector* G4QGSMFragmentation::FragmentString(const G4ExcitedString& theString)
117 { 82 {
118 << 83 // Can no longer modify Parameters for Fragmentation.
119 G4FragmentingString aString(theString); << 84 PastInitPhase=true;
120 SetMinimalStringMass(&aString); <<
121 <<
122 #ifdef debug_QGSMfragmentation <<
123 G4cout<<G4endl<<"QGSM StringFragm: String Ma <<
124 <<theString.Get4M <<
125 <<theString.Get4M <<
126 <<"-------------- <<
127 G4cout<<"String ends Direct "<<theString.Get <<
128 <<theString.Get <<
129 <<theString.Get <<
130 G4cout<<"Left mom "<<theString.GetLeftParto <<
131 G4cout<<"Right mom "<<theString.GetRightPart <<
132 G4cout<<"Check for Fragmentation "<<G4endl; <<
133 #endif <<
134 <<
135 // Can no longer modify Parameters for Fragm <<
136 PastInitPhase=true; <<
137 85
138 // Check if string has enough mass to fragme << 86 // check if string has enough mass to fragment...
139 G4KineticTrackVector * LeftVector=NULL; << 87 G4KineticTrackVector * LeftVector=LightFragmentationTest(&theString);
140 << 88 if ( LeftVector != 0 ) return LeftVector;
141 if ( !IsItFragmentable(&aString) ) { <<
142 LeftVector=ProduceOneHadron(&theString); <<
143 <<
144 #ifdef debug_QGSMfragmentation <<
145 if ( LeftVector != 0 ) G4cout<<"Non fragm <<
146 #endif <<
147 <<
148 if ( LeftVector == nullptr ) LeftVector = <<
149 return LeftVector; <<
150 } <<
151 <<
152 #ifdef debug_QGSMfragmentation <<
153 G4cout<<"The string will be fragmented. "<<G <<
154 #endif <<
155 89
156 LeftVector = new G4KineticTrackVector; << 90 LeftVector = new G4KineticTrackVector;
157 G4KineticTrackVector * RightVector=new G4Kin << 91 G4KineticTrackVector * RightVector=new G4KineticTrackVector;
158 <<
159 G4ExcitedString *theStringInCMS=CopyExcited( <<
160 G4LorentzRotation toCms=theStringInCMS->Tran <<
161 <<
162 G4bool success=false, inner_sucess=true; <<
163 G4int attempt=0; <<
164 while ( !success && attempt++ < StringLoopIn <<
165 { <<
166 #ifdef debug_QGSMfragmentation <<
167 G4cout<<"Loop_toFrag "<<theStr <<
168 <<theStr <<
169 <<theStr <<
170 #endif <<
171 92
>> 93 // this should work but its only a semi deep copy. %GF G4ExcitedString theStringInCMS(theString);
>> 94 G4ExcitedString *theStringInCMS=CPExcited(theString);
>> 95 G4LorentzRotation toCms=theStringInCMS->TransformToAlignedCms();
>> 96
>> 97 G4bool success=false, inner_sucess=true;
>> 98 G4int attempt=0;
>> 99 while ( !success && attempt++ < StringLoopInterrupt )
>> 100 {
172 G4FragmentingString *currentString=new G4F 101 G4FragmentingString *currentString=new G4FragmentingString(*theStringInCMS);
173 102
174 std::for_each(LeftVector->begin(), LeftVec 103 std::for_each(LeftVector->begin(), LeftVector->end(), DeleteKineticTrack());
175 LeftVector->clear(); 104 LeftVector->clear();
176 std::for_each(RightVector->begin(), RightV 105 std::for_each(RightVector->begin(), RightVector->end(), DeleteKineticTrack());
177 RightVector->clear(); 106 RightVector->clear();
178 107
179 inner_sucess=true; // set false on failur 108 inner_sucess=true; // set false on failure..
180 const G4int maxNumberOfLoops = << 109 while (! StopFragmenting(currentString) )
181 G4int loopCounter = -1; <<
182 while (! StopFragmenting(currentString) && <<
183 { // Split current string into hadron + n 110 { // Split current string into hadron + new string
184 <<
185 #ifdef debug_QGSMfragm <<
186 G4cout<<"The string ca <<
187 #endif <<
188 G4FragmentingString *newString=0; // us 111 G4FragmentingString *newString=0; // used as output from SplitUp...
189 G4KineticTrack * Hadron=Splitup(currentS 112 G4KineticTrack * Hadron=Splitup(currentString,newString);
190 << 113 if ( Hadron != 0 && IsFragmentable(newString))
191 if ( Hadron != 0 ) <<
192 { 114 {
193 #ifdef debug_QGSMfr <<
194 G4cout<<"Hadron pro <<
195 #endif <<
196 // To close the pro <<
197 if ( currentString->GetDecayDirection 115 if ( currentString->GetDecayDirection() > 0 )
198 LeftVector->push_back(Hadron); 116 LeftVector->push_back(Hadron);
199 else 117 else
200 RightVector->push_back(Hadron); 118 RightVector->push_back(Hadron);
201 <<
202 delete currentString; 119 delete currentString;
203 currentString=newString; 120 currentString=newString;
204 <<
205 } else { 121 } else {
206 << 122 // abandon ... start from the beginning
207 #ifdef debug_QGSMfr <<
208 G4cout<<"abandon .. <<
209 #endif <<
210 <<
211 // Abandon ... start from the beginni <<
212 if (newString) delete newString; 123 if (newString) delete newString;
>> 124 if (Hadron) delete Hadron;
213 inner_sucess=false; 125 inner_sucess=false;
214 break; 126 break;
215 } 127 }
216 } << 128 }
217 if ( loopCounter >= maxNumberO <<
218 inner_sucess=false; <<
219 } <<
220 <<
221 // Split current string into 2 final Hadro 129 // Split current string into 2 final Hadrons
222 #ifdef debug_QGSMfragmentation <<
223 if( inner_sucess ) { <<
224 G4cout<<"Split remaining str <<
225 } else { <<
226 G4cout<<" New attempt to fragment string <<
227 } <<
228 #endif <<
229 // To the close production of <<
230 if ( inner_sucess && 130 if ( inner_sucess &&
231 SplitLast(currentString,LeftVector, R 131 SplitLast(currentString,LeftVector, RightVector) )
232 { 132 {
233 success=true; 133 success=true;
234 } 134 }
235 delete currentString; 135 delete currentString;
236 } 136 }
237 137
238 delete theStringInCMS; 138 delete theStringInCMS;
239 139
240 if ( ! success ) 140 if ( ! success )
241 { 141 {
242 std::for_each(LeftVector->begin(), LeftVec 142 std::for_each(LeftVector->begin(), LeftVector->end(), DeleteKineticTrack());
243 LeftVector->clear(); 143 LeftVector->clear();
244 std::for_each(RightVector->begin(), RightV 144 std::for_each(RightVector->begin(), RightVector->end(), DeleteKineticTrack());
245 delete RightVector; 145 delete RightVector;
246 return LeftVector; 146 return LeftVector;
247 } 147 }
248 148
249 // Join Left- and RightVector into LeftVecto 149 // Join Left- and RightVector into LeftVector in correct order.
250 while(!RightVector->empty()) /* Loop checki << 150 while(!RightVector->empty())
251 { 151 {
252 LeftVector->push_back(RightVector->back( 152 LeftVector->push_back(RightVector->back());
253 RightVector->erase(RightVector->end()-1) 153 RightVector->erase(RightVector->end()-1);
254 } 154 }
255 delete RightVector; 155 delete RightVector;
256 156
257 CalculateHadronTimePosition(theString.Get4Mo 157 CalculateHadronTimePosition(theString.Get4Momentum().mag(), LeftVector);
258 158
259 G4LorentzRotation toObserverFrame(toCms.inve 159 G4LorentzRotation toObserverFrame(toCms.inverse());
260 160
261 for (size_t C1 = 0; C1 < LeftVector->size(); << 161 for(size_t C1 = 0; C1 < LeftVector->size(); C1++)
262 { 162 {
263 G4KineticTrack* Hadron = LeftVector->oper 163 G4KineticTrack* Hadron = LeftVector->operator[](C1);
264 G4LorentzVector Momentum = Hadron->Get4Mo 164 G4LorentzVector Momentum = Hadron->Get4Momentum();
265 Momentum = toObserverFrame*Momentum; 165 Momentum = toObserverFrame*Momentum;
266 Hadron->Set4Momentum(Momentum); 166 Hadron->Set4Momentum(Momentum);
267 G4LorentzVector Coordinate(Hadron->GetPos 167 G4LorentzVector Coordinate(Hadron->GetPosition(), Hadron->GetFormationTime());
268 Momentum = toObserverFrame*Coordinate; 168 Momentum = toObserverFrame*Coordinate;
269 Hadron->SetFormationTime(Momentum.e()); 169 Hadron->SetFormationTime(Momentum.e());
270 G4ThreeVector aPosition(Momentum.vect()); 170 G4ThreeVector aPosition(Momentum.vect());
271 Hadron->SetPosition(theString.GetPosition 171 Hadron->SetPosition(theString.GetPosition()+aPosition);
272 } 172 }
273 return LeftVector; 173 return LeftVector;
274 } << 174
275 175
276 //-------------------------------------------- <<
277 176
278 G4bool G4QGSMFragmentation::IsItFragmentable(c <<
279 { <<
280 return sqr( MinimalStringMass + MassCut ) < <<
281 } 177 }
282 178
283 //-------------------------------------------- 179 //----------------------------------------------------------------------------------------------------------
284 180
285 G4bool G4QGSMFragmentation::StopFragmenting(co << 181 G4double G4QGSMFragmentation::GetLightConeZ(G4double zmin, G4double zmax, G4int PartonEncoding, G4ParticleDefinition* , G4double , G4double )
286 { << 182 {
287 SetMinimalStringMass(string); << 183 G4double z;
288 if ( MinimalStringMass < 0.0 ) return << 184 G4double theA(0), d1, d2, yf;
289 << 185 G4int absCode = std::abs( PartonEncoding );
290 G4double smass = string->Mass(); << 186 if (absCode < 10)
291 G4double x = (string->IsAFourQuarkString()) << 187 {
292 : 0.66e-6*(smass - MinimalStringMass)*(sma << 188 if(absCode == 1 || absCode == 2) theA = arho;
293 << 189 else if(absCode == 3) theA = aphi;
294 G4bool res = true; << 190 else throw G4HadronicException(__FILE__, __LINE__, "Unknown PDGencoding in G4QGSMFragmentation::G4LightConeZ");
295 if(x > 0.0) { <<
296 res = (x < 200.) ? (G4UniformRand() <<
297 } <<
298 return res; <<
299 } <<
300 <<
301 //-------------------------------------------- <<
302 <<
303 G4KineticTrack * G4QGSMFragmentation::Splitup( <<
304 G4FragmentingStri <<
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 <<
365 delete newString; newString=0; <<
366 <<
367 G4KineticTrack * Hadron =0; <<
368 if ( HadronMomentum != 0 ) { <<
369 <<
370 #ifdef debug_QGSMfragmentation <<
371 G4cout<<"The attempt was successful <<
372 #endif <<
373 G4ThreeVector Pos; <<
374 Hadron = new G4KineticTrack(HadronDefinit <<
375 <<
376 newString=new G4FragmentingString(*st <<
377 <<
378 delete HadronMomentum; <<
379 } <<
380 else <<
381 { <<
382 #ifdef debug_QGSMfragmentation <<
383 G4cout<<"The attempt was not successf <<
384 #endif <<
385 } <<
386 <<
387 #ifdef debug_VStringDecay <<
388 G4cout<<"End SplitUP (G4VLongitudinalSt <<
389 #endif <<
390 <<
391 return Hadron; <<
392 } <<
393 <<
394 //-------------------------------------------- <<
395 <<
396 G4ParticleDefinition *G4QGSMFragmentation::DiQ <<
397 <<
398 { <<
399 //... can Diquark break or not? <<
400 if (G4UniformRand() < DiquarkBreakProb ) / <<
401 { <<
402 G4int stableQuarkEncoding = decay->GetPD <<
403 G4int decayQuarkEncoding = (decay->GetPD <<
404 <<
405 if (G4UniformRand() < 0.5) <<
406 { <<
407 G4int Swap = stableQuarkEncoding; <<
408 stableQuarkEncoding = decayQuarkEncod <<
409 decayQuarkEncoding = Swap; <<
410 } <<
411 <<
412 G4int IsParticle=(decayQuarkEncoding>0) <<
413 <<
414 G4double StrSup=GetStrangeSuppress(); <<
415 SetStrangenessSuppression((1.0 - 0.07)/2 <<
416 pDefPair QuarkPair = CreatePartonPair(Is <<
417 SetStrangenessSuppression(StrSup); <<
418 <<
419 //... Build new Diquark <<
420 G4int QuarkEncoding=QuarkPair.second->Ge <<
421 G4int i10 = std::max(std::abs(QuarkEnco <<
422 G4int i20 = std::min(std::abs(QuarkEnco <<
423 G4int spin = (i10 != i20 && G4UniformRan <<
424 G4int NewDecayEncoding = -1*IsParticle*( <<
425 <<
426 created = FindParticle(NewDecayEncoding) <<
427 G4ParticleDefinition * decayQuark=FindPa <<
428 G4ParticleDefinition * had=hadronizer->B <<
429 <<
430 DecayQuark = decay->GetPDGEncoding(); <<
431 NewQuark = NewDecayEncoding; <<
432 <<
433 return had; <<
434 <<
435 } else { //... Diquark does not break <<
436 <<
437 G4int IsParticle=(decay->GetPDGEncoding( <<
438 G4double StrSup=GetStrangeSuppress(); / <<
439 SetStrangenessSuppression((1.0 - 0.07)/2 <<
440 pDefPair QuarkPair = CreatePartonPair(Is <<
441 SetStrangenessSuppression(StrSup); <<
442 <<
443 created = QuarkPair.second; <<
444 191
445 DecayQuark = decay->GetPDGEncoding(); << 192 do
446 NewQuark = created->GetPDGEncoding(); << 193 {
>> 194 z = zmin + G4UniformRand() * (zmax - zmin);
>> 195 d1 = (1. - z);
>> 196 d2 = (alft - theA);
>> 197 yf = std::pow(d1, d2);
>> 198 }
>> 199 while (G4UniformRand() > yf);
>> 200 }
>> 201 else
>> 202 {
>> 203 if(absCode == 1103 || absCode == 2101 ||
>> 204 absCode == 2203 || absCode == 2103)
>> 205 {
>> 206 d2 = (alft - (2.*an - arho));
>> 207 }
>> 208 else if(absCode == 3101 || absCode == 3103 ||
>> 209 absCode == 3201 || absCode == 3203)
>> 210 {
>> 211 d2 = (alft - (2.*ala - arho));
>> 212 }
>> 213 else
>> 214 {
>> 215 d2 = (alft - (2.*aksi - arho));
>> 216 }
447 217
448 G4ParticleDefinition * had=hadronizer->B << 218 do
449 return had; << 219 {
450 } << 220 z = zmin + G4UniformRand() * (zmax - zmin);
>> 221 d1 = (1. - z);
>> 222 yf = std::pow(d1, d2);
>> 223 }
>> 224 while (G4UniformRand() > yf);
>> 225 }
>> 226 return z;
451 } 227 }
452 <<
453 //-------------------------------------------- 228 //-----------------------------------------------------------------------------------------
454 229
455 G4LorentzVector * G4QGSMFragmentation::SplitEa 230 G4LorentzVector * G4QGSMFragmentation::SplitEandP(G4ParticleDefinition * pHadron,
456 << 231 G4FragmentingString * string)
457 <<
458 { 232 {
459 G4double HadronMass = pHadron->GetPDGMa 233 G4double HadronMass = pHadron->GetPDGMass();
460 234
461 SetMinimalStringMass(NewString); <<
462 <<
463 if ( MinimalStringMass < 0.0 ) return n <<
464 <<
465 #ifdef debug_QGSMfragmentation <<
466 G4cout<<"G4QGSMFragmentation::SplitEand <<
467 G4cout<<"String 4 mom, String M "<<stri <<
468 G4cout<<"HadM MinimalStringMassLeft Str <<
469 <<string->Mass()<<" "<<HadronMass <<
470 #endif <<
471 <<
472 if (HadronMass + MinimalStringMass > st <<
473 { <<
474 #ifdef debug_QGSMfragmentation <<
475 G4cout<<"Mass of the string is not su <<
476 #endif <<
477 return 0; <<
478 } // have to start all over! <<
479 <<
480 // calculate and assign hadron transver 235 // calculate and assign hadron transverse momentum component HadronPx andHadronPy
481 G4double StringMT2 = string->MassT2(); << 236 G4ThreeVector thePt;
482 G4double StringMT = std::sqrt(StringMT << 237 thePt=SampleQuarkPt();
483 << 238 G4ThreeVector HadronPt = thePt +string->DecayPt();
484 G4LorentzVector String4Momentum = strin << 239 HadronPt.setZ(0);
485 String4Momentum.setPz(0.); <<
486 G4ThreeVector StringPt = String4Momentu <<
487 <<
488 G4ThreeVector HadronPt , RemSysPt; <<
489 G4double HadronMassT2, ResidualMas <<
490 <<
491 // Mt distribution is implemented <<
492 G4double HadronMt, Pt, Pt2, phi; <<
493 <<
494 //... sample Pt of the hadron <<
495 G4int attempt=0; <<
496 do <<
497 { <<
498 attempt++; if (attempt > StringLoopIn <<
499 <<
500 HadronMt = HadronMass - 200.0*G4Log(G <<
501 Pt2 = sqr(HadronMt)-sqr(HadronMass); <<
502 phi = 2.*pi*G4UniformRand(); <<
503 G4ThreeVector SampleQuarkPtw= G4Three <<
504 HadronPt =SampleQuarkPtw + string->D <<
505 HadronPt.setZ(0); <<
506 RemSysPt = StringPt - HadronPt; <<
507 <<
508 HadronMassT2 = sqr(HadronMass) + Hadr <<
509 ResidualMassT2=sqr(MinimalStringMass) <<
510 <<
511 } while (std::sqrt(HadronMassT2) + std: <<
512 <<
513 //... sample z to define hadron longit 240 //... sample z to define hadron longitudinal momentum and energy
514 //... but first check the available pha 241 //... but first check the available phase space
>> 242 G4double DecayQuarkMass2 = sqr(string->GetDecayParton()->GetPDGMass());
>> 243 G4double HadronMass2T = sqr(HadronMass) + HadronPt.mag2();
>> 244 if (DecayQuarkMass2 + HadronMass2T >= SmoothParam*(string->Mass2()) )
>> 245 return 0; // have to start all over!
515 246
516 G4double Pz2 = (sqr(StringMT2 - HadronM << 247 //... then compute allowed z region z_min <= z <= z_max
517 4*HadronMassT2 * ResidualMassT2)/4./ << 248
518 << 249 G4double zMin = HadronMass2T/(string->Mass2());
519 if ( Pz2 < 0 ) {return 0;} // << 250 G4double zMax = 1. - DecayQuarkMass2/(string->Mass2());
520 <<
521 //... then compute allowed z region z_ <<
522 <<
523 G4double Pz = std::sqrt(Pz2); <<
524 G4double zMin = (std::sqrt(HadronMassT2 <<
525 G4double zMax = (std::sqrt(HadronMassT2 <<
526 <<
527 if (zMin >= zMax) return 0; // have 251 if (zMin >= zMax) return 0; // have to start all over!
528 252
529 G4double z = GetLightConeZ(zMin, zMax, 253 G4double z = GetLightConeZ(zMin, zMax,
530 string->GetDecayParton() << 254 string->GetDecayParton()->GetPDGEncoding(), pHadron,
531 HadronPt.x(), HadronPt.y << 255 HadronPt.x(), HadronPt.y());
532 << 256
533 //... now compute hadron longitudinal m 257 //... now compute hadron longitudinal momentum and energy
534 // longitudinal hadron momentum compone 258 // longitudinal hadron momentum component HadronPz
535 259
536 HadronPt.setZ( 0.5* string->GetDecayDir << 260 HadronPt.setZ(0.5* string->GetDecayDirection() *
537 (z * string->LightConeDecay() - << 261 (z * string->LightConeDecay() -
538 HadronMassT2/(z * string->LightCone << 262 HadronMass2T/(z * string->LightConeDecay())));
539 G4double HadronE = 0.5* (z * string->L << 263 G4double HadronE = 0.5* (z * string->LightConeDecay() +
540 HadronMassT2/(z * string->LightConeDe << 264 HadronMass2T/(z * string->LightConeDecay()));
541 265
542 G4LorentzVector * a4Momentum= new G4Lor 266 G4LorentzVector * a4Momentum= new G4LorentzVector(HadronPt,HadronE);
543 267
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 <<
551 return a4Momentum; 268 return a4Momentum;
552 } 269 }
553 270
554 //-------------------------------------------- <<
555 <<
556 G4double G4QGSMFragmentation::GetLightConeZ(G4 <<
557 G4 <<
558 { <<
559 G4double lambda = 2.0*(sqr(ptx)+sqr(pty))/sq <<
560 <<
561 #ifdef debug_QGSMfragmentation <<
562 G4cout<<"GetLightConeZ zmin zmax Parton pHad <<
563 <<" "/*<< pHadron->GetParticleName() * <<
564 #endif <<
565 <<
566 G4double z(0.); <<
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 } <<
622 271
623 //-------------------------------------------- 272 //-----------------------------------------------------------------------------------------
624 273
625 G4bool G4QGSMFragmentation::SplitLast(G4Fragme 274 G4bool G4QGSMFragmentation::SplitLast(G4FragmentingString * string,
626 G4KineticTrackVector * Lef << 275 G4KineticTrackVector * LeftVector,
627 G4KineticTrackVector * Right << 276 G4KineticTrackVector * RightVector)
628 { 277 {
629 //... perform last cluster decay 278 //... perform last cluster decay
630 <<
631 G4ThreeVector ClusterVel =string->Get4Mome 279 G4ThreeVector ClusterVel =string->Get4Momentum().boostVector();
632 G4double ResidualMass =string->Mass(); << 280 G4double ResidualMass =string->Mass();
633 << 281 G4double ClusterMassCut = ClusterMass;
634 #ifdef debug_QGSMfragmentation <<
635 G4cout<<"Split last----------------------- <<
636 G4cout<<"StrMass "<<ResidualMass<<" q's " <<
637 <<string->GetLeftParton()->GetPartic <<
638 <<string->GetRightParton()->GetParti <<
639 #endif <<
640 <<
641 G4int cClusterInterrupt = 0; 282 G4int cClusterInterrupt = 0;
642 G4ParticleDefinition *LeftHadron = nullptr << 283 G4ParticleDefinition * LeftHadron, * RightHadron;
643 G4ParticleDefinition *RightHadron = nullpt <<
644 const G4int maxNumberOfLoops = 1000; <<
645 G4int loopCounter = 0; <<
646 <<
647 G4double LeftHadronMass(0.); G4double Righ <<
648 do 284 do
649 { 285 {
650 if (cClusterInterrupt++ >= ClusterLoop << 286 if (cClusterInterrupt++ >= ClusterLoopInterrupt)
651 LeftHadronMass = -MaxMass; RightHadron << 287 {
652 << 288 return false;
653 G4ParticleDefinition * quark = nullptr; << 289 }
>> 290 G4ParticleDefinition * quark = NULL;
654 string->SetLeftPartonStable(); // to query q 291 string->SetLeftPartonStable(); // to query quark contents..
655 <<
656 if (string->DecayIsQuark() && string->Stable 292 if (string->DecayIsQuark() && string->StableIsQuark() )
657 { 293 {
658 //... there are quarks on cluster ends << 294 //... there are quarks on cluster ends
659 << 295 LeftHadron= QuarkSplitup(string->GetLeftParton(), quark);
660 G4int IsParticle=(string->GetLeftParton()- << 296 } else {
661 // if we have a quark, we need << 297 //... there is a Diquark on cluster ends
662 << 298 G4int IsParticle;
663 pDefPair QuarkPair = CreatePartonPair(IsPa << 299 if ( string->StableIsQuark() ) {
664 quark = QuarkPair.second; << 300 IsParticle=(string->GetLeftParton()->GetPDGEncoding()>0) ? -1 : +1;
665 << 301 } else {
666 LeftHadron= hadronizer->Build(QuarkPair.fi << 302 IsParticle=(string->GetLeftParton()->GetPDGEncoding()>0) ? +1 : -1;
667 if ( LeftHadron == NULL ) continue; << 303 }
668 RightHadron = hadronizer->Build(stri << 304 pDefPair QuarkPair = CreatePartonPair(IsParticle,false); // no diquarks wanted
669 if ( RightHadron == NULL ) continue; << 305 quark = QuarkPair.second;
670 } else if( (!string->DecayIsQuark() && stri << 306 LeftHadron=hadronizer->Build(QuarkPair.first, string->GetLeftParton());
671 ( string->DecayIsQuark() && !stri << 307 }
672 //... there is a Diquark on one of cluster << 308 RightHadron = hadronizer->Build(string->GetRightParton(), quark);
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 309
712 if ( loopCounter >= maxNumberOfLoops ) { << 310 //... repeat procedure, if mass of cluster is too low to produce hadrons
713 return false; << 311 //... ClusterMassCut = 0.15*GeV model parameter
714 } << 312 if ( quark->GetParticleSubType()== "quark" ) {ClusterMassCut = 0.;}
>> 313 else {ClusterMassCut = ClusterMass;}
>> 314 }
>> 315 while (ResidualMass <= LeftHadron->GetPDGMass() + RightHadron->GetPDGMass() + ClusterMassCut);
715 316
716 //... compute hadron momenta and energies 317 //... compute hadron momenta and energies
717 G4LorentzVector LeftMom, RightMom; 318 G4LorentzVector LeftMom, RightMom;
718 G4ThreeVector Pos; 319 G4ThreeVector Pos;
719 Sample4Momentum(&LeftMom , LeftHadron->Get << 320 Sample4Momentum(&LeftMom, LeftHadron->GetPDGMass(), &RightMom, RightHadron->GetPDGMass(), ResidualMass);
720 &RightMom, RightHadron->Ge <<
721 LeftMom.boost(ClusterVel); 321 LeftMom.boost(ClusterVel);
722 RightMom.boost(ClusterVel); 322 RightMom.boost(ClusterVel);
723 <<
724 #ifdef debug_QGSMfragmentation <<
725 G4cout<<LeftHadron->GetParticleName()<<" " <<
726 G4cout<<"Left Hadrom P M "<<LeftMom<<" "< <<
727 G4cout<<"Right Hadrom P M "<<RightMom<<" " <<
728 #endif <<
729 <<
730 LeftVector->push_back(new G4KineticTrack(L 323 LeftVector->push_back(new G4KineticTrack(LeftHadron, 0, Pos, LeftMom));
731 RightVector->push_back(new G4KineticTrack( 324 RightVector->push_back(new G4KineticTrack(RightHadron, 0, Pos, RightMom));
732 325
733 return true; 326 return true;
>> 327
734 } 328 }
735 329
736 //-------------------------------------------- 330 //----------------------------------------------------------------------------------------------------------
737 331
738 void G4QGSMFragmentation::Sample4Momentum(G4Lo << 332 G4bool G4QGSMFragmentation::IsFragmentable(const G4FragmentingString * const string)
739 G4Lo <<
740 { 333 {
741 #ifdef debug_QGSMfragmentation << 334 return sqr(FragmentationMass(string)+MassCut) <
742 G4cout<<"Sample4Momentum Last------------- << 335 string->Mass2();
743 G4cout<<" StrMass "<<InitialMass<<" Mass1 << 336 }
744 G4cout<<" SumMass "<<Mass+AntiMass<<G4end << 337
745 #endif << 338 //----------------------------------------------------------------------------------------------------------
>> 339
>> 340 G4bool G4QGSMFragmentation::StopFragmenting(const G4FragmentingString * const string)
>> 341 {
>> 342 return
>> 343 sqr(FragmentationMass(string,&G4HadronBuilder::BuildHighSpin)+MassCut) >
>> 344 string->Get4Momentum().mag2();
>> 345 }
746 346
>> 347 //----------------------------------------------------------------------------------------------------------
>> 348
>> 349 //----------------------------------------------------------------------------------------------------------
>> 350
>> 351 void G4QGSMFragmentation::Sample4Momentum(G4LorentzVector* Mom, G4double Mass, G4LorentzVector* AntiMom, G4double AntiMass, G4double InitialMass)
>> 352 {
747 G4double r_val = sqr(InitialMass*InitialMa 353 G4double r_val = sqr(InitialMass*InitialMass - Mass*Mass - AntiMass*AntiMass) - sqr(2.*Mass*AntiMass);
748 G4double Pabs = (r_val > 0.)? std::sqrt(r_ 354 G4double Pabs = (r_val > 0.)? std::sqrt(r_val)/(2.*InitialMass) : 0;
749 355
750 //... sample unit vector 356 //... sample unit vector
751 G4double pz = 1. - 2.*G4UniformRand(); 357 G4double pz = 1. - 2.*G4UniformRand();
752 G4double st = std::sqrt(1. - pz * pz)* 358 G4double st = std::sqrt(1. - pz * pz)*Pabs;
753 G4double phi = 2.*pi*G4UniformRand(); 359 G4double phi = 2.*pi*G4UniformRand();
754 G4double px = st*std::cos(phi); 360 G4double px = st*std::cos(phi);
755 G4double py = st*std::sin(phi); 361 G4double py = st*std::sin(phi);
756 pz *= Pabs; 362 pz *= Pabs;
757 363
758 Mom->setPx(px); Mom->setPy(py); Mom->setPz 364 Mom->setPx(px); Mom->setPy(py); Mom->setPz(pz);
759 Mom->setE(std::sqrt(Pabs*Pabs + Mass*Mass) 365 Mom->setE(std::sqrt(Pabs*Pabs + Mass*Mass));
760 366
761 AntiMom->setPx(-px); AntiMom->setPy(-py); 367 AntiMom->setPx(-px); AntiMom->setPy(-py); AntiMom->setPz(-pz);
762 AntiMom->setE (std::sqrt(Pabs*Pabs + AntiM 368 AntiMom->setE (std::sqrt(Pabs*Pabs + AntiMass*AntiMass));
763 } << 369 }
764 <<
765 //-------------------------------------------- <<
766 <<
767 void G4QGSMFragmentation::SetFFq2q() // q-> q <<
768 { <<
769 for (G4int i=0; i < 5; i++) { <<
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 } <<
777 <<
778 //-------------------------------------------- <<
779 <<
780 void G4QGSMFragmentation::SetFFq2qq() // q-> <<
781 { <<
782 for (G4int i=0; i < 5; i++) { <<
783 FFq2qq[i][ 0][0] = 0.0 ; FFq2qq[i][ 0][1] <<
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 } <<
800 <<
801 //-------------------------------------------- <<
802 <<
803 void G4QGSMFragmentation::SetFFqq2q() // q1q2 <<
804 { <<
805 for (G4int i=0; i < 15; i++) { <<
806 FFqq2q[i][0][0] = 2.0*(arho - an); FFqq2q[ <<
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 <<
814 //-------------------------------------------- <<
815 <<
816 void G4QGSMFragmentation::SetFFqq2qq() // q1( <<
817 { <<
818 for (G4int i=0; i < 15; i++) { <<
819 FFqq2qq[i][0][0] = 0. ; FFqq2qq[i][0][1] <<
820 FFqq2qq[i][1][0] = 0. ; FFqq2qq[i][1][1] <<
821 FFqq2qq[i][2][0] = 0. ; FFqq2qq[i][2][1] <<
822 FFqq2qq[i][3][0] = 0. ; FFqq2qq[i][3][1] <<
823 FFqq2qq[i][4][0] = 0. ; FFqq2qq[i][4][1] <<
824 } <<
825 } <<
826 370
>> 371
>> 372 //*********************************************************************************************
827 373