Changeset 1315 for trunk/source/processes/hadronic/models/neutron_hp/src
- Timestamp:
- Jun 18, 2010, 11:42:07 AM (14 years ago)
- Location:
- trunk/source/processes/hadronic/models/neutron_hp/src
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/source/processes/hadronic/models/neutron_hp/src/G4NeutronHPEnAngCorrelation.cc
r962 r1315 28 28 // A prototype of the low energy neutron transport model. 29 29 // 30 // 100413 Fix bug in incidence energy by T. Koi 31 // 30 32 #include "G4NeutronHPEnAngCorrelation.hh" 31 33 #include "G4LorentzRotation.hh" 32 34 #include "G4LorentzVector.hh" 35 #include "G4RotationMatrix.hh" 33 36 34 37 G4ReactionProduct * G4NeutronHPEnAngCorrelation::SampleOne(G4double anEnergy) … … 67 70 { 68 71 // simplify and double check @ 69 G4ThreeVector the3Neutron = theNeutron.GetMomentum(); 72 G4ThreeVector the3Neutron = theNeutron.GetMomentum(); //theNeutron has value in LAB 70 73 G4double nEnergy = theNeutron.GetTotalEnergy(); 71 G4ThreeVector the3Target = theTarget.GetMomentum(); 74 G4ThreeVector the3Target = theTarget.GetMomentum(); //theTarget has value in LAB 72 75 G4double tEnergy = theTarget.GetTotalEnergy(); 73 76 G4double totE = nEnergy+tEnergy; … … 80 83 G4ReactionProduct aNeutron; 81 84 aNeutron.Lorentz(theNeutron, theCMS); 82 anEnergy = aNeutron.GetKineticEnergy(); 85 //TKDB 100413 86 //ENDF-6 Formats Manual ENDF-102 87 //CHAPTER 6. FILE 6: PRODUCT ENERGY-ANGLE DISTRIBUTIONS 88 //LCT Reference system for secondary energy and angle (incident energy is always given in the LAB system) 89 //anEnergy = aNeutron.GetKineticEnergy(); 90 anEnergy = theNeutron.GetKineticEnergy(); //should be same argumment of "anEnergy" 83 91 84 92 G4LorentzVector Ptmp (aNeutron.GetMomentum(), aNeutron.GetTotalEnergy()); 93 85 94 toZ.rotateZ(-1*Ptmp.phi()); 86 95 toZ.rotateY(-1*Ptmp.theta()); 87 88 96 } 89 97 theTotalMeanEnergy=0; 90 G4LorentzRotation toLab(toZ.inverse()); 98 G4LorentzRotation toLab(toZ.inverse()); //toLab only change axis NOT to LAB system 91 99 for(i=0; i<nProducts; i++) 92 100 { … … 110 118 it->operator[](ii)->SetMomentum(pTmp1.vect()); 111 119 it->operator[](ii)->SetTotalEnergy(pTmp1.e()); 112 if(frameFlag==1) // target rest 120 if(frameFlag==1) // target rest //TK 100413 should be LAB? 113 121 { 114 it->operator[](ii)->Lorentz(*(it->operator[](ii)), -1.*theTarget); 122 it->operator[](ii)->Lorentz(*(it->operator[](ii)), -1.*theTarget); //TK 100413 Is this really need? 115 123 } 116 124 else if(frameFlag==2) // CMS -
trunk/source/processes/hadronic/models/neutron_hp/src/G4NeutronHPInelasticCompFS.cc
r1055 r1315 38 38 // 090514 Fix bug in IC electron emission case 39 39 // Contribution from Chao Zhang (Chao.Zhang@usd.edu) and Dongming Mei(Dongming.Mei@usd.edu) 40 // 100406 "nothingWasKnownOnHadron=1" then sample mu isotropic in CM 41 // add two_body_reaction 40 42 // 41 43 #include "G4NeutronHPInelasticCompFS.hh" … … 436 438 if(nothingWasKnownOnHadron) 437 439 { 440 // TKDB 100405 441 // In this case, hadron should be isotropic in CM 442 // mu and p should be correlated 443 // 444 G4double totalPhotonEnergy = 0.0; 445 if ( thePhotons != 0 ) 446 { 447 unsigned int nPhotons = thePhotons->size(); 448 unsigned int i0; 449 for ( i0=0; i0<nPhotons; i0++) 450 { 451 //thePhotons has energies at LAB system 452 totalPhotonEnergy += thePhotons->operator[](i0)->GetTotalEnergy(); 453 } 454 } 455 456 //isotropic distribution in CM 457 G4double mu = 1.0 - 2 * G4UniformRand(); 458 459 // need momentums in target rest frame; 460 G4LorentzVector target_in_LAB ( theTarget.GetMomentum() , theTarget.GetTotalEnergy() ); 461 G4ThreeVector boostToTargetRest = -target_in_LAB.boostVector(); 462 G4LorentzVector proj_in_LAB = incidentParticle->Get4Momentum(); 463 464 G4DynamicParticle* proj = new G4DynamicParticle( G4Neutron::Neutron() , proj_in_LAB.boost( boostToTargetRest ) ); 465 G4DynamicParticle* targ = new G4DynamicParticle( G4ParticleTable::GetParticleTable()->GetIon ( (G4int)theBaseZ , (G4int)theBaseA , totalPhotonEnergy ) , G4ThreeVector(0) ); 466 G4DynamicParticle* hadron = new G4DynamicParticle( aHadron.GetDefinition() , G4ThreeVector(0) ); // will be fill momentum 467 468 two_body_reaction ( proj , targ , hadron , mu ); 469 470 G4LorentzVector hadron_in_trag_rest = hadron->Get4Momentum(); 471 G4LorentzVector hadron_in_LAB = hadron_in_trag_rest.boost ( -boostToTargetRest ); 472 aHadron.SetMomentum( hadron_in_LAB.v() ); 473 aHadron.SetKineticEnergy ( hadron_in_LAB.e() - hadron_in_LAB.m() ); 474 475 delete proj; 476 delete targ; 477 delete hadron; 478 479 //TKDB 100405 480 /* 438 481 G4double totalPhotonEnergy = 0; 439 482 if(thePhotons!=0) … … 463 506 464 507 aHadron.SetMomentum( Vector*p ); 508 */ 465 509 466 510 } … … 602 646 theResult.SetStatusChange(stopAndKill); 603 647 } 648 649 650 651 #include "G4RotationMatrix.hh" 652 void G4NeutronHPInelasticCompFS::two_body_reaction ( G4DynamicParticle* proj, G4DynamicParticle* targ, G4DynamicParticle* hadron, G4double mu ) 653 { 654 655 // Target rest flame 656 // 4vector in targ rest frame; 657 // targ could have excitation energy (photon energy will be emiited) tricky but,,, 658 659 G4LorentzVector before = proj->Get4Momentum() + targ->Get4Momentum(); 660 661 G4ThreeVector p3_proj = proj->GetMomentum(); 662 G4ThreeVector d = p3_proj.unit(); 663 G4RotationMatrix rot; 664 G4RotationMatrix rot1; 665 rot1.setPhi( pi/2 + d.phi() ); 666 G4RotationMatrix rot2; 667 rot2.setTheta( d.theta() ); 668 rot=rot2*rot1; 669 proj->SetMomentum( rot*p3_proj ); 670 671 // Now proj only has pz component; 672 673 // mu in CM system 674 675 //Valid only for neutron incidence 676 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) ); 677 678 G4double Q = proj->GetDefinition()->GetPDGMass() + targ->GetDefinition()->GetPDGMass() 679 - ( hadron->GetDefinition()->GetPDGMass() + residual->GetDefinition()->GetPDGMass() ); 680 681 // Non Relativistic Case 682 G4double A = targ->GetDefinition()->GetPDGMass() / proj->GetDefinition()->GetPDGMass(); 683 G4double AA = hadron->GetDefinition()->GetPDGMass() / proj->GetDefinition()->GetPDGMass(); 684 G4double E1 = proj->GetKineticEnergy(); 685 G4double beta = std::sqrt ( A*(A+1-AA)/AA*(1+(1+A)/A*Q/E1) ); 686 G4double gamma = AA/(A+1-AA)*beta; 687 G4double E3 = AA/std::pow((1+A),2)*(beta*beta+1+2*beta*mu)*E1; 688 G4double omega3 = (1+beta*mu)/std::sqrt(beta*beta+1+2*beta*mu); 689 690 G4double E4 = (A+1-AA)/std::pow((1+A),2)*(gamma*gamma+1-2*gamma*mu)*E1; 691 G4double omega4 = (1-gamma*mu)/std::sqrt(gamma*gamma+1-2*gamma*mu); 692 693 hadron->SetKineticEnergy ( E3 ); 694 695 G4double M = hadron->GetDefinition()->GetPDGMass(); 696 G4double pmag = std::sqrt ((E3+M)*(E3+M)-M*M) ; 697 G4ThreeVector p ( 0 , pmag*sqrt(1-omega3*omega3), pmag*omega3 ); 698 699 G4double M4 = residual->GetDefinition()->GetPDGMass(); 700 G4double pmag4 = std::sqrt ((E4+M4)*(E4+M4)-M4*M4) ; 701 G4ThreeVector p4 ( 0 , -pmag4*sqrt(1-omega4*omega4), pmag4*omega4 ); 702 703 // Rotate to orginal target rest flame. 704 p *= rot.inverse(); 705 hadron->SetMomentum( p ); 706 // Now hadron had 4 momentum in target rest flame 707 708 // TypeA 709 p4 *= rot.inverse(); 710 residual->SetMomentum ( p4 ); 711 712 //TypeB1 713 //residual->Set4Momentum ( p4_residual ); 714 //TypeB2 715 //residual->SetMomentum ( p4_residual.v() ); 716 717 // Type A make difference in Momenutum 718 // Type B1 make difference in Mass of residual 719 // Type B2 make difference in total energy. 720 721 delete residual; 722 723 }
Note: See TracChangeset
for help on using the changeset viewer.