00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046 #include "G4NeutronHPInelasticCompFS.hh"
00047 #include "G4PhysicalConstants.hh"
00048 #include "G4SystemOfUnits.hh"
00049 #include "G4Nucleus.hh"
00050 #include "G4NucleiProperties.hh"
00051 #include "G4He3.hh"
00052 #include "G4Alpha.hh"
00053 #include "G4Electron.hh"
00054 #include "G4NeutronHPDataUsed.hh"
00055 #include "G4ParticleTable.hh"
00056
00057 void G4NeutronHPInelasticCompFS::InitGammas(G4double AR, G4double ZR)
00058 {
00059
00060
00061
00062
00063
00064
00065 std::ostringstream ost;
00066 ost <<gammaPath<<"z"<<ZR<<".a"<<AR;
00067 G4String aName = ost.str();
00068 std::ifstream from(aName, std::ios::in);
00069
00070 if(!from) return;
00071
00072 std::ifstream theGammaData(aName, std::ios::in);
00073
00074 theGammas.Init(theGammaData);
00075
00076 }
00077
00078 void G4NeutronHPInelasticCompFS::Init (G4double A, G4double Z, G4int M, G4String & dirName, G4String & aFSType)
00079 {
00080 gammaPath = "/Inelastic/Gammas/";
00081 if(!getenv("G4NEUTRONHPDATA"))
00082 throw G4HadronicException(__FILE__, __LINE__, "Please setenv G4NEUTRONHPDATA to point to the neutron cross-section files.");
00083 G4String tBase = getenv("G4NEUTRONHPDATA");
00084 gammaPath = tBase+gammaPath;
00085 G4String tString = dirName;
00086 G4bool dbool;
00087 G4NeutronHPDataUsed aFile = theNames.GetName(static_cast<G4int>(A), static_cast<G4int>(Z), M, tString, aFSType, dbool);
00088 G4String filename = aFile.GetName();
00089 SetAZMs( A, Z, M, aFile );
00090
00091
00092
00093
00094
00095 if ( !dbool || ( Z<2.5 && ( std::abs(theNDLDataZ - Z)>0.0001 || std::abs(theNDLDataA - A)>0.0001)) )
00096 {
00097 if(getenv("NeutronHPNamesLogging")) G4cout << "Skipped = "<< filename <<" "<<A<<" "<<Z<<G4endl;
00098 hasAnyData = false;
00099 hasFSData = false;
00100 hasXsec = false;
00101 return;
00102 }
00103
00104
00105 std::ifstream theData(filename, std::ios::in);
00106 if(!theData)
00107 {
00108 hasAnyData = false;
00109 hasFSData = false;
00110 hasXsec = false;
00111 theData.close();
00112 return;
00113 }
00114
00115 G4int infoType, dataType, dummy;
00116 G4int sfType, it;
00117 hasFSData = false;
00118 while (theData >> infoType)
00119 {
00120 hasFSData = true;
00121 theData >> dataType;
00122 theData >> sfType >> dummy;
00123 it = 50;
00124 if(sfType>=600||(sfType<100&&sfType>=50)) it = sfType%50;
00125 if(dataType==3)
00126 {
00127
00128
00129
00130 G4double dqi;
00131 G4int ilr;
00132 theData >> dqi >> ilr;
00133
00134 QI[ it ] = dqi*eV;
00135 LR[ it ] = ilr;
00136 theXsection[it] = new G4NeutronHPVector;
00137 G4int total;
00138 theData >> total;
00139 theXsection[it]->Init(theData, total, eV);
00140
00141 }
00142 else if(dataType==4)
00143 {
00144 theAngularDistribution[it] = new G4NeutronHPAngular;
00145 theAngularDistribution[it]->Init(theData);
00146 }
00147 else if(dataType==5)
00148 {
00149 theEnergyDistribution[it] = new G4NeutronHPEnergyDistribution;
00150 theEnergyDistribution[it]->Init(theData);
00151 }
00152 else if(dataType==6)
00153 {
00154 theEnergyAngData[it] = new G4NeutronHPEnAngCorrelation;
00155 theEnergyAngData[it]->Init(theData);
00156 }
00157 else if(dataType==12)
00158 {
00159 theFinalStatePhotons[it] = new G4NeutronHPPhotonDist;
00160 theFinalStatePhotons[it]->InitMean(theData);
00161 }
00162 else if(dataType==13)
00163 {
00164 theFinalStatePhotons[it] = new G4NeutronHPPhotonDist;
00165 theFinalStatePhotons[it]->InitPartials(theData);
00166 }
00167 else if(dataType==14)
00168 {
00169 theFinalStatePhotons[it]->InitAngular(theData);
00170 }
00171 else if(dataType==15)
00172 {
00173 theFinalStatePhotons[it]->InitEnergies(theData);
00174 }
00175 else
00176 {
00177 throw G4HadronicException(__FILE__, __LINE__, "Data-type unknown to G4NeutronHPInelasticCompFS");
00178 }
00179 }
00180 theData.close();
00181 }
00182
00183 G4int G4NeutronHPInelasticCompFS::SelectExitChannel(G4double eKinetic)
00184 {
00185
00186
00187 G4double running[50];
00188 running[0] = 0;
00189 unsigned int i;
00190 for(i=0; i<50; i++)
00191 {
00192 if(i!=0) running[i]=running[i-1];
00193 if(theXsection[i] != 0)
00194 {
00195 running[i] += std::max(0., theXsection[i]->GetXsec(eKinetic));
00196 }
00197 }
00198 G4double random = G4UniformRand();
00199 G4double sum = running[49];
00200 G4int it = 50;
00201 if(0!=sum)
00202 {
00203 G4int i0;
00204 for(i0=0; i0<50; i0++)
00205 {
00206 it = i0;
00207 if(random < running[i0]/sum) break;
00208 }
00209 }
00210
00211 return it;
00212 }
00213
00214
00215
00216 void G4NeutronHPInelasticCompFS::CompositeApply(const G4HadProjectile & theTrack, G4ParticleDefinition * aDefinition)
00217 {
00218
00219
00220 theResult.Clear();
00221 G4double eKinetic = theTrack.GetKineticEnergy();
00222 const G4HadProjectile *incidentParticle = &theTrack;
00223 G4ReactionProduct theNeutron( const_cast<G4ParticleDefinition *>(incidentParticle->GetDefinition()) );
00224 theNeutron.SetMomentum( incidentParticle->Get4Momentum().vect() );
00225 theNeutron.SetKineticEnergy( eKinetic );
00226
00227
00228 G4int i;
00229 for(i=0; i<50; i++)
00230 { if(theXsection[i] != 0) { break; } }
00231
00232 G4double targetMass=0;
00233 G4double eps = 0.0001;
00234 targetMass = ( G4NucleiProperties::GetNuclearMass(static_cast<G4int>(theBaseA+eps), static_cast<G4int>(theBaseZ+eps))) /
00235 G4Neutron::Neutron()->GetPDGMass();
00236
00237
00238
00239
00240
00241
00242 G4Nucleus aNucleus;
00243 G4ReactionProduct theTarget;
00244 G4ThreeVector neuVelo = (1./incidentParticle->GetDefinition()->GetPDGMass())*theNeutron.GetMomentum();
00245 theTarget = aNucleus.GetBiasedThermalNucleus( targetMass, neuVelo, theTrack.GetMaterial()->GetTemperature());
00246
00247
00248 G4double residualMass=0;
00249 G4double residualZ = theBaseZ - aDefinition->GetPDGCharge();
00250 G4double residualA = theBaseA - aDefinition->GetBaryonNumber()+1;
00251 residualMass = ( G4NucleiProperties::GetNuclearMass(static_cast<G4int>(residualA+eps), static_cast<G4int>(residualZ+eps)) ) /
00252 G4Neutron::Neutron()->GetPDGMass();
00253
00254
00255 G4ReactionProduct boosted;
00256 boosted.Lorentz(theNeutron, theTarget);
00257 eKinetic = boosted.GetKineticEnergy();
00258
00259
00260
00261 G4int it = SelectExitChannel( eKinetic );
00262
00263
00264 InitDistributionInitialState(theNeutron, theTarget, it);
00265
00266 G4ReactionProductVector * thePhotons = 0;
00267 G4ReactionProductVector * theParticles = 0;
00268 G4ReactionProduct aHadron;
00269 aHadron.SetDefinition(aDefinition);
00270 G4double availableEnergy = theNeutron.GetKineticEnergy() + theNeutron.GetMass() - aHadron.GetMass() +
00271 (targetMass - residualMass)*G4Neutron::Neutron()->GetPDGMass();
00272
00273 if ( availableEnergy < 0 )
00274 {
00275
00276 availableEnergy = 0;
00277 }
00278 G4int nothingWasKnownOnHadron = 0;
00279 G4int dummy;
00280 G4double eGamm = 0;
00281 G4int iLevel=it-1;
00282
00283
00284 if( 50 == it )
00285 {
00286
00287
00288 iLevel=-1;
00289 aHadron.SetKineticEnergy(availableEnergy*residualMass*G4Neutron::Neutron()->GetPDGMass()/
00290 (aHadron.GetMass()+residualMass*G4Neutron::Neutron()->GetPDGMass()));
00291
00292
00293
00294
00295
00296
00297 G4double p2 = ( aHadron.GetTotalEnergy()*aHadron.GetTotalEnergy() - aHadron.GetMass()*aHadron.GetMass() );
00298 G4double p = 0.0;
00299 if ( p2 > 0.0 ) p = std::sqrt( p );
00300
00301 aHadron.SetMomentum(theNeutron.GetMomentum()*(1./theNeutron.GetTotalMomentum())*p );
00302
00303 }
00304 else
00305 {
00306 while( iLevel!=-1 && theGammas.GetLevel(iLevel) == 0 ) { iLevel--; }
00307 }
00308
00309
00310 if ( theAngularDistribution[it] != 0 )
00311 {
00312 if(theEnergyDistribution[it]!=0)
00313 {
00314
00315
00316
00317
00318
00319
00320
00321 G4double dqi = 0.0;
00322 if ( QI[it] < 0 || 849 < QI[it] ) dqi = QI[it];
00323 G4double MaxEne=eKinetic+dqi;
00324 G4double eSecN;
00325 do{
00326 eSecN=theEnergyDistribution[it]->Sample(eKinetic, dummy);
00327 }while(eSecN>MaxEne);
00328 aHadron.SetKineticEnergy(eSecN);
00329
00330 eGamm = eKinetic-eSecN;
00331 for(iLevel=theGammas.GetNumberOfLevels()-1; iLevel>=0; iLevel--)
00332 {
00333 if(theGammas.GetLevelEnergy(iLevel)<eGamm) break;
00334 }
00335 G4double random = 2*G4UniformRand();
00336 iLevel+=G4int(random);
00337 if(iLevel>theGammas.GetNumberOfLevels()-1)iLevel = theGammas.GetNumberOfLevels()-1;
00338 }
00339 else
00340 {
00341 G4double eExcitation = 0;
00342 if(iLevel>=0) eExcitation = theGammas.GetLevel(iLevel)->GetLevelEnergy();
00343 while (eKinetic-eExcitation < 0 && iLevel>0)
00344 {
00345 iLevel--;
00346 eExcitation = theGammas.GetLevel(iLevel)->GetLevelEnergy();
00347 }
00348
00349
00350 G4bool useQI=false;
00351 G4double dqi = QI[it];
00352 if ( dqi < 0 || 849 < dqi ) useQI = true;
00353
00354 if ( useQI )
00355 {
00356
00357 eExcitation = -QI[it];
00358
00359 iLevel = 0;
00360 G4bool find = false;
00361 G4int imaxEx = 0;
00362 while( !theGammas.GetLevel(iLevel+1) == 0 )
00363 {
00364 G4double maxEx = 0.0;
00365 if ( maxEx < theGammas.GetLevel(iLevel)->GetLevelEnergy() )
00366 {
00367 maxEx = theGammas.GetLevel(iLevel)->GetLevelEnergy();
00368 imaxEx = iLevel;
00369 }
00370 if ( eExcitation < theGammas.GetLevel(iLevel)->GetLevelEnergy() )
00371 {
00372 find = true;
00373 iLevel--;
00374
00375 if ( iLevel == -1 ) iLevel = 0;
00376 break;
00377 }
00378 iLevel++;
00379 }
00380
00381 if ( !find ) iLevel = imaxEx;
00382 }
00383
00384
00385 if(getenv("InelasticCompFSLogging") && eKinetic-eExcitation < 0)
00386 {
00387 throw G4HadronicException(__FILE__, __LINE__, "SEVERE: InelasticCompFS: Consistency of data not good enough, please file report");
00388 }
00389 if(eKinetic-eExcitation < 0) eExcitation = 0;
00390 if(iLevel!= -1) aHadron.SetKineticEnergy(eKinetic - eExcitation);
00391
00392 }
00393 theAngularDistribution[it]->SampleAndUpdate(aHadron);
00394
00395 if( theFinalStatePhotons[it] == 0 )
00396 {
00397
00398
00399 thePhotons = theGammas.GetDecayGammas(iLevel);
00400 eGamm -= theGammas.GetLevelEnergy(iLevel);
00401 if(eGamm>0)
00402 {
00403 G4ReactionProduct * theRestEnergy = new G4ReactionProduct;
00404 theRestEnergy->SetDefinition(G4Gamma::Gamma());
00405 theRestEnergy->SetKineticEnergy(eGamm);
00406 G4double costh = 2.*G4UniformRand()-1.;
00407 G4double phi = twopi*G4UniformRand();
00408 theRestEnergy->SetMomentum(eGamm*std::sin(std::acos(costh))*std::cos(phi),
00409 eGamm*std::sin(std::acos(costh))*std::sin(phi),
00410 eGamm*costh);
00411 if(thePhotons == 0) { thePhotons = new G4ReactionProductVector; }
00412 thePhotons->push_back(theRestEnergy);
00413 }
00414 }
00415 }
00416 else if(theEnergyAngData[it] != 0)
00417 {
00418 theParticles = theEnergyAngData[it]->Sample(eKinetic);
00419 }
00420 else
00421 {
00422
00423 nothingWasKnownOnHadron = 1;
00424 }
00425
00426
00427
00428
00429
00430
00431
00432 if ( theFinalStatePhotons[it] != 0 )
00433 {
00434
00435
00436 G4ReactionProduct boosted_tmp;
00437 boosted_tmp.Lorentz(theNeutron, theTarget);
00438 G4double anEnergy = boosted_tmp.GetKineticEnergy();
00439 thePhotons = theFinalStatePhotons[it]->GetPhotons(anEnergy);
00440 G4double aBaseEnergy = theFinalStatePhotons[it]->GetLevelEnergy();
00441 G4double testEnergy = 0;
00442 if(thePhotons!=0 && thePhotons->size()!=0)
00443 { aBaseEnergy-=thePhotons->operator[](0)->GetTotalEnergy(); }
00444 if(theFinalStatePhotons[it]->NeedsCascade())
00445 {
00446 while(aBaseEnergy>0.01*keV)
00447 {
00448
00449 G4bool foundMatchingLevel = false;
00450 G4int closest = 2;
00451 G4double deltaEold = -1;
00452 for(G4int j=1; j<it; j++)
00453 {
00454 if(theFinalStatePhotons[j]!=0)
00455 {
00456 testEnergy = theFinalStatePhotons[j]->GetLevelEnergy();
00457 }
00458 else
00459 {
00460 testEnergy = 0;
00461 }
00462 G4double deltaE = std::abs(testEnergy-aBaseEnergy);
00463 if(deltaE<0.1*keV)
00464 {
00465 G4ReactionProductVector * theNext =
00466 theFinalStatePhotons[j]->GetPhotons(anEnergy);
00467 thePhotons->push_back(theNext->operator[](0));
00468 aBaseEnergy = testEnergy-theNext->operator[](0)->GetTotalEnergy();
00469 delete theNext;
00470 foundMatchingLevel = true;
00471 break;
00472 }
00473 if(theFinalStatePhotons[j]!=0 && ( deltaE<deltaEold||deltaEold<0.) )
00474 {
00475 closest = j;
00476 deltaEold = deltaE;
00477 }
00478 }
00479 if(!foundMatchingLevel)
00480 {
00481 G4ReactionProductVector * theNext =
00482 theFinalStatePhotons[closest]->GetPhotons(anEnergy);
00483 thePhotons->push_back(theNext->operator[](0));
00484 aBaseEnergy = aBaseEnergy-theNext->operator[](0)->GetTotalEnergy();
00485 delete theNext;
00486 }
00487 }
00488 }
00489 }
00490 unsigned int i0;
00491 if(thePhotons!=0)
00492 {
00493 for(i0=0; i0<thePhotons->size(); i0++)
00494 {
00495
00496 thePhotons->operator[](i0)->Lorentz(*(thePhotons->operator[](i0)), -1.*theTarget);
00497 }
00498 }
00499
00500 if(nothingWasKnownOnHadron)
00501 {
00502
00503
00504
00505
00506 G4double totalPhotonEnergy = 0.0;
00507 if ( thePhotons != 0 )
00508 {
00509 unsigned int nPhotons = thePhotons->size();
00510 unsigned int ii0;
00511 for ( ii0=0; ii0<nPhotons; ii0++)
00512 {
00513
00514 totalPhotonEnergy += thePhotons->operator[](ii0)->GetTotalEnergy();
00515 }
00516 }
00517
00518
00519 G4double mu = 1.0 - 2 * G4UniformRand();
00520
00521
00522 G4LorentzVector target_in_LAB ( theTarget.GetMomentum() , theTarget.GetTotalEnergy() );
00523 G4ThreeVector boostToTargetRest = -target_in_LAB.boostVector();
00524 G4LorentzVector proj_in_LAB = incidentParticle->Get4Momentum();
00525
00526 G4DynamicParticle* proj = new G4DynamicParticle( G4Neutron::Neutron() , proj_in_LAB.boost( boostToTargetRest ) );
00527 G4DynamicParticle* targ = new G4DynamicParticle( G4ParticleTable::GetParticleTable()->GetIon ( (G4int)theBaseZ , (G4int)theBaseA , totalPhotonEnergy ) , G4ThreeVector(0) );
00528 G4DynamicParticle* hadron = new G4DynamicParticle( aHadron.GetDefinition() , G4ThreeVector(0) );
00529
00530 two_body_reaction ( proj , targ , hadron , mu );
00531
00532 G4LorentzVector hadron_in_trag_rest = hadron->Get4Momentum();
00533 G4LorentzVector hadron_in_LAB = hadron_in_trag_rest.boost ( -boostToTargetRest );
00534 aHadron.SetMomentum( hadron_in_LAB.v() );
00535 aHadron.SetKineticEnergy ( hadron_in_LAB.e() - hadron_in_LAB.m() );
00536
00537 delete proj;
00538 delete targ;
00539 delete hadron;
00540
00541
00542
00543
00544
00545
00546
00547
00548
00549
00550
00551
00552
00553
00554
00555
00556
00557
00558
00559
00560
00561
00562
00563
00564
00565
00566
00567
00568
00569
00570
00571
00572 }
00573
00574
00575
00576
00577
00578
00579
00580
00581
00582 G4int nSecondaries = 2;
00583 G4bool needsSeparateRecoil = false;
00584 G4int totalBaryonNumber = 0;
00585 G4int totalCharge = 0;
00586 G4ThreeVector totalMomentum(0);
00587 if(theParticles != 0)
00588 {
00589 nSecondaries = theParticles->size();
00590 G4ParticleDefinition * aDef;
00591 unsigned int ii0;
00592 for(ii0=0; ii0<theParticles->size(); ii0++)
00593 {
00594 aDef = theParticles->operator[](ii0)->GetDefinition();
00595 totalBaryonNumber+=aDef->GetBaryonNumber();
00596 totalCharge+=G4int(aDef->GetPDGCharge()+eps);
00597 totalMomentum += theParticles->operator[](ii0)->GetMomentum();
00598 }
00599 if(totalBaryonNumber!=G4int(theBaseA+eps+incidentParticle->GetDefinition()->GetBaryonNumber()))
00600 {
00601 needsSeparateRecoil = true;
00602 nSecondaries++;
00603 residualA = G4int(theBaseA+eps+incidentParticle->GetDefinition()->GetBaryonNumber()
00604 -totalBaryonNumber);
00605 residualZ = G4int(theBaseZ+eps+incidentParticle->GetDefinition()->GetPDGCharge()
00606 -totalCharge);
00607 }
00608 }
00609
00610 G4int nPhotons = 0;
00611 if(thePhotons!=0) { nPhotons = thePhotons->size(); }
00612 nSecondaries += nPhotons;
00613
00614 G4DynamicParticle * theSec;
00615
00616 if( theParticles==0 )
00617 {
00618 theSec = new G4DynamicParticle;
00619 theSec->SetDefinition(aHadron.GetDefinition());
00620 theSec->SetMomentum(aHadron.GetMomentum());
00621 theResult.AddSecondary(theSec);
00622
00623 aHadron.Lorentz(aHadron, theTarget);
00624 G4ReactionProduct theResidual;
00625 theResidual.SetDefinition(G4ParticleTable::GetParticleTable()
00626 ->GetIon(static_cast<G4int>(residualZ), static_cast<G4int>(residualA), 0));
00627 theResidual.SetKineticEnergy(aHadron.GetKineticEnergy()*aHadron.GetMass()/theResidual.GetMass());
00628
00629
00630
00631 G4ThreeVector incidentNeutronMomentum = theNeutron.GetMomentum();
00632 theResidual.SetMomentum(incidentNeutronMomentum - aHadron.GetMomentum());
00633
00634 theResidual.Lorentz(theResidual, -1.*theTarget);
00635 G4ThreeVector totalPhotonMomentum(0,0,0);
00636 if(thePhotons!=0)
00637 {
00638 for(i=0; i<nPhotons; i++)
00639 {
00640 totalPhotonMomentum += thePhotons->operator[](i)->GetMomentum();
00641 }
00642 }
00643 theSec = new G4DynamicParticle;
00644 theSec->SetDefinition(theResidual.GetDefinition());
00645 theSec->SetMomentum(theResidual.GetMomentum()-totalPhotonMomentum);
00646 theResult.AddSecondary(theSec);
00647 }
00648 else
00649 {
00650 for(i0=0; i0<theParticles->size(); i0++)
00651 {
00652 theSec = new G4DynamicParticle;
00653 theSec->SetDefinition(theParticles->operator[](i0)->GetDefinition());
00654 theSec->SetMomentum(theParticles->operator[](i0)->GetMomentum());
00655 theResult.AddSecondary(theSec);
00656 delete theParticles->operator[](i0);
00657 }
00658 delete theParticles;
00659 if(needsSeparateRecoil && residualZ!=0)
00660 {
00661 G4ReactionProduct theResidual;
00662 theResidual.SetDefinition(G4ParticleTable::GetParticleTable()
00663 ->GetIon(static_cast<G4int>(residualZ), static_cast<G4int>(residualA), 0));
00664 G4double resiualKineticEnergy = theResidual.GetMass()*theResidual.GetMass();
00665 resiualKineticEnergy += totalMomentum*totalMomentum;
00666 resiualKineticEnergy = std::sqrt(resiualKineticEnergy) - theResidual.GetMass();
00667
00668 theResidual.SetKineticEnergy(resiualKineticEnergy);
00669
00670
00671
00672
00673
00674
00675 theResidual.SetMomentum( theNeutron.GetMomentum() + theTarget.GetMomentum() - totalMomentum );
00676
00677 theSec = new G4DynamicParticle;
00678 theSec->SetDefinition(theResidual.GetDefinition());
00679 theSec->SetMomentum(theResidual.GetMomentum());
00680 theResult.AddSecondary(theSec);
00681 }
00682 }
00683 if(thePhotons!=0)
00684 {
00685 for(i=0; i<nPhotons; i++)
00686 {
00687 theSec = new G4DynamicParticle;
00688
00689
00690 theSec->SetDefinition( thePhotons->operator[](i)->GetDefinition() );
00691
00692 theSec->SetMomentum(thePhotons->operator[](i)->GetMomentum());
00693 theResult.AddSecondary(theSec);
00694 delete thePhotons->operator[](i);
00695 }
00696
00697 delete thePhotons;
00698 }
00699
00700
00701 G4ParticleDefinition* targ_pd = G4ParticleTable::GetParticleTable()->GetIon ( (G4int)theBaseZ , (G4int)theBaseA , 0.0 );
00702 G4LorentzVector targ_4p_lab ( theTarget.GetMomentum() , std::sqrt( targ_pd->GetPDGMass()*targ_pd->GetPDGMass() + theTarget.GetMomentum().mag2() ) );
00703 G4LorentzVector proj_4p_lab = theTrack.Get4Momentum();
00704 G4LorentzVector init_4p_lab = proj_4p_lab + targ_4p_lab;
00705 adjust_final_state ( init_4p_lab );
00706
00707
00708 theResult.SetStatusChange(stopAndKill);
00709 }
00710
00711
00712
00713 #include "G4RotationMatrix.hh"
00714 void G4NeutronHPInelasticCompFS::two_body_reaction ( G4DynamicParticle* proj, G4DynamicParticle* targ, G4DynamicParticle* hadron, G4double mu )
00715 {
00716
00717
00718
00719
00720
00721 G4LorentzVector before = proj->Get4Momentum() + targ->Get4Momentum();
00722
00723 G4ThreeVector p3_proj = proj->GetMomentum();
00724 G4ThreeVector d = p3_proj.unit();
00725 G4RotationMatrix rot;
00726 G4RotationMatrix rot1;
00727 rot1.setPhi( pi/2 + d.phi() );
00728 G4RotationMatrix rot2;
00729 rot2.setTheta( d.theta() );
00730 rot=rot2*rot1;
00731 proj->SetMomentum( rot*p3_proj );
00732
00733
00734
00735
00736
00737
00738 G4DynamicParticle* residual = new G4DynamicParticle ( G4ParticleTable::GetParticleTable()->GetIon ( (G4int)( targ->GetDefinition()->GetPDGCharge() - hadron->GetDefinition()->GetPDGCharge() ) , (G4int)(targ->GetDefinition()->GetBaryonNumber() - hadron->GetDefinition()->GetBaryonNumber()+1) , 0 ) , G4ThreeVector(0) );
00739
00740 G4double Q = proj->GetDefinition()->GetPDGMass() + targ->GetDefinition()->GetPDGMass()
00741 - ( hadron->GetDefinition()->GetPDGMass() + residual->GetDefinition()->GetPDGMass() );
00742
00743
00744 G4double A = targ->GetDefinition()->GetPDGMass() / proj->GetDefinition()->GetPDGMass();
00745 G4double AA = hadron->GetDefinition()->GetPDGMass() / proj->GetDefinition()->GetPDGMass();
00746 G4double E1 = proj->GetKineticEnergy();
00747
00748
00749
00750
00751 if ( ( 1 + (1+A)/A*Q/E1 ) < 0 )
00752 {
00753
00754 Q = -( A/(1+A)*E1 ) + 1.0e-6*eV;
00755 }
00756
00757 G4double beta = std::sqrt ( A*(A+1-AA)/AA*( 1 + (1+A)/A*Q/E1 ) );
00758 G4double gamma = AA/(A+1-AA)*beta;
00759 G4double E3 = AA/std::pow((1+A),2)*(beta*beta+1+2*beta*mu)*E1;
00760 G4double omega3 = (1+beta*mu)/std::sqrt(beta*beta+1+2*beta*mu);
00761 if ( omega3 > 1.0 ) omega3 = 1.0;
00762
00763 G4double E4 = (A+1-AA)/std::pow((1+A),2)*(gamma*gamma+1-2*gamma*mu)*E1;
00764 G4double omega4 = (1-gamma*mu)/std::sqrt(gamma*gamma+1-2*gamma*mu);
00765 if ( omega4 > 1.0 ) omega4 = 1.0;
00766
00767 hadron->SetKineticEnergy ( E3 );
00768
00769 G4double M = hadron->GetDefinition()->GetPDGMass();
00770 G4double pmag = std::sqrt ((E3+M)*(E3+M)-M*M) ;
00771 G4ThreeVector p ( 0 , pmag*std::sqrt(1-omega3*omega3), pmag*omega3 );
00772
00773 G4double M4 = residual->GetDefinition()->GetPDGMass();
00774 G4double pmag4 = std::sqrt ((E4+M4)*(E4+M4)-M4*M4) ;
00775 G4ThreeVector p4 ( 0 , -pmag4*std::sqrt(1-omega4*omega4), pmag4*omega4 );
00776
00777
00778 p *= rot.inverse();
00779 hadron->SetMomentum( p );
00780
00781
00782
00783 p4 *= rot.inverse();
00784 residual->SetMomentum ( p4 );
00785
00786
00787
00788
00789
00790
00791
00792
00793
00794
00795 delete residual;
00796
00797 }