Context Navigation

source: trunk/source/processes/hadronic/models/chiral_inv_phase_space/body/src/G4QSplitter.cc @ 1055

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

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

Download in other formats: