ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ParticleHPChannelList.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4ParticleHPChannelList.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 //
35 #include "G4Element.hh"
36 #include "G4HadFinalState.hh"
37 #include "G4HadProjectile.hh"
39 
41 
43  :theProjectile(projectile)
44  {
45  nChannels = n;
47  allChannelsCreated = false;
48  theInitCount = 0;
49  theElement = NULL;
50  }
51 
52 #include "G4Neutron.hh"
54  {
55  nChannels = 0;
56  theChannels = 0;
57  allChannelsCreated = false;
58  theInitCount = 0;
59  theElement = NULL;
61  }
62 
64  {
65  if(theChannels!=0)
66  {
67  for(G4int i=0;i<nChannels; i++)
68  {
69  delete theChannels[i];
70  }
71  delete [] theChannels;
72  }
73  }
74 
76  #include "G4ParticleHPManager.hh"
78  {
79  G4ParticleHPThermalBoost aThermalE;
80  G4int i, ii;
81  // decide on the isotope
82  G4int numberOfIsos(0);
83  for(ii=0; ii<nChannels; ii++)
84  {
85  numberOfIsos = theChannels[ii]->GetNiso();
86  if(numberOfIsos!=0) break;
87  }
88  G4double * running= new G4double [numberOfIsos];
89  running[0] = 0;
90  for(i=0;i<numberOfIsos; i++)
91  {
92  if(i!=0) running[i] = running[i-1];
93  for(ii=0; ii<nChannels; ii++)
94  {
95  if(theChannels[ii]->HasAnyData(i))
96  {
97  running[i] +=theChannels[ii]->GetWeightedXsec(aThermalE.GetThermalEnergy(aTrack,
98  theChannels[ii]->GetN(i),
99  theChannels[ii]->GetZ(i),
100  aTrack.GetMaterial()->GetTemperature()),
101  i);
102  }
103  }
104  }
105  G4int isotope=nChannels-1;
106  G4double random=G4UniformRand();
107  for(i=0;i<numberOfIsos; i++)
108  {
109  isotope = i;
110  //if(random<running[i]/running[numberOfIsos-1]) break;
111  if(running[numberOfIsos-1] != 0) if(random<running[i]/running[numberOfIsos-1]) break;
112  }
113  delete [] running;
114 
115  // decide on the channel
116  running = new G4double[nChannels];
117  running[0]=0;
118  G4int targA=-1; // For production of unChanged
119  G4int targZ=-1;
120  for(i=0; i<nChannels; i++)
121  {
122  if(i!=0) running[i] = running[i-1];
123  if(theChannels[i]->HasAnyData(isotope))
124  {
125  targA=(G4int)theChannels[i]->GetN(isotope); //Will be simply used the last valid value
126  targZ=(G4int)theChannels[i]->GetZ(isotope);
127  running[i] += theChannels[i]->GetFSCrossSection(aThermalE.GetThermalEnergy(aTrack,
128  targA,
129  targZ,
130  aTrack.GetMaterial()->GetTemperature()),
131  isotope);
132  targA=(G4int)theChannels[i]->GetN(isotope); //Will be simply used the last valid value
133  targZ=(G4int)theChannels[i]->GetZ(isotope);
134  // G4cout << " G4ParticleHPChannelList " << i << " targA " << targA << " targZ " << targZ << " running " << running[i] << G4endl;
135  }
136  }
137 
138  //TK120607
139  if ( running[nChannels-1] == 0 )
140  {
141  //This happened usually by the miss match between the cross section data and model
142  if ( targA == -1 && targZ == -1 ) {
143  throw G4HadronicException(__FILE__, __LINE__, "ParticleHP model encounter lethal discrepancy with cross section data");
144  }
145 
146  //TK121106
147  G4cout << "Warning from NeutronHP: could not find proper reaction channel. This may cause by inconsistency between cross section and model. Unchanged final states are returned." << G4endl;
148  unChanged.Clear();
149 
150  //For Ep Check create unchanged final state including rest target
151  G4ParticleDefinition* targ_pd = G4IonTable::GetIonTable()->GetIon ( targZ , targA , 0.0 );
152  G4DynamicParticle* targ_dp = new G4DynamicParticle( targ_pd , G4ThreeVector(1,0,0), 0.0 );
155  unChanged.AddSecondary(targ_dp);
156  //TK121106
159  delete [] running;
160  return &unChanged;
161  }
162  //TK120607
163 
164 
165  G4int lChan=0;
166  random=G4UniformRand();
167  for(i=0; i<nChannels; i++)
168  {
169  lChan = i;
170  if(running[nChannels-1] != 0) if(random<running[i]/running[nChannels-1]) break;
171  }
172  delete [] running;
173 #ifdef G4PHPDEBUG
174  if( std::getenv("G4ParticleHPDebug") ) G4cout << " G4ParticleHPChannelList SELECTED ISOTOPE " << isotope << " SELECTED CHANNEL " << lChan << G4endl;
175 #endif
176  return theChannels[lChan]->ApplyYourself(aTrack, isotope);
177  }
178 
179 void G4ParticleHPChannelList::Init(G4Element * anElement, const G4String & dirName, G4ParticleDefinition* projectile )
180  {
181  theDir = dirName;
182 // G4cout << theDir << G4endl;
183  theElement = anElement;
184 // G4cout << theElement << G4endl;
185  theProjectile = projectile;
186  }
187 
189  const G4String & aName )
190  {
191  if(!allChannelsCreated)
192  {
193  if(nChannels!=0)
194  {
195  G4ParticleHPChannel ** theBuffer = new G4ParticleHPChannel * [nChannels+1];
196  G4int i;
197  for(i=0; i<nChannels; i++)
198  {
199  theBuffer[i] = theChannels[i];
200  }
201  delete [] theChannels;
202  theChannels = theBuffer;
203  }
204  else
205  {
207  }
208  G4String name;
209  name = aName+"/";
212  // theChannels[nChannels]->SetProjectile( theProjectile );
213  nChannels++;
214  }
215 
216  //110527TKDB Unnessary codes, Detected by gcc4.6 compiler
217  //G4bool result;
218  //result = theChannels[theInitCount]->Register(theFS);
220  theInitCount++;
221  }
222 
224 
225  G4cout<<"================================================================"<<G4endl;
226  G4cout<<" Element: "<<theElement->GetName()<<G4endl;
227  G4cout<<" Number of channels: "<<nChannels<<G4endl;
228  G4cout<<" Projectile: "<<theProjectile->GetParticleName()<<G4endl;
229  G4cout<<" Directory name: "<<theDir<<G4endl;
230  for(int i=0;i<nChannels;i++){
232  G4cout<<"----------------------------------------------------------------"<<G4endl;
233  theChannels[i]->DumpInfo();
234  G4cout<<"----------------------------------------------------------------"<<G4endl;
235  }
236  }
237  G4cout<<"================================================================"<<G4endl;
238 
239 }