| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // this code implementation is the intellectua 26 // this code implementation is the intellectual property of
27 // neutron_hp -- source file 27 // neutron_hp -- source file
28 // J.P. Wellisch, Nov-1996 28 // J.P. Wellisch, Nov-1996
29 // A prototype of the low energy neutron trans 29 // A prototype of the low energy neutron transport model.
30 // 30 //
31 // By copying, distributing or modifying the P 31 // By copying, distributing or modifying the Program (or any work
32 // based on the Program) you indicate your acc 32 // based on the Program) you indicate your acceptance of this statement,
33 // and all its terms. 33 // and all its terms.
34 // 34 //
>> 35 // $Id$
35 // 36 //
36 // 070523 bug fix for G4FPE_DEBUG on by A. How 37 // 070523 bug fix for G4FPE_DEBUG on by A. Howard (and T. Koi)
37 // 081203 limit maximum trial for creating fin 38 // 081203 limit maximum trial for creating final states add protection for 1H isotope case by T. Koi
38 // 39 //
39 // P. Arce, June-2014 Conversion neutron_hp to 40 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
40 // V. Ivanchenko, July-2023 Basic revision of <<
41 // 41 //
42 #include "G4ParticleHPInelastic.hh" 42 #include "G4ParticleHPInelastic.hh"
>> 43 #include "G4SystemOfUnits.hh"
>> 44 #include "G4ParticleHPManager.hh"
>> 45 #include "G4Threading.hh"
>> 46
>> 47 G4ParticleHPInelastic::G4ParticleHPInelastic(G4ParticleDefinition* projectile, const char* name )
>> 48 :G4HadronicInteraction(name)
>> 49 ,theInelastic(NULL)
>> 50 ,numEle(0)
>> 51 ,theProjectile(projectile)
>> 52 {
>> 53 const char* dataDirVariable;
>> 54 if( theProjectile == G4Neutron::Neutron() ) {
>> 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
>> 220 G4ParticleHPInelastic::~G4ParticleHPInelastic()
>> 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
>> 235 G4HadFinalState * G4ParticleHPInelastic::ApplyYourself(const G4HadProjectile& aTrack, G4Nucleus& aNucleus )
>> 236 {
>> 237 G4ParticleHPManager::GetInstance()->OpenReactionWhiteBoard();
>> 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
>> 300 G4ParticleHPManager::GetInstance()->CloseReactionWhiteBoard();
>> 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 */
>> 392 G4int G4ParticleHPInelastic::GetVerboseLevel() const
>> 393 {
>> 394 return G4ParticleHPManager::GetInstance()->GetVerboseLevel();
>> 395 }
>> 396 void G4ParticleHPInelastic::SetVerboseLevel( G4int newValue )
>> 397 {
>> 398 G4ParticleHPManager::GetInstance()->SetVerboseLevel(newValue);
>> 399 }
43 400
44 #include "G4HadronicParameters.hh" <<
45 #include "G4ParticleHP2AInelasticFS.hh" 401 #include "G4ParticleHP2AInelasticFS.hh"
46 #include "G4ParticleHP2N2AInelasticFS.hh" 402 #include "G4ParticleHP2N2AInelasticFS.hh"
47 #include "G4ParticleHP2NAInelasticFS.hh" 403 #include "G4ParticleHP2NAInelasticFS.hh"
48 #include "G4ParticleHP2NDInelasticFS.hh" 404 #include "G4ParticleHP2NDInelasticFS.hh"
49 #include "G4ParticleHP2NInelasticFS.hh" 405 #include "G4ParticleHP2NInelasticFS.hh"
50 #include "G4ParticleHP2NPInelasticFS.hh" 406 #include "G4ParticleHP2NPInelasticFS.hh"
51 #include "G4ParticleHP2PInelasticFS.hh" 407 #include "G4ParticleHP2PInelasticFS.hh"
52 #include "G4ParticleHP3AInelasticFS.hh" 408 #include "G4ParticleHP3AInelasticFS.hh"
53 #include "G4ParticleHP3NAInelasticFS.hh" 409 #include "G4ParticleHP3NAInelasticFS.hh"
54 #include "G4ParticleHP3NInelasticFS.hh" 410 #include "G4ParticleHP3NInelasticFS.hh"
55 #include "G4ParticleHP3NPInelasticFS.hh" 411 #include "G4ParticleHP3NPInelasticFS.hh"
56 #include "G4ParticleHP4NInelasticFS.hh" 412 #include "G4ParticleHP4NInelasticFS.hh"
57 #include "G4ParticleHPAInelasticFS.hh" 413 #include "G4ParticleHPAInelasticFS.hh"
58 #include "G4ParticleHPD2AInelasticFS.hh" 414 #include "G4ParticleHPD2AInelasticFS.hh"
59 #include "G4ParticleHPDAInelasticFS.hh" 415 #include "G4ParticleHPDAInelasticFS.hh"
60 #include "G4ParticleHPDInelasticFS.hh" 416 #include "G4ParticleHPDInelasticFS.hh"
61 #include "G4ParticleHPHe3InelasticFS.hh" 417 #include "G4ParticleHPHe3InelasticFS.hh"
62 #include "G4ParticleHPManager.hh" <<
63 #include "G4ParticleHPN2AInelasticFS.hh" 418 #include "G4ParticleHPN2AInelasticFS.hh"
64 #include "G4ParticleHPN2PInelasticFS.hh" 419 #include "G4ParticleHPN2PInelasticFS.hh"
65 #include "G4ParticleHPN3AInelasticFS.hh" 420 #include "G4ParticleHPN3AInelasticFS.hh"
66 #include "G4ParticleHPNAInelasticFS.hh" 421 #include "G4ParticleHPNAInelasticFS.hh"
67 #include "G4ParticleHPND2AInelasticFS.hh" 422 #include "G4ParticleHPND2AInelasticFS.hh"
68 #include "G4ParticleHPNDInelasticFS.hh" 423 #include "G4ParticleHPNDInelasticFS.hh"
69 #include "G4ParticleHPNHe3InelasticFS.hh" 424 #include "G4ParticleHPNHe3InelasticFS.hh"
70 #include "G4ParticleHPNInelasticFS.hh" 425 #include "G4ParticleHPNInelasticFS.hh"
71 #include "G4ParticleHPNPAInelasticFS.hh" 426 #include "G4ParticleHPNPAInelasticFS.hh"
72 #include "G4ParticleHPNPInelasticFS.hh" 427 #include "G4ParticleHPNPInelasticFS.hh"
73 #include "G4ParticleHPNT2AInelasticFS.hh" 428 #include "G4ParticleHPNT2AInelasticFS.hh"
74 #include "G4ParticleHPNTInelasticFS.hh" 429 #include "G4ParticleHPNTInelasticFS.hh"
75 #include "G4ParticleHPNXInelasticFS.hh" 430 #include "G4ParticleHPNXInelasticFS.hh"
76 #include "G4ParticleHPPAInelasticFS.hh" 431 #include "G4ParticleHPPAInelasticFS.hh"
77 #include "G4ParticleHPPDInelasticFS.hh" 432 #include "G4ParticleHPPDInelasticFS.hh"
78 #include "G4ParticleHPPInelasticFS.hh" 433 #include "G4ParticleHPPInelasticFS.hh"
79 #include "G4ParticleHPPTInelasticFS.hh" 434 #include "G4ParticleHPPTInelasticFS.hh"
80 #include "G4ParticleHPT2AInelasticFS.hh" 435 #include "G4ParticleHPT2AInelasticFS.hh"
81 #include "G4ParticleHPTInelasticFS.hh" 436 #include "G4ParticleHPTInelasticFS.hh"
82 #include "G4ParticleHPThermalBoost.hh" <<
83 #include "G4SystemOfUnits.hh" <<
84 #include "G4AutoLock.hh" <<
85 <<
86 G4bool G4ParticleHPInelastic::fLock[] = {true, <<
87 std::vector<G4ParticleHPChannelList*>* <<
88 G4ParticleHPInelastic::theInelastic[] = {nullp <<
89 <<
90 namespace <<
91 { <<
92 G4Mutex theHPInelastic = G4MUTEX_INITIALIZER <<
93 } <<
94 <<
95 G4ParticleHPInelastic::G4ParticleHPInelastic(G <<
96 : G4HadronicInteraction(name), theProjectile <<
97 { <<
98 fManager = G4ParticleHPManager::GetInstance( <<
99 dirName = fManager->GetParticleHPPath(thePro <<
100 indexP = fManager->GetPHPIndex(theProjectile <<
101 <<
102 #ifdef G4VERBOSE <<
103 if (fManager->GetVerboseLevel() > 1) <<
104 G4cout << "@@@ G4ParticleHPInelastic insta <<
105 << p->GetParticleName() << " indexP <<
106 << "/n data directory " << dirName << <<
107 #endif <<
108 } <<
109 <<
110 G4ParticleHPInelastic::~G4ParticleHPInelastic( <<
111 { <<
112 // Vector is shared, only one delete <<
113 if (isFirst) { <<
114 ClearData(); <<
115 } <<
116 } <<
117 <<
118 void G4ParticleHPInelastic::ClearData() <<
119 { <<
120 if (theInelastic[indexP] != nullptr) { <<
121 for (auto const& p : *(theInelastic[indexP <<
122 delete p; <<
123 } <<
124 delete theInelastic[indexP]; <<
125 theInelastic[indexP] = nullptr; <<
126 } <<
127 } <<
128 <<
129 G4HadFinalState* G4ParticleHPInelastic::ApplyY <<
130 <<
131 { <<
132 G4ParticleHPManager::GetInstance()->OpenReac <<
133 const G4Material* theMaterial = aTrack.GetMa <<
134 auto n = (G4int)theMaterial->GetNumberOfElem <<
135 auto elm = theMaterial->GetElement(0); <<
136 std::size_t index = elm->GetIndex(); <<
137 G4int it = 0; <<
138 /* <<
139 G4cout << "G4ParticleHPInelastic::ApplyYours <<
140 << " indexP=" << indexP << " " <<
141 << aTrack.GetDefinition()->GetParticl <<
142 */ <<
143 if (n != 1) { <<
144 auto xSec = new G4double[n]; <<
145 G4double sum = 0; <<
146 G4int i; <<
147 const G4double* NumAtomsPerVolume = theMat <<
148 G4double rWeight; <<
149 G4double xs; <<
150 G4ParticleHPThermalBoost aThermalE; <<
151 for (i = 0; i < n; ++i) { <<
152 elm = theMaterial->GetElement(i); <<
153 index = elm->GetIndex(); <<
154 /* <<
155 G4cout << "i=" << i << " index=" << ind <<
156 << " " << (*(theInelastic[indexP]))[in <<
157 */ <<
158 rWeight = NumAtomsPerVolume[i]; <<
159 if (aTrack.GetDefinition() == G4Neutron: <<
160 xs = (*(theInelastic[indexP]))[index]- <<
161 theMaterial->GetTemperature <<
162 } <<
163 else { <<
164 xs = (*(theInelastic[indexP]))[index]- <<
165 } <<
166 xs *= rWeight; <<
167 sum += xs; <<
168 xSec[i] = sum; <<
169 #ifdef G4VERBOSE <<
170 if (fManager->GetDEBUG()) <<
171 G4cout << " G4ParticleHPInelastic XSEC <<
172 #endif <<
173 } <<
174 sum *= G4UniformRand(); <<
175 for (it = 0; it < n; ++it) { <<
176 elm = theMaterial->GetElement(it); <<
177 index = elm->GetIndex(); <<
178 if (sum <= xSec[it]) break; <<
179 } <<
180 delete[] xSec; <<
181 } <<
182 437
183 #ifdef G4VERBOSE << 438 void G4ParticleHPInelastic::BuildPhysicsTable(const G4ParticleDefinition& projectile) {
184 if (fManager->GetDEBUG()) <<
185 G4cout << " G4ParticleHPInelastic: Elem it <<
186 << elm->GetName() << " index=" << i <<
187 << " from material " << theMaterial->GetN <<
188 << G4endl; <<
189 #endif <<
190 <<
191 G4HadFinalState* result = <<
192 (*(theInelastic[indexP]))[index]->ApplyYou <<
193 439
194 aNucleus.SetParameters(fManager->GetReaction << 440 G4ParticleHPManager* hpmanager = G4ParticleHPManager::GetInstance();
195 fManager->GetReaction <<
196 <<
197 const G4Element* target_element = (*G4Elemen <<
198 const G4Isotope* target_isotope = nullptr; <<
199 auto iele = (G4int)target_element->GetNumber <<
200 for (G4int j = 0; j != iele; ++j) { <<
201 target_isotope = target_element->GetIsotop <<
202 if (target_isotope->GetN() == fManager->Ge <<
203 break; <<
204 } <<
205 aNucleus.SetIsotope(target_isotope); <<
206 441
207 G4ParticleHPManager::GetInstance()->CloseRea << 442 theInelastic = hpmanager->GetInelasticFinalStates( &projectile );
208 443
209 return result; << 444 if ( G4Threading::IsMasterThread() ) {
210 } <<
211 445
212 const std::pair<G4double, G4double> G4Particle << 446 if ( theInelastic == NULL ) theInelastic = new std::vector<G4ParticleHPChannelList*>;
213 { <<
214 // max energy non-conservation is mass of he <<
215 return std::pair<G4double, G4double>(10.0 * <<
216 } <<
217 447
218 void G4ParticleHPInelastic::BuildPhysicsTable( << 448 if ( numEle == (G4int)G4Element::GetNumberOfElements() ) return;
219 { <<
220 if (fLock[indexP]) { <<
221 G4AutoLock l(&theHPInelastic); <<
222 if (fLock[indexP]) { <<
223 isFirst = true; <<
224 fLock[indexP] = false; <<
225 } <<
226 l.unlock(); <<
227 } <<
228 449
229 G4int nelm = (G4int)G4Element::GetNumberOfEl << 450 if ( theInelastic->size() == G4Element::GetNumberOfElements() ) {
230 G4int n0 = numEle; << 451 numEle = G4Element::GetNumberOfElements();
231 numEle = nelm; << 452 return;
232 if (!isFirst || nelm == n0) { return; } << 453 }
233 454
234 // extra elements should be initialized << 455 const char* dataDirVariable;
235 G4AutoLock l(&theHPInelastic); << 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;
236 478
237 if (nullptr == theInelastic[indexP]) { << 479 G4String tString = "/Inelastic";
238 theInelastic[indexP] = new std::vector<G4P << 480 dirName = dirName + tString;
239 } <<
240 481
241 if (fManager->GetVerboseLevel() > 0 && isFir << 482 G4cout << "@@@ G4ParticleHPInelastic instantiated for particle " << projectile.GetParticleName() << " data directory variable is " << dataDirVariable << " pointing to " << dirName << G4endl;
242 fManager->DumpSetting(); <<
243 G4cout << "@@@ G4ParticleHPInelastic insta <<
244 << theProjectile->GetParticleName() << "/ <<
245 << dirName << G4endl; <<
246 } <<
247 483
248 auto table = G4Element::GetElementTable(); << 484 for (G4int i = numEle ; i < (G4int)G4Element::GetNumberOfElements(); i++)
249 for (G4int i = n0; i < nelm; ++i) { << 485 {
250 auto clist = new G4ParticleHPChannelList(3 << 486 theInelastic->push_back( new G4ParticleHPChannelList );
251 theInelastic[indexP]->push_back(clist); << 487 ((*theInelastic)[i])->Init((*(G4Element::GetElementTable()))[i], dirName, const_cast<G4ParticleDefinition*>(&projectile));
252 clist->Init((*table)[i], dirName, theProje << 488 G4int itry = 0;
253 clist->Register(new G4ParticleHPNInelastic << 489 do
254 clist->Register(new G4ParticleHPNXInelasti << 490 {
255 clist->Register(new G4ParticleHP2NDInelast << 491 ((*theInelastic)[i])->Register( new G4ParticleHPNInelasticFS , "F01"); // has
256 clist->Register(new G4ParticleHP2NInelasti << 492 ((*theInelastic)[i])->Register( new G4ParticleHPNXInelasticFS , "F02");
257 clist->Register(new G4ParticleHP3NInelasti << 493 ((*theInelastic)[i])->Register( new G4ParticleHP2NDInelasticFS , "F03");
258 clist->Register(new G4ParticleHPNAInelasti << 494 ((*theInelastic)[i])->Register( new G4ParticleHP2NInelasticFS , "F04"); // has, E Done
259 clist->Register(new G4ParticleHPN3AInelast << 495 ((*theInelastic)[i])->Register( new G4ParticleHP3NInelasticFS , "F05"); // has, E Done
260 clist->Register(new G4ParticleHP2NAInelast << 496 ((*theInelastic)[i])->Register( new G4ParticleHPNAInelasticFS , "F06");
261 clist->Register(new G4ParticleHP3NAInelast << 497 ((*theInelastic)[i])->Register( new G4ParticleHPN3AInelasticFS , "F07");
262 clist->Register(new G4ParticleHPNPInelasti << 498 ((*theInelastic)[i])->Register( new G4ParticleHP2NAInelasticFS , "F08");
263 clist->Register(new G4ParticleHPN2AInelast << 499 ((*theInelastic)[i])->Register( new G4ParticleHP3NAInelasticFS , "F09");
264 clist->Register(new G4ParticleHP2N2AInelas << 500 ((*theInelastic)[i])->Register( new G4ParticleHPNPInelasticFS , "F10");
265 clist->Register(new G4ParticleHPNDInelasti << 501 ((*theInelastic)[i])->Register( new G4ParticleHPN2AInelasticFS , "F11");
266 clist->Register(new G4ParticleHPNTInelasti << 502 ((*theInelastic)[i])->Register( new G4ParticleHP2N2AInelasticFS , "F12");
267 clist->Register(new G4ParticleHPNHe3Inelas << 503 ((*theInelastic)[i])->Register( new G4ParticleHPNDInelasticFS , "F13");
268 clist->Register(new G4ParticleHPND2AInelas << 504 ((*theInelastic)[i])->Register( new G4ParticleHPNTInelasticFS , "F14");
269 clist->Register(new G4ParticleHPNT2AInelas << 505 ((*theInelastic)[i])->Register( new G4ParticleHPNHe3InelasticFS , "F15");
270 clist->Register(new G4ParticleHP4NInelasti << 506 ((*theInelastic)[i])->Register( new G4ParticleHPND2AInelasticFS , "F16");
271 clist->Register(new G4ParticleHP2NPInelast << 507 ((*theInelastic)[i])->Register( new G4ParticleHPNT2AInelasticFS , "F17");
272 clist->Register(new G4ParticleHP3NPInelast << 508 ((*theInelastic)[i])->Register( new G4ParticleHP4NInelasticFS , "F18"); // has, E Done
273 clist->Register(new G4ParticleHPN2PInelast << 509 ((*theInelastic)[i])->Register( new G4ParticleHP2NPInelasticFS , "F19");
274 clist->Register(new G4ParticleHPNPAInelast << 510 ((*theInelastic)[i])->Register( new G4ParticleHP3NPInelasticFS , "F20");
275 clist->Register(new G4ParticleHPPInelastic << 511 ((*theInelastic)[i])->Register( new G4ParticleHPN2PInelasticFS , "F21");
276 clist->Register(new G4ParticleHPDInelastic << 512 ((*theInelastic)[i])->Register( new G4ParticleHPNPAInelasticFS , "F22");
277 clist->Register(new G4ParticleHPTInelastic << 513 ((*theInelastic)[i])->Register( new G4ParticleHPPInelasticFS , "F23");
278 clist->Register(new G4ParticleHPHe3Inelast << 514 ((*theInelastic)[i])->Register( new G4ParticleHPDInelasticFS , "F24");
279 clist->Register(new G4ParticleHPAInelastic << 515 ((*theInelastic)[i])->Register( new G4ParticleHPTInelasticFS , "F25");
280 clist->Register(new G4ParticleHP2AInelasti << 516 ((*theInelastic)[i])->Register( new G4ParticleHPHe3InelasticFS , "F26");
281 clist->Register(new G4ParticleHP3AInelasti << 517 ((*theInelastic)[i])->Register( new G4ParticleHPAInelasticFS , "F27");
282 clist->Register(new G4ParticleHP2PInelasti << 518 ((*theInelastic)[i])->Register( new G4ParticleHP2AInelasticFS , "F28");
283 clist->Register(new G4ParticleHPPAInelasti << 519 ((*theInelastic)[i])->Register( new G4ParticleHP3AInelasticFS , "F29");
284 clist->Register(new G4ParticleHPD2AInelast << 520 ((*theInelastic)[i])->Register( new G4ParticleHP2PInelasticFS , "F30");
285 clist->Register(new G4ParticleHPT2AInelast << 521 ((*theInelastic)[i])->Register( new G4ParticleHPPAInelasticFS , "F31");
286 clist->Register(new G4ParticleHPPDInelasti << 522 ((*theInelastic)[i])->Register( new G4ParticleHPD2AInelasticFS , "F32");
287 clist->Register(new G4ParticleHPPTInelasti << 523 ((*theInelastic)[i])->Register( new G4ParticleHPT2AInelasticFS , "F33");
288 clist->Register(new G4ParticleHPDAInelasti << 524 ((*theInelastic)[i])->Register( new G4ParticleHPPDInelasticFS , "F34");
289 #ifdef G4VERBOSE << 525 ((*theInelastic)[i])->Register( new G4ParticleHPPTInelasticFS , "F35");
290 if (fManager->GetVerboseLevel() > 1) { << 526 ((*theInelastic)[i])->Register( new G4ParticleHPDAInelasticFS , "F36");
291 G4cout << "ParticleHP::Inelastic for " << 527 ((*theInelastic)[i])->RestartRegistration();
292 << theProjectile->GetParticleName() << << 528 itry++;
293 << (*table)[i]->GetName() << G4endl; << 529 }
294 } << 530 while( !((*theInelastic)[i])->HasDataInAnyFinalState() && itry < 6 ); // Loop checking, 11.05.2015, T. Koi
295 #endif << 531 // 6 is corresponding to the value(5) of G4ParticleHPChannel. TK
296 } << 532
297 fManager->RegisterInelasticFinalStates( &pro << 533 if ( itry == 6 ) {
298 l.unlock(); << 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 */
>> 546 if ( G4ParticleHPManager::GetInstance()->GetVerboseLevel() > 1 ) {
>> 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 }
>> 559 numEle = G4Element::GetNumberOfElements();
299 } 560 }
300 561
301 void G4ParticleHPInelastic::ModelDescription(s 562 void G4ParticleHPInelastic::ModelDescription(std::ostream& outFile) const
302 { 563 {
303 outFile << "High Precision (HP) model for in << 564 outFile << "Extension of High Precision model for inelastic reaction of proton, deuteron, triton, He3 and alpha below 20MeV\n";
304 << theProjectile->GetParticleName() << " b <<
305 } 565 }
>> 566
306 567