Geant4 10.7.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4ParticleHPInelastic.cc
Go to the documentation of this file.
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// this code implementation is the intellectual property of
27// neutron_hp -- source file
28// J.P. Wellisch, Nov-1996
29// A prototype of the low energy neutron transport model.
30//
31// By copying, distributing or modifying the Program (or any work
32// based on the Program) you indicate your acceptance of this statement,
33// and all its terms.
34//
35//
36// 070523 bug fix for G4FPE_DEBUG on by A. Howard (and T. Koi)
37// 081203 limit maximum trial for creating final states add protection for 1H isotope case by T. Koi
38//
39// P. Arce, June-2014 Conversion neutron_hp to particle_hp
40//
41#include "G4ParticleHPInelastic.hh"
42#include "G4SystemOfUnits.hh"
43#include "G4ParticleHPManager.hh"
44#include "G4HadronicParameters.hh"
45#include "G4Threading.hh"
46
47G4ParticleHPInelastic::G4ParticleHPInelastic(G4ParticleDefinition* projectile, const char* name )
48 :G4HadronicInteraction(name)
49 ,theInelastic(NULL)
50 ,numEle(0)
51 ,theProjectile(projectile)
52{
53 G4String baseName;
54 if ( std::getenv("G4PARTICLEHPDATA") ) {
55 baseName = std::getenv( "G4PARTICLEHPDATA" );
56 }
57 //const char* dataDirVariable;
58 G4String particleName;
59 if ( theProjectile == G4Neutron::Neutron() ) {
60 dataDirVariable = "G4NEUTRONHPDATA";
61 } else if( theProjectile == G4Proton::Proton() ) {
62 dataDirVariable = "G4PROTONHPDATA";
63 particleName = "Proton";
64 } else if( theProjectile == G4Deuteron::Deuteron() ) {
65 dataDirVariable = "G4DEUTERONHPDATA";
66 particleName = "Deuteron";
67 } else if( theProjectile == G4Triton::Triton() ) {
68 dataDirVariable = "G4TRITONHPDATA";
69 particleName = "Triton";
70 } else if( theProjectile == G4He3::He3() ) {
71 dataDirVariable = "G4HE3HPDATA";
72 particleName = "He3";
73 } else if( theProjectile == G4Alpha::Alpha() ) {
74 dataDirVariable = "G4ALPHAHPDATA";
75 particleName = "Alpha";
76 } else {
77 G4String message("G4ParticleHPInelastic may only be called for neutron, proton, deuteron, triton, He3 or alpha, while it is called for " + theProjectile->GetParticleName());
78 throw G4HadronicException(__FILE__, __LINE__,message.c_str());
79 }
80
81 SetMinEnergy( 0.0 );
82 SetMaxEnergy( 20.*MeV );
83
84 //G4cout << " entering G4ParticleHPInelastic constructor"<<G4endl;
85 if ( !std::getenv("G4PARTICLEHPDATA") && !std::getenv(dataDirVariable) ) {
86 G4String message("Please setenv G4PARTICLEHPDATA (recommended) or, at least setenv " +
87 G4String(dataDirVariable) + " to point to the " + theProjectile->GetParticleName() + " cross-section files." );
88 throw G4HadronicException(__FILE__, __LINE__,message.c_str());
89 }
90 if ( std::getenv(dataDirVariable) ) {
91 dirName = std::getenv(dataDirVariable);
92 } else {
93 dirName = baseName + "/" + particleName;
94 }
95
96 #ifdef G4VERBOSE
97 if ( G4HadronicParameters::Instance()->GetVerboseLevel() > 0 ) G4cout << dirName << G4endl;
98 #endif
99
100 G4String tString = "/Inelastic";
101 dirName = dirName + tString;
102 //numEle = G4Element::GetNumberOfElements();
103
104 #ifdef G4VERBOSE
105 if ( G4HadronicParameters::Instance()->GetVerboseLevel() > 0 )
106 G4cout << "@@@ G4ParticleHPInelastic instantiated for particle " << theProjectile->GetParticleName() << " data directory variable is " << dataDirVariable << " pointing to " << dirName << G4endl;
107 #endif
108
109/*
110 theInelastic = new G4ParticleHPChannelList[numEle];
111 for (G4int i=0; i<numEle; i++)
112 {
113 theInelastic[i].Init((*(G4Element::GetElementTable()))[i], dirName);
114 G4int itry = 0;
115 do
116 {
117 theInelastic[i].Register(&theNFS, "F01"); // has
118 theInelastic[i].Register(&theNXFS, "F02");
119 theInelastic[i].Register(&the2NDFS, "F03");
120 theInelastic[i].Register(&the2NFS, "F04"); // has, E Done
121 theInelastic[i].Register(&the3NFS, "F05"); // has, E Done
122 theInelastic[i].Register(&theNAFS, "F06");
123 theInelastic[i].Register(&theN3AFS, "F07");
124 theInelastic[i].Register(&the2NAFS, "F08");
125 theInelastic[i].Register(&the3NAFS, "F09");
126 theInelastic[i].Register(&theNPFS, "F10");
127 theInelastic[i].Register(&theN2AFS, "F11");
128 theInelastic[i].Register(&the2N2AFS, "F12");
129 theInelastic[i].Register(&theNDFS, "F13");
130 theInelastic[i].Register(&theNTFS, "F14");
131 theInelastic[i].Register(&theNHe3FS, "F15");
132 theInelastic[i].Register(&theND2AFS, "F16");
133 theInelastic[i].Register(&theNT2AFS, "F17");
134 theInelastic[i].Register(&the4NFS, "F18"); // has, E Done
135 theInelastic[i].Register(&the2NPFS, "F19");
136 theInelastic[i].Register(&the3NPFS, "F20");
137 theInelastic[i].Register(&theN2PFS, "F21");
138 theInelastic[i].Register(&theNPAFS, "F22");
139 theInelastic[i].Register(&thePFS, "F23");
140 theInelastic[i].Register(&theDFS, "F24");
141 theInelastic[i].Register(&theTFS, "F25");
142 theInelastic[i].Register(&theHe3FS, "F26");
143 theInelastic[i].Register(&theAFS, "F27");
144 theInelastic[i].Register(&the2AFS, "F28");
145 theInelastic[i].Register(&the3AFS, "F29");
146 theInelastic[i].Register(&the2PFS, "F30");
147 theInelastic[i].Register(&thePAFS, "F31");
148 theInelastic[i].Register(&theD2AFS, "F32");
149 theInelastic[i].Register(&theT2AFS, "F33");
150 theInelastic[i].Register(&thePDFS, "F34");
151 theInelastic[i].Register(&thePTFS, "F35");
152 theInelastic[i].Register(&theDAFS, "F36");
153 theInelastic[i].RestartRegistration();
154 itry++;
155 }
156 //while(!theInelastic[i].HasDataInAnyFinalState());
157 while( !theInelastic[i].HasDataInAnyFinalState() && itry < 6 );
158 // 6 is corresponding to the value(5) of G4ParticleHPChannel. TK
159
160 if ( itry == 6 )
161 {
162 // No Final State at all.
163 G4bool exceptional = false;
164 if ( (*(G4Element::GetElementTable()))[i]->GetNumberOfIsotopes() == 1 )
165 {
166 if ( (*(G4Element::GetElementTable()))[i]->GetIsotope( 0 )->GetZ() == 1 && (*(G4Element::GetElementTable()))[i]->GetIsotope( 0 )->GetN() == 1 ) exceptional = true; //1H
167 }
168 if ( !exceptional ) throw G4HadronicException(__FILE__, __LINE__, "Channel: Do not know what to do with this element");
169 }
170 }
171*/
172/*
173 for (G4int i=0; i<numEle; i++)
174 {
175 theInelastic.push_back( new G4ParticleHPChannelList );
176 (*theInelastic[i]).Init((*(G4Element::GetElementTable()))[i], dirName, theProjectile);
177 G4int itry = 0;
178 do
179 {
180 (*theInelastic[i]).Register(&theNFS, "F01"); // has
181 (*theInelastic[i]).Register(&theNXFS, "F02");
182 (*theInelastic[i]).Register(&the2NDFS, "F03");
183 (*theInelastic[i]).Register(&the2NFS, "F04"); // has, E Done
184 (*theInelastic[i]).Register(&the3NFS, "F05"); // has, E Done
185 (*theInelastic[i]).Register(&theNAFS, "F06");
186 (*theInelastic[i]).Register(&theN3AFS, "F07");
187 (*theInelastic[i]).Register(&the2NAFS, "F08");
188 (*theInelastic[i]).Register(&the3NAFS, "F09");
189 (*theInelastic[i]).Register(&theNPFS, "F10");
190 (*theInelastic[i]).Register(&theN2AFS, "F11");
191 (*theInelastic[i]).Register(&the2N2AFS, "F12");
192 (*theInelastic[i]).Register(&theNDFS, "F13");
193 (*theInelastic[i]).Register(&theNTFS, "F14");
194 (*theInelastic[i]).Register(&theNHe3FS, "F15");
195 (*theInelastic[i]).Register(&theND2AFS, "F16");
196 (*theInelastic[i]).Register(&theNT2AFS, "F17");
197 (*theInelastic[i]).Register(&the4NFS, "F18"); // has, E Done
198 (*theInelastic[i]).Register(&the2NPFS, "F19");
199 (*theInelastic[i]).Register(&the3NPFS, "F20");
200 (*theInelastic[i]).Register(&theN2PFS, "F21");
201 (*theInelastic[i]).Register(&theNPAFS, "F22");
202 (*theInelastic[i]).Register(&thePFS, "F23");
203 (*theInelastic[i]).Register(&theDFS, "F24");
204 (*theInelastic[i]).Register(&theTFS, "F25");
205 (*theInelastic[i]).Register(&theHe3FS, "F26");
206 (*theInelastic[i]).Register(&theAFS, "F27");
207 (*theInelastic[i]).Register(&the2AFS, "F28");
208 (*theInelastic[i]).Register(&the3AFS, "F29");
209 (*theInelastic[i]).Register(&the2PFS, "F30");
210 (*theInelastic[i]).Register(&thePAFS, "F31");
211 (*theInelastic[i]).Register(&theD2AFS, "F32");
212 (*theInelastic[i]).Register(&theT2AFS, "F33");
213 (*theInelastic[i]).Register(&thePDFS, "F34");
214 (*theInelastic[i]).Register(&thePTFS, "F35");
215 (*theInelastic[i]).Register(&theDAFS, "F36");
216 (*theInelastic[i]).RestartRegistration();
217 itry++;
218 }
219 while( !(*theInelastic[i]).HasDataInAnyFinalState() && itry < 6 );
220 // 6 is corresponding to the value(5) of G4ParticleHPChannel. TK
221
222 if ( itry == 6 )
223 {
224 // No Final State at all.
225 G4bool exceptional = false;
226 if ( (*(G4Element::GetElementTable()))[i]->GetNumberOfIsotopes() == 1 )
227 {
228 if ( (*(G4Element::GetElementTable()))[i]->GetIsotope( 0 )->GetZ() == 1 && (*(G4Element::GetElementTable()))[i]->GetIsotope( 0 )->GetN() == 1 ) exceptional = true; //1H
229 }
230 if ( !exceptional ) {
231 G4cerr << " ELEMENT Z " << (*(G4Element::GetElementTable()))[i]->GetIsotope( 0 )->GetZ() << " N " << (*(G4Element::GetElementTable()))[i]->GetIsotope( 0 )->GetN() << G4endl; //1H
232throw G4HadronicException(__FILE__, __LINE__, "Channel: Do not know what to do with this element");
233 }
234 }
235
236 }
237*/
238
239 }
240
241 G4ParticleHPInelastic::~G4ParticleHPInelastic()
242 {
243// delete [] theInelastic;
244 //Vector is shared, only master deletes
245 if ( !G4Threading::IsWorkerThread() ) {
246 if ( theInelastic != NULL ) {
247 for ( std::vector<G4ParticleHPChannelList*>::iterator
248 it = theInelastic->begin() ; it != theInelastic->end() ; it++ ) {
249 delete *it;
250 }
251 theInelastic->clear();
252 }
253 }
254 }
255
256 #include "G4ParticleHPThermalBoost.hh"
257
258 G4HadFinalState * G4ParticleHPInelastic::ApplyYourself(const G4HadProjectile& aTrack, G4Nucleus& aNucleus )
259 {
260 G4ParticleHPManager::GetInstance()->OpenReactionWhiteBoard();
261 const G4Material * theMaterial = aTrack.GetMaterial();
262 G4int n = theMaterial->GetNumberOfElements();
263 G4int index = theMaterial->GetElement(0)->GetIndex();
264 G4int it=0;
265 if(n!=1)
266 {
267 G4double* xSec = new G4double[n];
268 G4double sum=0;
269 G4int i;
270 const G4double * NumAtomsPerVolume = theMaterial->GetVecNbOfAtomsPerVolume();
271 G4double rWeight;
272 G4ParticleHPThermalBoost aThermalE;
273 for (i=0; i<n; i++)
274 {
275 index = theMaterial->GetElement(i)->GetIndex();
276 rWeight = NumAtomsPerVolume[i];
277 if ( aTrack.GetDefinition() == G4Neutron::Neutron() ) {
278 xSec[i] = ((*theInelastic)[index])->GetXsec(aThermalE.GetThermalEnergy(aTrack,
279 theMaterial->GetElement(i),
280 theMaterial->GetTemperature()));
281 } else {
282 xSec[i] = ((*theInelastic)[index])->GetXsec(aTrack.GetKineticEnergy());
283 }
284 xSec[i] *= rWeight;
285 sum+=xSec[i];
286#ifdef G4PHPDEBUG
287 #ifdef G4VERBOSE
288 if( std::getenv("G4ParticleHPDebug") && G4HadronicParameters::Instance()->GetVerboseLevel() > 0 )
289 G4cout << " G4ParticleHPInelastic XSEC ELEM " << i << " = " << xSec[i] << G4endl;
290 #endif
291#endif
292
293 }
294 G4double random = G4UniformRand();
295 G4double running = 0;
296 for (i=0; i<n; i++)
297 {
298 running += xSec[i];
299 index = theMaterial->GetElement(i)->GetIndex();
300 it = i;
301 //if(random<=running/sum) break;
302 if( sum == 0 || random<=running/sum) break;
303 }
304 delete [] xSec;
305 }
306
307#ifdef G4PHPDEBUG
308 #ifdef G4VERBOSE
309 if( std::getenv("G4ParticleHPDebug") && G4HadronicParameters::Instance()->GetVerboseLevel() > 0 )
310 G4cout << " G4ParticleHPInelastic SELECTED ELEM " << it << " = " << theMaterial->GetElement(it)->GetName() << " FROM MATERIAL " << theMaterial->GetName() << G4endl;
311 #endif
312#endif
313 //return theInelastic[index].ApplyYourself(theMaterial->GetElement(it), aTrack);
314 G4HadFinalState* result = ((*theInelastic)[index])->ApplyYourself(theMaterial->GetElement(it), aTrack);
315 //
316 aNucleus.SetParameters(G4ParticleHPManager::GetInstance()->GetReactionWhiteBoard()->GetTargA(),G4ParticleHPManager::GetInstance()->GetReactionWhiteBoard()->GetTargZ());
317 const G4Element* target_element = (*G4Element::GetElementTable())[index];
318 const G4Isotope* target_isotope=NULL;
319 G4int iele = target_element->GetNumberOfIsotopes();
320 for ( G4int j = 0 ; j != iele ; j++ ) {
321 target_isotope=target_element->GetIsotope( j );
322 if ( target_isotope->GetN() == G4ParticleHPManager::GetInstance()->GetReactionWhiteBoard()->GetTargA() ) break;
323 }
324 //G4cout << "Target Material of this reaction is " << theMaterial->GetName() << G4endl;
325 //G4cout << "Target Element of this reaction is " << target_element->GetName() << G4endl;
326 //G4cout << "Target Isotope of this reaction is " << target_isotope->GetName() << G4endl;
327 aNucleus.SetIsotope( target_isotope );
328
329 G4ParticleHPManager::GetInstance()->CloseReactionWhiteBoard();
330
331 //GDEB
332 if( std::getenv("G4PHPTEST") ) {
333 G4HadSecondary* seco = result->GetSecondary(0);
334 if(seco) {
335 G4ThreeVector secoMom = seco->GetParticle()->GetMomentum();
336 #ifdef G4VERBOSE
337 if ( G4HadronicParameters::Instance()->GetVerboseLevel() > 0 )
338 G4cout << " G4ParticleHPinelastic COS THETA " << std::cos(secoMom.theta()) <<" " << secoMom << G4endl;
339 #endif
340 }
341 }
342
343 return result;
344 }
345
346const std::pair<G4double, G4double> G4ParticleHPInelastic::GetFatalEnergyCheckLevels() const
347{
348 // max energy non-conservation is mass of heavy nucleus
349 // This should be same to the hadron default value
350 return std::pair<G4double, G4double>(10*perCent,DBL_MAX);
351}
352
353/*
354void G4ParticleHPInelastic::addChannelForNewElement()
355{
356 for ( G4int i = numEle ; i < (G4int)G4Element::GetNumberOfElements() ; i++ )
357 {
358 G4cout << "G4ParticleHPInelastic Prepairing Data for the new element of " << (*(G4Element::GetElementTable()))[i]->GetName() << G4endl;
359
360 theInelastic.push_back( new G4ParticleHPChannelList );
361 (*theInelastic[i]).Init((*(G4Element::GetElementTable()))[i], dirName, theProjectile);
362 G4int itry = 0;
363 do
364 {
365 (*theInelastic[i]).Register(&theNFS, "F01"); // has
366 (*theInelastic[i]).Register(&theNXFS, "F02");
367 (*theInelastic[i]).Register(&the2NDFS, "F03");
368 (*theInelastic[i]).Register(&the2NFS, "F04"); // has, E Done
369 (*theInelastic[i]).Register(&the3NFS, "F05"); // has, E Done
370 (*theInelastic[i]).Register(&theNAFS, "F06");
371 (*theInelastic[i]).Register(&theN3AFS, "F07");
372 (*theInelastic[i]).Register(&the2NAFS, "F08");
373 (*theInelastic[i]).Register(&the3NAFS, "F09");
374 (*theInelastic[i]).Register(&theNPFS, "F10");
375 (*theInelastic[i]).Register(&theN2AFS, "F11");
376 (*theInelastic[i]).Register(&the2N2AFS, "F12");
377 (*theInelastic[i]).Register(&theNDFS, "F13");
378 (*theInelastic[i]).Register(&theNTFS, "F14");
379 (*theInelastic[i]).Register(&theNHe3FS, "F15");
380 (*theInelastic[i]).Register(&theND2AFS, "F16");
381 (*theInelastic[i]).Register(&theNT2AFS, "F17");
382 (*theInelastic[i]).Register(&the4NFS, "F18"); // has, E Done
383 (*theInelastic[i]).Register(&the2NPFS, "F19");
384 (*theInelastic[i]).Register(&the3NPFS, "F20");
385 (*theInelastic[i]).Register(&theN2PFS, "F21");
386 (*theInelastic[i]).Register(&theNPAFS, "F22");
387 (*theInelastic[i]).Register(&thePFS, "F23");
388 (*theInelastic[i]).Register(&theDFS, "F24");
389 (*theInelastic[i]).Register(&theTFS, "F25");
390 (*theInelastic[i]).Register(&theHe3FS, "F26");
391 (*theInelastic[i]).Register(&theAFS, "F27");
392 (*theInelastic[i]).Register(&the2AFS, "F28");
393 (*theInelastic[i]).Register(&the3AFS, "F29");
394 (*theInelastic[i]).Register(&the2PFS, "F30");
395 (*theInelastic[i]).Register(&thePAFS, "F31");
396 (*theInelastic[i]).Register(&theD2AFS, "F32");
397 (*theInelastic[i]).Register(&theT2AFS, "F33");
398 (*theInelastic[i]).Register(&thePDFS, "F34");
399 (*theInelastic[i]).Register(&thePTFS, "F35");
400 (*theInelastic[i]).Register(&theDAFS, "F36");
401 (*theInelastic[i]).RestartRegistration();
402 itry++;
403 }
404 while( !(*theInelastic[i]).HasDataInAnyFinalState() && itry < 6 );
405 // 6 is corresponding to the value(5) of G4ParticleHPChannel. TK
406
407 if ( itry == 6 )
408 {
409 // No Final State at all.
410 G4bool exceptional = false;
411 if ( (*(G4Element::GetElementTable()))[i]->GetNumberOfIsotopes() == 1 )
412 {
413 if ( (*(G4Element::GetElementTable()))[i]->GetIsotope( 0 )->GetZ() == 1 && (*(G4Element::GetElementTable()))[i]->GetIsotope( 0 )->GetN() == 1 ) exceptional = true; //1H
414 }
415 if ( !exceptional ) throw G4HadronicException(__FILE__, __LINE__, "Channel: Do not know what to do with this element");
416 }
417 }
418
419 numEle = (G4int)G4Element::GetNumberOfElements();
420}
421*/
422G4int G4ParticleHPInelastic::GetVerboseLevel() const
423{
424 return G4ParticleHPManager::GetInstance()->GetVerboseLevel();
425}
426void G4ParticleHPInelastic::SetVerboseLevel( G4int newValue )
427{
428 G4ParticleHPManager::GetInstance()->SetVerboseLevel(newValue);
429}
430
431#include "G4ParticleHP2AInelasticFS.hh"
432#include "G4ParticleHP2N2AInelasticFS.hh"
433#include "G4ParticleHP2NAInelasticFS.hh"
434#include "G4ParticleHP2NDInelasticFS.hh"
435#include "G4ParticleHP2NInelasticFS.hh"
436#include "G4ParticleHP2NPInelasticFS.hh"
437#include "G4ParticleHP2PInelasticFS.hh"
438#include "G4ParticleHP3AInelasticFS.hh"
439#include "G4ParticleHP3NAInelasticFS.hh"
440#include "G4ParticleHP3NInelasticFS.hh"
441#include "G4ParticleHP3NPInelasticFS.hh"
442#include "G4ParticleHP4NInelasticFS.hh"
443#include "G4ParticleHPAInelasticFS.hh"
444#include "G4ParticleHPD2AInelasticFS.hh"
445#include "G4ParticleHPDAInelasticFS.hh"
446#include "G4ParticleHPDInelasticFS.hh"
447#include "G4ParticleHPHe3InelasticFS.hh"
448#include "G4ParticleHPN2AInelasticFS.hh"
449#include "G4ParticleHPN2PInelasticFS.hh"
450#include "G4ParticleHPN3AInelasticFS.hh"
451#include "G4ParticleHPNAInelasticFS.hh"
452#include "G4ParticleHPND2AInelasticFS.hh"
453#include "G4ParticleHPNDInelasticFS.hh"
454#include "G4ParticleHPNHe3InelasticFS.hh"
455#include "G4ParticleHPNInelasticFS.hh"
456#include "G4ParticleHPNPAInelasticFS.hh"
457#include "G4ParticleHPNPInelasticFS.hh"
458#include "G4ParticleHPNT2AInelasticFS.hh"
459#include "G4ParticleHPNTInelasticFS.hh"
460#include "G4ParticleHPNXInelasticFS.hh"
461#include "G4ParticleHPPAInelasticFS.hh"
462#include "G4ParticleHPPDInelasticFS.hh"
463#include "G4ParticleHPPInelasticFS.hh"
464#include "G4ParticleHPPTInelasticFS.hh"
465#include "G4ParticleHPT2AInelasticFS.hh"
466#include "G4ParticleHPTInelasticFS.hh"
467
468void G4ParticleHPInelastic::BuildPhysicsTable(const G4ParticleDefinition& projectile) {
469
470 G4ParticleHPManager* hpmanager = G4ParticleHPManager::GetInstance();
471
472 theInelastic = hpmanager->GetInelasticFinalStates( &projectile );
473
474 if ( G4Threading::IsMasterThread() ) {
475
476 if ( theInelastic == NULL ) theInelastic = new std::vector<G4ParticleHPChannelList*>;
477
478 if ( numEle == (G4int)G4Element::GetNumberOfElements() ) return;
479
480 if ( theInelastic->size() == G4Element::GetNumberOfElements() ) {
481 numEle = G4Element::GetNumberOfElements();
482 return;
483 }
484/*
485 const char* dataDirVariable;
486 if( &projectile == G4Neutron::Neutron() ) {
487 dataDirVariable = "G4NEUTRONHPDATA";
488 }else if( &projectile == G4Proton::Proton() ) {
489 dataDirVariable = "G4PROTONHPDATA";
490 }else if( &projectile == G4Deuteron::Deuteron() ) {
491 dataDirVariable = "G4DEUTERONHPDATA";
492 }else if( &projectile == G4Triton::Triton() ) {
493 dataDirVariable = "G4TRITONHPDATA";
494 }else if( &projectile == G4He3::He3() ) {
495 dataDirVariable = "G4HE3HPDATA";
496 }else if( &projectile == G4Alpha::Alpha() ) {
497 dataDirVariable = "G4ALPHAHPDATA";
498 } else {
499 G4String message("G4ParticleHPInelastic may only be called for neutron, proton, deuteron, triton, He3 or alpha, while it is called for " + projectile.GetParticleName());
500 throw G4HadronicException(__FILE__, __LINE__,message.c_str());
501 }
502 if(!std::getenv(dataDirVariable)){
503 G4String message("Please set the environment variable " + G4String(dataDirVariable) + " to point to the " + projectile.GetParticleName() + " cross-section files.");
504 throw G4HadronicException(__FILE__, __LINE__,message.c_str());
505 }
506 dirName = std::getenv(dataDirVariable);
507 G4cout << dirName << G4endl;
508
509 G4String tString = "/Inelastic";
510 dirName = dirName + tString;
511
512*/
513 #ifdef G4VERBOSE
514 if ( G4HadronicParameters::Instance()->GetVerboseLevel() > 0 ) {
515 hpmanager->DumpSetting();
516 G4cout << "@@@ G4ParticleHPInelastic instantiated for particle " << projectile.GetParticleName() << " data directory variable is " << dataDirVariable << " pointing to " << dirName << G4endl;
517 }
518 #endif
519 for (G4int i = numEle ; i < (G4int)G4Element::GetNumberOfElements(); i++)
520 {
521 theInelastic->push_back( new G4ParticleHPChannelList );
522 ((*theInelastic)[i])->Init((*(G4Element::GetElementTable()))[i], dirName, const_cast<G4ParticleDefinition*>(&projectile));
523 G4int itry = 0;
524 do
525 {
526 ((*theInelastic)[i])->Register( new G4ParticleHPNInelasticFS , "F01"); // has
527 ((*theInelastic)[i])->Register( new G4ParticleHPNXInelasticFS , "F02");
528 ((*theInelastic)[i])->Register( new G4ParticleHP2NDInelasticFS , "F03");
529 ((*theInelastic)[i])->Register( new G4ParticleHP2NInelasticFS , "F04"); // has, E Done
530 ((*theInelastic)[i])->Register( new G4ParticleHP3NInelasticFS , "F05"); // has, E Done
531 ((*theInelastic)[i])->Register( new G4ParticleHPNAInelasticFS , "F06");
532 ((*theInelastic)[i])->Register( new G4ParticleHPN3AInelasticFS , "F07");
533 ((*theInelastic)[i])->Register( new G4ParticleHP2NAInelasticFS , "F08");
534 ((*theInelastic)[i])->Register( new G4ParticleHP3NAInelasticFS , "F09");
535 ((*theInelastic)[i])->Register( new G4ParticleHPNPInelasticFS , "F10");
536 ((*theInelastic)[i])->Register( new G4ParticleHPN2AInelasticFS , "F11");
537 ((*theInelastic)[i])->Register( new G4ParticleHP2N2AInelasticFS , "F12");
538 ((*theInelastic)[i])->Register( new G4ParticleHPNDInelasticFS , "F13");
539 ((*theInelastic)[i])->Register( new G4ParticleHPNTInelasticFS , "F14");
540 ((*theInelastic)[i])->Register( new G4ParticleHPNHe3InelasticFS , "F15");
541 ((*theInelastic)[i])->Register( new G4ParticleHPND2AInelasticFS , "F16");
542 ((*theInelastic)[i])->Register( new G4ParticleHPNT2AInelasticFS , "F17");
543 ((*theInelastic)[i])->Register( new G4ParticleHP4NInelasticFS , "F18"); // has, E Done
544 ((*theInelastic)[i])->Register( new G4ParticleHP2NPInelasticFS , "F19");
545 ((*theInelastic)[i])->Register( new G4ParticleHP3NPInelasticFS , "F20");
546 ((*theInelastic)[i])->Register( new G4ParticleHPN2PInelasticFS , "F21");
547 ((*theInelastic)[i])->Register( new G4ParticleHPNPAInelasticFS , "F22");
548 ((*theInelastic)[i])->Register( new G4ParticleHPPInelasticFS , "F23");
549 ((*theInelastic)[i])->Register( new G4ParticleHPDInelasticFS , "F24");
550 ((*theInelastic)[i])->Register( new G4ParticleHPTInelasticFS , "F25");
551 ((*theInelastic)[i])->Register( new G4ParticleHPHe3InelasticFS , "F26");
552 ((*theInelastic)[i])->Register( new G4ParticleHPAInelasticFS , "F27");
553 ((*theInelastic)[i])->Register( new G4ParticleHP2AInelasticFS , "F28");
554 ((*theInelastic)[i])->Register( new G4ParticleHP3AInelasticFS , "F29");
555 ((*theInelastic)[i])->Register( new G4ParticleHP2PInelasticFS , "F30");
556 ((*theInelastic)[i])->Register( new G4ParticleHPPAInelasticFS , "F31");
557 ((*theInelastic)[i])->Register( new G4ParticleHPD2AInelasticFS , "F32");
558 ((*theInelastic)[i])->Register( new G4ParticleHPT2AInelasticFS , "F33");
559 ((*theInelastic)[i])->Register( new G4ParticleHPPDInelasticFS , "F34");
560 ((*theInelastic)[i])->Register( new G4ParticleHPPTInelasticFS , "F35");
561 ((*theInelastic)[i])->Register( new G4ParticleHPDAInelasticFS , "F36");
562 ((*theInelastic)[i])->RestartRegistration();
563 itry++;
564 }
565 while( !((*theInelastic)[i])->HasDataInAnyFinalState() && itry < 6 ); // Loop checking, 11.05.2015, T. Koi
566 // 6 is corresponding to the value(5) of G4ParticleHPChannel. TK
567
568 if ( itry == 6 ) {
569 // No Final State at all.
570 /*
571 G4bool exceptional = false;
572 if ( (*(G4Element::GetElementTable()))[i]->GetNumberOfIsotopes() == 1 )
573 {
574 if ( (*(G4Element::GetElementTable()))[i]->GetIsotope( 0 )->GetZ() == 1 && (*(G4Element::GetElementTable()))[i]->GetIsotope( 0 )->GetN() == 1 ) exceptional = true; //1H
575 }
576 if ( !exceptional ) {
577 G4cerr << " ELEMENT Z " << (*(G4Element::GetElementTable()))[i]->GetIsotope( 0 )->GetZ() << " N " << (*(G4Element::GetElementTable()))[i]->GetIsotope( 0 )->GetN() << G4endl; //1H
578 throw G4HadronicException(__FILE__, __LINE__, "Channel: Do not know what to do with this element");
579 }
580 */
581 #ifdef G4VERBOSE
582 if ( G4ParticleHPManager::GetInstance()->GetVerboseLevel() > 1 &&
583 G4HadronicParameters::Instance()->GetVerboseLevel() > 0 ) {
584 G4cout << "ParticleHP::Inelastic for " << projectile.GetParticleName() << ". Do not know what to do with element of \"" << (*(G4Element::GetElementTable()))[i]->GetName() << "\"." << G4endl;
585 G4cout << "The components of the element are" << G4endl;
586 //G4cout << "TKDB dataDirVariable = " << dataDirVariable << G4endl;
587 for ( G4int ii = 0 ; ii < (G4int)( (*(G4Element::GetElementTable()))[i]->GetNumberOfIsotopes() ) ; ii++ ) {
588 G4cout << " Z: " << (*(G4Element::GetElementTable()))[i]->GetIsotope( ii )->GetZ() << ", A: " << (*(G4Element::GetElementTable()))[i]->GetIsotope( ii )->GetN() << G4endl;
589 }
590 G4cout << "No possible final state data of the element is found in " << dataDirVariable << "." << G4endl;
591 }
592 #endif
593 }
594 }
595 hpmanager->RegisterInelasticFinalStates( &projectile , theInelastic );
596 }
597 numEle = G4Element::GetNumberOfElements();
598}
599
600void G4ParticleHPInelastic::ModelDescription(std::ostream& outFile) const
601{
602 outFile << "Extension of High Precision model for inelastic reaction of proton, deuteron, triton, He3 and alpha below 20MeV\n";
603}
604
G4double
double G4double
Definition: G4Types.hh:83
G4int
int G4int
Definition: G4Types.hh:85
G4endl
#define G4endl
Definition: G4ios.hh:57
G4cout
G4GLOB_DLL std::ostream G4cout
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:52
CLHEP::Hep3Vector::theta
double theta() const
G4Alpha::Alpha
static G4Alpha * Alpha()
Definition: G4Alpha.cc:88
G4Deuteron::Deuteron
static G4Deuteron * Deuteron()
Definition: G4Deuteron.cc:93
G4DynamicParticle::GetMomentum
G4ThreeVector GetMomentum() const
G4Element::GetElementTable
static G4ElementTable * GetElementTable()
Definition: G4Element.cc:397
G4Element::GetNumberOfElements
static size_t GetNumberOfElements()
Definition: G4Element.cc:404
G4Element::GetIsotope
const G4Isotope * GetIsotope(G4int iso) const
Definition: G4Element.hh:169
G4Element::GetIndex
size_t GetIndex() const
Definition: G4Element.hh:181
G4Element::GetNumberOfIsotopes
size_t GetNumberOfIsotopes() const
Definition: G4Element.hh:158
G4Element::GetName
const G4String & GetName() const
Definition: G4Element.hh:126
G4HadFinalState::GetSecondary
G4HadSecondary * GetSecondary(size_t i)
Definition: G4HadFinalState.cc:79
G4HadProjectile::GetMaterial
const G4Material * GetMaterial() const
Definition: G4HadProjectile.hh:81
G4HadProjectile::GetDefinition
const G4ParticleDefinition * GetDefinition() const
Definition: G4HadProjectile.hh:86
G4HadProjectile::GetKineticEnergy
G4double GetKineticEnergy() const
Definition: G4HadProjectile.hh:116
G4HadSecondary::GetParticle
G4DynamicParticle * GetParticle()
Definition: G4HadSecondary.hh:42
G4HadronicInteraction::SetMinEnergy
void SetMinEnergy(G4double anEnergy)
Definition: G4HadronicInteraction.hh:89
G4HadronicInteraction::SetMaxEnergy
void SetMaxEnergy(const G4double anEnergy)
Definition: G4HadronicInteraction.hh:102
G4HadronicParameters::Instance
static G4HadronicParameters * Instance()
Definition: G4HadronicParameters.cc:52
G4He3::He3
static G4He3 * He3()
Definition: G4He3.cc:93
G4Isotope::GetN
G4int GetN() const
Definition: G4Isotope.hh:93
G4Material::GetTemperature
G4double GetTemperature() const
Definition: G4Material.hh:180
G4Material::GetElement
const G4Element * GetElement(G4int iel) const
Definition: G4Material.hh:200
G4Material::GetNumberOfElements
size_t GetNumberOfElements() const
Definition: G4Material.hh:184
G4Material::GetVecNbOfAtomsPerVolume
const G4double * GetVecNbOfAtomsPerVolume() const
Definition: G4Material.hh:204
G4Material::GetName
const G4String & GetName() const
Definition: G4Material.hh:175
G4Neutron::Neutron
static G4Neutron * Neutron()
Definition: G4Neutron.cc:103
G4Nucleus::SetParameters
void SetParameters(const G4double A, const G4double Z)
Definition: G4Nucleus.cc:212
G4Nucleus::SetIsotope
void SetIsotope(const G4Isotope *iso)
Definition: G4Nucleus.hh:122
G4ParticleDefinition::GetParticleName
const G4String & GetParticleName() const
Definition: G4ParticleDefinition.hh:108
G4ParticleHPInelastic::G4ParticleHPInelastic
G4ParticleHPInelastic(G4ParticleDefinition *projectile=G4Neutron::Neutron(), const char *name="NeutronHPInelastic")
Definition: G4ParticleHPInelastic.cc:47
G4ParticleHPInelastic::ApplyYourself
G4HadFinalState * ApplyYourself(const G4HadProjectile &aTrack, G4Nucleus &aTargetNucleus)
Definition: G4ParticleHPInelastic.cc:258
G4ParticleHPInelastic::BuildPhysicsTable
void BuildPhysicsTable(const G4ParticleDefinition &)
Definition: G4ParticleHPInelastic.cc:468
G4ParticleHPInelastic::GetFatalEnergyCheckLevels
virtual const std::pair< G4double, G4double > GetFatalEnergyCheckLevels() const
Definition: G4ParticleHPInelastic.cc:346
G4ParticleHPInelastic::theInelastic
std::vector< G4ParticleHPChannelList * > * theInelastic
Definition: G4ParticleHPInelastic.hh:109
G4ParticleHPInelastic::ModelDescription
virtual void ModelDescription(std::ostream &outFile) const
Definition: G4ParticleHPInelastic.cc:600
G4ParticleHPManager::GetInelasticFinalStates
std::vector< G4ParticleHPChannelList * > * GetInelasticFinalStates(const G4ParticleDefinition *)
Definition: G4ParticleHPManager.cc:309
G4ParticleHPManager::RegisterInelasticFinalStates
void RegisterInelasticFinalStates(const G4ParticleDefinition *, std::vector< G4ParticleHPChannelList * > *)
Definition: G4ParticleHPManager.cc:317
G4ParticleHPManager::GetInstance
static G4ParticleHPManager * GetInstance()
Definition: G4ParticleHPManager.cc:140
G4ParticleHPManager::GetReactionWhiteBoard
G4ParticleHPReactionWhiteBoard * GetReactionWhiteBoard()
Definition: G4ParticleHPManager.cc:155
G4ParticleHPThermalBoost::GetThermalEnergy
G4double GetThermalEnergy(const G4HadProjectile &aP, const G4Element *anE, G4double aT)
Definition: G4ParticleHPThermalBoost.hh:45
G4Proton::Proton
static G4Proton * Proton()
Definition: G4Proton.cc:92
G4Triton::Triton
static G4Triton * Triton()
Definition: G4Triton.cc:94
G4Threading::IsWorkerThread
G4bool IsWorkerThread()
Definition: G4Threading.cc:123
G4Threading::IsMasterThread
G4bool IsMasterThread()
Definition: G4Threading.cc:124
DBL_MAX
#define DBL_MAX
Definition: templates.hh:62