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