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