Ignore:
Timestamp:
Jun 18, 2010, 11:42:07 AM (14 years ago)
Author:
garnier
Message:

update geant4-09-04-beta-cand-01 interfaces-V09-03-09 vis-V09-03-08

Location:
trunk/source/processes/hadronic/models/neutron_hp
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • trunk/source/processes/hadronic/models/neutron_hp/History

    r1055 r1315  
    1515     ---------------------------------------------------------------
    1616
    17 14 May 2009  Tatsumi Koi  (hadr-hpn-V09-02-00)
     1725 April 2010  Tatsumi Koi (hadr-hpn-V09-03-03)
     18- Fix compiler warning
     19        G4NeutronHPInelasticCompFS.cc
     20
     2113 April 2010  Tatsumi Koi (hadr-hpn-V09-03-02)
     22- Fix bug about incidence energy
     23        G4NeutronHPEnAngCorrelation.cc
     24
     256 April 2010  Tatsumi Koi (hadr-hpn-V09-03-01)
     26- "nothingWasKnownOnHadron=1" then sample mu isotropic in CM
     27   mu and p are correlated
     28        G4NeutronHPInelasticCompFS.hh
     29        G4NeutronHPInelasticCompFS.cc
     30
     3124 March 2010  Dennis Wright (hadr-hpn-V09-03-00)
     32-------------------------------------------------
     33- move all code in hadr-hpn-V09-02-01 to HEAD and tag as
     34  hadr-hpn-V09-03-00
     35
     3618 November 2009 Tatsumi Koi (hadr-hpn-V09-02-01)
     37-------------------------------------------------
     38- Add Ignore and Enable On Flight Doppler Broadening methods 
     39        G4NeutronHPElasticData.cc
     40        G4NeutronHPInelasticData.cc
     41        G4NeutronHPCaptureData.cc
     42
     4322 September 2009  Tatsumi Koi  (hadr-hpn-V09-02-00)
     44- Add safty for Floating Point Exception caused by 0 kinetic energy neutron
     45        G4NeutronHPCaptureData.cc
     46        G4NeutronHPElasticData.cc
     47        G4NeutronHPInelasticData.cc
     48        G4NeutronHPFissionData.cc
     49        G4NeutronHPorLCaptureData.cc
     50        G4NeutronHPorLEInelasticData.cc
     51        G4NeutronHPorLElasticData.cc
     52        G4NeutronHPorLFissionData.cc
     53
     5414 May 2009  Tatsumi Koi  (hadr-hpn-V09-01-29)
    1855- Fix bug in IC electron emissions: Contribution from Chao Zhang (Chao.Zhang@usd.edu) and Dongming Mei(Dongming.Mei@usd.edu)
    1956        G4NeutronHPPhotonDist.cc

  • trunk/source/processes/hadronic/models/neutron_hp/include/G4NeutronHPCaptureData.hh

    r1228 r1315  
    2525//
    2626//
    27 // $Id: G4NeutronHPCaptureData.hh,v 1.11 2008/04/28 19:07:54 tkoi Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     27// $Id: G4NeutronHPCaptureData.hh,v 1.12 2009/11/18 23:16:57 tkoi Exp $
     28// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    2929//
    3030// 080417 Add IsZAApplicable method (return false) by T. Koi
    3131// 080428 Add bool onFlightDB by T. Koi
     32// 091118 Add Ignore and Enable On Flight Doppler Broadening methods by T. Koi
    3233//
    3334#ifndef G4NeutronHPCaptureData_h
    … …  
    6768
    6869   void DumpPhysicsTable(const G4ParticleDefinition&);
     70
     71      void IgnoreOnFlightDopplerBroadening(){ onFlightDB = false; };
     72      void EnableOnFlightDopplerBroadening(){ onFlightDB = true; };
    6973   
    7074   private:

  • trunk/source/processes/hadronic/models/neutron_hp/include/G4NeutronHPElasticData.hh

    r1228 r1315  
    2525//
    2626//
    27 // $Id: G4NeutronHPElasticData.hh,v 1.11 2008/04/28 19:07:53 tkoi Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     27// $Id: G4NeutronHPElasticData.hh,v 1.12 2009/11/18 23:16:57 tkoi Exp $
     28// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    2929//
    3030// 080417 Add IsZAApplicable method (return false) by T. Koi
    3131// 080428 Add bool onFlightDB by T. Koi
     32// 091118 Add Ignore and Enable On Flight Doppler Broadening methods by T. Koi
    3233//
    3334#ifndef G4NeutronHPElasticData_h
    … …  
    6768
    6869   void DumpPhysicsTable(const G4ParticleDefinition&);
     70
     71      void IgnoreOnFlightDopplerBroadening(){ onFlightDB = false; };
     72      void EnableOnFlightDopplerBroadening(){ onFlightDB = true; };
    6973   
    7074   private:

  • trunk/source/processes/hadronic/models/neutron_hp/include/G4NeutronHPInelasticCompFS.hh

    r1228 r1315  
    2525//
    2626//
    27 // $Id: G4NeutronHPInelasticCompFS.hh,v 1.13 2007/06/06 12:45:13 ahoward Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     27// $Id: G4NeutronHPInelasticCompFS.hh,v 1.14 2010/04/06 19:06:11 tkoi Exp $
     28// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    2929//
    3030#ifndef G4NeutronHPInelasticCompFS_h
    … …  
    109109  G4double theCurrentA;
    110110  G4double theCurrentZ;
     111
     112   private:
     113      //                       proj                 targ                 had                  mu of had   
     114      void two_body_reaction ( G4DynamicParticle* , G4DynamicParticle* , G4DynamicParticle* , G4double mu );
     115
    111116};
    112117#endif

  • trunk/source/processes/hadronic/models/neutron_hp/include/G4NeutronHPInelasticData.hh

    r1228 r1315  
    2525//
    2626//
    27 // $Id: G4NeutronHPInelasticData.hh,v 1.11 2008/04/28 19:07:54 tkoi Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     27// $Id: G4NeutronHPInelasticData.hh,v 1.12 2009/11/18 23:16:57 tkoi Exp $
     28// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    2929//
    3030// 080417 Add IsZAApplicable method (return false) by T. Koi
    3131// 080428 Add bool onFlightDB by T. Koi
     32// 091118 Add Ignore and Enable On Flight Doppler Broadening methods by T. Koi
    3233//
    3334#ifndef G4NeutronHPInelasticData_h
    … …  
    6768
    6869   void DumpPhysicsTable(const G4ParticleDefinition&);
     70
     71      void IgnoreOnFlightDopplerBroadening(){ onFlightDB = false; };
     72      void EnableOnFlightDopplerBroadening(){ onFlightDB = true; };
    6973   
    7074   private:

  • trunk/source/processes/hadronic/models/neutron_hp/src/G4NeutronHPEnAngCorrelation.cc

    r962 r1315  
    2828// A prototype of the low energy neutron transport model.
    2929//
     30// 100413 Fix bug in incidence energy by T. Koi 
     31//
    3032#include "G4NeutronHPEnAngCorrelation.hh"
    3133#include "G4LorentzRotation.hh"
    3234#include "G4LorentzVector.hh"
     35#include "G4RotationMatrix.hh"
    3336
    3437G4ReactionProduct * G4NeutronHPEnAngCorrelation::SampleOne(G4double anEnergy)
    … …  
    6770  {
    6871    // simplify and double check @
    69     G4ThreeVector the3Neutron = theNeutron.GetMomentum();
     72    G4ThreeVector the3Neutron = theNeutron.GetMomentum(); //theNeutron has value in LAB
    7073    G4double nEnergy = theNeutron.GetTotalEnergy();
    71     G4ThreeVector the3Target = theTarget.GetMomentum();
     74    G4ThreeVector the3Target = theTarget.GetMomentum();  //theTarget has value in LAB
    7275    G4double tEnergy = theTarget.GetTotalEnergy();
    7376    G4double totE = nEnergy+tEnergy;
    … …  
    8083    G4ReactionProduct aNeutron;
    8184    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"
    8391
    8492    G4LorentzVector Ptmp (aNeutron.GetMomentum(), aNeutron.GetTotalEnergy());
     93
    8594    toZ.rotateZ(-1*Ptmp.phi());
    8695    toZ.rotateY(-1*Ptmp.theta());
    87 
    8896  }
    8997  theTotalMeanEnergy=0;
    90   G4LorentzRotation toLab(toZ.inverse());
     98  G4LorentzRotation toLab(toZ.inverse()); //toLab only change axis NOT to LAB system
    9199  for(i=0; i<nProducts; i++)
    92100  {
    … …  
    110118        it->operator[](ii)->SetMomentum(pTmp1.vect());
    111119        it->operator[](ii)->SetTotalEnergy(pTmp1.e());
    112         if(frameFlag==1) // target rest
     120        if(frameFlag==1) // target rest //TK 100413 should be LAB?
    113121        {
    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?
    115123        }
    116124        else if(frameFlag==2) // CMS

  • trunk/source/processes/hadronic/models/neutron_hp/src/G4NeutronHPInelasticCompFS.cc

    r1055 r1315  
    3838// 090514 Fix bug in IC electron emission case
    3939//        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
    4042//
    4143#include "G4NeutronHPInelasticCompFS.hh"
    … …  
    436438    if(nothingWasKnownOnHadron)
    437439    {
     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/*
    438481      G4double totalPhotonEnergy = 0;
    439482      if(thePhotons!=0)
    … …  
    463506
    464507      aHadron.SetMomentum( Vector*p );
     508*/
    465509
    466510    }
    … …  
    602646    theResult.SetStatusChange(stopAndKill);
    603647}
     648
     649
     650
     651#include "G4RotationMatrix.hh"
     652void 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.