ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ParticleHPCapture.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4ParticleHPCapture.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 // 070523 bug fix for G4FPE_DEBUG on by A. Howard ( and T. Koi)
31 //
32 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
33 //
34 #include "G4ParticleHPCapture.hh"
35 #include "G4ParticleHPManager.hh"
36 #include "G4SystemOfUnits.hh"
37 #include "G4ParticleHPCaptureFS.hh"
39 #include "G4ParticleTable.hh"
40 #include "G4IonTable.hh"
41 #include "G4Threading.hh"
42 
44  :G4HadronicInteraction("NeutronHPCapture")
45  ,theCapture(NULL)
46  ,numEle(0)
47  {
48  SetMinEnergy( 0.0 );
49  SetMaxEnergy( 20.*MeV );
50 /*
51 // G4cout << "Capture : start of construction!!!!!!!!"<<G4endl;
52  if(!std::getenv("G4NEUTRONHPDATA"))
53  throw G4HadronicException(__FILE__, __LINE__, "Please setenv G4NEUTRONHPDATA to point to the neutron cross-section files.");
54  dirName = std::getenv("G4NEUTRONHPDATA");
55  G4String tString = "/Capture";
56  dirName = dirName + tString;
57  numEle = G4Element::GetNumberOfElements();
58 // G4cout << "+++++++++++++++++++++++++++++++++++++++++++++++++"<<G4endl;
59 // G4cout <<"Disname="<<dirName<<" numEle="<<numEle<<G4endl;
60  //theCapture = new G4ParticleHPChannel[numEle];
61 // G4cout <<"G4ParticleHPChannel constructed"<<G4endl;
62  G4ParticleHPCaptureFS * theFS = new G4ParticleHPCaptureFS;
63  //for (G4int i=0; i<numEle; i++)
64  //{
65 // // G4cout << "initializing theCapture "<<i<<" "<< numEle<<G4endl;
66  // theCapture[i].Init((*(G4Element::GetElementTable()))[i], dirName);
67  // theCapture[i].Register(theFS);
68  //}
69  for ( G4int i = 0 ; i < numEle ; i++ )
70  {
71  theCapture.push_back( new G4ParticleHPChannel );
72  (*theCapture[i]).Init((*(G4Element::GetElementTable()))[i], dirName);
73  (*theCapture[i]).Register(theFS);
74  }
75  delete theFS;
76 // G4cout << "-------------------------------------------------"<<G4endl;
77 // G4cout << "Leaving G4ParticleHPCapture::G4ParticleHPCapture"<<G4endl;
78 */
79  }
80 
82  {
83  //delete [] theCapture;
84  //vector is shared, only master deletes
85  if ( ! G4Threading::IsWorkerThread() ) {
86  if ( theCapture != NULL ) {
87  for ( std::vector<G4ParticleHPChannel*>::iterator
88  ite = theCapture->begin() ; ite != theCapture->end() ; ite++ ) {
89  delete *ite;
90  }
91  theCapture->clear();
92  }
93  }
94  }
95 
98  {
99 
100  //if ( numEle < (G4int)G4Element::GetNumberOfElements() ) addChannelForNewElement();
101 
103  if(std::getenv("NeutronHPCapture")) G4cout <<" ####### G4ParticleHPCapture called"<<G4endl;
104  const G4Material * theMaterial = aTrack.GetMaterial();
105  G4int n = theMaterial->GetNumberOfElements();
106  G4int index = theMaterial->GetElement(0)->GetIndex();
107  if(n!=1)
108  {
109  G4double* xSec = new G4double[n];
110  G4double sum=0;
111  G4int i;
112  const G4double * NumAtomsPerVolume = theMaterial->GetVecNbOfAtomsPerVolume();
113  G4double rWeight;
114  G4ParticleHPThermalBoost aThermalE;
115  for (i=0; i<n; i++)
116  {
117  index = theMaterial->GetElement(i)->GetIndex();
118  rWeight = NumAtomsPerVolume[i];
119  //xSec[i] = theCapture[index].GetXsec(aThermalE.GetThermalEnergy(aTrack,
120  xSec[i] = ((*theCapture)[index])->GetXsec(aThermalE.GetThermalEnergy(aTrack,
121  theMaterial->GetElement(i),
122  theMaterial->GetTemperature()));
123  xSec[i] *= rWeight;
124  sum+=xSec[i];
125  }
126  G4double random = G4UniformRand();
127  G4double running = 0;
128  for (i=0; i<n; i++)
129  {
130  running += xSec[i];
131  index = theMaterial->GetElement(i)->GetIndex();
132  //if(random<=running/sum) break;
133  if( sum == 0 || random <= running/sum ) break;
134  }
135  if(i==n) i=std::max(0, n-1);
136  delete [] xSec;
137  }
138 
139  //return theCapture[index].ApplyYourself(aTrack);
140  //G4HadFinalState* result = theCapture[index].ApplyYourself(aTrack);
141  G4HadFinalState* result = ((*theCapture)[index])->ApplyYourself(aTrack);
142 
143  //Overwrite target parameters
144  aNucleus.SetParameters(G4ParticleHPManager::GetInstance()->GetReactionWhiteBoard()->GetTargA(),G4ParticleHPManager::GetInstance()->GetReactionWhiteBoard()->GetTargZ());
145  const G4Element* target_element = (*G4Element::GetElementTable())[index];
146  const G4Isotope* target_isotope=NULL;
147  G4int iele = target_element->GetNumberOfIsotopes();
148  for ( G4int j = 0 ; j != iele ; j++ ) {
149  target_isotope=target_element->GetIsotope( j );
150  if ( target_isotope->GetN() == G4ParticleHPManager::GetInstance()->GetReactionWhiteBoard()->GetTargA() ) break;
151  }
152  //G4cout << "Target Material of this reaction is " << theMaterial->GetName() << G4endl;
153  //G4cout << "Target Element of this reaction is " << target_element->GetName() << G4endl;
154  //G4cout << "Target Isotope of this reaction is " << target_isotope->GetName() << G4endl;
155  aNucleus.SetIsotope( target_isotope );
156 
158  return result;
159  }
160 
161 const std::pair<G4double, G4double> G4ParticleHPCapture::GetFatalEnergyCheckLevels() const
162 {
163  //return std::pair<G4double, G4double>(10*perCent,10*GeV);
164  return std::pair<G4double, G4double>(10*perCent,DBL_MAX);
165 }
166 
167 /*
168 void G4ParticleHPCapture::addChannelForNewElement()
169 {
170  G4ParticleHPCaptureFS* theFS = new G4ParticleHPCaptureFS;
171  for ( G4int i = numEle ; i < (G4int)G4Element::GetNumberOfElements() ; i++ )
172  {
173  G4cout << "G4ParticleHPCapture Prepairing Data for the new element of " << (*(G4Element::GetElementTable()))[i]->GetName() << G4endl;
174  theCapture.push_back( new G4ParticleHPChannel );
175  (*theCapture[i]).Init((*(G4Element::GetElementTable()))[i], dirName);
176  (*theCapture[i]).Register(theFS);
177  }
178  delete theFS;
179  numEle = (G4int)G4Element::GetNumberOfElements();
180 }
181 */
182 
184 {
186 }
188 {
190 }
191 
193 {
194 
196 
197  theCapture = hpmanager->GetCaptureFinalStates();
198 
199  if ( G4Threading::IsMasterThread() ) {
200 
201  if ( theCapture == NULL ) theCapture = new std::vector<G4ParticleHPChannel*>;
202 
203  if ( numEle == (G4int)G4Element::GetNumberOfElements() ) return;
204 
205  if ( theCapture->size() == G4Element::GetNumberOfElements() ) {
207  return;
208  }
209 
210  if ( !std::getenv("G4NEUTRONHPDATA") )
211  throw G4HadronicException(__FILE__, __LINE__, "Please setenv G4NEUTRONHPDATA to point to the neutron cross-section files.");
212  dirName = std::getenv("G4NEUTRONHPDATA");
213  G4String tString = "/Capture";
214  dirName = dirName + tString;
215 
217  for ( G4int i = numEle ; i < (G4int)G4Element::GetNumberOfElements() ; i++ )
218  {
219  theCapture->push_back( new G4ParticleHPChannel );
220  ((*theCapture)[i])->Init((*(G4Element::GetElementTable()))[i], dirName);
221  ((*theCapture)[i])->Register(theFS);
222  }
223  delete theFS;
225  }
227 }
228 
229 void G4ParticleHPCapture::ModelDescription(std::ostream& outFile) const
230 {
231  outFile << "High Precision model based on Evaluated Nuclear Data Files (ENDF) for radiative capture reaction of neutrons below 20MeV\n";
232 }