Context Navigation

source: trunk/source/processes/hadronic/models/chiral_inv_phase_space/cross_sections/src/G4QDiffractionRatio.cc @ 1199

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

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

Download in other formats: