Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
G4QuasiFreeRatios.cc
Go to the documentation of this file.
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$
28//
29//
30// G4 Physics class: G4QuasiFreeRatios for N+A elastic cross sections
31// Created: M.V. Kossov, CERN/ITEP(Moscow), 10-OCT-01
32// The last update: M.V. Kossov, CERN/ITEP (Moscow) 28-Oct-11
33//
34// ----------------------------------------------------------------------
35// Short description: Provides percentage of quasi-free and quasi-elastic
36// reactions in the inelastic reactions.
37// ----------------------------------------------------------------------
38
39//#define debug
40//#define pdebug
41//#define ppdebug
42//#define nandebug
43
44#include "G4QuasiFreeRatios.hh"
45#include "G4SystemOfUnits.hh"
46
47// initialisation of statics
48std::vector<G4double*> G4QuasiFreeRatios::vT; // Vector of pointers to LinTable in C++ heap
49std::vector<G4double*> G4QuasiFreeRatios::vL; // Vector of pointers to LogTable in C++ heap
50
51G4QuasiFreeRatios::G4QuasiFreeRatios()
52{
53#ifdef pdebug
54 G4cout<<"***^^^*** G4QuasiFreeRatios singletone is created ***^^^***"<<G4endl;
55#endif
56 QFreeScat=G4QFreeScattering::GetPointer();
57}
58
59G4QuasiFreeRatios::~G4QuasiFreeRatios()
60{
61 std::vector<G4double*>::iterator pos;
62 for(pos=vT.begin(); pos<vT.end(); ++pos) delete [] *pos;
63 vT.clear();
64 for(pos=vL.begin(); pos<vL.end(); ++pos) delete [] *pos;
65 vL.clear();
66}
67
68// Returns Pointer to the G4VQCrossSection class
69G4QuasiFreeRatios* G4QuasiFreeRatios::GetPointer()
70{
71 static G4QuasiFreeRatios theRatios; // *** Static body of the QEl Cross Section ***
72 return &theRatios;
73}
74
75// Calculation of pair(QuasiFree/Inelastic,QuasiElastic/QuasiFree)
76std::pair<G4double,G4double> G4QuasiFreeRatios::GetRatios(G4double pIU, G4int pPDG,
77 G4int tgZ, G4int tgN)
78{
79#ifdef pdebug
80 G4cout<<">>>IN>>>G4QFRat::GetQF:P="<<pIU<<",pPDG="<<pPDG<<",Z="<<tgZ<<",N="<<tgN<<G4endl;
81#endif
82 G4double R=0.;
83 G4double QF2In=1.; // Prototype of QuasiFree/Inel ratio for hN_tot
84 G4int tgA=tgZ+tgN;
85 if(tgA<2) return std::make_pair(QF2In,R); // No quasi-elastic on the only nucleon
86 std::pair<G4double,G4double> ElTot=GetElTot(pIU, pPDG, tgZ, tgN); // mean (IU)
87 //if( ( (pPDG>999 && pIU<227.) || pIU<27.) && tgA>1) R=1.; // @@ TMP to accelerate @lowE
88 if(pPDG>999 && pIU<227. && tgZ+tgN>1) R=1.; // To accelerate @lowE
89 else if(ElTot.second>0.)
90 {
91 R=ElTot.first/ElTot.second; // El/Total ratio (does not depend on units
92 QF2In=GetQF2IN_Ratio(ElTot.second/millibarn, tgZ+tgN); // QuasiFree/Inelastic ratio
93 }
94#ifdef pdebug
95 G4cout<<">>>OUT>>>G4QuasiFreeRatio::GetQF2IN_Ratio: QF2In="<<QF2In<<", R="<<R<<G4endl;
96#endif
97 return std::make_pair(QF2In,R);
98}
99
100// Calculatio QasiFree/Inelastic Ratio as a function of total hN cross-section (mb) and A
101G4double G4QuasiFreeRatios::GetQF2IN_Ratio(G4double s_value, G4int A)
102{
103 static const G4int nps=150; // Number of steps in the R(s) LinTable
104 static const G4int mps=nps+1; // Number of elements in the R(s) LinTable
105 static const G4double sma=150.; // The first LinTabEl(s=0)=1., s>sma -> LogTable
106 static const G4double ds=sma/nps; // Step of the linear Table
107 static const G4int nls=100; // Number of steps in the R(lns) LogTable
108 static const G4int mls=nls+1; // Number of elements in the R(lns) LogTable
109 static const G4double lsi=5.; // The min ln(s) LogTabEl(s=148.4 < sma=150.)
110 static const G4double lsa=9.; // The max ln(s) LogTabEl(s=148.4 - 8103. mb)
111 static const G4double mi=std::exp(lsi);// The min s of LogTabEl(~ 148.4 mb)
112 static const G4double max_s=std::exp(lsa);// The max s of LogTabEl(~ 8103. mb)
113 static const G4double dl=(lsa-lsi)/nls;// A step Log-Length of the Logarithmic Table
114 static const G4double edl=std::exp(dl);// Multiplication step of the Logarithmic Table
115 static const G4double toler=.01; // The tolarence mb defining the same cross-section
116 static G4double lastS=0.; // The last sigma value for which R was calculated
117 static G4double lastR=0.; // The last ratio R which was calculated
118 // Local Associative Data Base:
119 static std::vector<G4int> vA; // Vector of calculated A
120 static std::vector<G4double> vH; // Vector of max s initialized in the LinTable
121 static std::vector<G4int> vN; // Vector of topBin number initialized in LinTable
122 static std::vector<G4double> vM; // Vector of rel max ln(s) initialized in LogTable
123 static std::vector<G4int> vK; // Vector of topBin number initialized in LogTable
124 // Last values of the Associative Data Base:
125 static G4int lastA=0; // theLast of calculated A
126 static G4double lastH=0.; // theLast of max s initialized in the LinTable
127 static G4int lastN=0; // theLast of topBin number initialized in LinTable
128 static G4double lastM=0.; // theLast of rel max ln(s) initialized in LogTable
129 static G4int lastK=0; // theLast of topBin number initialized in LogTable
130 static G4double* lastT=0; // theLast of pointer to LinTable in the C++ heap
131 static G4double* lastL=0; // theLast of pointer to LogTable in the C++ heap
132 // LogTable is created only if necessary. The ratio R(s>8100 mb) = 0 for any nuclei
133#ifdef pdebug
134 G4cout<<"+++G4QuasiFreeRatio::GetQF2IN_Ratio:A="<<A<<", s="<<s_value<<G4endl;
135#endif
136 if(s_value<toler || A<2) return 1.;
137 if(s_value>max_s) return 0.;
138 if(A>238)
139 {
140 G4cout<<"-Warning-G4QuasiFreeRatio::GetQF2IN_Ratio:A="<<A<<">238, return zero"<<G4endl;
141 return 0.;
142 }
143 G4int nDB=vA.size(); // A number of nuclei already initialized in AMDB
144 // if(nDB && lastA==A && std::fabs(s_value-lastS)<toler) return lastR;
145 if(nDB && lastA==A && s_value==lastS) return lastR; // VI do not use tolerance
146 G4bool found=false;
147 G4int i=-1;
148 if(nDB) for (i=0; i<nDB; i++) if(A==vA[i]) // Sirch for this A in AMDB
149 {
150 found=true; // The A value is found
151 break;
152 }
153#ifdef pdebug
154 G4cout<<"+++G4QuasiFreeRatio::GetQF2IN_Ratio: nDB="<<nDB<<", found="<<found<<G4endl;
155#endif
156 if(!nDB || !found) // Create new line in the AMDB
157 {
158 lastA = A;
159#ifdef pdebug
160 G4cout<<"G4QuasiFreeRatios::GetQF2IN_Ratio: NewT, A="<<A<<", nDB="<<nDB<<G4endl;
161#endif
162 lastT = new G4double[mps]; // Create the Linear Table
163 lastN = static_cast<int>(s_value/ds)+1; // MaxBin to be initialized in Lin Table
164 if(lastN>nps) // The Lin Table should be completely initialized
165 {
166 lastN=nps;
167 lastH=sma;
168 }
169 else lastH = lastN*ds; // Calculate max initialized s for Lin Table
170 G4double sv=0;
171 lastT[0]=1.;
172 for(G4int j=1; j<=lastN; j++) // Calculate Lin Table values
173 {
174 sv+=ds;
175 lastT[j]=CalcQF2IN_Ratio(sv,A);
176 }
177 lastL=new G4double[mls]; // Create the Logarithmic Table
178 if(s_value>sma) // Initialize the Log Table (at least a part of it)
179 {
180#ifdef pdebug
181 G4cout<<"G4QuasiFreeRatios::GetQF2IN_Ratio: NewL, A="<<A<<", nDB="<<nDB<<G4endl;
182#endif
183 G4double ls=std::log(s_value);
184 lastK = static_cast<int>((ls-lsi)/dl)+1; // MaxBin to be initialized in Log Table
185 if(lastK>nls) // The Log Table should be completely initialized
186 {
187 lastK=nls;
188 lastM=lsa-lsi; // lastM is a ln(s)-difference (not a s-difference)
189 }
190 else lastM = lastK*dl; // Calculate max initialized ln(s)-lsi for LogTab
191 sv=mi; // Min s (not ln(s)) for the Log Table
192 for(G4int j=0; j<=lastK; j++) // Calculate Log Table values
193 {
194 lastL[j]=CalcQF2IN_Ratio(sv,A);
195 if(j!=lastK) sv*=edl;
196 }
197 }
198 else // Log Tab is not initialized for the firs call
199 {
200 lastK = 0;
201 lastM = 0.;
202 }
203 i++; // Make a new record to AMDB and position on it
204 vA.push_back(lastA);
205 vH.push_back(lastH);
206 vN.push_back(lastN);
207 vM.push_back(lastM);
208 vK.push_back(lastK);
209 vT.push_back(lastT);
210 vL.push_back(lastL);
211 }
212 else // The A value was found in AMDB
213 {
214 lastA=vA[i];
215 lastH=vH[i];
216 lastN=vN[i];
217 lastM=vM[i];
218 lastK=vK[i];
219 lastT=vT[i];
220 lastL=vL[i];
221#ifdef pdebug
222 G4cout<<"G4QuasiFreeRatios::GetQF2IN_Ratio: Found, s="<<s_value<<", lastM="<<lastM<<G4endl;
223#endif
224 if(s_value>lastH) // At least Lin Table must be updated
225 {
226 G4int nextN=lastN+1; // The next bin to be initialized in Lin Table (p)
227#ifdef pdebug
228 G4cout<<"G4QuasiFreeRatios::GetQF2IN_Ratio: lastN="<<lastN<<" ?< nps="<<nps<<G4endl;
229#endif
230 if(lastN < nps) // The Lin Table is not completely initialized
231 {
232 lastN = static_cast<int>(s_value/ds)+1;// MaxBin to be initialized in Lin Table
233 G4double sv=lastH;
234 if(lastN>nps)
235 {
236 lastN=nps;
237 lastH=sma;
238 }
239 else lastH = lastN*ds; // Calculate max initialized s for LinTab
240 for(G4int j=nextN; j<=lastN; j++)// Calculate LogTab values
241 {
242 sv+=ds;
243 lastT[j]=CalcQF2IN_Ratio(sv,A);
244 }
245#ifdef pdebug
246 G4cout<<"G4QuasiFreeRatios::GetQF2IN_Ratio: End of LinTab update"<<G4endl;
247#endif
248 } // End of LinTab update
249#ifdef pdebug
250 G4cout<<"G4QFRatios::GetQF2IN_Ratio: lN="<<lastN<<", nN="<<nextN<<", i="<<i<<G4endl;
251#endif
252 if(lastN >= nextN) // The Lin Table update really took place
253 {
254 vH[i]=lastH;
255 vN[i]=lastN;
256 }
257 G4int nextK=lastK+1; // Possible next element in LogTable to be updated
258 if(!lastK) nextK=0; // The Log Table has not been filled at all
259#ifdef pdebug
260 G4cout<<"G4QFRat::GetQF2IN_Ratio: sma="<<sma<<", lastK="<<lastK<<" < "<<nls<<G4endl;
261#endif
262 if(s_value>sma && lastK<nls) // LogTab must be updated (not filled completely)
263 {
264 G4double sv=std::exp(lastM+lsi); // Define starting s-poit (lastM will be changed)
265 G4double ls=std::log(s_value);
266 lastK = static_cast<int>((ls-lsi)/dl)+1; // MaxBin to be initialized in Log Table
267 if(lastK>nls)
268 {
269 lastK=nls;
270 lastM=lsa-lsi;
271 }
272 else lastM = lastK*dl; // Calculate max initialized ln(s)-lsi for LogTab
273#ifdef pdebug
274 G4cout<<"G4QFRat::GetQF2IN_Ratio: nK="<<nextK<<", lK="<<lastK<<", sv="<<sv<<G4endl;
275#endif
276 for(G4int j=nextK; j<=lastK; j++)// Calculate LogTab values
277 {
278 sv*=edl;
279#ifdef pdebug
280 G4cout<<"G4QFRat::GetQF2IN_Ratio: j="<<j<<", sv="<<sv<<", A="<<A<<G4endl;
281#endif
282 lastL[j]=CalcQF2IN_Ratio(sv,A);
283 }
284#ifdef pdebug
285 G4cout<<"G4QuasiFreeRatios::GetQF2IN_Ratio: End of LinTab update"<<G4endl;
286#endif
287 } // End of LogTab update
288#ifdef pdebug
289 G4cout<<"G4QFRatios::GetQF2IN_Ratio: lK="<<lastK<<", nK="<<nextK<<", i="<<i<<G4endl;
290#endif
291 if(lastK >= nextK) // The Lin Table update really took place
292 {
293 vM[i]=lastM;
294 vK[i]=lastK;
295 }
296 }
297 }
298#ifdef pdebug
299 G4cout<<"G4QuasiFreeRatios::GetQF2IN_Ratio: BeforeTab s="<<s_value<<", sma="<<sma<<G4endl;
300#endif
301 // Now one can use tabeles to calculate the value
302 if(s_value<sma) // Use linear table
303 {
304 G4int n=static_cast<int>(s_value/ds); // Low edge number of the bin
305 G4double d=s_value-n*ds; // Linear shift
306 G4double v=lastT[n]; // Base
307 lastR=v+d*(lastT[n+1]-v)/ds; // Result
308 }
309 else // Use log table
310 {
311 G4double ls=std::log(s_value)-lsi; // ln(s)-l_min
312 G4int n=static_cast<int>(ls/dl); // Low edge number of the bin
313 G4double d=ls-n*dl; // Log shift
314 G4double v=lastL[n]; // Base
315 lastR=v+d*(lastL[n+1]-v)/dl; // Result
316 }
317 if(lastR<0.) lastR=0.;
318 if(lastR>1.) lastR=1.;
319#ifdef pdebug
320 G4cout<<"G4QuasiFreeRatios::GetQF2IN_Ratio: BeforeRet lastR="<<lastR<<G4endl;
321#endif
322 return lastR;
323} // End of CalcQF2IN_Ratio
324
325// Calculatio QasiFree/Inelastic Ratio as a function of total hN cross-section and A
326G4double G4QuasiFreeRatios::CalcQF2IN_Ratio(G4double s_value, G4int A)
327{
328 static const G4double C=1.246;
329 G4double s2=s_value*s_value;
330 G4double s4=s2*s2;
331 G4double ss=std::sqrt(std::sqrt(s_value));
332 G4double P=7.48e-5*s2/(1.+8.77e12/s4/s4/s2);
333 G4double E=.2644+.016/(1.+std::exp((29.54-s_value)/2.49));
334 G4double F=ss*.1526*std::exp(-s2*ss*.0000859);
335 return C*std::exp(-E*std::pow(G4double(A-1.),F))/std::pow(G4double(A),P);
336} // End of CalcQF2IN_Ratio
337
338// For hadron PDG with momentum Mom (GeV/c) on N (p/n) calculate <sig_el,sig_tot> pair (mb)
339std::pair<G4double,G4double> G4QuasiFreeRatios::GetElTotXS(G4double p, G4int PDG, G4bool F)
340{
341 G4int ind=0; // Prototype of the reaction index
342 G4bool kfl=true; // Flag of K0/aK0 oscillation
343 G4bool kf=false;
344 if(PDG==130||PDG==310)
345 {
346 kf=true;
347 if(G4UniformRand()>.5) kfl=false;
348 }
349 if ( (PDG == 2212 && F) || (PDG == 2112 && !F) ) ind=0; // pp/nn
350 else if ( (PDG == 2112 && F) || (PDG == 2212 && !F) ) ind=1; // np/pn
351 else if ( (PDG == -211 && F) || (PDG == 211 && !F) ) ind=2; // pimp/pipn
352 else if ( (PDG == 211 && F) || (PDG == -211 && !F) ) ind=3; // pipp/pimn
353 else if ( PDG == -321 || PDG == -311 || (kf && !kfl) ) ind=4; // KmN/K0N
354 else if ( PDG == 321 || PDG == 311 || (kf && kfl) ) ind=5; // KpN/aK0N
355 else if ( PDG > 3000 && PDG < 3335) ind=6; // @@ for all hyperons - take Lambda
356 else if ( PDG > -3335 && PDG < -2000) ind=7; // @@ for all anti-baryons (anti-p/anti-n)
357 else {
358 G4cout<<"*Error*G4QuasiFreeRatios::CalcElTotXS: PDG="<<PDG
359 <<", while it is defined only for p,n,hyperons,anti-baryons,pi,K/antiK"<<G4endl;
360 G4Exception("G4QuasiFreeRatio::CalcElTotXS:","22",FatalException,"CHIPScrash");
361 }
362 return QFreeScat->CalcElTot(p,ind);
363}
364
365// (Mean Elastic and Mean Total) Cross-Sections (mb) for PDG+(Z,N) at P=p[GeV/c]
366std::pair<G4double,G4double> G4QuasiFreeRatios::GetElTot(G4double pIU, G4int hPDG,
367 G4int Z, G4int N)
368{
369 G4double pGeV=pIU/gigaelectronvolt;
370#ifdef pdebug
371 G4cout<<"-->G4QuasiFreeR::GetElTot: P="<<pIU<<",pPDG="<<hPDG<<",Z="<<Z<<",N="<<N<<G4endl;
372#endif
373 if(Z<1 && N<1)
374 {
375 G4cout<<"-Warning-G4QuasiFreeRatio::GetElTot:Z="<<Z<<",N="<<N<<", return zero"<<G4endl;
376 return std::make_pair(0.,0.);
377 }
378 std::pair<G4double,G4double> hp=QFreeScat->FetchElTot(pGeV, hPDG, true);
379 std::pair<G4double,G4double> hn=QFreeScat->FetchElTot(pGeV, hPDG, false);
380#ifdef pdebug
381 G4cout<<"-OUT->G4QFRat::GetElTot: hp("<<hp.first<<","<<hp.second<<"), hn("<<hn.first<<","
382 <<hn.second<<")"<<G4endl;
383#endif
384 G4double A=(Z+N)/millibarn; // To make the result in independent units(IU)
385 return std::make_pair((Z*hp.first+N*hn.first)/A,(Z*hp.second+N*hn.second)/A);
386} // End of GetElTot
387
388// (Mean Elastic and Mean Total) Cross-Sections (mb) for PDG+(Z,N) at P=p[GeV/c]
389std::pair<G4double,G4double> G4QuasiFreeRatios::GetChExFactor(G4double pIU, G4int hPDG,
390 G4int Z, G4int N)
391{
392 G4double pGeV=pIU/gigaelectronvolt;
393 G4double resP=0.;
394 G4double resN=0.;
395 if(Z<1 && N<1)
396 {
397 G4cout<<"-Warning-G4QuasiFreeRatio::GetChExF:Z="<<Z<<",N="<<N<<", return zero"<<G4endl;
398 return std::make_pair(resP,resN);
399 }
400 G4double A=Z+N;
401 G4double pf=0.; // Possibility to interact with a proton
402 G4double nf=0.; // Possibility to interact with a neutron
403 if (hPDG==-211||hPDG==-321||hPDG==3112||hPDG==3212||hPDG==3312||hPDG==3334) pf=Z/(A+N);
404 else if(hPDG==211||hPDG==321||hPDG==3222||hPDG==3212||hPDG==3322) nf=N/(A+Z);
405 else if(hPDG==-311||hPDG==311||hPDG==130||hPDG==310)
406 {
407 G4double dA=A+A;
408 pf=Z/(dA+N+N);
409 nf=N/(dA+Z+Z);
410 }
411 G4double mult=1.; // Factor of increasing multiplicity ( ? @@)
412 if(pGeV>.5)
413 {
414 mult=1./(1.+std::log(pGeV+pGeV))/pGeV;
415 if(mult>1.) mult=1.;
416 }
417 if(pf)
418 {
419 std::pair<G4double,G4double> hp=QFreeScat->FetchElTot(pGeV, hPDG, true);
420 resP=pf*(hp.second/hp.first-1.)*mult;
421 }
422 if(nf)
423 {
424 std::pair<G4double,G4double> hn=QFreeScat->FetchElTot(pGeV, hPDG, false);
425 resN=nf*(hn.second/hn.first-1.)*mult;
426 }
427 return std::make_pair(resP,resN);
428} // End of GetChExFactor
429
430// scatter (pPDG,p4M) on a virtual nucleon (NPDG,N4M), result: final pair(newN4M,newp4M)
431// if(newN4M.e()==0.) - below threshold, XS=0, no scattering of the progectile happened
432std::pair<G4LorentzVector,G4LorentzVector> G4QuasiFreeRatios::Scatter(G4int NPDG,
433 G4LorentzVector N4M, G4int pPDG, G4LorentzVector p4M)
434{
435 static const G4double mNeut= G4QPDGCode(2112).GetMass();
436 static const G4double mProt= G4QPDGCode(2212).GetMass();
437 static const G4double mDeut= G4QPDGCode(2112).GetNuclMass(1,1,0);// Mass of deuteron
438 static const G4double mTrit= G4QPDGCode(2112).GetNuclMass(1,2,0);// Mass of tritium
439 static const G4double mHel3= G4QPDGCode(2112).GetNuclMass(2,1,0);// Mass of Helium3
440 static const G4double mAlph= G4QPDGCode(2112).GetNuclMass(2,2,0);// Mass of alpha
441 static const G4double two3d= 2./3.;
442 static const G4double two3d2= std::pow(2.,two3d); // The slope reductions for fragments
443 static const G4double two3d3= std::pow(3.,two3d);
444 static const G4double two3d4= std::pow(4.,two3d);
445 G4LorentzVector pr4M=p4M/megaelectronvolt; // Convert 4-momenta in MeV (keep p4M)
446 N4M/=megaelectronvolt;
447 G4LorentzVector tot4M=N4M+p4M;
448#ifdef ppdebug
449 G4cerr<<"->G4QFR::Scat:p4M="<<pr4M<<",N4M="<<N4M<<",t4M="<<tot4M<<",NPDG="<<NPDG<<G4endl;
450#endif
451 G4double mT=mNeut;
452 G4int Z=0;
453 G4int N=1;
454 if(NPDG==2212||NPDG==90001000)
455 {
456 mT=mProt;
457 Z=1;
458 N=0;
459 }
460 else if(NPDG==90001001)
461 {
462 mT=mDeut;
463 Z=1;
464 N=1;
465 }
466 else if(NPDG==90002001)
467 {
468 mT=mHel3;
469 Z=2;
470 N=1;
471 }
472 else if(NPDG==90001002)
473 {
474 mT=mTrit;
475 Z=1;
476 N=2;
477 }
478 else if(NPDG==90002002)
479 {
480 mT=mAlph;
481 Z=2;
482 N=2;
483 }
484 else if(NPDG!=2112&&NPDG!=90000001)
485 {
486 G4cout<<"Error:G4QuasiFreeRatios::Scatter:NPDG="<<NPDG<<" is not 2212 or 2112"<<G4endl;
487 G4Exception("G4QuasiFreeRatios::Scatter:","21",FatalException,"CHIPScomplain");
488 //return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M);// Use this if not exception
489 }
490 G4double mT2=mT*mT;
491 G4double mP2=pr4M.m2();
492 G4double E=(tot4M.m2()-mT2-mP2)/(mT+mT);
493#ifdef pdebug
494 G4cerr<<"G4QFR::Scat:qM="<<mT<<",qM2="<<mT2<<",pM2="<<mP2<<",totM2="<<tot4M.m2()<<G4endl;
495#endif
496 G4double E2=E*E;
497 if(E<0. || E2<mP2)
498 {
499#ifdef ppdebug
500 G4cerr<<"-Warning-G4QFR::Scat:*Negative Energy*E="<<E<<",E2="<<E2<<"<M2="<<mP2<<G4endl;
501#endif
502 return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); // Do Nothing Action
503 }
504 G4double P=std::sqrt(E2-mP2); // Momentum in pseudo laboratory system
505 G4VQCrossSection* PCSmanager=G4QProtonElasticCrossSection::GetPointer();
506 G4VQCrossSection* NCSmanager=G4QNeutronElasticCrossSection::GetPointer();
507#ifdef ppdebug
508 G4cout<<"G4QFR::Scatter: Before XS, P="<<P<<", Z="<<Z<<", N="<<N<<", PDG="<<pPDG<<G4endl;
509#endif
510 // @@ Temporary NN t-dependence for all hadrons
511 if(pPDG>3400 || pPDG<-3400) G4cout<<"-Warning-G4QuasiFree::Scatter: pPDG="<<pPDG<<G4endl;
512 G4int PDG=2212; // *TMP* instead of pPDG
513 if(pPDG==2112||pPDG==-211||pPDG==-321) PDG=2112; // *TMP* instead of pPDG
514 if(!Z && N==1) // Change for Quasi-Elastic on neutron
515 {
516 Z=1;
517 N=0;
518 if (PDG==2212) PDG=2112;
519 else if(PDG==2112) PDG=2212;
520 }
521 G4double xSec=0.; // Prototype of Recalculated Cross Section *TMP*
522 if(PDG==2212) xSec=PCSmanager->GetCrossSection(false, P, Z, N, PDG); // P CrossSect *TMP*
523 else xSec=NCSmanager->GetCrossSection(false, P, Z, N, PDG); // N CrossSect *TMP*
524#ifdef ppdebug
525 G4cout<<"G4QuasiFreeRat::Scatter: pPDG="<<pPDG<<",P="<<P<<",CS="<<xSec/millibarn<<G4endl;
526#endif
527#ifdef nandebug
528 if(xSec>0. || xSec<0. || xSec==0);
529 else G4cout<<"*Warning*G4QuasiFreeRatios::Scatter: xSec="<<xSec/millibarn<<G4endl;
530#endif
531 // @@ check a possibility to separate p, n, or alpha (!)
532 if(xSec <= 0.) // The cross-section iz 0 -> Do Nothing
533 {
534#ifdef ppdebug
535 G4cerr<<"-Warning-G4QFR::Scat:**Zero XS**PDG="<<pPDG<<",NPDG="<<NPDG<<",P="<<P<<G4endl;
536#endif
537 return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); //Do Nothing Action
538 }
539 G4double mint=0.; // Prototype of functional rand -t (MeV^2) *TMP*
540 if(PDG==2212) mint=PCSmanager->GetExchangeT(Z,N,PDG);// P functional rand -t(MeV^2) *TMP*
541 else mint=NCSmanager->GetExchangeT(Z,N,PDG);// N functional rand -t(MeV^2) *TMP*
542 if (mT==mDeut) mint/=two3d2;
543 else if(mT==mTrit || mT==mHel3) mint/=two3d3;
544 else if(mT==mAlph) mint/=two3d4;
545 G4double maxt=0.; // Prototype of max possible -t
546 if(PDG==2212) maxt=PCSmanager->GetHMaxT(); // max possible -t
547 else maxt=NCSmanager->GetHMaxT(); // max possible -t
548#ifdef ppdebug
549 G4cout<<"G4QFR::Scat:PDG="<<PDG<<",P="<<P<<",X="<<xSec<<",-t="<<mint<<"<"<<maxt<<", Z="
550 <<Z<<",N="<<N<<G4endl;
551#endif
552#ifdef nandebug
553 if(mint>-.0000001);
554 else G4cout<<"*Warning*G4QFR::Scat: -t="<<mint<<G4endl;
555#endif
556 G4double cost=1.-(mint+mint)/maxt; // cos(theta) in CMS
557#ifdef ppdebug
558 G4cout<<"G4QFR::Scat:-t="<<mint<<"<"<<maxt<<", cost="<<cost<<", Z="<<Z<<",N="<<N<<G4endl;
559#endif
560 if(cost>1. || cost<-1. || !(cost>-1. || cost<=1.))
561 {
562 if (cost>1.) cost=1.;
563 else if(cost<-1.) cost=-1.;
564 else
565 {
566 G4double tm=0.;
567 if(PDG==2212) tm=PCSmanager->GetHMaxT();
568 else tm=NCSmanager->GetHMaxT();
569 G4cerr<<"G4QuasiFreeRatio::Scat:*NAN* cost="<<cost<<",-t="<<mint<<",tm="<<tm<<G4endl;
570 return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); // Do Nothing Action
571 }
572 }
573 G4LorentzVector reco4M=G4LorentzVector(0.,0.,0.,mT); // 4mom of the recoil nucleon
574 G4LorentzVector dir4M=tot4M-G4LorentzVector(0.,0.,0.,(tot4M.e()-mT)*.01);
575 if(!G4QHadron(tot4M).RelDecayIn2(pr4M, reco4M, dir4M, cost, cost))
576 {
577 G4cerr<<"G4QFR::Scat:t="<<tot4M<<tot4M.m()<<",mT="<<mT<<",mP="<<std::sqrt(mP2)<<G4endl;
578 //G4Exception("G4QFR::Scat:","009",FatalException,"Decay of ElasticComp");
579 return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); // Do Nothing Action
580 }
581#ifdef ppdebug
582 G4cout<<"G4QFR::Scat:p4M="<<pr4M<<"+r4M="<<reco4M<<",dr="<<dir4M<<",t4M="<<tot4M<<G4endl;
583#endif
584 return std::make_pair(reco4M*megaelectronvolt,pr4M*megaelectronvolt); // Result
585} // End of Scatter
586
587// scatter (pPDG,p4M) on a virtual nucleon (NPDG,N4M), result: final pair(newN4M,newp4M)
588// if(newN4M.e()==0.) - below threshold, XS=0, no scattering of the progectile happened
589// User should himself change the charge (PDG) (e.g. pn->np, pi+n->pi0p, pi-p->pi0n etc.)
590std::pair<G4LorentzVector,G4LorentzVector> G4QuasiFreeRatios::ChExer(G4int NPDG,
591 G4LorentzVector N4M, G4int pPDG, G4LorentzVector p4M)
592{
593 static const G4double mNeut= G4QPDGCode(2112).GetMass();
594 static const G4double mProt= G4QPDGCode(2212).GetMass();
595 G4LorentzVector pr4M=p4M/megaelectronvolt; // Convert 4-momenta in MeV(keep p4M)
596 N4M/=megaelectronvolt;
597 G4LorentzVector tot4M=N4M+p4M;
598 G4int Z=0;
599 G4int N=1;
600 G4int sPDG=0; // PDG code of the scattered hadron
601 G4double mS=0.; // proto of mass of scattered hadron
602 G4double mT=mProt; // mass of the recoil nucleon
603 G4cout<<"-Warning-G4QFR::ChEx: Impossible for Omega-"<<G4endl; // Omega-
604 if(NPDG==2212)
605 {
606 mT=mNeut;
607 Z=1;
608 N=0;
609 if(pPDG==-211) sPDG=111; // pi+ -> pi0
610 else if(pPDG==-321)
611 {
612 sPDG=310; // K+ -> K0S
613 if(G4UniformRand()>.5) sPDG=130; // K+ -> K0L
614 }
615 else if(pPDG==-311||pPDG==311||pPDG==130||pPDG==310) sPDG=321; // K0 -> K+ (?)
616 else if(pPDG==3112) sPDG=3212; // Sigma- -> Sigma0
617 else if(pPDG==3212) sPDG=3222; // Sigma0 -> Sigma+
618 else if(pPDG==3312) sPDG=3322; // Xi- -> Xi0
619 }
620 else if(NPDG==2112) // Default
621 {
622 if(pPDG==211) sPDG=111; // pi+ -> pi0
623 else if(pPDG==321)
624 {
625 sPDG=310; // K+ -> K0S
626 if(G4UniformRand()>.5) sPDG=130; // K+ -> K0L
627 }
628 else if(pPDG==-311||pPDG==311||pPDG==130||pPDG==310) sPDG=-321; // K0 -> K- (?)
629 else if(pPDG==3222) sPDG=3212; // Sigma+ -> Sigma0
630 else if(pPDG==3212) sPDG=3112; // Sigma0 -> Sigma-
631 else if(pPDG==3322) sPDG=3312; // Xi0 -> Xi-
632 }
633 else
634 {
635 G4cout<<"Error:G4QuasiFreeRatios::ChExer: NPDG="<<NPDG<<" is not 2212 or 2112"<<G4endl;
636 G4Exception("G4QuasiFreeRatios::ChExer:","21",FatalException,"CHIPS complain");
637 //return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M);// Use this if not exception
638 }
639 if(sPDG) mS=mNeut;
640 else
641 {
642 G4cout<<"Error:G4QuasiFreeRatios::ChExer: BAD pPDG="<<pPDG<<", NPDG="<<NPDG<<G4endl;
643 G4Exception("G4QuasiFreeRatios::ChExer:","21",FatalException,"CHIPS complain");
644 //return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M);// Use this if not exception
645 }
646 G4double mT2=mT*mT;
647 G4double mS2=mS*mS;
648 G4double E=(tot4M.m2()-mT2-mS2)/(mT+mT);
649 G4double E2=E*E;
650 if(E<0. || E2<mS2)
651 {
652#ifdef pdebug
653 G4cerr<<"-Warning-G4QFR::ChEx:*Negative Energy*E="<<E<<",E2="<<E2<<"<M2="<<mS2<<G4endl;
654#endif
655 return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); // Do Nothing Action
656 }
657 G4double P=std::sqrt(E2-mS2); // Momentum in pseudo laboratory system
658 G4VQCrossSection* PCSmanager=G4QProtonElasticCrossSection::GetPointer();
659 G4VQCrossSection* NCSmanager=G4QNeutronElasticCrossSection::GetPointer();
660#ifdef debug
661 G4cout<<"G4QFR::ChExer: Before XS, P="<<P<<", Z="<<Z<<", N="<<N<<", PDG="<<pPDG<<G4endl;
662#endif
663 // @@ Temporary NN t-dependence for all hadrons
664 G4int PDG=2212; // *TMP* instead of pPDG
665 if(pPDG==2112||pPDG==-211||pPDG==-321) PDG=2112; // *TMP* instead of pPDG
666 if(!Z && N==1) // Change for Quasi-Elastic on neutron
667 {
668 Z=1;
669 N=0;
670 if (PDG==2212) PDG=2112;
671 else if(PDG==2112) PDG=2212;
672 }
673 G4double xSec=0.; // Prototype of Recalculated Cross Section *TMP*
674 if(PDG==2212) xSec=PCSmanager->GetCrossSection(false, P, Z, N, PDG); // P CrossSect *TMP*
675 else xSec=NCSmanager->GetCrossSection(false, P, Z, N, PDG); // N CrossSect *TMP*
676#ifdef debug
677 G4cout<<"G4QuasiFreeRat::ChExer: pPDG="<<pPDG<<",P="<<P<<", CS="<<xSec/millibarn<<G4endl;
678#endif
679#ifdef nandebug
680 if(xSec>0. || xSec<0. || xSec==0);
681 else G4cout<<"*Warning*G4QuasiFreeRatios::ChExer: xSec="<<xSec/millibarn<<G4endl;
682#endif
683 // @@ check a possibility to separate p, n, or alpha (!)
684 if(xSec <= 0.) // The cross-section iz 0 -> Do Nothing
685 {
686#ifdef pdebug
687 G4cerr<<"-Warning-G4QFR::ChEx:**Zero XS**PDG="<<pPDG<<",NPDG="<<NPDG<<",P="<<P<<G4endl;
688#endif
689 return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); //Do Nothing Action
690 }
691 G4double mint=0.; // Prototype of functional rand -t (MeV^2) *TMP*
692 if(PDG==2212) mint=PCSmanager->GetExchangeT(Z,N,PDG);// P functional rand -t(MeV^2) *TMP*
693 else mint=NCSmanager->GetExchangeT(Z,N,PDG);// N functional rand -t(MeV^2) *TMP*
694#ifdef pdebug
695 G4cout<<"G4QFR::ChEx:PDG="<<pPDG<<", P="<<P<<", CS="<<xSec<<", -t="<<mint<<G4endl;
696#endif
697#ifdef nandebug
698 if(mint>-.0000001);
699 else G4cout<<"*Warning*G4QFR::ChExer: -t="<<mint<<G4endl;
700#endif
701 G4double maxt=0.; // Prototype of max possible -t
702 if(PDG==2212) maxt=PCSmanager->GetHMaxT(); // max possible -t
703 else maxt=NCSmanager->GetHMaxT(); // max possible -t
704 G4double cost=1.-mint/maxt; // cos(theta) in CMS
705#ifdef pdebug
706 G4cout<<"G4QuasiFfreeRatio::ChExer: -t="<<mint<<", maxt="<<maxt<<", cost="<<cost<<G4endl;
707#endif
708 if(cost>1. || cost<-1. || !(cost>-1. || cost<=1.))
709 {
710 if (cost>1.) cost=1.;
711 else if(cost<-1.) cost=-1.;
712 else
713 {
714 G4cerr<<"G4QuasiFreeRatio::ChExer:*NAN* c="<<cost<<",t="<<mint<<",tm="<<maxt<<G4endl;
715 return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); // Do Nothing Action
716 }
717 }
718 G4LorentzVector reco4M=G4LorentzVector(0.,0.,0.,mT); // 4mom of the recoil nucleon
719 pr4M=G4LorentzVector(0.,0.,0.,mS); // 4mom of the scattered hadron
720 G4LorentzVector dir4M=tot4M-G4LorentzVector(0.,0.,0.,(tot4M.e()-mT)*.01);
721 if(!G4QHadron(tot4M).RelDecayIn2(pr4M, reco4M, dir4M, cost, cost))
722 {
723 G4cerr<<"G4QFR::ChEx:t="<<tot4M<<tot4M.m()<<",mT="<<mT<<",mP="<<mS<<G4endl;
724 //G4Exception("G4QFR::ChExer:","009",FatalException,"Decay of ElasticComp");
725 return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); // Do Nothing Action
726 }
727#ifdef debug
728 G4cout<<"G4QFR::ChEx:p4M="<<p4M<<"+r4M="<<reco4M<<"="<<p4M+reco4M<<"="<<tot4M<<G4endl;
729#endif
730 return std::make_pair(reco4M*megaelectronvolt,pr4M*megaelectronvolt); // Result
731} // End of ChExer
732
733// Calculate ChEx/El ratio (p is in independent units, (Z,N) is target, pPDG is projectile)
734G4double G4QuasiFreeRatios::ChExElCoef(G4double p, G4int Z, G4int N, G4int pPDG)
735{
736 p/=MeV; // Converted from independent units
737 G4double A=Z+N;
738 if(A<1.5) return 0.;
739 G4double C=0.;
740 if (pPDG==2212) C=N/(A+Z);
741 else if(pPDG==2112) C=Z/(A+N);
742 else G4cout<<"*Warning*G4QCohChrgExchange::ChExElCoef: wrong PDG="<<pPDG<<G4endl;
743 C*=C; // Coherent processes squares the amplitude
744 // @@ This is true only for nucleons: other projectiles must be treated differently
745 G4double sp=std::sqrt(p);
746 G4double p2=p*p;
747 G4double p4=p2*p2;
748 G4double dl1=std::log(p)-5.;
749 G4double T=(6.75+.14*dl1*dl1+13./p)/(1.+.14/p4)+.6/(p4+.00013);
750 G4double U=(6.25+8.33e-5/p4/p)*(p*sp+.34)/p2/p;
751 G4double R=U/T;
752 return C*R*R;
753}
FatalException
@ FatalException
Definition: G4ExceptionSeverity.hh:60
G4LorentzVector
CLHEP::HepLorentzVector G4LorentzVector
Definition: G4LorentzVector.hh:35
G4double
double G4double
Definition: G4Types.hh:64
G4int
int G4int
Definition: G4Types.hh:66
G4bool
bool G4bool
Definition: G4Types.hh:67
G4endl
#define G4endl
Definition: G4ios.hh:52
G4cerr
G4DLLIMPORT std::ostream G4cerr
G4cout
G4DLLIMPORT std::ostream G4cout
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:53
G4QFreeScattering::FetchElTot
std::pair< G4double, G4double > FetchElTot(G4double pGeV, G4int PDG, G4bool F)
Definition: G4QFreeScattering.cc:373
G4QFreeScattering::CalcElTot
std::pair< G4double, G4double > CalcElTot(G4double pGeV, G4int Index)
Definition: G4QFreeScattering.cc:71
G4QFreeScattering::GetPointer
static G4QFreeScattering * GetPointer()
Definition: G4QFreeScattering.cc:64
G4QPDGCode::GetMass
G4double GetMass()
Definition: G4QPDGCode.cc:693
G4QPDGCode::GetNuclMass
G4double GetNuclMass(G4int Z, G4int N, G4int S)
Definition: G4QPDGCode.cc:766
G4QProtonElasticCrossSection::GetPointer
static G4VQCrossSection * GetPointer()
Definition: G4QProtonElasticCrossSection.cc:143
G4QuasiFreeRatios::ChExer
std::pair< G4LorentzVector, G4LorentzVector > ChExer(G4int NPDG, G4LorentzVector N4M, G4int pPDG, G4LorentzVector p4M)
Definition: G4QuasiFreeRatios.cc:590
G4QuasiFreeRatios::GetChExFactor
std::pair< G4double, G4double > GetChExFactor(G4double pIU, G4int pPDG, G4int Z, G4int N)
Definition: G4QuasiFreeRatios.cc:389
G4QuasiFreeRatios::ChExElCoef
G4double ChExElCoef(G4double p, G4int Z, G4int N, G4int pPDG)
Definition: G4QuasiFreeRatios.cc:734
G4QuasiFreeRatios::Scatter
std::pair< G4LorentzVector, G4LorentzVector > Scatter(G4int NPDG, G4LorentzVector N4M, G4int pPDG, G4LorentzVector p4M)
Definition: G4QuasiFreeRatios.cc:432
G4QuasiFreeRatios::GetElTot
std::pair< G4double, G4double > GetElTot(G4double pIU, G4int hPDG, G4int Z, G4int N)
Definition: G4QuasiFreeRatios.cc:366
G4QuasiFreeRatios::GetPointer
static G4QuasiFreeRatios * GetPointer()
Definition: G4QuasiFreeRatios.cc:69
G4QuasiFreeRatios::GetElTotXS
std::pair< G4double, G4double > GetElTotXS(G4double Mom, G4int PDG, G4bool F)
Definition: G4QuasiFreeRatios.cc:339
G4QuasiFreeRatios::GetRatios
std::pair< G4double, G4double > GetRatios(G4double pIU, G4int prPDG, G4int tgZ, G4int tgN)
Definition: G4QuasiFreeRatios.cc:76
G4VQCrossSection::GetExchangeT
virtual G4double GetExchangeT(G4int tZ, G4int tN, G4int pPDG)
Definition: G4VQCrossSection.cc:70
G4VQCrossSection::GetCrossSection
virtual G4double GetCrossSection(G4bool, G4double, G4int, G4int, G4int pPDG=0)
Definition: G4VQCrossSection.hh:99
G4VQCrossSection::GetHMaxT
virtual G4double GetHMaxT()
Definition: G4VQCrossSection.cc:72
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41