ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ParticleHPKallbachMannSyst.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4ParticleHPKallbachMannSyst.cc
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // neutron_hp -- source file
27 // J.P. Wellisch, Nov-1996
28 // A prototype of the low energy neutron transport model.
29 //
30 // 080801 Protect div0 error, when theCompundFraction is 1 by T. Koi
31 //
32 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
33 //
34 // June-2019 - E. Mendoza --> perform some corrections
35 
37 #include "G4SystemOfUnits.hh"
38 #include "Randomize.hh"
39 #include "G4Exp.hh"
40 #include "G4Log.hh"
41 #include "G4Pow.hh"
42 #include "G4HadronicException.hh"
43 
45 {
46  G4double result;
47 
48  G4double zero = GetKallbachZero(anEnergy);
49  if(zero>1) zero=1.;
50  if(zero<-1)zero=-1.;
51  G4double max = Kallbach(zero, anEnergy);
52  G4double upper = Kallbach(1., anEnergy);
53  G4double lower = Kallbach(-1., anEnergy);
54  if(upper>max) max=upper;
55  if(lower>max) max=lower;
56  G4double value, random;
57 
58  G4int icounter=0;
59  G4int icounter_max=1024;
60  do
61  {
62  icounter++;
63  if ( icounter > icounter_max ) {
64  G4cout << "Loop-counter exceeded the threshold value at " << __LINE__ << "th line of " << __FILE__ << "." << G4endl;
65  break;
66  }
67  result = 2.*G4UniformRand()-1;
68  value = Kallbach(result, anEnergy)/max;
69  random = G4UniformRand();
70  }
71  while(random>value); // Loop checking, 11.05.2015, T. Koi
72 
73  return result;
74 }
75 
77 {
78  // Kallbach-Mann systematics without normalization.
79  G4double result;
80  G4double theX = A(anEnergy)*cosTh;
81  result = 0.5*(G4Exp( theX)*(1+theCompoundFraction)
82  +G4Exp(-theX)*(1-theCompoundFraction));
83  return result;
84 }
85 
87 {
88  G4double result;
89  //delta 2.0e-16 in not good.
90  //delta 4.0e-16 is OK
91  //safety factor of 2
92  G4double delta = 8.0e-16;
93  if ( std::abs (theCompoundFraction - 1 ) < delta ) {
95  }
96  result = 0.5 * (1./A(anEnergy)) * G4Log((1-theCompoundFraction)/(1+theCompoundFraction));
97  return result;
98 }
99 
101 {
102  G4double result;
103  G4double C1 = 0.04/MeV;
104  G4double C2 = 1.8E-6/(MeV*MeV*MeV);
105  G4double C3 = 6.7E-7/(MeV*MeV*MeV*MeV);
106 
109  G4int Nc = Ac - theTargetZ-theProjectileZ;
110  G4int AA = theTargetA;
111  G4int ZA = theTargetZ;
112  G4double ea = epsa+SeparationEnergy(Ac, Nc, AA, ZA,theProjectileA,theProjectileZ);
113  G4double Et1 = 130*MeV;
114  G4double R1 = std::min(ea, Et1);
115  // theProductEnergy is still in CMS!!!
117  G4int AB = theResidualA;
118  G4int ZB = theResidualZ;
119  G4double eb = epsb+SeparationEnergy(Ac, Nc, AB, ZB,theProductA, theProductZ);
120  G4double X1 = R1*eb/ea;
121  G4double Et3 = 41*MeV;
122  G4double R3 = std::min(ea, Et3);
123  G4double X3 = R3*eb/ea;
124 
125  G4double Ma=1;
126  G4double mb=1;
127  if(theProjectileA==1 || (theProjectileZ==1 && theProjectileA==2)){Ma=1;}//neutron,proton,deuteron
128  else if(theProjectileA==4 && theProjectileZ==2){Ma=0;}//alpha
129  else if(theProjectileA==3 && (theProjectileZ==1 || theProjectileZ==2)){Ma=0.5;}//tritum,He3 : set intermediate value
130  else
131  {
132  throw G4HadronicException(__FILE__, __LINE__, "Severe error in the sampling of Kallbach-Mann Systematics");
133  }
134  if(theProductA==1 && theProductZ==0){mb=1./2.;}//neutron
135  else if(theProductA==4 && theProductZ==2){mb=2;}//alpha
136  else{mb=1;}
137 
138  result = C1*X1 + C2*G4Pow::GetInstance()->powN(X1, 3) + C3*Ma*mb*G4Pow::GetInstance()->powN(X3, 4);
139  return result;
140 
141 }
142 
144 {
145  G4double result;
146  G4int NA = AA-ZA;
147  G4int Zc = Ac-Nc;
148  result = 15.68*(Ac-AA);
149  result += -28.07*((Nc-Zc)*(Nc-Zc)/(G4double)Ac - (NA-ZA)*(NA-ZA)/(G4double)AA);
150  result += -18.56*(G4Pow::GetInstance()->A23(G4double(Ac)) - G4Pow::GetInstance()->A23(G4double(AA)));
151  result += 33.22*((Nc-Zc)*(Nc-Zc)/G4Pow::GetInstance()->powA(G4double(Ac), 4./3.) - (NA-ZA)*(NA-ZA)/G4Pow::GetInstance()->powA(G4double(AA), 4./3.));
152  result += -0.717*(Zc*Zc/G4Pow::GetInstance()->A13(G4double(Ac))-ZA*ZA/G4Pow::GetInstance()->A13(G4double(AA)));
153  result += 1.211*(Zc*Zc/(G4double)Ac-ZA*ZA/(G4double)AA);
154  G4double totalBinding(0);
155  if(Zbinding==0&&Abinding==1) totalBinding=0;
156  if(Zbinding==1&&Abinding==1) totalBinding=0;
157  if(Zbinding==1&&Abinding==2) totalBinding=2.224596;
158  if(Zbinding==1&&Abinding==3) totalBinding=8.481798;
159  if(Zbinding==2&&Abinding==3) totalBinding=7.718043;
160  if(Zbinding==2&&Abinding==4) totalBinding=28.29566;
161  result += -totalBinding;
162  result *= MeV;
163  return result;
164 }