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