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