ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4NeutronRadCapture.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4NeutronRadCapture.cc
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 //
28 // Physics model class G4NeutronRadCapture
29 // Created: 31 August 2009
30 // Author V.Ivanchenko
31 //
32 // Modified:
33 // 09.09.2010 V.Ivanchenko added usage of G4PhotonEvaporation
34 //
35 
36 #include "G4NeutronRadCapture.hh"
37 #include "G4SystemOfUnits.hh"
38 #include "G4ParticleDefinition.hh"
39 #include "G4Fragment.hh"
40 #include "G4FragmentVector.hh"
41 #include "G4NucleiProperties.hh"
42 #include "G4VEvaporationChannel.hh"
43 #include "G4PhotonEvaporation.hh"
44 #include "G4DynamicParticle.hh"
45 #include "G4ParticleTable.hh"
46 #include "G4IonTable.hh"
47 #include "G4Electron.hh"
48 #include "G4Deuteron.hh"
49 #include "G4Triton.hh"
50 #include "G4He3.hh"
51 #include "G4Alpha.hh"
52 #include "G4RandomDirection.hh"
53 #include "G4HadronicParameters.hh"
54 
56  : G4HadronicInteraction("nRadCapture"),
57  photonEvaporation(nullptr),lab4mom(0.,0.,0.,0.)
58 {
61  SetMinEnergy( 0.0*CLHEP::GeV );
63 
65  icID = -1;
66 
68 }
69 
71 {
72  delete photonEvaporation;
73 }
74 
76 {
77  if(photonEvaporation != nullptr) { return; }
78  G4DeexPrecoParameters* param =
81  icID = param->GetInternalConversionID();
82 
86 }
87 
89  const G4HadProjectile& aTrack, G4Nucleus& theNucleus)
90 {
93 
94  G4int A = theNucleus.GetA_asInt();
95  G4int Z = theNucleus.GetZ_asInt();
96 
97  G4double time = aTrack.GetGlobalTime();
98 
99  // Create initial state
101  lab4mom += aTrack.Get4Momentum();
102 
103  G4double M = lab4mom.mag();
104  ++A;
106  //G4cout << "Capture start: Z= " << Z << " A= " << A
107  // << " LabM= " << M << " Mcompound= " << mass << G4endl;
108 
109  // simplified method of 1 gamma emission
110  if(A <= 4) {
111 
113 
114  if(M - mass <= lowestEnergyLimit) {
115  return &theParticleChange;
116  }
117 
118  if (verboseLevel > 1) {
119  G4cout << "G4NeutronRadCapture::DoIt: Eini(MeV)="
120  << aTrack.GetKineticEnergy()/MeV << " Eexc(MeV)= "
121  << (M - mass)/MeV
122  << " Z= " << Z << " A= " << A << G4endl;
123  }
124  G4double e1 = (M - mass)*(M + mass)/(2*M);
125  G4LorentzVector lv2(e1*G4RandomDirection(),e1);
126  lv2.boost(bst);
127  G4HadSecondary* news =
129  news->SetTime(time);
131  delete news;
132 
133  const G4ParticleDefinition* theDef = 0;
134 
135  lab4mom -= lv2;
136  if (Z == 1 && A == 2) {theDef = G4Deuteron::Deuteron();}
137  else if (Z == 1 && A == 3) {theDef = G4Triton::Triton();}
138  else if (Z == 2 && A == 3) {theDef = G4He3::He3();}
139  else if (Z == 2 && A == 4) {theDef = G4Alpha::Alpha();}
140  else { theDef = theTableOfIons->GetIon(Z,A,0.0,noFloat,0); }
141 
142  if (verboseLevel > 1) {
143  G4cout << "Gamma 4-mom: " << lv2 << " "
144  << theDef->GetParticleName() << " " << lab4mom << G4endl;
145  }
146  if(theDef) {
147  news = new G4HadSecondary(new G4DynamicParticle(theDef, lab4mom));
148  news->SetTime(time);
150  delete news;
151  }
152 
153  // Use photon evaporation
154  } else {
155 
156  // protection against wrong kinematic
157  if(M < mass) {
158  G4double etot = std::max(mass, lab4mom.e());
159  G4double ptot = std::sqrt((etot - mass)*(etot + mass));
161  lab4mom.set(v.x()*ptot,v.y()*ptot,v.z()*ptot,etot);
162  }
163 
164  G4Fragment* aFragment = new G4Fragment(A, Z, lab4mom);
165 
166  if (verboseLevel > 1) {
167  G4cout << "G4NeutronRadCapture::ApplyYourself initial G4Fragmet:"
168  << G4endl;
169  G4cout << aFragment << G4endl;
170  }
171 
172  //
173  // Sample final state
174  //
176  if(!fv) { fv = new G4FragmentVector(); }
177  fv->push_back(aFragment);
178  size_t n = fv->size();
179 
180  if (verboseLevel > 1) {
181  G4cout << "G4NeutronRadCapture: " << n << " final particle icID= " << icID << G4endl;
182  }
183  for(size_t i=0; i<n; ++i) {
184 
185  G4Fragment* f = (*fv)[i];
186  G4double etot = f->GetMomentum().e();
187 
188  Z = f->GetZ_asInt();
189  A = f->GetA_asInt();
190 
191  const G4ParticleDefinition* theDef;
192  if(0 == Z && 0 == A) {theDef = f->GetParticleDefinition();}
193  else if (Z == 1 && A == 2) {theDef = G4Deuteron::Deuteron();}
194  else if (Z == 1 && A == 3) {theDef = G4Triton::Triton();}
195  else if (Z == 2 && A == 3) {theDef = G4He3::He3();}
196  else if (Z == 2 && A == 4) {theDef = G4Alpha::Alpha();}
197  else {
198  G4double eexc = f->GetExcitationEnergy();
199  if(eexc <= minExcitation) { eexc = 0.0; }
200  theDef = theTableOfIons->GetIon(Z, A, eexc, noFloat, 0);
201  /*
202  G4cout << "### NC Find ion Z= " << Z << " A= " << A
203  << " Eexc(MeV)= " << eexc/MeV << " "
204  << theDef << G4endl;
205  */
206  }
207  G4double ekin = std::max(0.0,etot - theDef->GetPDGMass());
208  if (verboseLevel > 1) {
209  G4cout << i << ". " << theDef->GetParticleName()
210  << " Ekin(MeV)= " << etot/MeV
211  << " p: " << f->GetMomentum().vect()
212  << G4endl;
213  }
214  G4HadSecondary* news = new G4HadSecondary(
215  new G4DynamicParticle(theDef,
216  f->GetMomentum().vect().unit(),
217  ekin));
218  G4double timeF = f->GetCreationTime();
219  if(timeF < 0.0) { timeF = 0.0; }
220  news->SetTime(time + timeF);
221  if(theDef == electron) { news->SetCreatorModelType(icID); }
223  delete news;
224  delete f;
225  }
226  delete fv;
227  }
228  //G4cout << "Capture done" << G4endl;
229  return &theParticleChange;
230 }
231