Message: Re: New G4NeutroHPPhotonDist.cc G4NeutronHPVectro.cc&hh Not Logged In (login)
 Next-in-Thread Next-in-Thread
 Next-in-Forum Next-in-Forum

None Re: New G4NeutroHPPhotonDist.cc G4NeutronHPVectro.cc&hh 

Keywords: G4NeutronHPPhotonDist G4NeutronHPVector
Forum: Hadronic Processes
Re: None New G4NeutroHPPhotonDist.cc G4NeutronHPVectro.cc&hh (Enrico Billi)
Re: None Re: New G4NeutroHPPhotonDist.cc G4NeutronHPVectro.cc&hh (Sylvia Studeny)
Re: None Re: New G4NeutroHPPhotonDist.cc G4NeutronHPVectro.cc&hh (Sylvia Studeny)
Date: 07 Mar, 2007
From: Enrico Billi <billi@bo.infn.it>

Hi,

i did new corrections for this reason: with last development a gamma ray may have got less than 10 keV, in general a gamma has energy more than 10keV. For example i obtained these cascade from 197Au:

Mult : 4
Total Energy: 6.5625
0 : 6.3988554
1 : 0.0064916037        ***** Strange Gamma!
2 : 0.089561192
3 : 0.067591793
	***
Mult : 6
Total Energy: 6.5625
0 : 3.306503
1 : 0.68129768
2 : 1.6118733
3 : 0.25667473
4 : 0.00083944088       ***** Wrong Gamma!
5 : 0.70531183
	***

I did some cuts in maximum and minium (tipicaly 35KeV, first excited state of Au) energy to get for each gamma so i obtained this result:

Mult: 7
Total Energy: 6.5625 0 : 3.1568805 1 : 0.093424295 2 : 0.48907086 3 : 1.0061807 4 : 0.036964945 ****** Little energy but correct! 5 : 0.10232627 6 : 1.6776524

####################################################################
####################################################################
G4NeutronHPPhotonDist.cc
####################################################################
####################################################################

// ******************************************************************** // * DISCLAIMER * // * * // * The following disclaimer summarizes all the specific disclaimers * // * of contributors to this software. The specific disclaimers,which * // * govern, are listed with their locations in: * // * http://cern.ch/geant4/license * // * * // * 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. * // * * // * This code implementation is the intellectual property of the * // * GEANT4 collaboration. * // * By copying, distributing or modifying the Program (or any work * // * based on the Program) you indicate your acceptance of this * // * statement, and all its terms. * // ******************************************************************** // // neutron_hp -- source file // J.P. Wellisch, Nov-1996 // A prototype of the low energy neutron transport model. // // there is a lot of unused (and undebugged) code in this file. Kept for the moment just in case. @@

#include "G4NeutronHPPhotonDist.hh"
#include "G4NeutronHPLegendreStore.hh"
#include "G4Electron.hh"
#include "G4Poisson.hh"

G4bool G4NeutronHPPhotonDist::InitMean(std::ifstream & aDataFile)
{
  G4bool result = true;
  if(aDataFile >> repFlag)
  {
    aDataFile >> targetMass;
    if(repFlag==1)
    {
    // multiplicities
      aDataFile >> nDiscrete;
      disType = new G4int[nDiscrete];
      energy = new G4double[nDiscrete];
      actualMult = new G4int[nDiscrete];
      theYield = new G4NeutronHPVector[nDiscrete];
      for (G4int i=0; i<nDiscrete; i++)
      {
	aDataFile >> disType[i]>>energy[i];
	energy[i]*=eV;
	theYield[i].Init(aDataFile, eV);
      }
    }
    else if(repFlag == 2)
    {
       aDataFile >> theInternalConversionFlag;
       aDataFile >> theBaseEnergy;
       theBaseEnergy*=eV;
       aDataFile >> theInternalConversionFlag;
       aDataFile >> nGammaEnergies;
       theLevelEnergies = new G4double[nGammaEnergies];
       theTransitionProbabilities = new G4double[nGammaEnergies];
       if(theInternalConversionFlag == 2) thePhotonTransitionFraction = new G4double[nGammaEnergies];
       for(G4int  ii=0; ii<nGammaEnergies; ii++)
       {
	 if(theInternalConversionFlag == 1)
	 {
           aDataFile >> theLevelEnergies[ii] >> theTransitionProbabilities[ii];
	   theLevelEnergies[ii]*=eV;
	 }
	 else if(theInternalConversionFlag == 2)
	 {
           aDataFile >> theLevelEnergies[ii] >> theTransitionProbabilities[ii] >> thePhotonTransitionFraction[ii];
	   theLevelEnergies[ii]*=eV;
	 }
	 else
	 {
           throw G4HadronicException(__FILE__, __LINE__, "G4NeutronHPPhotonDist: Unknown conversion flag");
	 }
      }
       // Note, that this is equivalent to using the 'Gamma' classes.
      // throw G4HadronicException(__FILE__, __LINE__, "G4NeutronHPPhotonDist: Transition probability array not sampled for the moment.");
    }
    else
    {
      G4cout << "Data representation in G4NeutronHPPhotonDist: "<<repFlag<<G4endl;
      throw G4HadronicException(__FILE__, __LINE__, "G4NeutronHPPhotonDist: This data representation is not implemented.");
    }
  }
  else
  {
    result = false;
  }
  return result;
}

void G4NeutronHPPhotonDist::InitAngular(std::ifstream & aDataFile)
{
  G4int i, ii;
  //angular distributions
  aDataFile >> isoFlag;
  if (isoFlag != 1)
  {
    aDataFile >> tabulationType >> nDiscrete2 >> nIso;
    theShells = new G4double[nDiscrete2];
    theGammas = new G4double[nDiscrete2];
    for (i=0; i< nIso; i++) // isotropic photons
    {
        aDataFile >> theGammas[i] >> theShells[i];
        theGammas[i]*=eV;
        theShells[i]*=eV;
    }
    nNeu = new G4int [nDiscrete2-nIso];
    if(tabulationType==1)theLegendre=new G4NeutronHPLegendreTable *[nDiscrete2-nIso];
    if(tabulationType==2)theAngular =new G4NeutronHPAngularP *[nDiscrete2-nIso];
    for(i=nIso; i< nDiscrete2; i++)
    {
      if(tabulationType==1) 
      {
        aDataFile >> theGammas[i] >> theShells[i] >> nNeu[i-nIso];
        theGammas[i]*=eV;
        theShells[i]*=eV;
        theLegendre[i-nIso]=new G4NeutronHPLegendreTable[nNeu[i-nIso]];
        theLegendreManager.Init(aDataFile); 
        for (ii=0; ii<nNeu[i-nIso]; ii++)
        {
          theLegendre[i-nIso][ii].Init(aDataFile);
        }
      }
      else if(tabulationType==2)
      {
        aDataFile >> theGammas[i] >> theShells[i] >> nNeu[i-nIso];
        theGammas[i]*=eV;
        theShells[i]*=eV;
        theAngular[i-nIso]=new G4NeutronHPAngularP[nNeu[i-nIso]];
        for (ii=0; ii<nNeu[i-nIso]; ii++)
        {
          theAngular[i-nIso][ii].Init(aDataFile);
        }
      }
      else
      {
        G4cout << "tabulation type: tabulationType"<<G4endl;
        throw G4HadronicException(__FILE__, __LINE__, "cannot deal with this tabulation type for angular distributions.");
      }
    }
  }
}

void G4NeutronHPPhotonDist::InitEnergies(std::ifstream & aDataFile)
{
  G4int i, energyDistributionsNeeded = 0;
  for (i=0; i<nDiscrete; i++)
  {
    if( disType[i]==1) energyDistributionsNeeded =1;
  }
  if(!energyDistributionsNeeded) return;
  aDataFile >>  nPartials;
  distribution = new G4int[nPartials];
  probs = new G4NeutronHPVector[nPartials];
  partials = new G4NeutronHPPartial * [nPartials];
  G4int nen;
  G4int dummy;
  for (i=0; i<nPartials; i++)
  {
    aDataFile >> dummy;
    probs[i].Init(aDataFile, eV);
    aDataFile >> nen;
    partials[i] = new G4NeutronHPPartial(nen);
    partials[i]->InitInterpolation(aDataFile);
    partials[i]->Init(aDataFile);
  }
}

void G4NeutronHPPhotonDist::InitPartials(std::ifstream & aDataFile)
{
  aDataFile >> nDiscrete >> targetMass;
  if(nDiscrete != 1)
  {
    theTotalXsec.Init(aDataFile, eV);
  }
  G4int i;
  theGammas = new G4double[nDiscrete];
  theShells = new G4double[nDiscrete];
  isPrimary = new G4int[nDiscrete];
  disType = new G4int[nDiscrete];
  thePartialXsec = new G4NeutronHPVector[nDiscrete];
  for(i=0; i<nDiscrete; i++)
  {
    aDataFile>>theGammas[i]>>theShells[i]>>isPrimary[i]>>disType[i];
    theGammas[i]*=eV;
    theShells[i]*=eV;
    thePartialXsec[i].Init(aDataFile, eV);
  }  
}

G4ReactionProductVector * G4NeutronHPPhotonDist::GetPhotons(G4double anEnergy)
{
  // the partial cross-section case is not in this yet. @@@@
  G4int i, ii, iii;
  G4int nSecondaries = 0;
  G4ReactionProductVector * thePhotons = new G4ReactionProductVector;
  if(repFlag==1)
  {
    G4double current=0;
    for(i=0; i<nDiscrete; i++)
    {
      current = theYield[i].GetY(anEnergy);
      actualMult[i] = G4Poisson(current); // max cut-off still missing @@@
      if(nDiscrete==1&&current<1.0001) 
      {
        actualMult[i] = static_cast<G4int>(current);
        if(current<1) 
        {
          actualMult[i] = 0;
          if(G4UniformRand()<current) actualMult[i] = 1;
        }
      }
      nSecondaries += actualMult[i];
    }
    for(i=0;i<nSecondaries;i++)
    {
      G4ReactionProduct * theOne = new G4ReactionProduct;
      theOne->SetDefinition(G4Gamma::Gamma());
      thePhotons->push_back(theOne);
    }
    G4int count=0;
    G4double totalCascadeEnergy = 0.;
    G4double lastCascadeEnergy = 0.;
    G4double eGamm = 0;
    G4int maxEnergyIndex = 0;
    G4double gapEnergy = 60.*keV;

    for(i=0; i<nDiscrete; i++)
    { 
      for(ii=0; ii< actualMult[i]; ii++)
      { 
	if(disType[i]==1) // continuum
	{
          G4double  sum=0, run=0;
          for(iii=0; iii<nPartials; iii++) sum+=probs[iii].GetY(anEnergy);
          G4double random = G4UniformRand();
          G4int theP = 0;
          for(iii=0; iii<nPartials; iii++)
          {
            run+=probs[iii].GetY(anEnergy);
            theP = iii;
            if(random<run/sum) break;
          }
          if(theP==nPartials) theP=nPartials-1; // das sortiert J aus.
          sum=0; 
          G4NeutronHPVector * temp;
          temp = partials[theP]->GetY(anEnergy); //@@@ look at, seems fishy
          // Looking for TotalCascdeEnergy or LastMaxEnergy
          if (ii == 0)
          {
            maxEnergyIndex = temp->GetVectorLength()-1;
            totalCascadeEnergy = temp->GetX(maxEnergyIndex);
            lastCascadeEnergy = totalCascadeEnergy;
	    //G4cout << "\t***" << G4endl;
	    //G4cout << "Mult: " << actualMult[i] << "\nTotal Energy: " << totalCascadeEnergy << G4endl;
          }
          lastCascadeEnergy -= eGamm;
          if (ii != actualMult[i]-1) eGamm = temp->Sample(lastCascadeEnergy-(actualMult[i]-1-ii)*gapEnergy);
	  else eGamm = lastCascadeEnergy;
          //G4cout << ii << " : " << eGamm << G4endl; 
	  thePhotons->operator[](count)->SetKineticEnergy(eGamm);
          delete temp;
	}
	else // discrete
	{
          thePhotons->operator[](count)->SetKineticEnergy(energy[i]);
	}
	count++;
	if(count > nSecondaries)  throw G4HadronicException(__FILE__, __LINE__, "G4NeutronHPPhotonDist::GetPhotons inconsistancy");
      }
    }
    // now do the angular distributions...
    if( isoFlag == 1)
    {
      for (i=0; i< nSecondaries; i++)
      {
	G4double costheta = 2.*G4UniformRand()-1;
	G4double theta = std::acos(costheta);
	G4double phi = twopi*G4UniformRand();
	G4double sinth = std::sin(theta);
	G4double en = thePhotons->operator[](i)->GetTotalEnergy();
	G4ThreeVector temp(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
	thePhotons->operator[](i)->SetMomentum( temp ) ;
  //      G4cout << "Isotropic distribution in PhotonDist"<<temp<<G4endl;
      }
    }
    else
    {
      for(i=0; i<nSecondaries; i++)
      { 
	G4double currentEnergy = thePhotons->operator[](i)->GetTotalEnergy();
	for(ii=0; ii<nDiscrete2; ii++) 
	{
          if (std::abs(currentEnergy-theGammas[ii])<0.1*keV) break;
	}
	if(ii==nDiscrete2) ii--; // fix for what seems an (file12 vs file 14) inconsistancy found in the ENDF 7N14 data. @@
	if(ii<nIso)
	{
          // isotropic distribution
          G4double theta = pi*G4UniformRand();
          G4double phi = twopi*G4UniformRand();
          G4double sinth = std::sin(theta);
          G4double en = thePhotons->operator[](i)->GetTotalEnergy();
          G4ThreeVector tempVector(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
          thePhotons->operator[](i)->SetMomentum( tempVector ) ;
	}
	else if(tabulationType==1)
	{
          // legendre polynomials
          G4int it(0);
          for (iii=0; iii<nNeu[ii-nIso]; iii++) // find the neutron energy
          {
            it = iii;
	    if(theLegendre[ii-nIso][iii].GetEnergy()>anEnergy)
              break;
          }
          G4NeutronHPLegendreStore aStore(2);
          aStore.SetCoeff(1, &(theLegendre[ii-nIso][it]));  
          aStore.SetCoeff(0, &(theLegendre[ii-nIso][it-1])); 
          G4double cosTh = aStore.SampleMax(anEnergy);
          G4double theta = std::acos(cosTh);
          G4double phi = twopi*G4UniformRand();
          G4double sinth = std::sin(theta);
          G4double en = thePhotons->operator[](i)->GetTotalEnergy();
          G4ThreeVector tempVector(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
          thePhotons->operator[](i)->SetMomentum( tempVector ) ;
	}
	else
	{
          // tabulation of probabilities.
          G4int it(0);
          for (iii=0; iii<nNeu[ii-nIso]; iii++) // find the neutron energy
          {
            it = iii;
	    if(theAngular[ii-nIso][iii].GetEnergy()>anEnergy)
              break;
          }
          G4double costh = theAngular[ii-nIso][it].GetCosTh(); // no interpolation yet @@
          G4double theta = std::acos(costh);
          G4double phi = twopi*G4UniformRand();
          G4double sinth = std::sin(theta);
          G4double en = thePhotons->operator[](i)->GetTotalEnergy();
          G4ThreeVector tmpVector(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*costh );
          thePhotons->operator[](i)->SetMomentum( tmpVector ) ;
	}
      }  
    } 
  }
  else if(repFlag == 2)
  {
    G4double * running = new G4double[nGammaEnergies];
    running[0]=theTransitionProbabilities[0];
    G4int i;
    for(i=1; i<nGammaEnergies; i++)
    {
      running[i]=running[i-1]+theTransitionProbabilities[i];
    }
    G4double random = G4UniformRand();
    G4int it=0;
    for(i=0; i<nGammaEnergies; i++)
    {
      it = i;
      if(random < running[i]/running[nGammaEnergies-1]) break;
    }
    delete [] running;
    G4double totalEnergy = theBaseEnergy - theLevelEnergies[it];
    G4ReactionProduct * theOne = new G4ReactionProduct;
    theOne->SetDefinition(G4Gamma::Gamma());
    random = G4UniformRand();
    if(theInternalConversionFlag==2 && random>thePhotonTransitionFraction[it])
    {
      theOne->SetDefinition(G4Electron::Electron());
    }
    theOne->SetTotalEnergy(totalEnergy);
    if( isoFlag == 1)
    {
      G4double costheta = 2.*G4UniformRand()-1;
      G4double theta = std::acos(costheta);
      G4double phi = twopi*G4UniformRand();
      G4double sinth = std::sin(theta);
      G4double en = theOne->GetTotalEnergy();
      G4ThreeVector temp(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
      theOne->SetMomentum( temp ) ;
    }
    else
    {
      G4double currentEnergy = theOne->GetTotalEnergy();
      for(ii=0; ii<nDiscrete2; ii++) 
      {
        if (std::abs(currentEnergy-theGammas[ii])<0.1*keV) break;
      }
      if(ii==nDiscrete2) ii--; // fix for what seems an (file12 vs file 14) inconsistancy found in the ENDF 7N14 data. @@
      if(ii<nIso)
      {
        // isotropic distribution
        G4double theta = pi*G4UniformRand();
        G4double phi = twopi*G4UniformRand();
        G4double sinth = std::sin(theta);
        G4double en = theOne->GetTotalEnergy();
        G4ThreeVector tempVector(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
        theOne->SetMomentum( tempVector ) ;
      }
      else if(tabulationType==1)
      {
        // legendre polynomials
        G4int it(0);
        for (iii=0; iii<nNeu[ii-nIso]; iii++) // find the neutron energy
        {
          it = iii;
	  if(theLegendre[ii-nIso][iii].GetEnergy()>anEnergy)
            break;
        }
        G4NeutronHPLegendreStore aStore(2);
        aStore.SetCoeff(1, &(theLegendre[ii-nIso][it]));  
        aStore.SetCoeff(0, &(theLegendre[ii-nIso][it-1])); 
        G4double cosTh = aStore.SampleMax(anEnergy);
        G4double theta = std::acos(cosTh);
        G4double phi = twopi*G4UniformRand();
        G4double sinth = std::sin(theta);
        G4double en = theOne->GetTotalEnergy();
        G4ThreeVector tempVector(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
        theOne->SetMomentum( tempVector ) ;
      }
      else
      {
        // tabulation of probabilities.
        G4int it(0);
        for (iii=0; iii<nNeu[ii-nIso]; iii++) // find the neutron energy
        {
          it = iii;
	  if(theAngular[ii-nIso][iii].GetEnergy()>anEnergy)
            break;
        }
        G4double costh = theAngular[ii-nIso][it].GetCosTh(); // no interpolation yet @@
        G4double theta = std::acos(costh);
        G4double phi = twopi*G4UniformRand();
        G4double sinth = std::sin(theta);
        G4double en = theOne->GetTotalEnergy();
        G4ThreeVector tmpVector(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*costh );
        theOne->SetMomentum( tmpVector ) ;
      }
    }
    thePhotons->push_back(theOne);
  }
  else
  {
    delete thePhotons;
    thePhotons = NULL; // no gamma data available; some work needed @@@@@@@
  }    
  return thePhotons;
}

####################################################################
####################################################################
G4NeutronHPPhotonDist.hh
####################################################################
####################################################################

//
// ********************************************************************
// * 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.          *
// ********************************************************************
//
//
// $Id: G4NeutronHPPhotonDist.hh,v 1.10 2006/06/29 20:49:17 gunter Exp $
// GEANT4 tag $Name: geant4-08-01-patch-01 $
//
 // Hadronic Process: Very Low Energy Neutron X-Sections
 // original by H.P. Wellisch, TRIUMF, 14-Feb-97

#ifndef G4NeutronHPPhotonDist_h
#define G4NeutronHPPhotonDist_h 1
#include "globals.hh"
#include <fstream>
#include "G4ios.hh"
#include "globals.hh"
#include "G4NeutronHPVector.hh"
#include "G4NeutronHPLegendreTable.hh"
#include "G4NeutronHPAngularP.hh"
#include "G4NeutronHPPartial.hh"
#include "G4NeutronHPFastLegendre.hh"
#include "G4NeutronHPInterpolator.hh"
#include "G4ReactionProductVector.hh"
#include "G4ReactionProduct.hh"
#include "G4Gamma.hh"
#include "G4InterpolationManager.hh"

class G4NeutronHPPhotonDist { public:

  G4NeutronHPPhotonDist()
  {
     disType = NULL;
     energy = NULL;
     theYield = NULL;
     thePartialXsec = NULL;
     isPrimary = NULL;
     theShells = NULL;
     theGammas = NULL;
     nNeu = NULL;
     theLegendre = NULL;
     theAngular = NULL;
     distribution = NULL;
     probs = NULL;
     partials = NULL;
     actualMult = NULL;

     theLevelEnergies = NULL;
     theTransitionProbabilities = NULL;
     thePhotonTransitionFraction = NULL;
  }

  ~G4NeutronHPPhotonDist()
  {
     if(disType != NULL) delete [] disType;
     if(energy != NULL) delete [] energy;
     if(theYield != NULL) delete [] theYield;
     if(thePartialXsec != NULL) delete [] thePartialXsec;
     if(isPrimary != NULL) delete [] isPrimary;
     if(theShells != NULL) delete [] theShells;
     if(theGammas != NULL) delete [] theGammas;
     if(nNeu != NULL) delete [] nNeu;
     if(theLegendre != NULL) delete [] theLegendre;
     if(theAngular != NULL) delete [] theAngular;
     if(distribution != NULL) delete [] distribution;
     if(probs != NULL) delete [] probs;
     if(partials != NULL) delete [] partials;
     if(actualMult != NULL) delete [] actualMult;

     if(theLevelEnergies != NULL) delete theLevelEnergies;
     if(theTransitionProbabilities != NULL) delete theTransitionProbabilities;
     if(thePhotonTransitionFraction != NULL) delete thePhotonTransitionFraction;
  }

  G4bool InitMean(std::ifstream & aDataFile);

  void InitAngular(std::ifstream & aDataFile);

  void InitEnergies(std::ifstream & aDataFile);

  void InitPartials(std::ifstream & aDataFile);

  G4ReactionProductVector * GetPhotons(G4double anEnergy);

  inline G4double GetTargetMass() {return targetMass;}

  inline G4bool NeedsCascade() {return repFlag==2;}

  inline G4double GetLevelEnergy() {return theBaseEnergy;}

private:

   G4int repFlag;  //representation as multiplicities or transition probability arrays.
   G4double targetMass;

   G4int nDiscrete;  //number of discrete photons 
   G4int * disType;  // discrete, or continuum photons
   G4double * energy;  // photon energies
   G4NeutronHPVector * theYield; // multiplicity as a function of neutron energy.
   G4NeutronHPVector theTotalXsec;
   G4NeutronHPVector * thePartialXsec;
   G4int * isPrimary;

   G4int isoFlag; // isotropic or not?
   G4int tabulationType;
   G4int nDiscrete2;
   G4int nIso;
   G4double * theShells;
   G4double * theGammas;
   G4int * nNeu;
   G4InterpolationManager theLegendreManager;
   G4NeutronHPLegendreTable ** theLegendre;
   G4NeutronHPAngularP ** theAngular;

   G4int * distribution; // not used for the moment.                                 
   G4int nPartials;
   G4NeutronHPVector *  probs; // probabilities for the partial distributions.
   G4NeutronHPPartial ** partials; // the partials, parallel to the above

   G4int * actualMult;

    // for transition prob arrays start
   G4int theInternalConversionFlag;
   G4int nGammaEnergies;
   G4double theBaseEnergy;
   G4double * theLevelEnergies;
   G4double * theTransitionProbabilities;
   G4double * thePhotonTransitionFraction;
    // for transition prob arrays end

   G4NeutronHPFastLegendre theLegend; // fast look-up for leg-integrals
   G4NeutronHPInterpolator theInt; // interpolation
};

#endif

####################################################################
####################################################################
G4NeutronHPVector.cc
####################################################################
####################################################################

// // ******************************************************************** // * DISCLAIMER * // * * // * The following disclaimer summarizes all the specific disclaimers * // * of contributors to this software. The specific disclaimers,which * // * govern, are listed with their locations in: * // * http://cern.ch/geant4/license * // * * // * 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. * // * * // * This code implementation is the intellectual property of the * // * GEANT4 collaboration. * // * By copying, distributing or modifying the Program (or any work * // * based on the Program) you indicate your acceptance of this * // * statement, and all its terms. * // ******************************************************************** // // neutron_hp -- source file // J.P. Wellisch, Nov-1996 // A prototype of the low energy neutron transport model. #include "G4NeutronHPVector.hh"

  // if the ranges do not match, constant extrapolation is used.
  G4NeutronHPVector & operator + (G4NeutronHPVector & left, G4NeutronHPVector & right)
  {
    G4NeutronHPVector * result = new G4NeutronHPVector;
    G4int j=0;
    G4double x;
    G4double y;
    G4int running = 0;
    for(G4int i=0; i<left.GetVectorLength(); i++)
    {
      while(j<right.GetVectorLength())
      {
        if(right.GetX(j)<left.GetX(i)*1.001)
        {
          x = right.GetX(j);
          y = right.GetY(j)+left.GetY(x);
          result->SetData(running++, x, y);
          j++;
        }
        else if(std::abs((right.GetX(j)-left.GetX(i))/(left.GetX(i)+right.GetX(j)))>0.001)
        {
          x = left.GetX(i);
          y = left.GetY(i)+right.GetY(x);
          result->SetData(running++, x, y);
          break;
        }
        else
        {
          break;
        }
      }
      if(j==right.GetVectorLength())
      {
        x = left.GetX(i);
        y = left.GetY(i)+right.GetY(x);
        result->SetData(running++, x, y);     
      }
    }
    result->ThinOut(0.02);
    return *result;
  }

  G4NeutronHPVector::G4NeutronHPVector()
  {
    theData = new G4NeutronHPDataPoint[20]; 
    nPoints=20;
    nEntries=0;
    Verbose=0;
    theIntegral=NULL;
    totalIntegral=-1;
    isFreed = 0;
    maxValue = -DBL_MAX;
    the15percentBorderCash = -DBL_MAX;
    the50percentBorderCash = -DBL_MAX;
    label = -DBL_MAX;

  }

  G4NeutronHPVector::G4NeutronHPVector(G4int n)
  {
    nPoints=std::max(n, 20);
    theData = new G4NeutronHPDataPoint[nPoints]; 
    nEntries=0;
    Verbose=0;
    theIntegral=NULL;
    totalIntegral=-1;
    isFreed = 0;
    maxValue = -DBL_MAX;
    the15percentBorderCash = -DBL_MAX;
    the50percentBorderCash = -DBL_MAX;
  }

  G4NeutronHPVector::~G4NeutronHPVector()
  {
//    if(Verbose==1)G4cout <<"G4NeutronHPVector::~G4NeutronHPVector"<<G4endl;
    if(theData!=NULL)     
    {
      delete [] theData;
    }
//    if(Verbose==1)G4cout <<"Vector: delete theData"<<G4endl;
    if(theIntegral!=NULL) delete [] theIntegral;
//    if(Verbose==1)G4cout <<"Vector: delete theIntegral"<<G4endl;
    isFreed = 1;
  }

  G4NeutronHPVector & G4NeutronHPVector::  
  operator = (const G4NeutronHPVector & right)
  {
    if(&right == this) return *this;

    G4int i;

    totalIntegral = right.totalIntegral;
    if(right.theIntegral!=NULL) theIntegral = new G4double[right.nEntries];
    for(i=0; i<right.nEntries; i++)
    {
      SetPoint(i, right.GetPoint(i)); // copy theData
      if(right.theIntegral!=NULL) theIntegral[i] = right.theIntegral[i];
    }
    theManager = right.theManager; 
    label = right.label;

    Verbose = right.Verbose;
    the15percentBorderCash = right.the15percentBorderCash;
    the50percentBorderCash = right.the50percentBorderCash;
    theHash = right.theHash;
   return *this;
  }

  G4double G4NeutronHPVector::GetXsec(G4double e) 
  {
    if(nEntries == 0) return 0;
    if(!theHash.Prepared()) Hash();
    G4int min = theHash.GetMinIndex(e);
    G4int i;
    for(i=min ; i<nEntries; i++)
    {
      if(theData[i].GetX()>e) break;
    }
    G4int low = i-1;
    G4int high = i;
    if(i==0)
    {
      low = 0;
      high = 1;
    }
    else if(i==nEntries)
    {
      low = nEntries-2;
      high = nEntries-1;
    }
    G4double y;
    if(e<theData[nEntries-1].GetX()) 
    {
      // Protect against doubled-up x values
      if( (theData[high].GetX()-theData[low].GetX())/theData[high].GetX() < 0.000001)
      {
        y = theData[low].GetY();
      }
      else
      {
        y = theInt.Interpolate(theManager.GetScheme(high), e, 
                               theData[low].GetX(), theData[high].GetX(),
		  	       theData[low].GetY(), theData[high].GetY());
      }
    }
    else
    {
      y=theData[nEntries-1].GetY();
    }
    return y;
  }

  void G4NeutronHPVector::Dump()
  {
    G4cout << nEntries<<G4endl;
    for(G4int i=0; i<nEntries; i++)
    {
      G4cout << theData[i].GetX()<<" ";
      G4cout << theData[i].GetY()<<" ";
//      if (i!=1&&i==5*(i/5)) G4cout << G4endl;
      G4cout << G4endl;
    }
    G4cout << G4endl;
  }

  void G4NeutronHPVector::Check(G4int i)
  {
    if(i>nEntries) throw G4HadronicException(__FILE__, __LINE__, "Skipped some index numbers in G4NeutronHPVector");
    if(i==nPoints)
    {
      nPoints = static_cast<G4int>(1.2*nPoints);
      G4NeutronHPDataPoint * buff = new G4NeutronHPDataPoint[nPoints];
      for (G4int j=0; j<nEntries; j++) buff[j] = theData[j];
      delete [] theData;
      theData = buff;
    }
    if(i==nEntries) nEntries=i+1;
  }

  void G4NeutronHPVector::
  Merge(G4InterpolationScheme aScheme, G4double aValue, 
        G4NeutronHPVector * active, G4NeutronHPVector * passive)
  { 
    // interpolate between labels according to aScheme, cut at aValue, 
    // continue in unknown areas by substraction of the last difference.

    CleanUp();
    G4int s = 0, n=0, m=0;
    G4NeutronHPVector * tmp;
    G4int a = s, p = n, t;
    while ( a<active->GetVectorLength() )
    {
      if(active->GetEnergy(a) <= passive->GetEnergy(p))
      {
        G4double xa  = active->GetEnergy(a);
        G4double yy = theInt.Interpolate(aScheme, aValue, active->GetLabel(), passive->GetLabel(),
                                                          active->GetXsec(a), passive->GetXsec(xa));
        SetData(m, xa, yy);
        theManager.AppendScheme(m, active->GetScheme(a));
        m++;
        a++;
        G4double xp = passive->GetEnergy(p);
        if( std::abs(std::abs(xp-xa)/xa)<0.0000001&&a<active->GetVectorLength() ) 
        {
          p++;
          tmp = active; t=a;
          active = passive; a=p;
          passive = tmp; p=t;
        }
      } else {
        tmp = active; t=a;
        active = passive; a=p;
        passive = tmp; p=t;
      }
    }

    G4double deltaX = passive->GetXsec(GetEnergy(m-1)) - GetXsec(m-1);
    while (p!=passive->GetVectorLength()&&passive->GetEnergy(p)<=aValue)
    {
      G4double anX;
      anX = passive->GetXsec(p)-deltaX;
      if(anX>0)
      {
        if(std::abs(GetEnergy(m-1)-passive->GetEnergy(p))/passive->GetEnergy(p)>0.0000001)
        {
          SetData(m, passive->GetEnergy(p), anX);
          theManager.AppendScheme(m++, passive->GetScheme(p));
        }
      }
      p++;
    }
    // Rebuild the Hash;
    if(theHash.Prepared()) 
    {
      ReHash();
    }
  }

  void G4NeutronHPVector::ThinOut(G4double precision)
  {
    // anything in there?
    if(GetVectorLength()==0) return;
    // make the new vector
    G4NeutronHPDataPoint * aBuff = new G4NeutronHPDataPoint[nPoints];
    G4double x, x1, x2, y, y1, y2;
    G4int count = 0, current = 2, start = 1;

    // First element always goes and is never tested.
    aBuff[0] = theData[0];

    // Find the rest
    while(current < GetVectorLength())
    {
      x1=aBuff[count].GetX();
      y1=aBuff[count].GetY();
      x2=theData[current].GetX();
      y2=theData[current].GetY();
      for(G4int j=start; j<current; j++)
      {
	x = theData[j].GetX();
	if(x1-x2 == 0) y = (y2+y1)/2.;
	else y = theInt.Lin(x, x1, x2, y1, y2);
	if (std::abs(y-theData[j].GetY())>precision*y)
	{
	  aBuff[++count] = theData[current-1]; // for this one, everything was fine
          start = current; // the next candidate
	  break;
	}
      }
      current++ ;
    }
    // The last one also always goes, and is never tested.
    aBuff[++count] = theData[GetVectorLength()-1];
    delete [] theData;
    theData = aBuff;
    nEntries = count+1;

    // Rebuild the Hash;
    if(theHash.Prepared()) 
    {
      ReHash();
    }
  }

  G4bool G4NeutronHPVector::IsBlocked(G4double aX)
  {
    G4bool result = false;
    std::vector<G4double>::iterator i;
    for(i=theBlocked.begin(); i!=theBlocked.end(); i++)
    {
      G4double aBlock = *i;
      if(std::abs(aX-aBlock) < 0.1*MeV)
      {
        result = true;
	theBlocked.erase(i);
	break;
      }
    }
    return result;
  }
  G4double G4NeutronHPVector::Sample() // Samples X according to distribution Y
  {
    G4double result;
    G4int j;
    for(j=0; j<GetVectorLength(); j++)
    {
      if(GetY(j)<0) SetY(j, 0);
    }

    if(theBuffered.size() !=0 && G4UniformRand()<0.5)
    {
      result = theBuffered[0];
      theBuffered.erase(theBuffered.begin());
      if(result < GetX(GetVectorLength()-1) ) return result;
    }
    if(GetVectorLength()==1)
    {
      result = theData[0].GetX();
    }
    else
    {
      if(theIntegral==NULL) IntegrateAndNormalise();
      do
      {
        G4double value, test, baseline;
        baseline = theData[GetVectorLength()-1].GetX()-theData[0].GetX();
        G4double rand;
        do
        {
          value = baseline*G4UniformRand();
          value += theData[0].GetX();
          test = GetY(value)/maxValue;
          rand = G4UniformRand();
        }
        while(test<rand);
        result = value;
      }
      while(IsBlocked(result));
    }
    return result;
  }

  G4double G4NeutronHPVector::Sample(G4double anEnergy) // Samples X according to distribution Y
  {
    G4double result;
    G4int j, k;

    for(k=0; k<GetVectorLength(); k++)
    {
	if (GetX(k)>anEnergy)
	{
	    theIntegral = NULL;
	    break;
	}
    }

    for(j=0; j<GetVectorLength(); j++)
    {
      if(GetY(j)<0 || j>=k) SetY(j, 0);
    }

    if(theBuffered.size() !=0 && G4UniformRand()<0.5)
    {
      result = theBuffered[0];
      theBuffered.erase(theBuffered.begin());
      if(result < GetX(GetVectorLength()-1) ) return result;
    }
    if(GetVectorLength()==1)
    {
      result = theData[0].GetX();
    }
    else
    {
      if(theIntegral==NULL) IntegrateAndNormalise();
      do
      {
        G4double value, test, baseline;
        baseline = anEnergy-35*keV;
        G4double rand;
        do
        {
          value = baseline*G4UniformRand();
          value += 35*keV;
          test = GetY(value)/maxValue;
          rand = G4UniformRand();
        }
        while(test<rand);
        result = value;
      }
      while(IsBlocked(result));
    }
    return result;
  }

  G4double G4NeutronHPVector::Get15percentBorder()
  {
    if(the15percentBorderCash>-DBL_MAX/2.) return the15percentBorderCash;
    G4double result;
    if(GetVectorLength()==1)
    {
      result = theData[0].GetX();
      the15percentBorderCash = result;
    }
    else
    {
      if(theIntegral==NULL) IntegrateAndNormalise();
      G4int i;
      result = theData[GetVectorLength()-1].GetX();
      for(i=0;i<GetVectorLength();i++)
      {
	if(theIntegral[i]/theIntegral[GetVectorLength()-1]>0.15) 
	{
	  result = theData[std::min(i+1, GetVectorLength()-1)].GetX();
          the15percentBorderCash = result;
	  break;
	}
      }
      the15percentBorderCash = result;
    }
    return result;
  }

  G4double G4NeutronHPVector::Get50percentBorder()
  {
    if(the50percentBorderCash>-DBL_MAX/2.) return the50percentBorderCash;
    G4double result;
    if(GetVectorLength()==1)
    {
      result = theData[0].GetX();
      the50percentBorderCash = result;
    }
    else
    {
      if(theIntegral==NULL) IntegrateAndNormalise();
      G4int i;
      G4double x = 0.5;
      result = theData[GetVectorLength()-1].GetX();
      for(i=0;i<GetVectorLength();i++)
      {
	if(theIntegral[i]/theIntegral[GetVectorLength()-1]>x) 
	{
	  G4int it;
	  it = i;
	  if(it == GetVectorLength()-1)
	  {
	    result = theData[GetVectorLength()-1].GetX();
	  }
	  else
	  {
	    G4double x1, x2, y1, y2;
	    x1 = theIntegral[i-1]/theIntegral[GetVectorLength()-1];
	    x2 = theIntegral[i]/theIntegral[GetVectorLength()-1];
	    y1 = theData[i-1].GetX();
	    y2 = theData[i].GetX();
	    result = theLin.Lin(x, x1, x2, y1, y2);
	  }
          the50percentBorderCash = result;
	  break;
	}
      }
      the50percentBorderCash = result;
    }
    return result;
  }

####################################################################
####################################################################
G4NeutronHPVector.hh
####################################################################
####################################################################

// // ******************************************************************** // * DISCLAIMER * // * * // * The following disclaimer summarizes all the specific disclaimers * // * of contributors to this software. The specific disclaimers,which * // * govern, are listed with their locations in: * // * http://cern.ch/geant4/license * // * * // * 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. * // * * // * This code implementation is the intellectual property of the * // * GEANT4 collaboration. * // * By copying, distributing or modifying the Program (or any work * // * based on the Program) you indicate your acceptance of this * // * statement, and all its terms. * // ******************************************************************** //

#ifndef G4NeutronHPVector_h
#define G4NeutronHPVector_h 1

#include "G4NeutronHPDataPoint.hh"
#include "G4PhysicsVector.hh"
#include "G4NeutronHPInterpolator.hh"
#include "Randomize.hh"
#include "G4ios.hh"
#include <fstream>
#include "G4InterpolationManager.hh"
#include "G4NeutronHPInterpolator.hh"
#include "G4NeutronHPHash.hh"
#include <cmath>
#include <vector>

class G4NeutronHPVector
{
  friend G4NeutronHPVector & operator + (G4NeutronHPVector & left, 
                                         G4NeutronHPVector & right);

  public:

  G4NeutronHPVector();

  G4NeutronHPVector(G4int n);

  ~G4NeutronHPVector();

  G4NeutronHPVector & operator = (const G4NeutronHPVector & right);

  inline void SetVerbose(G4int ff)
  {
    Verbose = ff;
  }

  inline void Times(G4double factor)
  {
    G4int i;
    for(i=0; i<nEntries; i++)
    {
      theData[i].SetY(theData[i].GetY()*factor);
    }
    if(theIntegral!=NULL)
    {
      theIntegral[i] *= factor;
    }
  }

  inline void SetPoint(G4int i, const G4NeutronHPDataPoint & it)
  {
    G4double x = it.GetX();
    G4double y = it.GetY();
    SetData(i, x, y);
  }

  inline void SetData(G4int i, G4double x, G4double y) 
  { 
//    G4cout <<"G4NeutronHPVector::SetData called"<<nPoints<<" "<<nEntries<<G4endl;
    Check(i);
    if(y>maxValue) maxValue=y;
    theData[i].SetData(x, y);
  }
  inline void SetX(G4int i, G4double e)
  {
    Check(i);
    theData[i].SetX(e);
  }
  inline void SetEnergy(G4int i, G4double e)
  {
    Check(i);
    theData[i].SetX(e);
  }
  inline void SetY(G4int i, G4double x)
  {
    Check(i);
    if(x>maxValue) maxValue=x;
    theData[i].SetY(x);
  }
  inline void SetXsec(G4int i, G4double x)
  {
    Check(i);
    if(x>maxValue) maxValue=x;
    theData[i].SetY(x);
  }
  inline G4double GetEnergy(G4int i) const { return theData[i].GetX(); }
  inline G4double GetXsec(G4int i) { return theData[i].GetY(); }
  inline G4double GetX(G4int i) const 
  { 
    if (i<0) i=0;
    if(i>=GetVectorLength()) i=GetVectorLength()-1;
    return theData[i].GetX();
  }
  inline const G4NeutronHPDataPoint & GetPoint(G4int i) const { return theData[i]; }

  void Hash() 
  {
    G4int i;
    G4double x, y;
    for(i=0 ; i<nEntries; i++)
    {
      if(0 == (i+1)%10)
      {
        x = GetX(i);
	y = GetY(i);
	theHash.SetData(i, x, y);
      }
    }
  }

  void ReHash()
  {
    theHash.Clear();
    Hash();
  }

  G4double GetXsec(G4double e);
  G4double GetXsec(G4double e, G4int min)
  {
    G4int i;
    for(i=min ; i<nEntries; i++)
    {
      if(theData[i].GetX()>e) break;
    }
    G4int low = i-1;
    G4int high = i;
    if(i==0)
    {
      low = 0;
      high = 1;
    }
    else if(i==nEntries)
    {
      low = nEntries-2;
      high = nEntries-1;
    }
    G4double y;
    if(e<theData[nEntries-1].GetX()) 
    {
      // Protect against doubled-up x values
      if( (theData[high].GetX()-theData[low].GetX())/theData[high].GetX() < 0.000001)
      {
        y = theData[low].GetY();
      }
      else
      {
        y = theInt.Interpolate(theManager.GetScheme(high), e, 
                               theData[low].GetX(), theData[high].GetX(),
		  	       theData[low].GetY(), theData[high].GetY());
      }
    }
    else
    {
      y=theData[nEntries-1].GetY();
    }
    return y;
  }

  inline G4double GetY(G4double x)  {return GetXsec(x);}
  inline G4int GetVectorLength() const {return nEntries;}

  inline G4double GetY(G4int i)
  { 
    if (i<0) i=0;
    if(i>=GetVectorLength()) i=GetVectorLength()-1;
    return theData[i].GetY(); 
  }

  inline G4double GetY(G4int i) const
  {
    if (i<0) i=0;
    if(i>=GetVectorLength()) i=GetVectorLength()-1;
    return theData[i].GetY(); 
  }
  void Dump();

  inline void InitInterpolation(std::ifstream & aDataFile)
  {
    theManager.Init(aDataFile);
  }

  void Init(std::ifstream & aDataFile, G4int total, G4double ux=1., G4double uy=1.)
  {
    G4double x,y;
    for (G4int i=0;i<total;i++)
    {
      aDataFile >> x >> y;
      x*=ux;
      y*=uy;
      SetData(i,x,y);
      if(0 == nEntries%10)
      {
        theHash.SetData(nEntries-1, x, y);
      }
    }
  }

  void Init(std::ifstream & aDataFile,G4double ux=1., G4double uy=1.)
  {
    G4int total;
    aDataFile >> total;
    if(theData!=NULL) delete [] theData;
    theData = new G4NeutronHPDataPoint[total]; 
    nPoints=total;
    nEntries=0;    
    theManager.Init(aDataFile);
    Init(aDataFile, total, ux, uy);
  }

  void ThinOut(G4double precision);

  inline void SetLabel(G4double aLabel)
  {
    label = aLabel;
  }

  inline G4double GetLabel()
  {
    return label;
  }

  inline void CleanUp()
  {
    nEntries=0;   
    theManager.CleanUp();
    maxValue = -DBL_MAX;
    theHash.Clear();
  }

  // merges the vectors active and passive into *this
  inline void Merge(G4NeutronHPVector * active, G4NeutronHPVector * passive)
  {
    CleanUp();
    G4int s = 0, n=0, m=0;
    G4NeutronHPVector * tmp;
    G4int a = s, p = n, t;
    while (a<active->GetVectorLength()&&p<passive->GetVectorLength())
    {
      if(active->GetEnergy(a) <= passive->GetEnergy(p))
      {
        G4double xa  = active->GetEnergy(a);
        G4double yy = active->GetXsec(a);
        SetData(m, xa, yy);
        theManager.AppendScheme(m, active->GetScheme(a));
        m++;
        a++;
        G4double xp = passive->GetEnergy(p);
        if( std::abs(std::abs(xp-xa)/xa)<0.001 ) p++;
      } else {
        tmp = active; 
        t=a;
        active = passive; 
        a=p;
        passive = tmp; 
        p=t;
      }
    }
    while (a!=active->GetVectorLength())
    {
      SetData(m, active->GetEnergy(a), active->GetXsec(a));
      theManager.AppendScheme(m++, active->GetScheme(a));
      a++;
    }
    while (p!=passive->GetVectorLength())
    {
      if(std::abs(GetEnergy(m)-passive->GetEnergy(p))/passive->GetEnergy(p)>0.001)
      {
        SetData(m, passive->GetEnergy(p), passive->GetXsec(p));
        theManager.AppendScheme(m++, active->GetScheme(p));
      }
      p++;
    }
  }    

  void Merge(G4InterpolationScheme aScheme, G4double aValue, 
             G4NeutronHPVector * active, G4NeutronHPVector * passive);

  G4double SampleLin() // Samples X according to distribution Y, linear int
  {
    G4double result;
    if(theIntegral==NULL) IntegrateAndNormalise();
    if(GetVectorLength()==1)
    {
      result = theData[0].GetX();
    }
    else
    {
      G4int i;
      G4double rand = G4UniformRand();

// this was replaced // for(i=1;i<GetVectorLength();i++) // { // if(rand<theIntegral[i]/theIntegral[GetVectorLength()-1]) break; // }

// by this (begin)
      for(i=GetVectorLength()-1; i>=0 ;i--)
      {
	if(rand>theIntegral[i]/theIntegral[GetVectorLength()-1]) break;
      }
      if(i!=GetVectorLength()-1) i++;
// until this (end)

      G4double x1, x2, y1, y2;
      y1 = theData[i-1].GetX();
      x1 = theIntegral[i-1];
      y2 = theData[i].GetX();
      x2 = theIntegral[i];
      if(std::abs((y2-y1)/y2)<0.0000001) // not really necessary, since the case is excluded by construction
      {
	y1 = theData[i-2].GetX();
	x1 = theIntegral[i-2];
      }
      result = theLin.Lin(rand, x1, x2, y1, y2);
    }
    return result;
  }

  G4double Sample(); // Samples X according to distribution Y
  G4double Sample(G4double anEnergy);

  G4double * Debug()
  {
    return theIntegral;
  }

  inline void IntegrateAndNormalise()
  {
    G4int i;
    if(theIntegral!=NULL) return;
    theIntegral = new G4double[nEntries];
    if(nEntries == 1)
    {
      theIntegral[0] = 1;
      return;
    }
    theIntegral[0] = 0;
    G4double sum = 0;
    for(i=1;i<GetVectorLength();i++)
    {
      if(std::abs((theData[i].GetX()-theData[i-1].GetX())/theData[i].GetX())>0.0000001)
      {
        sum+= 0.5*(theData[i].GetY()+theData[i-1].GetY())*
                  (theData[i].GetX()-theData[i-1].GetX());
      }
      theIntegral[i] = sum;
    }
    G4double total = theIntegral[GetVectorLength()-1];
    for(i=1;i<GetVectorLength();i++)
    {
      theIntegral[i]/=total;
    }
  }

  inline void Integrate() 
  {
    G4int i;
    if(nEntries == 1)
    {
      totalIntegral = 0;
      return;
    }
    G4double sum = 0;
    for(i=1;i<GetVectorLength();i++)
    {
      if(std::abs((theData[i].GetX()-theData[i-1].GetX())/theData[i].GetX())>0.0000001)
      {
        G4double x1 = theData[i-1].GetX();
        G4double x2 = theData[i].GetX();
        G4double y1 = theData[i-1].GetY();
        G4double y2 = theData[i].GetY();
        G4InterpolationScheme aScheme = theManager.GetScheme(i);
        if(aScheme==LINLIN||aScheme==CLINLIN||aScheme==ULINLIN)
        {
          sum+= 0.5*(y2+y1)*(x2-x1);
        }
        else if(aScheme==LINLOG||aScheme==CLINLOG||aScheme==ULINLOG)
        {
          G4double a = y1;
          G4double b = (y2-y1)/(std::log(x2)-std::log(x1));
          sum+= (a-b)*(x2-x1) + b*(x2*std::log(x2)-x1*std::log(x1));
        }
        else if(aScheme==LOGLIN||aScheme==CLOGLIN||aScheme==ULOGLIN)
        {
          G4double a = std::log(y1);
          G4double b = (std::log(y2)-std::log(y1))/(x2-x1);
          sum += (std::exp(a)/b)*(std::exp(b*x2)-std::exp(b*x1));
        }
        else if(aScheme==HISTO||aScheme==CHISTO||aScheme==UHISTO)
        {
          sum+= y1*(x2-x1);
        }
        else if(aScheme==LOGLOG||aScheme==CLOGLOG||aScheme==ULOGLOG)
        {
          G4double a = std::log(y1);
          G4double b = (std::log(y2)-std::log(y1))/(std::log(x2)-std::log(x1));
          sum += (std::exp(a)/(b+1))*(std::pow(x2,b+1)-std::pow(x1,b+1));
        }
        else
        {
          throw G4HadronicException(__FILE__, __LINE__, "Unknown interpolation scheme in G4NeutronHPVector::Integrate");
        }

      }
    }
    totalIntegral = sum;
  }

  inline G4double GetIntegral() // linear interpolation; use with care
  { 
    if(totalIntegral<-0.5) Integrate();
    return totalIntegral; 
  }

  inline void SetInterpolationManager(const G4InterpolationManager & aManager)
  {
    theManager = aManager;
  }

  inline const G4InterpolationManager & GetInterpolationManager() const
  {
    return theManager;
  }

  inline void SetInterpolationManager(G4InterpolationManager & aMan)
  {
    theManager = aMan;
  }

  inline void SetScheme(G4int aPoint, const G4InterpolationScheme & aScheme)
  {
    theManager.AppendScheme(aPoint, aScheme);
  }

  inline G4InterpolationScheme GetScheme(G4int anIndex)
  {
    return theManager.GetScheme(anIndex);
  }

  G4double GetMeanX()
  {
    G4double result;
    G4double running = 0;
    G4double weighted = 0;
    for(G4int i=1; i<nEntries; i++)
    {
      running += theInt.GetBinIntegral(theManager.GetScheme(i-1),
                           theData[i-1].GetX(), theData[i].GetX(),
                           theData[i-1].GetY(), theData[i].GetY());
      weighted += theInt.GetWeightedBinIntegral(theManager.GetScheme(i-1),
                           theData[i-1].GetX(), theData[i].GetX(),
                           theData[i-1].GetY(), theData[i].GetY());
    }  
    result = weighted / running;  
    return result;
  }

  void Block(G4double aX)
  {
    theBlocked.push_back(aX);
  }

  void Buffer(G4double aX)
  {
    theBuffered.push_back(aX);
  }

  std::vector<G4double> GetBlocked() {return theBlocked;}
  std::vector<G4double> GetBuffered() {return theBuffered;}

  void SetBlocked(const std::vector<G4double> &aBlocked) {theBlocked = aBlocked;}
  void SetBuffered(const std::vector<G4double> &aBuffer) {theBuffered = aBuffer;}

  G4double Get15percentBorder();
  G4double Get50percentBorder();

  private:

  void Check(G4int i);

  G4bool IsBlocked(G4double aX);

  private:

  G4NeutronHPInterpolator theLin;

  private:

  G4double totalIntegral;

  G4NeutronHPDataPoint * theData; // the data
  G4InterpolationManager theManager; // knows how to interpolate the data.
  G4double * theIntegral;
  G4int nEntries;
  G4int nPoints;
  G4double label;

  G4NeutronHPInterpolator theInt;
  G4int Verbose;
  // debug only
  G4int isFreed;

  G4NeutronHPHash theHash;
  G4double maxValue;

  std::vector<G4double> theBlocked;
  std::vector<G4double> theBuffered;
  G4double the15percentBorderCash;
  G4double the50percentBorderCash;

};

#endif

Inline Depth:
 1 1
 All All
Outline Depth:
 1 1
 2 2
 All All
Add message: (add)

 Add Message Add Message
to: "Re: New G4NeutroHPPhotonDist.cc G4NeutronHPVectro.cc&hh"

 Subscribe Subscribe

This site runs SLAC HyperNews version 1.11-slac-98, derived from the original HyperNews


[ Geant 4 Home | Geant 4 HyperNews | Search | Request New Forum | Feedback ]