Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4QFreeScattering.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 // $Id$
28 //
29 //
30 // G4 Physics class: G4QFreeScattering for quasi-free scattering
31 // Created: M.V. Kossov, CERN/ITEP(Moscow), 26-OCT-11
32 // The last update: M.V. Kossov, CERN/ITEP (Moscow) 29-Oct-11
33 //
34 // ----------------------------------------------------------------------
35 // Short description: Provides quasi-free scattering on nuclear nucleons.
36 // ----------------------------------------------------------------------
37 
38 //#define debug
39 //#define pdebug
40 //#define ppdebug
41 //#define nandebug
42 
43 #include "G4QFreeScattering.hh"
44 #include "G4SystemOfUnits.hh"
45 
46 // initialisation of statics
47 std::vector<std::pair<G4double,G4double>*> G4QFreeScattering::vX; // ETPointers to LogTable
48 
49 G4QFreeScattering::G4QFreeScattering()
50 {
51 #ifdef pdebug
52  G4cout<<"***^^^*** G4QFreeScattering singletone is created ***^^^***"<<G4endl;
53 #endif
54 }
55 
56 G4QFreeScattering::~G4QFreeScattering()
57 {
58  std::vector<std::pair<G4double,G4double>*>::iterator pos2;
59  for(pos2=vX.begin(); pos2<vX.end(); pos2++) delete [] * pos2;
60  vX.clear();
61 }
62 
63 // Returns Pointer to the G4VQCrossSection class
64 G4QFreeScattering* G4QFreeScattering::GetPointer()
65 {
66  static G4QFreeScattering theQFS; // *** Static body of the Quasi-Free Scattering ***
67  return &theQFS;
68 }
69 
70 // Calculatio pair(hN_el,hN_tot) (mb): p in GeV/c, index(PDG,F) (see FetchElTot)
71 std::pair<G4double,G4double> G4QFreeScattering::CalcElTot(G4double p, G4int I)
72 {
73  // ---------> Each parameter set can have not more than nPoints=128 parameters
74  static const G4double lmi=3.5; // min of (lnP-lmi)^2 parabola
75  static const G4double pbe=.0557; // elastic (lnP-lmi)^2 parabola coefficient
76  static const G4double pbt=.3; // total (lnP-lmi)^2 parabola coefficient
77  static const G4double pmi=.1; // Below that fast LE calculation is made
78  static const G4double pma=1000.; // Above that fast HE calculation is made
79  G4double El=0.; // prototype of the elastic hN cross-section
80  G4double To=0.; // prototype of the total hN cross-section
81  if(p<=0.)
82  {
83  //G4cout<<"-Warning-G4QFreeScattering::CalcElTot: p="<<p<<" is 0 or negative"<<G4endl;
84  return std::make_pair(El,To);
85  }
86  if (!I) // pp/nn
87  {
88 #ifdef debug
89  G4cout<<"G4QFreeScatter::CalcElTot:I=0, p="<<p<<", pmi="<<pmi<<", pma="<<pma<<G4endl;
90 #endif
91  if(p<pmi)
92  {
93  G4double p2=p*p;
94  El=1./(.00012+p2*.2);
95  To=El;
96 #ifdef debug
97  G4cout<<"G4QFreeScatter::CalcElTot:I=0i, El="<<El<<", To="<<To<<", p2="<<p2<<G4endl;
98 #endif
99  }
100  else if(p>pma)
101  {
102  G4double lp=std::log(p)-lmi;
103  G4double lp2=lp*lp;
104  El=pbe*lp2+6.72;
105  To=pbt*lp2+38.2;
106 #ifdef debug
107  G4cout<<"G4QFreeScat::CalcElTot:I=0a, El="<<El<<", To="<<To<<", lp2="<<lp2<<G4endl;
108 #endif
109  }
110  else
111  {
112  G4double p2=p*p;
113  G4double LE=1./(.00012+p2*.2);
114  G4double lp=std::log(p)-lmi;
115  G4double lp2=lp*lp;
116  G4double rp2=1./p2;
117  El=LE+(pbe*lp2+6.72+32.6/p)/(1.+rp2/p);
118  To=LE+(pbt*lp2+38.2+52.7*rp2)/(1.+2.72*rp2*rp2);
119 #ifdef debug
120  G4cout<<"G4QFreeScat::CalcElTot:0,E="<<El<<",T="<<To<<",s="<<p2<<",l="<<lp2<<G4endl;
121 #endif
122  }
123  }
124  else if(I==1) // np/pn
125  {
126  if(p<pmi)
127  {
128  G4double p2=p*p;
129  El=1./(.00012+p2*(.051+.1*p2));
130  To=El;
131  }
132  else if(p>pma)
133  {
134  G4double lp=std::log(p)-lmi;
135  G4double lp2=lp*lp;
136  El=pbe*lp2+6.72;
137  To=pbt*lp2+38.2;
138  }
139  else
140  {
141  G4double p2=p*p;
142  G4double LE=1./(.00012+p2*(.051+.1*p2));
143  G4double lp=std::log(p)-lmi;
144  G4double lp2=lp*lp;
145  G4double rp2=1./p2;
146  El=LE+(pbe*lp2+6.72+30./p)/(1.+.49*rp2/p);
147  To=LE+(pbt*lp2+38.2)/(1.+.54*rp2*rp2);
148  }
149  }
150  else if(I==2) // pimp/pipn
151  {
152  G4double lp=std::log(p);
153  if(p<pmi)
154  {
155  G4double lr=lp+1.27;
156  El=1.53/(lr*lr+.0676);
157  To=El*3;
158  }
159  else if(p>pma)
160  {
161  G4double ld=lp-lmi;
162  G4double ld2=ld*ld;
163  G4double sp=std::sqrt(p);
164  El=pbe*ld2+2.4+7./sp;
165  To=pbt*ld2+22.3+12./sp;
166  }
167  else
168  {
169  G4double lr=lp+1.27; // p1
170  G4double LE=1.53/(lr*lr+.0676); // p2, p3
171  G4double ld=lp-lmi; // p4 (lmi=3.5)
172  G4double ld2=ld*ld;
173  G4double p2=p*p;
174  G4double p4=p2*p2;
175  G4double sp=std::sqrt(p);
176  G4double lm=lp+.36; // p5
177  G4double md=lm*lm+.04; // p6
178  G4double lh=lp-.017; // p7
179  G4double hd=lh*lh+.0025; // p8
180  El=LE+(pbe*ld2+2.4+7./sp)/(1.+.7/p4)+.6/md+.05/hd;//p9(pbe=.0557),p10,p11,p12,p13,p14
181  To=LE*3+(pbt*ld2+22.3+12./sp)/(1.+.4/p4)+1./md+.06/hd;
182  }
183  }
184  else if(I==3) // pipp/pimn
185  {
186  G4double lp=std::log(p);
187  if(p<pmi)
188  {
189  G4double lr=lp+1.27;
190  G4double lr2=lr*lr;
191  El=13./(lr2+lr2*lr2+.0676);
192  To=El;
193  }
194  else if(p>pma)
195  {
196  G4double ld=lp-lmi;
197  G4double ld2=ld*ld;
198  G4double sp=std::sqrt(p);
199  El=pbe*ld2+2.4+6./sp;
200  To=pbt*ld2+22.3+5./sp;
201  }
202  else
203  {
204  G4double lr=lp+1.27; // p1
205  G4double lr2=lr*lr;
206  G4double LE=13./(lr2+lr2*lr2+.0676); // p2, p3
207  G4double ld=lp-lmi; // p4 (lmi=3.5)
208  G4double ld2=ld*ld;
209  G4double p2=p*p;
210  G4double p4=p2*p2;
211  G4double sp=std::sqrt(p);
212  G4double lm=lp-.32; // p5
213  G4double md=lm*lm+.0576; // p6
214  El=LE+(pbe*ld2+2.4+6./sp)/(1.+3./p4)+.7/md; // p7(pbe=.0557), p8, p9, p10, p11
215  To=LE+(pbt*ld2+22.3+5./sp)/(1.+1./p4)+.8/md;
216  }
217  }
218  else if(I==4) // Kmp/Kmn/K0p/K0n
219  {
220 
221  if(p<pmi)
222  {
223  G4double psp=p*std::sqrt(p);
224  El=5.2/psp;
225  To=14./psp;
226  }
227  else if(p>pma)
228  {
229  G4double ld=std::log(p)-lmi;
230  G4double ld2=ld*ld;
231  El=pbe*ld2+2.23;
232  To=pbt*ld2+19.5;
233  }
234  else
235  {
236  G4double ld=std::log(p)-lmi;
237  G4double ld2=ld*ld;
238  G4double sp=std::sqrt(p);
239  G4double psp=p*sp;
240  G4double p2=p*p;
241  G4double p4=p2*p2;
242  G4double lm=p-.39;
243  G4double md=lm*lm+.000156;
244  G4double lh=p-1.;
245  G4double hd=lh*lh+.0156;
246  El=5.2/psp+(pbe*ld2+2.23)/(1.-.7/sp+.075/p4)+.004/md+.15/hd;
247  To=14./psp+(pbt*ld2+19.5)/(1.-.21/sp+.52/p4)+.006/md+.30/hd;
248  }
249  }
250  else if(I==5) // Kpp/Kpn/aKp/aKn
251  {
252  if(p<pmi)
253  {
254  G4double lr=p-.38;
255  G4double lm=p-1.;
256  G4double md=lm*lm+.372;
257  El=.7/(lr*lr+.0676)+2./md;
258  To=El+.6/md;
259  }
260  else if(p>pma)
261  {
262  G4double ld=std::log(p)-lmi;
263  G4double ld2=ld*ld;
264  El=pbe*ld2+2.23;
265  To=pbt*ld2+19.5;
266  }
267  else
268  {
269  G4double ld=std::log(p)-lmi;
270  G4double ld2=ld*ld;
271  G4double lr=p-.38;
272  G4double LE=.7/(lr*lr+.0676);
273  G4double sp=std::sqrt(p);
274  G4double p2=p*p;
275  G4double p4=p2*p2;
276  G4double lm=p-1.;
277  G4double md=lm*lm+.372;
278  El=LE+(pbe*ld2+2.23)/(1.-.7/sp+.1/p4)+2./md;
279  To=LE+(pbt*ld2+19.5)/(1.+.46/sp+1.6/p4)+2.6/md;
280  }
281  }
282  else if(I==6) // hyperon-N
283  {
284  if(p<pmi)
285  {
286  G4double p2=p*p;
287  El=1./(.002+p2*(.12+p2));
288  To=El;
289  }
290  else if(p>pma)
291  {
292  G4double lp=std::log(p)-lmi;
293  G4double lp2=lp*lp;
294  G4double sp=std::sqrt(p);
295  El=(pbe*lp2+6.72)/(1.+2./sp);
296  To=(pbt*lp2+38.2+900./sp)/(1.+27./sp);
297  }
298  else
299  {
300  G4double p2=p*p;
301  G4double LE=1./(.002+p2*(.12+p2));
302  G4double lp=std::log(p)-lmi;
303  G4double lp2=lp*lp;
304  G4double p4=p2*p2;
305  G4double sp=std::sqrt(p);
306  El=LE+(pbe*lp2+6.72+99./p2)/(1.+2./sp+2./p4);
307  To=LE+(pbt*lp2+38.2+900./sp)/(1.+27./sp+3./p4);
308  }
309  }
310  else if(I==7) // antibaryon-N
311  {
312  if(p>pma)
313  {
314  G4double lp=std::log(p)-lmi;
315  G4double lp2=lp*lp;
316  El=pbe*lp2+6.72;
317  To=pbt*lp2+38.2;
318  }
319  else
320  {
321  G4double ye=std::pow(p,1.25);
322  G4double yt=std::pow(p,.35);
323  G4double lp=std::log(p)-lmi;
324  G4double lp2=lp*lp;
325  El=80./(ye+1.)+pbe*lp2+6.72;
326  To=(80./yt+.3)/yt+pbt*lp2+38.2;
327  }
328  }
329  else
330  {
331  G4cout<<"*Error*G4QFreeScattering::CalcElTot:ind="<<I<<" is not defined (0-7)"<<G4endl;
332  G4Exception("G4QFreeScattering::CalcElTot:","23",FatalException,"CHIPScrash");
333  }
334  if(El>To) El=To;
335  return std::make_pair(El,To);
336 } // End of CalcElTot
337 
338 // For hadron PDG with momentum Mom (GeV/c) on N (p/n) calculate <sig_el,sig_tot> pair (mb)
339 // *** Only for external use. For simulation it is better to use GetElTot(...) [AMDB] ***
340 std::pair<G4double,G4double> G4QFreeScattering::GetElTotXS(G4double p, G4int PDG, G4bool F)
341 {
342  G4int ind=0; // Prototype of the reaction index
343  G4bool kfl=true; // Flag of K0/aK0 oscillation
344  G4bool kf=false;
345  if(PDG==90001000) PDG=2212;
346  if(PDG==90000001) PDG=2112;
347  if(PDG==91000000) PDG=3122;
348  if(PDG==130 || PDG==310)
349  {
350  kf=true;
351  if(G4UniformRand()>.5) kfl=false;
352  }
353  if ( (PDG == 2212 && F) || (PDG == 2112 && !F) ) ind=0; // pp/nn
354  else if ( (PDG == 2112 && F) || (PDG == 2212 && !F) ) ind=1; // np/pn
355  else if ( (PDG == -211 && F) || (PDG == 211 && !F) ) ind=2; // pimp/pipn
356  else if ( (PDG == 211 && F) || (PDG == -211 && !F) ) ind=3; // pipp/pimn
357  else if ( PDG == -321 || PDG == -311 || (kf && !kfl) ) ind=4; // KmN/K0N
358  else if ( PDG == 321 || PDG == 311 || (kf && kfl) ) ind=5; // KpN/aK0N
359  else if ( PDG > 3000 && PDG < 3335) ind=6; // @@ for all hyperons - take Lambda
360  else if ( PDG > -3335 && PDG < -2000) ind=7; // @@ for all anti-baryons (anti-p/anti-n)
361  else {
362  // VI changed Fatal to Warning, because this method cannot do better
363  // source of problem in classes which make this call
364  ind = 0;
365  G4cout<<"*Error*G4QFreeScattering::CalcElTotXS: PDG="<<PDG
366  <<", while it is defined only for p,n,hyperons,anti-baryons,pi,K/antiK"<<G4endl;
367  G4Exception("G4QFreeScattering::CalcElTotXS:","22",JustWarning,"CHIPS_crash");
368  }
369  return CalcElTot(p,ind); // This is a slow direct method, better use FetchElTot
370 }
371 
372 // Calculatio pair(hN_el,hN_tot)(mb): p in GeV/c, F=true -> N=proton, F=false -> N=neutron
373 std::pair<G4double,G4double> G4QFreeScattering::FetchElTot(G4double p, G4int PDG, G4bool F)
374 {
375  static const G4int nlp=300; // Number of steps in the S(lnp) logTable(5% step)
376  static const G4int mlp=nlp+1; // Number of elements in the S(lnp) logTable
377  static const G4double lpi=-5.; // The min ln(p) logTabEl(p=6.7 MeV/c - 22. TeV/c)
378  static const G4double lpa=10.; // The max ln(p) logTabEl(p=6.7 MeV/c - 22. TeV/c)
379  static const G4double mi=std::exp(lpi);// The min p of logTabEl(~ 6.7 MeV/c)
380  static const G4double ma=std::exp(lpa);// The max p of logTabEl(~ 22. TeV)
381  static const G4double dl=(lpa-lpi)/nlp;// Step of the logarithmic Table
382  static const G4double edl=std::exp(dl);// Multiplication step of the logarithmic Table
383  //static const G4double toler=.001; // Relative Tolarence defining "theSameMomentum"
384  static G4double lastP=0.; // The last momentum for which XS was calculated
385  static G4int lastH=0; // The last projPDG for which XS was calculated
386  static G4bool lastF=true; // The last nucleon for which XS was calculated
387  static std::pair<G4double,G4double> lastR(0.,0.); // The last result
388  // Local Associative Data Base:
389  static std::vector<G4int> vI; // Vector of index for which XS was calculated
390  static std::vector<G4double> vM; // Vector of rel max ln(p) initialized in LogTable
391  static std::vector<G4int> vK; // Vector of topBin number initialized in LogTable
392  // Last values of the Associative Data Base:
393  static G4int lastI=0; // The Last index for which XS was calculated
394  static G4double lastM=0.; // The Last rel max ln(p) initialized in LogTable
395  static G4int lastK=0; // The Last topBin number initialized in LogTable
396  static std::pair<G4double,G4double>* lastX=0; // The Last ETPointers to LogTable in heap
397  // LogTable is created only if necessary. The ratio R(s>8100 mb) = 0 for any nuclei
398  G4int nDB=vI.size(); // A number of hadrons already initialized in AMDB
399  if(PDG==90001000) PDG=2212;
400  if(PDG==90000001) PDG=2112;
401  if(PDG==91000000) PDG=3122;
402 #ifdef pdebug
403  G4cout<<"G4QFreeScatter::FetchElTot:p="<<p<<",PDG="<<PDG<<",F="<<F<<",nDB="<<nDB<<G4endl;
404 #endif
405  if(nDB && lastH==PDG && lastF==F && p>0. && p==lastP) return lastR;// VI don't use toler.
406  // if(nDB && lastH==PDG && lastF==F && p>0. && std::fabs(p-lastP)/p<toler) return lastR;
407  lastH=PDG;
408  lastF=F;
409  G4int ind=-1; // Prototipe of the index of the PDG/F combination
410  // i=0: pp(nn), i=1: np(pn), i=2: pimp(pipn), i=3: pipp(pimn), i=4: Kmp(Kmn,K0n,K0p),
411  // i=5: Kpp(Kpn,aK0n,aK0p), i=6: Hp(Hn), i=7: app(apn,ann,anp)
412  G4bool kfl=true; // Flag of K0/aK0 oscillation
413  G4bool kf=false;
414  if(PDG==130 || PDG==310)
415  {
416  kf=true;
417  if(G4UniformRand()>.5) kfl=false;
418  }
419  if ( (PDG == 2212 && F) || (PDG == 2112 && !F) ) ind=0; // pp/nn
420  else if ( (PDG == 2112 && F) || (PDG == 2212 && !F) ) ind=1; // np/pn
421  else if ( (PDG == -211 && F) || (PDG == 211 && !F) ) ind=2; // pimp/pipn
422  else if ( (PDG == 211 && F) || (PDG == -211 && !F) ) ind=3; // pipp/pimn
423  else if ( PDG == -321 || PDG == -311 || (kf && !kfl) ) ind=4; // KmN/K0N
424  else if ( PDG == 321 || PDG == 311 || (kf && kfl) ) ind=5; // KpN/aK0N
425  else if ( PDG > 3000 && PDG < 3335) ind=6; // @@ for all hyperons - take Lambda
426  else if ( PDG > -3335 && PDG < -2000) ind=7; // @@ for all anti-baryons (anti-p/anti-n)
427  else
428  {
429  // VI changed Fatal to Warning, because this method cannot do better
430  // source of problem in classes which make this call
431  ind = 0;
432  G4cout<<"*Error*G4QFreeScattering::FetchElTot: PDG="<<PDG
433  <<", while it is defined only for p,n,hyperons,anti-baryons,pi,K/antiK"<<G4endl;
434  G4Exception("G4QFreeScattering::FetchElTot:","22",JustWarning,"CHIPS problem");
435  }
436  if(nDB && lastI==ind && p>0. && p==lastP) return lastR; // VI do not use toler
437  // if(nDB && lastI==ind && p>0. && std::fabs(p-lastP)/p<toler) return lastR;
438  if(p<=mi || p>=ma) return CalcElTot(p,ind); // @@ Slow calculation ! (Warning?)
439  G4bool found=false;
440  G4int i=-1;
441  if(nDB) for (i=0; i<nDB; ++i) if(ind==vI[i]) // Sirch for this index in AMDB
442  {
443  found=true; // The index is found
444  break;
445  }
446  G4double lp=std::log(p);
447 #ifdef pdebug
448  G4cout<<"G4QFreeScat::FetchElTot: I="<<ind<<", i="<<i<<",fd="<<found<<",lp="<<lp<<G4endl;
449 #endif
450  if(!nDB || !found) // Create new line in the AMDB
451  {
452 #ifdef pdebug
453  G4cout<<"G4QFreeScattering::FetchElTot: NewX, ind="<<ind<<", nDB="<<nDB<<G4endl;
454 #endif
455  lastX = new std::pair<G4double,G4double>[mlp]; // Create logarithmic Table for ElTot
456  lastI = ind; // Remember the initialized inex
457  lastK = static_cast<int>((lp-lpi)/dl)+1; // MaxBin to be initialized in LogTaB
458  if(lastK>nlp)
459  {
460  lastK=nlp;
461  lastM=lpa-lpi;
462  }
463  else lastM = lastK*dl; // Calculate max initialized ln(p)-lpi for LogTab
464  G4double pv=mi;
465  for(G4int j=0; j<=lastK; ++j) // Calculate LogTab values
466  {
467  lastX[j]=CalcElTot(pv,ind);
468 #ifdef pdebug
469  G4cout<<"G4QFreeScat::FetchElTot:I,j="<<j<<",pv="<<pv<<",E="<<lastX[j].first<<",T="
470  <<lastX[j].second<<G4endl;
471 #endif
472  if(j!=lastK) pv*=edl;
473  }
474  i++; // Make a new record to AMDB and position on it
475  vI.push_back(lastI);
476  vM.push_back(lastM);
477  vK.push_back(lastK);
478  vX.push_back(lastX);
479  }
480  else // The A value was found in AMDB
481  {
482  lastI=vI[i];
483  lastM=vM[i];
484  lastK=vK[i];
485  lastX=vX[i];
486  G4int nextK=lastK+1;
487  G4double lpM=lastM+lpi;
488 #ifdef pdebug
489  G4cout<<"G4QFreeScatt::FetchElTo:M="<<lpM<<",l="<<lp<<",K="<<lastK<<",n="<<nlp<<G4endl;
490 #endif
491  if(lp>lpM && lastK<nlp) // LogTab must be updated
492  {
493  lastK = static_cast<int>((lp-lpi)/dl)+1; // MaxBin to be initialized in LogTab
494 #ifdef pdebug
495  G4cout<<"G4QFreeScat::FetET: K="<<lastK<<",lp="<<lp<<",li="<<lpi<<",dl="<<dl<<G4endl;
496 #endif
497  if(lastK>nlp)
498  {
499  lastK=nlp;
500  lastM=lpa-lpi;
501  }
502  else lastM = lastK*dl; // Calculate max initialized ln(p)-lpi for LogTab
503  G4double pv=std::exp(lpM); // momentum of the last calculated beam
504  for(G4int j=nextK; j<=lastK; ++j) // Calculate LogTab values
505  {
506  pv*=edl;
507  lastX[j]=CalcElTot(pv,ind);
508 #ifdef pdebug
509  G4cout<<"G4QFreeScat::FetchElTot: U:j="<<j<<",p="<<pv<<",E="<<lastX[j].first<<",T="
510  <<lastX[j].second<<G4endl;
511 #endif
512  }
513  } // End of LogTab update
514  if(lastK >= nextK) // The AMDB was apdated
515  {
516  vM[i]=lastM;
517  vK[i]=lastK;
518  }
519  }
520  // Now one can use tabeles to calculate the value
521  G4double dlp=lp-lpi; // Shifted log(p) value
522  G4int n=static_cast<int>(dlp/dl); // Low edge number of the bin
523  G4double d=dlp-n*dl; // Log shift
524  G4double e=lastX[n].first; // E-Base
525  lastR.first=e+d*(lastX[n+1].first-e)/dl; // E-Result
526  if(lastR.first<0.) lastR.first = 0.;
527  G4double t=lastX[n].second; // T-Base
528  lastR.second=t+d*(lastX[n+1].second-t)/dl; // T-Result
529  if(lastR.second<0.) lastR.second= 0.;
530 #ifdef pdebug
531  G4cout<<"=O=>G4QFreeScat::FetchElTot:1st="<<lastR.first<<", 2nd="<<lastR.second<<G4endl;
532 #endif
533  if(lastR.first>lastR.second) lastR.first = lastR.second;
534  return lastR;
535 } // End of FetchElTot
536 
537 // (Mean Elastic and Mean Total over (Z,N)) Cross-Sections (mb) for PDG+(Z,N) at P=p[GeV/c]
538 std::pair<G4double,G4double> G4QFreeScattering::GetElTotMean(G4double pIU, G4int hPDG,
539  G4int Z, G4int N)
540 {
541  G4double pGeV=pIU/gigaelectronvolt;
542  if(hPDG==90001000) hPDG=2212;
543  if(hPDG==90000001) hPDG=2112;
544  if(hPDG==91000000) hPDG=3122;
545 #ifdef pdebug
546  G4cout<<"->G4QFreeSc::GetElTotMean: P="<<pIU<<",pPDG="<<hPDG<<",Z="<<Z<<",N="<<N<<G4endl;
547 #endif
548  if(Z<1 && N<1)
549  {
550  G4cout<<"-Warning-G4QFreeScat::GetElTotMean: Z="<<Z<<",N="<<N<<", return zero"<<G4endl;
551  return std::make_pair(0.,0.);
552  }
553  std::pair<G4double,G4double> hp=FetchElTot(pGeV, hPDG, true);
554  std::pair<G4double,G4double> hn=FetchElTot(pGeV, hPDG, false);
555 #ifdef pdebug
556  G4cout<<"-OUT->G4QFreeScat::GetElTotMean: hp("<<hp.first<<","<<hp.second<<"), hn("
557  <<hn.first<<","<<hn.second<<")"<<G4endl;
558 #endif
559  G4double A=(Z+N)/millibarn; // To make the result in independent units(IU)
560  return std::make_pair((Z*hp.first+N*hn.first)/A,(Z*hp.second+N*hn.second)/A);
561 } // End of GetElTotMean
562 
563 // scatter (pPDG,p4M) on a virtual nucleon (NPDG,N4M), result: final pair(newN4M,newp4M)
564 // if(newN4M.e()==0.) - below threshold, XS=0, no scattering of the progectile happened
565 std::pair<G4LorentzVector,G4LorentzVector> G4QFreeScattering::Scatter(G4int NPDG,
566  G4LorentzVector N4M, G4int pPDG, G4LorentzVector p4M)
567 {
568  static const G4double mNeut= G4QPDGCode(2112).GetMass();
569  static const G4double mProt= G4QPDGCode(2212).GetMass();
570  p4M/=megaelectronvolt; // Convert 4-momenta in MeV (CHIPS units)
571  N4M/=megaelectronvolt;
572  G4LorentzVector tot4M=N4M+p4M;
573  if(pPDG==90001000) pPDG=2212;
574  if(pPDG==90000001) pPDG=2112;
575  if(pPDG==91000000) pPDG=3122;
576 #ifdef ppdebug
577  G4cout<<"->G4QFR::Scat:p4M="<<p4M<<",N4M="<<N4M<<",t4M="<<tot4M<<",NPDG="<<NPDG<<G4endl;
578 #endif
579  G4double mT=mNeut;
580 #ifdef ppdebug
581  G4int Z=0;
582  G4int N=1;
583 #endif
584  if(NPDG==2212 || NPDG==90001000)
585  {
586  mT=mProt;
587 #ifdef ppdebug
588  Z=1;
589  N=0;
590 #endif
591  }
592  else if(NPDG!=2112 && NPDG!=90000001)
593  {
594  G4cout<<"Error:G4QFreeScattering::Scatter:NPDG="<<NPDG<<" is not 2212 or 2112"<<G4endl;
595  G4Exception("G4QFreeScattering::Scatter:","21",FatalException,"CHIPScomplain");
596  //return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); // Use this if not exception
597  }
598  G4double mT2=mT*mT; // a squared mass of the free scattered nuclead cluster (FSNC)
599  G4double mP2=p4M.m2(); // a projectile squared mass
600  G4double E=(tot4M.m2()-mT2-mP2)/(mT+mT); // a projectile energy in the CMS of FSNC
601 #ifdef pdebug
602  G4cout<<"G4QFS::Scat:qM="<<mT<<",qM2="<<mT2<<",pM2="<<mP2<<",totM2="<<tot4M.m2()<<G4endl;
603 #endif
604  G4double E2=E*E;
605  if( E < 0. || E2 < mP2)
606  {
607 #ifdef ppdebug
608  G4cout<<"-Warning-G4QFS::Scat:*Negative Energy*E="<<E<<",E2="<<E2<<"<M2="<<mP2<<G4endl;
609 #endif
610  return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); // Do Nothing Action
611  }
612  G4double pP2=E2-mP2; // Squared Momentum in pseudo laboratory system (final particles)
613  G4double P=std::sqrt(pP2); // Momentum in pseudo laboratory system (final particles) ?
614 #ifdef ppdebug
615  G4cout<<"G4QFreeS::Scatter: Before XS, P="<<P<<",Z="<<Z<<",N="<<N<<",PDG="<<pPDG<<G4endl;
616 #endif
617  // @@ Temporary NN t-dependence for all hadrons
618  if(pPDG>3400 || pPDG<-3400) G4cout<<"-Warning-G4QFreeScat::Scatter: pPDG="<<pPDG<<G4endl;
619  // @@ check a possibility to separate p, n, or alpha (!)
620  G4double dmT=mT+mT;
621  G4double s_value=dmT*E2+mP2+mT2; // Mondelstam s (?)
622  G4double maxt=dmT*dmT*pP2/s_value; // max possible |t|
623  G4double mint=0.; // Prototype of functional rand -t (MeV^2)
624  if(P < 14.) mint=maxt*G4UniformRand(); // S-wave
625  else // Calculate slopes (no cash !)
626  {
627  G4double p4=pP2*pP2/1000000.; // @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
628  G4double theB1=(7.2+4.32/p4/(p4+12./P))/(1.+2.5/p4); // p4 in GeV, P in MeV
629  mint=-std::log(G4UniformRand())/theB1; // t-chan only
630  }
631 #ifdef ppdebug
632  G4cout<<"G4QFS::Scat:PDG="<<pPDG<<",P="<<P<<",-t="<<mint<<"<"<<maxt<<", Z="<<Z<<",N="<<N
633  <<G4endl;
634 #endif
635 #ifdef nandebug
636  if(mint>-.0000001);
637  else G4cout<<"*Warning*G4QFreeScattering::Scatter: -t="<<mint<<G4endl;
638 #endif
639  G4double cost=1.-(mint+mint)/maxt; // cos(theta) in CMS
640 #ifdef ppdebug
641  G4cout<<"G4QFS::Scat:-t="<<mint<<"<"<<maxt<<", cost="<<cost<<", Z="<<Z<<",N="<<N<<G4endl;
642 #endif
643  if(cost>1. || cost<-1. || !(cost>-1. || cost<=1.))
644  {
645  if (cost > 1.) cost = 1.;
646  else if(cost <-1.) cost =-1.;
647  else
648  {
649  G4cerr<<"G4QFreeScatter::Scat:*NAN* cost="<<cost<<",-t="<<mint<<",tm="<<maxt<<G4endl;
650  return std::make_pair(G4LorentzVector(0.,0.,0.,0.), p4M); // Do Nothing Action
651  }
652  }
653  G4LorentzVector reco4M=G4LorentzVector(0.,0.,0.,mT); // 4mom of the recoil nucleon
654  G4LorentzVector dir4M=tot4M-G4LorentzVector(0.,0.,0.,(tot4M.e()-mT)*.01);
655  if(!G4QHadron(tot4M).RelDecayIn2(p4M, reco4M, dir4M, cost, cost))
656  {
657  G4cerr<<"G4QFS::Scat:t="<<tot4M<<tot4M.m()<<",mT="<<mT<<",mP="<<std::sqrt(mP2)<<G4endl;
658  //G4Exception("G4QFS::Scat:","009",FatalException,"Decay of ElasticComp");
659  return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); // Do Nothing Action
660  }
661 #ifdef ppdebug
662  G4cout<<"G4QFS::Scat:p4M="<<p4M<<"+r4M="<<reco4M<<",dr="<<dir4M<<",t4M="<<tot4M<<G4endl;
663 #endif
664  return std::make_pair( reco4M*megaelectronvolt, p4M*megaelectronvolt ); // The Result
665 } // End of Scatter
666 
667 G4QHadronVector* G4QFreeScattering::InElF(G4int NPDG, G4LorentzVector N4M,
668  G4int pPDG, G4LorentzVector p4M)
669 {
670  static const G4double mPi0 = G4QPDGCode(111).GetMass();
671  static const G4double mPi = G4QPDGCode(211).GetMass(); // Pi+ (Pi-: -211)
672  static const G4double mK = G4QPDGCode(321).GetMass(); // K+ (K- : -321)
673  static const G4double mK0 = G4QPDGCode(311).GetMass(); // aK0 (K0 : -311)
674  static const G4double mNeut= G4QPDGCode(2112).GetMass();
675  static const G4double mProt= G4QPDGCode(2212).GetMass();
676  static const G4double mLamb= G4QPDGCode(3122).GetMass();
677  //static const G4double mSigZ= G4QPDGCode(3212).GetMass();
678  static const G4double mSigM= G4QPDGCode(3112).GetMass();
679  static const G4double mSigP= G4QPDGCode(3222).GetMass();
680  static const G4double mXiM = G4QPDGCode(3312).GetMass();
681  static const G4double mXiZ = G4QPDGCode(3322).GetMass();
682  //static const G4double mOmM = G4QPDGCode(3334).GetMass();
683  static const G4double third =1./3.;
684  static const G4double twothd =2./3.;
685 
686  p4M/=megaelectronvolt; // Convert 4-momenta in MeV (CHIPS units)
687  N4M/=megaelectronvolt;
688  G4LorentzVector c4M=N4M+p4M;
689  if(pPDG==90001000) pPDG=2212;
690  if(pPDG==90000001) pPDG=2112;
691  if(pPDG==91000000) pPDG=3122;
692 #ifdef ppdebug
693  G4cout<<"->G4QFR::InElF: p4M="<<p4M<<",N4M="<<N4M<<", c4M="<<c4M<<",NPDG="<<NPDG<<G4endl;
694 #endif
695  G4double mT=mNeut;
696  G4int Z=0;
697  //G4int N=1;
698  if(NPDG==2212 || NPDG==90001000)
699  {
700  NPDG=2212;
701  mT=mProt;
702  Z=1;
703  //N=0;
704  }
705  else if(NPDG!=2112 && NPDG!=90000001)
706  {
707  G4cout<<"Error:G4QFreeScattering::InElF:NPDG="<<NPDG<<" is not 2212 or 2112"<<G4endl;
708  G4Exception("G4QFreeScattering::InElF:","21",FatalException,"CHIPS complain");
709  }
710  else NPDG=2112;
711  G4double mC2=c4M.m2();
712  G4double mC=0.;
713  if( mC2 < -0.0000001 )
714  {
715 #ifdef ppdebug
716  G4cout<<"-Warning-G4QFS::InElF: Negative compoundMass ="<<mC2<<", c4M="<<c4M<<G4endl;
717 #endif
718  return 0; // Do Nothing Action
719  }
720  else if ( mC2 > 0.0000001) mC=std::sqrt(mC2);
721 #ifdef ppdebug
722  G4cout<<"G4QFS::InElF: mC ="<<mC<<G4endl;
723 #endif
724  G4double mP=0.;
725  G4double mP2=p4M.m2(); // a projectile squared mass
726  if( mP2 < -0.0000001 )
727  {
728 #ifdef ppdebug
729  G4cout<<"-Warning-G4QFS::InElF: Negative projectileMass ="<<mC2<<", c4M="<<c4M<<G4endl;
730 #endif
731  return 0; // Do Nothing Action
732  }
733  else if ( mP2 > 0.0000001) mP=std::sqrt(mP2);
734 #ifdef pdebug
735  G4cout<<"G4QFS::InElF:mT("<<mT<<")+mP("<<mP<<")+mPi0="<<mP+mT+mPi0<<"<? mC="<<mC<<G4endl;
736 #endif
737  if(pPDG > 3334 || pPDG < -321) // Annihilation/Inelastic of anti-barions not implemented
738  {
739  G4cout<<"-Warning-G4QFreeScat::InElF: pPDG="<<pPDG<<G4endl;
740  return 0; // Do Nothing Action
741  }
742  if(pPDG==130 || pPDG==310)
743  {
744  if(G4UniformRand()>.5) pPDG = 311; // K0
745  else pPDG =-311; // aK0
746  }
747  G4int mPDG=111; // Additional newtral pion by default
748  G4double mM=mPi0; // Default Pi0 mass
749  G4double r=G4UniformRand(); // A random number to split pi0/pi
750  if (pPDG == 2212) // p-N
751  {
752  if(Z) // p-p
753  {
754  if(r<twothd) // -> n + Pi+ + p
755  {
756  pPDG=2112; // n
757  mP=mNeut;
758  mPDG= 211; // Pi+
759  mM=mPi;
760  }
761  }
762  else // p-n
763  {
764  if(r<third) // -> n + Pi+ + n
765  {
766  pPDG=2112; // n
767  mP=mNeut;
768  mPDG= 211; // Pi+
769  mM=mPi;
770  }
771  else if(r<twothd) // -> p + Pi- + p
772  {
773  mPDG=-211; // Pi-
774  mM=mPi;
775  NPDG=2212; // p
776  mT=mProt;
777  }
778  }
779  }
780  else if (pPDG == 2112) // n-N
781  {
782  if(Z) // n-p
783  {
784  if(r<third) // -> n + Pi+ + n
785  {
786  mPDG= 211; // Pi+
787  mM=mPi;
788  NPDG=2112; // n
789  mT=mNeut;
790  }
791  else if(r<twothd) // -> p + Pi- + p
792  {
793  pPDG=2212; // p
794  mP=mProt;
795  mPDG=-211; // Pi-
796  mM=mPi;
797  }
798  }
799  else // n-n
800  {
801  if(r<twothd) // -> p + Pi- + n
802  {
803  pPDG=2212; // p
804  mP=mProt;
805  mPDG= -211; // Pi-
806  mM=mPi;
807  }
808  }
809  }
810  else if (pPDG == 211) // pip-N
811  {
812  if(Z) // pip-p
813  {
814  if(r<0.5) // -> Pi+ + Pi+ + n
815  {
816  mPDG= 211; // Pi+
817  mM=mPi;
818  NPDG=2112; // n
819  mT=mNeut;
820  }
821  }
822  else // pip-n
823  {
824  if(r<third) // -> Pi0 + Pi0 + p
825  {
826  pPDG= 111; // Pi0
827  mP=mPi0;
828  NPDG=2212; // p
829  mT=mProt;
830  }
831  else if(r<twothd) // -> Pi+ + Pi- + p
832  {
833  mPDG=-211; // Pi-
834  mM=mPi;
835  NPDG=2212; // p
836  mT=mProt;
837  }
838  }
839  }
840  else if (pPDG == -211) // pim-N
841  {
842  if(Z) // pim-p
843  {
844  if(r<third) // -> Pi0 + Pi0 + n
845  {
846  pPDG= 111; // Pi0
847  mP=mPi0;
848  NPDG=2112; // n
849  mT=mNeut;
850  }
851  else if(r<twothd) // -> Pi- + Pi+ + n
852  {
853  mPDG= 211; // Pi+
854  mM=mPi;
855  NPDG=2112; // n
856  mT=mNeut;
857  }
858  }
859  else // pim-n
860  {
861  if(r<0.5) // -> Pi- + Pi- + p
862  {
863  mPDG=-211; // Pi-
864  mM=mPi;
865  NPDG=2212; // p
866  mT=mProt;
867  }
868  }
869  }
870  else if (pPDG == -321) // Km-N
871  {
872  if(Z) // Km-p
873  {
874  if(r<0.25) // K0 + Pi0 + n
875  {
876  pPDG=-311; // K0
877  mP=mK0;
878  NPDG=2112; // n
879  mT=mNeut;
880  }
881  else if(r<0.5) // -> K- + Pi+ + n
882  {
883  mPDG= 211; // Pi+
884  mM=mPi;
885  NPDG=2112; // n
886  mT=mNeut;
887  }
888  else if(r<0.75) // -> K0 + Pi- + p
889  {
890  pPDG=-311; // K0
891  mP=mK0;
892  mPDG=-211; // Pi-
893  mM=mPi;
894  }
895  }
896  else // Km-n
897  {
898  if(r<third) // -> K- + Pi- + p
899  {
900  mPDG=-211; // Pi-
901  mM=mPi;
902  NPDG=2212; // p
903  mT=mProt;
904  }
905  else if(r<twothd) // -> K0 + Pi- + n
906  {
907  pPDG=-311; // K0
908  mP=mK0;
909  mPDG=-211; // Pi-
910  mM=mPi;
911  }
912  }
913  }
914  else if (pPDG == -311) // K0-N
915  {
916  if(Z) // K0-p
917  {
918  if(r<third) // K- + Pi+ + p
919  {
920  pPDG=-321; // K-
921  mP=mK;
922  NPDG=2212; // p
923  mT=mProt;
924  }
925  else if(r<twothd) // -> K0 + Pi+ + n
926  {
927  mPDG= 211; // Pi+
928  mM=mPi;
929  NPDG=2112; // n
930  mT=mNeut;
931  }
932  }
933  else // K0-n
934  {
935  if(r<0.25) // -> K- + Pi+ + n
936  {
937  pPDG=-321; // K-
938  mP=mK;
939  mPDG= 211; // Pi+
940  mM=mPi;
941  }
942  else if(r<0.5) // -> K- + Pi0 + p
943  {
944  pPDG=-321; // K-
945  mP=mK;
946  NPDG=2212; // p
947  mT=mProt;
948  }
949  else if(r<0.75) // -> K0 + Pi- + p
950  {
951  mPDG=-211; // Pi-
952  mM=mPi;
953  NPDG=2212; // p
954  mT=mProt;
955  }
956  }
957  }
958  else if (pPDG == 321) // Kp-N
959  {
960  if(Z) // Kp-p
961  {
962  if(r<third) // -> K+ + Pi+ + n
963  {
964  mPDG= 211; // Pi+
965  mM=mPi;
966  NPDG=2112; // n
967  mT=mNeut;
968  }
969  else if(r<twothd) // -> aK0 + Pi+ + p
970  {
971  pPDG= 311; // aK0
972  mP=mK0;
973  mPDG= 211; // Pi+
974  mM=mPi;
975  }
976  }
977  else // Kp-n
978  {
979  if(r<0.25) // -> aK0 + Pi+ + n
980  {
981  pPDG= 311; // aK0
982  mP=mK0;
983  mPDG= 211; // Pi+
984  mM=mPi;
985  }
986  else if(r<0.5) // -> Kp + Pi- + p
987  {
988  mPDG=-211; // Pi-
989  mM=mPi;
990  NPDG=2212; // p
991  mT=mProt;
992  }
993  else if(r<0.75) // -> aK0 + Pi0 + p
994  {
995  pPDG= 311; // aK0
996  mP=mK0;
997  NPDG=2212; // p
998  mT=mProt;
999  }
1000  }
1001  }
1002  else if (pPDG == 311) // aK0-N
1003  {
1004  if(Z) // aK0-p
1005  {
1006  if(r<0.25) // -> K+ + Pi- + p
1007  {
1008  pPDG= 321; // K+
1009  mP=mK;
1010  mPDG=-211; // Pi-
1011  mM=mPi;
1012  }
1013  else if(r<0.5) // -> aK0 + Pi+ + n
1014  {
1015  mPDG= 211; // Pi+
1016  mM=mPi;
1017  NPDG=2112; // n
1018  mT=mNeut;
1019  }
1020  else if(r<0.75) // -> K+ + Pi0 + n
1021  {
1022  pPDG= 321; // K+
1023  mP=mK;
1024  NPDG=2112; // n
1025  mT=mNeut;
1026  }
1027  }
1028  else // aK0-n
1029  {
1030  if(r<third) // -> aK0 + Pi- + p
1031  {
1032  mPDG=-211; // Pi-
1033  mM=mPi;
1034  NPDG=2212; // p
1035  mT=mProt;
1036  }
1037  else if(r<twothd) // -> K+ + Pi- + n
1038  {
1039  pPDG= 321; // K+
1040  mP=mK;
1041  mPDG=-211; // Pi-
1042  mM=mPi;
1043  }
1044  }
1045  }
1046  else if (pPDG == 3122 || pPDG== 3212) // Lambda/Sigma0-N
1047  {
1048  if(pPDG == 3212)
1049  {
1050  pPDG=3122;
1051  mP=mLamb;
1052  }
1053  if(Z) // L/S0-p
1054  {
1055  if(r<0.2) // -> SigP + Pi0 + n
1056  {
1057  pPDG=3222; // SigP
1058  mP=mSigP;
1059  NPDG=2112; // n
1060  mT=mNeut;
1061  }
1062  else if(r<0.4) // -> SigP + Pi- + p
1063  {
1064  pPDG=3222; // SigP
1065  mP=mSigP;
1066  mPDG=-211; // Pi-
1067  mM=mPi;
1068  }
1069  else if(r<0.6) // -> SigM + Pi+ + p
1070  {
1071  pPDG=3112; // SigM
1072  mP=mSigM;
1073  mPDG= 211; // Pi+
1074  mM=mPi;
1075  }
1076  else if(r<0.8) // -> Lamb + Pi+ + n
1077  {
1078  mPDG= 211; // Pi+
1079  mM=mPi;
1080  NPDG=2112; // n
1081  mT=mNeut;
1082  }
1083  }
1084  else // L/S0-n
1085  {
1086  if(r<0.2) // -> SigM + Pi0 + p
1087  {
1088  pPDG=3112; // SigM
1089  mP=mSigM;
1090  NPDG=2212; // p
1091  mT=mProt;
1092  }
1093  else if(r<0.4) // -> SigP + Pi- + n
1094  {
1095  pPDG=3222; // SigP
1096  mP=mSigP;
1097  mPDG=-211; // Pi-
1098  mM=mPi;
1099  }
1100  else if(r<0.6) // -> SigM + Pi+ + n
1101  {
1102  pPDG=3112; // SigM
1103  mP=mSigM;
1104  mPDG= 211; // Pi+
1105  mM=mPi;
1106  }
1107  else if(r<0.8) // -> Lamb + Pi- + p
1108  {
1109  mPDG=-211; // Pi-
1110  mM=mPi;
1111  NPDG=2212; // p
1112  mT=mProt;
1113  }
1114  }
1115  }
1116  else if (pPDG == 3112) // Sigma- -N
1117  {
1118  if(Z) // SigM-p
1119  {
1120  if(r<0.25) // -> Lamb + Pi0 + n
1121  {
1122  pPDG=3122; // Lamb
1123  mP=mLamb;
1124  NPDG=2112; // n
1125  mT=mNeut;
1126  }
1127  else if(r<0.5) // -> Lamb + Pi- + p
1128  {
1129  pPDG=3122; // Lamb
1130  mP=mLamb;
1131  mPDG=-211; // Pi-
1132  mM=mPi;
1133  }
1134  else if(r<0.75) // -> SigM + Pi+ + n
1135  {
1136  mPDG= 211; // Pi+
1137  mM=mPi;
1138  NPDG=2112; // n
1139  mT=mNeut;
1140  }
1141  }
1142  else // SigM-n
1143  {
1144  if(r<third) // -> Lamb + Pi- + n
1145  {
1146  pPDG=3122; // Lamb
1147  mP=mLamb;
1148  mPDG=-211; // Pi-
1149  mM=mPi;
1150  }
1151  else if(r<twothd) // -> SigM + Pi- + p
1152  {
1153  mPDG=-211; // Pi-
1154  mM=mPi;
1155  NPDG=2212; // p
1156  mT=mProt;
1157  }
1158  }
1159  }
1160  else if (pPDG == 3222) // Sigma+ -N
1161  {
1162  if(Z) // SigP-p
1163  {
1164  if(r<third) // -> Lamb + Pi+ + p
1165  {
1166  pPDG=3122; // Lamb
1167  mP=mLamb;
1168  mPDG= 211; // Pi+
1169  mM=mPi;
1170  }
1171  else if(r<twothd) // -> SigP + Pi+ + n
1172  {
1173  mPDG= 211; // Pi+
1174  mM=mPi;
1175  NPDG=2112; // n
1176  mT=mNeut;
1177  }
1178  }
1179  else // SigP-n
1180  {
1181  if(r<0.25) // -> Lamb + Pi0 + p
1182  {
1183  pPDG=3122; // Lamb
1184  mP=mLamb;
1185  NPDG=2212; // p
1186  mT=mProt;
1187  }
1188  else if(r<0.5) // -> Lamb + Pi+ + n
1189  {
1190  pPDG=3122; // Lamb
1191  mP=mLamb;
1192  mPDG= 211; // Pi+
1193  mM=mPi;
1194  }
1195  else if(r<0.75) // -> SigP + Pi- + p
1196  {
1197  mPDG=-211; // Pi-
1198  mM=mPi;
1199  NPDG=2212; // p
1200  mT=mProt;
1201  }
1202  }
1203  }
1204  else if (pPDG == 3312) // Xi- -N
1205  {
1206  if(Z) // XiM-p
1207  {
1208  if(r<0.25) // -> Xi0 + Pi0 + n
1209  {
1210  pPDG=3322; // Xi0
1211  mP=mXiZ;
1212  NPDG=2112; // n
1213  mT=mNeut;
1214  }
1215  else if(r<0.5) // -> Xi0 + Pi- + p
1216  {
1217  pPDG=3322; // Xi0
1218  mP=mXiZ;
1219  mPDG=-211; // Pi-
1220  mM=mPi;
1221  }
1222  else if(r<0.75) // -> XiM + Pi+ + n
1223  {
1224  mPDG= 211; // Pi+
1225  mM=mPi;
1226  NPDG=2112; // n
1227  mT=mNeut;
1228  }
1229  }
1230  else // XiM-n
1231  {
1232  if(r<third) // -> Xi0 + Pi- + n
1233  {
1234  pPDG=3322; // Xi0
1235  mP=mXiZ;
1236  mPDG=-211; // Pi-
1237  mM=mPi;
1238  }
1239  else if(r<twothd) // -> XiM + Pi- + p
1240  {
1241  mPDG=-211; // Pi-
1242  mM=mPi;
1243  NPDG=2212; // p
1244  mT=mProt;
1245  }
1246  }
1247  }
1248  else if (pPDG == 3322) // Xi0 -N
1249  {
1250  if(Z) // Xi0-p
1251  {
1252  if(r<third) // -> Xi- + Pi+ + p
1253  {
1254  pPDG=3312; // Xi-
1255  mP=mXiM;
1256  mPDG= 211; // Pi+
1257  mM=mPi;
1258  }
1259  else if(r<twothd) // -> Xi0 + Pi+ + n
1260  {
1261  mPDG= 211; // Pi+
1262  mM=mPi;
1263  NPDG=2112; // n
1264  mT=mNeut;
1265  }
1266  }
1267  else // Xi0-n
1268  {
1269  if(r<0.25) // -> Xi- + Pi0 + p
1270  {
1271  pPDG=3312; // Xi-
1272  mP=mXiM;
1273  NPDG=2212; // p
1274  mT=mProt;
1275  }
1276  else if(r<0.5) // -> Xi- + Pi+ + n
1277  {
1278  pPDG=3312; // Xi-
1279  mP=mXiM;
1280  mPDG= 211; // Pi+
1281  mM=mPi;
1282  }
1283  else if(r<0.75) // -> Xi0 + Pi- + p
1284  {
1285  mPDG=-211; // Pi-
1286  mM=mPi;
1287  NPDG=2212; // p
1288  mT=mProt;
1289  }
1290  }
1291  }
1292  else if (pPDG == 3334) // Om- -N
1293  {
1294  if(Z) // OmM-p
1295  {
1296  if(r<0.5) // -> Om- + Pi+ + n
1297  {
1298  mPDG= 211; // Pi+
1299  mM=mPi;
1300  NPDG=2112; // n
1301  mT=mNeut;
1302  }
1303  }
1304  else // OmM-n
1305  {
1306  if(r<0.5) // -> Om- + Pi- + p
1307  {
1308  mPDG=-211; // Pi-
1309  mM=mPi;
1310  NPDG=2212; // p
1311  mT=mProt;
1312  }
1313  }
1314  }
1315  else
1316  {
1317  G4cout<<"*Error*G4QFreeScattering::InElF: PDG="<<pPDG
1318  <<", while it is defined only for p,n,hyperons(not Omega),pi,K/antiK"<<G4endl;
1319  G4Exception("G4QFreeScattering::InElF:","22",FatalException,"CHIPS_crash");
1320  }
1321  if (mC-mP-mM-mT <-0.000001) return 0;
1322  G4QHadronVector* TripQH = new G4QHadronVector; // Proto of the Result
1323  G4LorentzVector m4M(0.,0.,0.,mM);
1324  if(mC-mP-mM-mT < 0.000001) // Equal share
1325  {
1326  p4M=(mP/mC)*c4M;
1327  m4M=(mM/mC)*c4M;
1328  N4M=(mT/mC)*c4M;
1329  }
1330  else
1331  {
1332  p4M=G4LorentzVector(0.,0.,0.,mP);
1333  N4M=G4LorentzVector(0.,0.,0.,mT);
1334  if(!G4QHadron(c4M).DecayIn3(p4M,m4M,N4M))
1335  {
1336  G4ExceptionDescription ed;
1337  ed << "DecayIn3, TotM=" << mC /* << " <? " << mT + mP + mM */ << G4endl;
1338  G4Exception("G4QFreeScattering::InElF()","HAD_CHPS_0027", FatalException, ed);
1339  }
1340  G4QHadron* h1 = new G4QHadron(pPDG, p4M);
1341  TripQH->push_back(h1); // (delete equivalent, responsibility of users)
1342 #ifdef debug
1343  G4cout << "G4QFreeScat::InElF: H1=" << pPDG << p4M << G4endl;
1344 #endif
1345  G4QHadron* h2 = new G4QHadron(mPDG, m4M);
1346  TripQH->push_back(h2); // (delete equivalent, responsibility of users)
1347 #ifdef debug
1348  G4cout << "G4QFreeScat::InElF: H2=" << mPDG << m4M << G4endl;
1349 #endif
1350  G4QHadron* h3 = new G4QHadron(NPDG, N4M);
1351  TripQH->push_back(h3); // (delete equivalent, responsibility of users)
1352 #ifdef debug
1353  G4cout << "G4QFreeScat::InElF: H3=" << NPDG << N4M << G4endl;
1354 #endif
1355  }
1356  return TripQH; // The Result
1357 } // End of Scatter
Generated on Sat May 25 2013 14:33:46 for Geant4 by   doxygen 1.8.4