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