Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4QuasiElRatios.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: G4QuasiElRatios.cc 96538 2016-04-20 12:39:18Z gcosmo $
28 //
29 //
30 // G4 Physics class: G4QuasiElRatios 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) 15-Oct-06
33 //
34 // ----------------------------------------------------------------------
35 // This class has been extracted from the CHIPS model.
36 // All the dependencies on CHIPS classes have been removed.
37 // Short description: Provides percentage of quasi-free and quasi-elastic
38 // reactions in the inelastic reactions.
39 // ----------------------------------------------------------------------
40 
41 
42 #include "G4QuasiElRatios.hh"
43 #include "G4PhysicalConstants.hh"
44 #include "G4SystemOfUnits.hh"
45 #include "G4Proton.hh"
46 #include "G4Neutron.hh"
47 #include "G4Deuteron.hh"
48 #include "G4Triton.hh"
49 #include "G4He3.hh"
50 #include "G4Alpha.hh"
51 #include "G4ThreeVector.hh"
53 #include "G4Pow.hh"
54 #include "G4Log.hh"
55 #include "G4Exp.hh"
56 
57 namespace {
58  const G4int nps=150; // Number of steps in the R(s) LinTable
59  const G4int mps=nps+1; // Number of elements in the R(s) LinTable
60  const G4double sma=150.; // The first LinTabEl(s=0)=1., s>sma -> logTab
61  const G4double ds=sma/nps; // Step of the linear Table
62  const G4int nls=100; // Number of steps in the R(lns) logTable
63  const G4int mls=nls+1; // Number of elements in the R(lns) logTable
64  const G4double lsi=5.; // The min ln(s) logTabEl(s=148.4 < sma=150.)
65  const G4double lsa=9.; // The max ln(s) logTabEl(s=148.4 - 8103. mb)
66  const G4double mi=G4Exp(lsi);// The min s of logTabEl(~ 148.4 mb)
67  const G4double min_s=G4Exp(lsa);// The max s of logTabEl(~ 8103. mb)
68  const G4double dls=(lsa-lsi)/nls;// Step of the logarithmic Table
69  const G4double edls=G4Exp(dls);// Multiplication step of the logarithmic Table
70  const G4double toler=.01; // The tolarence mb defining the same cross-section
71  const G4double C=1.246;
72  const G4double lmi=3.5; // min of (lnP-lmi)^2 parabola
73  const G4double pbe=.0557; // elastic (lnP-lmi)^2 parabola coefficient
74  const G4double pbt=.3; // total (lnP-lmi)^2 parabola coefficient
75  const G4double pmi=.1; // Below that fast LE calculation is made
76  const G4double pma=1000.; // Above that fast HE calculation is made
77  const G4int nlp=300; // Number of steps in the S(lnp) logTable(5% step)
78  const G4int mlp=nlp+1; // Number of elements in the S(lnp) logTable
79  const G4double lpi=-5.; // The min ln(p) logTabEl(p=6.7 MeV/c - 22. TeV/c)
80  const G4double lpa=10.; // The max ln(p) logTabEl(p=6.7 MeV/c - 22. TeV/c)
81  const G4double mip=G4Exp(lpi);// The min p of logTabEl(~ 6.7 MeV/c)
82  const G4double map=G4Exp(lpa);// The max p of logTabEl(~ 22. TeV)
83  const G4double dlp=(lpa-lpi)/nlp;// Step of the logarithmic Table
84  const G4double edlp=G4Exp(dlp);// Multiplication step of the logarithmic Table
85 }
86 
87 
89 {
90  vT = new std::vector<G4double*>;
91  vL = new std::vector<G4double*>;
92  vX = new std::vector<std::pair<G4double,G4double>*>;
93 
94  lastSRatio=0.; // The last sigma value for which R was calculated
95  lastRRatio=0.; // The last ratio R which was calculated
96  lastARatio=0; // theLast of calculated A
97  lastHRatio=0.; // theLast of max s initialized in the LinTable
98  lastNRatio=0; // theLast of topBin number initialized in LinTable
99  lastMRatio=0.; // theLast of rel max ln(s) initialized in LogTable
100  lastKRatio=0; // theLast of topBin number initialized in LogTable
101  lastTRatio=0; // theLast of pointer to LinTable in the C++ heap
102  lastLRatio=0; // theLast of pointer to LogTable in the C++ heap
103  lastPtot=0.; // The last momentum for which XS was calculated
104  lastHtot=0; // The last projPDG for which XS was calculated
105  lastFtot=true; // The last nucleon for which XS was calculated
106  lastItot=0; // The Last index for which XS was calculated
107  lastMtot=0.; // The Last rel max ln(p) initialized in LogTable
108  lastKtot=0; // The Last topBin number initialized in LogTable
109  lastXtot = new std::pair<G4double,G4double>; // The Last ETPointers to LogTable in heap
110 
112 
114 }
115 
117 {
118  std::vector<G4double*>::iterator pos;
119  for(pos=vT->begin(); pos<vT->end(); pos++)
120  { delete [] *pos; }
121  vT->clear();
122  delete vT;
123 
124  for(pos=vL->begin(); pos<vL->end(); pos++)
125  { delete [] *pos; }
126  vL->clear();
127  delete vL;
128 
129  std::vector<std::pair<G4double,G4double>*>::iterator pos2;
130  for(pos2=vX->begin(); pos2<vX->end(); pos2++)
131  { delete [] *pos2; }
132  vX->clear();
133  delete vX;
134 
135 }
136 
137 // Calculation of pair(QuasiFree/Inelastic,QuasiElastic/QuasiFree)
138 std::pair<G4double,G4double> G4QuasiElRatios::GetRatios(G4double pIU, G4int pPDG,
139  G4int tgZ, G4int tgN)
140 {
141  G4double R=0.;
142  G4double QF2In=1.; // Prototype of QuasiFree/Inel ratio for hN_tot
143  G4int tgA=tgZ+tgN;
144  if(tgA<2) return std::make_pair(QF2In,R); // No quasi-elastic on the only nucleon
145  std::pair<G4double,G4double> ElTot=GetElTot(pIU, pPDG, tgZ, tgN); // mean hN El&Tot(IU)
146  //if( ( (pPDG>999 && pIU<227.) || pIU<27.) && tgA>1) R=1.; // @@ TMP to accelerate @lowE
147  if(pPDG>999 && pIU<227. && tgZ+tgN>1) R=1.; // To accelerate @lowE
148  else if(ElTot.second>0.)
149  {
150  R=ElTot.first/ElTot.second; // El/Total ratio (does not depend on units
151  QF2In=GetQF2IN_Ratio(ElTot.second/millibarn, tgZ+tgN); // QuasiFree/Inelastic ratio
152  }
153  return std::make_pair(QF2In,R);
154 }
155 
156 // Calculatio QasiFree/Inelastic Ratio as a function of total hN cross-section (mb) and A
157 G4double G4QuasiElRatios::GetQF2IN_Ratio(G4double m_s, G4int A)
158 {
159  // LogTable is created only if necessary. The ratio R(s>8100 mb) = 0 for any nuclei
160  if(m_s<toler || A<2) return 1.;
161  if(m_s>min_s) return 0.;
162 
163  if(A>238)
164  {
165  G4cout<<"-Warning-G4QuasiElRatio::GetQF2IN_Ratio:A="<<A<<">238, return zero"<<G4endl;
166  return 0.;
167  }
168  G4int nDB=vARatio.size(); // A number of nuclei already initialized in AMDB
169  if(nDB && lastARatio==A && m_s==lastSRatio) return lastRRatio; // VI do not use tolerance
170  G4bool found=false;
171  G4int i=-1;
172  if(nDB) for (i=0; i<nDB; i++) if(A==vARatio[i]) // Search for this A in AMDB
173  {
174  found=true; // The A value is found
175  break;
176  }
177  if(!nDB || !found) // Create new line in the AMDB
178  {
179  lastARatio = A;
180  lastTRatio = new G4double[mps]; // Create the linear Table
181  lastNRatio = static_cast<int>(m_s/ds)+1; // MaxBin to be initialized
182  if(lastNRatio>nps)
183  {
184  lastNRatio=nps;
185  lastHRatio=sma;
186  }
187  else lastHRatio = lastNRatio*ds; // Calculate max initialized s for LinTab
188  G4double sv=0;
189  lastTRatio[0]=1.;
190  for(G4int j=1; j<=lastNRatio; j++) // Calculate LogTab values
191  {
192  sv+=ds;
193  lastTRatio[j]=CalcQF2IN_Ratio(sv,A);
194  }
195  lastLRatio=new G4double[mls]; // Create the logarithmic Table
196  if(m_s>sma) // Initialize the logarithmic Table
197  {
198  G4double ls=G4Log(m_s);
199  lastKRatio = static_cast<int>((ls-lsi)/dls)+1; // MaxBin to be initialized in LogTaB
200  if(lastKRatio>nls)
201  {
202  lastKRatio=nls;
203  lastMRatio=lsa-lsi;
204  }
205  else lastMRatio = lastKRatio*dls; // Calculate max initialized ln(s)-lsi for LogTab
206  sv=mi;
207  for(G4int j=0; j<=lastKRatio; j++) // Calculate LogTab values
208  {
209  lastLRatio[j]=CalcQF2IN_Ratio(sv,A);
210  if(j!=lastKRatio) sv*=edls;
211  }
212  }
213  else // LogTab is not initialized
214  {
215  lastKRatio = 0;
216  lastMRatio = 0.;
217  }
218  i++; // Make a new record to AMDB and position on it
219  vARatio.push_back(lastARatio);
220  vHRatio.push_back(lastHRatio);
221  vNRatio.push_back(lastNRatio);
222  vMRatio.push_back(lastMRatio);
223  vKRatio.push_back(lastKRatio);
224  vT->push_back(lastTRatio);
225  vL->push_back(lastLRatio);
226  }
227  else // The A value was found in AMDB
228  {
229  lastARatio=vARatio[i];
230  lastHRatio=vHRatio[i];
231  lastNRatio=vNRatio[i];
232  lastMRatio=vMRatio[i];
233  lastKRatio=vKRatio[i];
234  lastTRatio=(*vT)[i];
235  lastLRatio=(*vL)[i];
236  if(m_s>lastHRatio) // At least LinTab must be updated
237  {
238  G4int nextN=lastNRatio+1; // The next bin to be initialized
239  if(lastNRatio<nps)
240  {
241  G4double sv=lastHRatio; // bug fix by WP
242 
243  lastNRatio = static_cast<int>(m_s/ds)+1;// MaxBin to be initialized
244  if(lastNRatio>nps)
245  {
246  lastNRatio=nps;
247  lastHRatio=sma;
248  }
249  else lastHRatio = lastNRatio*ds; // Calculate max initialized s for LinTab
250 
251  for(G4int j=nextN; j<=lastNRatio; j++)// Calculate LogTab values
252  {
253  sv+=ds;
254  lastTRatio[j]=CalcQF2IN_Ratio(sv,A);
255  }
256  } // End of LinTab update
257  if(lastNRatio>=nextN)
258  {
259  vHRatio[i]=lastHRatio;
260  vNRatio[i]=lastNRatio;
261  }
262  G4int nextK=lastKRatio+1;
263  if(!lastKRatio) nextK=0;
264  if(m_s>sma && lastKRatio<nls) // LogTab must be updated
265  {
266  G4double sv=G4Exp(lastMRatio+lsi); // Define starting poit (lastM will be changed)
267  G4double ls=G4Log(m_s);
268  lastKRatio = static_cast<int>((ls-lsi)/dls)+1; // MaxBin to be initialized in LogTaB
269  if(lastKRatio>nls)
270  {
271  lastKRatio=nls;
272  lastMRatio=lsa-lsi;
273  }
274  else lastMRatio = lastKRatio*dls; // Calculate max initialized ln(s)-lsi for LogTab
275  for(G4int j=nextK; j<=lastKRatio; j++)// Calculate LogTab values
276  {
277  sv*=edls;
278  lastLRatio[j]=CalcQF2IN_Ratio(sv,A);
279  }
280  } // End of LogTab update
281  if(lastKRatio>=nextK)
282  {
283  vMRatio[i]=lastMRatio;
284  vKRatio[i]=lastKRatio;
285  }
286  }
287  }
288  // Now one can use tabeles to calculate the value
289  if(m_s<sma) // Use linear table
290  {
291  G4int n=static_cast<int>(m_s/ds); // Low edge number of the bin
292  G4double d=m_s-n*ds; // Linear shift
293  G4double v=lastTRatio[n]; // Base
294  lastRRatio=v+d*(lastTRatio[n+1]-v)/ds; // Result
295  }
296  else // Use log table
297  {
298  G4double ls=G4Log(m_s)-lsi; // ln(s)-l_min
299  G4int n=static_cast<int>(ls/dls); // Low edge number of the bin
300  G4double d=ls-n*dls; // Log shift
301  G4double v=lastLRatio[n]; // Base
302  lastRRatio=v+d*(lastLRatio[n+1]-v)/dls; // Result
303  }
304  if(lastRRatio<0.) lastRRatio=0.;
305  if(lastRRatio>1.) lastRRatio=1.;
306  return lastRRatio;
307 } // End of CalcQF2IN_Ratio
308 
309 // Calculatio QasiFree/Inelastic Ratio as a function of total hN cross-section and A
310 G4double G4QuasiElRatios::CalcQF2IN_Ratio(G4double m_s, G4int A)
311 {
312  G4double s2=m_s*m_s;
313  G4double s4=s2*s2;
314  G4double ss=std::sqrt(std::sqrt(m_s));
315  G4double P=7.48e-5*s2/(1.+8.77e12/s4/s4/s2);
316  G4double E=.2644+.016/(1.+G4Exp((29.54-m_s)/2.49));
317  G4double F=ss*.1526*G4Exp(-s2*ss*.0000859);
318  return C*G4Exp(-E*G4Pow::GetInstance()->powA(G4double(A-1.),F))/G4Pow::GetInstance()->powA(G4double(A),P);
319 } // End of CalcQF2IN_Ratio
320 
321 // Calculatio pair(hN_el,hN_tot) (mb): p in GeV/c, index(PDG,F) (see FetchElTot)
322 std::pair<G4double,G4double> G4QuasiElRatios::CalcElTot(G4double p, G4int I)
323 {
324  // ---------> Each parameter set can have not more than nPoints=128 parameters
325 
326  G4double El=0.; // prototype of the elastic hN cross-section
327  G4double To=0.; // prototype of the total hN cross-section
328  if(p<=0.)
329  {
330  G4cout<<"-Warning-G4QuasiElRatios::CalcElTot: p="<<p<<" is zero or negative"<<G4endl;
331  return std::make_pair(El,To);
332  }
333  if (!I) // pp/nn
334  {
335  if(p<pmi)
336  {
337  G4double p2=p*p;
338  El=1./(.00012+p2*.2);
339  To=El;
340  }
341  else if(p>pma)
342  {
343  G4double lp=G4Log(p)-lmi;
344  G4double lp2=lp*lp;
345  El=pbe*lp2+6.72;
346  To=pbt*lp2+38.2;
347  }
348  else
349  {
350  G4double p2=p*p;
351  G4double LE=1./(.00012+p2*.2);
352  G4double lp=G4Log(p)-lmi;
353  G4double lp2=lp*lp;
354  G4double rp2=1./p2;
355  El=LE+(pbe*lp2+6.72+32.6/p)/(1.+rp2/p);
356  To=LE+(pbt*lp2+38.2+52.7*rp2)/(1.+2.72*rp2*rp2);
357  }
358  }
359  else if(I==1) // np/pn
360  {
361  if(p<pmi)
362  {
363  G4double p2=p*p;
364  El=1./(.00012+p2*(.051+.1*p2));
365  To=El;
366  }
367  else if(p>pma)
368  {
369  G4double lp=G4Log(p)-lmi;
370  G4double lp2=lp*lp;
371  El=pbe*lp2+6.72;
372  To=pbt*lp2+38.2;
373  }
374  else
375  {
376  G4double p2=p*p;
377  G4double LE=1./(.00012+p2*(.051+.1*p2));
378  G4double lp=G4Log(p)-lmi;
379  G4double lp2=lp*lp;
380  G4double rp2=1./p2;
381  El=LE+(pbe*lp2+6.72+30./p)/(1.+.49*rp2/p);
382  To=LE+(pbt*lp2+38.2)/(1.+.54*rp2*rp2);
383  }
384  }
385  else if(I==2) // pimp/pipn
386  {
387  G4double lp=G4Log(p);
388  if(p<pmi)
389  {
390  G4double lr=lp+1.27;
391  El=1.53/(lr*lr+.0676);
392  To=El*3;
393  }
394  else if(p>pma)
395  {
396  G4double ld=lp-lmi;
397  G4double ld2=ld*ld;
398  G4double sp=std::sqrt(p);
399  El=pbe*ld2+2.4+7./sp;
400  To=pbt*ld2+22.3+12./sp;
401  }
402  else
403  {
404  G4double lr=lp+1.27; // p1
405  G4double LE=1.53/(lr*lr+.0676); // p2, p3
406  G4double ld=lp-lmi; // p4 (lmi=3.5)
407  G4double ld2=ld*ld;
408  G4double p2=p*p;
409  G4double p4=p2*p2;
410  G4double sp=std::sqrt(p);
411  G4double lm=lp+.36; // p5
412  G4double md=lm*lm+.04; // p6
413  G4double lh=lp-.017; // p7
414  G4double hd=lh*lh+.0025; // p8
415  El=LE+(pbe*ld2+2.4+7./sp)/(1.+.7/p4)+.6/md+.05/hd;//p9(pbe=.0557),p10,p11,p12,p13,p14
416  To=LE*3+(pbt*ld2+22.3+12./sp)/(1.+.4/p4)+1./md+.06/hd;
417  }
418  }
419  else if(I==3) // pipp/pimn
420  {
421  G4double lp=G4Log(p);
422  if(p<pmi)
423  {
424  G4double lr=lp+1.27;
425  G4double lr2=lr*lr;
426  El=13./(lr2+lr2*lr2+.0676);
427  To=El;
428  }
429  else if(p>pma)
430  {
431  G4double ld=lp-lmi;
432  G4double ld2=ld*ld;
433  G4double sp=std::sqrt(p);
434  El=pbe*ld2+2.4+6./sp;
435  To=pbt*ld2+22.3+5./sp;
436  }
437  else
438  {
439  G4double lr=lp+1.27; // p1
440  G4double lr2=lr*lr;
441  G4double LE=13./(lr2+lr2*lr2+.0676); // p2, p3
442  G4double ld=lp-lmi; // p4 (lmi=3.5)
443  G4double ld2=ld*ld;
444  G4double p2=p*p;
445  G4double p4=p2*p2;
446  G4double sp=std::sqrt(p);
447  G4double lm=lp-.32; // p5
448  G4double md=lm*lm+.0576; // p6
449  El=LE+(pbe*ld2+2.4+6./sp)/(1.+3./p4)+.7/md; // p7(pbe=.0557), p8, p9, p10, p11
450  To=LE+(pbt*ld2+22.3+5./sp)/(1.+1./p4)+.8/md;
451  }
452  }
453  else if(I==4) // Kmp/Kmn/K0p/K0n
454  {
455 
456  if(p<pmi)
457  {
458  G4double psp=p*std::sqrt(p);
459  El=5.2/psp;
460  To=14./psp;
461  }
462  else if(p>pma)
463  {
464  G4double ld=G4Log(p)-lmi;
465  G4double ld2=ld*ld;
466  El=pbe*ld2+2.23;
467  To=pbt*ld2+19.5;
468  }
469  else
470  {
471  G4double ld=G4Log(p)-lmi;
472  G4double ld2=ld*ld;
473  G4double sp=std::sqrt(p);
474  G4double psp=p*sp;
475  G4double p2=p*p;
476  G4double p4=p2*p2;
477  G4double lm=p-.39;
478  G4double md=lm*lm+.000156;
479  G4double lh=p-1.;
480  G4double hd=lh*lh+.0156;
481  El=5.2/psp+(pbe*ld2+2.23)/(1.-.7/sp+.075/p4)+.004/md+.15/hd;
482  To=14./psp+(pbt*ld2+19.5)/(1.-.21/sp+.52/p4)+.006/md+.30/hd;
483  }
484  }
485  else if(I==5) // Kpp/Kpn/aKp/aKn
486  {
487  if(p<pmi)
488  {
489  G4double lr=p-.38;
490  G4double lm=p-1.;
491  G4double md=lm*lm+.372;
492  El=.7/(lr*lr+.0676)+2./md;
493  To=El+.6/md;
494  }
495  else if(p>pma)
496  {
497  G4double ld=G4Log(p)-lmi;
498  G4double ld2=ld*ld;
499  El=pbe*ld2+2.23;
500  To=pbt*ld2+19.5;
501  }
502  else
503  {
504  G4double ld=G4Log(p)-lmi;
505  G4double ld2=ld*ld;
506  G4double lr=p-.38;
507  G4double LE=.7/(lr*lr+.0676);
508  G4double sp=std::sqrt(p);
509  G4double p2=p*p;
510  G4double p4=p2*p2;
511  G4double lm=p-1.;
512  G4double md=lm*lm+.372;
513  El=LE+(pbe*ld2+2.23)/(1.-.7/sp+.1/p4)+2./md;
514  To=LE+(pbt*ld2+19.5)/(1.+.46/sp+1.6/p4)+2.6/md;
515  }
516  }
517  else if(I==6) // hyperon-N
518  {
519  if(p<pmi)
520  {
521  G4double p2=p*p;
522  El=1./(.002+p2*(.12+p2));
523  To=El;
524  }
525  else if(p>pma)
526  {
527  G4double lp=G4Log(p)-lmi;
528  G4double lp2=lp*lp;
529  G4double sp=std::sqrt(p);
530  El=(pbe*lp2+6.72)/(1.+2./sp);
531  To=(pbt*lp2+38.2+900./sp)/(1.+27./sp);
532  }
533  else
534  {
535  G4double p2=p*p;
536  G4double LE=1./(.002+p2*(.12+p2));
537  G4double lp=G4Log(p)-lmi;
538  G4double lp2=lp*lp;
539  G4double p4=p2*p2;
540  G4double sp=std::sqrt(p);
541  El=LE+(pbe*lp2+6.72+99./p2)/(1.+2./sp+2./p4);
542  To=LE+(pbt*lp2+38.2+900./sp)/(1.+27./sp+3./p4);
543  }
544  }
545  else if(I==7) // antibaryon-N
546  {
547  if(p>pma)
548  {
549  G4double lp=G4Log(p)-lmi;
550  G4double lp2=lp*lp;
551  El=pbe*lp2+6.72;
552  To=pbt*lp2+38.2;
553  }
554  else
555  {
556  G4double ye=G4Pow::GetInstance()->powA(p,1.25);
557  G4double yt=G4Pow::GetInstance()->powA(p,.35);
558  G4double lp=G4Log(p)-lmi;
559  G4double lp2=lp*lp;
560  El=80./(ye+1.)+pbe*lp2+6.72;
561  To=(80./yt+.3)/yt+pbt*lp2+38.2;
562  }
563  }
564  else
565  {
566  G4cout<<"*Error*G4QuasiElRatios::CalcElTot:ind="<<I<<" is not defined (0-7)"<<G4endl;
567  G4Exception("G4QuasiElRatios::CalcElTot:","23",FatalException,"QEcrash");
568  }
569  if(El>To) El=To;
570  return std::make_pair(El,To);
571 } // End of CalcElTot
572 
573 // For hadron PDG with momentum Mom (GeV/c) on N (p/n) calculate <sig_el,sig_tot> pair (mb)
574 std::pair<G4double,G4double> G4QuasiElRatios::GetElTotXS(G4double p, G4int PDG, G4bool F)
575 {
576  G4int ind=0; // Prototype of the reaction index
577  G4bool kfl=true; // Flag of K0/aK0 oscillation
578  G4bool kf=false;
579  if(PDG==130||PDG==310)
580  {
581  kf=true;
582  if(G4UniformRand()>.5) kfl=false;
583  }
584  if ( (PDG == 2212 && F) || (PDG == 2112 && !F) ) ind=0; // pp/nn
585  else if ( (PDG == 2112 && F) || (PDG == 2212 && !F) ) ind=1; // np/pn
586  else if ( (PDG == -211 && F) || (PDG == 211 && !F) ) ind=2; // pimp/pipn
587  else if ( (PDG == 211 && F) || (PDG == -211 && !F) ) ind=3; // pipp/pimn
588  else if ( PDG == -321 || PDG == -311 || (kf && !kfl) ) ind=4; // KmN/K0N
589  else if ( PDG == 321 || PDG == 311 || (kf && kfl) ) ind=5; // KpN/aK0N
590  else if ( PDG > 3000 && PDG < 3335) ind=6; // @@ for all hyperons - take Lambda
591  else if ( PDG > -3335 && PDG < -2000) ind=7; // @@ for all anti-baryons (anti-p/anti-n)
592  else {
593  G4cout<<"*Error*G4QuasiElRatios::CalcElTotXS: PDG="<<PDG
594  <<", while it is defined only for p,n,hyperons,anti-baryons,pi,K/antiK"<<G4endl;
595  G4Exception("G4QuasiElRatio::CalcElTotXS:","22",FatalException,"QEcrash");
596  }
597  return CalcElTot(p,ind);
598 }
599 
600 // Calculatio pair(hN_el,hN_tot)(mb): p in GeV/c, F=true -> N=proton, F=false -> N=neutron
601 std::pair<G4double,G4double> G4QuasiElRatios::FetchElTot(G4double p, G4int PDG, G4bool F)
602 {
603  // LogTable is created only if necessary. The ratio R(s>8100 mb) = 0 for any nuclei
604  G4int nDB=vItot.size(); // A number of hadrons already initialized in AMDB
605  if(nDB && lastHtot==PDG && lastFtot==F && p>0. && p==lastPtot) return lastRtot;// VI don't use toler.
606  // if(nDB && lastH==PDG && lastF==F && p>0. && std::fabs(p-lastP)/p<toler) return lastR;
607  lastHtot=PDG;
608  lastFtot=F;
609  G4int ind=-1; // Prototipe of the index of the PDG/F combination
610  // i=0: pp(nn), i=1: np(pn), i=2: pimp(pipn), i=3: pipp(pimn), i=4: Kmp(Kmn,K0n,K0p),
611  // i=5: Kpp(Kpn,aK0n,aK0p), i=6: Hp(Hn), i=7: app(apn,ann,anp)
612  G4bool kfl=true; // Flag of K0/aK0 oscillation
613  G4bool kf=false;
614  if(PDG==130||PDG==310)
615  {
616  kf=true;
617  if(G4UniformRand()>.5) kfl=false;
618  }
619  if ( (PDG == 2212 && F) || (PDG == 2112 && !F) ) ind=0; // pp/nn
620  else if ( (PDG == 2112 && F) || (PDG == 2212 && !F) ) ind=1; // np/pn
621  else if ( (PDG == -211 && F) || (PDG == 211 && !F) ) ind=2; // pimp/pipn
622  else if ( (PDG == 211 && F) || (PDG == -211 && !F) ) ind=3; // pipp/pimn
623  else if ( PDG == -321 || PDG == -311 || (kf && !kfl) ) ind=4; // KmN/K0N
624  else if ( PDG == 321 || PDG == 311 || (kf && kfl) ) ind=5; // KpN/aK0N
625  else if ( PDG > 3000 && PDG < 3335) ind=6; // @@ for all hyperons - take Lambda
626  else if ( PDG > -3335 && PDG < -2000) ind=7; // @@ for all anti-baryons (anti-p/anti-n)
627  else {
628  G4cout<<"*Error*G4QuasiElRatios::FetchElTot: PDG="<<PDG
629  <<", while it is defined only for p,n,hyperons,anti-baryons,pi,K/antiK"<<G4endl;
630  G4Exception("G4QuasiELRatio::FetchElTot:","22",FatalException,"QECrash");
631  }
632  if(nDB && lastItot==ind && p>0. && p==lastPtot) return lastRtot; // VI do not use toler
633  // if(nDB && lastI==ind && p>0. && std::fabs(p-lastP)/p<toler) return lastR;
634  if(p<=mip || p>=map) return CalcElTot(p,ind); // @@ Slow calculation ! (Warning?)
635  G4bool found=false;
636  G4int i=-1;
637  if(nDB) for (i=0; i<nDB; i++) if(ind==vItot[i]) // Sirch for this index in AMDB
638  {
639  found=true; // The index is found
640  break;
641  }
642  G4double lp=G4Log(p);
643  if(!nDB || !found) // Create new line in the AMDB
644  {
645  lastXtot = new std::pair<G4double,G4double>[mlp]; // Create logarithmic Table for ElTot
646  lastItot = ind; // Remember the initialized inex
647  lastKtot = static_cast<int>((lp-lpi)/dlp)+1; // MaxBin to be initialized in LogTaB
648  if(lastKtot>nlp)
649  {
650  lastKtot=nlp;
651  lastMtot=lpa-lpi;
652  }
653  else lastMtot = lastKtot*dlp; // Calculate max initialized ln(p)-lpi for LogTab
654  G4double pv=mip;
655  for(G4int j=0; j<=lastKtot; j++) // Calculate LogTab values
656  {
657  lastXtot[j]=CalcElTot(pv,ind);
658  if(j!=lastKtot) pv*=edlp;
659  }
660  i++; // Make a new record to AMDB and position on it
661  vItot.push_back(lastItot);
662  vMtot.push_back(lastMtot);
663  vKtot.push_back(lastKtot);
664  vX->push_back(lastXtot);
665  }
666  else // The A value was found in AMDB
667  {
668  lastItot=vItot[i];
669  lastMtot=vMtot[i];
670  lastKtot=vKtot[i];
671  lastXtot=(*vX)[i];
672  G4int nextK=lastKtot+1;
673  G4double lpM=lastMtot+lpi;
674  if(lp>lpM && lastKtot<nlp) // LogTab must be updated
675  {
676  lastKtot = static_cast<int>((lp-lpi)/dlp)+1; // MaxBin to be initialized in LogTab
677  if(lastKtot>nlp)
678  {
679  lastKtot=nlp;
680  lastMtot=lpa-lpi;
681  }
682  else lastMtot = lastKtot*dlp; // Calculate max initialized ln(p)-lpi for LogTab
683  G4double pv=G4Exp(lpM); // momentum of the last calculated beam
684  for(G4int j=nextK; j<=lastKtot; j++)// Calculate LogTab values
685  {
686  pv*=edlp;
687  lastXtot[j]=CalcElTot(pv,ind);
688  }
689  } // End of LogTab update
690  if(lastKtot>=nextK) // The AMDB was apdated
691  {
692  vMtot[i]=lastMtot;
693  vKtot[i]=lastKtot;
694  }
695  }
696  // Now one can use tabeles to calculate the value
697  G4double dlpp=lp-lpi; // Shifted log(p) value
698  G4int n=static_cast<int>(dlpp/dlp); // Low edge number of the bin
699  G4double d=dlpp-n*dlp; // Log shift
700  G4double e=lastXtot[n].first; // E-Base
701  lastRtot.first=e+d*(lastXtot[n+1].first-e)/dlp; // E-Result
702  if(lastRtot.first<0.) lastRtot.first = 0.;
703  G4double t=lastXtot[n].second; // T-Base
704  lastRtot.second=t+d*(lastXtot[n+1].second-t)/dlp; // T-Result
705  if(lastRtot.second<0.) lastRtot.second= 0.;
706  if(lastRtot.first>lastRtot.second) lastRtot.first = lastRtot.second;
707  return lastRtot;
708 } // End of FetchElTot
709 
710 // (Mean Elastic and Mean Total) Cross-Sections (mb) for PDG+(Z,N) at P=p[GeV/c]
711 std::pair<G4double,G4double> G4QuasiElRatios::GetElTot(G4double pIU, G4int hPDG,
712  G4int Z, G4int N)
713 {
714  G4double pGeV=pIU/gigaelectronvolt;
715  if(Z<1 && N<1)
716  {
717  G4cout<<"-Warning-G4QuasiElRatio::GetElTot:Z="<<Z<<",N="<<N<<", return zero"<<G4endl;
718  return std::make_pair(0.,0.);
719  }
720  std::pair<G4double,G4double> hp=FetchElTot(pGeV, hPDG, true);
721  std::pair<G4double,G4double> hn=FetchElTot(pGeV, hPDG, false);
722  G4double A=(Z+N)/millibarn; // To make the result in independent units(IU)
723  return std::make_pair((Z*hp.first+N*hn.first)/A,(Z*hp.second+N*hn.second)/A);
724 } // End of GetElTot
725 
726 // (Mean Elastic and Mean Total) Cross-Sections (mb) for PDG+(Z,N) at P=p[GeV/c]
727 std::pair<G4double,G4double> G4QuasiElRatios::GetChExFactor(G4double pIU, G4int hPDG,
728  G4int Z, G4int N)
729 {
730  G4double pGeV=pIU/gigaelectronvolt;
731  G4double resP=0.;
732  G4double resN=0.;
733  if(Z<1 && N<1)
734  {
735  G4cout<<"-Warning-G4QuasiElRatio::GetChExF:Z="<<Z<<",N="<<N<<", return zero"<<G4endl;
736  return std::make_pair(resP,resN);
737  }
738  G4double A=Z+N;
739  G4double pf=0.; // Possibility to interact with a proton
740  G4double nf=0.; // Possibility to interact with a neutron
741  if (hPDG==-211||hPDG==-321||hPDG==3112||hPDG==3212||hPDG==3312) pf=Z/(A+N);
742  else if(hPDG==211||hPDG==321||hPDG==3222||hPDG==3212||hPDG==3322) nf=N/(A+Z);
743  else if(hPDG==-311||hPDG==311||hPDG==130||hPDG==310)
744  {
745  G4double dA=A+A;
746  pf=Z/(dA+N+N);
747  nf=N/(dA+Z+Z);
748  }
749  G4double mult=1.; // Factor of increasing multiplicity ( ? @@)
750  if(pGeV>.5)
751  {
752  mult=1./(1.+G4Log(pGeV+pGeV))/pGeV;
753  if(mult>1.) mult=1.;
754  }
755  if(pf)
756  {
757  std::pair<G4double,G4double> hp=FetchElTot(pGeV, hPDG, true);
758  resP=pf*(hp.second/hp.first-1.)*mult;
759  }
760  if(nf)
761  {
762  std::pair<G4double,G4double> hn=FetchElTot(pGeV, hPDG, false);
763  resN=nf*(hn.second/hn.first-1.)*mult;
764  }
765  return std::make_pair(resP,resN);
766 } // End of GetChExFactor
767 
768 // scatter (pPDG,p4M) on a virtual nucleon (NPDG,N4M), result: final pair(newN4M,newp4M)
769 // if(newN4M.e()==0.) - below threshold, XS=0, no scattering of the progectile happened
770 std::pair<G4LorentzVector,G4LorentzVector> G4QuasiElRatios::Scatter(G4int NPDG,
771  G4LorentzVector N4M, G4int pPDG, G4LorentzVector p4M)
772 {
773  static const G4double mNeut= G4Neutron::Neutron()->GetPDGMass();
774  static const G4double mProt= G4Proton::Proton()->GetPDGMass();
775  static const G4double mDeut= G4Deuteron::Deuteron()->GetPDGMass();
776  static const G4double mTrit= G4Triton::Triton()->GetPDGMass();
777  static const G4double mHel3= G4He3::He3()->GetPDGMass();
778  static const G4double mAlph= G4Alpha::Alpha()->GetPDGMass();
779 
780  G4LorentzVector pr4M=p4M/megaelectronvolt; // Convert 4-momenta in MeV (keep p4M)
781  N4M/=megaelectronvolt;
782  G4LorentzVector tot4M=N4M+p4M;
783  G4double mT=mNeut;
784  G4int Z=0;
785  G4int N=1;
786  if(NPDG==2212||NPDG==90001000)
787  {
788  mT=mProt;
789  Z=1;
790  N=0;
791  }
792  else if(NPDG==90001001)
793  {
794  mT=mDeut;
795  Z=1;
796  N=1;
797  }
798  else if(NPDG==90002001)
799  {
800  mT=mHel3;
801  Z=2;
802  N=1;
803  }
804  else if(NPDG==90001002)
805  {
806  mT=mTrit;
807  Z=1;
808  N=2;
809  }
810  else if(NPDG==90002002)
811  {
812  mT=mAlph;
813  Z=2;
814  N=2;
815  }
816  else if(NPDG!=2112&&NPDG!=90000001)
817  {
818  G4cout<<"Error:G4QuasiElRatios::Scatter:NPDG="<<NPDG<<" is not 2212 or 2112"<<G4endl;
819  G4Exception("G4QuasiElRatios::Scatter:","21",FatalException,"QEcomplain");
820  //return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M);// Use this if not exception
821  }
822  G4double mT2=mT*mT;
823  G4double mP2=pr4M.m2();
824  G4double E=(tot4M.m2()-mT2-mP2)/(mT+mT);
825  G4double E2=E*E;
826  if(E<0. || E2<mP2)
827  {
828  return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); // Do Nothing Action
829  }
830  G4double P=std::sqrt(E2-mP2); // Momentum in pseudo laboratory system
831  // @@ Temporary NN t-dependence for all hadrons
832  if(pPDG>3400 || pPDG<-3400) G4cout<<"-Warning-G4QE::Scatter: pPDG="<<pPDG<<G4endl;
833  G4int PDG=2212; // *TMP* instead of pPDG
834  if(pPDG==2112||pPDG==-211||pPDG==-321) PDG=2112; // *TMP* instead of pPDG
835  if(!Z && N==1) // Change for Quasi-Elastic on neutron
836  {
837  Z=1;
838  N=0;
839  if (PDG==2212) PDG=2112;
840  else if(PDG==2112) PDG=2212;
841  }
842  G4double xSec=0.; // Prototype of Recalculated Cross Section *TMP*
843  if(PDG==2212) xSec=PCSmanager->GetChipsCrossSection(P, Z, N, PDG); // P CrossSect *TMP*
844  else xSec=NCSmanager->GetChipsCrossSection(P, Z, N, PDG); // N CrossSect *TMP*
845  // @@ check a possibility to separate p, n, or alpha (!)
846  if(xSec <= 0.) // The cross-section iz 0 -> Do Nothing
847  {
848  return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); //Do Nothing Action
849  }
850  G4double mint=0.; // Prototype of functional rand -t (MeV^2) *TMP*
851  if(PDG==2212) mint=PCSmanager->GetExchangeT(Z,N,PDG);// P functional rand -t(MeV^2) *TMP*
852  else mint=NCSmanager->GetExchangeT(Z,N,PDG);// N functional rand -t(MeV^2) *TMP*
853  G4double maxt=0.; // Prototype of max possible -t
854  if(PDG==2212) maxt=PCSmanager->GetHMaxT(); // max possible -t
855  else maxt=NCSmanager->GetHMaxT(); // max possible -t
856  G4double cost=1.-(mint+mint)/maxt; // cos(theta) in CMS
857  if(cost>1. || cost<-1. || !(cost>-1. || cost<=1.))
858  {
859  if (cost>1.) cost=1.;
860  else if(cost<-1.) cost=-1.;
861  else
862  {
863  G4double tm=0.;
864  if(PDG==2212) tm=PCSmanager->GetHMaxT();
865  else tm=NCSmanager->GetHMaxT();
866  G4cerr<<"G4QuasiFreeRatio::Scat:*NAN* cost="<<cost<<",-t="<<mint<<",tm="<<tm<<G4endl;
867  return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); // Do Nothing Action
868  }
869  }
870  G4LorentzVector reco4M=G4LorentzVector(0.,0.,0.,mT); // 4mom of the recoil nucleon
871  G4LorentzVector dir4M=tot4M-G4LorentzVector(0.,0.,0.,(tot4M.e()-mT)*.01);
872  if(!RelDecayIn2(tot4M, pr4M, reco4M, dir4M, cost, cost))
873  {
874  G4cerr<<"G4QFR::Scat:t="<<tot4M<<tot4M.m()<<",mT="<<mT<<",mP="<<std::sqrt(mP2)<<G4endl;
875  //G4Exception("G4QFR::Scat:","009",FatalException,"Decay of ElasticComp");
876  return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); // Do Nothing Action
877  }
878  return std::make_pair(reco4M*megaelectronvolt,pr4M*megaelectronvolt); // Result
879 } // End of Scatter
880 
881 // scatter (pPDG,p4M) on a virtual nucleon (NPDG,N4M), result: final pair(newN4M,newp4M)
882 // if(newN4M.e()==0.) - below threshold, XS=0, no scattering of the progectile happened
883 // User should himself change the charge (PDG) (e.g. pn->np, pi+n->pi0p, pi-p->pi0n etc.)
884 std::pair<G4LorentzVector,G4LorentzVector> G4QuasiElRatios::ChExer(G4int NPDG,
885  G4LorentzVector N4M, G4int pPDG, G4LorentzVector p4M)
886 {
887  static const G4double mNeut= G4Neutron::Neutron()->GetPDGMass();
888  static const G4double mProt= G4Proton::Proton()->GetPDGMass();
889  G4LorentzVector pr4M=p4M/megaelectronvolt; // Convert 4-momenta in MeV(keep p4M)
890  N4M/=megaelectronvolt;
891  G4LorentzVector tot4M=N4M+p4M;
892  G4int Z=0;
893  G4int N=1;
894  G4int sPDG=0; // PDG code of the scattered hadron
895  G4double mS=0.; // proto of mass of scattered hadron
896  G4double mT=mProt; // mass of the recoil nucleon
897  if(NPDG==2212)
898  {
899  mT=mNeut;
900  Z=1;
901  N=0;
902  if(pPDG==-211) sPDG=111; // pi+ -> pi0
903  else if(pPDG==-321)
904  {
905  sPDG=310; // K+ -> K0S
906  if(G4UniformRand()>.5) sPDG=130; // K+ -> K0L
907  }
908  else if(pPDG==-311||pPDG==311||pPDG==130||pPDG==310) sPDG=321; // K0 -> K+ (?)
909  else if(pPDG==3112) sPDG=3212; // Sigma- -> Sigma0
910  else if(pPDG==3212) sPDG=3222; // Sigma0 -> Sigma+
911  else if(pPDG==3312) sPDG=3322; // Xi- -> Xi0
912  }
913  else if(NPDG==2112) // Default
914  {
915  if(pPDG==211) sPDG=111; // pi+ -> pi0
916  else if(pPDG==321)
917  {
918  sPDG=310; // K+ -> K0S
919  if(G4UniformRand()>.5) sPDG=130; // K+ -> K0L
920  }
921  else if(pPDG==-311||pPDG==311||pPDG==130||pPDG==310) sPDG=-321; // K0 -> K- (?)
922  else if(pPDG==3222) sPDG=3212; // Sigma+ -> Sigma0
923  else if(pPDG==3212) sPDG=3112; // Sigma0 -> Sigma-
924  else if(pPDG==3322) sPDG=3312; // Xi0 -> Xi-
925  }
926  else
927  {
928  G4cout<<"Error:G4QuasiElRatios::ChExer: NPDG="<<NPDG<<" is not 2212 or 2112"<<G4endl;
929  G4Exception("G4QuasiElRatios::ChExer:","21",FatalException,"QE complain");
930  //return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M);// Use this if not exception
931  }
932  if(sPDG) mS=mNeut;
933  else
934  {
935  G4cout<<"Error:G4QuasiElRatios::ChExer: BAD pPDG="<<pPDG<<", NPDG="<<NPDG<<G4endl;
936  G4Exception("G4QuasiElRatios::ChExer:","21",FatalException,"QE complain");
937  //return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M);// Use this if not exception
938  }
939  G4double mT2=mT*mT;
940  G4double mS2=mS*mS;
941  G4double E=(tot4M.m2()-mT2-mS2)/(mT+mT);
942  G4double E2=E*E;
943  if(E<0. || E2<mS2)
944  {
945  return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); // Do Nothing Action
946  }
947  G4double P=std::sqrt(E2-mS2); // Momentum in pseudo laboratory system
948  // @@ Temporary NN t-dependence for all hadrons
949  G4int PDG=2212; // *TMP* instead of pPDG
950  if(pPDG==2112||pPDG==-211||pPDG==-321) PDG=2112; // *TMP* instead of pPDG
951  if(!Z && N==1) // Change for Quasi-Elastic on neutron
952  {
953  Z=1;
954  N=0;
955  if (PDG==2212) PDG=2112;
956  else if(PDG==2112) PDG=2212;
957  }
958  G4double xSec=0.; // Prototype of Recalculated Cross Section *TMP*
959  if(PDG==2212) xSec=PCSmanager->GetChipsCrossSection(P, Z, N, PDG); // P CrossSect *TMP*
960  else xSec=NCSmanager->GetChipsCrossSection(P, Z, N, PDG); // N CrossSect *TMP*
961  // @@ check a possibility to separate p, n, or alpha (!)
962  if(xSec <= 0.) // The cross-section iz 0 -> Do Nothing
963  {
964  return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); //Do Nothing Action
965  }
966  G4double mint=0.; // Prototype of functional rand -t (MeV^2) *TMP*
967  if(PDG==2212) mint=PCSmanager->GetExchangeT(Z,N,PDG);// P functional rand -t(MeV^2) *TMP*
968  else mint=NCSmanager->GetExchangeT(Z,N,PDG);// N functional rand -t(MeV^2) *TMP*
969  G4double maxt=0.; // Prototype of max possible -t
970  if(PDG==2212) maxt=PCSmanager->GetHMaxT(); // max possible -t
971  else maxt=NCSmanager->GetHMaxT(); // max possible -t
972  G4double cost=1.-mint/maxt; // cos(theta) in CMS
973  if(cost>1. || cost<-1. || !(cost>-1. || cost<=1.))
974  {
975  if (cost>1.) cost=1.;
976  else if(cost<-1.) cost=-1.;
977  else
978  {
979  G4cerr<<"G4QuasiFreeRatio::ChExer:*NAN* c="<<cost<<",t="<<mint<<",tm="<<maxt<<G4endl;
980  return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); // Do Nothing Action
981  }
982  }
983  G4LorentzVector reco4M=G4LorentzVector(0.,0.,0.,mT); // 4mom of the recoil nucleon
984  pr4M=G4LorentzVector(0.,0.,0.,mS); // 4mom of the scattered hadron
985  G4LorentzVector dir4M=tot4M-G4LorentzVector(0.,0.,0.,(tot4M.e()-mT)*.01);
986  if(!RelDecayIn2(tot4M, pr4M, reco4M, dir4M, cost, cost))
987  {
988  G4cerr<<"G4QFR::ChEx:t="<<tot4M<<tot4M.m()<<",mT="<<mT<<",mP="<<mS<<G4endl;
989  //G4Exception("G4QFR::ChExer:","009",FatalException,"Decay of ElasticComp");
990  return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); // Do Nothing Action
991  }
992  return std::make_pair(reco4M*megaelectronvolt,pr4M*megaelectronvolt); // Result
993 } // End of ChExer
994 
995 // Calculate ChEx/El ratio (p is in independent units, (Z,N) is target, pPDG is projectile)
997 {
998 
999  p/=MeV; // Converted from independent units
1000  G4double A=Z+N;
1001  if(A<1.5) return 0.;
1002  G4double Cex=0.;
1003  if (pPDG==2212) Cex=N/(A+Z);
1004  else if(pPDG==2112) Cex=Z/(A+N);
1005  else G4cout<<"*Warning*G4CohChrgExchange::ChExElCoef: wrong PDG="<<pPDG<<G4endl;
1006  Cex*=Cex; // Coherent processes squares the amplitude
1007  // @@ This is true only for nucleons: other projectiles must be treated differently
1008  G4double sp=std::sqrt(p);
1009  G4double p2=p*p;
1010  G4double p4=p2*p2;
1011  G4double dl1=G4Log(p)-5.;
1012  G4double T=(6.75+.14*dl1*dl1+13./p)/(1.+.14/p4)+.6/(p4+.00013);
1013  G4double U=(6.25+8.33e-5/p4/p)*(p*sp+.34)/p2/p;
1014  G4double R=U/T;
1015  return Cex*R*R;
1016 }
1017 
1018 // Decay of Hadron In2Particles f&s, f is in respect to the direction of HadronMomentumDir
1020  G4LorentzVector& dir, G4double maxCost, G4double minCost)
1021 {
1022  G4double fM2 = f4Mom.m2();
1023  G4double fM = std::sqrt(fM2); // Mass of the 1st Hadron
1024  G4double sM2 = s4Mom.m2();
1025  G4double sM = std::sqrt(sM2); // Mass of the 2nd Hadron
1026  G4double iM2 = theMomentum.m2();
1027  G4double iM = std::sqrt(iM2); // Mass of the decaying hadron
1028  G4double vP = theMomentum.rho(); // Momentum of the decaying hadron
1029  G4double dE = theMomentum.e(); // Energy of the decaying hadron
1030  if(dE<vP)
1031  {
1032  G4cerr<<"***G4QHad::RelDecIn2: Tachionic 4-mom="<<theMomentum<<", E-p="<<dE-vP<<G4endl;
1033  G4double accuracy=.000001*vP;
1034  G4double emodif=std::fabs(dE-vP);
1035  //if(emodif<accuracy)
1036  //{
1037  G4cerr<<"G4QHadron::RelDecIn2: *Boost* E-p shift is corrected to "<<emodif<<G4endl;
1038  theMomentum.setE(vP+emodif+.01*accuracy);
1039  //}
1040  }
1041  G4ThreeVector ltb = theMomentum.boostVector();// Boost vector for backward Lorentz Trans.
1042  G4ThreeVector ltf = -ltb; // Boost vector for forward Lorentz Trans.
1043  G4LorentzVector cdir = dir; // A copy to make a transformation to CMS
1044  cdir.boost(ltf); // Direction transpormed to CMS of the Momentum
1045  G4ThreeVector vdir = cdir.vect(); // 3-Vector of the direction-particle
1046  G4ThreeVector vx(0.,0.,1.); // Ort in the direction of the reference particle
1047  G4ThreeVector vy(0.,1.,0.); // First ort orthogonal to the direction
1048  G4ThreeVector vz(1.,0.,0.); // Second ort orthoganal to the direction
1049  if(vdir.mag2() > 0.) // the refference particle isn't at rest in CMS
1050  {
1051  vx = vdir.unit(); // Ort in the direction of the reference particle
1052  G4ThreeVector vv= vx.orthogonal(); // Not normed orthogonal vector (!)
1053  vy = vv.unit(); // First ort orthogonal to the direction
1054  vz = vx.cross(vy); // Second ort orthoganal to the direction
1055  }
1056  if(maxCost> 1.) maxCost= 1.;
1057  if(minCost<-1.) minCost=-1.;
1058  if(maxCost<-1.) maxCost=-1.;
1059  if(minCost> 1.) minCost= 1.;
1060  if(minCost> maxCost) minCost=maxCost;
1061  if(std::fabs(iM-fM-sM)<.00000001)
1062  {
1063  G4double fR=fM/iM;
1064  G4double sR=sM/iM;
1065  f4Mom=fR*theMomentum;
1066  s4Mom=sR*theMomentum;
1067  return true;
1068  }
1069  else if (iM+.001<fM+sM || iM==0.)
1070  {//@@ Later on make a quark content check for the decay
1071  G4cerr<<"***G4QH::RelDecIn2: fM="<<fM<<"+sM="<<sM<<">iM="<<iM<<",d="<<iM-fM-sM<<G4endl;
1072  return false;
1073  }
1074  G4double d2 = iM2-fM2-sM2;
1075  G4double p2 = (d2*d2/4.-fM2*sM2)/iM2; // Decay momentum(^2) in CMS of Quasmon
1076  if(p2<0.)
1077  {
1078  p2=0.;
1079  }
1080  G4double p = std::sqrt(p2);
1081  G4double ct = maxCost;
1082  if(maxCost>minCost)
1083  {
1084  G4double dcost=maxCost-minCost;
1085  ct = minCost+dcost*G4UniformRand();
1086  }
1087  G4double phi= twopi*G4UniformRand(); // @@ Change 360.*deg to M_TWOPI (?)
1088  G4double pss=0.;
1089  if(std::fabs(ct)<1.) pss = p * std::sqrt(1.-ct*ct);
1090  else
1091  {
1092  if(ct>1.) ct=1.;
1093  if(ct<-1.) ct=-1.;
1094  }
1095  G4ThreeVector pVect=(pss*std::sin(phi))*vz+(pss*std::cos(phi))*vy+p*ct*vx;
1096 
1097  f4Mom.setVect(pVect);
1098  f4Mom.setE(std::sqrt(fM2+p2));
1099  s4Mom.setVect((-1)*pVect);
1100  s4Mom.setE(std::sqrt(sM2+p2));
1101 
1102  if(f4Mom.e()+.001<f4Mom.rho())G4cerr<<"*G4QH::RDIn2:*Boost* f4M="<<f4Mom<<",e-p="
1103  <<f4Mom.e()-f4Mom.rho()<<G4endl;
1104  f4Mom.boost(ltb); // Lor.Trans. of 1st hadron back to LS
1105  if(s4Mom.e()+.001<s4Mom.rho())G4cerr<<"*G4QH::RDIn2:*Boost* s4M="<<s4Mom<<",e-p="
1106  <<s4Mom.e()-s4Mom.rho()<<G4endl;
1107  s4Mom.boost(ltb); // Lor.Trans. of 2nd hadron back to LS
1108  return true;
1109 } // End of "RelDecayIn2"
1110 
1111 
1112 
1113 
1114 
1115 
G4VCrossSectionDataSet * GetCrossSectionDataSet(const G4String &name, G4bool warning=true)
std::pair< G4double, G4double > GetRatios(G4double pIU, G4int prPDG, G4int tgZ, G4int tgN)
static G4Pow * GetInstance()
Definition: G4Pow.cc:55
G4double powA(G4double A, G4double y) const
Definition: G4Pow.hh:259
Hep3Vector boostVector() const
virtual G4double GetChipsCrossSection(G4double momentum, G4int Z, G4int N, G4int pdg)
const char * p
Definition: xmltok.h:285
static const G4double d2
static constexpr double megaelectronvolt
Definition: G4SIunits.hh:208
std::pair< G4double, G4double > FetchElTot(G4double pGeV, G4int PDG, G4bool F)
double C(double temp)
int G4int
Definition: G4Types.hh:78
std::pair< G4double, G4double > GetElTot(G4double pIU, G4int hPDG, G4int Z, G4int N)
std::pair< G4LorentzVector, G4LorentzVector > Scatter(G4int NPDG, G4LorentzVector N4M, G4int pPDG, G4LorentzVector p4M)
static double P[]
static constexpr double twopi
Definition: G4SIunits.hh:76
std::pair< G4double, G4double > GetChExFactor(G4double pIU, G4int pPDG, G4int Z, G4int N)
static const G4double dE
Hep3Vector vect() const
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
double A(double temperature)
bool G4bool
Definition: G4Types.hh:79
HepLorentzVector & boost(double, double, double)
G4double GetExchangeT(G4int tZ, G4int tN, G4int pPDG)
static const char * Default_Name()
static G4CrossSectionDataSetRegistry * Instance()
double rho() const
static G4Triton * Triton()
Definition: G4Triton.cc:95
static G4Proton * Proton()
Definition: G4Proton.cc:93
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
const G4int n
static G4Deuteron * Deuteron()
Definition: G4Deuteron.cc:94
tuple v
Definition: test.py:18
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
std::pair< G4LorentzVector, G4LorentzVector > ChExer(G4int NPDG, G4LorentzVector N4M, G4int pPDG, G4LorentzVector p4M)
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
G4double ChExElCoef(G4double p, G4int Z, G4int N, G4int pPDG)
G4bool RelDecayIn2(G4LorentzVector &theMomentum, G4LorentzVector &f4Mom, G4LorentzVector &s4Mom, G4LorentzVector &dir, G4double maxCost=1., G4double minCost=-1.)
G4double GetPDGMass() const
static constexpr double gigaelectronvolt
Definition: G4SIunits.hh:209
double m2() const
Hep3Vector unit() const
Hep3Vector orthogonal() const
double mag2() const
static const char * Default_Name()
#define G4endl
Definition: G4ios.hh:61
**D E S C R I P T I O N
Definition: HEPEvtcom.cc:77
static constexpr double MeV
Definition: G4SIunits.hh:214
static G4Alpha * Alpha()
Definition: G4Alpha.cc:89
Hep3Vector cross(const Hep3Vector &) const
std::pair< G4double, G4double > GetElTotXS(G4double Mom, G4int PDG, G4bool F)
Definition: Evaluator.cc:66
void setVect(const Hep3Vector &)
double G4double
Definition: G4Types.hh:76
static G4He3 * He3()
Definition: G4He3.cc:94
static const G4double pos
static constexpr double millibarn
Definition: G4SIunits.hh:106
virtual G4double GetChipsCrossSection(G4double momentum, G4int Z, G4int N, G4int pdg)
G4double GetExchangeT(G4int tZ, G4int tN, G4int pPDG)
G4GLOB_DLL std::ostream G4cerr
CLHEP::HepLorentzVector G4LorentzVector