Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
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