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