| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer << 3 // * DISCLAIMER *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th << 5 // * The following disclaimer summarizes all the specific disclaimers *
6 // * the Geant4 Collaboration. It is provided << 6 // * of contributors to this software. The specific disclaimers,which *
7 // * conditions of the Geant4 Software License << 7 // * govern, are listed with their locations in: *
8 // * LICENSE and available at http://cern.ch/ << 8 // * http://cern.ch/geant4/license *
9 // * include a list of copyright holders. <<
10 // * 9 // * *
11 // * Neither the authors of this software syst 10 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 11 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 12 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 13 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file << 14 // * use. *
16 // * for the full disclaimer and the limitatio <<
17 // * 15 // * *
18 // * This code implementation is the result << 16 // * This code implementation is the intellectual property of the *
19 // * technical work of the GEANT4 collaboratio << 17 // * GEANT4 collaboration. *
20 // * By using, copying, modifying or distri << 18 // * By copying, distributing or modifying the Program (or any work *
21 // * any work based on the software) you ag << 19 // * based on the Program) you indicate your acceptance of this *
22 // * use in resulting scientific publicati << 20 // * statement, and all its terms. *
23 // * acceptance of all terms of the Geant4 Sof <<
24 // ******************************************* 21 // ********************************************************************
25 // 22 //
26 // 23 //
>> 24 // $Id: G4QGSMFragmentation.cc,v 1.4 2005/06/04 13:47:01 jwellisc Exp $
>> 25 // GEANT4 tag $Name: geant4-08-00 $
27 // 26 //
28 // ------------------------------------------- 27 // -----------------------------------------------------------------------------
29 // GEANT 4 class implementation file 28 // GEANT 4 class implementation file
30 // 29 //
31 // History: first implementation, Maxim K 30 // History: first implementation, Maxim Komogorov, 10-Jul-1998
32 // ------------------------------------------- 31 // -----------------------------------------------------------------------------
33 #include "G4QGSMFragmentation.hh" 32 #include "G4QGSMFragmentation.hh"
34 #include "G4PhysicalConstants.hh" <<
35 #include "G4SystemOfUnits.hh" <<
36 #include "Randomize.hh" 33 #include "Randomize.hh"
37 #include "G4ios.hh" 34 #include "G4ios.hh"
38 #include "G4FragmentingString.hh" <<
39 #include "G4DiQuarks.hh" <<
40 #include "G4Quarks.hh" <<
41 #include "G4HadronicParameters.hh" <<
42 #include "G4Pow.hh" <<
43 <<
44 //#define debug_QGSMfragmentation <<
45 35
46 // Class G4QGSMFragmentation 36 // Class G4QGSMFragmentation
47 //******************************************** 37 //****************************************************************************************
48 38
49 G4QGSMFragmentation::G4QGSMFragmentation() << 39 G4QGSMFragmentation::G4QGSMFragmentation() :
50 { << 40 arho(0.5), aphi(0.), an(-0.5), ala(-0.75), aksi(-1.), alft(0.5)
51 SigmaQT = 0.45 * GeV; << 41 {
52 << 42 }
53 MassCut = 0.35*GeV; <<
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 43
94 SetFFq2q(); << 44 G4QGSMFragmentation::G4QGSMFragmentation(const G4QGSMFragmentation &) : G4VLongitudinalStringDecay(),
95 SetFFq2qq(); << 45 arho(0.5), aphi(0.), an(-0.5), ala(-0.75), aksi(-1.), alft(0.5)
96 SetFFqq2q(); << 46 {
97 SetFFqq2qq(); << 47 }
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 48
111 G4QGSMFragmentation::~G4QGSMFragmentation() 49 G4QGSMFragmentation::~G4QGSMFragmentation()
112 {} <<
113 <<
114 //-------------------------------------------- <<
115 <<
116 G4KineticTrackVector* G4QGSMFragmentation::Fra <<
117 { <<
118 <<
119 G4FragmentingString aString(theString); <<
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 <<
138 // Check if string has enough mass to fragme <<
139 G4KineticTrackVector * LeftVector=NULL; <<
140 <<
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 <<
156 LeftVector = new G4KineticTrackVector; <<
157 G4KineticTrackVector * RightVector=new G4Kin <<
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 <<
172 G4FragmentingString *currentString=new G4F <<
173 <<
174 std::for_each(LeftVector->begin(), LeftVec <<
175 LeftVector->clear(); <<
176 std::for_each(RightVector->begin(), RightV <<
177 RightVector->clear(); <<
178 <<
179 inner_sucess=true; // set false on failur <<
180 const G4int maxNumberOfLoops = <<
181 G4int loopCounter = -1; <<
182 while (! StopFragmenting(currentString) && <<
183 { // Split current string into hadron + n <<
184 <<
185 #ifdef debug_QGSMfragm <<
186 G4cout<<"The string ca <<
187 #endif <<
188 G4FragmentingString *newString=0; // us <<
189 G4KineticTrack * Hadron=Splitup(currentS <<
190 <<
191 if ( Hadron != 0 ) <<
192 { <<
193 #ifdef debug_QGSMfr <<
194 G4cout<<"Hadron pro <<
195 #endif <<
196 // To close the pro <<
197 if ( currentString->GetDecayDirection <<
198 LeftVector->push_back(Hadron); <<
199 else <<
200 RightVector->push_back(Hadron); <<
201 <<
202 delete currentString; <<
203 currentString=newString; <<
204 <<
205 } else { <<
206 <<
207 #ifdef debug_QGSMfr <<
208 G4cout<<"abandon .. <<
209 #endif <<
210 <<
211 // Abandon ... start from the beginni <<
212 if (newString) delete newString; <<
213 inner_sucess=false; <<
214 break; <<
215 } <<
216 } <<
217 if ( loopCounter >= maxNumberO <<
218 inner_sucess=false; <<
219 } <<
220 <<
221 // Split current string into 2 final Hadro <<
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 && <<
231 SplitLast(currentString,LeftVector, R <<
232 { <<
233 success=true; <<
234 } <<
235 delete currentString; <<
236 } <<
237 <<
238 delete theStringInCMS; <<
239 <<
240 if ( ! success ) <<
241 { <<
242 std::for_each(LeftVector->begin(), LeftVec <<
243 LeftVector->clear(); <<
244 std::for_each(RightVector->begin(), RightV <<
245 delete RightVector; <<
246 return LeftVector; <<
247 } <<
248 <<
249 // Join Left- and RightVector into LeftVecto <<
250 while(!RightVector->empty()) /* Loop checki <<
251 { <<
252 LeftVector->push_back(RightVector->back( <<
253 RightVector->erase(RightVector->end()-1) <<
254 } <<
255 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 <<
276 //-------------------------------------------- <<
277 <<
278 G4bool G4QGSMFragmentation::IsItFragmentable(c <<
279 { <<
280 return sqr( MinimalStringMass + MassCut ) < <<
281 } <<
282 <<
283 //-------------------------------------------- <<
284 <<
285 G4bool G4QGSMFragmentation::StopFragmenting(co <<
286 { <<
287 SetMinimalStringMass(string); <<
288 if ( MinimalStringMass < 0.0 ) return <<
289 <<
290 G4double smass = string->Mass(); <<
291 G4double x = (string->IsAFourQuarkString()) <<
292 : 0.66e-6*(smass - MinimalStringMass)*(sma <<
293 <<
294 G4bool res = true; <<
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 { 50 {
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 <<
445 DecayQuark = decay->GetPDGEncoding(); <<
446 NewQuark = created->GetPDGEncoding(); <<
447 <<
448 G4ParticleDefinition * had=hadronizer->B <<
449 return had; <<
450 } 51 }
451 } <<
452 <<
453 //-------------------------------------------- <<
454 <<
455 G4LorentzVector * G4QGSMFragmentation::SplitEa <<
456 <<
457 <<
458 { <<
459 G4double HadronMass = pHadron->GetPDGMa <<
460 <<
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 <<
481 G4double StringMT2 = string->MassT2(); <<
482 G4double StringMT = std::sqrt(StringMT <<
483 <<
484 G4LorentzVector String4Momentum = strin <<
485 String4Momentum.setPz(0.); <<
486 G4ThreeVector StringPt = String4Momentu <<
487 <<
488 G4ThreeVector HadronPt , RemSysPt; <<
489 G4double HadronMassT2, ResidualMas <<
490 52
491 // Mt distribution is implemented << 53 //****************************************************************************************
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 <<
514 //... but first check the available pha <<
515 <<
516 G4double Pz2 = (sqr(StringMT2 - HadronM <<
517 4*HadronMassT2 * ResidualMassT2)/4./ <<
518 <<
519 if ( Pz2 < 0 ) {return 0;} // <<
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 <<
528 <<
529 G4double z = GetLightConeZ(zMin, zMax, <<
530 string->GetDecayParton() <<
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 54
544 #ifdef debug_QGSMfragmentation << 55 const G4QGSMFragmentation & G4QGSMFragmentation::operator=(const G4QGSMFragmentation &)
545 G4cout<<"string->GetDecayDirection() st << 56 {
546 <<string->GetDecayDirection()<<" << 57 throw G4HadronicException(__FILE__, __LINE__, "G4QGSMFragmentation::operator= meant to not be accessable");
547 G4cout<<"HadronPt,HadronE "<<HadronPt<< << 58 return *this;
548 G4cout<<"Out of QGSM SplitEandP "<<G4en << 59 }
549 #endif <<
550 60
551 return a4Momentum; << 61 int G4QGSMFragmentation::operator==(const G4QGSMFragmentation &right) const
552 } << 62 {
>> 63 return !memcmp(this, &right, sizeof(G4QGSMFragmentation));
>> 64 }
553 65
554 //-------------------------------------------- << 66 int G4QGSMFragmentation::operator!=(const G4QGSMFragmentation &right) const
>> 67 {
>> 68 return memcmp(this, &right, sizeof(G4QGSMFragmentation));
>> 69 }
>> 70
>> 71 //****************************************************************************************
555 72
556 G4double G4QGSMFragmentation::GetLightConeZ(G4 << 73 G4double G4QGSMFragmentation::GetLightConeZ(G4double zmin, G4double zmax, G4int PartonEncoding, G4ParticleDefinition* , G4double , G4double )
557 G4 <<
558 { 74 {
559 G4double lambda = 2.0*(sqr(ptx)+sqr(pty))/sq << 75 G4double z;
560 << 76 G4double theA(0), d1, d2, yf;
561 #ifdef debug_QGSMfragmentation << 77 G4int absCode = std::abs( PartonEncoding );
562 G4cout<<"GetLightConeZ zmin zmax Parton pHad << 78 if (absCode < 10)
563 <<" "/*<< pHadron->GetParticleName() * << 79 {
564 #endif << 80 if(absCode == 1 || absCode == 2) theA = arho;
565 << 81 else if(absCode == 3) theA = aphi;
566 G4double z(0.); << 82 else throw G4HadronicException(__FILE__, __LINE__, "Unknown PDGencoding in G4QGSMFragmentation::G4LightConeZ");
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 <<
623 //-------------------------------------------- <<
624 <<
625 G4bool G4QGSMFragmentation::SplitLast(G4Fragme <<
626 G4KineticTrackVector * Lef <<
627 G4KineticTrackVector * Right <<
628 { <<
629 //... perform last cluster decay <<
630 83
631 G4ThreeVector ClusterVel =string->Get4Mome << 84 do
632 G4double ResidualMass =string->Mass(); <<
633 <<
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; <<
642 G4ParticleDefinition *LeftHadron = nullptr <<
643 G4ParticleDefinition *RightHadron = nullpt <<
644 const G4int maxNumberOfLoops = 1000; <<
645 G4int loopCounter = 0; <<
646 <<
647 G4double LeftHadronMass(0.); G4double Righ <<
648 do <<
649 { 85 {
650 if (cClusterInterrupt++ >= ClusterLoop << 86 z = zmin + G4UniformRand() * (zmax - zmin);
651 LeftHadronMass = -MaxMass; RightHadron << 87 d1 = (1. - z);
652 << 88 d2 = (alft - theA);
653 G4ParticleDefinition * quark = nullptr; << 89 yf = std::pow(d1, d2);
654 string->SetLeftPartonStable(); // to query q << 90 }
655 << 91 while (G4UniformRand() > yf);
656 if (string->DecayIsQuark() && string->Stable << 92 }
657 { << 93 else
658 //... there are quarks on cluster ends << 94 {
659 << 95 if(absCode == 1103 || absCode == 2101 ||
660 G4int IsParticle=(string->GetLeftParton()- << 96 absCode == 2203 || absCode == 2103)
661 // if we have a quark, we need << 97 {
662 << 98 d2 = (alft - (2.*an - arho));
663 pDefPair QuarkPair = CreatePartonPair(IsPa << 99 }
664 quark = QuarkPair.second; << 100 else if(absCode == 3101 || absCode == 3103 ||
665 << 101 absCode == 3201 || absCode == 3203)
666 LeftHadron= hadronizer->Build(QuarkPair.fi << 102 {
667 if ( LeftHadron == NULL ) continue; << 103 d2 = (alft - (2.*ala - arho));
668 RightHadron = hadronizer->Build(stri << 104 }
669 if ( RightHadron == NULL ) continue; << 105 else
670 } else if( (!string->DecayIsQuark() && stri << 106 {
671 ( string->DecayIsQuark() && !stri << 107 d2 = (alft - (2.*aksi - arho));
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; <<
714 } 108 }
715 109
716 //... compute hadron momenta and energies << 110 do
717 G4LorentzVector LeftMom, RightMom; << 111 {
718 G4ThreeVector Pos; << 112 z = zmin + G4UniformRand() * (zmax - zmin);
719 Sample4Momentum(&LeftMom , LeftHadron->Get << 113 d1 = (1. - z);
720 &RightMom, RightHadron->Ge << 114 yf = std::pow(d1, d2);
721 LeftMom.boost(ClusterVel); << 115 }
722 RightMom.boost(ClusterVel); << 116 while (G4UniformRand() > yf);
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 <<
731 RightVector->push_back(new G4KineticTrack( <<
732 <<
733 return true; <<
734 } <<
735 <<
736 //-------------------------------------------- <<
737 <<
738 void G4QGSMFragmentation::Sample4Momentum(G4Lo <<
739 G4Lo <<
740 { <<
741 #ifdef debug_QGSMfragmentation <<
742 G4cout<<"Sample4Momentum Last------------- <<
743 G4cout<<" StrMass "<<InitialMass<<" Mass1 <<
744 G4cout<<" SumMass "<<Mass+AntiMass<<G4end <<
745 #endif <<
746 <<
747 G4double r_val = sqr(InitialMass*InitialMa <<
748 G4double Pabs = (r_val > 0.)? std::sqrt(r_ <<
749 <<
750 //... sample unit vector <<
751 G4double pz = 1. - 2.*G4UniformRand(); <<
752 G4double st = std::sqrt(1. - pz * pz)* <<
753 G4double phi = 2.*pi*G4UniformRand(); <<
754 G4double px = st*std::cos(phi); <<
755 G4double py = st*std::sin(phi); <<
756 pz *= Pabs; <<
757 <<
758 Mom->setPx(px); Mom->setPy(py); Mom->setPz <<
759 Mom->setE(std::sqrt(Pabs*Pabs + Mass*Mass) <<
760 <<
761 AntiMom->setPx(-px); AntiMom->setPy(-py); <<
762 AntiMom->setE (std::sqrt(Pabs*Pabs + AntiM <<
763 } <<
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 } 117 }
>> 118 return z;
825 } 119 }
826 << 120
>> 121 //*********************************************************************************************
827 122