source: trunk/source/processes/hadronic/models/chiral_inv_phase_space/cross_sections/src/G4QuasiFreeRatios.cc @ 1315

Last change on this file since 1315 was 1315, checked in by garnier, 14 years ago

update geant4-09-04-beta-cand-01 interfaces-V09-03-09 vis-V09-03-08
File size: 45.4 KB
Line 
1//
2// ********************************************************************
3// * License and Disclaimer                                           *
4// *                                                                  *
5// * The  Geant4 software  is  copyright of the Copyright Holders  of *
6// * the Geant4 Collaboration.  It is provided  under  the terms  and *
7// * conditions of the Geant4 Software License,  included in the file *
8// * LICENSE and available at  http://cern.ch/geant4/license .  These *
9// * include a list of copyright holders.                             *
10// *                                                                  *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work  make  any representation or  warranty, express or implied, *
14// * regarding  this  software system or assume any liability for its *
15// * use.  Please see the license in the file  LICENSE  and URL above *
16// * for the full disclaimer and the limitation of liability.         *
17// *                                                                  *
18// * This  code  implementation is the result of  the  scientific and *
19// * technical work of the GEANT4 collaboration.                      *
20// * By using,  copying,  modifying or  distributing the software (or *
21// * any work based  on the software)  you  agree  to acknowledge its *
22// * use  in  resulting  scientific  publications,  and indicate your *
23// * acceptance of all terms of the Geant4 Software license.          *
24// ********************************************************************
25//
26//
27// $Id: G4QuasiFreeRatios.cc,v 1.3 2010/01/22 17:02:48 mkossov Exp $
28// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
29//
30//
31// G4 Physics class: G4QuasiFreeRatios for N+A elastic cross sections
32// Created: M.V. Kossov, CERN/ITEP(Moscow), 10-OCT-01
33// The last update: M.V. Kossov, CERN/ITEP (Moscow) 15-Oct-06
34//
35//=======================================================================
36// Short description: Provides percentage of quasi-free and quasi-elastic
37// reactions in the inelastic reactions.
38// ----------------------------------------------------------------------
39
40//#define debug
41//#define pdebug
42//#define ppdebug
43//#define nandebug
44
45#include "G4QuasiFreeRatios.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
50std::vector<std::pair<G4double,G4double>*> G4QuasiFreeRatios::vX; // ETPointers to LogTable
51
52G4QuasiFreeRatios::G4QuasiFreeRatios()
53{
54#ifdef pdebug
55  G4cout<<"***^^^*** G4QuasiFreeRatios singletone is created ***^^^***"<<G4endl;
56#endif
57}
58
59G4QuasiFreeRatios::~G4QuasiFreeRatios()
60{
61  std::vector<G4double*>::iterator pos;
62  for(pos=vT.begin(); pos<vT.end(); pos++)
63  { delete [] *pos; }
64  vT.clear();
65  for(pos=vL.begin(); pos<vL.end(); pos++)
66  { delete [] *pos; }
67  vL.clear();
68
69  std::vector<std::pair<G4double,G4double>*>::iterator pos2;
70  for(pos2=vX.begin(); pos2<vX.end(); pos2++)
71  { delete [] *pos2; }
72  vX.clear();
73}
74
75// Returns Pointer to the G4VQCrossSection class
76G4QuasiFreeRatios* G4QuasiFreeRatios::GetPointer()
77{
78  static G4QuasiFreeRatios theRatios;   // *** Static body of the QEl Cross Section ***
79  return &theRatios;
80}
81
82// Calculation of pair(QuasiFree/Inelastic,QuasiElastic/QuasiFree)
83std::pair<G4double,G4double> G4QuasiFreeRatios::GetRatios(G4double pIU, G4int pPDG,
84                                                           G4int tgZ,    G4int tgN)
85{
86#ifdef pdebug
87  G4cout<<">>>IN>>>G4QFRat::GetQF:P="<<pIU<<",pPDG="<<pPDG<<",Z="<<tgZ<<",N="<<tgN<<G4endl;
88#endif
89  G4double R=0.;
90  G4double QF2In=1.;                        // Prototype of QuasiFree/Inel ratio for hN_tot
91  G4int tgA=tgZ+tgN;
92  if(tgA<2) return std::make_pair(QF2In,R); // No quasi-elastic on the only nucleon
93  std::pair<G4double,G4double> ElTot=GetElTot(pIU, pPDG, tgZ, tgN); // mean hN El&Tot(IU)
94  //if( ( (pPDG>999 && pIU<227.) || pIU<27.) && tgA>1) R=1.; // @@ TMP to accelerate @lowE
95  if(pPDG>999 && pIU<227. && tgZ+tgN>1) R=1.;                // To accelerate @lowE
96  else if(ElTot.second>0.)
97  {
98    R=ElTot.first/ElTot.second;             // El/Total ratio (does not depend on units
99    QF2In=GetQF2IN_Ratio(ElTot.second/millibarn, tgZ+tgN);   // QuasiFree/Inelastic ratio
100  }
101#ifdef pdebug
102  G4cout<<">>>OUT>>>G4QuasiFreeRatio::GetQF2IN_Ratio: QF2In="<<QF2In<<", R="<<R<<G4endl;
103#endif
104  return std::make_pair(QF2In,R);
105}
106
107// Calculatio QasiFree/Inelastic Ratio as a function of total hN cross-section (mb) and A
108G4double G4QuasiFreeRatios::GetQF2IN_Ratio(G4double s, G4int A)
109{
110  static const G4int    nps=150;        // Number of steps in the R(s) LinTable
111  static const G4int    mps=nps+1;      // Number of elements in the R(s) LinTable
112  static const G4double sma=150.;       // The first LinTabEl(s=0)=1., s>sma -> logTab
113  static const G4double ds=sma/nps;     // Step of the linear Table
114  static const G4int    nls=100;        // Number of steps in the R(lns) logTable
115  static const G4int    mls=nls+1;      // Number of elements in the R(lns) logTable
116  static const G4double lsi=5.;         // The min ln(s) logTabEl(s=148.4 < sma=150.)
117  static const G4double lsa=9.;         // The max ln(s) logTabEl(s=148.4 - 8103. mb)
118  static const G4double mi=std::exp(lsi);// The min s of logTabEl(~ 148.4 mb)
119  static const G4double ms=std::exp(lsa);// The max s of logTabEl(~ 8103. mb)
120  static const G4double dl=(lsa-lsi)/nls;// Step of the logarithmic Table
121  static const G4double edl=std::exp(dl);// Multiplication step of the logarithmic Table
122  static const G4double toler=.01;      // The tolarence mb defining the same cross-section
123  static G4double lastS=0.;             // The last sigma value for which R was calculated
124  static G4double lastR=0.;             // The last ratio R which was calculated
125  // Local Associative Data Base:
126  static std::vector<G4int>     vA;     // Vector of calculated A
127  static std::vector<G4double>  vH;     // Vector of max s initialized in the LinTable
128  static std::vector<G4int>     vN;     // Vector of topBin number initialized in LinTable
129  static std::vector<G4double>  vM;     // Vector of rel max ln(s) initialized in LogTable
130  static std::vector<G4int>     vK;     // Vector of topBin number initialized in LogTable
131  // Last values of the Associative Data Base:
132  static G4int     lastA=0;             // theLast of calculated A
133  static G4double  lastH=0.;            // theLast of max s initialized in the LinTable
134  static G4int     lastN=0;             // theLast of topBin number initialized in LinTable
135  static G4double  lastM=0.;            // theLast of rel max ln(s) initialized in LogTable
136  static G4int     lastK=0;             // theLast of topBin number initialized in LogTable
137  static G4double* lastT=0;             // theLast of pointer to LinTable in the C++ heap
138  static G4double* lastL=0;             // theLast of pointer to LogTable in the C++ heap
139  // LogTable is created only if necessary. The ratio R(s>8100 mb) = 0 for any nuclei
140#ifdef pdebug
141  G4cout<<"+++G4QuasiFreeRatio::GetQF2IN_Ratio:A="<<A<<", s="<<s<<G4endl;
142#endif
143  if(s<toler || A<2) return 1.;
144  if(s>ms) return 0.;
145  if(A>238)
146  {
147    G4cout<<"-Warning-G4QuasiFreeRatio::GetQF2IN_Ratio:A="<<A<<">238, return zero"<<G4endl;
148    return 0.;
149  }
150  G4int nDB=vA.size();                  // A number of nuclei already initialized in AMDB
151  //  if(nDB && lastA==A && std::fabs(s-lastS)<toler) return lastR;
152  if(nDB && lastA==A && s==lastS) return lastR;  // VI do not use tolerance
153  G4bool found=false;
154  G4int i=-1;
155  if(nDB) for (i=0; i<nDB; i++) if(A==vA[i]) // Sirch for this A in AMDB
156  {
157    found=true;                         // The A value is found
158    break;
159  }
160#ifdef pdebug
161  G4cout<<"+++G4QuasiFreeRatio::GetQF2IN_Ratio: nDB="<<nDB<<", found="<<found<<G4endl;
162#endif
163  if(!nDB || !found)                    // Create new line in the AMDB
164  {
165    lastA = A;
166#ifdef pdebug
167    G4cout<<"G4QuasiFreeRatios::GetQF2IN_Ratio: NewT, A="<<A<<", nDB="<<nDB<<G4endl;
168#endif
169    lastT = new G4double[mps];          // Create the linear Table
170    lastN = static_cast<int>(s/ds)+1;   // MaxBin to be initialized
171    if(lastN>nps)
172    {
173      lastN=nps;
174      lastH=sma;
175    }
176    else lastH = lastN*ds;              // Calculate max initialized s for LinTab
177    G4double sv=0;
178    lastT[0]=1.;
179    for(G4int j=1; j<=lastN; j++)       // Calculate LogTab values
180    {
181      sv+=ds;
182      lastT[j]=CalcQF2IN_Ratio(sv,A);
183    }
184    lastL=new G4double[mls];          // Create the logarithmic Table
185    if(s>sma)                           // Initialize the logarithmic Table
186    {
187#ifdef pdebug
188      G4cout<<"G4QuasiFreeRatios::GetQF2IN_Ratio: NewL, A="<<A<<", nDB="<<nDB<<G4endl;
189#endif
190      G4double ls=std::log(s);
191      lastK = static_cast<int>((ls-lsi)/dl)+1; // MaxBin to be initialized in LogTaB
192      if(lastK>nls)
193      {
194        lastK=nls;
195        lastM=lsa-lsi;
196      }
197      else lastM = lastK*dl;            // Calculate max initialized ln(s)-lsi for LogTab
198      sv=mi;
199      for(G4int j=0; j<=lastK; j++)     // Calculate LogTab values
200      {
201        lastL[j]=CalcQF2IN_Ratio(sv,A);
202        if(j!=lastK) sv*=edl;
203      }
204    }
205    else                                // LogTab is not initialized
206    {
207      lastK = 0;
208      lastM = 0.;
209    }
210    i++;                                // Make a new record to AMDB and position on it
211    vA.push_back(lastA);
212    vH.push_back(lastH);
213    vN.push_back(lastN);
214    vM.push_back(lastM);
215    vK.push_back(lastK);
216    vT.push_back(lastT);
217    vL.push_back(lastL);
218  }
219  else                                  // The A value was found in AMDB
220  {
221    lastA=vA[i];
222    lastH=vH[i];
223    lastN=vN[i];
224    lastM=vM[i];
225    lastK=vK[i];
226    lastT=vT[i];
227    lastL=vL[i];
228#ifdef pdebug
229    G4cout<<"G4QuasiFreeRatios::GetQF2IN_Ratio: Found, s="<<s<<", lastM="<<lastM<<G4endl;
230#endif
231    if(s>lastM)                          // At least LinTab must be updated
232    {
233      G4int nextN=lastN+1;               // The next bin to be initialized
234#ifdef pdebug
235      G4cout<<"G4QuasiFreeRatios::GetQF2IN_Ratio: lastN="<<lastN<<" ?< nps="<<nps<<G4endl;
236#endif
237      if(lastN<nps)
238      {
239        lastN = static_cast<int>(s/ds)+1;// MaxBin to be initialized
240        if(lastN>nps)
241        {
242          lastN=nps;
243          lastH=sma;
244        }
245        else lastH = lastN*ds;           // Calculate max initialized s for LinTab
246        G4double sv=lastM;
247        for(G4int j=nextN; j<=lastN; j++)// Calculate LogTab values
248        {
249          sv+=ds;
250          lastT[j]=CalcQF2IN_Ratio(sv,A);
251        }
252#ifdef pdebug
253        G4cout<<"G4QuasiFreeRatios::GetQF2IN_Ratio: End of LinTab update"<<G4endl;
254#endif
255      } // End of LinTab update
256#ifdef pdebug
257      G4cout<<"G4QFRatios::GetQF2IN_Ratio: lN="<<lastN<<", nN="<<nextN<<", i="<<i<<G4endl;
258#endif
259      if(lastN>=nextN)
260      {
261        vH[i]=lastH;
262        vN[i]=lastN;
263      }
264      G4int nextK=lastK+1;
265      if(!lastK) nextK=0;
266#ifdef pdebug
267      G4cout<<"G4QFRat::GetQF2IN_Ratio: sma="<<sma<<", lastK="<<lastK<<" < "<<nls<<G4endl;
268#endif
269      if(s>sma && lastK<nls)             // LogTab must be updated
270      {
271        G4double sv=std::exp(lastM+lsi); // Define starting poit (lastM will be changed)
272        G4double ls=std::log(s);
273        lastK = static_cast<int>((ls-lsi)/dl)+1; // MaxBin to be initialized in LogTaB
274        if(lastK>nls)
275        {
276          lastK=nls;
277          lastM=lsa-lsi;
278        }
279        else lastM = lastK*dl;           // Calculate max initialized ln(s)-lsi for LogTab
280#ifdef pdebug
281        G4cout<<"G4QFRat::GetQF2IN_Ratio: nK="<<nextK<<", lK="<<lastK<<", sv="<<sv<<G4endl;
282#endif
283        for(G4int j=nextK; j<=lastK; j++)// Calculate LogTab values
284        {
285          sv*=edl;
286#ifdef pdebug
287          G4cout<<"G4QFRat::GetQF2IN_Ratio: j="<<j<<", sv="<<sv<<", A="<<A<<G4endl;
288#endif
289          lastL[j]=CalcQF2IN_Ratio(sv,A);
290        }
291#ifdef pdebug
292        G4cout<<"G4QuasiFreeRatios::GetQF2IN_Ratio: End of LinTab update"<<G4endl;
293#endif
294      } // End of LogTab update
295#ifdef pdebug
296      G4cout<<"G4QFRatios::GetQF2IN_Ratio: lK="<<lastK<<", nK="<<nextK<<", i="<<i<<G4endl;
297#endif
298      if(lastK>=nextK)
299      {
300        vM[i]=lastM;
301        vK[i]=lastK;
302      }
303    }
304  }
305#ifdef pdebug
306  G4cout<<"G4QuasiFreeRatios::GetQF2IN_Ratio: BeforeTab s="<<s<<", sma="<<sma<<G4endl;
307#endif
308  // Now one can use tabeles to calculate the value
309  if(s<sma)                             // Use linear table
310  {
311    G4int n=static_cast<int>(s/ds);     // Low edge number of the bin
312    G4double d=s-n*ds;                  // Linear shift
313    G4double v=lastT[n];                // Base
314    lastR=v+d*(lastT[n+1]-v)/ds;        // Result
315  }
316  else                                  // Use log table
317  {
318    G4double ls=std::log(s)-lsi;        // ln(s)-l_min
319    G4int n=static_cast<int>(ls/dl);    // Low edge number of the bin
320    G4double d=ls-n*dl;                 // Log shift
321    G4double v=lastL[n];                // Base
322    lastR=v+d*(lastL[n+1]-v)/dl;        // Result
323  }
324  if(lastR<0.) lastR=0.;
325  if(lastR>1.) lastR=1.;
326#ifdef pdebug
327  G4cout<<"G4QuasiFreeRatios::GetQF2IN_Ratio: BeforeRet lastR="<<lastR<<G4endl;
328#endif
329  return lastR;
330} // End of CalcQF2IN_Ratio
331
332// Calculatio QasiFree/Inelastic Ratio as a function of total hN cross-section and A
333G4double G4QuasiFreeRatios::CalcQF2IN_Ratio(G4double s, G4int A)
334{
335  static const G4double C=1.246;
336  G4double s2=s*s;
337  G4double s4=s2*s2;
338  G4double ss=std::sqrt(std::sqrt(s));
339  G4double P=7.48e-5*s2/(1.+8.77e12/s4/s4/s2);
340  G4double E=.2644+.016/(1.+std::exp((29.54-s)/2.49));
341  G4double F=ss*.1526*std::exp(-s2*ss*.0000859);
342  return C*std::exp(-E*std::pow(G4double(A-1.),F))/std::pow(G4double(A),P);
343} // End of CalcQF2IN_Ratio
344
345// Calculatio pair(hN_el,hN_tot) (mb): p in GeV/c, index(PDG,F) (see FetchElTot)
346std::pair<G4double,G4double> G4QuasiFreeRatios::CalcElTot(G4double p, G4int I)
347{
348  // ---------> Each parameter set can have not more than nPoints=128 parameters
349  static const G4double lmi=3.5;       // min of (lnP-lmi)^2 parabola
350  static const G4double pbe=.0557;     // elastic (lnP-lmi)^2 parabola coefficient
351  static const G4double pbt=.3;        // total (lnP-lmi)^2 parabola coefficient
352  static const G4double pmi=.1;        // Below that fast LE calculation is made
353  static const G4double pma=1000.;     // Above that fast HE calculation is made
354  // ======================================================================================
355  G4double El=0.;                      // prototype of the elastic hN cross-section
356  G4double To=0.;                      // prototype of the total hN cross-section
357  if(p<=0.)
358  {
359    G4cout<<"-Warning-G4QuasiFreeRatios::CalcElTot: p="<<p<<" is zero or negative"<<G4endl;
360    return std::make_pair(El,To);
361  }
362  if     (!I)                          // pp/nn
363  {
364#ifdef debug
365    G4cout<<"G4QuasiFreeR::CalcElTot:I=0, p="<<p<<", pmi="<<pmi<<", pma="<<pma<<G4endl;
366#endif
367    if(p<pmi)
368    {
369      G4double p2=p*p;
370      El=1./(.00012+p2*.2);
371      To=El;
372#ifdef debug
373      G4cout<<"G4QuasiFreeR::CalcElTot:I=0i, El="<<El<<", To="<<To<<", p2="<<p2<<G4endl;
374#endif
375    }
376    else if(p>pma)
377    {
378      G4double lp=std::log(p)-lmi;
379      G4double lp2=lp*lp;
380      El=pbe*lp2+6.72;
381      To=pbt*lp2+38.2;
382#ifdef debug
383      G4cout<<"G4QuasiFreeR::CalcElTot:I=0a, El="<<El<<", To="<<To<<", lp2="<<lp2<<G4endl;
384#endif
385    }
386    else
387    {
388      G4double p2=p*p;
389      G4double LE=1./(.00012+p2*.2);
390      G4double lp=std::log(p)-lmi;
391      G4double lp2=lp*lp;
392      G4double rp2=1./p2;
393      El=LE+(pbe*lp2+6.72+32.6/p)/(1.+rp2/p);
394      To=LE+(pbt*lp2+38.2+52.7*rp2)/(1.+2.72*rp2*rp2);
395#ifdef debug
396      G4cout<<"G4QuasiFreeR::CalcElTot:0,E="<<El<<",T="<<To<<",s="<<p2<<",l="<<lp2<<G4endl;
397#endif
398    }
399  }
400  else if(I==1)                        // np/pn
401  {
402    if(p<pmi)
403    {
404      G4double p2=p*p;
405      El=1./(.00012+p2*(.051+.1*p2));
406      To=El;
407    }
408    else if(p>pma)
409    {
410      G4double lp=std::log(p)-lmi;
411      G4double lp2=lp*lp;
412      El=pbe*lp2+6.72;
413      To=pbt*lp2+38.2;
414    }
415    else
416    {
417      G4double p2=p*p;
418      G4double LE=1./(.00012+p2*(.051+.1*p2));
419      G4double lp=std::log(p)-lmi;
420      G4double lp2=lp*lp;
421      G4double rp2=1./p2;
422      El=LE+(pbe*lp2+6.72+30./p)/(1.+.49*rp2/p);
423      To=LE+(pbt*lp2+38.2)/(1.+.54*rp2*rp2);
424    }
425  }
426  else if(I==2)                        // pimp/pipn
427  {
428    G4double lp=std::log(p);
429    if(p<pmi)
430    {
431      G4double lr=lp+1.27;
432      El=1.53/(lr*lr+.0676);
433      To=El*3;
434    }
435    else if(p>pma)
436    {
437      G4double ld=lp-lmi;
438      G4double ld2=ld*ld;
439      G4double sp=std::sqrt(p);
440      El=pbe*ld2+2.4+7./sp;
441      To=pbt*ld2+22.3+12./sp;
442    }
443    else
444    {
445      G4double lr=lp+1.27;                    // p1
446      G4double LE=1.53/(lr*lr+.0676);         // p2, p3       
447      G4double ld=lp-lmi;                     // p4 (lmi=3.5)
448      G4double ld2=ld*ld;
449      G4double p2=p*p;
450      G4double p4=p2*p2;
451      G4double sp=std::sqrt(p);
452      G4double lm=lp+.36;                     // p5
453      G4double md=lm*lm+.04;                  // p6
454      G4double lh=lp-.017;                    // p7
455      G4double hd=lh*lh+.0025;                // p8
456      El=LE+(pbe*ld2+2.4+7./sp)/(1.+.7/p4)+.6/md+.05/hd;//p9(pbe=.0557),p10,p11,p12,p13,p14
457      To=LE*3+(pbt*ld2+22.3+12./sp)/(1.+.4/p4)+1./md+.06/hd;
458    }
459  }
460  else if(I==3)                        // pipp/pimn
461  {
462    G4double lp=std::log(p);
463    if(p<pmi)
464    {
465      G4double lr=lp+1.27;
466      G4double lr2=lr*lr;
467      El=13./(lr2+lr2*lr2+.0676);
468      To=El;
469    }
470    else if(p>pma)
471    {
472      G4double ld=lp-lmi;
473      G4double ld2=ld*ld;
474      G4double sp=std::sqrt(p);
475      El=pbe*ld2+2.4+6./sp;
476      To=pbt*ld2+22.3+5./sp;
477    }
478    else
479    {
480      G4double lr=lp+1.27;                   // p1
481      G4double lr2=lr*lr;
482      G4double LE=13./(lr2+lr2*lr2+.0676);   // p2, p3
483      G4double ld=lp-lmi;                    // p4 (lmi=3.5)
484      G4double ld2=ld*ld;
485      G4double p2=p*p;
486      G4double p4=p2*p2;
487      G4double sp=std::sqrt(p);
488      G4double lm=lp-.32;                    // p5
489      G4double md=lm*lm+.0576;               // p6
490      El=LE+(pbe*ld2+2.4+6./sp)/(1.+3./p4)+.7/md; // p7(pbe=.0557), p8, p9, p10, p11
491      To=LE+(pbt*ld2+22.3+5./sp)/(1.+1./p4)+.8/md;
492    }
493  }
494  else if(I==4)                        // Kmp/Kmn/K0p/K0n
495  {
496
497    if(p<pmi)
498    {
499      G4double psp=p*std::sqrt(p);
500      El=5.2/psp;
501      To=14./psp;
502    }
503    else if(p>pma)
504    {
505      G4double ld=std::log(p)-lmi;
506      G4double ld2=ld*ld;
507      El=pbe*ld2+2.23;
508      To=pbt*ld2+19.5;
509    }
510    else
511    {
512      G4double ld=std::log(p)-lmi;
513      G4double ld2=ld*ld;
514      G4double sp=std::sqrt(p);
515      G4double psp=p*sp;
516      G4double p2=p*p;
517      G4double p4=p2*p2;
518      G4double lm=p-.39;
519      G4double md=lm*lm+.000156;
520      G4double lh=p-1.;
521      G4double hd=lh*lh+.0156;
522      El=5.2/psp+(pbe*ld2+2.23)/(1.-.7/sp+.075/p4)+.004/md+.15/hd;
523      To=14./psp+(pbt*ld2+19.5)/(1.-.21/sp+.52/p4)+.006/md+.30/hd;
524    }
525  }
526  else if(I==5)                        // Kpp/Kpn/aKp/aKn
527  {
528    if(p<pmi)
529    {
530      G4double lr=p-.38;
531      G4double lm=p-1.;
532      G4double md=lm*lm+.372;   
533      El=.7/(lr*lr+.0676)+2./md;
534      To=El+.6/md;
535    }
536    else if(p>pma)
537    {
538      G4double ld=std::log(p)-lmi;
539      G4double ld2=ld*ld;
540      El=pbe*ld2+2.23;
541      To=pbt*ld2+19.5;
542    }
543    else
544    {
545      G4double ld=std::log(p)-lmi;
546      G4double ld2=ld*ld;
547      G4double lr=p-.38;
548      G4double LE=.7/(lr*lr+.0676);
549      G4double sp=std::sqrt(p);
550      G4double p2=p*p;
551      G4double p4=p2*p2;
552      G4double lm=p-1.;
553      G4double md=lm*lm+.372;
554      El=LE+(pbe*ld2+2.23)/(1.-.7/sp+.1/p4)+2./md;
555      To=LE+(pbt*ld2+19.5)/(1.+.46/sp+1.6/p4)+2.6/md;
556    }
557  }
558  else if(I==6)                        // hyperon-N
559  {
560    if(p<pmi)
561    {
562      G4double p2=p*p;
563      El=1./(.002+p2*(.12+p2));
564      To=El;
565    }
566    else if(p>pma)
567    {
568      G4double lp=std::log(p)-lmi;
569      G4double lp2=lp*lp;
570      G4double sp=std::sqrt(p);
571      El=(pbe*lp2+6.72)/(1.+2./sp);
572      To=(pbt*lp2+38.2+900./sp)/(1.+27./sp);
573    }
574    else
575    {
576      G4double p2=p*p;
577      G4double LE=1./(.002+p2*(.12+p2));
578      G4double lp=std::log(p)-lmi;
579      G4double lp2=lp*lp;
580      G4double p4=p2*p2;
581      G4double sp=std::sqrt(p);
582      El=LE+(pbe*lp2+6.72+99./p2)/(1.+2./sp+2./p4);
583      To=LE+(pbt*lp2+38.2+900./sp)/(1.+27./sp+3./p4);
584    }
585  }
586  else if(I==7)                        // antibaryon-N
587  {
588    if(p>pma)
589    {
590      G4double lp=std::log(p)-lmi;
591      G4double lp2=lp*lp;
592      El=pbe*lp2+6.72;
593      To=pbt*lp2+38.2;
594    }
595    else
596    {
597      G4double ye=std::pow(p,1.25);
598      G4double yt=std::pow(p,.35);
599      G4double lp=std::log(p)-lmi;
600      G4double lp2=lp*lp;
601      El=80./(ye+1.)+pbe*lp2+6.72;
602      To=(80./yt+.3)/yt+pbt*lp2+38.2;
603    }
604  }
605  else
606  {
607    G4cout<<"*Error*G4QuasiFreeRatios::CalcElTot:ind="<<I<<" is not defined (0-7)"<<G4endl;
608    G4Exception("G4QuasiFreeRatios::CalcElTot:","23",FatalException,"CHIPScrash");
609  }
610  if(El>To) El=To;
611  return std::make_pair(El,To);
612} // End of CalcElTot
613
614// For hadron PDG with momentum Mom (GeV/c) on N (p/n) calculate <sig_el,sig_tot> pair (mb)
615std::pair<G4double,G4double> G4QuasiFreeRatios::GetElTotXS(G4double p, G4int PDG, G4bool F)
616{
617  G4int ind=0;                                 // Prototype of the reaction index
618  G4bool kfl=true;                             // Flag of K0/aK0 oscillation
619  G4bool kf=false;
620  if(PDG==130||PDG==310)
621  {
622    kf=true;
623    if(G4UniformRand()>.5) kfl=false;
624  }
625  if      ( (PDG == 2212 && F) || (PDG == 2112 && !F) ) ind=0; // pp/nn
626  else if ( (PDG == 2112 && F) || (PDG == 2212 && !F) ) ind=1; // np/pn
627  else if ( (PDG == -211 && F) || (PDG == 211 && !F) ) ind=2; // pimp/pipn
628  else if ( (PDG == 211 && F) || (PDG == -211 && !F) ) ind=3; // pipp/pimn
629  else if ( PDG == -321 || PDG == -311 || (kf && !kfl) ) ind=4; // KmN/K0N
630  else if ( PDG == 321 || PDG == 311 || (kf && kfl) ) ind=5; // KpN/aK0N
631  else if ( PDG >  3000 && PDG <  3335) ind=6; // @@ for all hyperons - take Lambda
632  else if ( PDG > -3335 && PDG < -2000) ind=7; // @@ for all anti-baryons (anti-p/anti-n)
633  else {
634    G4cout<<"*Error*G4QuasiFreeRatios::CalcElTotXS: PDG="<<PDG
635          <<", while it is defined only for p,n,hyperons,anti-baryons,pi,K/antiK"<<G4endl;
636    G4Exception("G4QuasiFreeRatio::CalcElTotXS:","22",FatalException,"CHIPScrash");
637  }
638  return CalcElTot(p,ind);
639}
640
641// Calculatio pair(hN_el,hN_tot)(mb): p in GeV/c, F=true -> N=proton, F=false -> N=neutron
642std::pair<G4double,G4double> G4QuasiFreeRatios::FetchElTot(G4double p, G4int PDG, G4bool F)
643{
644  static const G4int    nlp=300;         // Number of steps in the S(lnp) logTable(5% step)
645  static const G4int    mlp=nlp+1;       // Number of elements in the S(lnp) logTable
646  static const G4double lpi=-5.;         // The min ln(p) logTabEl(p=6.7 MeV/c - 22. TeV/c)
647  static const G4double lpa=10.;         // The max ln(p) logTabEl(p=6.7 MeV/c - 22. TeV/c)
648  static const G4double mi=std::exp(lpi);// The min p of logTabEl(~ 6.7 MeV/c)
649  static const G4double ma=std::exp(lpa);// The max p of logTabEl(~ 22. TeV)
650  static const G4double dl=(lpa-lpi)/nlp;// Step of the logarithmic Table
651  static const G4double edl=std::exp(dl);// Multiplication step of the logarithmic Table
652  //static const G4double toler=.001;      // Relative Tolarence defining "theSameMomentum"
653  static G4double lastP=0.;              // The last momentum for which XS was calculated
654  static G4int    lastH=0;               // The last projPDG for which XS was calculated
655  static G4bool   lastF=true;            // The last nucleon for which XS was calculated
656  static std::pair<G4double,G4double> lastR=std::make_pair(0.,0.); // The last result
657  // Local Associative Data Base:
658  static std::vector<G4int>     vI;      // Vector of index for which XS was calculated
659  static std::vector<G4double>  vM;      // Vector of rel max ln(p) initialized in LogTable
660  static std::vector<G4int>     vK;      // Vector of topBin number initialized in LogTable
661  // Last values of the Associative Data Base:
662  static G4int     lastI=0;              // The Last index for which XS was calculated
663  static G4double  lastM=0.;             // The Last rel max ln(p) initialized in LogTable
664  static G4int     lastK=0;             // The Last topBin number initialized in LogTable
665  static std::pair<G4double,G4double>* lastX=0; // The Last ETPointers to LogTable in heap
666  // LogTable is created only if necessary. The ratio R(s>8100 mb) = 0 for any nuclei
667  G4int nDB=vI.size();                   // A number of hadrons already initialized in AMDB
668#ifdef pdebug
669  G4cout<<"G4QuasiFreeR::FetchElTot:p="<<p<<",PDG="<<PDG<<",F="<<F<<",nDB="<<nDB<<G4endl;
670#endif
671  if(nDB && lastH==PDG && lastF==F && p>0. && p==lastP) return lastR;// VI don't use toler.
672  //  if(nDB && lastH==PDG && lastF==F && p>0. && std::fabs(p-lastP)/p<toler) return lastR;
673  lastH=PDG;
674  lastF=F;
675  G4int ind=-1;                          // Prototipe of the index of the PDG/F combination
676  // i=0: pp(nn), i=1: np(pn), i=2: pimp(pipn), i=3: pipp(pimn), i=4: Kmp(Kmn,K0n,K0p),
677  // i=5: Kpp(Kpn,aK0n,aK0p), i=6: Hp(Hn), i=7: app(apn,ann,anp)
678  G4bool kfl=true;                             // Flag of K0/aK0 oscillation
679  G4bool kf=false;
680  if(PDG==130||PDG==310)
681  {
682    kf=true;
683    if(G4UniformRand()>.5) kfl=false;
684  }
685  if      ( (PDG == 2212 && F) || (PDG == 2112 && !F) ) ind=0; // pp/nn
686  else if ( (PDG == 2112 && F) || (PDG == 2212 && !F) ) ind=1; // np/pn
687  else if ( (PDG == -211 && F) || (PDG == 211 && !F) ) ind=2; // pimp/pipn
688  else if ( (PDG == 211 && F) || (PDG == -211 && !F) ) ind=3; // pipp/pimn
689  else if ( PDG == -321 || PDG == -311 || (kf && !kfl) ) ind=4; // KmN/K0N
690  else if ( PDG == 321 || PDG == 311 || (kf && kfl) ) ind=5; // KpN/aK0N
691  else if ( PDG >  3000 && PDG <  3335) ind=6; // @@ for all hyperons - take Lambda
692  else if ( PDG > -3335 && PDG < -2000) ind=7; // @@ for all anti-baryons (anti-p/anti-n)
693  else {
694    G4cout<<"*Error*G4QuasiFreeRatios::FetchElTot: PDG="<<PDG
695          <<", while it is defined only for p,n,hyperons,anti-baryons,pi,K/antiK"<<G4endl;
696    G4Exception("G4QuasiFreeRatio::FetchElTot:","22",FatalException,"CHIPScrash");
697  }
698  if(nDB && lastI==ind && p>0. && p==lastP) return lastR;  // VI do not use toler
699  //  if(nDB && lastI==ind && p>0. && std::fabs(p-lastP)/p<toler) return lastR;
700  if(p<=mi || p>=ma) return CalcElTot(p,ind);   // @@ Slow calculation ! (Warning?)
701  G4bool found=false;
702  G4int i=-1;
703  if(nDB) for (i=0; i<nDB; i++) if(ind==vI[i])  // Sirch for this index in AMDB
704  {
705    found=true;                                 // The index is found
706    break;
707  }
708  G4double lp=std::log(p);
709#ifdef pdebug
710  G4cout<<"G4QuasiFreeR::FetchElTot:I="<<ind<<",i="<<i<<",fd="<<found<<",lp="<<lp<<G4endl;
711#endif
712  if(!nDB || !found)                            // Create new line in the AMDB
713  {
714#ifdef pdebug
715    G4cout<<"G4QuasiFreeRatios::FetchElTot: NewX, ind="<<ind<<", nDB="<<nDB<<G4endl;
716#endif
717    lastX = new std::pair<G4double,G4double>[mlp]; // Create logarithmic Table for ElTot
718    lastI = ind;                                // Remember the initialized inex
719    lastK = static_cast<int>((lp-lpi)/dl)+1;    // MaxBin to be initialized in LogTaB
720    if(lastK>nlp)
721    {
722      lastK=nlp;
723      lastM=lpa-lpi;
724    }
725    else lastM = lastK*dl;               // Calculate max initialized ln(p)-lpi for LogTab
726    G4double pv=mi;
727    for(G4int j=0; j<=lastK; j++)        // Calculate LogTab values
728    {
729      lastX[j]=CalcElTot(pv,ind);
730#ifdef pdebug
731      G4cout<<"G4QuasiFreeR::FetchElTot:I,j="<<j<<",pv="<<pv<<",E="<<lastX[j].first<<",T="
732            <<lastX[j].second<<G4endl;
733#endif
734      if(j!=lastK) pv*=edl;
735    }
736    i++;                                 // Make a new record to AMDB and position on it
737    vI.push_back(lastI);
738    vM.push_back(lastM);
739    vK.push_back(lastK);
740    vX.push_back(lastX);
741  }
742  else                                   // The A value was found in AMDB
743  {
744    lastI=vI[i];
745    lastM=vM[i];
746    lastK=vK[i];
747    lastX=vX[i];
748    G4int nextK=lastK+1;
749    G4double lpM=lastM+lpi;
750#ifdef pdebug
751    G4cout<<"G4QuasiFreeR::FetchElTo:M="<<lpM<<",l="<<lp<<",K="<<lastK<<",n="<<nlp<<G4endl;
752#endif
753    if(lp>lpM && lastK<nlp)              // LogTab must be updated
754    {
755      lastK = static_cast<int>((lp-lpi)/dl)+1; // MaxBin to be initialized in LogTab
756#ifdef pdebug
757      G4cout<<"G4QuasiFreeR::FetET:K="<<lastK<<",lp="<<lp<<",li="<<lpi<<",dl="<<dl<<G4endl;
758#endif
759      if(lastK>nlp)
760      {
761        lastK=nlp;
762        lastM=lpa-lpi;
763      }
764      else lastM = lastK*dl;           // Calculate max initialized ln(p)-lpi for LogTab
765      G4double pv=std::exp(lpM);       // momentum of the last calculated beam
766      for(G4int j=nextK; j<=lastK; j++)// Calculate LogTab values
767      {
768        pv*=edl;
769        lastX[j]=CalcElTot(pv,ind);
770#ifdef pdebug
771        G4cout<<"G4QuasiFreeR::FetchElTot:U:j="<<j<<",p="<<pv<<",E="<<lastX[j].first<<",T="
772              <<lastX[j].second<<G4endl;
773#endif
774      }
775    } // End of LogTab update
776    if(lastK>=nextK)                   // The AMDB was apdated
777    {
778      vM[i]=lastM;
779      vK[i]=lastK;
780    }
781  }
782  // Now one can use tabeles to calculate the value
783  G4double dlp=lp-lpi;                       // Shifted log(p) value
784  G4int n=static_cast<int>(dlp/dl);          // Low edge number of the bin
785  G4double d=dlp-n*dl;                       // Log shift
786  G4double e=lastX[n].first;                 // E-Base
787  lastR.first=e+d*(lastX[n+1].first-e)/dl;   // E-Result
788  if(lastR.first<0.)  lastR.first = 0.;
789  G4double t=lastX[n].second;                // T-Base
790  lastR.second=t+d*(lastX[n+1].second-t)/dl; // T-Result
791  if(lastR.second<0.) lastR.second= 0.;
792#ifdef pdebug
793  G4cout<<"=O=>G4QuasiFreeR::FetchElTot:1st="<<lastR.first<<", 2nd="<<lastR.second<<G4endl;
794#endif
795  if(lastR.first>lastR.second) lastR.first = lastR.second;
796  return lastR;
797} // End of FetchElTot
798
799// (Mean Elastic and Mean Total) Cross-Sections (mb) for PDG+(Z,N) at P=p[GeV/c]
800std::pair<G4double,G4double> G4QuasiFreeRatios::GetElTot(G4double pIU, G4int hPDG,
801                                                         G4int Z,       G4int N)
802{
803  G4double pGeV=pIU/gigaelectronvolt;
804#ifdef pdebug
805  G4cout<<"-->G4QuasiFreeR::GetElTot: P="<<pIU<<",pPDG="<<hPDG<<",Z="<<Z<<",N="<<N<<G4endl;
806#endif
807  if(Z<1 && N<1)
808  {
809    G4cout<<"-Warning-G4QuasiFreeRatio::GetElTot:Z="<<Z<<",N="<<N<<", return zero"<<G4endl;
810    return std::make_pair(0.,0.);
811  }
812  std::pair<G4double,G4double> hp=FetchElTot(pGeV, hPDG, true);
813  std::pair<G4double,G4double> hn=FetchElTot(pGeV, hPDG, false);
814#ifdef pdebug
815  G4cout<<"-OUT->G4QFRat::GetElTot: hp("<<hp.first<<","<<hp.second<<"), hn("<<hn.first<<","
816        <<hn.second<<")"<<G4endl;
817#endif
818  G4double A=(Z+N)/millibarn;                // To make the result in independent units(IU)
819  return std::make_pair((Z*hp.first+N*hn.first)/A,(Z*hp.second+N*hn.second)/A);
820} // End of GetElTot
821
822// (Mean Elastic and Mean Total) Cross-Sections (mb) for PDG+(Z,N) at P=p[GeV/c]
823std::pair<G4double,G4double> G4QuasiFreeRatios::GetChExFactor(G4double pIU, G4int hPDG,
824                                                              G4int Z, G4int N)
825{
826  G4double pGeV=pIU/gigaelectronvolt;
827  G4double resP=0.;
828  G4double resN=0.;
829  if(Z<1 && N<1)
830  {
831    G4cout<<"-Warning-G4QuasiFreeRatio::GetChExF:Z="<<Z<<",N="<<N<<", return zero"<<G4endl;
832    return std::make_pair(resP,resN);
833  }
834  G4double A=Z+N;
835  G4double pf=0.;                              // Possibility to interact with a proton
836  G4double nf=0.;                              // Possibility to interact with a neutron
837  if   (hPDG==-211||hPDG==-321||hPDG==3112||hPDG==3212||hPDG==3312) pf=Z/(A+N);
838  else if(hPDG==211||hPDG==321||hPDG==3222||hPDG==3212||hPDG==3322) nf=N/(A+Z);
839  else if(hPDG==-311||hPDG==311||hPDG==130||hPDG==310)
840  {
841    G4double dA=A+A;
842    pf=Z/(dA+N+N);
843    nf=N/(dA+Z+Z);
844  }
845  G4double mult=1.;  // Factor of increasing multiplicity ( ? @@)
846  if(pGeV>.5)
847  {
848    mult=1./(1.+std::log(pGeV+pGeV))/pGeV;
849    if(mult>1.) mult=1.;
850  }
851  if(pf)
852  {
853    std::pair<G4double,G4double> hp=FetchElTot(pGeV, hPDG, true);
854    resP=pf*(hp.second/hp.first-1.)*mult;
855  }
856  if(nf)
857  {
858    std::pair<G4double,G4double> hn=FetchElTot(pGeV, hPDG, false);
859    resN=nf*(hn.second/hn.first-1.)*mult;
860  }
861  return std::make_pair(resP,resN);
862} // End of GetChExFactor
863
864// scatter (pPDG,p4M) on a virtual nucleon (NPDG,N4M), result: final pair(newN4M,newp4M)
865// if(newN4M.e()==0.) - below threshold, XS=0, no scattering of the progectile happened
866std::pair<G4LorentzVector,G4LorentzVector> G4QuasiFreeRatios::Scatter(G4int NPDG,
867                                     G4LorentzVector N4M, G4int pPDG, G4LorentzVector p4M)
868{
869  static const G4double mNeut= G4QPDGCode(2112).GetMass();
870  static const G4double mProt= G4QPDGCode(2212).GetMass();
871  static const G4double mDeut= G4QPDGCode(2112).GetNuclMass(1,1,0);// Mass of deuteron
872  static const G4double mTrit= G4QPDGCode(2112).GetNuclMass(1,2,0);// Mass of tritium
873  static const G4double mHel3= G4QPDGCode(2112).GetNuclMass(2,1,0);// Mass of Helium3
874  static const G4double mAlph= G4QPDGCode(2112).GetNuclMass(2,2,0);// Mass of alpha
875  G4LorentzVector pr4M=p4M/megaelectronvolt;   // Convert 4-momenta in MeV (keep p4M)
876  N4M/=megaelectronvolt;
877  G4LorentzVector tot4M=N4M+p4M;
878#ifdef ppdebug
879  G4cerr<<"->G4QFR::Scat:p4M="<<pr4M<<",N4M="<<N4M<<",t4M="<<tot4M<<",NPDG="<<NPDG<<G4endl;
880#endif
881  G4double mT=mNeut;
882  G4int Z=0;
883  G4int N=1;
884  if(NPDG==2212||NPDG==90001000)
885  {
886    mT=mProt;
887    Z=1;
888    N=0;
889  }
890  else if(NPDG==90001001)
891  {
892    mT=mDeut;
893    Z=1;
894    N=1;
895  }
896  else if(NPDG==90002001)
897  {
898    mT=mHel3;
899    Z=2;
900    N=1;
901  }
902  else if(NPDG==90001002)
903  {
904    mT=mTrit;
905    Z=1;
906    N=2;
907  }
908  else if(NPDG==90002002)
909  {
910    mT=mAlph;
911    Z=2;
912    N=2;
913  }
914  else if(NPDG!=2112&&NPDG!=90000001)
915  {
916    G4cout<<"Error:G4QuasiFreeRatios::Scatter:NPDG="<<NPDG<<" is not 2212 or 2112"<<G4endl;
917    G4Exception("G4QuasiFreeRatios::Scatter:","21",FatalException,"CHIPScomplain");
918    //return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M);// Use this if not exception
919  }
920  G4double mT2=mT*mT;
921  G4double mP2=pr4M.m2();
922  G4double E=(tot4M.m2()-mT2-mP2)/(mT+mT);
923#ifdef pdebug
924  G4cerr<<"G4QFR::Scat:qM="<<mT<<",qM2="<<mT2<<",pM2="<<mP2<<",totM2="<<tot4M.m2()<<G4endl;
925#endif
926  G4double E2=E*E;
927  if(E<0. || E2<mP2)
928  {
929#ifdef ppdebug
930    G4cerr<<"-Warning-G4QFR::Scat:*Negative Energy*E="<<E<<",E2="<<E2<<"<M2="<<mP2<<G4endl;
931#endif
932    return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); // Do Nothing Action
933  }
934  G4double P=std::sqrt(E2-mP2);                   // Momentum in pseudo laboratory system
935  G4VQCrossSection* PCSmanager=G4QProtonElasticCrossSection::GetPointer();
936  G4VQCrossSection* NCSmanager=G4QNeutronElasticCrossSection::GetPointer();
937#ifdef ppdebug
938  G4cout<<"G4QFR::Scatter: Before XS, P="<<P<<", Z="<<Z<<", N="<<N<<", PDG="<<pPDG<<G4endl;
939#endif
940  // @@ Temporary NN t-dependence for all hadrons
941  if(pPDG>3400 || pPDG<-3400) G4cout<<"-Warning-G4QElast::Scatter: pPDG="<<pPDG<<G4endl;
942  G4int PDG=2212;                                                // *TMP* instead of pPDG
943  if(pPDG==2112||pPDG==-211||pPDG==-321) PDG=2112;               // *TMP* instead of pPDG
944  if(!Z && N==1)                 // Change for Quasi-Elastic on neutron
945  {
946    Z=1;
947    N=0;
948    if     (PDG==2212) PDG=2112;
949    else if(PDG==2112) PDG=2212;
950  }
951  G4double xSec=0.;                        // Prototype of Recalculated Cross Section *TMP*
952  if(PDG==2212) xSec=PCSmanager->GetCrossSection(false, P, Z, N, PDG); // P CrossSect *TMP*
953  else          xSec=NCSmanager->GetCrossSection(false, P, Z, N, PDG); // N CrossSect *TMP*
954#ifdef ppdebug
955  G4cout<<"G4QElast::Scatter:pPDG="<<pPDG<<",P="<<P<<",CS="<<xSec/millibarn<<G4endl;
956#endif
957#ifdef nandebug
958  if(xSec>0. || xSec<0. || xSec==0);
959  else  G4cout<<"*Warning*G4QElast::Scatter: xSec="<<xSec/millibarn<<G4endl;
960#endif
961  // @@ check a possibility to separate p, n, or alpha (!)
962  if(xSec <= 0.)                                    // The cross-section iz 0 -> Do Nothing
963  {
964#ifdef ppdebug
965    G4cerr<<"-Warning-G4QFR::Scat:**Zero XS**PDG="<<pPDG<<",NPDG="<<NPDG<<",P="<<P<<G4endl;
966#endif
967    return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); //Do Nothing Action
968  }
969  G4double mint=0.;                        // Prototype of functional rand -t (MeV^2) *TMP*
970  if(PDG==2212) mint=PCSmanager->GetExchangeT(Z,N,PDG);// P functional rand -t(MeV^2) *TMP*
971  else          mint=NCSmanager->GetExchangeT(Z,N,PDG);// N functional rand -t(MeV^2) *TMP*
972  G4double maxt=0.;                                    // Prototype of max possible -t
973  if(PDG==2212) maxt=PCSmanager->GetHMaxT();           // max possible -t
974  else          maxt=NCSmanager->GetHMaxT();           // max possible -t
975#ifdef ppdebug
976  G4cout<<"G4QFR::Scat:PDG="<<PDG<<",P="<<P<<",X="<<xSec<<",-t="<<mint<<"<"<<maxt<<", Z="
977        <<Z<<",N="<<N<<G4endl;
978#endif
979#ifdef nandebug
980  if(mint>-.0000001);
981  else  G4cout<<"*Warning*G4QFR::Scat: -t="<<mint<<G4endl;
982#endif
983  G4double cost=1.-(mint+mint)/maxt; // cos(theta) in CMS
984#ifdef ppdebug
985  G4cout<<"G4QFR::Scat:-t="<<mint<<"<"<<maxt<<", cost="<<cost<<", Z="<<Z<<",N="<<N<<G4endl;
986#endif
987  if(cost>1. || cost<-1. || !(cost>-1. || cost<=1.))
988  {
989    if     (cost>1.)  cost=1.;
990    else if(cost<-1.) cost=-1.;
991    else
992    {
993      G4double tm=0.;
994      if(PDG==2212) tm=PCSmanager->GetHMaxT();
995      else          tm=NCSmanager->GetHMaxT();
996      G4cerr<<"G4QuasiFreeRatio::Scat:*NAN* cost="<<cost<<",-t="<<mint<<",tm="<<tm<<G4endl;
997      return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); // Do Nothing Action
998    }
999  }
1000  G4LorentzVector reco4M=G4LorentzVector(0.,0.,0.,mT);      // 4mom of the recoil nucleon
1001  G4LorentzVector dir4M=tot4M-G4LorentzVector(0.,0.,0.,(tot4M.e()-mT)*.01);
1002  if(!G4QHadron(tot4M).RelDecayIn2(pr4M, reco4M, dir4M, cost, cost))
1003  {
1004    G4cerr<<"G4QFR::Scat:t="<<tot4M<<tot4M.m()<<",mT="<<mT<<",mP="<<std::sqrt(mP2)<<G4endl;
1005    //G4Exception("G4QFR::Scat:","009",FatalException,"Decay of ElasticComp");
1006    return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); // Do Nothing Action
1007  }
1008#ifdef ppdebug
1009  G4cout<<"G4QFR::Scat:p4M="<<pr4M<<"+r4M="<<reco4M<<",dr="<<dir4M<<",t4M="<<tot4M<<G4endl;
1010#endif
1011  return std::make_pair(reco4M*megaelectronvolt,pr4M*megaelectronvolt); // Result
1012} // End of Scatter
1013
1014// scatter (pPDG,p4M) on a virtual nucleon (NPDG,N4M), result: final pair(newN4M,newp4M)
1015// if(newN4M.e()==0.) - below threshold, XS=0, no scattering of the progectile happened
1016// User should himself change the charge (PDG) (e.g. pn->np, pi+n->pi0p, pi-p->pi0n etc.)
1017std::pair<G4LorentzVector,G4LorentzVector> G4QuasiFreeRatios::ChExer(G4int NPDG,
1018                                     G4LorentzVector N4M, G4int pPDG, G4LorentzVector p4M)
1019{
1020  static const G4double mNeut= G4QPDGCode(2112).GetMass();
1021  static const G4double mProt= G4QPDGCode(2212).GetMass();
1022  G4LorentzVector pr4M=p4M/megaelectronvolt;          // Convert 4-momenta in MeV(keep p4M)
1023  N4M/=megaelectronvolt;
1024  G4LorentzVector tot4M=N4M+p4M;
1025  G4int Z=0;
1026  G4int N=1;
1027  G4int RPDG=2212;                                     // PDG of the recoil nucleon
1028  G4int sPDG=0;                                        // PDG code of the scattered hadron
1029  G4double mS=0.;                                      // proto of mass of scattered hadron
1030  G4double mT=mProt;                                   // mass of the recoil nucleon
1031  if(NPDG==2212)
1032  {
1033    mT=mNeut;
1034    Z=1;
1035    N=0;
1036    RPDG=2112;                                         // PDG of the recoil nucleon
1037    if(pPDG==-211) sPDG=111;                           // pi+    -> pi0
1038    else if(pPDG==-321)
1039    {
1040      sPDG=310;                                        // K+     -> K0S
1041      if(G4UniformRand()>.5) sPDG=130;                 // K+     -> K0L
1042    }
1043    else if(pPDG==-311||pPDG==311||pPDG==130||pPDG==310) sPDG=321;  // K0     -> K+ (?)
1044    else if(pPDG==3112) sPDG=3212;                     // Sigma- -> Sigma0
1045    else if(pPDG==3212) sPDG=3222;                     // Sigma0 -> Sigma+
1046    else if(pPDG==3312) sPDG=3322;                     // Xi-    -> Xi0
1047  }
1048  else if(NPDG==2112) // Default
1049  {
1050    if(pPDG==211)  sPDG=111;                           // pi+    -> pi0
1051    else if(pPDG==321)
1052    {
1053      sPDG=310;                                        // K+     -> K0S
1054      if(G4UniformRand()>.5) sPDG=130;                 // K+     -> K0L
1055    }
1056    else if(pPDG==-311||pPDG==311||pPDG==130||pPDG==310) sPDG=-321; // K0     -> K- (?)
1057    else if(pPDG==3222) sPDG=3212;                     // Sigma+ -> Sigma0
1058    else if(pPDG==3212) sPDG=3112;                     // Sigma0 -> Sigma-
1059    else if(pPDG==3322) sPDG=3312;                     // Xi0    -> Xi-
1060  }
1061  else
1062  {
1063    G4cout<<"Error:G4QuasiFreeRatios::ChExer: NPDG="<<NPDG<<" is not 2212 or 2112"<<G4endl;
1064    G4Exception("G4QuasiFreeRatios::ChExer:","21",FatalException,"CHIPS complain");
1065    //return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M);// Use this if not exception
1066  }
1067  if(sPDG) mS=mNeut;
1068  else
1069  {
1070    G4cout<<"Error:G4QuasiFreeRatios::ChExer: BAD pPDG="<<pPDG<<", NPDG="<<NPDG<<G4endl;
1071    G4Exception("G4QuasiFreeRatios::ChExer:","21",FatalException,"CHIPS complain");
1072    //return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M);// Use this if not exception
1073  }
1074  G4double mT2=mT*mT;
1075  G4double mS2=mS*mS;
1076  G4double E=(tot4M.m2()-mT2-mS2)/(mT+mT);
1077  G4double E2=E*E;
1078  if(E<0. || E2<mS2)
1079  {
1080#ifdef pdebug
1081    G4cerr<<"-Warning-G4QFR::ChEx:*Negative Energy*E="<<E<<",E2="<<E2<<"<M2="<<mS2<<G4endl;
1082#endif
1083    return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); // Do Nothing Action
1084  }
1085  G4double P=std::sqrt(E2-mS2);                   // Momentum in pseudo laboratory system
1086  G4VQCrossSection* PCSmanager=G4QProtonElasticCrossSection::GetPointer();
1087  G4VQCrossSection* NCSmanager=G4QNeutronElasticCrossSection::GetPointer();
1088#ifdef debug
1089  G4cout<<"G4QFR::ChExer: Before XS, P="<<P<<", Z="<<Z<<", N="<<N<<", PDG="<<pPDG<<G4endl;
1090#endif
1091  // @@ Temporary NN t-dependence for all hadrons
1092  G4int PDG=2212;                                                // *TMP* instead of pPDG
1093  if(pPDG==2112||pPDG==-211||pPDG==-321) PDG=2112;               // *TMP* instead of pPDG
1094  if(!Z && N==1)                 // Change for Quasi-Elastic on neutron
1095  {
1096    Z=1;
1097    N=0;
1098    if     (PDG==2212) PDG=2112;
1099    else if(PDG==2112) PDG=2212;
1100  }
1101  G4double xSec=0.;                        // Prototype of Recalculated Cross Section *TMP*
1102  if(PDG==2212) xSec=PCSmanager->GetCrossSection(false, P, Z, N, PDG); // P CrossSect *TMP*
1103  else          xSec=NCSmanager->GetCrossSection(false, P, Z, N, PDG); // N CrossSect *TMP*
1104#ifdef debug
1105  G4cout<<"G4QElast::ChExer:pPDG="<<pPDG<<",P="<<P<<",CS="<<xSec/millibarn<<G4endl;
1106#endif
1107#ifdef nandebug
1108  if(xSec>0. || xSec<0. || xSec==0);
1109  else  G4cout<<"*Warning*G4QElast::ChExer: xSec="<<xSec/millibarn<<G4endl;
1110#endif
1111  // @@ check a possibility to separate p, n, or alpha (!)
1112  if(xSec <= 0.) // The cross-section iz 0 -> Do Nothing
1113  {
1114#ifdef pdebug
1115    G4cerr<<"-Warning-G4QFR::ChEx:**Zero XS**PDG="<<pPDG<<",NPDG="<<NPDG<<",P="<<P<<G4endl;
1116#endif
1117    return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); //Do Nothing Action
1118  }
1119  G4double mint=0.;                        // Prototype of functional rand -t (MeV^2) *TMP*
1120  if(PDG==2212) mint=PCSmanager->GetExchangeT(Z,N,PDG);// P functional rand -t(MeV^2) *TMP*
1121  else          mint=NCSmanager->GetExchangeT(Z,N,PDG);// N functional rand -t(MeV^2) *TMP*
1122#ifdef pdebug
1123  G4cout<<"G4QFR::ChEx:PDG="<<pPDG<<", P="<<P<<", CS="<<xSec<<", -t="<<mint<<G4endl;
1124#endif
1125#ifdef nandebug
1126  if(mint>-.0000001);
1127  else  G4cout<<"*Warning*G4QFR::ChExer: -t="<<mint<<G4endl;
1128#endif
1129  G4double maxt=0.;                                    // Prototype of max possible -t
1130  if(PDG==2212) maxt=PCSmanager->GetHMaxT();           // max possible -t
1131  else          maxt=NCSmanager->GetHMaxT();           // max possible -t
1132  G4double cost=1.-mint/maxt;                          // cos(theta) in CMS
1133#ifdef pdebug
1134  G4cout<<"G4QuasiFfreeRatio::ChExer: -t="<<mint<<", maxt="<<maxt<<", cost="<<cost<<G4endl;
1135#endif
1136  if(cost>1. || cost<-1. || !(cost>-1. || cost<=1.))
1137  {
1138    if     (cost>1.)  cost=1.;
1139    else if(cost<-1.) cost=-1.;
1140    else
1141    {
1142      G4cerr<<"G4QuasiFreeRatio::ChExer:*NAN* c="<<cost<<",t="<<mint<<",tm="<<maxt<<G4endl;
1143      return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); // Do Nothing Action
1144    }
1145  }
1146  G4LorentzVector reco4M=G4LorentzVector(0.,0.,0.,mT);      // 4mom of the recoil nucleon
1147  pr4M=G4LorentzVector(0.,0.,0.,mS);                        // 4mom of the scattered hadron
1148  G4LorentzVector dir4M=tot4M-G4LorentzVector(0.,0.,0.,(tot4M.e()-mT)*.01);
1149  if(!G4QHadron(tot4M).RelDecayIn2(pr4M, reco4M, dir4M, cost, cost))
1150  {
1151    G4cerr<<"G4QFR::ChEx:t="<<tot4M<<tot4M.m()<<",mT="<<mT<<",mP="<<mS<<G4endl;
1152    //G4Exception("G4QFR::ChExer:","009",FatalException,"Decay of ElasticComp");
1153    return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); // Do Nothing Action
1154  }
1155#ifdef debug
1156  G4cout<<"G4QFR::ChEx:p4M="<<p4M<<"+r4M="<<reco4M<<"="<<p4M+reco4M<<"="<<tot4M<<G4endl;
1157#endif
1158  return std::make_pair(reco4M*megaelectronvolt,pr4M*megaelectronvolt); // Result
1159} // End of ChExer
1160
1161// Calculate ChEx/El ratio (p is in independent units, (Z,N) is target, pPDG is projectile)
1162G4double G4QuasiFreeRatios::ChExElCoef(G4double p, G4int Z, G4int N, G4int pPDG) 
1163{
1164  p/=MeV;                                // Converted from independent units
1165  G4double A=Z+N;
1166  if(A<1.5) return 0.;
1167  G4double C=0.;
1168  if     (pPDG==2212) C=N/(A+Z);
1169  else if(pPDG==2112) C=Z/(A+N);
1170  else G4cout<<"*Warning*G4QCohChrgExchange::ChExElCoef: wrong PDG="<<pPDG<<G4endl;
1171  C*=C;                         // Coherent processes squares the amplitude
1172  // @@ This is true only for nucleons: other projectiles must be treated differently
1173  G4double sp=std::sqrt(p);
1174  G4double p2=p*p;           
1175  G4double p4=p2*p2;
1176  G4double dl1=std::log(p)-5.;
1177  G4double T=(6.75+.14*dl1*dl1+13./p)/(1.+.14/p4)+.6/(p4+.00013);
1178  G4double U=(6.25+8.33e-5/p4/p)*(p*sp+.34)/p2/p; 
1179  G4double R=U/T;
1180  return C*R*R;
1181}
Note: See TracBrowser for help on using the repository browser.