source: trunk/source/processes/hadronic/models/neutron_hp/src/G4NeutronHPJENDLHEData.cc @ 1358

//
// ********************************************************************
// * License and Disclaimer                                           *
// *                                                                  *
// * The  Geant4 software  is  copyright of the Copyright Holders  of *
// * the Geant4 Collaboration.  It is provided  under  the terms  and *
// * conditions of the Geant4 Software License,  included in the file *
// * LICENSE and available at  http://cern.ch/geant4/license .  These *
// * include a list of copyright holders.                             *
// *                                                                  *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work  make  any representation or  warranty, express or implied, *
// * regarding  this  software system or assume any liability for its *
// * use.  Please see the license in the file  LICENSE  and URL above *
// * for the full disclaimer and the limitation of liability.         *
// *                                                                  *
// * This  code  implementation is the result of  the  scientific and *
// * technical work of the GEANT4 collaboration.                      *
// * By using,  copying,  modifying or  distributing the software (or *
// * any work based  on the software)  you  agree  to acknowledge its *
// * use  in  resulting  scientific  publications,  and indicate your *
// * acceptance of all terms of the Geant4 Software license.          *
// ********************************************************************
//
// Class Description
// Cross-section data set for a high precision (based on JENDL_HE evaluated data
// libraries) description of elastic scattering 20 MeV ~ 3 GeV;
// Class Description - End

// 15-Nov-06 First Implementation is done by T. Koi (SLAC/SCCS)

#include "G4NeutronHPJENDLHEData.hh"
#include "G4LPhysicsFreeVector.hh"
#include "G4ElementTable.hh"
#include "G4NeutronHPData.hh"

G4bool G4NeutronHPJENDLHEData::IsApplicable(const G4DynamicParticle*aP, const G4Element* anE)
{

   G4bool result = true;
   G4double eKin = aP->GetKineticEnergy();
   //if(eKin>20*MeV||aP->GetDefinition()!=G4Neutron::Neutron()) result = false;
   if ( eKin < 20*MeV || 3*GeV < eKin || aP->GetDefinition()!=G4Neutron::Neutron() ) 
   {
      result = false;
   } 
// Element Check 
   else if ( !(vElement[ anE->GetIndex() ]) ) result = false;

   return result;

}



G4NeutronHPJENDLHEData::G4NeutronHPJENDLHEData()
{
   ;
}
   


G4NeutronHPJENDLHEData::G4NeutronHPJENDLHEData( G4String reaction , G4ParticleDefinition* pd )
{
   reactionName = reaction;
   BuildPhysicsTable( *pd );
}



G4NeutronHPJENDLHEData::~G4NeutronHPJENDLHEData()
{
   ; 
   //delete theCrossSections;
}
 


void G4NeutronHPJENDLHEData::BuildPhysicsTable( const G4ParticleDefinition& aP )
{

//   if ( &aP != G4Neutron::Neutron() ) 
//      throw G4HadronicException(__FILE__, __LINE__, "Attempt to use NeutronHP data for particles other than neutrons!!!");  
   particleName = aP.GetParticleName();

   G4String baseName = getenv( "G4NEUTRONHPDATA" );
   G4String dirName = baseName+"/JENDL_HE/"+particleName+"/"+reactionName ;
   G4String aFSType = "/CrossSection/";
   G4NeutronHPNames theNames; 

   G4String filename;

// Create JENDL_HE data 
// Create map element or isotope  

   size_t numberOfElements = G4Element::GetNumberOfElements();
   //theCrossSections = new G4PhysicsTable( numberOfElements );

   // make a PhysicsVector for each element

   static const G4ElementTable *theElementTable = G4Element::GetElementTable();
   vElement.clear();
   vElement.resize( numberOfElements );
   for ( size_t i = 0; i < numberOfElements; ++i )
   {

      G4Element* theElement = (*theElementTable)[i];
      vElement[i] = false;

      // isotope
      G4int nIso = (*theElementTable)[i]->GetNumberOfIsotopes();
      G4int Z = static_cast<G4int> ((*theElementTable)[i]->GetZ());
      if ( nIso!=0 )
      {
         G4bool found_at_least_one = false; 
         for ( G4int i1 = 0; i1 < nIso; i1++ )
         {
             G4int A = theElement->GetIsotope(i1)->GetN();

             if ( isThisNewIsotope( Z , A ) ) 
             { 

                std::stringstream ss; 
                ss << dirName << aFSType << Z << "_" << A << "_" << theNames.GetName( Z-1 );
                filename = ss.str();
                std::fstream file;
                file.open ( filename , std::fstream::in );
                G4int dummy;
                file >> dummy;
                if ( file.good() ) 
                {

                   //G4cout << "Found file for Z=" << Z << ", A=" << A << ", as " << filename << G4endl;
                   found_at_least_one = true;

                   // read the file
                   G4PhysicsVector* aPhysVec = readAFile ( &file );

                   //Regist 

                   registAPhysicsVector( Z , A , aPhysVec );

                }
                else 
                { 
                   //G4cout << "No file for "<< reactionType << " Z=" << Z << ", A=" << A << G4endl;
                }

                file.close();

             }
             else
             {
                found_at_least_one = TRUE;
             }
          }

          if ( found_at_least_one ) vElement[i] = true;

       }
       else
       {
          G4StableIsotopes theStableOnes;
          G4int first = theStableOnes.GetFirstIsotope( Z );
          G4bool found_at_least_one = FALSE; 
          for ( G4int i1 = 0; i1 < theStableOnes.GetNumberOfIsotopes( static_cast<G4int>(theElement->GetZ() ) ); i1++)
          {
             G4int A = theStableOnes.GetIsotopeNucleonCount( first+i1 );

             if ( isThisNewIsotope( Z , A ) ) 
             {

                std::stringstream ss; 
                ss << dirName << aFSType << Z << "_" << A << "_" << theNames.GetName( Z-1 );
                filename = ss.str();

                std::fstream file;
                file.open ( filename , std::fstream::in );
                G4int dummy;
                file >> dummy;
                if ( file.good() ) 
                {
                   //G4cout << "Found file for Z=" << Z << ", A=" << A << ", as " << filename << G4endl;
                   found_at_least_one = TRUE;
                   //Read the file

                   G4PhysicsVector* aPhysVec = readAFile ( &file );

                   //Regist the PhysicsVector
                   registAPhysicsVector( Z , A , aPhysVec );

                }
                else 
                { 
                   //G4cout << "No file for "<< reactionType << " Z=" << Z << ", A=" << A << G4endl;
                }

                file.close();
             }
             else
             {
                found_at_least_one = TRUE;
             }
          }

          if ( found_at_least_one ) vElement[i] = true;

       }

   }

}



void G4NeutronHPJENDLHEData::DumpPhysicsTable(const G4ParticleDefinition& aP)
{
  if(&aP!=G4Neutron::Neutron()) 
     throw G4HadronicException(__FILE__, __LINE__, "Attempt to use NeutronHP data for particles other than neutrons!!!");  
//  G4cout << "G4NeutronHPJENDLHEData::DumpPhysicsTable still to be implemented"<<G4endl;
}



G4double G4NeutronHPJENDLHEData::
GetCrossSection(const G4DynamicParticle* aP, const G4Element*anE, G4double )
//                                                                    aTemp  
{

   // Primary energy >20MeV
   // Thus
   // Not take account of Doppler broadening 
   // also
   // Not take account of Target thermal motions

   G4double result = 0;

   G4double ek = aP->GetKineticEnergy();

   G4int nIso = anE->GetNumberOfIsotopes();
   G4int Z = static_cast<G4int> ( anE->GetZ() );
   if ( nIso!=0 )
   {
      for ( G4int i1 = 0; i1 < nIso; i1++ )
      {

         G4int A = anE->GetIsotope(i1)->GetN();
         G4double frac = anE->GetRelativeAbundanceVector()[ i1 ];   // This case do NOT request "*perCent".    

         result += frac * getXSfromThisIsotope( Z , A , ek );
         //G4cout << reactionType << " XS in barn " << Z << " " << A << " " << frac << " " << getXSfromThisIsotope( Z , A , ek )/barn << G4endl; 

      }
   }
   else
   {

      G4StableIsotopes theStableOnes;
      G4int first = theStableOnes.GetFirstIsotope( Z );
      for ( G4int i1 = 0; i1 < theStableOnes.GetNumberOfIsotopes( static_cast<G4int>(anE->GetZ() ) ); i1++)
      {
         
         G4int A = theStableOnes.GetIsotopeNucleonCount( first+i1 );
         G4double frac = theStableOnes.GetAbundance( first+i1 )*perCent;  // This case request "*perCent". 

         result += frac * getXSfromThisIsotope( Z , A , ek );
         //G4cout << reactionType << " XS in barn " << Z << " " << A << " " << frac << " " << getXSfromThisIsotope( Z , A , ek )/barn << G4endl; 
          
      }
   }
   return result;

}



G4PhysicsVector* G4NeutronHPJENDLHEData::readAFile ( std::fstream* file )
{

   G4int dummy;
   G4int len;
   *file >> dummy;
   *file >> len;

   std::vector< G4double > v_e; 
   std::vector< G4double > v_xs; 

   for ( G4int i = 0 ; i < len ; i++ )
   {
      G4double e;
      G4double xs;

      *file >> e; 
      *file >> xs;
      // data are written in eV and barn.    
      v_e.push_back( e*eV );
      v_xs.push_back( xs*barn );
   }

   G4LPhysicsFreeVector* aPhysVec = new G4LPhysicsFreeVector( static_cast< size_t >( len ) , v_e.front() , v_e.back() );

   for ( G4int i = 0 ; i < len ; i++ )
   {
      aPhysVec->PutValues( static_cast< size_t >( i ) , v_e[ i ] , v_xs[ i ] );   
   }

   return aPhysVec;
}



G4bool G4NeutronHPJENDLHEData::isThisInMap( G4int z , G4int a )
{
   if ( mIsotope.find ( z ) == mIsotope.end() ) return false;
   if ( mIsotope.find ( z ) -> second->find ( a ) ==  mIsotope.find ( z ) -> second->end() ) return false;
   return true; 
}



void G4NeutronHPJENDLHEData::registAPhysicsVector( G4int Z , G4int A , G4PhysicsVector* aPhysVec )
{

    std::pair< G4int , G4PhysicsVector* > aPair = std::pair < G4int , G4PhysicsVector* > ( A , aPhysVec );  

    std::map < G4int , std::map< G4int , G4PhysicsVector* >* >::iterator itm; 
    itm = mIsotope.find ( Z );
    if ( itm !=  mIsotope.end() ) 
    { 
       itm->second->insert ( aPair ); 
    }  
    else
    {
       std::map< G4int , G4PhysicsVector* >* aMap = new std::map< G4int , G4PhysicsVector* >;
       aMap->insert ( aPair ); 
       mIsotope.insert( std::pair< G4int , std::map< G4int , G4PhysicsVector* >* > ( Z , aMap ) );
    }

}



G4double G4NeutronHPJENDLHEData::getXSfromThisIsotope( G4int Z , G4int A , G4double ek )
{

   G4double aXSection = 0.0;
   G4bool outOfRange;

   G4PhysicsVector* aPhysVec;
   if ( mIsotope.find ( Z )->second->find ( A ) != mIsotope.find ( Z )->second->end() )
   {
  
      aPhysVec = mIsotope.find ( Z )->second->find ( A )->second;
      aXSection = aPhysVec->GetValue( ek , outOfRange );

   } 
   else
   { 

      //Select closest one in the same Z
      std::map < G4int , G4PhysicsVector* >::iterator it; 
      G4int delta0 = 99; // no mean for 99 
      for ( it = mIsotope.find ( Z )->second->begin() ; it != mIsotope.find ( Z )->second->end() ; it++ )
      {
         G4int delta = std::abs( A - it->first );
         if ( delta < delta0 ) delta0 = delta;
      }

      // Randomize of selection larger or smaller than A
      if ( G4UniformRand() < 0.5 ) delta0 *= -1;
      G4int A1 = A + delta0;
      if ( mIsotope.find ( Z )->second->find ( A1 ) != mIsotope.find ( Z )->second->end() )
      {
         aPhysVec = mIsotope.find ( Z )->second->find ( A1 )->second;
      }
      else
      {
         A1 = A - delta0;
         aPhysVec = mIsotope.find ( Z )->second->find ( A1 )->second;
      }

      aXSection = aPhysVec->GetValue( ek , outOfRange );
      // X^(2/3) factor
      aXSection *= std::pow ( 1.0*A/ A1 , 2.0 / 3.0 );

   }

   return aXSection;
}