Skip to content
Snippets Groups Projects
Commit c126c0ed authored by BaM's avatar BaM
Browse files 

git-svn-id: svn+ssh://svn.in2p3.fr/class@328 0e7d625b-0364-4367-a6be-d5be4a48d228
parent 58e2418f
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 8 additions and 2178 deletions
source/branches/CLASSV3/src/EQM_LIN_PWR_MOX.cxx deleted 100644 → 0
+ 0
− 262
View file @ 58e2418f
#include "EQM_LIN_PWR_MOX.hxx"
#include "CLASSConstante.hxx"
#include <vector>
#include "StringLine.hxx"
#include "CLASSLogger.hxx"
#include "IsotopicVector.hxx"
EQM_LIN_PWR_MOX::EQM_LIN_PWR_MOX(string WeightPath):EquivalenceModel(new CLASSLogger("EQM_LIN_PWR_MOX.log"))
{
fWeightPath = WeightPath;
ifstream DataDB(fWeightPath.c_str()); // Open the File
if(!DataDB)
WARNING << "Can't open \"" << fWeightPath << "\"\n" << endl;
string line;
int start = 0; // First Get Fuel Parameter
getline(DataDB, line);
if( StringLine::NextWord(line, start, ' ') != "PARAM")
{
ERROR << " Bad Database file : " << fWeightPath << " Can't find the Parameter of the DataBase " << endl;
exit (1);
}
while(start < (int)line.size())
fFuelParameter.push_back(atof(StringLine::NextWord(line, start, ' ').c_str()));
INFO << fFuelParameter.size() << " have been read " << endl;
//-----------------------------------------------------------------------------//
//-----------------------------------------------------------------------------//
//-----------------------------------------------------------------------------//
//-----------------------------------------------------------------------------//
//-----------------------------------------------------------------------------//
// ADD ENrichment of the U reading !!!!!!!!!!!!!!!!!!!!!!!!!!!! //
//-----------------------------------------------------------------------------//
//-----------------------------------------------------------------------------//
//-----------------------------------------------------------------------------//
//-----------------------------------------------------------------------------//
ZAI U8(92,238,0);
ZAI U5(92,235,0);
double U5_enrich= 0.0025;
fFertileList = U5*U5_enrich + U8*(1-U5_enrich);
ZAI Pu8(94,238,0);
ZAI Pu9(94,239,0);
ZAI Pu0(94,240,0);
ZAI Pu1(94,241,0);
ZAI Pu2(94,242,0);
fFissileList = Pu8*1+Pu9*1+Pu0*1+Pu1*1+Pu2*1;
}
EQM_LIN_PWR_MOX::EQM_LIN_PWR_MOX(CLASSLogger* log, string WeightPath):EquivalenceModel(log)
{
fWeightPath = WeightPath;
ifstream DataDB(fWeightPath.c_str()); // Open the File
if(!DataDB)
WARNING << " Can't open \"" << fWeightPath << "\"\n" << endl;
string line;
int start = 0; // First Get Fuel Parameter
getline(DataDB, line);
if( StringLine::NextWord(line, start, ' ') != "PARAM")
{
ERROR << " Bad Database file : " << fWeightPath << " Can't find the Parameter of the DataBase"<< endl;
exit (1);
}
while(start < (int)line.size())
fFuelParameter.push_back(atof(StringLine::NextWord(line, start, ' ').c_str()));
INFO << fFuelParameter.size() << " have been read"<< endl;
//-----------------------------------------------------------------------------//
//-----------------------------------------------------------------------------//
//-----------------------------------------------------------------------------//
//-----------------------------------------------------------------------------//
//-----------------------------------------------------------------------------//
// ADD ENrichment of the U reading !!!!!!!!!!!!!!!!!!!!!!!!!!!! //
//-----------------------------------------------------------------------------//
//-----------------------------------------------------------------------------//
//-----------------------------------------------------------------------------//
//-----------------------------------------------------------------------------//
ZAI U8(92,238,0);
ZAI U5(92,235,0);
double U5_enrich= 0.0025;
fFertileList = U5*U5_enrich + U8*(1-U5_enrich);
ZAI Pu8(94,238,0);
ZAI Pu9(94,239,0);
ZAI Pu0(94,240,0);
ZAI Pu1(94,241,0);
ZAI Pu2(94,242,0);
fFissileList = Pu8*1+Pu9*1+Pu0*1+Pu1*1+Pu2*1;
}
EQM_LIN_PWR_MOX::~EQM_LIN_PWR_MOX()
{
}
//________________________________________________________________________
vector<double> EQM_LIN_PWR_MOX::BuildFuel(double BurnUp, double HMMass,vector<IsotopicVector> FissilArray, vector<IsotopicVector> FertilArray)
{
//-----------------------------------------------------------------------------//
//-----------------------------------------------------------------------------//
//-----------------------------------------------------------------------------//
//-----------------------------------------------------------------------------//
//-----------------------------------------------------------------------------//
// ADD ENrichment of the U check ++ Check Un seul fertile !!!! //
//-----------------------------------------------------------------------------//
//-----------------------------------------------------------------------------//
//-----------------------------------------------------------------------------//
//-----------------------------------------------------------------------------//
vector<double> lambda;
for(int i = 0; i < (int) (FissilArray.size() + FertilArray.size()); i++)
lambda.push_back(0);
double Na = 6.02214129e23; //N Avogadro
IsotopicVector FullUsedStock;
IsotopicVector stock;
bool FuelBuild = false;
if(FissilArray.size() == 0)
{
for(int i = 0; i < (int)lambda.size(); i++)
lambda[i] = -1;
FuelBuild = true;
}
int N_FissilStock_OnCheck = 0;
while(!FuelBuild)
{
double nPu_0 = 0;
double MPu_0 = 0;
{
map<ZAI ,double>::iterator it;
map<ZAI ,double> isotopicquantity = FullUsedStock.GetSpeciesComposition(94).GetIsotopicQuantity();
for( it = isotopicquantity.begin(); it != isotopicquantity.end(); it++ )
nPu_0 += (*it).second;
isotopicquantity = (FullUsedStock.GetSpeciesComposition(94) + ZAI(94,241,0)*FullUsedStock.GetZAIIsotopicQuantity(95,241,0)).GetIsotopicQuantity(); //Add the 241Am as 241Pu... the Pu is not old in the Eq Model but is in the FissileArray....;
for( it = isotopicquantity.begin(); it != isotopicquantity.end(); it++ )
MPu_0 += (*it).second*cZAIMass.fZAIMass.find( (*it).first )->second/Na*1e-6;
}
stock = FissilArray[N_FissilStock_OnCheck];
double nPu_1 = 0;
double MPu_1 = 0;
double Sum_AlphaI_nPuI = 0;
double Sum_AlphaI_nPuI0 = 0;
{
map<ZAI ,double>::iterator it;
map<ZAI ,double> isotopicquantity = stock.GetSpeciesComposition(94).GetIsotopicQuantity();
for( it = isotopicquantity.begin(); it != isotopicquantity.end(); it++ )
{
if ((*it).first.A() >= 238 && (*it).first.A() <= 242)
{
nPu_1 += (*it).second;
Sum_AlphaI_nPuI += fFuelParameter[(*it).first.A() -237]*(*it).second;
}
}
isotopicquantity = (stock.GetSpeciesComposition(94) + ZAI(94,241,0)*stock.GetZAIIsotopicQuantity(95,241,0)).GetIsotopicQuantity(); //Add the 241Am as 241Pu... the Pu is not old in the Eq Model but is in the FissileArray....
for( it = isotopicquantity.begin(); it != isotopicquantity.end(); it++ )
if ((*it).first.A() >= 238 && (*it).first.A() <= 242)
{
MPu_1 += (*it).second * (cZAIMass.fZAIMass.find( (*it).first )->second)/Na*1e-6;
}
isotopicquantity = FullUsedStock.GetSpeciesComposition(94).GetIsotopicQuantity();
for( it = isotopicquantity.begin(); it != isotopicquantity.end(); it++ )
if ((*it).first.A() >= 238 && (*it).first.A() <= 242)
{
Sum_AlphaI_nPuI0 += fFuelParameter[(*it).first.A() -237]*(*it).second;
}
}
double StockFactionToUse = 0;
double NT = HMMass*1e6 * Na / (cZAIMass.GetMass( ZAI(92,238,0) ) * 0.997
+ cZAIMass.GetMass( ZAI(92,235,0) ) * 0.003 );
double N1 = (BurnUp - fFuelParameter[6]) * NT;
double N2 = -Sum_AlphaI_nPuI0;
double N3 = -fFuelParameter[0] * Na / (cZAIMass.fZAIMass.find( ZAI(92,238,0) )->second*0.997
+ cZAIMass.fZAIMass.find( ZAI(92,235,0) )->second*0.003 )
* (HMMass*1e6 - MPu_0*1e6);
double D1 = Sum_AlphaI_nPuI;
double D2 = -fFuelParameter[0] * MPu_1*1e6 * Na / (cZAIMass.fZAIMass.find( ZAI(92,238,0) )->second*0.997
+ cZAIMass.fZAIMass.find( ZAI(92,235,0) )->second*0.003 ) ;
StockFactionToUse = (N1 + N2 + N3) / (D1 + D2);
if(StockFactionToUse < 0)
{
WARNING << "!!!FabricationPlant!!! Oups Bug in calculating stock fraction to use "<< endl;
lambda[N_FissilStock_OnCheck] = 0.;
N_FissilStock_OnCheck++;
FuelBuild = false;
}
else if( StockFactionToUse > 1 )
{
FullUsedStock += stock;
lambda[N_FissilStock_OnCheck] = 1;
N_FissilStock_OnCheck++;
FuelBuild = false;
}
else
{
lambda[N_FissilStock_OnCheck] = StockFactionToUse;
FuelBuild = true;
double U8_Quantity = (HMMass - (MPu_0+StockFactionToUse*MPu_1 ))/(cZAIMass.fZAIMass.find( ZAI(92,238,0) )->second*0.997 + cZAIMass.fZAIMass.find( ZAI(92,235,0) )->second*0.003 )*Na/1e-6;
lambda.back() = U8_Quantity / FertilArray[0].GetSumOfAll();
}
if( N_FissilStock_OnCheck == (int) FissilArray.size() ) // Check if the last Fissil stock has been tested... quit if so...
{
for(int i = 0; i < (int)lambda.size(); i++)
lambda[i] = -1;
FuelBuild = true;
}
}
return lambda;
}
source/branches/CLASSV3/src/EQM_MLP_PWR_MOX.cxx deleted 100644 → 0
+ 0
− 182
View file @ 58e2418f
#include "EquivalenceModel.hxx"
#include "EQM_MLP_PWR_MOX.hxx"
#include "CLASSLogger.hxx"
#include "StringLine.hxx"
#include <string>
#include <iostream>
#include <fstream>
#include <algorithm>
#include <cmath>
#include <cassert>
#include "TSystem.h"
#include "TMVA/Reader.h"
#include "TMVA/Tools.h"
#include "TMVA/MethodCuts.h"
//________________________________________________________________________
//
// EQM_MLP_MOX
//
// Equivalenve Model based on multi layer perceptron from TMVA (root cern)
// For REP MOX use
//
//________________________________________________________________________
EQM_MLP_MOX::EQM_MLP_MOX(string TMVAWeightPath):EquivalenceModel(new CLASSLogger("EQM_MLP_MOX.log"))
{
fTMVAWeightPath = TMVAWeightPath;
ZAI U8(92,238,0);
ZAI U5(92,235,0);
double U5_enrich= 0.0025;
ZAI Pu8(94,238,0);
ZAI Pu9(94,239,0);
ZAI Pu0(94,240,0);
ZAI Pu1(94,241,0);
ZAI Pu2(94,242,0);
fFissileList = Pu8*1+Pu9*1+Pu0*1+Pu1*1+Pu2*1;
fFertileList = U5*U5_enrich + U8*(1-U5_enrich);
INFO<<"__An equivalence model of PWR MOX has been define__"<<endl;
INFO<<"\tThis model is based on a multi layer perceptron"<<endl;
INFO<<"\t\tThe TMVA weight file is :"<<endl;
INFO<<"\t\t\t"<<fTMVAWeightPath<<endl;
}
//________________________________________________________________________
EQM_MLP_MOX::EQM_MLP_MOX(CLASSLogger* log, string TMVAWeightPath):EquivalenceModel(log)
{
fTMVAWeightPath = TMVAWeightPath;
ZAI U8(92,238,0);
ZAI U5(92,235,0);
double U5_enrich= 0.0025;
ZAI Pu8(94,238,0);
ZAI Pu9(94,239,0);
ZAI Pu0(94,240,0);
ZAI Pu1(94,241,0);
ZAI Pu2(94,242,0);
fFissileList = Pu8*1+Pu9*1+Pu0*1+Pu1*1+Pu2*1;
fFertileList = U5*U5_enrich + U8*(1-U5_enrich);
INFO<<"__An equivalence model of PWR MOX has been define__"<<endl;
INFO<<"\tThis model is based on a multi layer perceptron"<<endl;
INFO<<"\t\tThe TMVA weight file is :"<<endl;
INFO<<"\t\t\t"<<fTMVAWeightPath<<endl;
}
//________________________________________________________________________
TTree* EQM_MLP_MOX::CreateTMVAInputTree(IsotopicVector Fissil,IsotopicVector Fertil,double BurnUp)
{
TTree* InputTree = new TTree("EQTMP", "EQTMP");
float Pu8 = 0;
float Pu9 = 0;
float Pu10 = 0;
float Pu11 = 0;
float Pu12 = 0;
float Am1 = 0;
float U5_enrichment = 0;
float BU = 0;
InputTree->Branch( "Pu8" ,&Pu8 ,"Pu8/F" );
InputTree->Branch( "Pu9" ,&Pu9 ,"Pu9/F" );
InputTree->Branch( "Pu10" ,&Pu10 ,"Pu10/F" );
InputTree->Branch( "Pu11" ,&Pu11 ,"Pu11/F" );
InputTree->Branch( "Pu12" ,&Pu12 ,"Pu12/F" );
InputTree->Branch( "Am1" ,&Am1 ,"Am1/F" );
InputTree->Branch( "U5_enrichment" ,&U5_enrichment ,"U5_enrichment/F" );
InputTree->Branch( "BU" ,&BU ,"BU/F" );
float U8 = Fertil.GetZAIIsotopicQuantity(92,238,0);
float U5 = Fertil.GetZAIIsotopicQuantity(92,235,0);
float U4 = Fertil.GetZAIIsotopicQuantity(92,234,0);
float UTOT = U8 + U5 + U4;
Pu8 = Fissil.GetZAIIsotopicQuantity(94,238,0);
Pu9 = Fissil.GetZAIIsotopicQuantity(94,239,0);
Pu10 = Fissil.GetZAIIsotopicQuantity(94,240,0);
Pu11 = Fissil.GetZAIIsotopicQuantity(94,241,0);
Pu12 = Fissil.GetZAIIsotopicQuantity(94,242,0);
Am1 = Fissil.GetZAIIsotopicQuantity(95,241,0);
double TOTPU=(Pu8+Pu9+Pu10+Pu11+Pu12+Am1);
Pu8 = Pu8 / TOTPU;
Pu9 = Pu9 / TOTPU;
Pu10= Pu10 / TOTPU;
Pu11= Pu11 / TOTPU;
Pu12= Pu12 / TOTPU;
Am1 = Am1 / TOTPU;
U5_enrichment = U5 / UTOT;
BU=BurnUp;
//cout<<"Pu8 "<<Pu8 <<" Pu9 "<< Pu9 <<" Pu10 "<< Pu10 << " Pu11 "<< Pu11 <<" Pu12 "<<Pu12 <<" Am1 "<<Am1<<endl;
//cout<<"BU "<<BU<<" U5_enrichment "<<U5_enrichment<<endl;
if(Pu8 + Pu9 + Pu10 + Pu11 + Pu12 + Am1 > 1.00001 )//?????1.00001??? I don't know it! goes in condition if =1 !! may be float/double issue ...
{
ERROR << Pu8 << " " << Pu9 << " " << Pu10 << " " << Pu11 << " " << Pu12 << " " << Am1 << endl;
exit(0);
}
// All value are molar (!weight)
InputTree->Fill();
return InputTree;
}
//________________________________________________________________________
double EQM_MLP_MOX::ExecuteTMVA(TTree* theTree)
{
// --- Create the Reader object
TMVA::Reader *reader = new TMVA::Reader( "Silent" );
// Create a set of variables and declare them to the reader
// - the variable names MUST corresponds in name and type to those given in the weight file(s) used
Float_t Pu8,Pu9,Pu10,Pu11,Pu12,Am1,BU,U5_enrichment;
reader->AddVariable( "BU" ,&BU );
reader->AddVariable( "U5_enrichment",&U5_enrichment );
reader->AddVariable( "Pu8" ,&Pu8 );
reader->AddVariable( "Pu9" ,&Pu9 );
reader->AddVariable( "Pu10" ,&Pu10);
reader->AddVariable( "Pu11" ,&Pu11);
reader->AddVariable( "Pu12" ,&Pu12);
reader->AddVariable( "Am1" ,&Am1 );
// --- Book the MVA methods
// Book method MLP
TString methodName = "MLP method";
reader->BookMVA( methodName, fTMVAWeightPath );
//theTree->SetBranchAddress("teneur",&teneur);
theTree->SetBranchAddress( "BU" ,&BU );
theTree->SetBranchAddress( "U5_enrichment" ,&U5_enrichment ) ;
theTree->SetBranchAddress( "Pu8" ,&Pu8 );
theTree->SetBranchAddress( "Pu9" ,&Pu9 );
theTree->SetBranchAddress( "Pu10" ,&Pu10 );
theTree->SetBranchAddress( "Pu11" ,&Pu11 );
theTree->SetBranchAddress( "Pu12" ,&Pu12 );
theTree->SetBranchAddress( "Am1" ,&Am1 );
theTree->GetEntry(0);
Float_t val = (reader->EvaluateRegression( methodName ))[0];
delete reader;
delete theTree;
return (double)val; //retourne teneur
}
//________________________________________________________________________
double EQM_MLP_MOX::GetFissileMolarFraction(IsotopicVector Fissil,IsotopicVector Fertil,double BurnUp)
{DBGL
return ExecuteTMVA(CreateTMVAInputTree(Fissil,Fertil,BurnUp));
}
source/branches/CLASSV3/src/EQM_QUAD_PWR_MOX.cxx deleted 100644 → 0
+ 0
− 133
View file @ 58e2418f
#include "EQM_QUAD_PWR_MOX.hxx"
#include <vector>
#include "StringLine.hxx"
#include "CLASSLogger.hxx"
EQM_QUAD_PWR_MOX::EQM_QUAD_PWR_MOX(string WeightPath):EquivalenceModel(new CLASSLogger("EQM_QUAD_PWR_MOX.log"))
{
fWeightPath = WeightPath;
ifstream DataDB(fWeightPath.c_str()); // Open the File
if(!DataDB)
WARNING << " Can't open \"" << fWeightPath << "\"" << endl;
string line;
int start = 0; // First Get Fuel Parameter
getline(DataDB, line);
if( StringLine::NextWord(line, start, ' ') != "PARAM")
{
ERROR << "Bad Database file : " << fWeightPath << " Can't find the Parameter of the DataBase " << endl;
exit (1);
}
while(start < (int)line.size())
fFuelParameter.push_back(atof(StringLine::NextWord(line, start, ' ').c_str()));
INFO << fFuelParameter.size() << " have been read " << endl;
ZAI U8(92,238,0);
ZAI U5(92,235,0);
double U5_enrich= 0.0025;
fFertileList = U5*U5_enrich + U8*(1-U5_enrich);
ZAI Pu8(94,238,0);
ZAI Pu9(94,239,0);
ZAI Pu0(94,240,0);
ZAI Pu1(94,241,0);
ZAI Pu2(94,242,0);
fFissileList = Pu8*1+Pu9*1+Pu0*1+Pu1*1+Pu2*1;
}
EQM_QUAD_PWR_MOX::EQM_QUAD_PWR_MOX(CLASSLogger* log, string WeightPath):EquivalenceModel(log)
{
fWeightPath = WeightPath;
ifstream DataDB(fWeightPath.c_str()); // Open the File
if(!DataDB)
WARNING << " Can't open \"" << fWeightPath << "\"" << endl;
string line;
int start = 0; // First Get Fuel Parameter
getline(DataDB, line);
if( StringLine::NextWord(line, start, ' ') != "PARAM")
{
ERROR << "Bad Database file : " << fWeightPath << " Can't find the Parameter of the DataBase " << endl;
exit (1);
}
while(start < (int)line.size())
fFuelParameter.push_back(atof(StringLine::NextWord(line, start, ' ').c_str()));
INFO << fFuelParameter.size() << " have been read " << endl;
ZAI U8(92,238,0);
ZAI U5(92,235,0);
double U5_enrich= 0.0025;
fFertileList = U5*U5_enrich + U8*(1-U5_enrich);
ZAI Pu8(94,238,0);
ZAI Pu9(94,239,0);
ZAI Pu0(94,240,0);
ZAI Pu1(94,241,0);
ZAI Pu2(94,242,0);
fFissileList = Pu8*1+Pu9*1+Pu0*1+Pu1*1+Pu2*1;
}
EQM_QUAD_PWR_MOX::~EQM_QUAD_PWR_MOX()
{
}
double EQM_QUAD_PWR_MOX::GetFissileMolarFraction(IsotopicVector Fissile,IsotopicVector Fertile,double BurnUp)
{
ZAI ZAIList[6] = {ZAI(94,238,0), ZAI(94,239,0), ZAI(94,240,0), ZAI(94,241,0), ZAI(94,242,0), ZAI(95,241,0) };
vector<double> PuCompo;
double Sum = Fissile.GetSumOfAll();
for(int i = 0; i< 5; i++)
PuCompo.push_back( Fissile.GetZAIIsotopicQuantity(ZAIList[i])/Sum *100 );
PuCompo[2] += Fissile.GetZAIIsotopicQuantity(ZAIList[5])/Sum *100;
double A = 0;
if(PuCompo[0] <= PuCompo[2] && PuCompo[0] <= PuCompo[4] && PuCompo[1] + PuCompo[3] >= 40 && PuCompo[1] >0 )
{
int par = 0;
for(int j = 0 ; j < 5 ; j++)
{
A += fFuelParameter[par] * PuCompo[j]/100 ;
par++;
for(int i = j ; i < 5 ; i++)
{
A += fFuelParameter[par] *PuCompo[i]/100 *PuCompo[j]/100;
par++;
}
}
A += fFuelParameter[par];
}
return A;
}
source/branches/CLASSV3/src/IM_Matrix.cxx deleted 100644 → 0
+ 0
− 315
View file @ 58e2418f
//
// IM_Matrix.cxx
// CLASSSource
//
// Created by BaM on 04/05/2014.
// Copyright (c) 2014 BaM. All rights reserved.
//
#include "IM_Matrix.hxx"
#include "IsotopicVector.hxx"
#include "CLASSLogger.hxx"
#include "StringLine.hxx"
#include <TGraph.h>
#include <TString.h>
#include <sstream>
#include <string>
#include <iostream>
#include <fstream>
#include <algorithm>
#include <cmath>
using namespace std;
//________________________________________________________________________
IM_Matrix::IM_Matrix():IrradiationModel(new CLASSLogger("IM_Matrix.log")), DynamicalSystem()
{
fShorstestHalflife = 3600.*24*160.; //cut by default all nuclei with a shorter liftime than the Cm242 -> remain 33 actinides
}
IM_Matrix::IM_Matrix(CLASSLogger* log):IrradiationModel(log), DynamicalSystem()
{
fShorstestHalflife = 3600.*24*160.; //cut by default all nuclei with a shorter liftime than the Cm242 -> remain 33 actinides
}
//________________________________________________________________________
/* Evolution Calculation */
//________________________________________________________________________
EvolutionData IM_Matrix::GenerateEvolutionData(IsotopicVector isotopicvector, EvolutionData XSSet, double Power, double cycletime)
{
DBGL
if(fFastDecay.size() == 0)
{
BuildDecayMatrix();
fNVar = findex_inver.size();
}
string ReactorType;
vector< TMatrixT<double> > NMatrix ;// TMatrixT<double>(decayindex.size(),1))
{ // Filling the t=0 State;
map<ZAI, double > isotopicquantity = isotopicvector.GetIsotopicQuantity();
TMatrixT<double> N_0Matrix = TMatrixT<double>( findex.size(),1) ;
for(int i = 0; i < (int)findex.size(); i++)
N_0Matrix[i] = 0;
map<ZAI, double >::iterator it ;
for(int i = 0; i < (int)findex.size(); i++)
N_0Matrix[i] = 0;
for(it = isotopicquantity.begin(); it != isotopicquantity.end(); it++)
{
/// Need TO change with FP managment
map<ZAI, int >::iterator it2;
if( (*it).first.Z() < fZAIThreshold )
it2 = findex_inver.find( ZAI(-2,-2,-2) );
else it2 = findex_inver.find( (*it).first );
if(it2 == findex_inver.end() ) //If not in index should be TMP, can't be fast decay for new Fuel !!!
it2 = findex_inver.find( ZAI(-3,-3,-3) );
N_0Matrix[ (*it2).second ][0] = (*it).second ;
}
isotopicquantity.clear();
NMatrix.push_back(N_0Matrix);
N_0Matrix.Clear();
}
//-------------------------//
//--- Perform Evolution ---//
//-------------------------//
ReactorType = XSSet.GetReactorType();
double M_ref = XSSet.GetHeavyMetalMass();
double M = 0;
double Power_ref = XSSet.GetPower();
// Get the mass of the fuel to irradiate in order to perform the evolution at fixed burnup (the burnup step of the calculation will match the burnup step of the XSSet
{
map<ZAI, double >::iterator it ;
IsotopicVector IVtmp = isotopicvector.GetActinidesComposition();
map<ZAI, double >isotopicquantity = IVtmp.GetIsotopicQuantity();
for( it = isotopicquantity.begin(); it != isotopicquantity.end(); it++ )
M += isotopicvector.GetActinidesComposition().GetZAIIsotopicQuantity( (*it).first )*cZAIMass.fZAIMass.find( (*it).first )->second/AVOGADRO*1e-6;
isotopicquantity.clear();
}
int NStep = XSSet.GetFissionXS().begin()->second->GetN();
double* DBTimeStep = XSSet.GetFissionXS().begin()->second->GetX();
int InsideStep = 10;
double timevector[NStep];
timevector[0] = 0;
double Flux[NStep];
TMatrixT<double> FissionEnergy = TMatrixT<double>(findex.size(),1);
for(int i = 0; i < (int)findex.size(); i++)
FissionEnergy[i] = 0;
{
map< ZAI, int >::iterator it;
for(it = findex_inver.begin(); it != findex_inver.end(); it++)
{
map< ZAI, double >::iterator it2 = fFissionEnergy.find(it->first);
if(it2 == fFissionEnergy.end())
{
if(it->first.Z() > fZAIThreshold)
FissionEnergy[it->second][0] = 1.9679e6*it->first.A()-2.601e8; // //simple linear fit to known values ;extrapolation to unknown isotopes
else FissionEnergy[it->second][0] = 0;
}
else
FissionEnergy[it->second][0] = it2->second;
}
}
vector< TMatrixT<double> > FissionXSMatrix; // Store The Fisison XS Matrix
vector< TMatrixT<double> > CaptureXSMatrix; // Store The Capture XS Matrix
vector< TMatrixT<double> > n2nXSMatrix; // Store The n2N XS Matrix
DBGL
for(int i = 0; i < NStep-1; i++)
{
double TStepMax = ( (DBTimeStep[i+1]-DBTimeStep[i] ) ) * Power_ref/M_ref / Power*M ; // Get the next Time step
TMatrixT<double> BatemanMatrix = TMatrixT<double>(findex.size(),findex.size());
TMatrixT<double> BatemanReactionMatrix = TMatrixT<double>(findex.size(),findex.size());
TMatrixT<double> NEvolutionMatrix = TMatrixT<double>(findex.size(),1);
NEvolutionMatrix = NMatrix.back();
FissionXSMatrix.push_back(GetFissionXsMatrix(XSSet, DBTimeStep[i])); //Feel the fission reaction Matrix
CaptureXSMatrix.push_back(GetCaptureXsMatrix(XSSet, DBTimeStep[i])); //Feel the capture reaction Matrix
n2nXSMatrix.push_back(Getn2nXsMatrix(XSSet, DBTimeStep[i])); //Feel the (n,2n) reaction Matrix
// ---------------- Evolution
BatemanReactionMatrix = FissionXSMatrix[i];
BatemanReactionMatrix += CaptureXSMatrix[i];
BatemanReactionMatrix += n2nXSMatrix[i];
for(int k=0; k < InsideStep; k++)
{
double ESigmaN = 0;
for (int j = 0; j < (int)findex.size() ; j++)
ESigmaN -= FissionXSMatrix[i][j][j]*NEvolutionMatrix[j][0]*1.6e-19*FissionEnergy[j][0];
// Update Flux
double Flux_k = Power/ESigmaN;
if(k==0)
Flux[i]=Flux_k;
BatemanMatrix = BatemanReactionMatrix;
BatemanMatrix *= Flux_k;
BatemanMatrix += fDecayMatrix ;
BatemanMatrix *= TStepMax/InsideStep ;
TMatrixT<double> IdMatrix = TMatrixT<double>(findex.size(),findex.size());
for(int j = 0; j < (int)findex.size(); j++)
for(int k = 0; k < (int)findex.size(); k++)
{
if(k == j) IdMatrix[j][k] = 1;
else IdMatrix[j][k] = 0;
}
TMatrixT<double> BatemanMatrixDL = TMatrixT<double>(findex.size(),findex.size()); // Order 0 Term from the DL : Id
TMatrixT<double> BatemanMatrixDLTermN = TMatrixT<double>(findex.size(),findex.size()); // Addind it;
{
BatemanMatrix *= TStepMax ;
BatemanMatrixDLTermN = IdMatrix;
BatemanMatrixDL = BatemanMatrixDLTermN;
int j = 1;
double NormN;
do
{
TMatrixT<double> BatemanMatrixDLTermtmp = TMatrixT<double>(findex.size(),findex.size()); // Adding it;
BatemanMatrixDLTermtmp = BatemanMatrixDLTermN;
BatemanMatrixDLTermN.Mult(BatemanMatrixDLTermtmp, BatemanMatrix );
BatemanMatrixDLTermN *= 1./j;
BatemanMatrixDL += BatemanMatrixDLTermN;
NormN = 0;
for(int m = 0; m < (int)findex.size(); m++)
for(int n = 0; n < (int)findex.size(); n++)
NormN += BatemanMatrixDLTermN[m][n]*BatemanMatrixDLTermN[m][n];
j++;
} while ( NormN != 0 );
}
NEvolutionMatrix = BatemanMatrixDL * NEvolutionMatrix ;
}
NMatrix.push_back(NEvolutionMatrix);
timevector[i+1] = timevector[i] + TStepMax;
BatemanMatrix.Clear();
BatemanReactionMatrix.Clear();
NEvolutionMatrix.Clear();
}
DBGL
FissionXSMatrix.push_back(GetFissionXsMatrix(XSSet, DBTimeStep[NStep-1])); //Feel the reaction Matrix
CaptureXSMatrix.push_back(GetCaptureXsMatrix(XSSet, DBTimeStep[NStep-1])); //Feel the reaction Matrix
n2nXSMatrix.push_back(Getn2nXsMatrix(XSSet, DBTimeStep[NStep-1])); //Feel the reaction Matrix
EvolutionData GeneratedDB = EvolutionData(GetLog());
double ESigmaN = 0;
for (int j = 0; j < (int)findex.size() ; j++)
ESigmaN -= FissionXSMatrix.back()[j][j]*NMatrix.back()[j][0]*1.6e-19*FissionEnergy[j][0];
Flux[NStep-1] = Power/ESigmaN;
GeneratedDB.SetFlux( new TGraph(NStep, timevector, Flux) );
for(int i = 0; i < (int)findex.size(); i++)
{
double ZAIQuantity[NMatrix.size()];
double FissionXS[NStep];
double CaptureXS[NStep];
double n2nXS[NStep];
for(int j = 0; j < (int)NMatrix.size(); j++)
ZAIQuantity[j] = (NMatrix[j])[i][0];
for(int j = 0; j < NStep; j++)
{
FissionXS[j] = FissionXSMatrix[j][i][i];
CaptureXS[j] = CaptureXSMatrix[j][i][i];
n2nXS[j] = n2nXSMatrix[j][i][i];
}
GeneratedDB.NucleiInsert(pair<ZAI, TGraph*> (findex.find(i)->second, new TGraph(NMatrix.size(), timevector, ZAIQuantity)));
GeneratedDB.FissionXSInsert(pair<ZAI, TGraph*> (findex.find(i)->second, new TGraph(NStep, timevector, FissionXS)));
GeneratedDB.CaptureXSInsert(pair<ZAI, TGraph*> (findex.find(i)->second, new TGraph(NStep, timevector, CaptureXS)));
GeneratedDB.n2nXSInsert(pair<ZAI, TGraph*> (findex.find(i)->second, new TGraph(NStep, timevector, n2nXS)));
}
GeneratedDB.SetPower(Power );
GeneratedDB.SetReactorType(ReactorType );
GeneratedDB.SetCycleTime(cycletime);
for (int i = 0; i < (int) FissionXSMatrix.size(); i++)
{
FissionXSMatrix[i].Clear();
CaptureXSMatrix[i].Clear();
n2nXSMatrix[i].Clear();
}
FissionXSMatrix.clear();
CaptureXSMatrix.clear();
n2nXSMatrix.clear();
DBGL
return GeneratedDB;
}
source/branches/CLASSV3/src/IM_RK4.cxx deleted 100644 → 0
+ 0
− 408
View file @ 58e2418f
//
// IM_RK4.cxx
// CLASSSource
//
// Created by BaM on 04/05/2014.
// Copyright (c) 2014 BaM. All rights reserved.
//
#include "IM_RK4.hxx"
#include "IsotopicVector.hxx"
#include "CLASSConstante.hxx"
#include "CLASSLogger.hxx"
#include "StringLine.hxx"
#include <TGraph.h>
#include <TString.h>
#include <sstream>
#include <string>
#include <iostream>
#include <fstream>
#include <algorithm>
#include <cmath>
using namespace std;
//________________________________________________________________________
IM_RK4::IM_RK4():IrradiationModel(new CLASSLogger("IM_RK4.log")), DynamicalSystem()
{
fTheNucleiVector = 0;
fTheMatrix = 0;
SetForbidNegativeValue();
}
IM_RK4::IM_RK4(CLASSLogger* log):IrradiationModel(log), DynamicalSystem()
{
fTheNucleiVector = 0;
fTheMatrix = 0;
SetForbidNegativeValue();
}
//________________________________________________________________________
/* Evolution Calculation */
//________________________________________________________________________
EvolutionData IM_RK4::GenerateEvolutionData(IsotopicVector isotopicvector, EvolutionData XSSet, double Power, double cycletime)
{
DBGL
if(fFastDecay.size() == 0)
{
BuildDecayMatrix();
fNVar = findex_inver.size();
}
SetTheMatrixToZero();
SetTheNucleiVectorToZero();
string ReactorType;
vector< TMatrixT<double> > NMatrix ;// TMatrixT<double>(decayindex.size(),1))
{ // Filling the t=0 State;
map<ZAI, double > isotopicquantity = isotopicvector.GetIsotopicQuantity();
TMatrixT<double> N_0Matrix = TMatrixT<double>( findex.size(),1) ;
for(int i = 0; i < (int)findex.size(); i++)
N_0Matrix[i] = 0;
map<ZAI, double >::iterator it ;
for(int i = 0; i < (int)findex.size(); i++)
N_0Matrix[i] = 0;
for(it = isotopicquantity.begin(); it != isotopicquantity.end(); it++)
{
/// Need TO change with FP managment
map<ZAI, int >::iterator it2;
if( (*it).first.Z() < fZAIThreshold )
it2 = findex_inver.find( ZAI(-2,-2,-2) );
else it2 = findex_inver.find( (*it).first );
if(it2 == findex_inver.end() ) //If not in index should be TMP, can't be fast decay for new Fuel !!!
it2 = findex_inver.find( ZAI(-3,-3,-3) );
N_0Matrix[ (*it2).second ][0] = (*it).second ;
}
isotopicquantity.clear();
NMatrix.push_back(N_0Matrix);
N_0Matrix.Clear();
}
//-------------------------//
//--- Perform Evolution ---//
//-------------------------//
ReactorType = XSSet.GetReactorType();
double M_ref = XSSet.GetHeavyMetalMass();
double M = 0;
double Power_ref = XSSet.GetPower();
// Get the mass of the fuel to irradiate in order to perform the evolution at fixed burnup (the burnup step of the calculation will match the burnup step of the XSSet
{
map<ZAI, double >::iterator it ;
IsotopicVector IVtmp = isotopicvector.GetActinidesComposition();
map<ZAI, double >isotopicquantity = IVtmp.GetIsotopicQuantity();
for( it = isotopicquantity.begin(); it != isotopicquantity.end(); it++ )
M += isotopicvector.GetActinidesComposition().GetZAIIsotopicQuantity( (*it).first )*cZAIMass.fZAIMass.find( (*it).first )->second/AVOGADRO*1e-6;
isotopicquantity.clear();
}
int NStep = XSSet.GetFissionXS().begin()->second->GetN();
double* DBTimeStep = XSSet.GetFissionXS().begin()->second->GetX();
int InsideStep = 10;
double timevector[NStep];
timevector[0] = 0;
double Flux[NStep];
TMatrixT<double> FissionEnergy = TMatrixT<double>(findex.size(),1);
for(int i = 0; i < (int)findex.size(); i++)
FissionEnergy[i] = 0;
{
map< ZAI, int >::iterator it;
for(it = findex_inver.begin(); it != findex_inver.end(); it++)
{
map< ZAI, double >::iterator it2 = fFissionEnergy.find(it->first);
if(it2 == fFissionEnergy.end())
{
if(it->first.Z() > fZAIThreshold)
FissionEnergy[it->second][0] = 1.9679e6*it->first.A()-2.601e8; // //simple linear fit to known values ;extrapolation to unknown isotopes
else FissionEnergy[it->second][0] = 0;
}
else
FissionEnergy[it->second][0] = it2->second;
}
}
vector< TMatrixT<double> > FissionXSMatrix; // Store The Fisison XS Matrix
vector< TMatrixT<double> > CaptureXSMatrix; // Store The Capture XS Matrix
vector< TMatrixT<double> > n2nXSMatrix; // Store The n2N XS Matrix
DBGL
for(int i = 0; i < NStep-1; i++)
{
double TStepMax = ( (DBTimeStep[i+1]-DBTimeStep[i] ) ) * Power_ref/M_ref / Power*M ; // Get the next Time step
TMatrixT<double> BatemanMatrix = TMatrixT<double>(findex.size(),findex.size());
TMatrixT<double> BatemanReactionMatrix = TMatrixT<double>(findex.size(),findex.size());
TMatrixT<double> NEvolutionMatrix = TMatrixT<double>(findex.size(),1);
NEvolutionMatrix = NMatrix.back();
FissionXSMatrix.push_back(GetFissionXsMatrix(XSSet, DBTimeStep[i])); //Feel the fission reaction Matrix
CaptureXSMatrix.push_back(GetCaptureXsMatrix(XSSet, DBTimeStep[i])); //Feel the capture reaction Matrix
n2nXSMatrix.push_back(Getn2nXsMatrix(XSSet, DBTimeStep[i])); //Feel the (n,2n) reaction Matrix
// ---------------- Evolution
BatemanReactionMatrix = FissionXSMatrix[i];
BatemanReactionMatrix += CaptureXSMatrix[i];
BatemanReactionMatrix += n2nXSMatrix[i];
for(int k=0; k < InsideStep; k++)
{
double ESigmaN = 0;
for (int j = 0; j < (int)findex.size() ; j++)
ESigmaN -= FissionXSMatrix[i][j][j]*NEvolutionMatrix[j][0]*1.6e-19*FissionEnergy[j][0];
// Update Flux
double Flux_k = Power/ESigmaN;
if(k==0)
Flux[i]=Flux_k;
BatemanMatrix = BatemanReactionMatrix;
BatemanMatrix *= Flux_k;
BatemanMatrix += fDecayMatrix ;
SetTheMatrixToZero();
SetTheNucleiVectorToZero();
SetTheMatrix(BatemanMatrix);
SetTheNucleiVector(NEvolutionMatrix);
RungeKutta(fTheNucleiVector, timevector[i]+TStepMax/InsideStep*k, timevector[i]+TStepMax/InsideStep*(k+1), fNVar);
NEvolutionMatrix = GetTheNucleiVector();
}
NEvolutionMatrix = GetTheNucleiVector();
NMatrix.push_back(NEvolutionMatrix);
timevector[i+1] = timevector[i] + TStepMax;
BatemanMatrix.Clear();
BatemanReactionMatrix.Clear();
NEvolutionMatrix.Clear();
}
FissionXSMatrix.push_back(GetFissionXsMatrix(XSSet, DBTimeStep[NStep-1])); //Feel the reaction Matrix
CaptureXSMatrix.push_back(GetCaptureXsMatrix(XSSet, DBTimeStep[NStep-1])); //Feel the reaction Matrix
n2nXSMatrix.push_back(Getn2nXsMatrix(XSSet, DBTimeStep[NStep-1])); //Feel the reaction Matrix
EvolutionData GeneratedDB = EvolutionData(GetLog());
double ESigmaN = 0;
for (int j = 0; j < (int)findex.size() ; j++)
ESigmaN -= FissionXSMatrix.back()[j][j]*NMatrix.back()[j][0]*1.6e-19*FissionEnergy[j][0];
Flux[NStep-1] = Power/ESigmaN;
GeneratedDB.SetFlux( new TGraph(NStep, timevector, Flux) );
for(int i = 0; i < (int)findex.size(); i++)
{
double ZAIQuantity[NMatrix.size()];
double FissionXS[NStep];
double CaptureXS[NStep];
double n2nXS[NStep];
for(int j = 0; j < (int)NMatrix.size(); j++)
ZAIQuantity[j] = (NMatrix[j])[i][0];
for(int j = 0; j < NStep; j++)
{
FissionXS[j] = FissionXSMatrix[j][i][i];
CaptureXS[j] = CaptureXSMatrix[j][i][i];
n2nXS[j] = n2nXSMatrix[j][i][i];
}
GeneratedDB.NucleiInsert(pair<ZAI, TGraph*> (findex.find(i)->second, new TGraph(NMatrix.size(), timevector, ZAIQuantity)));
GeneratedDB.FissionXSInsert(pair<ZAI, TGraph*> (findex.find(i)->second, new TGraph(NStep, timevector, FissionXS)));
GeneratedDB.CaptureXSInsert(pair<ZAI, TGraph*> (findex.find(i)->second, new TGraph(NStep, timevector, CaptureXS)));
GeneratedDB.n2nXSInsert(pair<ZAI, TGraph*> (findex.find(i)->second, new TGraph(NStep, timevector, n2nXS)));
}
DBGL
GeneratedDB.SetPower(Power );
GeneratedDB.SetReactorType(ReactorType );
GeneratedDB.SetCycleTime(cycletime);
ResetTheMatrix();
ResetTheNucleiVector();
for (int i = 0; i < (int) FissionXSMatrix.size(); i++)
{
FissionXSMatrix[i].Clear();
CaptureXSMatrix[i].Clear();
n2nXSMatrix[i].Clear();
}
FissionXSMatrix.clear();
CaptureXSMatrix.clear();
n2nXSMatrix.clear();
DBGL
return GeneratedDB;
}
void IM_RK4::ResetTheMatrix()
{
if(fTheMatrix)
{
for(int i= 0; i<fNVar; i++)
delete [] fTheMatrix[i];
delete [] fTheMatrix;
}
fTheMatrix = 0;
}
void IM_RK4::SetTheMatrixToZero()
{
ResetTheMatrix();
fNVar = findex.size();
fTheMatrix = new double*[fNVar];
#pragma omp parallel for
for(int i= 0; i < fNVar; i++)
fTheMatrix[i] = new double[fNVar];
for(int i = 0; i < fNVar; i++)
for(int k = 0; k < fNVar; k++)
{
fTheMatrix[i][k]=0.0;
}
}
//________________________________________________________________________
void IM_RK4::ResetTheNucleiVector()
{
if(fTheNucleiVector)
delete [] fTheNucleiVector;
fTheNucleiVector = 0;
}
//________________________________________________________________________
void IM_RK4::SetTheNucleiVectorToZero()
{
ResetTheNucleiVector();
fTheNucleiVector = new double[fNVar];
#pragma omp parallel for
for(int i = 0; i < fNVar; i++)
fTheNucleiVector[i]=0.0;
}
//________________________________________________________________________
void IM_RK4::BuildEqns(double t, double *N, double *dNdt)
{
double sum=0;
// pragma omp parallel for reduction(+:sum)
for(int i = 0; i < fNVar; i++)
{
sum=0;
for(int k = 0; k < fNVar; k++)
{
sum += fTheMatrix[i][k]*N[k];
}
dNdt[i] = sum;
}
}
//________________________________________________________________________
void IM_RK4::SetTheMatrix(TMatrixT<double> BatemanMatrix)
{
for (int k = 0; k < (int)fNVar; k++)
for (int l = 0; l < (int)findex_inver.size(); l++)
fTheMatrix[l][k] = BatemanMatrix[l][k];
}
//________________________________________________________________________
TMatrixT<double> IM_RK4::GetTheMatrix()
{
TMatrixT<double> BatemanMatrix = TMatrixT<double>(findex.size(),findex.size());
for (int k = 0; k < (int)fNVar; k++)
for (int l = 0; l < (int)findex_inver.size(); l++)
BatemanMatrix[l][k] = fTheMatrix[l][k];
return BatemanMatrix;
}
//________________________________________________________________________
void IM_RK4::SetTheNucleiVector(TMatrixT<double> NEvolutionMatrix)
{
for (int k = 0; k < (int)fNVar; k++)
fTheNucleiVector[k] = NEvolutionMatrix[k][0];
}
//________________________________________________________________________
TMatrixT<double> IM_RK4::GetTheNucleiVector()
{
TMatrixT<double> NEvolutionMatrix = TMatrixT<double>(findex.size(),1);
for (int k = 0; k < (int)fNVar; k++)
NEvolutionMatrix[k][0] = fTheNucleiVector[k];
return NEvolutionMatrix;
}
source/branches/CLASSV3/src/XSM_CLOSEST.cxx deleted 100644 → 0
+ 0
− 366
View file @ 58e2418f
#include "XSModel.hxx"
#include "XSM_CLOSEST.hxx"
#include "CLASSLogger.hxx"
#include "StringLine.hxx"
#include <TGraph.h>
#include <TString.h>
#include "TSystem.h"
#include "TROOT.h"
#include <sstream>
#include <string>
#include <iostream>
#include <fstream>
#include <algorithm>
#include <cmath>
//________________________________________________________________________
//
// XSM_CLOSEST
//
//
//
//
//________________________________________________________________________
XSM_CLOSEST::XSM_CLOSEST(string DB_index_file, bool oldreadmethod ): XSModel(new CLASSLogger("XSM_CLOSEST.log"))
{
fOldReadMethod = oldreadmethod;
fDataBaseIndex = DB_index_file;
fDistanceType = 0;
fWeightedDistance = false;
fEvolutionDataInterpolation = false;
ReadDataBase();
// Warning
INFO << " A EvolutionData has been define :" << endl;
INFO << "\t His index is : \"" << DB_index_file << "\" " << endl;
INFO << "\t " << fFuelDataBank.size() << " EvolutionData have been read."<< endl << endl;
}
//________________________________________________________________________
XSM_CLOSEST::XSM_CLOSEST(CLASSLogger* Log,string DB_index_file, bool oldreadmethod ): XSModel(Log)
{
fOldReadMethod = oldreadmethod;
fDataBaseIndex = DB_index_file;
fDistanceType = 0;
fWeightedDistance = false;
fEvolutionDataInterpolation = false;
ReadDataBase();
// Warning
INFO << " A EvolutionData has been define :" << endl;
INFO << "\t His index is : \"" << DB_index_file << "\" " << endl;
INFO << "\t " << fFuelDataBank.size() << " EvolutionData have been read."<< endl << endl;
}
//________________________________________________________________________
XSM_CLOSEST::~XSM_CLOSEST()
{
for( int i = 0; i < (int)fFuelDataBank.size(); i++)
fFuelDataBank[i].DeleteEvolutionData();
fFuelDataBank.clear();
}
//________________________________________________________________________
void XSM_CLOSEST::ReadDataBase()
{
if(fFuelDataBank.size() != 0)
{
for( int i = 0; i < (int)fFuelDataBank.size(); i++)
fFuelDataBank[i].DeleteEvolutionData();
fFuelDataBank.clear();
}
ifstream DataDB(fDataBaseIndex.c_str()); // Open the File
if(!DataDB)
WARNING << " Can't open \"" << fDataBaseIndex << "\"\n" << endl;
vector<double> vTime;
vector<double> vTimeErr;
string line;
int start = 0;
// First Get Fuel Type
getline(DataDB, line);
if( StringLine::NextWord(line, start, ' ') != "TYPE")
{
ERROR << " Bad Database file : " << fDataBaseIndex << " Can't find the type of the DataBase"<< endl;
exit (1);
}
fFuelType = StringLine::NextWord(line, start, ' ');
//Then Get All the Database
while (!DataDB.eof())
{
getline(DataDB, line);
if(line != "")
{
EvolutionData evolutionproduct(GetLog(), line, fOldReadMethod);
fFuelDataBank.push_back(evolutionproduct);
}
}
}
//________________________________________________________________________
//________________________________________________________________________
//________________________________________________________________________
/* Distance Calculation */
//________________________________________________________________________
//________________________________________________________________________
//________________________________________________________________________
map<double, int> XSM_CLOSEST::GetDistancesTo(IsotopicVector isotopicvector, double t)
{
map<double, int> distances;
for( int i = 0; i < (int)fFuelDataBank.size(); i++)
{
pair<map<double, int>::iterator, bool> IResult;
IsotopicVector ActinidesCompositionAtT = fFuelDataBank[i].GetIsotopicVectorAt(t).GetActinidesComposition();
IsotopicVector IV_ActinidesComposition = isotopicvector.GetActinidesComposition();
double NormalisationFactor = Norme(IV_ActinidesComposition) / Norme( ActinidesCompositionAtT );
double distance = Distance( IV_ActinidesComposition,
ActinidesCompositionAtT / NormalisationFactor,
fDistanceType,
fDistanceParameter);
IResult = distances.insert( pair< double, int >(distance, i) );
}
return distances;
}
//________________________________________________________________________
EvolutionData XSM_CLOSEST::GetCrossSections(IsotopicVector isotopicvector, double t)
{
DBGL
double distance = 0;
int N_BestEvolutionData = 0;
if(fWeightedDistance)
{
DBGL
IsotopicVector ActinidesCompositionAtT = fFuelDataBank[0].GetIsotopicVectorAt(t).GetActinidesComposition();
IsotopicVector IV_ActinidesComposition = isotopicvector.GetActinidesComposition();
double NormalisationFactor = Norme( ActinidesCompositionAtT ) / Norme(IV_ActinidesComposition);
distance = Distance( IV_ActinidesComposition / NormalisationFactor,
fFuelDataBank[0]);
for( int i = 1; i < (int)fFuelDataBank.size(); i++)
{
double D = 0;
ActinidesCompositionAtT = fFuelDataBank[i].GetIsotopicVectorAt(t).GetActinidesComposition();
IV_ActinidesComposition = isotopicvector.GetActinidesComposition();
D = Distance( IV_ActinidesComposition / NormalisationFactor,
fFuelDataBank[i]);
if (D< distance)
{
distance = D;
N_BestEvolutionData = i;
}
}
DBGL
return fFuelDataBank[N_BestEvolutionData];
}
else if (fEvolutionDataInterpolation)
{
DBGL
map<double, int> distance_map = GetDistancesTo(isotopicvector, t);
map<double, int>::iterator it_distance;
int NClose = 64;
int Nstep = 0;
EvolutionData EvolInterpolate;
double SumOfDistance = 0;
for( it_distance = distance_map.begin(); it_distance != distance_map.end(); it_distance++)
{
double distance = (*it_distance).first;
int ED_Indice = (*it_distance).second;
if(distance == 0 )
{
EvolInterpolate = Multiply(1,fFuelDataBank[ED_Indice]);
return EvolInterpolate;
}
if(Nstep == 0)
EvolInterpolate = Multiply(1./distance, fFuelDataBank[ED_Indice]);
else
{
EvolutionData Evoltmp = EvolInterpolate;
EvolutionData Evoltmp2 = Multiply(1./distance, fFuelDataBank[ED_Indice]);
EvolInterpolate = Sum(Evoltmp, Evoltmp2);
Evoltmp.DeleteEvolutionData();
Evoltmp2.DeleteEvolutionData();
}
SumOfDistance += 1./distance;
Nstep++;
if(Nstep == NClose) break;
}
EvolutionData Evoltmp = EvolInterpolate;
EvolInterpolate.Clear();
EvolInterpolate = Multiply(1/SumOfDistance, Evoltmp);
Evoltmp.DeleteEvolutionData();
DBGL
return EvolInterpolate;
}
else
{
DBGL
IsotopicVector ActinidesCompositionAtT = fFuelDataBank[0].GetIsotopicVectorAt(t).GetActinidesComposition();
IsotopicVector IV_ActinidesComposition = isotopicvector.GetActinidesComposition();
double NormalisationFactor = Norme(IV_ActinidesComposition) / Norme( ActinidesCompositionAtT );
distance = Distance( IV_ActinidesComposition,
ActinidesCompositionAtT / NormalisationFactor,
fDistanceType,
fDistanceParameter);
for( int i = 1; i < (int)fFuelDataBank.size(); i++)
{
double D = 0;
ActinidesCompositionAtT = fFuelDataBank[i].GetIsotopicVectorAt(t).GetActinidesComposition();
IV_ActinidesComposition = isotopicvector.GetActinidesComposition();
D = Distance( IV_ActinidesComposition,
ActinidesCompositionAtT / NormalisationFactor,
fDistanceType,
fDistanceParameter);
if (D< distance)
{
distance = D;
N_BestEvolutionData = i;
}
}
DBGL
return fFuelDataBank[N_BestEvolutionData];
}
}
//________________________________________________________________________
void XSM_CLOSEST::CalculateDistanceParameter()
{
DBGL
if(fDistanceType!=1){
WARNING << " Distance Parameter will be calculate even if the distance type is not the good one. Any Distance Parameters given by the user will be overwriten" << endl;
}
fDistanceParameter.Clear();
//We calculate the weight for the distance calculation.
int NevolutionDatainFuelDataBank = 0;
for( int i = 0; i < (int)fFuelDataBank.size(); i++)
{
NevolutionDatainFuelDataBank++;
map<ZAI ,double>::iterator itit;
map<ZAI ,double> isovector = fFuelDataBank[i].GetIsotopicVectorAt(0).GetIsotopicQuantity();
for(itit=isovector.begin(); itit != isovector.end(); itit++) //Boucle sur ZAI
{
double TmpXS=0;
for( int i=1; i<4; i++ ) //Loop on Reactions 1==fission, 2==capture, 3==n2n
TmpXS+= fFuelDataBank[i].GetXSForAt(0, (*itit).first, i);
fDistanceParameter.Add((*itit).first,TmpXS);
}
}
fDistanceParameter.Multiply( (double)1.0/NevolutionDatainFuelDataBank );
if(GetLog())
{
INFO <<"!!INFO!! Distance Parameters "<<endl;
map<ZAI ,double >::iterator it2;
for(it2 = fDistanceParameter.GetIsotopicQuantity().begin();it2 != fDistanceParameter.GetIsotopicQuantity().end(); it2++)
{
INFO << (*it2).first.Z() << " ";
INFO << (*it2).first.A() << " ";
INFO << (*it2).first.I() << " ";
INFO << ": " << (*it2).second;
INFO << endl;
}
INFO << endl;
}
DBGL
}
//________________________________________________________________________
void XSM_CLOSEST::SetDistanceParameter(IsotopicVector DistanceParameter)
{
fDistanceParameter = DistanceParameter;
INFO <<"!!INFO!! Distance Parameters "<<endl;
map<ZAI ,double >::iterator it2;
for(it2 = fDistanceParameter.GetIsotopicQuantity().begin();it2 != fDistanceParameter.GetIsotopicQuantity().end(); it2++)
{
INFO << (*it2).first.Z() << " ";
INFO << (*it2).first.A() << " ";
INFO << (*it2).first.I() << " ";
INFO << ": " << (*it2).second;
INFO << endl;
}
INFO << endl;
}
//________________________________________________________________________
void XSM_CLOSEST::SetDistanceType(int DistanceType)
{
fDistanceType=DistanceType;
if(fDistanceType==1){
CalculateDistanceParameter();
}
else if(fDistanceType == 2 && Norme(fDistanceParameter)==0){
// This is so bad!! You will probably unsynchronize all the reactor....
WARNING << " Distance use weight defined by user for each isotope, but no weight have been given" << endl<<"Use SetDistanceParameter()"<<endl;
exit(1);
}
else if (fDistanceType != 0 && fDistanceType != 1 && fDistanceType != 2 ){
ERROR << " Distancetype defined by the user isn't recognized by the code"<<endl;
exit(1);
}
}
\ No newline at end of file
source/branches/CLASSV3/src/XSM_MLP.cxx deleted 100644 → 0
+ 0
− 509
View file @ 58e2418f
#include "XSModel.hxx"
#include "XSM_MLP.hxx"
#include "CLASSLogger.hxx"
#include "StringLine.hxx"
#include "TMVA/Reader.h"
#include "TMVA/Tools.h"
#include "TMVA/MethodCuts.h"
#include <TGraph.h>
#include <TString.h>
#include "TFile.h"
#include "TSystem.h"
#include "TROOT.h"
#include "TStopwatch.h"
#include <dirent.h>
#include <errno.h>
#include <sstream>
#include <string>
#include <iostream>
#include <fstream>
#include <algorithm>
#include <cmath>
//________________________________________________________________________
//
// XSM_MLP
//
//
//
//
//________________________________________________________________________
XSM_MLP::XSM_MLP(string TMVA_Weight_Directory,string InformationFile, bool IsTimeStep):XSModel(new CLASSLogger("XSM_MLP.log"))
{
fIsStepTime=IsTimeStep;
fTMVAWeightFolder = TMVA_Weight_Directory;
if(InformationFile=="")
fMLPInformationFile = TMVA_Weight_Directory+"/Data_Base_Info.nfo";
else
fMLPInformationFile=fTMVAWeightFolder+InformationFile;
GetMLPWeightFiles();
GetDataBaseInformation();
INFO<<"__A cross section interpolator using" <<endl;
INFO<<"Multi Layer Perceptron has been define__"<<endl;
INFO << " \t His TMVA folder is : \" " << fTMVAWeightFolder << "\"" << endl;
}
//________________________________________________________________________
XSM_MLP::XSM_MLP(CLASSLogger* Log,string TMVA_Weight_Directory,string InformationFile, bool IsTimeStep):XSModel(Log)
{
fIsStepTime=IsTimeStep;
fTMVAWeightFolder = TMVA_Weight_Directory;
if(InformationFile=="")
fMLPInformationFile = TMVA_Weight_Directory+"/Data_Base_Info.nfo";
else
fMLPInformationFile=fTMVAWeightFolder+InformationFile;
GetMLPWeightFiles();
GetDataBaseInformation();
INFO<<"__A cross section interpolator using" <<endl;
INFO<<"Multi Layer Perceptron has been define__"<<endl;
INFO << " \t His TMVA folder is : \" " << fTMVAWeightFolder << "\"" << endl;
}
//________________________________________________________________________
XSM_MLP::~XSM_MLP()
{
fMapOfTMVAVariableNames.clear();
}
//________________________________________________________________________
void XSM_MLP::GetDataBaseInformation()
{
ifstream FILE(fMLPInformationFile.c_str());
if(FILE.good())
{
while(!FILE.eof())
{
string line;
getline(FILE, line);
size_t foundRType = line.find("ReactorType :");
size_t foundFType = line.find("FuelType :");
size_t foundHM = line.find("Heavy Metal (t) :");
size_t foundPower = line.find("Thermal Power (W) :");
size_t foundTime = line.find("Time (s) :");
size_t foundZAI = line.find("Z A I Name (input MLP) :");
int pos=0;
if(foundRType != std::string::npos)
{ StringLine::NextWord(line,pos,':');
fDataBaseRType = atof( (StringLine::NextWord(line,pos,':')).c_str() );
}
if(foundFType != std::string::npos)
{ StringLine::NextWord(line,pos,':');
fDataBaseFType = atof( (StringLine::NextWord(line,pos,':')).c_str() );
}
if(foundHM != std::string::npos)
{ StringLine::NextWord(line,pos,':');
fDataBaseHMMass = atof( (StringLine::NextWord(line,pos,':')).c_str() );
}
if(foundPower !=std::string::npos)
{ StringLine::NextWord(line,pos,':');
fDataBasePower = atof( (StringLine::NextWord(line,pos,':') ).c_str() );
}
pos=0;
if(foundTime!=std::string::npos)
{
StringLine::NextWord(line,pos,':');
while( pos< (int)line.size() )
fMLP_Time.push_back( atof( (StringLine::NextWord(line,pos,' ')).c_str() ));
}
if(foundZAI != std::string::npos)
{
while(!FILE.eof())
{
int Z=-4;
int A=-4;
int I=-4;
string Name;
getline(FILE, line);
stringstream ssline;
ssline<<line;
ssline>>Z>>A>>I>>Name;
//cout<<Z<<" "<<A<<" "<<I<<" "<<Name<<endl;
fMapOfTMVAVariableNames.insert( pair<ZAI,string>(ZAI(Z,A,I),Name) );
}
}
}
}
else
{
ERROR << "Can't find/open file " << fMLPInformationFile << endl;
exit(0);
}
/********DEBUG*************************************/
INFO<<"\tMLP XS Data Base Information : "<<endl;
INFO<<"\t\tHeavy Metal (t) :"<<fDataBaseHMMass<<endl;
INFO<<"\t\tThermal Power (W) :"<<fDataBasePower<<endl;
INFO<<"\t\tTime (s) :"<<endl;
for (int i = 0; i < (int)fMLP_Time.size(); ++i)
INFO<<"\t\t\t"<<fMLP_Time[i]<<endl;
INFO<<"\t\tZ A I Name (input MLP) :"<<endl;
map<ZAI ,string >::iterator it;
for (it= fMapOfTMVAVariableNames.begin();it!=fMapOfTMVAVariableNames.end();it++)
INFO<<"\t\t\t"<< it->first.Z()<<" "<<it->first.A()<<" "<<it->second<<endl;
}
//________________________________________________________________________
void XSM_MLP::GetMLPWeightFiles()
{DBGL
/**********Get All TMVA weight files*******************/
//check if the folder containing weights exists
DIR* rep = NULL;
struct dirent* fichierLu = NULL;
rep = opendir(fTMVAWeightFolder.c_str());
if (rep == NULL)
{
perror("");
exit(1);
}
/**Save file names of TMVA weights*/
fWeightFiles.resize(0);
while ((fichierLu = readdir(rep)) != NULL)
{
string FileName= fichierLu->d_name ;
if(FileName != "." && FileName != ".." )
{
if(FileName[FileName.size()-3]=='x' && FileName[FileName.size()-2]=='m' && FileName[FileName.size()-1]=='l' )
fWeightFiles.push_back(FileName);
}
}
DBGL
}
//________________________________________________________________________
//________________________________________________________________________
//
// Time (MLP take time as parameter)
//________________________________________________________________________
//________________________________________________________________________
void XSM_MLP::ReadWeightFile(string Filename, int &Z, int &A, int &I, int &Reaction)
{
Z=-1;
A=-1;
I=-1;
Reaction=-1;
size_t found = Filename.find("XS");
string NameJOB;
NameJOB=Filename.substr(found);
int pos=0;
StringLine::NextWord(NameJOB, pos, '_');
Z = atof( (StringLine::NextWord(NameJOB,pos,'_') ).c_str() );
A = atof( (StringLine::NextWord(NameJOB,pos,'_') ).c_str() );
I = atof( (StringLine::NextWord(NameJOB,pos,'_') ).c_str() );
string sReaction = (StringLine::NextWord(NameJOB,pos,'_') ).c_str() ;
size_t foundext = sReaction.find(".weights.xml");
sReaction = sReaction.substr(0,foundext);
if(sReaction=="fis")
Reaction=0;
if(sReaction=="cap")
Reaction=1;
if(sReaction=="n2n")
Reaction=2;
if(Z<=0 || A<=0 || I<0 || Reaction==-1)
{
ERROR << " wrong TMVA weight format " << endl;
exit(0);
}
}
//________________________________________________________________________
TTree* XSM_MLP::CreateTMVAInputTree(IsotopicVector isotopicvector,int TimeStep)
{
/******Create Input data tree to be interpreted by TMVA::Reader***/
TTree* InputTree = new TTree("InTMP", "InTMP");
vector<float> InputTMVA;
for(int i = 0 ; i< (int)fMapOfTMVAVariableNames.size() ; i++)
InputTMVA.push_back(0);
float Time=0;
IsotopicVector IVInputTMVA;
map<ZAI ,string >::iterator it;
int j=0;
for( it = fMapOfTMVAVariableNames.begin() ; it!=fMapOfTMVAVariableNames.end() ; it++)
{
InputTree->Branch( ((*it).second).c_str() ,&InputTMVA[j], ((*it).second + "/F").c_str());
IVInputTMVA+= ((*it).first)*1;
j++;
}
if( !fIsStepTime)
InputTree->Branch( "Time" ,&Time ,"Time/F" );
IsotopicVector IVAccordingToUserInfoFile = isotopicvector.GetThisComposition(IVInputTMVA);
double Ntot = IVAccordingToUserInfoFile.GetTotalMass()*1e6/IVAccordingToUserInfoFile.MeanMolar()*6.02214129e23;
IVAccordingToUserInfoFile = IVAccordingToUserInfoFile/Ntot;
j=0;
map<ZAI ,string >::iterator it2;
for( it2 = fMapOfTMVAVariableNames.begin() ; it2!=fMapOfTMVAVariableNames.end() ; it2++)
{
InputTMVA[j] = IVAccordingToUserInfoFile.GetZAIIsotopicQuantity( (*it2).first ) ;
INFO<< (*it2).first.Z()<<" "<<(*it2).first.A()<<" "<<InputTMVA[j]<<endl;
j++;
}
Time=fMLP_Time[TimeStep];
InputTree->Fill();
return InputTree;
}
//________________________________________________________________________
double XSM_MLP::ExecuteTMVA(string WeightFile,TTree* InputTree)
{
// --- Create the Reader object
TMVA::Reader *reader = new TMVA::Reader( "Silent" );
// Create a set of variables and declare them to the reader
// - the variable names MUST corresponds in name and type to those given in the weight file(s) used
vector<float> InputTMVA;
for(int i = 0 ; i< (int)fMapOfTMVAVariableNames.size() ; i++)
InputTMVA.push_back(0);
Float_t Time;
map<ZAI ,string >::iterator it;
int j=0;
for( it = fMapOfTMVAVariableNames.begin() ; it!=fMapOfTMVAVariableNames.end() ; it++)
{ reader->AddVariable( ( (*it).second ).c_str(),&InputTMVA[j]);
j++;
}
if(!fIsStepTime)
reader->AddVariable( "Time" ,&Time);
// --- Book the MVA methods
string dir = fTMVAWeightFolder;
if(dir[dir.size()-1]!='/')
dir+="/";
// Book method MLP
TString methodName = "MLP method";
TString weightpath = dir + WeightFile ;
reader->BookMVA( methodName, weightpath );
map<ZAI ,string >::iterator it2;
j=0;
for( it2 = fMapOfTMVAVariableNames.begin() ; it2!=fMapOfTMVAVariableNames.end() ; it2++)
{
InputTree->SetBranchAddress(( (*it2).second ).c_str(),&InputTMVA[j]);
j++;
}
if(!fIsStepTime)
InputTree->SetBranchAddress( "Time" ,&Time );
InputTree->GetEntry(0);
Float_t val = (reader->EvaluateRegression( methodName ))[0];
delete reader;
return (double)val;
}
//________________________________________________________________________
EvolutionData XSM_MLP::GetCrossSectionsTime(IsotopicVector IV)
{DBGL
EvolutionData EvolutionDataFromMLP = EvolutionData();
map<ZAI,TGraph*> ExtrapolatedXS[3];
/*************DATA BASE INFO****************/
EvolutionDataFromMLP.SetReactorType(fDataBaseRType);
EvolutionDataFromMLP.SetFuelType(fDataBaseFType);
EvolutionDataFromMLP.SetPower(fDataBasePower);
EvolutionDataFromMLP.SetHeavyMetalMass(fDataBaseHMMass);
/************* The Cross sections***********/
for(int i=0;i<int(fWeightFiles.size());i++)
{
int Z=-2;
int A=-2;
int I=-2;
int Reaction=-2;
ReadWeightFile( fWeightFiles[i], Z, A, I, Reaction);
for(int TimeStep=0;TimeStep<int(fMLP_Time.size());TimeStep++)
{
TTree* InputTree = CreateTMVAInputTree(IV,TimeStep);
pair< map<ZAI, TGraph*>::iterator, bool> IResult;
IResult = ExtrapolatedXS[Reaction].insert(pair<ZAI ,TGraph* >(ZAI(Z,A,I), new TGraph() ) );
double XSValue = ExecuteTMVA(fWeightFiles[i],InputTree );
if(IResult.second )
{
(IResult.first)->second->SetPoint(0, (double)GetMLPTime()[TimeStep], XSValue );
}
else
{
(IResult.first)->second->SetPoint( (IResult.first)->second->GetN(), (double)GetMLPTime()[TimeStep], XSValue );
}
delete InputTree;
}
}
/**********Sorting TGraph*********/
for(int x=0;x<3;x++)
{ map<ZAI,TGraph*>::iterator it;
for(it = ExtrapolatedXS[x].begin(); it != ExtrapolatedXS[x].end(); it++)
it->second->Sort();
}
/**********Filling Matrices*/
EvolutionDataFromMLP.SetFissionXS(ExtrapolatedXS[0]);
EvolutionDataFromMLP.SetCaptureXS(ExtrapolatedXS[1]);
EvolutionDataFromMLP.Setn2nXS(ExtrapolatedXS[2]);
DBGL
return EvolutionDataFromMLP;
}
//________________________________________________________________________
//________________________________________________________________________
//
// Time step (1 MLP per time step)
//________________________________________________________________________
//________________________________________________________________________
void XSM_MLP::ReadWeightFileStep(string Filename, int &Z, int &A, int &I, int &Reaction, int &TimeStep)
{
Z=-1;
A=-1;
I=-1;
Reaction=-1;
size_t found = Filename.find("XS");
string NameJOB;
NameJOB=Filename.substr(found);
int pos=0;
StringLine::NextWord(NameJOB, pos, '_');
Z = atof( (StringLine::NextWord(NameJOB,pos,'_') ).c_str() );
A = atof( (StringLine::NextWord(NameJOB,pos,'_') ).c_str() );
I = atof( (StringLine::NextWord(NameJOB,pos,'_') ).c_str() );
string sReaction = (StringLine::NextWord(NameJOB,pos,'_') ).c_str() ;
if(sReaction=="fis")
Reaction=0;
if(sReaction=="cap")
Reaction=1;
if(sReaction=="n2n")
Reaction=2;
TimeStep = atof( (StringLine::NextWord(NameJOB,pos,'_') ).c_str() );
if(Z==-1 || A==-1 || I==-1 || Reaction==-1 || TimeStep==-1)
{
ERROR << " wrong TMVA weight format " << endl;
exit(0);
}
}
//________________________________________________________________________
EvolutionData XSM_MLP::GetCrossSectionsStep(IsotopicVector IV)
{DBGL
TTree* InputTree=CreateTMVAInputTree(IV);
//cout<<"=====Building Evolution Data From TMVA MLP====="<<endl;
EvolutionData EvolutionDataFromMLP = EvolutionData();
map<ZAI,TGraph*> ExtrapolatedXS[3];
/*************DATA BASE INFO****************/
EvolutionDataFromMLP.SetReactorType("PWR");
EvolutionDataFromMLP.SetFuelType("MOX");
EvolutionDataFromMLP.SetPower(fDataBasePower);
EvolutionDataFromMLP.SetHeavyMetalMass(fDataBaseHMMass);
/************* The Cross sections***********/
for(int i=0;i<int(fWeightFiles.size());i++)
{
int Z=-2;
int A=-2;
int I=-2;
int Reaction=-2;
int TimeStep=-2;
ReadWeightFileStep( fWeightFiles[i], Z, A, I, Reaction, TimeStep);
ZAI zaitmp = ZAI(Z,A,I);
pair< map<ZAI, TGraph*>::iterator, bool> IResult;
IResult = ExtrapolatedXS[Reaction].insert(pair<ZAI ,TGraph* >(ZAI(Z,A,I), new TGraph() ) );
if( IResult.second )
{
(IResult.first)->second->SetPoint(0, (double)GetMLPTime()[TimeStep], ExecuteTMVA(fWeightFiles[i],InputTree) );
}
else
{
(IResult.first)->second->SetPoint( (IResult.first)->second->GetN(), (double)GetMLPTime()[TimeStep], ExecuteTMVA(fWeightFiles[i],InputTree) );
}
}
/**********Sorting TGraph*********/
for(int x=0;x<3;x++)
{ map<ZAI,TGraph*>::iterator it;
for(it = ExtrapolatedXS[x].begin(); it != ExtrapolatedXS[x].end(); it++)
it->second->Sort();
}
/**********Filling Matrices*/
EvolutionDataFromMLP.SetFissionXS(ExtrapolatedXS[0]);
EvolutionDataFromMLP.SetCaptureXS(ExtrapolatedXS[1]);
EvolutionDataFromMLP.Setn2nXS(ExtrapolatedXS[2]);
delete InputTree;
DBGL
return EvolutionDataFromMLP;
}
//________________________________________________________________________
EvolutionData XSM_MLP::GetCrossSections(IsotopicVector IV ,double t)
{DBGL
if(t!=0)
WARNING << " Argument t has non effect here " << endl;
EvolutionData EV;
if(fIsStepTime)
EV=GetCrossSectionsStep(IV);
else
EV=GetCrossSectionsTime(IV);
DBGL
return EV;
}
//________________________________________________________________________
source/branches/CLASSV3/src/XSModel.cxx
+ 8
− 3
View file @ c126c0ed
......@@ -7,12 +7,17 @@
//
#include "XSModel.hxx"
#include "CLASSHeaders.hxx"
using namespace std;
XSModel::XSModel(): CLASSObject()
XSModel::XSModel(): CLASSObject(new CLASSLogger("XSModel.log"))
{
}
XSModel::XSModel(CLASSLogger* log): CLASSObject(log)
{
}
\ No newline at end of file
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment