Context Navigation

source: trunk/source/processes/hadronic/models/chiral_inv_phase_space/body/src/G4QuasmonString.cc @ 1071

Last change on this file since 1071 was 1055, checked in by garnier, 15 years ago
maj sur la beta de geant 4.9.3
File size: 41.9 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	// 1 2 3 4 5 6 7 8 9
27	//34567890123456789012345678901234567890123456789012345678901234567890123456789012345678901
28	//
29	//
30	// $Id: G4QuasmonString.cc,v 1.8 2009/02/23 09:49:24 mkossov Exp $
31	// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
32	//
33	// ---------------- G4QuasmonString ----------------
34	// by Mikhail Kossov, August 2000.
35	// class for Hadron-Hadron String Interaction used by the CHIPS Model
36	// -------------------------------------------------------------------
37	// Short description: The Quark-Gluon String can fragment in hadrons
38	// (QGS), small hadronic clusters (HERWIG) or in unlimited in mass
39	// Quasmons, defined by the G4QPartonPair. This class is made to
40	// simulate the string fragmentation in Quasmons, which hadronize in
41	// vacuuum or in nuclear environment.
42	// -------------------------------------------------------------------
43
44	//#define chdebug
45	//#define debug
46	//#define sdebug
47	//#define ppdebug
48	//#define cdebug
49	//#define cldebug
50	//#define edebug
51	//#define fdebug
52	//#define pdebug
53	//#define rdebug
54	//#define ffdebug
55	//#define pcdebug
56	//#define mudebug
57
58	#include "G4QuasmonString.hh"
59	#include <cmath>
60	using namespace std;
61
62
63	G4QuasmonString::G4QuasmonString(G4QHadron projHadron, const G4bool projEnvFlag,
64	const G4int targPDG, const G4bool targEnvFlag) :
65	theProjEnvFlag(projEnvFlag), theTargEnvFlag(targEnvFlag), theWeight(1.)
66	{
67	//static const G4double mPi0 = G4QPDGCode(111).GetMass(); // Pi0 mass
68	//static const G4double mPi0 = (1232.-938.272)/2; // between p&Delta (147.365) shift
69	static const G4double mPi0 = sqrt((1232.1232+938.272938.272)/2)-938.272; // Pi0 mass
70	//static const G4double mPi0 = 170.; // Pi0 mass
71	static const G4double FreeFra = 1.; // @@ Free fraction for n-n scattering
72	static const G4double interc = -1.; // @@ Pomeron intercept
73	//static const G4double Temperature = 180.; // Temperature as a parameter
74	//static const G4double Temp2 = Temperature*Temperature; // Squared Temperature
75	static G4Quasmon classG4Quasmon(G4QContent(1,2,0,0,0,0));
76	//static const G4QPDGCode pimQPDG(-211);
77	theWorld= G4QCHIPSWorld::Get(); // Get a pointer to the CHIPS World
78	#ifdef debug
79	G4cout<<"-->G4QuasStr::Constr:pPDG="<<projHadron.GetPDGCode()<<",tPDG="<<targPDG<<G4endl;
80	#endif
81	// Target initialization
82	G4QPDGCode targQPDG(targPDG);
83	totBaryNum= targQPDG.GetBaryNum(); // Prototype
84	totCharge = targQPDG.GetCharge(); // Prototype
85	theTargQC = targQPDG.GetQuarkContent();
86	G4double tM = targQPDG.GetMass();
87	theTarg4Mom = G4LorentzVector(0.,0.,0.,tM);
88	G4QHadron targHadron(targQPDG,theTarg4Mom); // Target Hadron
89
90	// Projectile initialization
91	theProj4Mom = projHadron.Get4Momentum();
92	G4ThreeVector projBoost = theProj4Mom.boostVector(); // Projectile BoostVector to LS
93	G4ThreeVector projRBoost= -projBoost; // Projevtile Boost vector to projectile CMS
94	G4QPDGCode projQPDG(projHadron.GetPDGCode());
95	theProjQC = projQPDG.GetQuarkContent();
96	totBaryNum+= projQPDG.GetBaryNum(); // Final control value
97	totCharge += projQPDG.GetCharge(); // Final control value
98	tot4Mom = theProj4Mom+theTarg4Mom; // Final control value
99
100	G4double projM = projQPDG.GetMass();
101	G4double projM2 = projM*projM;
102	G4double pM2 = theProj4Mom.m2();
103	G4double tM2 = tM*tM;
104	G4double dM2 = fabs(pM2-projM2);
105	if(dM2>1.) G4cout<<"-Wor-G4QuStr::Cons:dM2="<<dM2<<",M2="<<projM2<<",LVM2="<<pM2<<G4endl;
106	G4double pM=sqrt(pM2); // @@ do we need pM ? @@ (in print)
107
108	// === Print out of the input information at Creation time & tot 4-mom Calculation ======
109	#ifdef pdebug
110	G4cout<<"G4QuaString::Cons:PQC="<<theProjQC<<",TQC="<<theTargQC<<",P4Mom="<<theProj4Mom<<
111	theProj4Mom.m2()<<theProj4Mom.m()<<G4endl;
112	G4cout<<"G4QuaString::C: tC="<<totCharge<<",tB="<<totBaryNum<<",tot4M="<<tot4Mom<<G4endl;
113	#endif
114	//G4int nP=theWorld->GetQPEntries(); // A#of init'ed particles in CHIPS World (@@?)
115	//G4int nCl=nP-90; // A#of init'ed clusters in CHIPS World (@@?)
116	//#ifdef pdebug
117	//G4cout<<"G4QS:Const:Before QEX:n="<<nP<<G4endl;
118	//#endif
119	// @@@@@@ ===> Here the Quark Exchange Quasmon Creation must be added <=== @@@@@@
120	// @@ --- old ---
121	//G4int nQTarg=theTargQC.GetTot(); // a#of Quark-partons in the Target
122	//G4int nQProj=theProjQC.GetTot(); // a#of Quark-partons in the Progectile
123	// @@ --- new ---
124	G4LorentzVector tSum=theProj4Mom+theTarg4Mom;
125	G4double sumM=tM+pM;
126	G4double dn=2.17*(tSum.m2()/sumM/sumM-1.); // additional partons (gluons) - real
127	G4int np=static_cast<int>(dn); // additional partons (gluons) - basic int
128	if(G4UniformRand()>1.-dn+np) np++; // randomized int additional partons(gluons)
129	G4int nQTarg=theTargQC.GetTot()+np; // a#of Quark-partons in the Target
130	G4int nQProj=theProjQC.GetTot()+np; // a#of Quark-partons in the Progectile
131	// @@ --- new ---
132	//G4int mQTarg=theTargQC.GetTot(); // a#of Quark-partons in the Target
133	//G4int mQProj=theProjQC.GetTot(); // a#of Quark-partons in the Progectile
134	////G4double tothM=tot4Mom.m()/2.; // C.M. mass of the protons
135	////G4int nQCM=classG4Quasmon.CalculateNumberOfQPartons(tothM);
136	//G4double tothM2=tot4Mom.m2()/4.; // C.M. mass of the protons
137	//G4int nQCM=static_cast<int>((1.+sqrt(tothM2/Temp2+1.))/2.);
138	//G4int nQTarg=nQCM; // a#of Quark-partons in the Target
139	//if(nQTarg<mQTarg) nQTarg=mQTarg;
140	//G4int nQProj=nQCM; // a#of Quark-partons in the Progectile
141	//if(nQProj<mQProj) nQProj=mQProj;
142	// @@ --- end ---
143	#ifdef pdebug
144	G4cout<<"G4QuasmonString::Constr: nP="<<nQProj<<", nT="<<nQTarg<<G4endl;
145	#endif
146	// @@ Now projectile can be only meson or baryon @@ -- @@ Improve for clusters @@ --
147	G4LorentzVector pq4Mom(0.,0.,0.,0.); // Prototype of LV of quark of progectile
148	G4double rPMass=0.; // Prototype of the residProjMass (Meson case)
149	G4bool FreeFraF=false; // Prototype of the free exchange Flag
150	if(G4UniformRand()<FreeFra) FreeFraF=true; // Confirm the free exchange
151	if(nQProj<2) G4cout<<"*G4QuaStr::Con: nQProj="<<nQProj<<"<2 FatalError**"<<G4endl;
152	else if(nQProj>2) // ---> Baryon case (clusters are not implem.)
153	{
154	//if(nQProj>3)G4cout<<"-Wor-G4QS::Const:nQProj="<<nQProj<<">3 is not implem'd"<<G4endl;
155	G4double xP=0.;
156	if(FreeFraF) xP=RandomizeMomFractionFree(nQProj);
157	else
158	{
159	xP=RandomizeMomFractionString(nQProj);
160	theWeight*=pow(xP,interc);
161	G4cout<<"************G4QuasmonString::Constr: string xP="<<xP<<G4endl;
162	}
163	rPMass = sqrt(pM2*(1.-xP)); // Residual Projectile mass
164	#ifdef pdebug
165	G4cout<<"G4QS::C: nQProj="<<nQProj<<", xProj="<<xP<<", rPMass="<<rPMass<<G4endl;
166	#endif
167	}
168	G4LorentzVector pr4Mom(0.,0.,0.,rPMass); // Prototype of LV of the residual projectile
169	G4bool projFl=false;
170	G4bool targFl=false;
171	G4bool tmpBl=projHadron.DecayIn2(pq4Mom,pr4Mom);
172	if(tmpBl) projFl=true; // Proj decay is OK
173	//if(projHadron.DecayIn2(pq4Mom,pr4Mom)) projFl=true;
174	else G4cout<<"WarningG4QuasmonString::C:ProjDecIn2 rpM="<<rPMass<<",pM="<<pM<<G4endl;
175	#ifdef pdebug
176	G4cout<<"G4QuaStr::C:"<<projFl<<" split PROJ in R4M="<<pr4Mom<<" & Q4M="<<pq4Mom<<G4endl;
177	#endif
178	G4LorentzVector tq4Mom(0.,0.,0.,0.); // Prototype of LV of quark of the target
179	//if(nQTarg<3)G4cout<<"*G4QStr::Const: nQTarg="<<nQTarg<<"<3 FatalError**"<<G4endl;
180	//if(nQTarg>3)G4cout<<"-W-G4QS::Const: nQTarg="<<nQTarg<<">3 is not implemented"<<G4endl;
181	G4double xT=0.;
182	if(FreeFraF) xT=RandomizeMomFractionFree(nQTarg);
183	else
184	{
185	xT=RandomizeMomFractionString(nQTarg);
186	theWeight*=pow(xT,interc);
187	G4cout<<"************G4QuaString::Constr: string xT="<<xT<<G4endl;
188	}
189	G4double rTMass = sqrt(tM2*(1.-xT)); // Residual Target mass
190	#ifdef pdebug
191	G4cout<<"G4QuaString::C: nQTarg="<<nQTarg<<", xProj="<<xT<<", rTMass="<<rTMass<<G4endl;
192	#endif
193	G4LorentzVector tr4Mom(0.,0.,0.,rTMass); // Prototype of LV of the residual projectile
194	if(targHadron.DecayIn2(tq4Mom,tr4Mom)) targFl=true; // Targ decay is OK
195	else G4cout<<"WorningG4QuaStr::Constr:TargDecIn2 rtM="<<rTMass<<", tM="<<tM<<G4endl;
196	#ifdef pdebug
197	G4cout<<"G4QuaStr::C:"<<targFl<<" split TARG in R4M="<<tr4Mom<<" & Q4M="<<tq4Mom<<G4endl;
198	#endif
199	G4bool elasFl=false; // ByDefault avoid the elastic scattering
200	if (targFl && projFl) // --- @@ Now only for pp case @@ ---
201	{
202	G4LorentzVector newProj4M=pr4Mom+tq4Mom;
203	G4double nPM2=newProj4M.m2();
204	G4double npM=sqrt(nPM2);
205	G4bool pcorFl=false; // ByDefault no MassSh correction for newProj.
206	// @@ Just to try ---- Start
207	//G4LorentzVector newTarg4M=tr4Mom+pq4Mom;
208	//G4double nTM2=newTarg4M.m2();
209	//G4double ntM=sqrt(nTM2);
210	#ifdef pdebug
211	//G4cout<<"G4QStr::C:ntM="<<ntM<<" <? tM="<<tM<<"+mPi0="<<mPi0<<" = "<<tM+mPi0<<G4endl;
212	#endif
213	//if(ntM<tM+mPi0 && npM<pM+mPi0) elasFl=true; // Target&Project are underMinMass => El
214	// @@ Just to try ---- End
215	#ifdef pdebug
216	G4cout<<"G4QuaStr::C:npM="<<npM<<" <? pM="<<pM<<"+mPi0="<<mPi0<<" = "<<pM+mPi0<<G4endl;
217	#endif
218	if(npM<pM+mPi0) // The projectile is under min mass (@@ In Future put other cut @@)
219	{
220	G4double valk=tq4Mom.e(); // Momentum of the target quark
221	//#ifdef pdebug
222	G4double dvalk=valk-tq4Mom.rho();
223	if(fabs(dvalk)>.00001) G4cout<<"kpG4QuaStr::C:vk="<<valk<<",dvk="<<dvalk<<G4endl;
224	//#endif
225	G4double dmas=pM2-pr4Mom.m2(); // Difference of squared masses
226	G4double vale=pr4Mom.e();
227	G4ThreeVector pr=pr4Mom.v();
228	G4double valp=pr.mag();
229	G4ThreeVector upr=pr/valp; // Unit vector in the pr direction
230	G4double cost=-(dmas/(valk+valk)-vale)/valp; // Initial cos(theta)
231	if(fabs(cost)>1.)
232	{
233	#ifdef pdebug
234	G4cout<<"*p*>>>G4QuasmonString::Constr: cost="<<cost<<G4endl;
235	#endif
236	// Get max value of \|cos(theta)\| and change tq value to get the pM on mass shell
237	if(cost>0.) // --->>> cost>0.
238	{
239	valk=dmas/(vale-valp)/2; // Momentum is too big (must be reduced)
240	G4ThreeVector tqf=upr*valk; // final targetQuark 3-vector
241	tq4Mom.set(valk,tqf); // Fill new 4-mom for targetQuark
242	tr4Mom.set(tM-valk,-tqf); // Fill new 4-mom for targetResidual
243	newProj4M=tq4Mom+pr4Mom; // Final Projectile 4-mom in LS
244	}
245	else // --->>> cost<0.
246	{
247	G4double hvalk=dmas/(vale+valp); // Momentum's too small (must be increased, LIM)
248	if(hvalk>tM) // Momentum can not be increased to this value
249	{
250	#ifdef pdebug
251	G4cout<<"p-Cor>G4QuaStr::C: hvalk="<<hvalk<<" > tM="<<tM<<", dm="<<dmas
252	<<", e="<<vale<<", p="<<valp<<", ct="<<cost<<G4endl;
253	#endif
254	// Put the scatteredProjectile on the massShell, and rescatter the target quark
255	G4LorentzVector tmpProj4M=newProj4M+pq4Mom; // TempCriticalCompound for Project
256	newProj4M=G4LorentzVector(0.,0.,0.,pM); // New Project is on the mass shell
257	G4QHadron tmpHadron(tmpProj4M); // Create a Hadron for decay
258	pq4Mom=G4LorentzVector(0.,0.,0.,0.); // NewProjParton on lightMassShell
259	tmpBl=tmpHadron.DecayIn2(newProj4M,pq4Mom); // Decay the Compound in newProg+pQ
260	if(!tmpBl)
261	G4cout<<"G4QuaStr::Con: DecIn2 err "<<sqrt(tmpProj4M.m2())<<"<"<<pM<<G4endl;
262	}
263	else
264	{
265	#ifdef pdebug
266	G4cout<<"*---*>G4QuasmonStr::Constr: hvalk="<<hvalk<<" < tM="<<tM<<G4endl;
267	#endif
268	valk=hvalk/2; // Momentum is too small (must be reduced)
269	G4ThreeVector tqf=upr*(-valk); // Final targetQuark 3-vector
270	tq4Mom.set(valk,tqf); // Fill new 4-mom for targetQuark
271	tr4Mom.set(tM-valk,-tqf); // Fill new 4-mom for targetResidual
272	newProj4M=tq4Mom+pr4Mom;
273	}
274	}
275	if(fabs(newProj4M.m()-pM)>.0001)
276	G4cout<<"*G4QuaString::Constr: Change of mass "<<newProj4M.m()<<"#"<<pM<<G4endl;
277	}
278	else
279	{
280	// Turn the target quark-parton to the projectile residual in LS = CMSofTheTarget
281	#ifdef pdebug
282	G4cout<<"---p--->>>G4QuasmonString::Constr: cost="<<cost<<G4endl;
283	#endif
284	G4ThreeVector tq=tq4Mom.v(); // 3-mom of the targetQ (LS)
285	G4double tqlp=upr.dot(tq); // Projection of tq on the projectile
286	G4ThreeVector tql=upr*tqlp; // tq part along pr
287	G4ThreeVector tqt=tq-tql; // tq part perpendicular to pr
288	G4double cosi=tqlp/valk; // Initial cos(theta)
289	G4ThreeVector tqlf=tql*(cost/cosi); // final tq part along pr
290	G4double sini=sqrt(1.-cosi*cosi); // Initial sin(theta)
291	G4double sint=sqrt(1.-cost*cost); // Desired sin(theta)
292	G4ThreeVector tqpf=tqt*(sint/sini); // final tq part perpendicular pr
293	G4ThreeVector tqf=tqlf+tqpf; // final tq
294	tq4Mom.setV(tqf); // Change the momentum direction of targetQ
295	tr4Mom.setV(-tqf); // Change the momentum direction of targetR
296	if(fabs(tqf.mag()-valk)>.0001)
297	G4cout<<"*G4QuasmonString::Constr: \|k\|="<<tqf.mag()<<" # "<<valk<<G4endl;
298	newProj4M=tq4Mom+pr4Mom;
299	if(fabs(newProj4M.m()-pM)>.001)
300	G4cout<<"*G4QuasmonString::Constr: M="<<newProj4M.m()<<" # "<<pM<<G4endl;
301	}
302	#ifdef pdebug
303	G4cout<<"G4QuaStr::C: Proj under GS, newP4M="<<newProj4M<<", pq4M="<<pq4Mom<<G4endl;
304	#endif
305	pcorFl=true; // Projectile is on the GS mass shell
306	}
307	G4bool tcorFl=false; // ByDefault no MassSh correction for newTarg.
308	G4LorentzVector newTarg4M=tr4Mom+pq4Mom;
309	G4double nTM2=newTarg4M.m2();
310	G4double ntM=sqrt(nTM2);
311	//newTarg4M=tr4Mom+pq4Mom;
312	//nTM2=newTarg4M.m2();
313	//ntM=sqrt(nTM2);
314	#ifdef pdebug
315	G4cout<<"G4QuaStr::C:ntM="<<ntM<<" <? tM="<<tM<<"+mPi0="<<mPi0<<" = "<<tM+mPi0<<G4endl;
316	#endif
317	if(ntM<tM+mPi0 && !pcorFl) // The target is under min mass (@@ InFut put another cut@@)
318	{
319	// Turn the projectile quark-parton to the target rsidual in CMS of the Projectile
320	G4LorentzVector pqc4M=pq4Mom; // projectileQuark => projCM system <step1>
321	pqc4M.boost(projRBoost); // projectileQuark => projCM system <step2>
322	G4double valk=pqc4M.e(); // Momentum of the target quark in projCM
323	//#ifdef pdebug
324	G4double dvalk=valk-pqc4M.rho();
325	if(fabs(dvalk)>.00001) G4cout<<"*kt*G4QuaS::C:vk="<<valk<<",dvk="<<dvalk<<G4endl;
326	//#endif
327	G4double dmas=tM2-tr4Mom.m2(); // Difference of squared masses (targ - targRes)
328	G4LorentzVector trc4M=tr4Mom; // targetResidual => projCM system <step1>
329	trc4M.boost(projRBoost); // targetResidual => projCM system <step2>
330	G4double vale=trc4M.e(); // Energy of the targetResidual in projCM
331	G4ThreeVector tr=trc4M.v(); // 3-mom of the targetResidual in projCM
332	G4double valp=tr.mag(); // momentum of the targetResidual in projCM
333	if(fabs(dmas-tM2+trc4M.m2())>.1) G4cout<<"tG4QuaStr::C: trM2="<<tr4Mom.m2()<<"="
334	<<trc4M.m2()<<","<<valevale-valpvalp<<G4endl;
335	G4ThreeVector utr=tr/valp; // Unit vector in the tr direction in projCM
336	G4double cost=-(dmas/(valk+valk)-vale)/valp; // Initial cos(theta)
337	if(fabs(cost)>1.)
338	{
339	#ifdef pdebug
340	G4cout<<"*t*>>>G4QuasmonString::Constr: cost="<<cost<<G4endl;
341	#endif
342	// Get max value of \|cos(theta)\| and change pq value to get the tM on mass shell
343	if(cost>0.) // --->>> cost>0.
344	{
345	valk=dmas/(vale-valp)/2; // Momentum is too big (must be reduced)
346	G4ThreeVector pqf=utr*valk; // final projectileQuark 3-vector
347	G4LorentzVector pqc4M(valk,pqf); // Fill new 4-mom for projectileQuark in CM
348	pq4Mom=pqc4M; // <step1> // Fill new 4-mom for projectileQuark in LS
349	pq4Mom.boost(projBoost); // <step2> // Fill new 4-mom for projectileQuark in LS
350	G4LorentzVector prc4M(pM-valk,-pqf); // Fill new 4-mom for projectResidual in CM
351	pr4Mom=prc4M; // <step1> // Fill new 4-mom for projectResidual in LS
352	pr4Mom.boost(projBoost); // <step2> // Fill new 4-mom for projectResidual in LS
353	newTarg4M=pq4Mom+tr4Mom; // Final Target 4-mom in LS
354	}
355	else // --->>> cost<-1 => cost=-1
356	{
357	G4double hvalk=dmas/(vale+valp); // Momentum's too small (must be increased, LIM)
358	if(hvalk>pM) // Momentum can not be increased to this value
359	{
360	#ifdef pdebug
361	G4cout<<"t-Cor>G4QuasmonString::C: hvalk="<<hvalk<<" > pM="<<pM<<", dm="
362	<<dmas<<", e="<<vale<<", p="<<valp<<", ct="<<cost<<G4endl;
363	#endif
364	// Put the scatteredProjectile on the massShell, rescatter the targetQuark (LS)
365	G4LorentzVector tmpTarg4M=newTarg4M+tq4Mom; // TempCriticalCompound for Target
366	newTarg4M=G4LorentzVector(0.,0.,0.,tM); // New Target is on the mass shell
367	G4QHadron tmpHadron(tmpTarg4M); // Create a CompHadron for decay
368	tq4Mom=G4LorentzVector(0.,0.,0.,0.); // NewTargParton on lightMassShell
369	tmpBl=tmpHadron.DecayIn2(newTarg4M,tq4Mom); // Decay the Compound in newTarg+tQ
370	if(!tmpBl)
371	G4cout<<"G4QuaString::C: DecIn2-err "<<sqrt(tmpTarg4M.m2())<<"<"<<tM<<G4endl;
372	}
373	else
374	{
375	#ifdef pdebug
376	G4cout<<"*---*>G4QuaString::Constr: hvalk="<<hvalk<<" < pM="<<pM<<G4endl;
377	#endif
378	valk=hvalk/2; // Momentum is too small (mustBeIncreased)
379	G4ThreeVector pqf=utr*(-valk); // Final projectileQuark 3-vector in CM
380	G4LorentzVector pqc4M(valk,pqf); // Fill new 4-mom for projectQuark in CM
381	pq4Mom=pqc4M; // <step1> // Fill new 4-mom for projectQuark in LS
382	pq4Mom.boost(projBoost); // <step1> // Fill new 4-mom for projectQuark in LS
383	G4LorentzVector nTc4M=pqc4M+trc4M;
384	#ifdef pdebug
385	G4cout<<"G4QuasmonString::Constr:After Boost="<<projBoost<<G4endl;
386	G4cout<<"G4QuS::C:Aft:q="<<pqc4M<<pqc4M.m2()<<",t="<<nTc4M<<nTc4M.m2()<<G4endl;
387	G4cout<<"G4QuaStr::C:E="<<nTc4M.e()<<" = q="<<pqc4M.e()<<"="<<valk<<" + r="<<
388	trc4M.e()<<"="<<vale<<" = "<<pqc4M.e()+trc4M.e()<<",c="<<pq4Mom<<G4endl;
389	if(fabs(nTc4M.m()-tM)>.0001) G4cout<<"*G4QuS::C:"<<nTc4M.m()<<"!="<<tM<<G4endl;
390	G4double pp=pqf.dot(tr); // trc3Mpqc=-valkvalp
391	G4double ee=vale*valk;
392	G4cout<<"G4QS::C:tM2="<<tM2<<"="<<tMtM<<",pp="<<pp2<<"="<<-valkvalp2<<",d="
393	<<ee-pp<<"="<<dmas<<"="<<(tM2-trc4M.m2())<<",u="<<utr.dot(utr)<<G4endl;
394	G4double sen=nTc4M.e();
395	G4double smo=nTc4M.rho();
396	G4cout<<"G4QuaStr::C: qM2="<<pqc4M.m2()<<",rM2="<<trc4M.m2()<<",E="<<sen<<"="<<
397	vale+valk<<", P="<<smo<<"="<<valp-valk<<",M="<<sqrt(sensen-smosmo)<<G4endl;
398	#endif
399	G4LorentzVector prc4M(pM-valk,-pqf); // Fill new 4-mom for projectResid in CM
400	pr4Mom=prc4M; // <step1> // Fill new 4-mom for projectResid in LS
401	pr4Mom.boost(projBoost); // <step2> // Fill new 4-mom for projectResid in LS
402	newTarg4M=pq4Mom+tr4Mom;
403	}
404	}
405	if(fabs(newTarg4M.m()-tM)>.0001)
406	G4cout<<"****************G4QuasmonString::C: M="<<newTarg4M.m()<<"="<<pM<<G4endl;
407	}
408	else // -1<cos(theta)<1 - only rotate the projQ
409	{
410	// Turn the projectile quark-parton to the target residual in CMS of the Projectile
411	#ifdef pdebug
412	G4cout<<"---t--->>>G4QuasmonString::Constr: cost="<<cost<<G4endl;
413	#endif
414	G4LorentzVector prc4M=pr4Mom; // projectileQuark => projCM system <step1>
415	prc4M.boost(projRBoost); // projectileQuark => projCM system <step2>
416	G4ThreeVector pq=pqc4M.v(); // 3-vector of the projQ in projCM
417	G4double pqlp=utr.dot(pq); // Projection of pq on the target in projCM
418	G4ThreeVector pql=utr*pqlp; // pq part along pr
419	G4ThreeVector pqt=pq-pql; // pq part perpendicular to tr in projCM
420	G4double cosi=pqlp/valk; // Initial cos(theta)
421	G4ThreeVector pqlf=pql*(cost/cosi); // final pq part along tr
422	G4double sini=sqrt(1.-cosi*cosi); // Initial sin(theta)
423	G4double sint=sqrt(1.-cost*cost); // Desired sin(theta)
424	G4ThreeVector pqpf=pqt*(sint/sini); // final pq part perpendicular tr
425	G4ThreeVector pqf=pqlf+pqpf; // final pq
426	pqc4M.setV(pqf); // Change the momentumDirection of projQ(CM)
427	pq4Mom=pqc4M; // Fill new 4-mom for projQuark in LS<step1>
428	pq4Mom.boost(projBoost); // Fill new 4-mom for projQuark in LS<step2>
429	prc4M.setV(-pqf); // Change the momentumDirection of projR(CM)
430	#ifdef pdebug
431	G4LorentzVector nT4M=pqc4M+trc4M;
432	if(fabs(nT4M.m()-tM)>.001)G4cout<<"tG4QuaStr::C:M="<<nT4M.m()<<"="<<tM<<G4endl;
433	#endif
434	pr4Mom=prc4M; // <step1> // Fill new 4-mom for projResidual in LS
435	pr4Mom.boost(projBoost); // <step2> // Fill new 4-mom for projResidual in LS
436	if(fabs(pqf.mag()-valk)>.0001)
437	G4cout<<"*G4QuasmonString::Constr: \|k\|="<<pqf.mag()<<"="<<valk<<G4endl;
438	newTarg4M=pq4Mom+tr4Mom;
439	if(fabs(newTarg4M.m()-tM)>.001)
440	G4cout<<"*G4QuasmonString::Constr: M="<<newTarg4M.m()<<"="<<tM<<G4endl;
441	}
442	#ifdef pdebug
443	G4cout<<"*G4QuaStr::C: Targ under GS, newT4M="<<newTarg4M<<", tq4M="<<tq4Mom<<G4endl;
444	#endif
445	tcorFl=true; // Target is on the GS mass shell
446	newProj4M=pr4Mom+tq4Mom; // Recalculate the Projectile (!)
447	nPM2=newProj4M.m2();
448	npM=sqrt(nPM2);
449	#ifdef pdebug
450	G4cout<<"G4QuaStr::C:npM="<<npM<<" <? pM="<<pM<<"+mPi0="<<mPi0<<"="<<pM+mPi0<<G4endl;
451	#endif
452	if(npM<pM+mPi0) elasFl=true; // Force elastic scattering (@@InFut putAnotherCut@@)
453	}
454	else if(ntM<tM+mPi0) elasFl=true; // Force elastScattering (@@InFut put another cut@@)
455	if(elasFl) // Ellastic scattering happened
456	{
457	// ** Put the hadrons on the mass shell conserving the CMS scattering angle **
458	G4LorentzVector theTot4M=theProj4Mom+theTarg4Mom;// 4-momentum of CMS of "targ+proj"
459	G4ThreeVector cmsBoost = theTot4M.boostVector(); // CMS Boost Vector "CMS to LS"
460	G4ThreeVector cmsRBoost= -cmsBoost; // CMS Boost Vector "LS to CMS"
461	G4LorentzVector cmsProj4M=theProj4Mom; // LS projectile => CMS <step1>
462	cmsProj4M.boost(cmsRBoost); // LS projectile => CMS <step2>
463	G4LorentzVector cmsTarg4M=theTarg4Mom; // LS target => CMS <step1>
464	cmsTarg4M.boost(cmsRBoost); // LS target => CMS <step2>
465	G4double pcm=cmsProj4M.rho(); // CMS momentum for the elastic scattering
466	//#ifdef pdebug
467	if(fabs(cmsTarg4M.rho()-pcm) > 0.0001)
468	G4cout<<"-Worning-G4QuasmonString::Constr: P="<<cmsTarg4M.rho()<<"#"<<pcm<<G4endl;
469	//#endif
470	G4LorentzVector cmsNewPr4M=newProj4M; // LS finalProj => CMS <step1>
471	cmsNewPr4M.boost(cmsRBoost); // LS finalProj => CMS <step2>
472	G4ThreeVector puV=cmsNewPr4M.v()/cmsNewPr4M.rho(); // Direction of the projectile
473	//#ifdef pdebug
474	G4LorentzVector cmsNewTg4M=newTarg4M; // LS finalTarg => CMS <step1> @@ TMP
475	cmsNewTg4M.boost(cmsRBoost); // LS finalTarg => CMS <step2> @@ TMP
476	G4ThreeVector tuV=cmsNewTg4M.v()/cmsNewTg4M.rho(); // Direction of the projectile
477	if(1.+puV.dot(tuV) > 0.001)
478	G4cout<<"-Worning-G4QuasmonString::Constr: ct="<<puV.dot(tuV)<<G4endl;
479	//#endif
480	cmsProj4M.setV(puV*pcm);
481	newProj4M=cmsProj4M;
482	newProj4M.boost(cmsBoost); // CMS FinalProjectile => LS <step2>
483	G4QHadron* projH = new G4QHadron(projHadron); // Prototype of the Projectile Hadron
484	projH->Set4Momentum(newProj4M);
485	theQHadrons.push_back(projH); // Fill theProjectileHadron(delete equivalent)
486	#ifdef pdebug
487	G4cout<<"G4QuasmonString::Constr: Fill ElastProjH="<<projQPDG<<newProj4M<<G4endl;
488	#endif
489	cmsTarg4M.setV(puV*(-pcm));
490	newTarg4M=cmsTarg4M;
491	newTarg4M.boost(cmsBoost); // CMS FinalTarget => LS <step2>
492	G4QHadron* targH = new G4QHadron(targQPDG,newTarg4M); // Prototype of theTargetHadron
493	theQHadrons.push_back(targH); // Fill the Target Hadron (delete equivalent)
494	#ifdef pdebug
495	G4cout<<"G4QuasmonString::Constr: Fill ElastTargH="<<targQPDG<<newTarg4M<<G4endl;
496	#endif
497	}
498	else
499	{
500	// Inelastic scattering: one or both hadrons are excited (charge exchange is not in).
501	if(pcorFl) // Projectile is on the mass shell
502	{
503	G4QHadron* projH = new G4QHadron(projHadron); // Prototype of the Projectile Hadron
504	projH->Set4Momentum(newProj4M);
505	theQHadrons.push_back(projH); // Fill theProjectileHadron(delete equivalent)
506	G4Quasmon* targQ = new G4Quasmon(targQPDG.GetQuarkContent(),newTarg4M);
507	theQuasmons.push_back(targQ); // Insert Projectile Quasmon (delete equival.)
508	}
509	if(tcorFl) // Target is on the mass shell
510	{
511	G4Quasmon* projQ = new G4Quasmon(projHadron.GetQC(),newProj4M);
512	theQuasmons.push_back(projQ); // Insert Projectile Quasmon (delete equival.)
513	G4QHadron* targH = new G4QHadron(targQPDG,newTarg4M);//Prototype of theTargetHadron
514	theQHadrons.push_back(targH); // Fill the Target Hadron (delete equivalent)
515	}
516	else // Both are excited (two Quasmons only)
517	{
518	G4Quasmon* projQ = new G4Quasmon(projHadron.GetQC(),newProj4M);
519	theQuasmons.push_back(projQ); // Insert Projectile Quasmon (delete equival.)
520	G4Quasmon* targQ = new G4Quasmon(targQPDG.GetQuarkContent(),newTarg4M);
521	theQuasmons.push_back(targQ); // Insert Projectile Quasmon (delete equival.)
522	}
523	}
524	}
525	else G4cout<<"-Wor-G4QuasmonString::Constr:T="<<targFl<<" or P="<<projFl<<" err"<<G4endl;
526	#ifdef pdebug
527	G4cout<<"G4QuasmonString::Constructor: * End of Constructor *, elF="<<elasFl<<G4endl;
528	#endif
529	// This is the first Check of the control values -- @@ Must be remade @@ --
530	#ifdef chdebug
531	G4int finCharge=0;
532	G4int finBaryoN=0;
533	G4int nHad=theQHadrons.size();
534	if(nHad) for(G4int ih=0; ih<nHad; ih++)
535	{
536	finCharge+=theQHadrons[ih]->GetCharge();
537	finBaryoN+=theQHadrons[ih]->GetBaryonNumber();
538	}
539	G4int nQuas=theQuasmons.size();
540	if(nQuas) for(G4int iq=0; iq<nQuas; iq++)
541	{
542	finCharge+=theQuasmons[iq]->GetCharge();
543	finBaryoN+=theQuasmons[iq]->GetBaryonNumber();
544	}
545	G4cout<<"G4QuaStr::C:n="<<nHad<<",N="<<nQuas<<",C="<<finCharge<<",B="<<finBaryoN<<G4endl;
546	if(finCharge!=totCharge \|\| finBaryoN!=totBaryNum)
547	{
548	G4cerr<<"***G4QuasmonString::C: tC="<<totCharge<<",C="<<finCharge<<",tB="<<totBaryNum
549	<<",B="<<finBaryoN<<G4endl;
550	if(nHad) for(G4int h=0; h<nHad; h++)
551	{
552	G4QHadron* cH = theQHadrons[h];
553	G4cerr<<">G4QuaStr::C:h#"<<h<<",QC="<<cH->GetQC()<<",PDG="<<cH->GetPDGCode()<<G4endl;
554	}
555	if(nQuas) for(G4int q=0; q<nQuas; q++)
556	{
557	G4Quasmon* cQ = theQuasmons[q];
558	G4cerr<<">G4QuaStr::C:q#"<<q<<",C="<<cQ->GetCharge()<<",QCont="<<cQ->GetQC()<<G4endl;
559	}
560	}
561	#endif
562	} // End of the G4QuasmonString constructor
563
564	G4QuasmonString::G4QuasmonString(const G4QuasmonString &right)
565	{
566	// theQuasmons (Vector)
567	G4int nQ = right.theQuasmons.size();
568	if(nQ) for(G4int iq=0; iq<nQ; iq++)
569	{
570	G4Quasmon* curQ = new G4Quasmon(right.theQuasmons[iq]);
571	#ifdef fdebug
572	G4cout<<"G4QuaStr::CopyByVal:Q#"<<iq<<","<<curQ->GetQC()<<curQ->Get4Momentum()<<G4endl;
573	#endif
574	theQuasmons.push_back(curQ); // (delete equivalent)
575	}
576	theProjEnvFlag = right.theProjEnvFlag;
577	theTargEnvFlag = right.theTargEnvFlag;
578	theWeight = right.theWeight;
579	theProjQC = right.theProjQC;
580	theTargQC = right.theTargQC;
581	theProj4Mom = right.theProj4Mom;
582	theTarg4Mom = right.theTarg4Mom;
583
584	theWorld = right.theWorld;
585	tot4Mom = right.tot4Mom;
586	totCharge = right.totCharge;
587	totBaryNum = right.totBaryNum;
588	}
589
590	const G4QuasmonString& G4QuasmonString::operator=(const G4QuasmonString &right)
591	{// ========================================================================
592	if(this != &right) // Beware of self assignment
593	{
594	// theQuasmons (Vector)
595	G4int iQ = theQuasmons.size();
596	if(iQ) for(G4int jq=0; jq<iQ; jq++) delete theQuasmons[jq];
597	theQuasmons.clear();
598	G4int nQ = right.theQuasmons.size();
599	if(nQ) for(G4int iq=0; iq<nQ; iq++)
600	{
601	G4Quasmon* curQ = new G4Quasmon(right.theQuasmons[iq]);
602	#ifdef fdebug
603	G4cout<<"G4QuaStr::CpByVal:Q#"<<iq<<","<<curQ->GetQC()<<curQ->Get4Momentum()<<G4endl;
604	#endif
605	theQuasmons.push_back(curQ); // (delete equivalent)
606	}
607	theProjEnvFlag = right.theProjEnvFlag;
608	theTargEnvFlag = right.theTargEnvFlag;
609	theWeight = right.theWeight;
610	theProjQC = right.theProjQC;
611	theTargQC = right.theTargQC;
612	theProj4Mom = right.theProj4Mom;
613	theTarg4Mom = right.theTarg4Mom;
614
615	theWorld = right.theWorld;
616	tot4Mom = right.tot4Mom;
617	totCharge = right.totCharge;
618	totBaryNum = right.totBaryNum;
619	}
620	return *this;
621	}
622
623	G4QuasmonString::G4QuasmonString(G4QuasmonString* right)
624	{
625	// theQuasmons (Vector)
626	G4int nQ = right->theQuasmons.size();
627	if(nQ) for(G4int iq=0; iq<nQ; iq++)
628	{
629	G4Quasmon* curQ = new G4Quasmon(right->theQuasmons[iq]);
630	#ifdef fdebug
631	G4cout<<"G4QuaStr::CpByPoint:Q#"<<iq<<","<<curQ->GetQC()<<curQ->Get4Momentum()<<G4endl;
632	#endif
633	theQuasmons.push_back(curQ); // (delete equivalent)
634	}
635	theProjEnvFlag = right->theProjEnvFlag;
636	theTargEnvFlag = right->theTargEnvFlag;
637	theWeight = right->theWeight;
638	theProjQC = right->theProjQC;
639	theTargQC = right->theTargQC;
640	theProj4Mom = right->theProj4Mom;
641	theTarg4Mom = right->theTarg4Mom;
642
643	theWorld = right->theWorld;
644	tot4Mom = right->tot4Mom;
645	totCharge = right->totCharge;
646	totBaryNum = right->totBaryNum;
647	}
648
649	G4QuasmonString::~G4QuasmonString()
650	{
651	#ifdef debug
652	G4cout<<"~G4QuasmonString: before theQHadrons nH="<<theQHadrons.size()<<G4endl;
653	#endif
654	for_each(theQHadrons.begin(), theQHadrons.end(), DeleteQHadron());
655	#ifdef debug
656	G4cout<<"~G4QuasmonString: before theQuasmons nQ="<<theQuasmons.size()<<G4endl;
657	#endif
658	for_each(theQuasmons.begin(), theQuasmons.end(), DeleteQuasmon());
659	#ifdef debug
660	G4cout<<"~G4QuasmonString: === DONE ==="<<G4endl;
661	#endif
662	}
663
664	// ================== Initialization of Static Parameters ============
665	//G4double G4QuasmonString::StParName=0.; // Static Parameter (Example)
666	//G4bool G4QuasmonString::stFlag=false; // Static Flag (Example)
667
668	// Fill the private static parameters
669	//void G4QuasmonString::SetParameters(G4double aStPar, G4bool aStFlag)
670	//{// ====================================================================================
671	// StParName=aStPar; // (Example)
672	// stFlag=aStFlag; // (Example)
673	//}
674
675
676	//The public Hadronisation function with the Exception treatment (del respons. of User !)
677	G4QHadronVector* G4QuasmonString::Fragment()
678	{// ========================== -- @@ Must be changed @@ --
679	// Make the final check before filling the output -- @@ Must be changed @@ --
680	#ifdef chdebug
681	G4int fCharge=0;
682	G4int fBaryoN=0;
683	G4int nHad=theQHadrons.size();
684	if(nHad) for(G4int ih=0; ih<nHad; ih++)
685	{
686	fCharge+=theQHadrons[ih]->GetCharge();
687	fBaryoN+=theQHadrons[ih]->GetBaryonNumber();
688	}
689	G4int nQuas=theQuasmons.size();
690	if(nQuas)for(G4int iqs=0; iqs<nQuas; iqs++)
691	{
692	fCharge+=theQuasmons[iqs]->GetCharge();
693	fBaryoN+=theQuasmons[iqs]->GetBaryonNumber();
694	}
695	if(fCharge!=totCharge \|\| fBaryoN!=totBaryNum)
696	{
697	G4cerr<<"***G4QuasmonString::Fragment: tC="<<totCharge<<", C="<<fCharge<<", tB="
698	<<totBaryNum<<", B="<<fBaryoN<<G4endl;
699	if(nHad) for(G4int h=0; h<nHad; h++)
700	{
701	G4QHadron* cH = theQHadrons[h];
702	G4cerr<<":G4QuaSt::HQ:h#"<<h<<",QC="<<cH->GetQC()<<",PDG="<<cH->GetPDGCode()<<G4endl;
703	}
704	if(nQuas) for(G4int q=0; q<nQuas; q++)
705	{
706	G4Quasmon* cQ = theQuasmons[q];
707	G4cerr<<":G4QuaSt::HQ:q#"<<q<<",C="<<cQ->GetCharge()<<",QCont="<<cQ->GetQC()<<G4endl;
708	}
709	}
710	#endif
711	G4QHadronVector dummy; // Prototype of the output G4QuasmonVector to avoid wornings
712	G4QHadronVector* finalQHadrons = &dummy; // Prototype of the output G4QuasmonVector
713	G4int nH = theQHadrons.size();
714	if(nH)
715	{
716	for(G4int ih=0; ih<nH; ih++)
717	{
718	G4QHadron* curH = new G4QHadron(theQHadrons[ih]);
719	finalQHadrons->push_back(curH); // deleted after the "while LOOP"
720	}
721	}
722	for_each(theQuasmons.begin(),theQuasmons.end(),DeleteQuasmon()); //CleanUp Quasm's
723	theQuasmons.clear();
724	return finalQHadrons;
725	} // End of the Fragmentation member function
726
727	//The public extraction of the number of the created (in Constructor) G4QHadrons
728	G4int G4QuasmonString::GetNOfHadrons() {return theQHadrons.size();}
729
730	//The public extraction of the number of the created (in Constructor) G4Quasmons
731	G4int G4QuasmonString::GetNOfQuasmons() {return theQuasmons.size();}
732
733	//The public extraction of the projectile environment flag ("true" - exists)
734	G4bool G4QuasmonString::GetProjEnvFlag() {return theProjEnvFlag;}
735
736	//The public extraction of the target environment flag ("true" - exists)
737	G4bool G4QuasmonString::GetTargEnvFlag() {return theTargEnvFlag;}
738
739	//The public Quasmons duplication with delete responsibility of User (!)
740	G4QuasmonVector* G4QuasmonString::GetQuasmons()
741	{// ==============================
742	G4int nQ=theQuasmons.size();
743	#ifdef pdebug
744	G4cout<<"G4QuasmonString::GetQuasmons is called nQ="<<nQ<<G4endl;
745	#endif
746	G4QuasmonVector* quasmons = new G4QuasmonVector; // Intermediate
747	if(nQ) for(G4int iq=0; iq<nQ; iq++)
748	{
749	#ifdef pdebug
750	G4cout<<"G4QStr::GetQuasmons:Q#"<<iq<<",QPDG="<<theQuasmons[iq]->GetQPDG()<<",QQC="
751	<<theQuasmons[iq]->GetQC()<<",Q4M="<<theQuasmons[iq]->Get4Momentum()<<G4endl;
752	#endif
753	if(theQuasmons[iq]->GetStatus())
754	{
755	G4Quasmon* curQ = new G4Quasmon(theQuasmons[iq]);
756	quasmons->push_back(curQ); // (del. equiv. - user is responsibile)
757	}
758	}
759	#ifdef pdebug
760	G4cout<<"G4QuasmonString::GetQuasmons ===OUT==="<<G4endl;
761	#endif
762	return quasmons;
763	} // End of GetQuasmons
764
765	//The public Quasmons duplication with delete responsibility of User (!)
766	G4QHadronVector* G4QuasmonString::GetHadrons()
767	{// =============================
768	G4int nH=theQHadrons.size();
769	#ifdef pdebug
770	G4cout<<"G4QuasmonString::GetHadrons is called nH="<<nH<<G4endl;
771	#endif
772	G4QHadronVector* hadrons = new G4QHadronVector; // Intermediate
773	if(nH) for(G4int ih=0; ih<nH; ih++)
774	{
775	#ifdef pdebug
776	G4cout<<"G4QStr::GetHadrons: H#"<<ih<<",QPDG="<<theQHadrons[ih]->GetQPDG()<<",QQC="
777	<<theQHadrons[ih]->GetQC()<<",H_Mass="<<theQHadrons[ih]->GetMass()<<G4endl;
778	#endif
779	if(!theQHadrons[ih]->GetNFragments())
780	{
781	G4QHadron* curH = new G4QHadron(theQHadrons[ih]);
782	hadrons->push_back(curH); // (del. equiv. - user is responsibile)
783	}
784	}
785	#ifdef pdebug
786	G4cout<<"G4QuasmonString::GetHadrons ===OUT==="<<G4endl;
787	#endif
788	return hadrons;
789	} // End of GetHadrons
790
791	// Randomize the momentum fraction for the CHIPS of nPart free partons [x*(1-x)^(n-3)]
792	G4double G4QuasmonString::RandomizeMomFractionFree(G4int nPart)
793	{// ==============================================
794	// @@ TMP --- Begin ---
795	if(2>1)
796	{
797	if(nPart<2)
798	{
799	G4cerr<<"**G4QuasmonString::RandMomFractionString: n="<<nPart<<" < 2, retun 0"<<G4endl;
800	return 0.;
801	}
802	G4double x=0.5;
803	if(nPart==2) return x; // GS meson
804	G4double r=G4UniformRand();
805	if(r==0.) return 0.;
806	if(r==1.) return 1.;
807	if (nPart==3) x=r; // GS baryon
808	else if(nPart==4) x=1.-sqrt(r); // GS quaternion
809	else x=1.-pow(r,1./(nPart-2.)); // nPart>4
810	return x;
811	}
812	//if(nPart!=3) G4cerr<<">>>>>G4QuasmonString::RandMomFractionFree: n="<<nPart<<G4endl;
813	// @@ TMP --- End ---
814	if(nPart<2)
815	{
816	G4cerr<<"**G4QuasmonString::RandMomFractionFree: n="<<nPart<<" < 2, retun 0."<<G4endl;
817	return 0.;
818	}
819	G4double x=0.5;
820	if(nPart==2) return x; // GS meson
821	G4double r=G4UniformRand();
822	if(r==0.) return 0.;
823	if(r==1.) return 1.;
824	if(nPart==3) x=sqrt(r); // GS baryon
825	else if(nPart==4) // GS quaternion
826	{
827	if (r==0.5) x=0.5;
828	else if(r<0.5) x=sqrt(r+r)(.5+.1579(r-.5));
829	else x=1.-sqrt(2.-r-r)(.5+.1579(.5-r));
830	}
831	else
832	{
833	G4int n1=nPart-2;
834	G4double r1=n1;
835	G4double r2=r1-1.;
836	G4double rr=r2/r1;
837	G4double rp=pow(rr,n1);
838	G4double p2=rp+rp;
839	if (r==rr) x=p2;
840	else
841	{
842	if(r<rr)
843	{
844	G4double pr=0.;
845	G4double pra=0.;
846	if(nPart>8)
847	{
848	if(nPart>10)
849	{
850	if(nPart>11) // >11
851	{
852	pr=.614/pow((nPart+1.25),.75);
853	pra=.915/pow((nPart+6.7),1.75);
854	}
855	else // 11
856	{
857	pr=.09945;
858	pra=.00667;
859	}
860	}
861	else
862	{
863	if(nPart>9) // 10
864	{
865	pr=.1064;
866	pra=.00741;
867	}
868	else // 9
869	{
870	pr=.11425;
871	pra=.00828;
872	}
873	}
874	}
875	else
876	{
877	if(nPart>6)
878	{
879	if(nPart>7) // 8
880	{
881	pr=.12347;
882	pra=.00926;
883	}
884	else // 7
885	{
886	pr=.13405;
887	pra=.01027;
888	}
889	}
890	else
891	{
892	if(nPart>5) // 6
893	{
894	pr=.1454;
895	pra=.01112;
896	}
897	else // 5
898	{
899	pr=.15765;
900	pra=.00965;
901	}
902	}
903	}
904	x=pow((r/p2),(1.-rr+pra))(rr+pr(r-p2));
905	}
906	else
907	{
908	G4double sr=0.;
909	if(nPart>8)
910	{
911	if(nPart>10)
912	{
913	if(nPart>11) sr=.86/(nPart+1.05); // >11
914	else sr=.0774; // 11
915	}
916	else
917	{
918	if(nPart>9) sr=.0849; // 10
919	else sr=.0938; // 9
920	}
921	}
922	else
923	{
924	if(nPart>6)
925	{
926	if(nPart>7) sr=.1047; // 8
927	else sr=.1179; // 7
928	}
929	else
930	{
931	if(nPart>5) sr=.1339; // 6
932	else sr=.15135; // 5
933	}
934	}
935	x=1.-sqrt((1.-r)/(1.-p2))(1.-rr+sr(p2-r));
936	}
937	}
938	}
939	return x;
940	} // End of RandomizeMomFractionFree
941
942	// Randomize MomentumFraction for CHIPS of nPart partons for Pomeron exchange [(1-x)^(n-2)]
943	G4double G4QuasmonString::RandomizeMomFractionString(G4int nPart)
944	{// ================================================
945	if(nPart<2)
946	{
947	G4cerr<<"**G4QuasmonString::RandMomFractionString: n="<<nPart<<" < 2, retun 0"<<G4endl;
948	return 0.;
949	}
950	G4double x=0.5;
951	if(nPart==2) return x; // GS meson
952	G4double r=G4UniformRand();
953	if(r==0.) return 0.;
954	if(r==1.) return 1.;
955	if (nPart==3) x=r; // GS baryon
956	else if(nPart==4) x=1.-sqrt(r); // GS quaternion
957	else x=1.-pow(r,1./(nPart-2.)); // nPart>4
958	return x;
959	} // End of RandomizeMomFractionString

Note: See TracBrowser for help on using the repository browser.

Download in other formats: