Context Navigation

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

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

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

Download in other formats: