source: trunk/source/processes/hadronic/models/parton_string/hadronization/src/G4LundStringFragmentation.cc@ 1338

Last change on this file since 1338 was 1337, checked in by garnier, 15 years ago

tag geant4.9.4 beta 1 + modifs locales
File size: 24.5 KB
Line 
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26//
27// $Id: G4LundStringFragmentation.cc,v 1.20 2010/06/21 17:50:48 vuzhinsk Exp $
28// GEANT4 tag $Name: geant4-09-04-beta-01 $ 1.8
29//
30// -----------------------------------------------------------------------------
31// GEANT 4 class implementation file
32//
33// History: first implementation, Maxim Komogorov, 10-Jul-1998
34// -----------------------------------------------------------------------------
35#include "G4LundStringFragmentation.hh"
36#include "G4FragmentingString.hh"
37#include "G4DiQuarks.hh"
38#include "G4Quarks.hh"
39
40#include "Randomize.hh"
41
42// Class G4LundStringFragmentation
43//*************************************************************************************
44
45G4LundStringFragmentation::G4LundStringFragmentation()
46 {
47// ------ For estimation of a minimal string mass ---------------
48 Mass_of_light_quark =140.*MeV;
49 Mass_of_heavy_quark =500.*MeV;
50 Mass_of_string_junction=720.*MeV;
51// ------ An estimated minimal string mass ----------------------
52 MinimalStringMass = 0.;
53 MinimalStringMass2 = 0.;
54// ------ Minimal invariant mass used at a string fragmentation -
55 WminLUND = 0.7*GeV; // Uzhi 0.8 1.5
56// ------ Smooth parameter used at a string fragmentation for ---
57// ------ smearinr sharp mass cut-off ---------------------------
58 SmoothParam = 0.2;
59
60// SetStringTensionParameter(0.25);
61 SetStringTensionParameter(1.);
62 }
63
64// --------------------------------------------------------------
65G4LundStringFragmentation::G4LundStringFragmentation(const G4LundStringFragmentation &) : G4VLongitudinalStringDecay()
66 {
67 }
68
69G4LundStringFragmentation::~G4LundStringFragmentation()
70 {
71 }
72
73//*************************************************************************************
74
75const G4LundStringFragmentation & G4LundStringFragmentation::operator=(const G4LundStringFragmentation &)
76 {
77 throw G4HadronicException(__FILE__, __LINE__, "G4LundStringFragmentation::operator= meant to not be accessable");
78 return *this;
79 }
80
81int G4LundStringFragmentation::operator==(const G4LundStringFragmentation &right) const
82 {
83 return !memcmp(this, &right, sizeof(G4LundStringFragmentation));
84 }
85
86int G4LundStringFragmentation::operator!=(const G4LundStringFragmentation &right) const
87 {
88 return memcmp(this, &right, sizeof(G4LundStringFragmentation));
89 }
90
91//--------------------------------------------------------------------------------------
92void G4LundStringFragmentation::SetMinimalStringMass(const G4FragmentingString * const string)
93{
94 G4double EstimatedMass=0.;
95 G4int Number_of_quarks=0;
96
97 G4int Qleft =std::abs(string->GetLeftParton()->GetPDGEncoding());
98
99 if( Qleft > 1000)
100 {
101 Number_of_quarks+=2;
102 G4int q1=Qleft/1000;
103 if( q1 < 3) {EstimatedMass +=Mass_of_light_quark;}
104 if( q1 > 2) {EstimatedMass +=Mass_of_heavy_quark;}
105
106 G4int q2=(Qleft/100)%10;
107 if( q2 < 3) {EstimatedMass +=Mass_of_light_quark;}
108 if( q2 > 2) {EstimatedMass +=Mass_of_heavy_quark;}
109 EstimatedMass +=Mass_of_string_junction;
110 }
111 else
112 {
113 Number_of_quarks++;
114 if( Qleft < 3) {EstimatedMass +=Mass_of_light_quark;}
115 if( Qleft > 2) {EstimatedMass +=Mass_of_heavy_quark;}
116 }
117
118 G4int Qright=std::abs(string->GetRightParton()->GetPDGEncoding());
119
120 if( Qright > 1000)
121 {
122 Number_of_quarks+=2;
123 G4int q1=Qright/1000;
124 if( q1 < 3) {EstimatedMass +=Mass_of_light_quark;}
125 if( q1 > 2) {EstimatedMass +=Mass_of_heavy_quark;}
126
127 G4int q2=(Qright/100)%10;
128 if( q2 < 3) {EstimatedMass +=Mass_of_light_quark;}
129 if( q2 > 2) {EstimatedMass +=Mass_of_heavy_quark;}
130 EstimatedMass +=Mass_of_string_junction;
131 }
132 else
133 {
134 Number_of_quarks++;
135 if( Qright < 3) {EstimatedMass +=Mass_of_light_quark;}
136 if( Qright > 2) {EstimatedMass +=Mass_of_heavy_quark;}
137 }
138
139 if(Number_of_quarks==2){EstimatedMass +=100.*MeV;}
140 if(Number_of_quarks==3){EstimatedMass += 20.*MeV;}
141 if(Number_of_quarks==4){EstimatedMass -=2.*Mass_of_string_junction;
142 if(EstimatedMass <= 1600.*MeV){EstimatedMass-=200.*MeV;}
143 else {EstimatedMass+=100.*MeV;}
144 }
145 MinimalStringMass=EstimatedMass;
146 SetMinimalStringMass2(EstimatedMass);
147}
148
149//--------------------------------------------------------------------------------------
150void G4LundStringFragmentation::SetMinimalStringMass2(
151 const G4double aValue)
152{
153 MinimalStringMass2=aValue * aValue;
154}
155
156//--------------------------------------------------------------------------------------
157G4KineticTrackVector* G4LundStringFragmentation::FragmentString(
158 const G4ExcitedString& theString)
159{
160// Can no longer modify Parameters for Fragmentation.
161 PastInitPhase=true;
162
163// check if string has enough mass to fragment...
164
165 SetMassCut(160.*MeV); // For LightFragmentationTest it is required
166 // that no one pi-meson can be produced
167/*
168G4cout<<G4endl<<"G4LundStringFragmentation::"<<G4endl;
169G4cout<<"FragmentString Position"<<theString.GetPosition()/fermi<<" "<<
170theString.GetTimeOfCreation()/fermi<<G4endl;
171G4cout<<"FragmentString Momentum"<<theString.Get4Momentum()<<theString.Get4Momentum().mag()<<G4endl;
172*/
173 G4FragmentingString aString(theString);
174 SetMinimalStringMass(&aString);
175
176 if((MinimalStringMass+WminLUND)*(1.-SmoothParam) > theString.Get4Momentum().mag())
177 {SetMassCut(1000.*MeV);}
178// V.U. 20.06.10 in order to put un correspondence LightFragTest and MinStrMass
179
180 G4KineticTrackVector * LeftVector=LightFragmentationTest(&theString);
181
182 if ( LeftVector != 0 ) {
183// Uzhi insert 6.05.08 start
184 if(LeftVector->size() == 1){
185 // One hadron is saved in the interaction
186 LeftVector->operator[](0)->SetFormationTime(theString.GetTimeOfCreation());
187 LeftVector->operator[](0)->SetPosition(theString.GetPosition());
188
189/* // To set large formation time open *
190 LeftVector->operator[](0)->SetFormationTime(theString.GetTimeOfCreation()+100.*fermi);
191 LeftVector->operator[](0)->SetPosition(theString.GetPosition());
192 G4ThreeVector aPosition(theString.GetPosition().x(),
193 theString.GetPosition().y(),
194 theString.GetPosition().z()+100.*fermi);
195 LeftVector->operator[](0)->SetPosition(aPosition);
196*/
197 } else { // 2 hadrons created from qq-qqbar are stored
198
199 LeftVector->operator[](0)->SetFormationTime(theString.GetTimeOfCreation());
200 LeftVector->operator[](0)->SetPosition(theString.GetPosition());
201 LeftVector->operator[](1)->SetFormationTime(theString.GetTimeOfCreation());
202 LeftVector->operator[](1)->SetPosition(theString.GetPosition());
203 }
204 return LeftVector;
205 }
206
207//--------------------- The string can fragment -------------------------------
208//--------------- At least two particles can be produced ----------------------
209 LeftVector =new G4KineticTrackVector;
210 G4KineticTrackVector * RightVector=new G4KineticTrackVector;
211
212 G4ExcitedString *theStringInCMS=CPExcited(theString);
213 G4LorentzRotation toCms=theStringInCMS->TransformToAlignedCms();
214
215 G4bool success=false, inner_sucess=true;
216 G4int attempt=0;
217 while ( !success && attempt++ < StringLoopInterrupt )
218 { // If the string fragmentation do not be happend, repeat the fragmentation---
219 G4FragmentingString *currentString=new G4FragmentingString(*theStringInCMS);
220
221 // Cleaning up the previously produced hadrons ------------------------------
222 std::for_each(LeftVector->begin() , LeftVector->end() , DeleteKineticTrack());
223 LeftVector->clear();
224
225 std::for_each(RightVector->begin(), RightVector->end(), DeleteKineticTrack());
226 RightVector->clear();
227
228 // Main fragmentation loop until the string will not be able to fragment ----
229 inner_sucess=true; // set false on failure..
230 while (! StopFragmenting(currentString) )
231 { // Split current string into hadron + new string
232 G4FragmentingString *newString=0; // used as output from SplitUp...
233
234 G4KineticTrack * Hadron=Splitup(currentString,newString);
235
236 if ( Hadron != 0 ) // Store the hadron
237 {
238 if ( currentString->GetDecayDirection() > 0 )
239 LeftVector->push_back(Hadron);
240 else
241 RightVector->push_back(Hadron);
242 delete currentString;
243 currentString=newString;
244 }
245 };
246 // Split remaining string into 2 final Hadrons ------------------------
247 if ( inner_sucess &&
248 SplitLast(currentString,LeftVector, RightVector) )
249 {
250 success=true;
251 }
252 delete currentString;
253 } // End of the loop in attemps to fragment the string
254
255 delete theStringInCMS;
256
257 if ( ! success )
258 {
259 std::for_each(LeftVector->begin(), LeftVector->end(), DeleteKineticTrack());
260 LeftVector->clear();
261 std::for_each(RightVector->begin(), RightVector->end(), DeleteKineticTrack());
262 delete RightVector;
263 return LeftVector;
264 }
265
266 // Join Left- and RightVector into LeftVector in correct order.
267 while(!RightVector->empty())
268 {
269 LeftVector->push_back(RightVector->back());
270 RightVector->erase(RightVector->end()-1);
271 }
272 delete RightVector;
273
274 CalculateHadronTimePosition(theString.Get4Momentum().mag(), LeftVector);
275
276 G4LorentzRotation toObserverFrame(toCms.inverse());
277
278// LeftVector->operator[](0)->SetFormationTime(theString.GetTimeOfCreation());
279// LeftVector->operator[](0)->SetPosition(theString.GetPosition());
280
281 G4double TimeOftheStringCreation=theString.GetTimeOfCreation();
282 G4ThreeVector PositionOftheStringCreation(theString.GetPosition());
283
284 for(size_t C1 = 0; C1 < LeftVector->size(); C1++)
285 {
286 G4KineticTrack* Hadron = LeftVector->operator[](C1);
287 G4LorentzVector Momentum = Hadron->Get4Momentum();
288 Momentum = toObserverFrame*Momentum;
289 Hadron->Set4Momentum(Momentum);
290
291 G4LorentzVector Coordinate(Hadron->GetPosition(), Hadron->GetFormationTime());
292 Momentum = toObserverFrame*Coordinate;
293 Hadron->SetFormationTime(TimeOftheStringCreation+Momentum.e()
294 -fermi/c_light);
295 G4ThreeVector aPosition(Momentum.vect());
296// Hadron->SetPosition(theString.GetPosition()+aPosition);
297 Hadron->SetPosition(PositionOftheStringCreation+aPosition);
298 };
299
300 return LeftVector;
301
302}
303
304//----------------------------------------------------------------------------------
305G4bool G4LundStringFragmentation::IsFragmentable(const G4FragmentingString * const string)
306{
307 SetMinimalStringMass(string);
308 return sqr(MinimalStringMass + WminLUND) < string->Get4Momentum().mag2();
309}
310
311//----------------------------------------------------------------------------------------
312G4bool G4LundStringFragmentation::StopFragmenting(const G4FragmentingString * const string)
313{
314 SetMinimalStringMass(string);
315
316//G4cout<<"StopFragm MinMass "<<MinimalStringMass<<" String Mass "<<std::sqrt(string->Get4Momentum().mag2())<<G4endl;
317//G4cout<<"WminLUND "<<WminLUND<<" SmoothParam "<<SmoothParam<<" "<<string->Mass()<<G4endl;
318 return (MinimalStringMass + WminLUND)*
319 (1 + SmoothParam * (1.-2*G4UniformRand())) >
320 string->Mass();
321}
322
323//----------------------------------------------------------------------------------------
324G4bool G4LundStringFragmentation::SplitLast(G4FragmentingString * string,
325 G4KineticTrackVector * LeftVector,
326 G4KineticTrackVector * RightVector)
327{
328 //... perform last cluster decay
329
330 G4LorentzVector Str4Mom=string->Get4Momentum();
331
332 G4ThreeVector ClusterVel =string->Get4Momentum().boostVector();
333
334 G4double ResidualMass = string->Mass();
335
336 G4ParticleDefinition * LeftHadron, * RightHadron;
337
338 G4int cClusterInterrupt = 0;
339 do
340 {
341 if (cClusterInterrupt++ >= ClusterLoopInterrupt)
342 {
343 return false;
344 }
345 G4ParticleDefinition * quark = NULL;
346 string->SetLeftPartonStable(); // to query quark contents..
347
348 if (!string->FourQuarkString() )
349 {
350 // The string is q-qbar, or q-qq, or qbar-qqbar type
351 if (string->DecayIsQuark() && string->StableIsQuark() )
352 {
353 //... there are quarks on cluster ends
354 LeftHadron= QuarkSplitup(string->GetLeftParton(), quark);
355 } else
356 {
357 //... there is a Diquark on one of the cluster ends
358 G4int IsParticle;
359
360 if ( string->StableIsQuark() )
361 {
362 IsParticle=(string->GetLeftParton()->GetPDGEncoding()>0) ? -1 : +1;
363 } else
364 {
365 IsParticle=(string->GetLeftParton()->GetPDGEncoding()>0) ? +1 : -1;
366 }
367
368 pDefPair QuarkPair = CreatePartonPair(IsParticle,false); // no diquarks wanted
369 quark = QuarkPair.second;
370
371 LeftHadron=hadronizer->Build(QuarkPair.first, string->GetLeftParton());
372 }
373
374 RightHadron = hadronizer->Build(string->GetRightParton(), quark);
375 } else
376 {
377 // The string is qq-qqbar type. Diquarks are on the string ends
378 G4int LiftQuark1= string->GetLeftParton()->GetPDGEncoding()/1000;
379 G4int LiftQuark2=(string->GetLeftParton()->GetPDGEncoding()/100)%10;
380
381 G4int RightQuark1= string->GetRightParton()->GetPDGEncoding()/1000;
382 G4int RightQuark2=(string->GetRightParton()->GetPDGEncoding()/100)%10;
383
384 if(G4UniformRand()<0.5)
385 {
386 LeftHadron =hadronizer->Build(FindParticle( LiftQuark1),
387 FindParticle(RightQuark1));
388 RightHadron=hadronizer->Build(FindParticle( LiftQuark2),
389 FindParticle(RightQuark2));
390 } else
391 {
392 LeftHadron =hadronizer->Build(FindParticle( LiftQuark1),
393 FindParticle(RightQuark2));
394 RightHadron=hadronizer->Build(FindParticle( LiftQuark2),
395 FindParticle(RightQuark1));
396 }
397 }
398
399 //... repeat procedure, if mass of cluster is too low to produce hadrons
400 //... ClusterMassCut = 0.15*GeV model parameter
401 }
402 while (ResidualMass <= LeftHadron->GetPDGMass() + RightHadron->GetPDGMass());
403
404 //... compute hadron momenta and energies
405
406 G4LorentzVector LeftMom, RightMom;
407 G4ThreeVector Pos;
408
409 Sample4Momentum(&LeftMom, LeftHadron->GetPDGMass(),
410 &RightMom, RightHadron->GetPDGMass(),
411 ResidualMass);
412
413 LeftMom.boost(ClusterVel);
414 RightMom.boost(ClusterVel);
415
416 LeftVector->push_back(new G4KineticTrack(LeftHadron, 0, Pos, LeftMom));
417 RightVector->push_back(new G4KineticTrack(RightHadron, 0, Pos, RightMom));
418
419 return true;
420
421}
422
423//----------------------------------------------------------------------------------------------------------
424
425void G4LundStringFragmentation::Sample4Momentum(G4LorentzVector* Mom, G4double Mass, G4LorentzVector* AntiMom, G4double AntiMass, G4double InitialMass)
426 {
427// ------ Sampling of momenta of 2 last produced hadrons --------------------
428 G4ThreeVector Pt;
429 G4double MassMt2, AntiMassMt2;
430 G4double AvailablePz, AvailablePz2;
431
432 if(Mass > 930. || AntiMass > 930.) // If there is a baryon
433 {
434 // ----------------- Isotripic decay ------------------------------------
435 G4double r_val = sqr(InitialMass*InitialMass - Mass*Mass - AntiMass*AntiMass) -
436 sqr(2.*Mass*AntiMass);
437 G4double Pabs = (r_val > 0.)? std::sqrt(r_val)/(2.*InitialMass) : 0;
438
439 //... sample unit vector
440 G4double pz = 1. - 2.*G4UniformRand();
441 G4double st = std::sqrt(1. - pz * pz)*Pabs;
442 G4double phi = 2.*pi*G4UniformRand();
443 G4double px = st*std::cos(phi);
444 G4double py = st*std::sin(phi);
445 pz *= Pabs;
446
447 Mom->setPx(px); Mom->setPy(py); Mom->setPz(pz);
448 Mom->setE(std::sqrt(Pabs*Pabs + Mass*Mass));
449
450 AntiMom->setPx(-px); AntiMom->setPy(-py); AntiMom->setPz(-pz);
451 AntiMom->setE (std::sqrt(Pabs*Pabs + AntiMass*AntiMass));
452 }
453 else
454 {
455 do
456 {
457 // GF 22-May-09, limit sampled pt to allowed range
458
459 G4double termD = InitialMass*InitialMass -Mass*Mass - AntiMass*AntiMass;
460 G4double termab = 4*sqr(Mass*AntiMass);
461 G4double termN = 2*termD + 4*Mass*Mass + 4*AntiMass*AntiMass;
462 G4double pt2max=(termD*termD - termab )/ termN ;
463
464 Pt=SampleQuarkPt(std::sqrt(pt2max)); Pt.setZ(0); G4double Pt2=Pt.mag2();
465
466 MassMt2 = Mass * Mass + Pt2;
467 AntiMassMt2= AntiMass * AntiMass + Pt2;
468
469 AvailablePz2= sqr(InitialMass*InitialMass - MassMt2 - AntiMassMt2) -
470 4.*MassMt2*AntiMassMt2;
471 }
472 while(AvailablePz2 < 0.); // GF will occur only for numerical precision problem with limit in sampled pt
473
474 AvailablePz2 /=(4.*InitialMass*InitialMass);
475 AvailablePz = std::sqrt(AvailablePz2);
476
477 G4double Px=Pt.getX();
478 G4double Py=Pt.getY();
479
480 Mom->setPx(Px); Mom->setPy(Py); Mom->setPz(AvailablePz);
481 Mom->setE(std::sqrt(MassMt2+AvailablePz2));
482
483 AntiMom->setPx(-Px); AntiMom->setPy(-Py); AntiMom->setPz(-AvailablePz);
484 AntiMom->setE (std::sqrt(AntiMassMt2+AvailablePz2));
485 }
486 }
487
488//-----------------------------------------------------------------------------
489
490G4LorentzVector * G4LundStringFragmentation::SplitEandP(G4ParticleDefinition * pHadron,
491 G4FragmentingString * string, G4FragmentingString * newString)
492{
493 G4LorentzVector String4Momentum=string->Get4Momentum();
494 G4double StringMT2=string->Get4Momentum().mt2();
495
496 G4double HadronMass = pHadron->GetPDGMass();
497
498 SetMinimalStringMass(newString);
499 String4Momentum.setPz(0.);
500 G4ThreeVector StringPt=String4Momentum.vect();
501
502// calculate and assign hadron transverse momentum component HadronPx and HadronPy
503 G4ThreeVector thePt;
504 thePt=SampleQuarkPt();
505
506 G4ThreeVector HadronPt = thePt +string->DecayPt();
507 HadronPt.setZ(0);
508
509 G4ThreeVector RemSysPt = StringPt - HadronPt;
510
511 //... sample z to define hadron longitudinal momentum and energy
512 //... but first check the available phase space
513
514 G4double HadronMassT2 = sqr(HadronMass) + HadronPt.mag2();
515 G4double ResidualMassT2=sqr(MinimalStringMass) + RemSysPt.mag2();
516
517 G4double Pz2 = (sqr(StringMT2 - HadronMassT2 - ResidualMassT2) -
518 4*HadronMassT2 * ResidualMassT2)/4./StringMT2;
519
520 if(Pz2 < 0 ) {return 0;} // have to start all over!
521
522 //... then compute allowed z region z_min <= z <= z_max
523
524 G4double Pz = std::sqrt(Pz2);
525 G4double zMin = (std::sqrt(HadronMassT2+Pz2) - Pz)/std::sqrt(StringMT2);
526 G4double zMax = (std::sqrt(HadronMassT2+Pz2) + Pz)/std::sqrt(StringMT2);
527
528 if (zMin >= zMax) return 0; // have to start all over!
529
530 G4double z = GetLightConeZ(zMin, zMax,
531 string->GetDecayParton()->GetPDGEncoding(), pHadron,
532 HadronPt.x(), HadronPt.y());
533
534 //... now compute hadron longitudinal momentum and energy
535 // longitudinal hadron momentum component HadronPz
536
537 HadronPt.setZ(0.5* string->GetDecayDirection() *
538 (z * string->LightConeDecay() -
539 HadronMassT2/(z * string->LightConeDecay())));
540
541 G4double HadronE = 0.5* (z * string->LightConeDecay() +
542 HadronMassT2/(z * string->LightConeDecay()));
543
544 G4LorentzVector * a4Momentum= new G4LorentzVector(HadronPt,HadronE);
545
546 return a4Momentum;
547}
548
549//-----------------------------------------------------------------------------------------
550G4double G4LundStringFragmentation::GetLightConeZ(G4double zmin, G4double zmax,
551 G4int PDGEncodingOfDecayParton,
552 G4ParticleDefinition* pHadron,
553 G4double Px, G4double Py)
554{
555 G4double alund;
556
557// If blund get restored, you MUST adapt the calculation of zOfMaxyf.
558// const G4double blund = 1;
559
560 G4double z, yf;
561 G4double Mass = pHadron->GetPDGMass();
562// G4int HadronEncoding=pHadron->GetPDGEncoding();
563
564 G4double Mt2 = Px*Px + Py*Py + Mass*Mass;
565
566 if(std::abs(PDGEncodingOfDecayParton) < 1000)
567 {
568 // ---------------- Quark fragmentation ----------------------
569 alund=0.35/GeV/GeV; // Instead of 0.7 because kinks are not considered
570 G4double zOfMaxyf=alund*Mt2/(alund*Mt2 + 1.);
571 G4double maxYf=(1-zOfMaxyf)/zOfMaxyf * std::exp(-alund*Mt2/zOfMaxyf);
572
573 do
574 {
575 z = zmin + G4UniformRand()*(zmax-zmin);
576// yf = std::pow(1. - z, blund)/z*std::exp(-alund*Mt2/z);
577 yf = (1-z)/z * std::exp(-alund*Mt2/z);
578 }
579 while (G4UniformRand()*maxYf > yf);
580 }
581 else
582 {
583 // ---------------- Di-quark fragmentation ----------------------
584 alund=0.7/GeV/GeV; // 0.7 2.0
585 G4double zOfMaxyf=alund*Mt2/(alund*Mt2 + 1.);
586 G4double maxYf=(1-zOfMaxyf)/zOfMaxyf * std::exp(-alund*Mt2/zOfMaxyf);
587 do
588 {
589 z = zmin + G4UniformRand()*(zmax-zmin);
590// yf = std::pow(1. - z, blund)/z*std::exp(-alund*Mt2/z);
591 yf = (1-z)/z * std::exp(-alund*Mt2/z);
592 }
593 while (G4UniformRand()*maxYf > yf);
594 };
595
596 return z;
597 }
Note: See TracBrowser for help on using the repository browser.