Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/lend/src/G4LENDElastic.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  1 //
  2 // ********************************************************************
  3 // * License and Disclaimer                                           *
  4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.                             *
 10 // *                                                                  *
 11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                                                  *
 18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // ********************************************************************
 25 //
 26 
 27 #include "G4LENDElastic.hh"
 28 #include "G4Pow.hh"
 29 #include "G4PhysicalConstants.hh"
 30 #include "G4SystemOfUnits.hh"
 31 #include "G4Nucleus.hh"
 32 #include "G4IonTable.hh"
 33 
 34 //extern "C" double MyRNG(void*) { return drand48(); }
 35 //extern "C" double MyRNG(void*) { return  CLHEP::HepRandom::getTheEngine()->flat(); }
 36 
 37 G4HadFinalState * G4LENDElastic::ApplyYourself(const G4HadProjectile& aTrack, G4Nucleus& aTarg )
 38 {
 39 
 40    G4double temp = aTrack.GetMaterial()->GetTemperature();
 41 
 42    //G4int iZ = int ( aTarg.GetZ() );
 43    //G4int iA = int ( aTarg.GetN() );
 44    //migrate to integer A and Z (GetN_asInt returns number of neutrons in the nucleus since this) 
 45    G4int iZ = aTarg.GetZ_asInt();
 46    G4int iA = aTarg.GetA_asInt();
 47    G4int iM = 0;
 48    if ( aTarg.GetIsotope() != NULL ) {
 49       iM = aTarg.GetIsotope()->Getm();
 50    }
 51 
 52    G4double ke = aTrack.GetKineticEnergy();
 53 
 54    G4HadFinalState* theResult = &theParticleChange;
 55    theResult->Clear();
 56 
 57    G4GIDI_target* aTarget = get_target_from_map( lend_manager->GetNucleusEncoding( iZ , iA , iM ) );
 58    if ( aTarget == NULL ) return returnUnchanged( aTrack , theResult );
 59 
 60    G4double aMu = aTarget->getElasticFinalState( ke*MeV, temp, MyRNG , NULL );
 61 
 62    G4double phi = twopi*G4UniformRand();
 63    G4double theta = std::acos( aMu );
 64    //G4double sinth = std::sin( theta );
 65 
 66    G4ReactionProduct theNeutron( aTrack.GetDefinition() );
 67    theNeutron.SetMomentum( aTrack.Get4Momentum().vect() );
 68    theNeutron.SetKineticEnergy( ke );
 69 
 70    //G4ParticleDefinition* pd = G4IonTable::GetIonTable()->GetIon( iZ , iA , iM );
 71    //TK 170509 Fix for the case of excited isomer target 
 72    G4double EE = 0.0;
 73    if ( iM != 0 ) {
 74       G4LENDManager::GetInstance()->GetExcitationEnergyOfExcitedIsomer( iZ , iA , iM );
 75    }
 76    G4ParticleDefinition* target_pd = G4IonTable::GetIonTable()->GetIon( iZ , iA , EE );
 77    G4ReactionProduct theTarget( target_pd );
 78 
 79    G4double mass = target_pd->GetPDGMass();
 80 
 81 // add Thermal motion 
 82    G4double kT = k_Boltzmann*temp;
 83    G4ThreeVector v ( G4RandGauss::shoot() * std::sqrt( kT*mass ) 
 84                    , G4RandGauss::shoot() * std::sqrt( kT*mass ) 
 85                    , G4RandGauss::shoot() * std::sqrt( kT*mass ) );
 86    theTarget.SetMomentum( v );
 87 
 88      G4ThreeVector the3Neutron = theNeutron.GetMomentum();
 89      G4double nEnergy = theNeutron.GetTotalEnergy();
 90      G4ThreeVector the3Target = theTarget.GetMomentum();
 91      G4double tEnergy = theTarget.GetTotalEnergy();
 92      G4ReactionProduct theCMS;
 93      G4double totE = nEnergy+tEnergy;
 94      G4ThreeVector the3CMS = the3Target+the3Neutron;
 95      theCMS.SetMomentum(the3CMS);
 96      G4double cmsMom = std::sqrt(the3CMS*the3CMS);
 97      G4double sqrts = std::sqrt((totE-cmsMom)*(totE+cmsMom));
 98      theCMS.SetMass(sqrts);
 99      theCMS.SetTotalEnergy(totE);
100 
101        theNeutron.Lorentz(theNeutron, theCMS);
102        theTarget.Lorentz(theTarget, theCMS);
103        G4double en = theNeutron.GetTotalMomentum(); // already in CMS.
104        G4ThreeVector cms3Mom=theNeutron.GetMomentum(); // for neutron direction in CMS
105        G4double cms_theta=cms3Mom.theta();
106        G4double cms_phi=cms3Mom.phi();
107        G4ThreeVector tempVector;
108        tempVector.setX( std::cos(theta)*std::sin(cms_theta)*std::cos(cms_phi)
109                        +std::sin(theta)*std::cos(phi)*std::cos(cms_theta)*std::cos(cms_phi)
110                        -std::sin(theta)*std::sin(phi)*std::sin(cms_phi) );
111        tempVector.setY( std::cos(theta)*std::sin(cms_theta)*std::sin(cms_phi)
112                        +std::sin(theta)*std::cos(phi)*std::cos(cms_theta)*std::sin(cms_phi)
113                        +std::sin(theta)*std::sin(phi)*std::cos(cms_phi) );
114        tempVector.setZ( std::cos(theta)*std::cos(cms_theta)
115                        -std::sin(theta)*std::cos(phi)*std::sin(cms_theta) );
116        tempVector *= en;
117        theNeutron.SetMomentum(tempVector);
118        theTarget.SetMomentum(-tempVector);
119        G4double tP = theTarget.GetTotalMomentum();
120        G4double tM = theTarget.GetMass();
121        theTarget.SetTotalEnergy(std::sqrt((tP+tM)*(tP+tM)-2.*tP*tM));
122 
123 
124       theNeutron.Lorentz(theNeutron, -1.*theCMS);
125       theTarget.Lorentz(theTarget, -1.*theCMS);
126 
127 //110913 Add Protection for very low energy (1e-6eV) scattering 
128       if ( theNeutron.GetKineticEnergy() <= 0 )
129       {
130          theNeutron.SetTotalEnergy ( theNeutron.GetMass() * ( 1 + G4Pow::GetInstance()->powA( 10 , -15.65 ) ) );
131       }
132 
133       if ( theTarget.GetKineticEnergy() < 0 )
134       {
135          theTarget.SetTotalEnergy ( theTarget.GetMass() * ( 1 + G4Pow::GetInstance()->powA( 10 , -15.65 ) ) );
136       }
137 //110913 END
138 
139      theResult->SetEnergyChange(theNeutron.GetKineticEnergy());
140      theResult->SetMomentumChange(theNeutron.GetMomentum().unit());
141      G4DynamicParticle* theRecoil = new G4DynamicParticle;
142 
143      theRecoil->SetDefinition( target_pd );
144      theRecoil->SetMomentum( theTarget.GetMomentum() );
145      theResult->AddSecondary( theRecoil, secID );
146 
147    return theResult; 
148 
149 }
150 
151