Context Navigation

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

Last change on this file since 836 was 819, checked in by garnier, 16 years ago
import all except CVS
File size: 40.5 KB

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

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

Download in other formats: