Context Navigation

source: trunk/source/processes/hadronic/models/chiral_inv_phase_space/interface/src/G4QDiffractionRatio.cc @ 968

Last change on this file since 968 was 962, checked in by garnier, 15 years ago
update processes
File size: 63.5 KB

Line
1	//
2	// ********************************************************************
3	// * License and Disclaimer *
4	// * *
5	// * The Geant4 software is copyright of the Copyright Holders of *
6	// * the Geant4 Collaboration. It is provided under the terms and *
7	// * conditions of the Geant4 Software License, included in the file *
8	// * LICENSE and available at http://cern.ch/geant4/license . These *
9	// * include a list of copyright holders. *
10	// * *
11	// * Neither the authors of this software system, nor their employing *
12	// * institutes,nor the agencies providing financial support for this *
13	// * work make any representation or warranty, express or implied, *
14	// * regarding this software system or assume any liability for its *
15	// * use. Please see the license in the file LICENSE and URL above *
16	// * for the full disclaimer and the limitation of liability. *
17	// * *
18	// * This code implementation is the result of the scientific and *
19	// * technical work of the GEANT4 collaboration. *
20	// * By using, copying, modifying or distributing the software (or *
21	// * any work based on the software) you agree to acknowledge its *
22	// * use in resulting scientific publications, and indicate your *
23	// * acceptance of all terms of the Geant4 Software license. *
24	// ********************************************************************
25	//
26	//
27	// $Id: G4QDiffractionRatio.cc,v 1.9 2008/03/21 21:40:08 dennis Exp $
28	// GEANT4 tag $Name: geant4-09-02-ref-02 $
29	//
30	//
31	// G4 Physics class: G4QDiffractionRatio for N+A elastic cross sections
32	// Created: M.V. Kossov, CERN/ITEP(Moscow), 10-OCT-01
33	// The last update: M.V. Kossov, CERN/ITEP (Moscow) 15-Oct-06
34	//
35	//================================================================================
36
37	//#define debug
38	//#define pdebug
39	//#define fdebug
40	//#define nandebug
41
42	#include "G4QDiffractionRatio.hh"
43
44	// Returns Pointer to the G4VQCrossSection class
45	G4QDiffractionRatio* G4QDiffractionRatio::GetPointer()
46	{
47	static G4QDiffractionRatio theRatios; // * Static body of the Diffraction Ratio *
48	return &theRatios;
49	}
50
51	// Calculation of pair(QuasiFree/Inelastic,QuasiElastic/QuasiFree)
52	G4double G4QDiffractionRatio::GetRatio(G4double pIU, G4int pPDG, G4int tgZ, G4int tgN)
53	{
54	static const G4double mNeut= G4QPDGCode(2112).GetMass();
55	static const G4double mProt= G4QPDGCode(2212).GetMass();
56	static const G4double mN=.5*(mNeut+mProt);
57	static const G4double dmN=mN+mN;
58	static const G4double mN2=mN*mN;
59	// Table parameters
60	static const G4int nps=100; // Number of steps in the R(s) LinTable
61	static const G4int mps=nps+1; // Number of elements in the R(s) LinTable
62	static const G4double sma=6.; // The first LinTabEl(sqs=0)=1., sqs>sma -> logTab
63	static const G4double ds=sma/nps; // Step of the linear Table
64	static const G4int nls=150; // Number of steps in the R(lns) logTable
65	static const G4int mls=nls+1; // Number of elements in the R(lns) logTable
66	static const G4double lsi=1.79; // The min ln(sqs) logTabEl(sqs=5.99 < sma=6.)
67	static const G4double lsa=8.; // The max ln(sqs) logTabEl(sqs=5.99 - 2981 GeV)
68	static const G4double mi=std::exp(lsi);// The min s of logTabEl(~ 5.99 GeV)
69	static const G4double ms=std::exp(lsa);// The max s of logTabEl(~ 2981 GeV)
70	static const G4double dl=(lsa-lsi)/nls;// Step of the logarithmic Table
71	static const G4double edl=std::exp(dl);// Multiplication step of the logarithmic Table
72	static const G4double toler=.0001; // Tolarence (GeV) defining the same sqs
73	static G4double lastS=0.; // Last sqs value for which R was calculated
74	static G4double lastR=0.; // Last ratio R which was calculated
75	// Local Associative Data Base:
76	static std::vector<G4int> vA; // Vector of calculated A
77	static std::vector<G4double> vH; // Vector of max sqs initialized in the LinTable
78	static std::vector<G4int> vN; // Vector of topBin number initialized in LinTable
79	static std::vector<G4double> vM; // Vector of relMax ln(sqs) initialized in LogTable
80	static std::vector<G4int> vK; // Vector of topBin number initialized in LogTable
81	static std::vector<G4double*> vT; // Vector of pointers to LinTable in C++ heap
82	static std::vector<G4double*> vL; // Vector of pointers to LogTable in C++ heap
83	// Last values of the Associative Data Base:
84	static G4int lastPDG=0; // Last PDG for which R was calculated (now indep)
85	static G4int lastA=0; // theLast of calculated A
86	static G4double lastH=0.; // theLast of max sqs initialized in the LinTable
87	static G4int lastN=0; // theLast of topBin number initialized in LinTable
88	static G4double lastM=0.; // theLast of relMax ln(sqs) initialized in LogTab.
89	static G4int lastK=0; // theLast of topBin number initialized in LogTable
90	static G4double* lastT=0; // theLast of pointer to LinTable in the C++ heap
91	static G4double* lastL=0; // theLast of pointer to LogTable in the C++ heap
92	// LogTable is created only if necessary. R(sqs>2981GeV) calcul by formula for any nuclei
93	G4int A=tgN+tgZ;
94	if(pIU<toler \|\| A<1) return 1.; // Fake use of toler as non zero number
95	if(A>238)
96	{
97	G4cout<<"--Warning--G4QuasiFreeRatio::GetRatio: A="<<A<<">238, return zero"<<G4endl;
98	return 0.;
99	}
100	lastPDG=pPDG; // @@ at present ratio is PDG independent @@
101	// Calculate sqs
102	G4double pM=G4QPDGCode(pPDG).GetMass()*.001; // Projectile mass in GeV
103	G4double pM2=pM*pM;
104	G4double mom=pIU/gigaelectronvolt; // Projectile momentum in GeV
105	G4double s=std::sqrt(mN2+pM2+dmNstd::sqrt(pM2+mommom));
106	G4int nDB=vA.size(); // A number of nuclei already initialized in AMDB
107	// if(nDB && lastA==A && std::fabs(s-lastS)<toler) return lastR;
108	if(nDB && lastA==A && s==lastS) return lastR; // VI do not use toler
109	if(s>ms)
110	{
111	lastR=CalcDiff2Prod_Ratio(s,A); // @@ Probably user ought to be notified about bigS
112	return lastR;
113	}
114	G4bool found=false;
115	G4int i=-1;
116	if(nDB) for (i=0; i<nDB; i++) if(A==vA[i]) // Sirch for this A in AMDB
117	{
118	found=true; // The A value is found
119	break;
120	}
121	if(!nDB \|\| !found) // Create new line in the AMDB
122	{
123	lastA = A;
124	lastT = new G4double[mps]; // Create the linear Table
125	lastN = static_cast<int>(s/ds)+1; // MaxBin to be initialized
126	if(lastN>nps)
127	{
128	lastN=nps;
129	lastH=sma;
130	}
131	else lastH = lastN*ds; // Calculate max initialized s for LinTab
132	G4double sv=0;
133	lastT[0]=1.;
134	for(G4int j=1; j<=lastN; j++) // Calculate LinTab values
135	{
136	sv+=ds;
137	lastT[j]=CalcDiff2Prod_Ratio(sv,A);
138	}
139	if(s>sma) // Initialize the logarithmic Table
140	{
141	lastL=new G4double[mls]; // Create the logarithmic Table
142	G4double ls=std::log(s);
143	lastK = static_cast<int>((ls-lsi)/dl)+1; // MaxBin to be initialized in LogTaB
144	if(lastK>nls)
145	{
146	lastK=nls;
147	lastM=lsa-lsi;
148	}
149	else lastM = lastK*dl; // Calculate max initialized ln(s)-lsi for LogTab
150	sv=mi;
151	for(G4int j=0; j<=lastK; j++) // Calculate LogTab values
152	{
153	lastL[j]=CalcDiff2Prod_Ratio(sv,A);
154	if(j!=lastK) sv*=edl;
155	}
156	}
157	else // LogTab is not initialized
158	{
159	lastL = 0;
160	lastK = 0;
161	lastM = 0.;
162	}
163	i++; // Make a new record to AMDB and position on it
164	vA.push_back(lastA);
165	vH.push_back(lastH);
166	vN.push_back(lastN);
167	vM.push_back(lastM);
168	vK.push_back(lastK);
169	vT.push_back(lastT);
170	vL.push_back(lastL);
171	}
172	else // The A value was found in AMDB
173	{
174	lastA=vA[i];
175	lastH=vH[i];
176	lastN=vN[i];
177	lastM=vM[i];
178	lastK=vK[i];
179	lastT=vT[i];
180	lastL=vL[i];
181	if(s>lastM) // At least LinTab must be updated
182	{
183	G4int nextN=lastN+1; // The next bin to be initialized
184	if(lastN<nps)
185	{
186	lastN = static_cast<int>(s/ds)+1;// MaxBin to be initialized
187	if(lastN>nps)
188	{
189	lastN=nps;
190	lastH=sma;
191	}
192	else lastH = lastN*ds; // Calculate max initialized s for LinTab
193	G4double sv=lastM;
194	for(G4int j=nextN; j<=lastN; j++)// Calculate LogTab values
195	{
196	sv+=ds;
197	lastT[j]=CalcDiff2Prod_Ratio(sv,A);
198	}
199	} // End of LinTab update
200	if(lastN>=nextN)
201	{
202	vH[i]=lastH;
203	vN[i]=lastN;
204	}
205	G4int nextK=lastK+1;
206	if(s>sma && lastK<nls) // LogTab must be updated
207	{
208	G4double sv=std::exp(lastM+lsi); // Define starting poit (lastM will be changed)
209	G4double ls=std::log(s);
210	lastK = static_cast<int>((ls-lsi)/dl)+1; // MaxBin to be initialized in LogTaB
211	if(lastK>nls)
212	{
213	lastK=nls;
214	lastM=lsa-lsi;
215	}
216	else lastM = lastK*dl; // Calculate max initialized ln(s)-lsi for LogTab
217	for(G4int j=nextK; j<=lastK; j++)// Calculate LogTab values
218	{
219	sv*=edl;
220	lastL[j]=CalcDiff2Prod_Ratio(sv,A);
221	}
222	} // End of LogTab update
223	if(lastK>=nextK)
224	{
225	vM[i]=lastM;
226	vK[i]=lastK;
227	}
228	}
229	}
230	// Now one can use tabeles to calculate the value
231	if(s<sma) // Use linear table
232	{
233	G4int n=static_cast<int>(s/ds); // Low edge number of the bin
234	G4double d=s-n*ds; // Linear shift
235	G4double v=lastT[n]; // Base
236	lastR=v+d*(lastT[n+1]-v)/ds; // Result
237	}
238	else // Use log table
239	{
240	G4double ls=std::log(s)-lsi; // ln(s)-l_min
241	G4int n=static_cast<int>(ls/dl); // Low edge number of the bin
242	G4double d=ls-n*dl; // Log shift
243	G4double v=lastL[n]; // Base
244	lastR=v+d*(lastL[n+1]-v)/dl; // Result
245	}
246	if(lastR<0.) lastR=0.;
247	if(lastR>1.) lastR=1.;
248	return lastR;
249	} // End of CalcDiff2Prod_Ratio
250
251	// Calculate Diffraction/Production Ratio as a function of total sq(s)(hN) (in GeV), A=Z+N
252	G4double G4QDiffractionRatio::CalcDiff2Prod_Ratio(G4double s, G4int A)
253	{
254	static G4int mA=0;
255	static G4double S=.1; // s=SQRT(M_N^2+M_h^2+2E_hM_N)
256	static G4double R=0.; // Prototype of the result
257	static G4double p1=0.;
258	static G4double p2=0.;
259	static G4double p4=0.;
260	static G4double p5=0.;
261	static G4double p6=0.;
262	static G4double p7=0.;
263	if(s<=0. \|\| A<=1) return 0.;
264	if(A!=mA && A!=1)
265	{
266	mA=A;
267	G4double a=mA;
268	G4double sa=std::sqrt(a);
269	G4double a2=a*a;
270	G4double a3=a2*a;
271	G4double a4=a3*a;
272	G4double a5=a4*a;
273	G4double a6=a5*a;
274	G4double a7=a6*a;
275	G4double a8=a7*a;
276	G4double a11=a8*a3;
277	G4double a12=a8*a4;
278	G4double p=std::pow(a,0.37);
279	p1=(.023p+3.5/a3+2.1e6/a12+4.e-14a5)/(1.+7.6e-4asa+2.15e7/a11);
280	p2=(1.42std::pow(a,0.61)+1.6e5/a8+4.5e-8a4)/(1.+4.e-8*a4+1.2e4/a6);
281	G4double q=std::pow(a,0.7);
282	p4=(.036/q+.0009q)/(1.+6./a3+1.e-7a3);
283	p5=1.3std::pow(a,0.1168)/(1.+1.2e-8a3);
284	p6=.00046*(a+11830./a2);
285	p7=1./(1.+6.17/a2+.00406*a);
286	}
287	else if(A==1 && mA!=1)
288	{
289	p1=.0315;
290	p2=.73417;
291	p4=.01109;
292	p5=1.0972;
293	p6=.065787;
294	p7=.62976;
295	}
296	else if(std::fabs(s-S)/S<.0001) return R;
297	G4double s2=s*s;
298	G4double s4=s2*s2;
299	G4double dl=std::log(s)-p5;
300	R=1./(1.+1./(p1+p2/s4+p4dldl/(1.+p6*std::pow(s,p7))));
301	return R;
302	} // End of CalcQF2IN_Ratio
303
304
305	G4QHadronVector* G4QDiffractionRatio::TargFragment(G4int pPDG, G4LorentzVector p4M,
306	G4int tgZ, G4int tgN)
307	{
308	static const G4double pFm= 0.; // Fermi momentum in MeV (delta function)
309	//static const G4double pFm= 250.; // Fermi momentum in MeV (delta function)
310	static const G4double pFm2= pFm*pFm; // Squared Fermi momentum in MeV^2 (delta function)
311	static const G4double mPi0= G4QPDGCode(111).GetMass(); // pi0 mass (MeV =min diffraction)
312	//static const G4double mPi= G4QPDGCode(211).GetMass(); // pi+- mass (MeV)
313	static const G4double mNeut= G4QPDGCode(2112).GetMass();
314	static const G4double mNeut2=mNeut*mNeut;
315	static const G4double dmNeut=mNeut+mNeut;
316	static const G4double mProt= G4QPDGCode(2212).GetMass();
317	static const G4double mProt2=mProt*mProt;
318	static const G4double dmProt=mProt+mProt;
319	static const G4double maxDM=mProt*12.;
320	//static const G4double mLamb= G4QPDGCode(3122).GetMass();
321	//static const G4double mSigZ= G4QPDGCode(3212).GetMass();
322	//static const G4double mSigM= G4QPDGCode(3112).GetMass();
323	//static const G4double mSigP= G4QPDGCode(3222).GetMass();
324	//static const G4double eps=.003;
325	static const G4double third=1./3.;
326	//
327	G4LorentzVector pr4M=p4M/megaelectronvolt; // Convert 4-momenta in MeV (keep p4M)
328	// prepare the DONOTHING answer
329	G4QHadronVector* ResHV = new G4QHadronVector;// !! MUST BE DESTROYE/DDELETER by CALLER !!
330	G4QHadron* hadron = new G4QHadron(pPDG,p4M); // Hadron for not scattered projectile
331	ResHV->push_back(hadron); // It must be cleaned up for real scattering sec
332	// @@ diffraction is simulated as noncoherent (coherent is small)
333	G4int tgA=tgZ+tgN; // A of the target
334	G4int tPDG=90000000+tgZ*1000+tgN; // PDG code of the targetNucleus/recoilNucleus
335	G4double tgM=G4QPDGCode(tPDG).GetMass(); // Mass of the target nucleus
336	G4int rPDG=2112; // prototype of PDG code of the recoiled nucleon
337	if(tgA*G4UniformRand()>tgN) // Substitute by a proton
338	{
339	rPDG=2212; // PDG code of the recoiled QF nucleon
340	tPDG-=1000; // PDG code of the recoiled nucleus
341	}
342	else tPDG-=1; // PDG code of the recoiled nucleus
343	G4double tM=G4QPDGCode(tPDG).GetMass(); // Mass of the recoiled nucleus
344	G4double tE=std::sqrt(tM*tM+pFm2); // Free energy of the recoil nucleus
345	G4ThreeVector tP=pFm*G4RandomDirection();// 3-mom of the recoiled nucleus
346	G4LorentzVector t4M(tP,tE); // 4M of the recoil nucleus
347	G4LorentzVector tg4M(0.,0.,0.,tgM); // Full target 4-momentum
348	G4LorentzVector N4M=tg4M-t4M; // 4-mom of Quasi-free target nucleon
349	G4LorentzVector tot4M=N4M+p4M; // total momentum of quasi-free diffraction
350	G4double mT=mNeut; // Prototype of mass of QF nucleon
351	G4double mT2=mNeut2; // Squared mass of a free nucleon to be excited
352	G4double dmT=dmNeut; // Doubled mass
353	G4int Z=0; // Prototype of the isotope Z
354	G4int N=1; // Prototype of the Isotope N
355	if(rPDG==2212) // Correct it, if this is a proton
356	{
357	mT=mProt; // Prototype of mass of QF nucleon to be excited
358	mT2=mProt2; // Squared mass of the free nucleon
359	dmT=dmProt; // Doubled mass
360	Z=1; // Z of the isotope
361	N=0; // N of the Isotope
362	}
363	G4double mP2=pr4M.m2(); // Squared mass of the projectile
364	if(mP2<0.) mP2=0.; // Can be a problem for photon (m_min = 2*m_pi0)
365	G4double s=tot4M.m2(); // @@ Check <0 ...
366	G4double E=(s-mT2-mP2)/dmT; // Effective interactinEnergy (virt.nucl.target)
367	G4double E2=E*E;
368	if(E<0. \|\| E2<mP2)
369	{
370	#ifdef pdebug
371	G4cerr<<"_Warning-G4DifR::TFra:<NegativeEnergy>E="<<E<<",E2="<<E2<<"<M2="<<mP2<<G4endl;
372	#endif
373	return ResHV; // Do Nothing Action
374	}
375	G4double mP=std::sqrt(mP2); // Calculate mass of the projectile (to be exc.)
376	if(mP<.1)mP=mPi0; // For photons minDiffraction is gam+P->P+Pi0
377	G4double dmP=mP+mP; // Doubled mass of the projectile
378	G4double mMin=mP+mPi0; // Minimum diffractive mass
379	G4double tA=tgA; // Real A of the target
380	G4double sA=5./std::pow(tA,third); // Mass-screaning
381	//mMin+=mPi0+G4UniformRand()(mPsA+mPi0); // Experimental
382	mMin+=G4UniformRand()(mPsA+mPi0); // Experimental
383	G4double ss=std::sqrt(s); // CM compound mass (sqrt(s))
384	G4double mMax=ss-mP; // Maximum diffraction mass of the projectile
385	if(mMax>maxDM) mMax=maxDM; // Restriction to avoid too big masses
386	if(mMin>=mMax)
387	{
388	#ifdef pdebug
389	G4cerr<<"Warning-G4DifR::TFra:ZeroDiffractionMRange, mi="<<mMin<<", ma="<<mMax<<G4endl;
390	#endif
391	return ResHV; // Do Nothing Action
392	}
393	G4double R = G4UniformRand();
394	G4double mDif=std::exp(Rstd::log(mMax)+(1.-R)std::log(mMin)); // Low Mass Approximation
395	G4double mDif2=mDif*mDif;
396	G4double ds=s-mP2-mDif2;
397	G4double e=ds/dmP;
398	G4double P=std::sqrt(e*e-mDif2); // Momentum in pseudo laboratory system
399	G4VQCrossSection* CSmanager=G4QElasticCrossSection::GetPointer();
400	#ifdef debug
401	G4cout<<"G4QDiffR::TargFrag:Before XS, P="<<P<<",Z="<<Z<<",N="<<N<<",PDG="<<pPDG<<G4endl;
402	#endif
403	// @@ Temporary NN t-dependence for all hadrons
404	if(pPDG>3400 \|\| pPDG<-3400) G4cout<<"-Warning-G4QDifR::Fragment: pPDG="<<pPDG<<G4endl;
405	G4int PDG=2212; // TMP instead of pPDG
406	G4double xSec=CSmanager->GetCrossSection(false, P, tgZ, tgN, PDG);// Rec.CrossSect TMP
407	//G4double xSec=CSmanager->GetCrossSection(false, P, tgZ, tgN, pPDG); // Rec.CrossSect
408	#ifdef debug
409	G4cout<<"G4QDiffRat::TargFragment:pPDG="<<pPDG<<",P="<<P<<",CS="<<xSec/millibarn<<G4endl;
410	#endif
411	#ifdef nandebug
412	if(xSec>0. \|\| xSec<0. \|\| xSec==0);
413	else G4cout<<"*NAN*G4QDiffractionRatio::TargFragment:xSec="<<xSec/millibarn<<G4endl;
414	#endif
415	// @@ check a possibility to separate p, n, or alpha (!)
416	if(xSec <= 0.) // The cross-section iz 0 -> Do Nothing
417	{
418	#ifdef pdebug
419	G4cerr<<"-Warning-G4QDiffrRatio::TargFragment:Zero XSPDG="<<pPDG<<",P="<<P<<G4endl;
420	#endif
421	return ResHV; //Do Nothing Action
422	}
423	//G4double t=CSmanager->GetExchangeT(tgZ,tgN,pPDG); // functional randomized -t (MeV^2)
424	G4double maxt=(dsds-4mP2*mDif2)/s; // maximum possible -t
425	G4double tsl=140000.; // slope in MeV^2
426	G4double t=-std::log(G4UniformRand())*tsl;
427	#ifdef pdebug
428	G4cout<<"G4QDifR::TFra:ph="<<pPDG<<",P="<<P<<",X="<<xSec<<",t="<<t<<"<"<<maxt<<G4endl;
429	#endif
430	#ifdef nandebug
431	if(mint>-.0000001); // To make the Warning for NAN
432	else G4cout<<"******G4QDiffractionRatio::TargFragment: -t="<<mint<<G4endl;
433	#endif
434	G4double rt=t/maxt;
435	G4double cost=1.-rt-rt; // cos(theta) in CMS
436	#ifdef ppdebug
437	G4cout<<"G4QDiffraRatio::TargFragment: -t="<<t<<", maxt="<<maxt<<", cost="<<cost<<G4endl;
438	#endif
439	if(cost>1. \|\| cost<-1. \|\| !(cost>-1. \|\| cost<=1.))
440	{
441	if (cost>1.) cost=1.;
442	else if(cost<-1.) cost=-1.;
443	else
444	{
445	G4cerr<<"G4QDiffRat::TargFragm: NAN cost="<<cost<<",t="<<t<<",tmax="<<maxt<<G4endl;
446	return ResHV; // Do Nothing Action
447	}
448	}
449	G4LorentzVector r4M=G4LorentzVector(0.,0.,0.,mP); // 4mom of the leading nucleon
450	G4LorentzVector d4M=G4LorentzVector(0.,0.,0.,mDif); // 4mom of the diffract. Quasmon
451	G4LorentzVector dir4M=tot4M-G4LorentzVector(0.,0.,0.,(tot4M.e()-mT)*.01);
452	if(!G4QHadron(tot4M).RelDecayIn2(r4M, d4M, dir4M, cost, cost))
453	{
454	G4cerr<<"G4QDifR::TFr:M="<<tot4M.m()<<",T="<<mT<<",D="<<mDif<<",T+D="<<mT+mDif<<G4endl;
455	//G4Exception("G4QDifR::Fragm:","009",FatalException,"Decay of ElasticComp");
456	return ResHV; // Do Nothing Action
457	}
458	#ifdef debug
459	G4cout<<"G4QDifRat::TargFragm:d4M="<<d4M<<"+r4M="<<r4M<<"="<<d4M+r4M<<"="<<tot4M<<G4endl;
460	#endif
461	// Now everything is ready for fragmentation and DoNothing projHadron must be wiped out
462	ResHV->pop_back(); // Clean up pointer to the fake (doNothing) projectile
463	delete hadron; // Delete the fake (doNothing) projectile hadron
464	hadron = new G4QHadron(pPDG,r4M); // Hadron for the recoil nucleon
465	ResHV->push_back(hadron); // Fill the recoil nucleon
466	#ifdef debug
467	G4cout<<"G4QDiffractionRatio::TargFragm: Filled LeadingNuc="<<r4M<<pPDG<<G4endl;
468	#endif
469	G4QHadronVector* leadhs=new G4QHadronVector;// Quasmon Output G4QHadronVectorum --->---+
470	G4QContent dQC=G4QPDGCode(rPDG).GetQuarkContent(); // QuarkContent of quasiFreeNucleon \|
471	G4Quasmon* quasm = new G4Quasmon(dQC,d4M); // Quasmon=DiffractionExcitationQuasmon-+ \|
472	#ifdef debug
473	G4cout<<"G4QDiffRatio::TgFrag:tPDG="<<tPDG<<",rPDG="<<rPDG<<",d4M="<<d4M<<G4endl;//\| \|
474	#endif
475	G4QEnvironment* pan= new G4QEnvironment(G4QNucleus(tPDG));// --> DELETED --->---+ \| \|
476	pan->AddQuasmon(quasm); // Add diffractiveQuasmon to Environ.\| \| \|
477	#ifdef debug
478	G4cout<<"G4QDiffractionRatio::TargFragment: EnvPDG="<<tPDG<<G4endl; // \| \| \|
479	#endif
480	try // \| \| \|
481	{ // \| \| \|
482	delete leadhs; // \| \| \|
483	leadhs = pan->Fragment();// DESTROYED in the end of the LOOP work space \| \| \|
484	} // \| \| \|
485	catch (G4QException& error)// \| \| \|
486	{ // \| \| \|
487	//#ifdef pdebug
488	G4cerr<<"***G4QCollision::PostStepDoIt: G4QE Exception is catched"<<G4endl;// \| \| \|
489	//#endif
490	G4Exception("G4QCollision::PostStepDoIt:","27",FatalException,"CHIPSCrash");//\| \| \|
491	} // \| \| \|
492	delete pan; // Delete the Nuclear Environment <-<--+--+ \|
493	G4int qNH=leadhs->size(); // A#of collected hadrons from diff.frag. \|
494	if(qNH) for(G4int iq=0; iq<qNH; iq++) // Loop over hadrons to fill the result \|
495	{ // \|
496	G4QHadron* loh=(*leadhs)[iq]; // Pointer to the output hadron \|
497	ResHV->push_back(loh); // Fill in the result \|
498	} // \|
499	delete leadhs; // <----<----<----<----<----<----<----<----<----<----<----<----<----<---*
500	return ResHV; // Result
501	} // End of TargFragment
502
503
504	G4QHadronVector* G4QDiffractionRatio::ProjFragment(G4int pPDG, G4LorentzVector p4M,
505	G4int tgZ, G4int tgN)
506	{
507	static const G4double pFm= 250.; // Fermi momentum in MeV (delta function)
508	static const G4double pFm2= pFm*pFm; // Squared Fermi momentum in MeV^2 (delta function)
509	static const G4double mPi0= G4QPDGCode(111).GetMass(); // pi0 mass (MeV =min diffraction)
510	static const G4double mPi= G4QPDGCode(211).GetMass(); // pi+- mass (MeV)
511	static const G4double mNeut= G4QPDGCode(2112).GetMass();
512	static const G4double mNeut2=mNeut*mNeut;
513	static const G4double dmNeut=mNeut+mNeut;
514	static const G4double mProt= G4QPDGCode(2212).GetMass();
515	static const G4double mProt2=mProt*mProt;
516	static const G4double dmProt=mProt+mProt;
517	static const G4double maxDM=mProt*12.;
518	static const G4double mLamb= G4QPDGCode(3122).GetMass();
519	static const G4double mSigZ= G4QPDGCode(3212).GetMass();
520	static const G4double mSigM= G4QPDGCode(3112).GetMass();
521	static const G4double mSigP= G4QPDGCode(3222).GetMass();
522	static const G4double eps=.003;
523	//
524	G4LorentzVector pr4M=p4M/megaelectronvolt; // Convert 4-momenta in MeV (keep p4M)
525	// prepare the DONOTHING answer
526	G4QHadronVector* ResHV = new G4QHadronVector;// !! MUST BE DESTROYE/DDELETER by CALLER !!
527	G4QHadron* hadron = new G4QHadron(pPDG,p4M); // Hadron for not scattered projectile
528	ResHV->push_back(hadron); // It must be cleaned up for real scattering sec
529	// @@ diffraction is simulated as noncoherent (coherent is small)
530	G4int tgA=tgZ+tgN; // A of the target
531	G4int tPDG=90000000+tgZ*1000+tgN; // PDG code of the targetNucleus/recoilNucleus
532	G4double tgM=G4QPDGCode(tPDG).GetMass(); // Mass of the target nucleus
533	G4int rPDG=2112; // prototype of PDG code of the recoiled nucleon
534	if(tgA*G4UniformRand()>tgN) // Substitute by a proton
535	{
536	rPDG=2212; // PDG code of the recoiled QF nucleon
537	tPDG-=1000; // PDG code of the recoiled nucleus
538	}
539	else tPDG-=1; // PDG code of the recoiled nucleus
540	G4double tM=G4QPDGCode(tPDG).GetMass(); // Mass of the recoiled nucleus
541	G4double tE=std::sqrt(tM*tM+pFm2);
542	G4ThreeVector tP=pFm*G4RandomDirection();
543	G4LorentzVector t4M(tP,tE); // 4M of the recoil nucleus
544	G4LorentzVector tg4M(0.,0.,0.,tgM);
545	G4LorentzVector N4M=tg4M-t4M; // Quasi-free target nucleon
546	G4LorentzVector tot4M=N4M+p4M; // total momentum of quasi-free diffraction
547	G4double mT=mNeut;
548	G4double mT2=mNeut2; // Squared mass of the free nucleon spectator
549	G4double dmT=dmNeut;
550	G4int Z=0;
551	G4int N=1;
552	if(rPDG==2212)
553	{
554	mT=mProt;
555	mT2=mProt2;
556	dmT=dmProt;
557	Z=1;
558	N=0;
559	}
560	G4double mP2=pr4M.m2(); // Squared mass of the projectile
561	if(mP2<0.) mP2=0.; // A possible problem for photon (m_min = 2*m_pi0)
562	G4double s=tot4M.m2(); // @@ Check <0 ...
563	G4double E=(s-mT2-mP2)/dmT; // Effective interactin energy (virt. nucl. target)
564	G4double E2=E*E;
565	if(E<0. \|\| E2<mP2)
566	{
567	#ifdef pdebug
568	G4cerr<<"_Warning-G4DifR::PFra:<NegativeEnergy>E="<<E<<",E2="<<E2<<"<M2="<<mP2<<G4endl;
569	#endif
570	return ResHV; // Do Nothing Action
571	}
572	G4double mP=std::sqrt(mP2);
573	if(mP<.1)mP=mPi0; // For photons min diffraction is gamma+P->Pi0+Pi0
574	G4double mMin=mP+mPi0; // Minimum diffractive mass
575	G4double ss=std::sqrt(s); // CM compound mass (sqrt(s))
576	G4double mMax=ss-mT; // Maximum diffraction mass
577	if(mMax>maxDM) mMax=maxDM; // Restriction to avoid too big masses
578	if(mMin>=mMax)
579	{
580	#ifdef pdebug
581	G4cerr<<"Warning-G4DifR::PFra:ZeroDiffractionMRange, mi="<<mMin<<", ma="<<mMax<<G4endl;
582	#endif
583	return ResHV; // Do Nothing Action
584	}
585	G4double R = G4UniformRand();
586	G4double mDif=std::exp(Rstd::log(mMax)+(1.-R)std::log(mMin)); // LowMassApproximation
587	G4double mDif2=mDif*mDif;
588	G4double ds=s-mT2-mDif2;
589	G4double e=ds/dmT;
590	G4double P=std::sqrt(e*e-mDif2); // Momentum in pseudo laboratory system
591	G4VQCrossSection* CSmanager=G4QElasticCrossSection::GetPointer();
592	#ifdef debug
593	G4cout<<"G4QDiffR::PFra: Before XS, P="<<P<<", Z="<<Z<<", N="<<N<<", PDG="<<pPDG<<G4endl;
594	#endif
595	// @@ Temporary NN t-dependence for all hadrons
596	if(pPDG>3400 \|\| pPDG<-3400) G4cout<<"-Warning-G4QDifR::Fragment: pPDG="<<pPDG<<G4endl;
597	G4int PDG=2212; // TMP instead of pPDG
598	G4double xSec=CSmanager->GetCrossSection(false, P, tgZ, tgN, PDG);// Rec.CrossSect TMP
599	//G4double xSec=CSmanager->GetCrossSection(false, P, tgZ, tgN, pPDG); // Rec.CrossSect
600	#ifdef debug
601	G4cout<<"G4QDiffR::ProjFragment:pPDG="<<pPDG<<",P="<<P<<",CS="<<xSec/millibarn<<G4endl;
602	#endif
603	#ifdef nandebug
604	if(xSec>0. \|\| xSec<0. \|\| xSec==0);
605	else G4cout<<"*NAN*G4QDiffRatio::ProjFragment: xSec="<<xSec/millibarn<<G4endl;
606	#endif
607	// @@ check a possibility to separate p, n, or alpha (!)
608	if(xSec <= 0.) // The cross-section iz 0 -> Do Nothing
609	{
610	#ifdef pdebug
611	G4cerr<<"-Warning-G4QDiffRatio::ProjFragment:Zero XSPDG="<<pPDG<<",P="<<P<<G4endl;
612	#endif
613	return ResHV; //Do Nothing Action
614	}
615	G4double t=CSmanager->GetExchangeT(tgZ,tgN,pPDG); // functional randomized -t (MeV^2)
616	G4double maxt=(dsds-4mT2*mDif2)/s; // maximum possible -t
617	#ifdef pdebug
618	G4cout<<"G4QDifR::PFra:ph="<<pPDG<<",P="<<P<<",X="<<xSec<<",t="<<mint<<"<"<<maxt<<G4endl;
619	#endif
620	#ifdef nandebug
621	if(mint>-.0000001); // To make the Warning for NAN
622	else G4cout<<"******G4QDiffractionRatio::ProjFragment: -t="<<mint<<G4endl;
623	#endif
624	G4double rt=t/maxt;
625	G4double cost=1.-rt-rt; // cos(theta) in CMS
626	#ifdef ppdebug
627	G4cout<<"G4QDiffRatio::ProjFragment: -t="<<t<<", maxt="<<maxt<<", cost="<<cost<<G4endl;
628	#endif
629	if(cost>1. \|\| cost<-1. \|\| !(cost>-1. \|\| cost<=1.))
630	{
631	if (cost>1.) cost=1.;
632	else if(cost<-1.) cost=-1.;
633	else
634	{
635	G4cerr<<"G4QDiffRat::ProjFragm: NAN cost="<<cost<<",t="<<t<<",tmax="<<maxt<<G4endl;
636	return ResHV; // Do Nothing Action
637	}
638	}
639	G4LorentzVector r4M=G4LorentzVector(0.,0.,0.,mT); // 4mom of the recoil nucleon
640	G4LorentzVector d4M=G4LorentzVector(0.,0.,0.,mDif); // 4mom of the diffract. Quasmon
641	G4LorentzVector dir4M=tot4M-G4LorentzVector(0.,0.,0.,(tot4M.e()-mT)*.01);
642	if(!G4QHadron(tot4M).RelDecayIn2(d4M, r4M, dir4M, cost, cost))
643	{
644	G4cerr<<"G4QDifR::PFr:M="<<tot4M.m()<<",T="<<mT<<",D="<<mDif<<",T+D="<<mT+mDif<<G4endl;
645	//G4Exception("G4QDifR::Fragm:","009",FatalException,"Decay of ElasticComp");
646	return ResHV; // Do Nothing Action
647	}
648	#ifdef debug
649	G4cout<<"G4QDiffR::ProjFragm:d4M="<<d4M<<"+r4M="<<r4M<<"="<<d4M+r4M<<"="<<tot4M<<G4endl;
650	#endif
651	// Now everything is ready for fragmentation and DoNothing projHadron must be wiped out
652	ResHV->pop_back(); // Clean up pointer to the fake (doNothing) projectile
653	delete hadron; // Delete the fake (doNothing) projectile hadron
654	hadron = new G4QHadron(tPDG,t4M); // Hadron for the recoil neucleus
655	ResHV->push_back(hadron); // Fill the recoil nucleus
656	#ifdef debug
657	G4cout<<"G4QDiffractionRatio::ProjFragment: Filled RecNucleus="<<t4M<<tPDG<<G4endl;
658	#endif
659	hadron = new G4QHadron(rPDG,r4M); // Hadron for the recoil nucleon
660	ResHV->push_back(hadron); // Fill the recoil nucleon
661	#ifdef debug
662	G4cout<<"G4QDiffractionRatio::ProjFragment: Filled RecNucleon="<<r4M<<rPDG<<G4endl;
663	#endif
664	G4LorentzVector sum4M(0.,0.,0.,0.);
665	// Now the (pPdg,d4M) Quasmon must be fragmented
666	G4QHadronVector* leadhs = new G4QHadronVector;// Prototype of QuasmOutput G4QHadronVector
667	G4QContent dQC=G4QPDGCode(pPDG).GetQuarkContent(); // Quark Content of the projectile
668	G4Quasmon* pan= new G4Quasmon(dQC,d4M); // --->---->---->----->-----> DELETED -->---*
669	try // \|
670	{ // \|
671	G4QNucleus vac(90000000); // \|
672	leadhs=pan->Fragment(vac,1); // DELETED after it is copied to ResHV vector -->---+-*
673	} // \| \|
674	catch (G4QException& error) // \| \|
675	{ // \| \|
676	G4cerr<<"***G4QDiffractionRatio::ProjFragment: G4Quasmon Exception"<<G4endl; //\| \|
677	G4Exception("G4QDiffractionRatio::ProjFragment","72",FatalException,"*Quasmon");//\| \|
678	} // \| \|
679	delete pan; // Delete the Nuclear Environment <----<---* \|
680	G4int qNH=leadhs->size(); // A#of collected hadrons from diff.frag. \|
681	if(qNH) for(G4int iq=0; iq<qNH; iq++) // Loop over hadrons to fill the result \|
682	{ // \|
683	G4QHadron* loh=(*leadhs)[iq]; // Pointer to the output hadron \|
684	G4int nL=loh->GetStrangeness(); // A number of Lambdas in the Hypernucleus \|
685	G4int nB=loh->GetBaryonNumber(); // Total Baryon Number of the Hypernucleus \|
686	G4int nC = loh->GetCharge(); // Charge of the Hypernucleus \|
687	G4int oPDG = loh->GetPDGCode(); // Original CHIPS PDG Code of the hadron \|
688	//if((nC>nB \|\| nC<0) && nB>0 && nL>=0 && nL<=nB && oPDG>80000000) // Iso-nucleus \|
689	if(2>3) // Closed because "G4QDR::F:90002999,M=-7.768507e-04,B=2,S=0,C=3" is found \|
690	{
691	G4LorentzVector q4M = loh->Get4Momentum(); // Get 4-momentum of the Isonucleus \|
692	G4double qM=q4M.m(); // Real mass of the Isonucleus
693	#ifdef fdebug
694	G4cout<<"G4QDR::PF:"<<oPDG<<",M="<<qM<<",B="<<nB<<",S="<<nL<<",C="<<nC<<G4endl;// \|
695	#endif
696	G4int qPN=nC-nB; // Number of pions in the Isonucleus \|
697	G4int fPDG = 2212; // Prototype for nP+(Pi+) case \|
698	G4int sPDG = 211;
699	G4int tPDG = 3122; // @@ Sigma0 (?) \|
700	G4double fMass= mProt;
701	G4double sMass= mPi;
702	G4double tMass= mLamb; // @@ Sigma0 (?) \|
703	G4bool cont=true; // Continue flag \|
704	// ========= Negative state ======
705	if(nC<0) // ====== Only Pi- can help \|
706	{
707	if(nL&&nB==nL) // --- nLamb + k(Pi-) State --- \|
708	{
709	sPDG = -211;
710	if(-nC==nL && nL==1) // Only one Sigma- like (nB=1) \|
711	{
712	if(std::fabs(qM-mSigM)<eps)
713	{
714	loh->SetQPDG(G4QPDGCode(3112)); // This is Sigma- \|
715	cont=false; // Skip decay \|
716	}
717	else if(qM>mLamb+mPi) //(2) Sigma- => Lambda + Pi- decay \|
718	{
719	fPDG = 3122;
720	fMass= mLamb;
721	}
722	else if(qM>mSigM) //(2) Sigma+=>Sigma++gamma decay \|
723	{
724	fPDG = 3112;
725	fMass= mSigM;
726	sPDG = 22;
727	sMass= 0.;
728	}
729	else //(2) Sigma-=>Neutron+Pi- decay \|
730	{
731	fPDG = 2112;
732	fMass= mNeut;
733	}
734	qPN = 1; // #of (Pi+ or gamma)'s = 1 \|
735	}
736	else if(-nC==nL) //(2) a few Sigma- like \|
737	{
738	qPN = 1; // One separated Sigma- \|
739	fPDG = 3112;
740	sPDG = 3112;
741	sMass= mSigM;
742	nB--;
743	fMass= mSigM;
744	}
745	else if(-nC>nL) //(2) n(Sigma-)+m(Pi-) \|
746	{
747	qPN = -nC-nL; // #of Pi-'s \|
748	fPDG = 3112;
749	fMass= mSigM;
750	}
751	else //(2) n(Sigma-)+mLambda(-nC<nL) \|
752	{
753	nB += nC; // #of Lambda's \|
754	fPDG = 3122;
755	fMass= mLamb;
756	qPN = -nC; // #of Sigma+'s \|
757	sPDG = 3112;
758	sMass= mSigM;
759	}
760	nL = 0; // Only decays in two are above \|
761	}
762	else if(nL) // ->nLamb+mNeut+k*(Pi-) State (nL<nB) \|
763	{
764	nB -= nL; // #of neutrons \|
765	fPDG = 2112;
766	fMass= mNeut;
767	G4int nPin = -nC; // #of Pi-'s
768	if(nL==nPin) //(2) mNeut+nSigma- \|
769	{
770	qPN = nL; // #of Sigma- \|
771	sPDG = 3112;
772	sMass= mSigM;
773	nL = 0;
774	}
775	else if(nL>nPin) //(3) mP+n(Sigma+)+k*Lambda \|
776	{
777	nL-=nPin; // #of Lambdas \|
778	qPN = nPin; // #of Sigma+ \|
779	sPDG = 3112;
780	sMass= mSigM;
781	}
782	else //(3) mN+n(Sigma-)+k*(Pi-) (nL<nPin) \|
783	{
784	qPN = nPin-nL; // #of Pi- \|
785	sPDG = -211;
786	tPDG = 3112;
787	tMass= mSigM;
788	}
789	}
790	else //(2) nN+m(Pi-) (nL=0) \|
791	{
792	sPDG = -211;
793	qPN = -nC;
794	fPDG = 2112;
795	fMass= mNeut;
796	}
797	}
798	else if(!nC) // * Should not be here * \|
799	{
800	if(nL && nL<nB) //(2) nLamb+mN *Should not be here* \|
801	{
802	qPN = nL;
803	fPDG = 2112; // mN+nL case \|
804	sPDG = 3122;
805	sMass= mLamb;
806	nB -= nL;
807	fMass= mNeut;
808	nL = 0;
809	}
810	else if(nL>1 && nB==nL) //(2) mLamb(m>1) Should not be here* \|
811	{
812	qPN = 1;
813	fPDG = 3122;
814	sPDG = 3122;
815	sMass= mLamb;
816	nB--;
817	fMass= mLamb;
818	}
819	else if(!nL && nB>1) //(2) nNeut(n>1) Should not be here* \|
820	{
821	qPN = 1;
822	fPDG = 2112;
823	sPDG = 2112;
824	sMass= mNeut;
825	nB--;
826	fMass= mNeut;
827	}
828	else G4cout<<"?G4QDiffractionRatio::ProjFragment: (1) oPDG="<<oPDG<<G4endl;// \|
829	}
830	else if(nC>0) // nLamb+(mP)+(k*Pi+) \|
831	{
832	if(nL && nL+nC==nB) //(2) nLamb+mP *Should not be here* \|
833	{
834	qPN = nL;
835	nL = 0;
836	fPDG = 2212;
837	sPDG = 3122;
838	sMass= mLamb;
839	nB = nC;
840	fMass= mProt;
841	}
842	else if(nL && nC<nB-nL) //(3)nL+mP+kN Should not be here* \|
843	{
844	qPN = nC; // #of protons \|
845	fPDG = 2112; // mP+nL case \|
846	sPDG = 2212;
847	sMass= mProt;
848	nB -= nL+nC; // #of neutrons \|
849	fMass= mNeut;
850	}
851	else if(nL && nB==nL) // ---> nL+mPi+ State \|
852	{
853	if(nC==nL && nL==1) // Only one Sigma+ like State \|
854	{
855	if(std::fabs(qM-mSigP)<eps)
856	{
857	loh->SetQPDG(G4QPDGCode(3222)); // This is GS Sigma+ \|
858	cont=false; // Skip decay \|
859	}
860	else if(qM>mLamb+mPi) //(2) Sigma+=>Lambda+Pi+ decay \|
861	{
862	fPDG = 3122;
863	fMass= mLamb;
864	}
865	else if(qM>mNeut+mPi) //(2) Sigma+=>Neutron+Pi+ decay \|
866	{
867	fPDG = 2112;
868	fMass= mNeut;
869	}
870	else if(qM>mSigP) //(2) Sigma+=>Sigma++gamma decay \|
871	{
872	fPDG = 3222;
873	fMass= mSigP;
874	sPDG = 22;
875	sMass= 0.;
876	}
877	else //(2) Sigma+=>Proton+gamma decay \|
878	{
879	fPDG = 2212;
880	fMass= mProt;
881	sPDG = 22;
882	sMass= 0.;
883	}
884	qPN = 1; // #of (Pi+ or gamma)'s = 1 \|
885	}
886	else if(nC==nL) //(2) a few Sigma+ like hyperons \|
887	{
888	qPN = 1;
889	fPDG = 3222;
890	sPDG = 3222;
891	sMass= mSigP;
892	nB--;
893	fMass= mSigP;
894	}
895	else if(nC>nL) //(2) n(Sigma+)+m(Pi+) \|
896	{
897	qPN = nC-nL; // #of Pi+'s \|
898	fPDG = 3222;
899	nB = nL; // #of Sigma+'s \|
900	fMass= mSigP;
901	}
902	else //(2) n(Sigma+)+mLambda \|
903	{
904	nB -= nC; // #of Lambda's \|
905	fPDG = 3122;
906	fMass= mLamb;
907	qPN = nC; // #of Sigma+'s \|
908	sPDG = 3222;
909	sMass= mSigP;
910	}
911	nL = 0; // All above are decays in 2 \|
912	}
913	else if(nL && nC>nB-nL) // nLamb+mP+k*Pi+ \|
914	{
915	nB -= nL; // #of protons \|
916	G4int nPip = nC-nB; // #of Pi+'s \|
917	if(nL==nPip) //(2) mP+nSigma+ \|
918	{
919	qPN = nL; // #of Sigma+ \|
920	sPDG = 3222;
921	sMass= mSigP;
922	nL = 0;
923	}
924	else if(nL>nPip) //(3) mP+n(Sigma+)+k*Lambda \|
925	{
926	nL -= nPip; // #of Lambdas \|
927	qPN = nPip; // #of Sigma+ \|
928	sPDG = 3222;
929	sMass= mSigP;
930	}
931	else //(3) mP+n(Sigma+)+k*(Pi+) \|
932	{
933	qPN = nPip-nL; // #of Pi+ \|
934	tPDG = 3222;
935	tMass= mSigP;
936	}
937	}
938	if(nC<nB) //(2) nP+mN *Should not be here* \|
939	{
940	fPDG = 2112;
941	fMass= mNeut;
942	qPN = nC;
943	sPDG = 2212;
944	sMass= mProt;
945	}
946	else if(nB==nC && nC>1) //(2) mProt(m>1) Should not be here* \|
947	{
948	qPN = 1;
949	fPDG = 2212;
950	sPDG = 2212;
951	sMass= mProt;
952	nB--;
953	fMass= mProt;
954	}
955	else if(nC<=nB\|\|!nB) G4cout<<"?G4QDR::ProjFragm: (2) oPDG="<<oPDG<<G4endl; // \|
956	// !nL && nC>nB //(2) Default condition nP+m(Pi+) \|
957	}
958	if(cont) // Make a decay \|
959	{
960	G4double tfM=nB*fMass;
961	G4double tsM=qPN*sMass;
962	G4double ttM=0.;
963	if(nL) ttM=nL*tMass;
964	G4LorentzVector f4Mom(0.,0.,0.,tfM);
965	G4LorentzVector s4Mom(0.,0.,0.,tsM);
966	G4LorentzVector t4Mom(0.,0.,0.,ttM);
967	G4double sum=tfM+tsM+ttM;
968	if(std::fabs(qM-sum)<eps)
969	{
970	f4Mom=q4M*(tfM/sum);
971	s4Mom=q4M*(tsM/sum);
972	if(nL) t4Mom=q4M*(ttM/sum);
973	}
974	else if(!nL && (qM<sum \|\| !G4QHadron(q4M).DecayIn2(f4Mom, s4Mom))) // Error \|
975	{
976	//#ifdef fdebug
977	G4cout<<"**G4QDR::PrFragm:fPDG="<<fPDG<<""<<nB<<"(fM="<<fMass<<")+sPDG="<<sPDG
978	<<"*"<<qPN<<"(sM="<<sMass<<")"<<"="<<sum<<" > TM="<<qM<<q4M<<oPDG<<G4endl;
979	//#endif
980	throw G4QException("*G4QDiffractionRatio::ProjFragment: Bad decay in 2"); // \|
981	}
982	else if(nL && (qM<sum \|\| !G4QHadron(q4M).DecayIn3(f4Mom, s4Mom, t4Mom)))// Error\|
983	{
984	//#ifdef fdebug
985	G4cout<<"**G4DF::PrFrag: "<<fPDG<<""<<nB<<"("<<fMass<<")+"<<sPDG<<"*"<<qPN<<"("
986	<<sMass<<")+Lamb*"<<nL<<"="<<sum<<" > TotM="<<qM<<q4M<<oPDG<<G4endl;
987	//#endif
988	throw G4QException("*G4QDiffractionRatio::ProjFragment: Bad decay in 3"); // \|
989	}
990	#ifdef fdebug
991	G4cout<<"G4QDF::ProjFragm: DONE n="<<nB<<f4Mom<<fPDG<<", m="<<qPN<<s4Mom<<sPDG
992	<<", l="<<nL<<t4Mom<<G4endl;
993	#endif
994	G4bool notused=true;
995	if(nB) // There are baryons \|
996	{
997	f4Mom/=nB;
998	loh->Set4Momentum(f4Mom); // ! Update the Hadron ! \|
999	loh->SetQPDG(G4QPDGCode(fPDG)); // Baryons \|
1000	notused=false; // Loh was used \|
1001	if(nB>1) for(G4int ih=1; ih<nB; ih++) // Loop over the rest of baryons \|
1002	{
1003	G4QHadron* Hi = new G4QHadron(fPDG,f4Mom); // Create a Hadron for Baryon \|
1004	ResHV->push_back(Hi); // Fill in the additional nucleon \|
1005	#ifdef fdebug
1006	sum4M+=r4M; // Sum 4-momenta for the EnMom check \|
1007	G4cout<<"G4QDR::ProjFrag: additional Nucleon="<<f4Mom<<fPDG<<G4endl; // \|
1008	#endif
1009	}
1010	}
1011	if(qPN) // There are pions \|
1012	{
1013	s4Mom/=qPN;
1014	G4int min=0;
1015	if(notused)
1016	{
1017	loh->Set4Momentum(s4Mom); // ! Update the Hadron 4M ! \|
1018	loh->SetQPDG(G4QPDGCode(sPDG)); // Update PDG \|
1019	notused=false; // loh was used \|
1020	min=1; // start value \|
1021	}
1022	if(qPN>min) for(G4int ip=min; ip<qPN; ip++) // Loop over pions \|
1023	{
1024	G4QHadron* Hj = new G4QHadron(sPDG,s4Mom); // Create a Hadron for the meson \|
1025	ResHV->push_back(Hj); // Fill in the additional pion \|
1026	#ifdef fdebug
1027	sum4M+=r4M; // Sum 4-momenta for the EnMom check \|
1028	G4cout<<"G4QDR::ProjFragm: additional Pion="<<f4Mom<<fPDG<<G4endl; // \|
1029	#endif
1030	}
1031	}
1032	if(nL) // There are Hyperons \|
1033	{
1034	t4Mom/=nL;
1035	G4int min=0;
1036	if(notused)
1037	{
1038	loh->Set4Momentum(t4Mom); // ! Update the Hadron 4M ! \|
1039	loh->SetQPDG(G4QPDGCode(tPDG));// Update PDG \|
1040	notused=false; // loh was used \|
1041	min=1; //
1042	}
1043	if(nL>min) for(G4int il=min; il<nL; il++) // Loop over Hyperons \|
1044	{
1045	G4QHadron* Hk = new G4QHadron(tPDG,t4Mom); // Create a Hadron for Lambda \|
1046	ResHV->push_back(Hk); // Fill in the additional pion \|
1047	#ifdef fdebug
1048	sum4M+=r4M; // Sum 4-momenta for the EnMom check \|
1049	G4cout<<"G4QDR::ProjFragm: additional Hyperon="<<f4Mom<<fPDG<<G4endl; // \|
1050	#endif
1051	}
1052	}
1053	} // --> End of decay \|
1054	} // -> End of Iso-nuclear treatment \|
1055	else if( (nL > 0 && nB > 1) \|\| (nL < 0 && nB < -1) )
1056	{ // Hypernucleus is found \|
1057	G4bool anti=false; // Default=Nucleus (true=antinucleus \|
1058	if(nB<0) // Anti-nucleus \|
1059	{
1060	anti=true; // Flag of anti-hypernucleus \|
1061	nB=-nB; // Reverse the baryon number \|
1062	nC=-nC; // Reverse the charge \|
1063	nL=-nL; // Reverse the strangeness \|
1064	}
1065	G4int hPDG = 90000000+nL999999+nC999+nB; // CHIPS PDG Code for Hypernucleus \|
1066	G4int nSM=0; // A#0f unavoidable Sigma- \|
1067	G4int nSP=0; // A#0f unavoidable Sigma+ \|
1068	if(nC<0) // Negative hypernucleus \|
1069	{
1070	if(-nC<=nL) // Partial compensation by Sigma- \|
1071	{
1072	nSM=-nC; // Can be compensated by Sigma- \|
1073	nL+=nC; // Reduce the residual strangeness \|
1074	}
1075	else // All Charge is compensated by Sigma- \|
1076	{
1077	nSM=nL; // The maximum number of Sigma- \|
1078	nL=0; // Kill the residual strangeness \|
1079	}
1080	}
1081	else if(nC>nB-nL) // Extra positive hypernucleus \|
1082	{
1083	if(nC<=nB) // Partial compensation by Sigma+ \|
1084	{
1085	G4int dH=nB-nC; // Isotopic shift \|
1086	nSP=nL-dH; // Can be compensated by Sigma+ \|
1087	nL=dH; // Reduce the residual strangeness \|
1088	}
1089	else // All Charge is compensated by Sigma+ \|
1090	{
1091	nSP=nL; // The maximum number of Sigma+ \|
1092	nL=0; // Kill the residual strangeness \|
1093	}
1094	}
1095	G4LorentzVector r4M=loh->Get4Momentum(); // Real 4-momentum of the hypernucleus !
1096	G4double reM=r4M.m(); // Real mass of the hypernucleus \|
1097	#ifdef fdebug
1098	G4cout<<"G4QDiffRatio::PrFrag:oPDG=="<<oPDG<<",hPDG="<<hPDG<<",M="<<reM<<G4endl;//\|
1099	#endif
1100	G4int rlPDG=hPDG-nL1000000-nSP1000999-nSM*999001;// Subtract Lamb/Sig from Nucl.\|
1101	G4int sPDG=3122; // Prototype for the Hyperon PDG (Lambda)\|
1102	G4double MLa=mLamb; // Prototype for one Hyperon decay \|
1103	#ifdef fdebug
1104	G4cout<<"G4QDiffRatio::PrFrag:G4nS+="<<nSP<<",nS-="<<nSM<<",nL="<<nL<<G4endl;// \|
1105	#endif
1106	if(nSP\|\|nSM) // Sigma+/- improvement \|
1107	{
1108	if(nL) // By mistake Lambda improvement is found \|
1109	{
1110	G4cout<<"***G4QDR::PFr:HypN="<<hPDG<<": bothSigm&Lamb -> ImproveIt"<<G4endl;//\|
1111	//throw G4QException("*G4QDiffractionRatio::Fragment:BothLambda&SigmaInHN");//\|
1112	// @@ Correction, which does not conserv the charge !! (-> add decay in 3) \|
1113	if(nSP) nL+=nSP; // Convert Sigma+ to Lambda \|
1114	else nL+=nSM; // Convert Sigma- to Lambda \|
1115	}
1116	if(nSP) // Sibma+ should be decayed \|
1117	{
1118	nL=nSP; // #of decaying hyperons \|
1119	sPDG=3222; // PDG code of decaying hyperons \|
1120	MLa=mSigP; // Mass of decaying hyperons \|
1121	}
1122	else // Sibma+ should be decayed \|
1123	{
1124	nL=nSM; // #of decaying hyperons \|
1125	sPDG=3112; // PDG code of decaying hyperons \|
1126	MLa=mSigM; // Mass of decaying hyperons \|
1127	}
1128	}
1129	#ifdef fdebug
1130	G4cout<<"G4QDiffRat::ProjFrag:G4mS="<<MLa<<",sPDG="<<sPDG<<",nL="<<nL<<G4endl;//\|
1131	#endif
1132	if(nL>1) MLa*=nL; // Total mass of the decaying hyperons \|
1133	G4double rlM=G4QNucleus(rlPDG).GetMZNS();// Mass of the NonstrangeNucleus \|
1134	if(!nSP&&!nSM&&nL==1&&reM>rlM+mSigZ&&G4UniformRand()>.5) // Conv Lambda->Sigma0 \|
1135	{
1136	sPDG=3212; // PDG code of a decaying hyperon \|
1137	MLa=mSigZ; // Mass of the decaying hyperon \|
1138	}
1139	G4int rnPDG = hPDG-nL*999999; // Convert Lambdas to neutrons (for convInN) \|
1140	G4QNucleus rnN(rnPDG); // New nonstrange nucleus \|
1141	G4double rnM=rnN.GetMZNS(); // Mass of the new nonstrange nucleus \|
1142	// @@ In future take into account Iso-Hypernucleus (Add PI+,R & Pi-,R decays) \|
1143	if(rlPDG==90000000) // Multy Hyperon (HyperNuc of only hyperons) \|
1144	{
1145	if(nL>1) r4M=r4M/nL; // split the 4-mom for the MultyLambda \|
1146	for(G4int il=0; il<nL; il++) // loop over Lambdas \|
1147	{
1148	if(anti) sPDG=-sPDG; // For anti-nucleus case \|
1149	G4QHadron* theLam = new G4QHadron(sPDG,r4M); // Make NewHadr for the Hyperon \|
1150	ResHV->push_back(theLam); // Fill in the Lambda \|
1151	#ifdef fdebug
1152	sum4M+=r4M; // Sum 4-momenta for the EnMom check \|
1153	G4cout<<"G4QDR::ProjFrag: additional Lambda="<<r4M<<sPDG<<G4endl; // \|
1154	#endif
1155	}
1156	}
1157	else if(reM>rlM+MLa-eps) // Lambda (or Sigma) can be split \|
1158	{
1159	G4LorentzVector n4M(0.,0.,0.,rlM); // 4-mom of the residual nucleus \|
1160	G4LorentzVector h4M(0.,0.,0.,MLa); // 4-mom of the Hyperon \|
1161	G4double sum=rlM+MLa; // Safety sum \|
1162	if(std::fabs(reM-sum)<eps) // At rest in CMS \|
1163	{
1164	n4M=r4M*(rlM/sum); // Split tot 4-mom for resNuc \|
1165	h4M=r4M*(MLa/sum); // Split tot 4-mom for Hyperon \|
1166	}
1167	else if(reM<sum \|\| !G4QHadron(r4M).DecayIn2(n4M,h4M)) // Error in decay \|
1168	{
1169	G4cerr<<"**G4QDF::PF:HypN,M="<<reM<<"<A+nL="<<sum<<",d="<<sum-reM<<G4endl;//\|
1170	throw G4QException("***G4QDiffractionRatio::ProjFragment:HypernuclusDecay");//\|
1171	}
1172	#ifdef fdebug
1173	G4cout<<"G4QDR::PF:HypN="<<r4M<<"->A="<<rlPDG<<n4M<<",nL="<<nL<<h4M<<G4endl;//\|
1174	#endif
1175	loh->Set4Momentum(n4M); // ! Update the Hadron ! \|
1176	if(anti && rlPDG==90000001) rlPDG=-2112; // Convert to anti-neutron \|
1177	if(anti && rlPDG==90001000) rlPDG=-2212; // Convert to anti-proton \|
1178	loh->SetQPDG(G4QPDGCode(rlPDG)); // ConvertedHypernucleus to nonstrange(@anti)\|
1179	if(rlPDG==90000002) // Additional action with loH changed to 2n \|
1180	{
1181	G4LorentzVector newLV=n4M/2.; // Split 4-momentum \|
1182	loh->Set4Momentum(newLV); // Reupdate the hadron \|
1183	if(anti) loh->SetQPDG(G4QPDGCode(-2112)); // Make anti-neutron PDG \|
1184	else loh->SetQPDG(G4QPDGCode(2112)); // Make neutron PDG \|
1185	G4QHadron* secHadr = new G4QHadron(loh); // Duplicate the neutron \|
1186	ResHV->push_back(secHadr); // Fill in the additional neutron \|
1187	#ifdef fdebug
1188	sum4M+=r4M; // Sum 4-momenta for the EnMom check \|
1189	G4cout<<"G4QDR::ProgFrag: additional Neutron="<<r4M<<sPDG<<G4endl; // \|
1190	#endif
1191	}
1192	else if(rlPDG==90002000) // Additional action with loH change to 2p \|
1193	{
1194	G4LorentzVector newLV=n4M/2.; // Split 4-momentum \|
1195	loh->Set4Momentum(newLV); // Reupdate the hadron \|
1196	if(anti) loh->SetQPDG(G4QPDGCode(-2212)); // Make anti-neutron PDG \|
1197	else loh->SetQPDG(G4QPDGCode(2112)); // Make neutron PDG \|
1198	G4QHadron* secHadr = new G4QHadron(loh); // Duplicate the proton \|
1199	ResHV->push_back(secHadr); // Fill in the additional neutron \|
1200	#ifdef fdebug
1201	sum4M+=r4M; // Sum 4-momenta for the EnMom check \|
1202	G4cout<<"G4QDR::ProjFrag: additional Proton="<<r4M<<sPDG<<G4endl; // \|
1203	#endif
1204	}
1205	// @@(?) Add multybaryon decays if necessary (Now it anyhow is made later) \|
1206	#ifdef fdebug
1207	G4cout<<"*G4QDiffractionRatio::PrFrag:resNucPDG="<<loh->GetPDGCode()<<G4endl;// \|
1208	#endif
1209	if(nL>1) h4M=h4M/nL; // split the lambda's 4-mom if necessary \|
1210	for(G4int il=0; il<nL; il++) // A loop over excessive hyperons \|
1211	{
1212	if(anti) sPDG=-sPDG; // For anti-nucleus case \|
1213	G4QHadron* theLamb = new G4QHadron(sPDG,h4M); // Make NewHadr for the Hyperon \|
1214	ResHV->push_back(theLamb); // Fill in the additional neutron \|
1215	#ifdef fdebug
1216	sum4M+=r4M; // Sum 4-momenta for the EnMom check \|
1217	G4cout<<"G4QDR::ProjFrag: additional Hyperon="<<r4M<<sPDG<<G4endl; // \|
1218	#endif
1219	}
1220	}
1221	else if(reM>rnM+mPi0-eps&&!nSP&&!nSM)// Lambda->N only if Sigmas are absent \|
1222	{
1223	G4int nPi=static_cast<G4int>((reM-rnM)/mPi0); // Calc. pion multiplicity \|
1224	if (nPi>nL) nPi=nL; // Cut the pion multiplicity \|
1225	G4double npiM=nPi*mPi0; // Total pion mass \|
1226	G4LorentzVector n4M(0.,0.,0.,rnM); // Residual nucleus 4-momentum \|
1227	G4LorentzVector h4M(0.,0.,0.,npiM);// 4-momentum of pions \|
1228	G4double sum=rnM+npiM; // Safety sum \|
1229	if(std::fabs(reM-sum)<eps) // At rest \|
1230	{
1231	n4M=r4M*(rnM/sum); // The residual nucleus part \|
1232	h4M=r4M*(npiM/sum); // The pion part \|
1233	}
1234	else if(reM<sum \|\| !G4QHadron(r4M).DecayIn2(n4M,h4M)) // Error in decay \|
1235	{
1236	G4cerr<<"G4QDR::PF:HypN,M="<<reM<<"<A+nPi0="<<sum<<",d="<<sum-reM<<G4endl;//\|
1237	throw G4QException("***G4QDiffractionRatio::ProjFragment:HypernuclDecay"); // \|
1238	}
1239	loh->Set4Momentum(n4M); // ! Update the Hadron ! \|
1240	if(anti && rnPDG==90000001) rnPDG=-2112; // Convert to anti-neutron \|
1241	if(anti && rnPDG==90001000) rnPDG=-2212; // Convert to anti-proton \|
1242	loh->SetQPDG(G4QPDGCode(rnPDG)); // convert hyperNuc to nonstrangeNuc(@@anti) \|
1243	#ifdef fdebug
1244	G4cout<<"G4QDR::PF:R="<<r4M<<"->A="<<rnPDG<<n4M<<",nPi0="<<nPi<<h4M<<G4endl;//\|
1245	#endif
1246	if(nPi>1) h4M=h4M/nPi; // Split the 4-mom if necessary \|
1247	for(G4int ihn=0; ihn<nPi; ihn++) // A loop over additional pions \|
1248	{
1249	G4QHadron* thePion = new G4QHadron(111,h4M); // Make a New Hadr for the pi0 \|
1250	ResHV->push_back(thePion); // Fill in the Pion \|
1251	#ifdef fdebug
1252	sum4M+=r4M; // Sum 4-momenta for the EnMom check \|
1253	G4cout<<"G4QDR::ProjFrag: additional Pion="<<r4M<<sPDG<<G4endl; // \|
1254	#endif
1255	}
1256	if(rnPDG==90000002) // Additional action with loH change to 2n \|
1257	{
1258	G4LorentzVector newLV=n4M/2.; // Split 4-momentum \|
1259	loh->Set4Momentum(newLV); // Reupdate the hadron \|
1260	if(anti) loh->SetQPDG(G4QPDGCode(-2112)); // Make anti-neutron PDG \|
1261	else loh->SetQPDG(G4QPDGCode(2112)); // Make neutron PDG \|
1262	G4QHadron* secHadr = new G4QHadron(loh); // Duplicate the neutron \|
1263	ResHV->push_back(secHadr); // Fill in the additional neutron \|
1264	#ifdef fdebug
1265	sum4M+=r4M; // Sum 4-momenta for the EnMom check \|
1266	G4cout<<"G4QDR::ProjFrag: additional Neutron="<<r4M<<sPDG<<G4endl; // \|
1267	#endif
1268	}
1269	else if(rnPDG==90002000) // Additional action with loH change to 2p \|
1270	{
1271	G4LorentzVector newLV=n4M/2.; // Split 4-momentum \|
1272	loh->Set4Momentum(newLV); // Reupdate the hadron \|
1273	if(anti) loh->SetQPDG(G4QPDGCode(-2212)); // Make anti-neutron PDG \|
1274	else loh->SetQPDG(G4QPDGCode(2112)); // Make neutron PDG \|
1275	G4QHadron* secHadr = new G4QHadron(loh); // Duplicate the proton \|
1276	ResHV->push_back(secHadr); // Fill in the additional neutron \|
1277	#ifdef fdebug
1278	sum4M+=r4M; // Sum 4-momenta for the EnMom check \|
1279	G4cout<<"G4QDR::ProjFrag: additional Proton="<<r4M<<sPDG<<G4endl; // \|
1280	#endif
1281	}
1282	// @@ Add multybaryon decays if necessary \|
1283	}
1284	else // If this Excepton shows up (lowProbable appearance) => include gamma decay \|
1285	{
1286	G4double d=rlM+MLa-reM; // Hyperon Excessive energy \|
1287	G4cerr<<"G4QDR::PF:R="<<rlM<<",S+="<<nSP<<",S-="<<nSM<<",L="<<nL<<",d="<<d<<G4endl;
1288	d=rnM+mPi0-reM; // Pion Excessive energy \|
1289	G4cerr<<"G4QDR::PF:"<<oPDG<<","<<hPDG<<",M="<<reM<<"<"<<rnM+mPi0<<",d="<<d<<G4endl;
1290	throw G4QException("G4QDiffractionRatio::ProjFragment: Hypernuclear conver");// \|
1291	}
1292	} // => End of G4 Hypernuclear decay \|
1293	ResHV->push_back(loh); // Fill in the result \|
1294	#ifdef debug
1295	sum4M+=loh->Get4Momentum(); // Sum 4-momenta for the EnMom check \|
1296	G4cout<<"G4QDR::PrFra:#"<<iq<<","<<loh->Get4Momentum()<<loh->GetPDGCode()<<G4endl;//\|
1297	#endif
1298	} // \|
1299	delete leadhs; // <----<----<----<----<----<----<----<----<----<----<----<----<----<--*
1300	#ifdef debug
1301	G4cout<<"G4QDiffractionRatio::ProjFragment: End Sum="<<sum4M<<" =?= d4M="<<d4M<<G4endl;
1302	#endif
1303	return ResHV; // Result
1304	} // End of ProjFragment
1305
1306	// Calculates Single Diffraction Taarget Excitation Cross-Section (independent Units)
1307	G4double G4QDiffractionRatio::GetTargSingDiffXS(G4double pIU, G4int pPDG, G4int Z, G4int N)
1308	{
1309	G4double mom=pIU/gigaelectronvolt; // Projectile momentum in GeV
1310	if ( mom < 1. \|\| (pPDG != 2212 && pPDG != 2112) )
1311	G4cerr<<"G4QDiffractionRatio::GetTargSingDiffXS isn't applicable p="<<mom<<" GeV, PDG="
1312	<<pPDG<<G4endl;
1313	G4double A=Z+N; // A of the target
1314	//return 4.5std::pow(A,.364)millibarn; // Result
1315	return 3.7std::pow(A,.364)millibarn; // Result after mpi0 correction
1316
1317	} // End of ProjFragment

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

Download in other formats: