| 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 // 35 //
36 // 070523 bug fix for G4FPE_DEBUG on by A. How 36 // 070523 bug fix for G4FPE_DEBUG on by A. Howard (and T. Koi)
37 // 081203 limit maximum trial for creating fin 37 // 081203 limit maximum trial for creating final states add protection for 1H isotope case by T. Koi
38 // 38 //
39 // P. Arce, June-2014 Conversion neutron_hp to 39 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
40 // V. Ivanchenko, July-2023 Basic revision of <<
41 // 40 //
42 #include "G4ParticleHPInelastic.hh" 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 }
43 416
44 #include "G4HadronicParameters.hh" <<
45 #include "G4ParticleHP2AInelasticFS.hh" 417 #include "G4ParticleHP2AInelasticFS.hh"
46 #include "G4ParticleHP2N2AInelasticFS.hh" 418 #include "G4ParticleHP2N2AInelasticFS.hh"
47 #include "G4ParticleHP2NAInelasticFS.hh" 419 #include "G4ParticleHP2NAInelasticFS.hh"
48 #include "G4ParticleHP2NDInelasticFS.hh" 420 #include "G4ParticleHP2NDInelasticFS.hh"
49 #include "G4ParticleHP2NInelasticFS.hh" 421 #include "G4ParticleHP2NInelasticFS.hh"
50 #include "G4ParticleHP2NPInelasticFS.hh" 422 #include "G4ParticleHP2NPInelasticFS.hh"
51 #include "G4ParticleHP2PInelasticFS.hh" 423 #include "G4ParticleHP2PInelasticFS.hh"
52 #include "G4ParticleHP3AInelasticFS.hh" 424 #include "G4ParticleHP3AInelasticFS.hh"
53 #include "G4ParticleHP3NAInelasticFS.hh" 425 #include "G4ParticleHP3NAInelasticFS.hh"
54 #include "G4ParticleHP3NInelasticFS.hh" 426 #include "G4ParticleHP3NInelasticFS.hh"
55 #include "G4ParticleHP3NPInelasticFS.hh" 427 #include "G4ParticleHP3NPInelasticFS.hh"
56 #include "G4ParticleHP4NInelasticFS.hh" 428 #include "G4ParticleHP4NInelasticFS.hh"
57 #include "G4ParticleHPAInelasticFS.hh" 429 #include "G4ParticleHPAInelasticFS.hh"
58 #include "G4ParticleHPD2AInelasticFS.hh" 430 #include "G4ParticleHPD2AInelasticFS.hh"
59 #include "G4ParticleHPDAInelasticFS.hh" 431 #include "G4ParticleHPDAInelasticFS.hh"
60 #include "G4ParticleHPDInelasticFS.hh" 432 #include "G4ParticleHPDInelasticFS.hh"
61 #include "G4ParticleHPHe3InelasticFS.hh" 433 #include "G4ParticleHPHe3InelasticFS.hh"
62 #include "G4ParticleHPManager.hh" <<
63 #include "G4ParticleHPN2AInelasticFS.hh" 434 #include "G4ParticleHPN2AInelasticFS.hh"
64 #include "G4ParticleHPN2PInelasticFS.hh" 435 #include "G4ParticleHPN2PInelasticFS.hh"
65 #include "G4ParticleHPN3AInelasticFS.hh" 436 #include "G4ParticleHPN3AInelasticFS.hh"
66 #include "G4ParticleHPNAInelasticFS.hh" 437 #include "G4ParticleHPNAInelasticFS.hh"
67 #include "G4ParticleHPND2AInelasticFS.hh" 438 #include "G4ParticleHPND2AInelasticFS.hh"
68 #include "G4ParticleHPNDInelasticFS.hh" 439 #include "G4ParticleHPNDInelasticFS.hh"
69 #include "G4ParticleHPNHe3InelasticFS.hh" 440 #include "G4ParticleHPNHe3InelasticFS.hh"
70 #include "G4ParticleHPNInelasticFS.hh" 441 #include "G4ParticleHPNInelasticFS.hh"
71 #include "G4ParticleHPNPAInelasticFS.hh" 442 #include "G4ParticleHPNPAInelasticFS.hh"
72 #include "G4ParticleHPNPInelasticFS.hh" 443 #include "G4ParticleHPNPInelasticFS.hh"
73 #include "G4ParticleHPNT2AInelasticFS.hh" 444 #include "G4ParticleHPNT2AInelasticFS.hh"
74 #include "G4ParticleHPNTInelasticFS.hh" 445 #include "G4ParticleHPNTInelasticFS.hh"
75 #include "G4ParticleHPNXInelasticFS.hh" 446 #include "G4ParticleHPNXInelasticFS.hh"
76 #include "G4ParticleHPPAInelasticFS.hh" 447 #include "G4ParticleHPPAInelasticFS.hh"
77 #include "G4ParticleHPPDInelasticFS.hh" 448 #include "G4ParticleHPPDInelasticFS.hh"
78 #include "G4ParticleHPPInelasticFS.hh" 449 #include "G4ParticleHPPInelasticFS.hh"
79 #include "G4ParticleHPPTInelasticFS.hh" 450 #include "G4ParticleHPPTInelasticFS.hh"
80 #include "G4ParticleHPT2AInelasticFS.hh" 451 #include "G4ParticleHPT2AInelasticFS.hh"
81 #include "G4ParticleHPTInelasticFS.hh" 452 #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 453
183 #ifdef G4VERBOSE << 454 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 455
191 G4HadFinalState* result = << 456 G4ParticleHPManager* hpmanager = G4ParticleHPManager::GetInstance();
192 (*(theInelastic[indexP]))[index]->ApplyYou <<
193 457
194 aNucleus.SetParameters(fManager->GetReaction << 458 theInelastic = hpmanager->GetInelasticFinalStates( &projectile );
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 459
207 G4ParticleHPManager::GetInstance()->CloseRea << 460 if ( G4Threading::IsMasterThread() ) {
208 461
209 return result; << 462 if ( theInelastic == NULL ) theInelastic = new std::vector<G4ParticleHPChannelList*>;
210 } <<
211 463
212 const std::pair<G4double, G4double> G4Particle << 464 if ( numEle == (G4int)G4Element::GetNumberOfElements() ) return;
213 { <<
214 // max energy non-conservation is mass of he <<
215 return std::pair<G4double, G4double>(10.0 * <<
216 } <<
217 <<
218 void G4ParticleHPInelastic::BuildPhysicsTable( <<
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 465
229 G4int nelm = (G4int)G4Element::GetNumberOfEl << 466 if ( theInelastic->size() == G4Element::GetNumberOfElements() ) {
230 G4int n0 = numEle; << 467 numEle = G4Element::GetNumberOfElements();
231 numEle = nelm; << 468 return;
232 if (!isFirst || nelm == n0) { return; } << 469 }
233 << 470 /*
234 // extra elements should be initialized << 471 const char* dataDirVariable;
235 G4AutoLock l(&theHPInelastic); << 472 if( &projectile == G4Neutron::Neutron() ) {
236 << 473 dataDirVariable = "G4NEUTRONHPDATA";
237 if (nullptr == theInelastic[indexP]) { << 474 }else if( &projectile == G4Proton::Proton() ) {
238 theInelastic[indexP] = new std::vector<G4P << 475 dataDirVariable = "G4PROTONHPDATA";
239 } << 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;
240 494
241 if (fManager->GetVerboseLevel() > 0 && isFir << 495 G4String tString = "/Inelastic";
242 fManager->DumpSetting(); << 496 dirName = dirName + tString;
243 G4cout << "@@@ G4ParticleHPInelastic insta <<
244 << theProjectile->GetParticleName() << "/ <<
245 << dirName << G4endl; <<
246 } <<
247 497
248 auto table = G4Element::GetElementTable(); << 498 */
249 for (G4int i = n0; i < nelm; ++i) { << 499 G4cout << "@@@ G4ParticleHPInelastic instantiated for particle " << projectile.GetParticleName() << " data directory variable is " << dataDirVariable << " pointing to " << dirName << G4endl;
250 auto clist = new G4ParticleHPChannelList(3 << 500 for (G4int i = numEle ; i < (G4int)G4Element::GetNumberOfElements(); i++)
251 theInelastic[indexP]->push_back(clist); << 501 {
252 clist->Init((*table)[i], dirName, theProje << 502 theInelastic->push_back( new G4ParticleHPChannelList );
253 clist->Register(new G4ParticleHPNInelastic << 503 ((*theInelastic)[i])->Init((*(G4Element::GetElementTable()))[i], dirName, const_cast<G4ParticleDefinition*>(&projectile));
254 clist->Register(new G4ParticleHPNXInelasti << 504 G4int itry = 0;
255 clist->Register(new G4ParticleHP2NDInelast << 505 do
256 clist->Register(new G4ParticleHP2NInelasti << 506 {
257 clist->Register(new G4ParticleHP3NInelasti << 507 ((*theInelastic)[i])->Register( new G4ParticleHPNInelasticFS , "F01"); // has
258 clist->Register(new G4ParticleHPNAInelasti << 508 ((*theInelastic)[i])->Register( new G4ParticleHPNXInelasticFS , "F02");
259 clist->Register(new G4ParticleHPN3AInelast << 509 ((*theInelastic)[i])->Register( new G4ParticleHP2NDInelasticFS , "F03");
260 clist->Register(new G4ParticleHP2NAInelast << 510 ((*theInelastic)[i])->Register( new G4ParticleHP2NInelasticFS , "F04"); // has, E Done
261 clist->Register(new G4ParticleHP3NAInelast << 511 ((*theInelastic)[i])->Register( new G4ParticleHP3NInelasticFS , "F05"); // has, E Done
262 clist->Register(new G4ParticleHPNPInelasti << 512 ((*theInelastic)[i])->Register( new G4ParticleHPNAInelasticFS , "F06");
263 clist->Register(new G4ParticleHPN2AInelast << 513 ((*theInelastic)[i])->Register( new G4ParticleHPN3AInelasticFS , "F07");
264 clist->Register(new G4ParticleHP2N2AInelas << 514 ((*theInelastic)[i])->Register( new G4ParticleHP2NAInelasticFS , "F08");
265 clist->Register(new G4ParticleHPNDInelasti << 515 ((*theInelastic)[i])->Register( new G4ParticleHP3NAInelasticFS , "F09");
266 clist->Register(new G4ParticleHPNTInelasti << 516 ((*theInelastic)[i])->Register( new G4ParticleHPNPInelasticFS , "F10");
267 clist->Register(new G4ParticleHPNHe3Inelas << 517 ((*theInelastic)[i])->Register( new G4ParticleHPN2AInelasticFS , "F11");
268 clist->Register(new G4ParticleHPND2AInelas << 518 ((*theInelastic)[i])->Register( new G4ParticleHP2N2AInelasticFS , "F12");
269 clist->Register(new G4ParticleHPNT2AInelas << 519 ((*theInelastic)[i])->Register( new G4ParticleHPNDInelasticFS , "F13");
270 clist->Register(new G4ParticleHP4NInelasti << 520 ((*theInelastic)[i])->Register( new G4ParticleHPNTInelasticFS , "F14");
271 clist->Register(new G4ParticleHP2NPInelast << 521 ((*theInelastic)[i])->Register( new G4ParticleHPNHe3InelasticFS , "F15");
272 clist->Register(new G4ParticleHP3NPInelast << 522 ((*theInelastic)[i])->Register( new G4ParticleHPND2AInelasticFS , "F16");
273 clist->Register(new G4ParticleHPN2PInelast << 523 ((*theInelastic)[i])->Register( new G4ParticleHPNT2AInelasticFS , "F17");
274 clist->Register(new G4ParticleHPNPAInelast << 524 ((*theInelastic)[i])->Register( new G4ParticleHP4NInelasticFS , "F18"); // has, E Done
275 clist->Register(new G4ParticleHPPInelastic << 525 ((*theInelastic)[i])->Register( new G4ParticleHP2NPInelasticFS , "F19");
276 clist->Register(new G4ParticleHPDInelastic << 526 ((*theInelastic)[i])->Register( new G4ParticleHP3NPInelasticFS , "F20");
277 clist->Register(new G4ParticleHPTInelastic << 527 ((*theInelastic)[i])->Register( new G4ParticleHPN2PInelasticFS , "F21");
278 clist->Register(new G4ParticleHPHe3Inelast << 528 ((*theInelastic)[i])->Register( new G4ParticleHPNPAInelasticFS , "F22");
279 clist->Register(new G4ParticleHPAInelastic << 529 ((*theInelastic)[i])->Register( new G4ParticleHPPInelasticFS , "F23");
280 clist->Register(new G4ParticleHP2AInelasti << 530 ((*theInelastic)[i])->Register( new G4ParticleHPDInelasticFS , "F24");
281 clist->Register(new G4ParticleHP3AInelasti << 531 ((*theInelastic)[i])->Register( new G4ParticleHPTInelasticFS , "F25");
282 clist->Register(new G4ParticleHP2PInelasti << 532 ((*theInelastic)[i])->Register( new G4ParticleHPHe3InelasticFS , "F26");
283 clist->Register(new G4ParticleHPPAInelasti << 533 ((*theInelastic)[i])->Register( new G4ParticleHPAInelasticFS , "F27");
284 clist->Register(new G4ParticleHPD2AInelast << 534 ((*theInelastic)[i])->Register( new G4ParticleHP2AInelasticFS , "F28");
285 clist->Register(new G4ParticleHPT2AInelast << 535 ((*theInelastic)[i])->Register( new G4ParticleHP3AInelasticFS , "F29");
286 clist->Register(new G4ParticleHPPDInelasti << 536 ((*theInelastic)[i])->Register( new G4ParticleHP2PInelasticFS , "F30");
287 clist->Register(new G4ParticleHPPTInelasti << 537 ((*theInelastic)[i])->Register( new G4ParticleHPPAInelasticFS , "F31");
288 clist->Register(new G4ParticleHPDAInelasti << 538 ((*theInelastic)[i])->Register( new G4ParticleHPD2AInelasticFS , "F32");
289 #ifdef G4VERBOSE << 539 ((*theInelastic)[i])->Register( new G4ParticleHPT2AInelasticFS , "F33");
290 if (fManager->GetVerboseLevel() > 1) { << 540 ((*theInelastic)[i])->Register( new G4ParticleHPPDInelasticFS , "F34");
291 G4cout << "ParticleHP::Inelastic for " << 541 ((*theInelastic)[i])->Register( new G4ParticleHPPTInelasticFS , "F35");
292 << theProjectile->GetParticleName() << << 542 ((*theInelastic)[i])->Register( new G4ParticleHPDAInelasticFS , "F36");
293 << (*table)[i]->GetName() << G4endl; << 543 ((*theInelastic)[i])->RestartRegistration();
294 } << 544 itry++;
295 #endif << 545 }
296 } << 546 while( !((*theInelastic)[i])->HasDataInAnyFinalState() && itry < 6 ); // Loop checking, 11.05.2015, T. Koi
297 fManager->RegisterInelasticFinalStates( &pro << 547 // 6 is corresponding to the value(5) of G4ParticleHPChannel. TK
298 l.unlock(); << 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();
299 } 576 }
300 577
301 void G4ParticleHPInelastic::ModelDescription(s 578 void G4ParticleHPInelastic::ModelDescription(std::ostream& outFile) const
302 { 579 {
303 outFile << "High Precision (HP) model for in << 580 outFile << "Extension of High Precision model for inelastic reaction of proton, deuteron, triton, He3 and alpha below 20MeV\n";
304 << theProjectile->GetParticleName() << " b <<
305 } 581 }
>> 582
306 583