Context Navigation

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

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

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

Download in other formats: