| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 //
27 //
28 // -------------------------------------------
29 // GEANT 4 class implementation file
30 //
31 // History: first implementation, Maxim K
32 // -------------------------------------------
33 #include "G4QGSMFragmentation.hh"
34 #include "G4PhysicalConstants.hh"
35 #include "G4SystemOfUnits.hh"
36 #include "Randomize.hh"
37 #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
46 // Class G4QGSMFragmentation
47 //********************************************
48
49 G4QGSMFragmentation::G4QGSMFragmentation()
50 {
51 SigmaQT = 0.45 * GeV;
52
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
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
111 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 {
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 }
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
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
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
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;
552 }
553
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
623 //--------------------------------------------
624
625 G4bool G4QGSMFragmentation::SplitLast(G4Fragme
626 G4KineticTrackVector * Lef
627 G4KineticTrackVector * Right
628 {
629 //... perform last cluster decay
630
631 G4ThreeVector ClusterVel =string->Get4Mome
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 {
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;
714 }
715
716 //... compute hadron momenta and energies
717 G4LorentzVector LeftMom, RightMom;
718 G4ThreeVector Pos;
719 Sample4Momentum(&LeftMom , LeftHadron->Get
720 &RightMom, RightHadron->Ge
721 LeftMom.boost(ClusterVel);
722 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
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 }
825 }
826
827