Context Navigation

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

Last change on this file since 1337 was 1337, checked in by garnier, 14 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
45	G4LundStringFragmentation::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	// --------------------------------------------------------------
65	G4LundStringFragmentation::G4LundStringFragmentation(const G4LundStringFragmentation &) : G4VLongitudinalStringDecay()
66	{
67	}
68
69	G4LundStringFragmentation::~G4LundStringFragmentation()
70	{
71	}
72
73	//*************************************************************************************
74
75	const G4LundStringFragmentation & G4LundStringFragmentation::operator=(const G4LundStringFragmentation &)
76	{
77	throw G4HadronicException(__FILE__, __LINE__, "G4LundStringFragmentation::operator= meant to not be accessable");
78	return *this;
79	}
80
81	int G4LundStringFragmentation::operator==(const G4LundStringFragmentation &right) const
82	{
83	return !memcmp(this, &right, sizeof(G4LundStringFragmentation));
84	}
85
86	int G4LundStringFragmentation::operator!=(const G4LundStringFragmentation &right) const
87	{
88	return memcmp(this, &right, sizeof(G4LundStringFragmentation));
89	}
90
91	//--------------------------------------------------------------------------------------
92	void 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	//--------------------------------------------------------------------------------------
150	void G4LundStringFragmentation::SetMinimalStringMass2(
151	const G4double aValue)
152	{
153	MinimalStringMass2=aValue * aValue;
154	}
155
156	//--------------------------------------------------------------------------------------
157	G4KineticTrackVector* 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	/*
168	G4cout<<G4endl<<"G4LundStringFragmentation::"<<G4endl;
169	G4cout<<"FragmentString Position"<<theString.GetPosition()/fermi<<" "<<
170	theString.GetTimeOfCreation()/fermi<<G4endl;
171	G4cout<<"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	//----------------------------------------------------------------------------------
305	G4bool G4LundStringFragmentation::IsFragmentable(const G4FragmentingString * const string)
306	{
307	SetMinimalStringMass(string);
308	return sqr(MinimalStringMass + WminLUND) < string->Get4Momentum().mag2();
309	}
310
311	//----------------------------------------------------------------------------------------
312	G4bool 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	//----------------------------------------------------------------------------------------
324	G4bool 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
425	void 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(InitialMassInitialMass - MassMass - AntiMass*AntiMass) -
436	sqr(2.MassAntiMass);
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.piG4UniformRand();
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(PabsPabs + MassMass));
449
450	AntiMom->setPx(-px); AntiMom->setPy(-py); AntiMom->setPz(-pz);
451	AntiMom->setE (std::sqrt(PabsPabs + AntiMassAntiMass));
452	}
453	else
454	{
455	do
456	{
457	// GF 22-May-09, limit sampled pt to allowed range
458
459	G4double termD = InitialMassInitialMass -MassMass - AntiMass*AntiMass;
460	G4double termab = 4sqr(MassAntiMass);
461	G4double termN = 2termD + 4MassMass + 4AntiMass*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.MassMt2AntiMassMt2;
471	}
472	while(AvailablePz2 < 0.); // GF will occur only for numerical precision problem with limit in sampled pt
473
474	AvailablePz2 /=(4.InitialMassInitialMass);
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
490	G4LorentzVector * 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	4HadronMassT2 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	//-----------------------------------------------------------------------------------------
550	G4double 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 = PxPx + PyPy + 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=alundMt2/(alundMt2 + 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)/zstd::exp(-alundMt2/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=alundMt2/(alundMt2 + 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)/zstd::exp(-alundMt2/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.

Download in other formats: