diff --git a/source/Model/Equivalence/EQ_OneParameter.cxx b/source/Model/Equivalence/EQ_OneParameter.cxx
index 656f740188d915c02e6f167ccae7688702a768e2..2a314ed941e3385cb86930383286e094a6bde3d6 100644
--- a/source/Model/Equivalence/EQ_OneParameter.cxx
+++ b/source/Model/Equivalence/EQ_OneParameter.cxx
@@ -1,4 +1,3 @@
-#include "EquivalenceModel.hxx"
#include "EQ_OneParameter.hxx"
#include "StringLine.hxx"
#include "CLASSMethod.hxx"
@@ -19,7 +18,7 @@
//________________________________________________________________________
//________________________________________________________________________
-EQ_OneParameter::EQ_OneParameter(string TMVAXMLFilePath, string TMVANFOFilePath): EquivalenceModel(new CLASSLogger("EQ_OneParameter.log"))
+EQ_OneParameter::EQ_OneParameter(string TMVAXMLFilePath, string TMVANFOFilePath):EquivalenceModel(new CLASSLogger("EQ_OneParameter.log"))
{
fUseTMVAPredictor = true;
@@ -41,22 +40,22 @@ EQ_OneParameter::EQ_OneParameter(string TMVAXMLFilePath, string TMVANFOFilePath)
LoadKeyword(); // Load Key words defineds in NFO file
ReadNFO(); //Getting information from file NFO
-
+
//Check if any information is missing in NFO file
- if (fZAILimits.empty()) {ERROR << "Missing information for k_zail in : " << fTMVANFOFilePath << endl; exit(1);}
- if (fDBRType.empty()) {ERROR << "Missing information for k_reactor in : " << fTMVANFOFilePath << endl; exit(1);}
- if (fDBFType.empty()) {ERROR << "Missing information for k_fuel in : " << fTMVANFOFilePath << endl; exit(1);}
- if (!fSpecificPower) {ERROR << "Missing information for k_specpower in : " << fTMVANFOFilePath << endl; exit(1);}
- if (!fMaximalBU) {ERROR << "Missing information for k_maxburnup in : " << fTMVANFOFilePath << endl; exit(1);}
- if (fStreamListEqMMassFractionMin.empty() || fStreamListEqMMassFractionMax.empty()) { ERROR << "Missing information for k_massfractionmin and/or k_massfractionmax in : " << fTMVANFOFilePath << endl; exit(1);}
- if (fStreamList.empty()) { ERROR << "Missing information for k_list in : " << fTMVANFOFilePath << endl; exit(1); }
- if (fMapOfTMVAVariableNames.empty()) { ERROR << "Missing information for k_zainame in : " << fTMVANFOFilePath << endl; exit(1);}
- if (fTargetParameter.empty()) { ERROR << "Missing information for k_targetparameter in : " << fTMVANFOFilePath << endl; exit(1);}
- if (!fTargetParameterStDev) { ERROR << "Missing information for fTargetParameterStDev in : " << fTMVANFOFilePath << endl; exit(1);}
- if (fModelParameter.empty()) { ERROR << "Missing information for k_modelparameter in : " << fTMVANFOFilePath << endl; exit(1);}
- if (fBuffer.empty()) { ERROR << "Missing information for k_buffer in : " << fTMVANFOFilePath << endl; exit(1);}
- if (fPredictorType.empty()) { ERROR << "Missing information for k_predictortype in : " << fTMVANFOFilePath << endl; exit(1);}
- if (fOutput.empty()) { ERROR << "Missing information for k_output in : " << fTMVANFOFilePath << endl; exit(1);}
+ if(fZAILimits.empty()) {ERROR<<"Missing information for k_zail in : "<<fTMVANFOFilePath<<endl; exit(1);}
+ if(fDBRType.empty()) {ERROR<<"Missing information for k_reactor in : "<<fTMVANFOFilePath<<endl; exit(1);}
+ if(fDBFType.empty()) {ERROR<<"Missing information for k_fuel in : "<<fTMVANFOFilePath<<endl; exit(1);}
+ if(!fSpecificPower) {ERROR<<"Missing information for k_specpower in : "<<fTMVANFOFilePath<<endl; exit(1);}
+ if(!fMaximalBU) {ERROR<<"Missing information for k_maxburnup in : "<<fTMVANFOFilePath<<endl; exit(1);}
+ if(fStreamListEqMMassFractionMin.empty() || fStreamListEqMMassFractionMax.empty()) { ERROR<<"Missing information for k_massfractionmin and/or k_massfractionmax in : "<<fTMVANFOFilePath<<endl; exit(1);}
+ if(fStreamList.empty()) { ERROR<<"Missing information for k_list in : "<<fTMVANFOFilePath<<endl; exit(1); }
+ if(fMapOfTMVAVariableNames.empty()) { ERROR<<"Missing information for k_zainame in : "<<fTMVANFOFilePath<<endl; exit(1);}
+ if(fTargetParameter.empty()) { ERROR<<"Missing information for k_targetparameter in : "<<fTMVANFOFilePath<<endl; exit(1);}
+ if(!fTargetParameterStDev) { ERROR<<"Missing information for fTargetParameterStDev in : "<<fTMVANFOFilePath<<endl; exit(1);}
+ if(fModelParameter.empty()) { ERROR<<"Missing information for k_modelparameter in : "<<fTMVANFOFilePath<<endl; exit(1);}
+ if(fBuffer.empty()) { ERROR<<"Missing information for k_buffer in : "<<fTMVANFOFilePath<<endl; exit(1);}
+ if(fPredictorType.empty()) { ERROR<<"Missing information for k_predictortype in : "<<fTMVANFOFilePath<<endl; exit(1);}
+ if(fOutput.empty()) { ERROR<<"Missing information for k_output in : "<<fTMVANFOFilePath<<endl; exit(1);}
INFO << "__An equivalence model has been define__" << endl;
INFO << "\tThe TMVA weights file is :" << fTMVAXMLFilePath << endl;
@@ -65,12 +64,12 @@ EQ_OneParameter::EQ_OneParameter(string TMVAXMLFilePath, string TMVANFOFilePath)
}
//________________________________________________________________________
-EQ_OneParameter::EQ_OneParameter(CLASSLogger* log, string TMVAXMLFilePath, string TMVANFOFilePath): EquivalenceModel(log)
+EQ_OneParameter::EQ_OneParameter(CLASSLogger* log, string TMVAXMLFilePath, string TMVANFOFilePath):EquivalenceModel(log)
{
fUseTMVAPredictor = true;
fMaxIterration = 500;
- freaded = false;
+ freaded = false;
fPCMprecision = 10;
fTMVAXMLFilePath = TMVAXMLFilePath;
@@ -90,28 +89,27 @@ EQ_OneParameter::EQ_OneParameter(CLASSLogger* log, string TMVAXMLFilePath, strin
ReadNFO(); //Getting information from file NFO
//Check if any information is missing in NFO file
- if (fZAILimits.empty()) {ERROR << "Missing information for k_zail in : " << fTMVANFOFilePath << endl; exit(1);}
- if (fDBRType.empty()) {ERROR << "Missing information for k_reactor in : " << fTMVANFOFilePath << endl; exit(1);}
- if (fDBFType.empty()) {ERROR << "Missing information for k_fuel in : " << fTMVANFOFilePath << endl; exit(1);}
- if (!fSpecificPower) {ERROR << "Missing information for k_specpower in : " << fTMVANFOFilePath << endl; exit(1);}
- if (!fMaximalBU) {ERROR << "Missing information for k_maxburnup in : " << fTMVANFOFilePath << endl; exit(1);}
- if (fStreamListEqMMassFractionMin.empty() || fStreamListEqMMassFractionMax.empty()) { ERROR << "Missing information for k_massfractionmin and/or k_massfractionmax in : " << fTMVANFOFilePath << endl; exit(1);}
- if (fStreamList.empty()) { ERROR << "Missing information for k_list in : " << fTMVANFOFilePath << endl; exit(1); }
- if (fMapOfTMVAVariableNames.empty()) { ERROR << "Missing information for k_zainame in : " << fTMVANFOFilePath << endl; exit(1);}
- if (fTargetParameter.empty()) { ERROR << "Missing information for k_targetparameter in : " << fTMVANFOFilePath << endl; exit(1);}
- if (!fTargetParameterStDev) { ERROR << "Missing information for fTargetParameterStDev in : " << fTMVANFOFilePath << endl; exit(1);}
- if (fModelParameter.empty()) { ERROR << "Missing information for k_modelparameter in : " << fTMVANFOFilePath << endl; exit(1);}
- if (fBuffer.empty()) { ERROR << "Missing information for k_buffer in : " << fTMVANFOFilePath << endl; exit(1);}
- if (fPredictorType.empty()) { ERROR << "Missing information for k_predictortype in : " << fTMVANFOFilePath << endl; exit(1);}
- if (fOutput.empty()) { ERROR << "Missing information for k_output in : " << fTMVANFOFilePath << endl; exit(1);}
+ if(fZAILimits.empty()) {ERROR<<"Missing information for k_zail in : "<<fTMVANFOFilePath<<endl; exit(1);}
+ if(fDBRType.empty()) {ERROR<<"Missing information for k_reactor in : "<<fTMVANFOFilePath<<endl; exit(1);}
+ if(fDBFType.empty()) {ERROR<<"Missing information for k_fuel in : "<<fTMVANFOFilePath<<endl; exit(1);}
+ if(!fSpecificPower) {ERROR<<"Missing information for k_specpower in : "<<fTMVANFOFilePath<<endl; exit(1);}
+ if(!fMaximalBU) {ERROR<<"Missing information for k_maxburnup in : "<<fTMVANFOFilePath<<endl; exit(1);}
+ if(fStreamListEqMMassFractionMin.empty() || fStreamListEqMMassFractionMax.empty()) { ERROR<<"Missing information for k_massfractionmin and/or k_massfractionmax in : "<<fTMVANFOFilePath<<endl; exit(1);}
+ if(fStreamList.empty()) { ERROR<<"Missing information for k_list in : "<<fTMVANFOFilePath<<endl; exit(1); }
+ if(fMapOfTMVAVariableNames.empty()) { ERROR<<"Missing information for k_zainame in : "<<fTMVANFOFilePath<<endl; exit(1);}
+ if(fTargetParameter.empty()) { ERROR<<"Missing information for k_targetparameter in : "<<fTMVANFOFilePath<<endl; exit(1);}
+ if(!fTargetParameterStDev) { ERROR<<"Missing information for fTargetParameterStDev in : "<<fTMVANFOFilePath<<endl; exit(1);}
+ if(fModelParameter.empty()) { ERROR<<"Missing information for k_modelparameter in : "<<fTMVANFOFilePath<<endl; exit(1);}
+ if(fBuffer.empty()) { ERROR<<"Missing information for k_buffer in : "<<fTMVANFOFilePath<<endl; exit(1);}
+ if(fPredictorType.empty()) { ERROR<<"Missing information for k_predictortype in : "<<fTMVANFOFilePath<<endl; exit(1);}
+ if(fOutput.empty()) { ERROR<<"Missing information for k_output in : "<<fTMVANFOFilePath<<endl; exit(1);}
INFO << "__An equivalence model has been define__" << endl;
INFO << "\tThe TMVA weights file is :" << fTMVAXMLFilePath << endl;
INFO << "\tThe TMVA NFO file is :" << fTMVANFOFilePath << endl;
- PrintInfo();
-}
+ PrintInfo();}
//________________________________________________________________________
-EQ_OneParameter::EQ_OneParameter(string TMVANFOFilePath): EquivalenceModel(new CLASSLogger("EQ_OneParameter.log"))
+EQ_OneParameter::EQ_OneParameter(string TMVANFOFilePath):EquivalenceModel(new CLASSLogger("EQ_OneParameter.log"))
{
fUseTMVAPredictor = false;
@@ -127,23 +125,23 @@ EQ_OneParameter::EQ_OneParameter(string TMVANFOFilePath): EquivalenceModel(new C
LoadKeyword(); // Load Key words defineds in NFO file
ReadNFO(); //Getting information from file NFO
-
+
//Check if any information is missing in NFO file
- if (fZAILimits.empty()) {ERROR << "Missing information for k_zail in : " << fTMVANFOFilePath << endl; exit(1);}
- if (fDBRType.empty()) {ERROR << "Missing information for k_reactor in : " << fTMVANFOFilePath << endl; exit(1);}
- if (fDBFType.empty()) {ERROR << "Missing information for k_fuel in : " << fTMVANFOFilePath << endl; exit(1);}
- if (!fSpecificPower) {ERROR << "Missing information for k_specpower in : " << fTMVANFOFilePath << endl; exit(1);}
- if (!fMaximalBU) {ERROR << "Missing information for k_maxburnup in : " << fTMVANFOFilePath << endl; exit(1);}
- if (fStreamListEqMMassFractionMin.empty() || fStreamListEqMMassFractionMax.empty()) { ERROR << "Missing information for k_massfractionmin and/or k_massfractionmax in : " << fTMVANFOFilePath << endl; exit(1);}
- if (fStreamList.empty()) { ERROR << "Missing information for k_list in : " << fTMVANFOFilePath << endl; exit(1); }
- if (fBuffer.empty()) { ERROR << "Missing information for k_buffer in : " << fTMVANFOFilePath << endl; exit(1);}
+ if(fZAILimits.empty()) {ERROR<<"Missing information for k_zail in : "<<fTMVANFOFilePath<<endl; exit(1);}
+ if(fDBRType.empty()) {ERROR<<"Missing information for k_reactor in : "<<fTMVANFOFilePath<<endl; exit(1);}
+ if(fDBFType.empty()) {ERROR<<"Missing information for k_fuel in : "<<fTMVANFOFilePath<<endl; exit(1);}
+ if(!fSpecificPower) {ERROR<<"Missing information for k_specpower in : "<<fTMVANFOFilePath<<endl; exit(1);}
+ if(!fMaximalBU) {ERROR<<"Missing information for k_maxburnup in : "<<fTMVANFOFilePath<<endl; exit(1);}
+ if(fStreamListEqMMassFractionMin.empty() || fStreamListEqMMassFractionMax.empty()) { ERROR<<"Missing information for k_massfractionmin and/or k_massfractionmax in : "<<fTMVANFOFilePath<<endl; exit(1);}
+ if(fStreamList.empty()) { ERROR<<"Missing information for k_list in : "<<fTMVANFOFilePath<<endl; exit(1); }
+ if(fBuffer.empty()) { ERROR<<"Missing information for k_buffer in : "<<fTMVANFOFilePath<<endl; exit(1);}
INFO << "__An equivalence model without TMVA data has been define__" << endl;
INFO << "\tThe NFO file is :" << fTMVANFOFilePath << endl;
PrintInfo();
}
//________________________________________________________________________
-EQ_OneParameter::EQ_OneParameter(CLASSLogger* log, string TMVANFOFilePath): EquivalenceModel(log)
+EQ_OneParameter::EQ_OneParameter(CLASSLogger* log, string TMVANFOFilePath):EquivalenceModel(log)
{
fUseTMVAPredictor = false;
@@ -159,24 +157,23 @@ EQ_OneParameter::EQ_OneParameter(CLASSLogger* log, string TMVANFOFilePath): Equi
ReadNFO(); //Getting information from file NFO
//Check if any information is missing in NFO file
- if (fZAILimits.empty()) {ERROR << "Missing information for k_zail in : " << fTMVANFOFilePath << endl; exit(1);}
- if (fDBRType.empty()) {ERROR << "Missing information for k_reactor in : " << fTMVANFOFilePath << endl; exit(1);}
- if (fDBFType.empty()) {ERROR << "Missing information for k_fuel in : " << fTMVANFOFilePath << endl; exit(1);}
- if (!fSpecificPower) {ERROR << "Missing information for k_specpower in : " << fTMVANFOFilePath << endl; exit(1);}
- if (!fMaximalBU) {ERROR << "Missing information for k_maxburnup in : " << fTMVANFOFilePath << endl; exit(1);}
- if (fStreamListEqMMassFractionMin.empty() || fStreamListEqMMassFractionMax.empty()) { ERROR << "Missing information for k_massfractionmin and/or k_massfractionmax in : " << fTMVANFOFilePath << endl; exit(1);}
- if (fStreamList.empty()) { ERROR << "Missing information for k_list in : " << fTMVANFOFilePath << endl; exit(1); }
- if (fBuffer.empty()) { ERROR << "Missing information for k_buffer in : " << fTMVANFOFilePath << endl; exit(1);}
+ if(fZAILimits.empty()) {ERROR<<"Missing information for k_zail in : "<<fTMVANFOFilePath<<endl; exit(1);}
+ if(fDBRType.empty()) {ERROR<<"Missing information for k_reactor in : "<<fTMVANFOFilePath<<endl; exit(1);}
+ if(fDBFType.empty()) {ERROR<<"Missing information for k_fuel in : "<<fTMVANFOFilePath<<endl; exit(1);}
+ if(!fSpecificPower) {ERROR<<"Missing information for k_specpower in : "<<fTMVANFOFilePath<<endl; exit(1);}
+ if(!fMaximalBU) {ERROR<<"Missing information for k_maxburnup in : "<<fTMVANFOFilePath<<endl; exit(1);}
+ if(fStreamListEqMMassFractionMin.empty() || fStreamListEqMMassFractionMax.empty()) { ERROR<<"Missing information for k_massfractionmin and/or k_massfractionmax in : "<<fTMVANFOFilePath<<endl; exit(1);}
+ if(fStreamList.empty()) { ERROR<<"Missing information for k_list in : "<<fTMVANFOFilePath<<endl; exit(1); }
+ if(fBuffer.empty()) { ERROR<<"Missing information for k_buffer in : "<<fTMVANFOFilePath<<endl; exit(1);}
INFO << "__An equivalence model has been define__" << endl;
INFO << "\tThe TMVA weights file is :" << fTMVAXMLFilePath << endl;
INFO << "\tThe TMVA NFO file is :" << fTMVANFOFilePath << endl;
- PrintInfo();
-}
+ PrintInfo();}
//________________________________________________________________________
EQ_OneParameter::~EQ_OneParameter()
{
-
+
}
//________________________________________________________________________
IsotopicVector EQ_OneParameter::BuildFuelToTest(map < string, vector<double> >& lambda, map < string , vector <IsotopicVector> > const& StreamArray, double HMMass, map <string, bool> StreamListIsBuffer)
@@ -187,29 +184,29 @@ IsotopicVector EQ_OneParameter::BuildFuelToTest(map < string, vector<double> >&
map < string , bool >::iterator it_s_B;
//Find the buffer and set its lambda to 0
- string BufferDenomination = "";
- for ( it_s_B = StreamListIsBuffer.begin(); it_s_B != StreamListIsBuffer.end(); it_s_B++)
+ string BufferDenomination ="";
+ for( it_s_B = StreamListIsBuffer.begin(); it_s_B != StreamListIsBuffer.end(); it_s_B++)
{
- if ((*it_s_B ).second == true) { BufferDenomination = (*it_s_B).first; }
+ if((*it_s_B ).second==true){ BufferDenomination = (*it_s_B).first; }
}
- for (int i = 0; i < lambda[BufferDenomination].size(); i++)
+ for(int i = 0; i< lambda[BufferDenomination].size(); i++)
{
- lambda[BufferDenomination][i] = 0;
+ lambda[BufferDenomination][i]=0;
}
-
+
//Build an IV with all materials besides buffer to get the total mass of others materials
IsotopicVector IV;
- for ( it_s_vIV = StreamArray.begin(); it_s_vIV != StreamArray.end(); it_s_vIV++)
- {
- for (int i = 0; i < (int)StreamArray.at( it_s_vIV->first ).size(); i++)
+ for( it_s_vIV = StreamArray.begin(); it_s_vIV != StreamArray.end(); it_s_vIV++)
+ {
+ for(int i=0; i < (int)StreamArray.at( it_s_vIV->first ).size(); i++)
{
- IV += lambda[(*it_s_vIV).first][i] * StreamArray.at( it_s_vIV->first )[i];
+ IV += lambda[(*it_s_vIV).first][i] * StreamArray.at( it_s_vIV->first )[i];
}
}
//Calculate MassBuffer
- double MassBuffer = HMMass - IV.GetTotalMass() * 1e06;
+ double MassBuffer = HMMass - IV.GetTotalMass()*1e06;
//Set buffer lambda according to MassBuffer
ConvertMassToLambdaVector(BufferDenomination, lambda[BufferDenomination], MassBuffer, StreamArray.at(BufferDenomination));
@@ -217,15 +214,15 @@ IsotopicVector EQ_OneParameter::BuildFuelToTest(map < string, vector<double> >&
IV.Clear();
//Build fuel with all materials
- for ( it_s_vIV = StreamArray.begin(); it_s_vIV != StreamArray.end(); it_s_vIV++)
- {
- for (int i = 0; i < (int)StreamArray.at( it_s_vIV->first ).size(); i++)
+ for( it_s_vIV = StreamArray.begin(); it_s_vIV != StreamArray.end(); it_s_vIV++)
+ {
+ for(int i=0; i < (int)StreamArray.at( it_s_vIV->first ).size(); i++)
{
- IV += lambda[(*it_s_vIV).first][i] * StreamArray.at( it_s_vIV->first )[i];
+ IV += lambda[(*it_s_vIV).first][i] * StreamArray.at( it_s_vIV->first )[i];
}
}
DBGL
- return IV;
+ return IV;
}
@@ -245,36 +242,34 @@ map <string , vector<double> > EQ_OneParameter::BuildFuel(double BurnUp, double
map < int , string >::iterator it_i_s;
// Initialize lambda to 0 //
- for ( it_s_vIV = StreamArray.begin(); it_s_vIV != StreamArray.end(); it_s_vIV++)
- {
- for (int i = 0; i < (int)StreamArray[(*it_s_vIV).first].size(); i++)
+ for( it_s_vIV = StreamArray.begin(); it_s_vIV != StreamArray.end(); it_s_vIV++)
+ {
+ for(int i=0; i < (int)StreamArray[(*it_s_vIV).first].size(); i++)
{
lambda[(*it_s_vIV).first].push_back(0);
}
}
// Test if there is at least one stock available in each list, otherwise fuel is not built //
- bool BreakReturnLambda = true;
+ bool BreakReturnLambda = false;
for( it_s_vIV = StreamArray.begin(); it_s_vIV != StreamArray.end(); it_s_vIV++)
{
- if(StreamArray[(*it_s_vIV).first].size() != 0)
+ if(StreamArray[(*it_s_vIV).first].size() == 0)
{
- BreakReturnLambda &= false;
+ WARNING << " No stock available for stream : "<< (*it_s_vIV).first <<". Fuel not built." << endl;
+ SetLambdaToErrorCode(lambda[(*it_s_vIV).first]);
+ BreakReturnLambda = true;
}
}
- if(BreakReturnLambda) {
- WARNING << " No stock available for stream : "<< (*it_s_vIV).first <<". Fuel not built." << endl;
- SetLambdaToErrorCode(lambda[(*it_s_vIV).first]);
- return lambda;
- }
+ if(BreakReturnLambda) { return lambda;}
// Check if the targeted burn-up is lower than maximum burn-up of model //
- if (BurnUp > fMaximalBU)
+ if(BurnUp > fMaximalBU)
{
- ERROR << " Targeted burn-up is higher than maximum burn-up defined in NFO file..." << endl;
- ERROR << " Targeted burn-up : " << BurnUp << " GWd/t" << endl;
- ERROR << " Maximum burn-up : " << fMaximalBU << " GWd/t" << endl;
- exit(1);
+ ERROR << " Targeted burn-up is higher than maximum burn-up defined in NFO file..."<< endl;
+ ERROR << " Targeted burn-up : "<<BurnUp<<" GWd/t"<<endl;
+ ERROR << " Maximum burn-up : "<<fMaximalBU<<" GWd/t"<<endl;
+ exit(1);
}
// Fissile fraction calculation is needed.
@@ -283,242 +278,242 @@ map <string , vector<double> > EQ_OneParameter::BuildFuel(double BurnUp, double
// Check if EQ_OneParameter->SetTMVAXMLFilePath() and/or EQ_OneParameter->SetTMVANFOFilePath() have been defined
if (fTMVAXMLFilePath.empty() || fTMVANFOFilePath.empty())
{
- ERROR << " TMVA XML and/or NFO File path are not defined..." << endl;
- ERROR << " You have to use EQ_OneParameter->SetTMVAXMLFilePath() and/or EQ_OneParameter->SetTMVANFOFilePath() methods." << endl;
+ ERROR << " TMVA XML and/or NFO File path are not defined..."<< endl;
+ ERROR << " You have to use EQ_OneParameter->SetTMVAXMLFilePath() and/or EQ_OneParameter->SetTMVANFOFilePath() methods."<<endl;
exit(1);
}
-
+
double TargetParameterValue = 0;
-
- for ( it_s_D = fModelParameter.begin(); it_s_D != fModelParameter.end(); it_s_D++)
- {
- if (fModelParameter[(*it_s_D).first] == -1)
+
+ for( it_s_D = fModelParameter.begin(); it_s_D != fModelParameter.end(); it_s_D++)
+ {
+ if(fModelParameter[(*it_s_D).first] == -1)
{
- ERROR << "Model parameter ( " << fModelParameter[(*it_s_D).first] << " ) value is not defined in the input." << endl;
- ERROR << "Use EqM->SetModelParameter( \" " << (*it_s_D).first << " \", value) to define it." << endl;
- exit(1);
+ ERROR<< "Model parameter ( "<<fModelParameter[(*it_s_D).first] << " ) value is not defined in the input." <<endl;
+ ERROR<< "Use EqM->SetModelParameter( \" "<<(*it_s_D).first<<" \", value) to define it." <<endl;
+ exit(1);
}
}
-
- if (fTargetParameter == "BurnUpMax") {TargetParameterValue = BurnUp;}
- else if (fTargetParameter == "keffBOC") {TargetParameterValue = fModelParameter["keffBOC"];}
+
+ if(fTargetParameter=="BurnUpMax") {TargetParameterValue = BurnUp;}
+ else if (fTargetParameter=="keffBOC") {TargetParameterValue = fModelParameter["keffBOC"];}
else
{
- ERROR << "Target parameter defined in InformationFile ( " << fTargetParameter << " ) doesn't exist." << endl;
- ERROR << "Possible target parameters for the moment are : " << endl;
- ERROR << " - BurnUpMax - Used for PWR" << endl;
- ERROR << " - keffBOC - Used for SFR" << endl;
- exit(1);
+ ERROR<< "Target parameter defined in InformationFile ( "<<fTargetParameter<<" ) doesn't exist." <<endl;
+ ERROR<< "Possible target parameters for the moment are : "<< endl;
+ ERROR<< " - BurnUpMax - Used for PWR" <<endl;
+ ERROR<< " - keffBOC - Used for SFR" <<endl;
+ exit(1);
}
-
+
/// Search for the minimum and maximum fraction of each material in fuel ///
- map < string, double > StreamListMassFractionMin ;
- map < string, double > StreamListMassFractionMax ;
- for ( it_s_D = StreamListFPMassFractionMin.begin(); it_s_D != StreamListFPMassFractionMin.end(); it_s_D++)
- {
- if (StreamListFPMassFractionMin[(*it_s_D).first] < fStreamListEqMMassFractionMin[(*it_s_D).first]) // if limits FP are lower than limits EqM
+ map < string, double > StreamListMassFractionMin ;
+ map < string, double > StreamListMassFractionMax ;
+ for( it_s_D = StreamListFPMassFractionMin.begin(); it_s_D != StreamListFPMassFractionMin.end(); it_s_D++)
+ {
+ if(StreamListFPMassFractionMin[(*it_s_D).first] < fStreamListEqMMassFractionMin[(*it_s_D).first]) // if limits FP are lower than limits EqM
{
- ERROR << " User mass fraction min requirement is lower than the model mass fraction min for list : " << (*it_s_D).first << endl;
- ERROR << " User mass fraction min requirement : " << StreamListFPMassFractionMin[(*it_s_D).first] << endl;
- ERROR << " Model mass fraction min requirement : " << fStreamListEqMMassFractionMin[(*it_s_D).first] << endl;
+ ERROR << " User mass fraction min requirement is lower than the model mass fraction min for list : "<<(*it_s_D).first << endl;
+ ERROR << " User mass fraction min requirement : "<<StreamListFPMassFractionMin[(*it_s_D).first]<<endl;
+ ERROR << " Model mass fraction min requirement : "<<fStreamListEqMMassFractionMin[(*it_s_D).first]<<endl;
exit(1);
}
else
{
StreamListMassFractionMin[(*it_s_D).first] = StreamListFPMassFractionMin[(*it_s_D).first];
}
- }
-
- for ( it_s_D = StreamListFPMassFractionMax.begin(); it_s_D != StreamListFPMassFractionMax.end(); it_s_D++)
- {
- if (StreamListFPMassFractionMax[(*it_s_D).first] > fStreamListEqMMassFractionMax[(*it_s_D).first]) // if limits FP are higher than limits EqM
+ }
+
+ for( it_s_D = StreamListFPMassFractionMax.begin(); it_s_D != StreamListFPMassFractionMax.end(); it_s_D++)
+ {
+ if(StreamListFPMassFractionMax[(*it_s_D).first] > fStreamListEqMMassFractionMax[(*it_s_D).first]) // if limits FP are higher than limits EqM
{
- ERROR << " User mass fraction max requirement is higher than the model mass fraction max for list : " << (*it_s_D).first << endl;
- ERROR << " User mass fraction max requirement : " << StreamListFPMassFractionMax[(*it_s_D).first] << endl;
- ERROR << " Model mass fraction max requirement : " << fStreamListEqMMassFractionMax[(*it_s_D).first] << endl;
+ ERROR << " User mass fraction max requirement is higher than the model mass fraction max for list : "<<(*it_s_D).first << endl;
+ ERROR << " User mass fraction max requirement : "<<StreamListFPMassFractionMax[(*it_s_D).first]<<endl;
+ ERROR << " Model mass fraction max requirement : "<<fStreamListEqMMassFractionMax[(*it_s_D).first]<<endl;
exit(1);
}
else
{
StreamListMassFractionMax[(*it_s_D).first] = StreamListFPMassFractionMax[(*it_s_D).first];
}
-
+
}
-
+
//Calculate Total mass in stock for each stream and fill fTotalMassInStocks
StocksTotalMassCalculation(StreamArray);
-
+
// Check if there is enough material in stock to satisfy mass fraction min //
- BreakReturnLambda = false;
- for ( it_s_D = StreamListMassFractionMin.begin(); it_s_D != StreamListMassFractionMin.end(); it_s_D++)
+ BreakReturnLambda = false;
+ for( it_s_D = StreamListMassFractionMin.begin(); it_s_D != StreamListMassFractionMin.end(); it_s_D++)
{
- if (fTotalMassInStocks[(*it_s_D).first] < HMMass * StreamListMassFractionMin[(*it_s_D).first])
+ if(fTotalMassInStocks[(*it_s_D).first]< HMMass*StreamListMassFractionMin[(*it_s_D).first])
{
- WARNING << " Not enough material : " << (*it_s_D).first << " in stocks to reach the build fuel lower limit of " << StreamListMassFractionMin[(*it_s_D).first] << " reactor mass. Fuel not built." << endl;
+ WARNING << " Not enough material : "<< (*it_s_D).first << " in stocks to reach the build fuel lower limit of "<<StreamListMassFractionMin[(*it_s_D).first]<<" reactor mass. Fuel not built." << endl;
SetLambdaToErrorCode(lambda[(*it_s_D).first]);
- BreakReturnLambda = true;
+ BreakReturnLambda = true;
}
}
- if (BreakReturnLambda) { return lambda;}
-
+ if(BreakReturnLambda) { return lambda;}
+
// Check if there is enough material in stock to satisfy mass fraction max, if not mass fraction max is set to MassINStock/MassReactor//
- for ( it_s_D = StreamListMassFractionMax.begin(); it_s_D != StreamListMassFractionMax.end(); it_s_D++)
+ for( it_s_D = StreamListMassFractionMax.begin(); it_s_D != StreamListMassFractionMax.end(); it_s_D++)
{
- if (fTotalMassInStocks[(*it_s_D).first] < HMMass * StreamListMassFractionMax[(*it_s_D).first])
- {
- StreamListMassFractionMax[(*it_s_D).first] = fTotalMassInStocks[(*it_s_D).first] / HMMass;
- WARNING << " Not enough material : " << (*it_s_D).first << " in stocks to reach the build fuel higher limit of " << StreamListMassFractionMax[(*it_s_D).first] << " reactor mass. " << endl;
- WARNING << " Mass fraction max of material : " << (*it_s_D).first << " is set to MassInStock/HMMassReactor : " << StreamListMassFractionMax[(*it_s_D).first] << endl;
+ if(fTotalMassInStocks[(*it_s_D).first]< HMMass*StreamListMassFractionMax[(*it_s_D).first])
+ {
+ StreamListMassFractionMax[(*it_s_D).first] = fTotalMassInStocks[(*it_s_D).first]/HMMass;
+ WARNING << " Not enough material : "<< (*it_s_D).first << " in stocks to reach the build fuel higher limit of "<<StreamListMassFractionMax[(*it_s_D).first]<<" reactor mass. " << endl;
+ WARNING << " Mass fraction max of material : "<< (*it_s_D).first << " is set to MassInStock/HMMassReactor : "<< StreamListMassFractionMax[(*it_s_D).first]<< endl;
}
}
-
+
//Check if TargetParameter is inside [TargetParameterMin, TargetParameterMax] associated to fraction Min et Max//
-
- map < string , double > MassMin;
- map < string , double > MassMax;
-
- map < string , double > TargetParameterMin;
- map < string , double > TargetParameterMax;
-
+
+ map < string , double > MassMin;
+ map < string , double > MassMax;
+
+ map < string , double > TargetParameterMin;
+ map < string , double > TargetParameterMax;
+
IsotopicVector FuelToTest;
-
+
bool TargetParameterIncluded = false;
- for ( it_i_s = StreamListPriority.begin(); it_i_s != StreamListPriority.end(); it_i_s++)
- {
+ for( it_i_s = StreamListPriority.begin(); it_i_s != StreamListPriority.end(); it_i_s++)
+ {
//Calculate TargetParameterMin for each possibility : min1 ; max1 + min2 ; max1 + max2 + min3 ....
MassMin[(*it_i_s ).second] = HMMass * StreamListMassFractionMin[(*it_i_s).second];
ConvertMassToLambdaVector((*it_i_s ).second, lambda[(*it_i_s ).second], MassMin[(*it_i_s ).second], StreamArray[(*it_i_s ).second]);
FuelToTest = BuildFuelToTest(lambda, StreamArray, HMMass, StreamListIsBuffer);
- FuelToTest = FuelToTest / FuelToTest.GetSumOfAll();
+ FuelToTest = FuelToTest/FuelToTest.GetSumOfAll();
TargetParameterMin[(*it_i_s ).second] = CalculateTargetParameter(FuelToTest, fTargetParameter);
-
+
//Check is TargetParameterMin < TargetParameter
- if (TargetParameterMin[(*it_i_s ).second] > TargetParameterValue)
+ if(TargetParameterMin[(*it_i_s ).second]>TargetParameterValue)
{
- if ((*it_i_s).first == 1) //Minimum of first material is too high
- {
- WARNING << "CRITICAL ! Minimum parameter value associated to the first priority material ( " << (*it_i_s ).second << " ) is higher than targeted parameter." << endl;
- WARNING << "Targeted parameter : " << fTargetParameter << " = " << TargetParameterValue << endl;
- WARNING << "Minimum parameter value : " << TargetParameterMin[(*it_i_s ).second] << endl;
- WARNING << "Try to increase targeted parameter." << endl;
- SetLambdaToErrorCode(lambda[(*it_i_s).second]);
- return lambda;
- DBGL
+ if((*it_i_s).first ==1) //Minimum of first material is too high
+ {
+ WARNING << "CRITICAL ! Minimum parameter value associated to the first priority material ( "<<(*it_i_s ).second <<" ) is higher than targeted parameter."<< endl;
+ WARNING << "Targeted parameter : "<<fTargetParameter<<" = "<<TargetParameterValue<<endl;
+ WARNING << "Minimum parameter value : " <<TargetParameterMin[(*it_i_s ).second]<<endl;
+ WARNING << "Try to increase targeted parameter." <<endl;
+ SetLambdaToErrorCode(lambda[(*it_i_s).second]);
+ return lambda;
+ DBGL
}
- else if ((*it_i_s).first > 1) //TargetParameter is located between max n-1 and min n
+ else if ((*it_i_s).first >1) //TargetParameter is located between max n-1 and min n
{
- WARNING << "CRITICAL ! Targeted parameter value ( " << fTargetParameter << " ) is located between 2 materials. " << endl;
+ WARNING << "CRITICAL ! Targeted parameter value ( "<<fTargetParameter<<" ) is located between 2 materials. "<<endl;
it_i_s --;
- WARNING << fTargetParameter << " of max fraction of material : " << (*it_i_s).second << " ---> " << TargetParameterMax[(*it_i_s ).second] << endl;
+ WARNING << fTargetParameter <<" of max fraction of material : "<< (*it_i_s).second<<" ---> "<<TargetParameterMax[(*it_i_s ).second]<<endl;
it_i_s ++;
- WARNING << fTargetParameter << " of min fraction of material : " << (*it_i_s ).second << " ---> " << TargetParameterMin[(*it_i_s ).second] << endl;
- WARNING << "Targeted " << fTargetParameter << " : " << TargetParameterValue << endl;
- WARNING << "Try to decrease mimimum fraction of : " << (*it_i_s ).second << endl;
- SetLambdaToErrorCode(lambda[(*it_i_s).second]);
- return lambda;
+ WARNING << fTargetParameter<< " of min fraction of material : "<< (*it_i_s ).second<<" ---> "<<TargetParameterMin[(*it_i_s ).second]<<endl;
+ WARNING << "Targeted "<<fTargetParameter<<" : " <<TargetParameterValue<<endl;
+ WARNING << "Try to decrease mimimum fraction of : "<< (*it_i_s ).second<<endl;
+ SetLambdaToErrorCode(lambda[(*it_i_s).second]);
+ return lambda;
}
}
FuelToTest.Clear();
-
+
//Calculate TargetParameter max for each possibility : max1 ; max1 + max2 ; max1 + max2 + max3 ....
- MassMax[(*it_i_s ).second] = HMMass * StreamListMassFractionMax[(*it_i_s).second];
+ MassMax[(*it_i_s ).second] = HMMass * StreamListMassFractionMax[(*it_i_s).second];
ConvertMassToLambdaVector((*it_i_s ).second, lambda[(*it_i_s ).second], MassMax[(*it_i_s ).second], StreamArray[(*it_i_s ).second]);
FuelToTest = BuildFuelToTest(lambda, StreamArray, HMMass, StreamListIsBuffer);
- FuelToTest = FuelToTest / FuelToTest.GetSumOfAll();
+ FuelToTest = FuelToTest/FuelToTest.GetSumOfAll();
TargetParameterMax[(*it_i_s ).second] = CalculateTargetParameter(FuelToTest, fTargetParameter);
-
- if (TargetParameterMax[(*it_i_s ).second] >= TargetParameterValue)
+
+ if(TargetParameterMax[(*it_i_s ).second]>=TargetParameterValue)
{
- TargetParameterIncluded = true ;
+ TargetParameterIncluded = true ;
break;
}
}
-
+
//Check if target parameter increases monotously with the material mass
CheckTargetParameterConsistency(StreamListPriority, TargetParameterMin, TargetParameterMax);
-
- if (!TargetParameterIncluded)
+
+ if(!TargetParameterIncluded)
{
- WARNING << "CRITICAL ! Maximum reachable " << fTargetParameter << " is lower than targeted " << fTargetParameter << ". " << endl;
- WARNING << "Targeted " << fTargetParameter << " = " << TargetParameterValue << endl;
- WARNING << "Maximum reachable " << fTargetParameter << " : " << TargetParameterMax[(*--StreamListPriority.end()).second] << endl;
- WARNING << "Try to increase maximum fraction of materials, or decrease " << fTargetParameter << " ." << endl;
- SetLambdaToErrorCode(lambda[(*--StreamListPriority.end()).second]);
- return lambda;
+ WARNING << "CRITICAL ! Maximum reachable "<<fTargetParameter<<" is lower than targeted "<< fTargetParameter<<". "<< endl;
+ WARNING << "Targeted "<<fTargetParameter<<" = "<<TargetParameterValue<<endl;
+ WARNING << "Maximum reachable "<<fTargetParameter<<" : "<<TargetParameterMax[(*--StreamListPriority.end()).second]<<endl;
+ WARNING << "Try to increase maximum fraction of materials, or decrease "<< fTargetParameter<<" ." <<endl;
+ SetLambdaToErrorCode(lambda[(*--StreamListPriority.end()).second]);
+ return lambda;
}
-
- //Search the TargetParameterValue location in the mass damain //
+
+ //Search the TargetParameterValue location in the mass damain //
string MaterialToSearch = (*it_i_s ).second;
double CalculatedTargetParameter = TargetParameterMax[MaterialToSearch] ; //Algo start with maximum point
double MassToAdd = MassMax[MaterialToSearch]; //Algo start with maximum point
-
- double LastMassMinus = MassMin[MaterialToSearch]; //Used in bissection method
- double LastMassPlus = MassMax[MaterialToSearch]; //Used in bissection method
-
+
+ double LastMassMinus = MassMin[MaterialToSearch]; //Used in bissection method
+ double LastMassPlus = MassMax[MaterialToSearch]; //Used in bissection method
+
int count = 0;
-
+
FuelToTest.Clear();
-
- /*
- if (fDBFType == "MOX")
- {
- cout<<"------------------------------------------------------"<<endl;
- cout<<"START ALGO -> BU, Mass "<<BurnUp<<" "<<HMMass<<endl;
- cout<<"------------------------------------------------------"<<endl;
- double MassTest = MassMin[MaterialToSearch];
- cout<<MaterialToSearch<<" "<<MassMax[MaterialToSearch]<<" "<<MassMin[MaterialToSearch]<<" "<<endl;
- do
- {
- ConvertMassToLambdaVector(MaterialToSearch, lambda[MaterialToSearch], MassTest, StreamArray[MaterialToSearch]);
- FuelToTest = BuildFuelToTest(lambda, StreamArray, HMMass, StreamListIsBuffer);
- FuelToTest = FuelToTest/FuelToTest.GetSumOfAll();
- CalculatedTargetParameter = CalculateTargetParameter(FuelToTest, fTargetParameter);
-
- cout<<"Lambda vector : "<<MaterialToSearch<<" - "; for(int i=0; i < (int)lambda[MaterialToSearch].size(); i++) cout<<lambda[MaterialToSearch][i]<<" ";
- cout<<endl;
-
-
- MassTest += (MassMax[MaterialToSearch] - MassMin[MaterialToSearch])/100.;
-
- cout<<MassTest<<" "<<CalculatedTargetParameter<<endl;
-
- } while (MassTest <= MassMax[MaterialToSearch]);
- cout<<"------------------------------------------------------"<<endl;
- cout<<"STOP ALGO EXIT(1)..."<<endl; exit(1);
- cout<<"------------------------------------------------------"<<endl;
- }
- */
-
+
+ /*
+ if (fDBFType == "MOX")
+ {
+ cout<<"------------------------------------------------------"<<endl;
+ cout<<"START ALGO -> BU, Mass "<<BurnUp<<" "<<HMMass<<endl;
+ cout<<"------------------------------------------------------"<<endl;
+ double MassTest = MassMin[MaterialToSearch];
+ cout<<MaterialToSearch<<" "<<MassMax[MaterialToSearch]<<" "<<MassMin[MaterialToSearch]<<" "<<endl;
+ do
+ {
+ ConvertMassToLambdaVector(MaterialToSearch, lambda[MaterialToSearch], MassTest, StreamArray[MaterialToSearch]);
+ FuelToTest = BuildFuelToTest(lambda, StreamArray, HMMass, StreamListIsBuffer);
+ FuelToTest = FuelToTest/FuelToTest.GetSumOfAll();
+ CalculatedTargetParameter = CalculateTargetParameter(FuelToTest, fTargetParameter);
+
+ cout<<"Lambda vector : "<<MaterialToSearch<<" - "; for(int i=0; i < (int)lambda[MaterialToSearch].size(); i++) cout<<lambda[MaterialToSearch][i]<<" ";
+ cout<<endl;
+
+
+ MassTest += (MassMax[MaterialToSearch] - MassMin[MaterialToSearch])/100.;
+
+ cout<<MassTest<<" "<<CalculatedTargetParameter<<endl;
+
+ } while (MassTest <= MassMax[MaterialToSearch]);
+ cout<<"------------------------------------------------------"<<endl;
+ cout<<"STOP ALGO EXIT(1)..."<<endl; exit(1);
+ cout<<"------------------------------------------------------"<<endl;
+ }
+ */
+
do
{
- if (count > fMaxIterration)
+ if(count > fMaxIterration)
{
ERROR << "CRITICAL ! Can't manage to predict fissile content\nHint : Try to decrease the precision on the target parameter using :\nYourEQ_OneParameter->SetTargetParameterStDev(Precision); " << endl;
- ERROR << "Targeted " << fTargetParameter << " : " << TargetParameterValue << endl;
- ERROR << "Last calculated " << fTargetParameter << " : " << CalculatedTargetParameter << endl;
- ERROR << "Last Fresh fuel normalized composition : " << endl;
- ERROR << FuelToTest.sPrint() << endl;
+ ERROR << "Targeted "<<fTargetParameter<<" : "<<TargetParameterValue<<endl;
+ ERROR << "Last calculated "<<fTargetParameter<<" : "<<CalculatedTargetParameter<<endl;
+ ERROR << "Last Fresh fuel normalized composition : " <<endl;
+ ERROR << FuelToTest.sPrint()<<endl;
exit(1);
}
-
- if ( (CalculatedTargetParameter - TargetParameterValue) < 0 ) //Need to add more fissile material in fuel
+
+ if( (CalculatedTargetParameter - TargetParameterValue) < 0 ) //Need to add more fissile material in fuel
{
LastMassMinus = MassToAdd;
- MassToAdd = MassToAdd + fabs(LastMassPlus - MassToAdd) / 2.;
+ MassToAdd = MassToAdd + fabs(LastMassPlus - MassToAdd)/2.;
}
- else if ( (CalculatedTargetParameter - TargetParameterValue) > 0) //Need to add less fissile material in fuel
+ else if( (CalculatedTargetParameter - TargetParameterValue) > 0) //Need to add less fissile material in fuel
{
LastMassPlus = MassToAdd;
- MassToAdd = MassToAdd - fabs(LastMassMinus - MassToAdd) / 2.;
+ MassToAdd = MassToAdd - fabs(LastMassMinus - MassToAdd)/2.;
}
- ConvertMassToLambdaVector(MaterialToSearch, lambda[MaterialToSearch], MassToAdd, StreamArray[MaterialToSearch]);
+ ConvertMassToLambdaVector(MaterialToSearch, lambda[MaterialToSearch], MassToAdd, StreamArray[MaterialToSearch]);
FuelToTest = BuildFuelToTest(lambda, StreamArray, HMMass, StreamListIsBuffer);
- FuelToTest = FuelToTest / FuelToTest.GetSumOfAll();
+ FuelToTest = FuelToTest/FuelToTest.GetSumOfAll();
CalculatedTargetParameter = CalculateTargetParameter(FuelToTest, fTargetParameter);
-
+
count ++;
-
- } while (fabs(TargetParameterValue - CalculatedTargetParameter) > GetTargetParameterStDev()*TargetParameterValue);
+
+ }while(fabs(TargetParameterValue - CalculatedTargetParameter) > GetTargetParameterStDev()*TargetParameterValue);
}
// Fissile fraction is imposed by the FP
// No need to use algo
@@ -528,84 +523,84 @@ map <string , vector<double> > EQ_OneParameter::BuildFuel(double BurnUp, double
// Check if EQ_OneParameter->SetTMVANFOFilePath() have been defined
if (fTMVANFOFilePath.empty())
{
- ERROR << " TMVA NFO File path is not defined..." << endl;
- ERROR << " You have to use EQ_OneParameter->SetTMVANFOFilePath() methods." << endl;
+ ERROR << " TMVA NFO File path is not defined..."<< endl;
+ ERROR << " You have to use EQ_OneParameter->SetTMVANFOFilePath() methods."<<endl;
exit(1);
}
/// Search for the fraction of each material in fuel ///
- map < string, double > StreamListMassFraction;
- for ( it_s_D = StreamListFPMassFractionMin.begin(); it_s_D != StreamListFPMassFractionMin.end(); it_s_D++)
- {
- if (StreamListFPMassFractionMin[(*it_s_D).first] < fStreamListEqMMassFractionMin[(*it_s_D).first]) // if limits FP are lower than limits EqM
+ map < string, double > StreamListMassFraction;
+ for( it_s_D = StreamListFPMassFractionMin.begin(); it_s_D != StreamListFPMassFractionMin.end(); it_s_D++)
+ {
+ if(StreamListFPMassFractionMin[(*it_s_D).first] < fStreamListEqMMassFractionMin[(*it_s_D).first]) // if limits FP are lower than limits EqM
{
- ERROR << " User mass fraction requirement is lower than the model mass fraction min for list : " << (*it_s_D).first << endl;
- ERROR << " User mass fraction requirement : " << StreamListFPMassFractionMin[(*it_s_D).first] << endl;
- ERROR << " Model mass fraction min requirement : " << fStreamListEqMMassFractionMin[(*it_s_D).first] << endl;
+ ERROR << " User mass fraction requirement is lower than the model mass fraction min for list : "<<(*it_s_D).first << endl;
+ ERROR << " User mass fraction requirement : "<<StreamListFPMassFractionMin[(*it_s_D).first]<<endl;
+ ERROR << " Model mass fraction min requirement : "<<fStreamListEqMMassFractionMin[(*it_s_D).first]<<endl;
exit(1);
}
- else if (StreamListFPMassFractionMax[(*it_s_D).first] > fStreamListEqMMassFractionMax[(*it_s_D).first]) // if limits FP are higher than limits EqM
+ else if(StreamListFPMassFractionMax[(*it_s_D).first] > fStreamListEqMMassFractionMax[(*it_s_D).first]) // if limits FP are higher than limits EqM
{
- ERROR << " User mass fraction requirement is higher than the model mass fraction max for list : " << (*it_s_D).first << endl;
- ERROR << " User mass fraction requirement : " << StreamListFPMassFractionMax[(*it_s_D).first] << endl;
- ERROR << " Model mass fraction max requirement : " << fStreamListEqMMassFractionMax[(*it_s_D).first] << endl;
+ ERROR << " User mass fraction requirement is higher than the model mass fraction max for list : "<<(*it_s_D).first << endl;
+ ERROR << " User mass fraction requirement : "<<StreamListFPMassFractionMax[(*it_s_D).first]<<endl;
+ ERROR << " Model mass fraction max requirement : "<<fStreamListEqMMassFractionMax[(*it_s_D).first]<<endl;
exit(1);
}
else
{
StreamListMassFraction[(*it_s_D).first] = StreamListFPMassFractionMin[(*it_s_D).first]; // Because here, min = max
}
- }
+ }
//Calculate Total mass in stock for each stream and fill fTotalMassInStocks
StocksTotalMassCalculation(StreamArray);
// Check if there is enough material in stock to satisfy requested mass fraction //
- BreakReturnLambda = false;
- for ( it_s_D = StreamListMassFraction.begin(); it_s_D != StreamListMassFraction.end(); it_s_D++)
+ BreakReturnLambda = false;
+ for( it_s_D = StreamListMassFraction.begin(); it_s_D != StreamListMassFraction.end(); it_s_D++)
{
- if (fTotalMassInStocks[(*it_s_D).first] < HMMass * StreamListMassFraction[(*it_s_D).first])
+ if(fTotalMassInStocks[(*it_s_D).first]< HMMass*StreamListMassFraction[(*it_s_D).first])
{
- WARNING << " Not enough material : " << (*it_s_D).first << " in stocks to reach the build fuel limit of " << StreamListMassFraction[(*it_s_D).first] << " reactor mass. Fuel not built." << endl;
+ WARNING << " Not enough material : "<< (*it_s_D).first << " in stocks to reach the build fuel limit of "<<StreamListMassFraction[(*it_s_D).first]<<" reactor mass. Fuel not built." << endl;
SetLambdaToErrorCode(lambda[(*it_s_D).first]);
- BreakReturnLambda = true;
+ BreakReturnLambda = true;
}
}
- if (BreakReturnLambda) { return lambda;}
+ if(BreakReturnLambda) { return lambda;}
IsotopicVector FuelToTest;
// Build Fuel
- for ( it_i_s = StreamListPriority.begin(); it_i_s != StreamListPriority.end(); it_i_s++)
- {
+ for( it_i_s = StreamListPriority.begin(); it_i_s != StreamListPriority.end(); it_i_s++)
+ {
FuelToTest.Clear();
- ConvertMassToLambdaVector((*it_i_s).second, lambda[(*it_i_s).second], HMMass * StreamListMassFraction[(*it_i_s).second], StreamArray[(*it_i_s).second]);
+ ConvertMassToLambdaVector((*it_i_s).second, lambda[(*it_i_s).second], HMMass*StreamListMassFraction[(*it_i_s).second], StreamArray[(*it_i_s).second]);
FuelToTest = BuildFuelToTest(lambda, StreamArray, HMMass, StreamListIsBuffer);
}
}
- //Final builded fuel
+ //Final builded fuel
IsotopicVector IVStream;
- for ( it_s_vD = lambda.begin(); it_s_vD != lambda.end(); it_s_vD++)
+ for( it_s_vD = lambda.begin(); it_s_vD != lambda.end(); it_s_vD++)
{
- for (int i = 0; i < (int)lambda[(*it_s_vD).first].size(); i++)
+ for(int i=0; i<(int)lambda[(*it_s_vD).first].size(); i++)
{
- IVStream += lambda[(*it_s_vD).first][i] * StreamArray[(*it_s_vD).first][i];
+ IVStream +=lambda[(*it_s_vD).first][i] * StreamArray[(*it_s_vD).first][i];
}
- }
-
- //Check if BuildedFuel is in Model isotopic bounds
+ }
+
+ //Check if BuildedFuel is in Model isotopic bounds
(*this).isIVInDomain(IVStream);
+
+ for( it_s_vD = lambda.begin(); it_s_vD != lambda.end(); it_s_vD++)
+ {
+ DBGV( "Lambda vector : "<<(*it_s_vD).first );
- for ( it_s_vD = lambda.begin(); it_s_vD != lambda.end(); it_s_vD++)
- {
- DBGV( "Lambda vector : " << (*it_s_vD).first );
-
- for (int i = 0; i < (int)lambda[(*it_s_vD).first].size(); i++)
+ for(int i=0; i < (int)lambda[(*it_s_vD).first].size(); i++)
{
- DBGV(lambda[(*it_s_vD).first][i]);
+ DBGV(lambda[(*it_s_vD).first][i]);
}
}
-
+
DBGL
return lambda;
}
@@ -615,44 +610,57 @@ TTree* EQ_OneParameter::CreateTMVAInputTree(IsotopicVector TheFreshfuel, double
{
/******Create Input data tree to be interpreted by TMVA::Reader***/
TTree* InputTree = new TTree(fOutput.c_str(), fOutput.c_str());
-
+
vector<float> InputTMVA;
- for (int i = 0 ; i < (int)fMapOfTMVAVariableNames.size() ; i++)
+ for(int i = 0 ; i< (int)fListOfNonZaiTMVAVariables.size() ; i++)
+ InputTMVA.push_back(0);
+ for(int i = 0 ; i< (int)fMapOfTMVAVariableNames.size() ; i++)
InputTMVA.push_back(0);
-
+
float Time = 0;
-
+
IsotopicVector IVInputTMVA;
- map<ZAI , string >::iterator it_ZAI_s;
+ map<ZAI ,string >::iterator it_ZAI_s;
int j = 0;
+
+ for( j = 0; j<fListOfNonZaiTMVAVariables.size(); j++) {
+ InputTree->Branch( (fListOfNonZaiTMVAVariables[j].second).c_str(),
+ &InputTMVA[j], (fListOfNonZaiTMVAVariables[j].second + "/F").c_str());
+ }
- for ( it_ZAI_s = fMapOfTMVAVariableNames.begin() ; it_ZAI_s != fMapOfTMVAVariableNames.end() ; it_ZAI_s++)
+
+
+ for( it_ZAI_s = fMapOfTMVAVariableNames.begin() ; it_ZAI_s != fMapOfTMVAVariableNames.end() ; it_ZAI_s++)
{
InputTree->Branch( ((*it_ZAI_s).second).c_str(), &InputTMVA[j], ((*it_ZAI_s).second + "/F").c_str());
- IVInputTMVA += ((*it_ZAI_s).first) * 1;
+ IVInputTMVA+= ((*it_ZAI_s).first)*1;
j++;
}
-
- if (ThisTime != -1)
- InputTree->Branch("Time" , &Time , "Time/F");
-
+
+ if(ThisTime != -1)
+ InputTree->Branch("Time" ,&Time ,"Time/F");
+
IsotopicVector IVAccordingToUserInfoFile = TheFreshfuel.GetThisComposition(IVInputTMVA);
double Ntot = IVAccordingToUserInfoFile.GetSumOfAll();
- IVAccordingToUserInfoFile = IVAccordingToUserInfoFile / Ntot;
-
+ IVAccordingToUserInfoFile = IVAccordingToUserInfoFile/Ntot;
+
j = 0;
+
+ for( j = 0; j<fListOfNonZaiTMVAVariables.size(); j++) {
+ InputTMVA[j] =fListOfNonZaiTMVAVariables[j].first;
+ }
- for ( it_ZAI_s = fMapOfTMVAVariableNames.begin() ; it_ZAI_s != fMapOfTMVAVariableNames.end() ; it_ZAI_s++)
+ for( it_ZAI_s = fMapOfTMVAVariableNames.begin() ; it_ZAI_s != fMapOfTMVAVariableNames.end() ; it_ZAI_s++)
{
InputTMVA[j] = IVAccordingToUserInfoFile.GetZAIIsotopicQuantity( (*it_ZAI_s).first ) ;
j++;
}
-
+
Time = ThisTime;
InputTree->Fill();
-
+
return InputTree;
-
+
}
//________________________________________________________________________
void EQ_OneParameter::CheckTargetParameterConsistency(map < int , string > StreamListPriority, map < string , double > TargetParameterMin, map < string , double > TargetParameterMax)
@@ -660,31 +668,31 @@ void EQ_OneParameter::CheckTargetParameterConsistency(map < int , string > Strea
map < int , string >::iterator it_i_s;
//Loop on priority order to check if target parameter increases monotously with the material mass
- for ( it_i_s = StreamListPriority.begin(); it_i_s != StreamListPriority.end(); it_i_s++)
- {
+ for( it_i_s = StreamListPriority.begin(); it_i_s != StreamListPriority.end(); it_i_s++)
+ {
double TargetParameterUp = -1.0; //to be sure BUMin is > to BUmax even if BUmin is zero
double TargetParameterDown = 0.0;
- if (TargetParameterMin.find((*it_i_s).second) == TargetParameterMin.end())
+ if(TargetParameterMin.find((*it_i_s).second) == TargetParameterMin.end())
{
- break; //if material is not in map, break the loop
+ break; //if material is not in map, break the loop
}
TargetParameterDown = TargetParameterMin[(*it_i_s).second];
if (TargetParameterDown < 0.0 )
{
- ERROR << "Target parameter evolution should always be positive." << endl;
- ERROR << "TargetParameterDown = " << TargetParameterDown << " is negative " << endl;
- ERROR << "Check the evolution..." << endl;
+ ERROR<< "Target parameter evolution should always be positive." <<endl;
+ ERROR<< "TargetParameterDown = "<< TargetParameterDown<<" is negative "<<endl;
+ ERROR<< "Check the evolution..." <<endl;
exit(1);
- }
+ }
TargetParameterUp = TargetParameterMax[(*it_i_s).second];
if (TargetParameterDown > TargetParameterUp )
- {
- ERROR << "Target parameter evolution as a function of material mass is not monotonous." << endl;
- ERROR << "TargetParameterDown = " << TargetParameterDown << " is greater than TargetParameterUp = " << TargetParameterUp << endl;
- ERROR << "Check the evolution..." << endl;
+ {
+ ERROR<< "Target parameter evolution as a function of material mass is not monotonous." <<endl;
+ ERROR<< "TargetParameterDown = "<< TargetParameterDown<<" is greater than TargetParameterUp = "<< TargetParameterUp<<endl;
+ ERROR<< "Check the evolution..." <<endl;
exit(1);
}
}
@@ -692,14 +700,14 @@ void EQ_OneParameter::CheckTargetParameterConsistency(map < int , string > Strea
//________________________________________________________________________
double EQ_OneParameter::CalculateTargetParameter(IsotopicVector TheFuel, string TargetParameterName)
{
- double ParameterToCalculate = 0;
- if (TargetParameterName == "BurnUpMax") ParameterToCalculate = CalculateBurnUpMax(TheFuel, fModelParameter);
- else if (TargetParameterName == "keffBOC") ParameterToCalculate = CalculateKeffAtBOC(TheFuel);
+ double ParameterToCalculate = 0;
+ if(TargetParameterName=="BurnUpMax") ParameterToCalculate = CalculateBurnUpMax(TheFuel, fModelParameter);
+ else if(TargetParameterName=="keffBOC") ParameterToCalculate = CalculateKeffAtBOC(TheFuel);
else
{
- ERROR << "Target parameter defined in InformationFile ( " << TargetParameterName << " ) doesn't exist" << endl;
- ERROR << "Possible target parameters for the moment are : BurnUpMax and keffBOC." << endl;
- exit(1);
+ ERROR<< "Target parameter defined in InformationFile ( "<<TargetParameterName<<" ) doesn't exist" <<endl;
+ ERROR<< "Possible target parameters for the moment are : BurnUpMax and keffBOC." <<endl;
+ exit(1);
}
return ParameterToCalculate ;
@@ -718,143 +726,158 @@ double EQ_OneParameter::CalculateBurnUpMax(IsotopicVector TheFuel, map<string, d
double OldFinalTimeMinus = 0;
double MaximumBU = fMaximalBU;
double MinimumBU = 0 ;
- double TheFinalTime = BurnupToSecond((MaximumBU - MinimumBU) / 2.);
+ double TheFinalTime = BurnupToSecond((MaximumBU-MinimumBU)/2.);
double OldFinalTimePlus = BurnupToSecond(MaximumBU);
double k_av = 0; //average kinf
double OldPredictedk_av = 0;
-
+
CLASSReader * reader = new CLASSReader( fMapOfTMVAVariableNames );
reader->AddVariable( "Time" );
reader->BookMVA( "MLP method" , fTMVAXMLFilePath );
-
- for (int b = 0; b < NumberOfBatch; b++)
+
+ for(int b = 0;b<NumberOfBatch;b++)
{
- float TheTime = (b + 1) * TheFinalTime / NumberOfBatch;
+ float TheTime = (b+1)*TheFinalTime/NumberOfBatch;
- TTree* InputTree = CreateTMVAInputTree(TheFuel, TheTime);
+ TTree* InputTree = CreateTMVAInputTree(TheFuel,TheTime);
reader->SetInputData( InputTree );
-
+
OldPredictedk_av += reader->EvaluateRegression( "MLP method" )[0];
-
+
delete InputTree;
}
OldPredictedk_av /= NumberOfBatch;
-
+
//Algorithm control
int count = 0;
int MaximumLoopCount = 500;
do
{
- if (count > MaximumLoopCount )
+ if(count > MaximumLoopCount )
{
ERROR << "CRITICAL ! Can't manage to predict burnup\nHint : Try to increase the precision on k effective using :\n YourEQM_MLP_Kinf->SetPCMPrecision(pcm); with pcm the precision in pcm (default 10) REDUCE IT\n If this message still appear mail to leniau@subatech.in2p3.fr\nor nicolas.thiolliere@subatech.in2p3.fr " << endl;
exit(1);
}
-
- if ( (OldPredictedk_av - KThreshold) > 0) //The burnup can be increased
+
+ if( (OldPredictedk_av-KThreshold) > 0) //The burnup can be increased
{
OldFinalTimeMinus = TheFinalTime;
- TheFinalTime = TheFinalTime + fabs(OldFinalTimePlus - TheFinalTime) / 2.;
-
- if (SecondToBurnup(TheFinalTime) >= (MaximumBU - MaximumBU * GetTargetParameterStDev() ) )
- { delete reader; return MaximumBU; }
+ TheFinalTime = TheFinalTime + fabs(OldFinalTimePlus - TheFinalTime)/2.;
+
+ if(SecondToBurnup(TheFinalTime) >= (MaximumBU-MaximumBU*GetTargetParameterStDev() ) )
+ { delete reader; return MaximumBU; }
}
-
- else if ( (OldPredictedk_av - KThreshold) < 0) //The burnup is too high
+
+ else if( (OldPredictedk_av-KThreshold) < 0)//The burnup is too high
{
OldFinalTimePlus = TheFinalTime;
- TheFinalTime = TheFinalTime - fabs(OldFinalTimeMinus - TheFinalTime) / 2.;
- if ( SecondToBurnup(TheFinalTime) < (MaximumBU - MinimumBU) / 2.*GetTargetParameterStDev() )
- { delete reader; return 0; }
+ TheFinalTime = TheFinalTime - fabs(OldFinalTimeMinus-TheFinalTime)/2.;
+ if( SecondToBurnup(TheFinalTime) < (MaximumBU-MinimumBU)/2.*GetTargetParameterStDev() )
+ { delete reader; return 0; }
}
-
+
k_av = 0;
- for (int b = 0; b < NumberOfBatch; b++)
+ for(int b = 0;b<NumberOfBatch;b++)
{
- float TheTime = (b + 1) * TheFinalTime / NumberOfBatch;
- TTree* InputTree = CreateTMVAInputTree(TheFuel, TheTime);
+ float TheTime = (b+1)*TheFinalTime/NumberOfBatch;
+ TTree* InputTree = CreateTMVAInputTree(TheFuel,TheTime);
reader->SetInputData( InputTree );
-
+
k_av += reader->EvaluateRegression("MLP method")[0];
-
+
delete InputTree;
}
- k_av /= NumberOfBatch;
+ k_av/= NumberOfBatch;
//cout<<SecondToBurnup(TheFinalTime)<<" ";
OldPredictedk_av = k_av;
count++;
-//std::clog << "-> " << k_av << "\t\t(" << count << ") \t [" << TheFinalTime << "]" << "\t" << OldPredictedk_av-KThreshold << "\t" << GetPCMPrecision() << std::endl;
- } while ( fabs(OldPredictedk_av - KThreshold) > GetPCMPrecision() ) ;
-
+//std::clog << "-> " << k_av << "\t\t(" << count << ") \t [" << TheFinalTime << "]" << "\t" << OldPredictedk_av-KThreshold << "\t" << GetPCMPrecision() << std::endl;
+ } while( fabs(OldPredictedk_av-KThreshold) > GetPCMPrecision() ) ;
+
delete reader;
//cout<<endl;
return SecondToBurnup(TheFinalTime);
}
//________________________________________________________________________
-double EQ_OneParameter::CalculateKeffAtBOC(IsotopicVector FreshFuel)
-{
- CLASSReader * reader = new CLASSReader( fMapOfTMVAVariableNames );
- reader->BookMVA( "MLP method" , fTMVAXMLFilePath );
-
- TTree* InputTree = CreateTMVAInputTree(FreshFuel, -1) ;
- reader->SetInputData( InputTree );
-
- double keff = reader->EvaluateRegression( "MLP method" )[0];
-
- delete InputTree;
-
- return keff;
-}
+double EQ_OneParameter::CalculateKeffAtBOC(IsotopicVector FreshFuel)
+{
+ double keff(-1);
+ double TimeOfInterest(-1);
+ CLASSReader* reader;
+ if(fListOfNonZaiTMVAVariables.size()>0){
+ vector<string> VectorOfAllTMVAVariableNames;
+ for(unsigned int j = 0; j<fListOfNonZaiTMVAVariables.size(); j++) {
+ VectorOfAllTMVAVariableNames.push_back(fListOfNonZaiTMVAVariables[j].second);
+ }
+ for(map<ZAI,string>::iterator it = fMapOfTMVAVariableNames.begin(); it != fMapOfTMVAVariableNames.end(); it++) {
+ VectorOfAllTMVAVariableNames.push_back(it->second);
+ }
+ reader = new CLASSReader( VectorOfAllTMVAVariableNames );
+ reader->AddVariable("Time");
+ reader->BookMVA( "MLP method" , fTMVAXMLFilePath );
+ TimeOfInterest=0;//0 because BOC
+ }
+ else{
+ reader = new CLASSReader( fMapOfTMVAVariableNames );
+ reader->BookMVA( "MLP method" , fTMVAXMLFilePath );
+ TimeOfInterest=-1;
+ }
+ TTree* InputTree = CreateTMVAInputTree(FreshFuel,TimeOfInterest) ;
+ reader->SetInputData( InputTree );
+ keff = reader->EvaluateRegression( "MLP method" )[0];
+ delete InputTree;
+ return keff;
+}
//________________________________________________________________________
void EQ_OneParameter::ReadNFO()
{
DBGL
ifstream NFO(fTMVANFOFilePath.c_str());
-
- if (!NFO)
+
+ if(!NFO)
{
ERROR << "Can't find/open file " << fTMVANFOFilePath << endl;
exit(0);
}
-
+
do
{
string line;
- getline(NFO, line);
-
+ getline(NFO,line);
+
EQ_OneParameter::ReadLine(line);
-
- } while (!NFO.eof());
-
+
+ } while(!NFO.eof());
+
DBGL
}
//________________________________________________________________________
void EQ_OneParameter::ReadLine(string line)
{
DBGL
-
+
if (!freaded)
{
int pos = 0;
string keyword = tlc(StringLine::NextWord(line, pos, ' '));
-
+
map<string, EQOP_MthPtr>::iterator it = fKeyword.find(keyword);
-
- if (it != fKeyword.end())
+
+ if(it != fKeyword.end())
(this->*(it->second))( line );
-
+
freaded = true;
ReadLine(line);
-
+
}
-
+
freaded = false;
-
+
DBGL
}
//________________________________________________________________________
-void EQ_OneParameter::LoadKeyword()
+void EQ_OneParameter::LoadKeyword()
{
DBGL
fKeyword.insert( pair<string, EQOP_MthPtr>( "k_zail", & EQ_OneParameter::ReadZAIlimits) );
@@ -865,13 +888,13 @@ void EQ_OneParameter::LoadKeyword()
fKeyword.insert( pair<string, EQOP_MthPtr>( "k_list", & EQ_OneParameter::ReadList) );
fKeyword.insert( pair<string, EQOP_MthPtr>( "k_specpower", & EQ_OneParameter::ReadSpecificPower) );
if (fUseTMVAPredictor) fKeyword.insert( pair<string, EQOP_MthPtr>( "k_zainame", & EQ_OneParameter::ReadZAIName) );
- fKeyword.insert( pair<string, EQOP_MthPtr>( "k_maxburnup", & EQ_OneParameter::ReadMaxBurnUp) );
+ fKeyword.insert( pair<string, EQOP_MthPtr>( "k_maxburnup", & EQ_OneParameter::ReadMaxBurnUp) );
if (fUseTMVAPredictor) fKeyword.insert( pair<string, EQOP_MthPtr>( "k_targetparameter", & EQ_OneParameter::ReadTargetParameter) );
if (fUseTMVAPredictor) fKeyword.insert( pair<string, EQOP_MthPtr>( "k_predictortype", & EQ_OneParameter::ReadPredictorType) );
if (fUseTMVAPredictor) fKeyword.insert( pair<string, EQOP_MthPtr>( "k_output", & EQ_OneParameter::ReadOutput) );
- fKeyword.insert( pair<string, EQOP_MthPtr>( "k_buffer", & EQ_OneParameter::ReadBuffer) );
- if (fUseTMVAPredictor) fKeyword.insert( pair<string, EQOP_MthPtr>( "k_modelparameter", & EQ_OneParameter::ReadModelParameter) );
- if (fUseTMVAPredictor) fKeyword.insert( pair<string, EQOP_MthPtr>( "k_targetparameterstdev", & EQ_OneParameter::ReadTargetParameterStDev) );
+ fKeyword.insert( pair<string, EQOP_MthPtr>( "k_buffer", & EQ_OneParameter::ReadBuffer) );
+ if (fUseTMVAPredictor) fKeyword.insert( pair<string, EQOP_MthPtr>( "k_modelparameter", & EQ_OneParameter::ReadModelParameter) );
+ if (fUseTMVAPredictor) fKeyword.insert( pair<string, EQOP_MthPtr>( "k_targetparameterstdev", & EQ_OneParameter::ReadTargetParameterStDev) );
DBGL
}
@@ -881,16 +904,16 @@ void EQ_OneParameter::ReadType(const string &line)
DBGL
int pos = 0;
string keyword = tlc(StringLine::NextWord(line, pos, ' '));
- if ( keyword != "k_fuel" && keyword != "k_reactor" ) // Check the keyword
+ if( keyword != "k_fuel" && keyword != "k_reactor" ) // Check the keyword
{
ERROR << " Bad keyword : " << keyword << " Not found !" << endl;
exit(1);
}
- if ( keyword == "k_fuel" )
+ if( keyword == "k_fuel" )
fDBFType = StringLine::NextWord(line, pos, ' ');
- else if ( keyword == "k_reactor" )
+ else if( keyword == "k_reactor" )
fDBRType = StringLine::NextWord(line, pos, ' ');
-
+
DBGL
}
//________________________________________________________________________
@@ -899,26 +922,26 @@ void EQ_OneParameter::ReadZAIlimits(const string &line)
DBGL
int pos = 0;
string keyword = tlc(StringLine::NextWord(line, pos, ' '));
- if ( keyword != "k_zail" ) // Check the keyword
+ if( keyword != "k_zail" ) // Check the keyword
{
ERROR << " Bad keyword : \"k_zail\" not found !" << endl;
exit(1);
}
-
+
int Z = atoi(StringLine::NextWord(line, pos, ' ').c_str());
int A = atoi(StringLine::NextWord(line, pos, ' ').c_str());
int I = atoi(StringLine::NextWord(line, pos, ' ').c_str());
-
+
double downLimit = atof(StringLine::NextWord(line, pos, ' ').c_str());
double upLimit = atof(StringLine::NextWord(line, pos, ' ').c_str());
-
+
if (upLimit < downLimit)
{
double tmp = upLimit;
upLimit = downLimit;
downLimit = tmp;
}
- fZAILimits.insert(pair<ZAI, pair<double, double> >(ZAI(Z, A, I), pair<double, double>(downLimit, upLimit)));
+ fZAILimits.insert(pair<ZAI, pair<double, double> >(ZAI(Z,A,I), pair<double,double>(downLimit, upLimit)));
DBGL
}
//________________________________________________________________________
@@ -927,18 +950,18 @@ void EQ_OneParameter::ReadList(const string &line)
DBGL
int pos = 0;
string keyword = tlc(StringLine::NextWord(line, pos, ' '));
- if ( keyword != "k_list" ) // Check the keyword
+ if( keyword != "k_list" ) // Check the keyword
{
ERROR << " Bad keyword : \"k_list\" not found !" << endl;
exit(1);
}
- string ListName = StringLine::NextWord(line, pos, ' ');
+ string ListName= StringLine::NextWord(line, pos, ' ');
int Z = atoi(StringLine::NextWord(line, pos, ' ').c_str());
int A = atoi(StringLine::NextWord(line, pos, ' ').c_str());
int I = atoi(StringLine::NextWord(line, pos, ' ').c_str());
double Q = atof(StringLine::NextWord(line, pos, ' ').c_str());
fStreamList[ListName].Add(Z, A, I, Q);
-
+
DBGL
}
//________________________________________________________________________
@@ -947,12 +970,12 @@ void EQ_OneParameter::ReadEqMinFraction(const string &line)
DBGL
int pos = 0;
string keyword = tlc(StringLine::NextWord(line, pos, ' '));
- if ( keyword != "k_massfractionmin" ) // Check the keyword
+ if( keyword != "k_massfractionmin" ) // Check the keyword
{
ERROR << " Bad keyword : \"k_massfractionmin\" not found !" << endl;
exit(1);
}
- string ListName = StringLine::NextWord(line, pos, ' ');
+ string ListName= StringLine::NextWord(line, pos, ' ');
double Q = atof(StringLine::NextWord(line, pos, ' ').c_str());
fStreamListEqMMassFractionMin[ListName] = Q;
@@ -965,12 +988,12 @@ void EQ_OneParameter::ReadEqMaxFraction(const string &line)
DBGL
int pos = 0;
string keyword = tlc(StringLine::NextWord(line, pos, ' '));
- if ( keyword != "k_massfractionmax" ) // Check the keyword
+ if( keyword != "k_massfractionmax" ) // Check the keyword
{
ERROR << " Bad keyword : \"k_massfractionmax\" not found !" << endl;
exit(1);
}
- string ListName = StringLine::NextWord(line, pos, ' ');
+ string ListName= StringLine::NextWord(line, pos, ' ');
double Q = atof(StringLine::NextWord(line, pos, ' ').c_str());
fStreamListEqMMassFractionMax[ListName] = Q;
@@ -983,38 +1006,38 @@ void EQ_OneParameter::ReadSpecificPower(const string &line)
DBGL
int pos = 0;
string keyword = tlc(StringLine::NextWord(line, pos, ' '));
- if ( keyword != "k_specpower") // Check the keyword
+ if( keyword != "k_specpower") // Check the keyword
{
ERROR << " Bad keyword : \"k_specpower\" Not found !" << endl;
exit(1);
}
-
+
fSpecificPower = atof(StringLine::NextWord(line, pos, ' ').c_str());
-
+
DBGL
}
//________________________________________________________________________
void EQ_OneParameter::ReadZAIName(const string &line)
{
DBGL
-
+
int pos = 0;
string keyword = tlc(StringLine::NextWord(line, pos, ' '));
- if ( keyword != "k_zainame" ) // Check the keyword
+ if( keyword != "k_zainame" ) // Check the keyword
{
ERROR << " Bad keyword : \"k_zainame\" not found !" << endl;
exit(1);
}
-
+
int Z = atoi(StringLine::NextWord(line, pos, ' ').c_str());
int A = atoi(StringLine::NextWord(line, pos, ' ').c_str());
int I = atoi(StringLine::NextWord(line, pos, ' ').c_str());
-
+
string name = StringLine::NextWord(line, pos, ' ');
+
+ fMapOfTMVAVariableNames.insert( pair<ZAI,string>( ZAI(Z, A, I), name ) );
- fMapOfTMVAVariableNames.insert( pair<ZAI, string>( ZAI(Z, A, I), name ) );
-
- DBGL
+ DBGL
}
//________________________________________________________________________
void EQ_OneParameter::ReadMaxBurnUp(const string &line)
@@ -1022,12 +1045,12 @@ void EQ_OneParameter::ReadMaxBurnUp(const string &line)
DBGL
int pos = 0;
string keyword = tlc(StringLine::NextWord(line, pos, ' '));
- if ( keyword != "k_maxburnup" ) // Check the keyword
+ if( keyword != "k_maxburnup" ) // Check the keyword
{
ERROR << " Bad keyword : \"k_maxburnup\" not found !" << endl;
exit(1);
}
-
+
fMaximalBU = atof(StringLine::NextWord(line, pos, ' ').c_str());
DBGL
@@ -1039,12 +1062,12 @@ void EQ_OneParameter::ReadTargetParameter(const string &line)
DBGL
int pos = 0;
string keyword = tlc(StringLine::NextWord(line, pos, ' '));
- if ( keyword != "k_targetparameter" ) // Check the keyword
+ if( keyword != "k_targetparameter" ) // Check the keyword
{
ERROR << " Bad keyword : \"k_targetparameter\" not found !" << endl;
exit(1);
}
-
+
fTargetParameter = StringLine::NextWord(line, pos, ' ');
DBGL
@@ -1054,132 +1077,163 @@ void EQ_OneParameter::ReadTargetParameter(const string &line)
void EQ_OneParameter::ReadPredictorType(const string &line)
{
DBGL
-
+
int pos = 0;
string keyword = tlc(StringLine::NextWord(line, pos, ' '));
- if ( keyword != "k_predictortype" ) // Check the keyword
+ if( keyword != "k_predictortype" ) // Check the keyword
{
ERROR << " Bad keyword : \"k_predictortype\" not found !" << endl;
exit(1);
}
-
+
fPredictorType = StringLine::NextWord(line, pos, ' ');
-
- DBGL
+
+ DBGL
}
//________________________________________________________________________
void EQ_OneParameter::ReadOutput(const string &line)
{
DBGL
-
+
int pos = 0;
string keyword = tlc(StringLine::NextWord(line, pos, ' '));
- if ( keyword != "k_output" ) // Check the keyword
+ if( keyword != "k_output" ) // Check the keyword
{
ERROR << " Bad keyword : \"k_output\" not found !" << endl;
exit(1);
}
-
+
fOutput = StringLine::NextWord(line, pos, ' ');
-
- DBGL
+
+ DBGL
}
//________________________________________________________________________
void EQ_OneParameter::ReadBuffer(const string &line)
{
DBGL
-
+
int pos = 0;
string keyword = tlc(StringLine::NextWord(line, pos, ' '));
- if ( keyword != "k_buffer" ) // Check the keyword
+ if( keyword != "k_buffer" ) // Check the keyword
{
ERROR << " Bad keyword : \"k_buffer\" not found !" << endl;
exit(1);
}
-
+
fBuffer = StringLine::NextWord(line, pos, ' ');
-
- DBGL
+
+ DBGL
}
//________________________________________________________________________
void EQ_OneParameter::ReadModelParameter(const string &line)
{
DBGL
-
+
int pos = 0;
string keyword = tlc(StringLine::NextWord(line, pos, ' '));
- if ( keyword != "k_modelparameter" ) // Check the keyword
+ if( keyword != "k_modelparameter" ) // Check the keyword
{
ERROR << " Bad keyword : \"k_modelparameter\" not found !" << endl;
exit(1);
}
-
+
keyword = StringLine::NextWord(line, pos, ' ');
fModelParameter[keyword] = -1;
-
- DBGL
+
+ DBGL
+}
+//________________________________________________________________________
+void EQ_OneParameter::ReadNonZaiTMVAVariables(const string &line)
+{
+ DBGL
+
+ int pos = 0;
+ string keyword = tlc(StringLine::NextWord(line, pos, ' '));
+ if( keyword != "k_nonZAIforTMVA" ) // Check the keyword
+ {
+ ERROR << " Bad keyword : \"k_nonZAIforTMVA\" not found !" << endl;
+ exit(1);
+ }
+
+ keyword = StringLine::NextWord(line, pos, ' ');
+ fListOfNonZaiTMVAVariables.push_back(make_pair(-1.0,keyword));
+
+ DBGL
+}
+//________________________________________________________________________
+void EQ_OneParameter::SetNonZaiTMVAVariable(string snZP, double dnZP)
+{
+ DBGL
+ for(unsigned int j=0;j<fListOfNonZaiTMVAVariables.size();j++){
+ if(fListOfNonZaiTMVAVariables[j].second==snZP){
+ fListOfNonZaiTMVAVariables[j].first=dnZP;
+ return;
+ }
+ }
+ fListOfNonZaiTMVAVariables.push_back(make_pair(dnZP,snZP));
+
+ DBGL
}
-
//________________________________________________________________________
void EQ_OneParameter::ReadTargetParameterStDev(const string &line)
{
DBGL
-
+
int pos = 0;
string keyword = tlc(StringLine::NextWord(line, pos, ' '));
- if ( keyword != "k_targetparameterstdev" ) // Check the keyword
+ if( keyword != "k_targetparameterstdev" ) // Check the keyword
{
ERROR << " Bad keyword : \"k_targetparameterstdev\" not found !" << endl;
exit(1);
}
-
+
fTargetParameterStDev = atof(StringLine::NextWord(line, pos, ' ').c_str());;
-
- DBGL
+
+ DBGL
}
//________________________________________________________________________
void EQ_OneParameter::PrintInfo()
{
- INFO << "Reactor Type : " << fDBRType << endl;
- INFO << "Fuel Type : " << fDBFType << endl;
- INFO << "Specific Power [W/g]: " << fSpecificPower << endl;
-
+ INFO << "Reactor Type : "<< fDBRType << endl;
+ INFO << "Fuel Type : "<< fDBFType << endl;
+ INFO << "Specific Power [W/g]: "<< fSpecificPower << endl;
+
map < string , IsotopicVector >::iterator it_s_IV;
map < string , double >::iterator it_s_D;
- for ( it_s_IV = fStreamList.begin(); it_s_IV != fStreamList.end(); it_s_IV++)
- {
- INFO << (* it_s_IV).first << " (Z A I) :" << endl;
- map<ZAI , double >::iterator it1;
- map<ZAI , double > fMap1 = fStreamList[(* it_s_IV).first].GetIsotopicQuantity();
- for (it1 = fMap1.begin() ; it1 != fMap1.end() ; it1++)
- INFO << (*it1).first.Z() << " " << (*it1).first.A() << " " << (*it1).first.I() << endl;
+ for( it_s_IV = fStreamList.begin(); it_s_IV != fStreamList.end(); it_s_IV++)
+ {
+ INFO <<(* it_s_IV).first<<" (Z A I) :" << endl;
+ map<ZAI ,double >::iterator it1;
+ map<ZAI ,double > fMap1 = fStreamList[(* it_s_IV).first].GetIsotopicQuantity();
+ for(it1 = fMap1.begin() ; it1 != fMap1.end() ; it1++)
+ INFO << (*it1).first.Z() <<" "<< (*it1).first.A() <<" "<< (*it1).first.I() << endl;
}
- INFO << "Minimum fraction in the fuel for each material : " << endl;
- for ( it_s_D = fStreamListEqMMassFractionMin.begin(); it_s_D != fStreamListEqMMassFractionMin.end(); it_s_D++)
+ INFO<<"Minimum fraction in the fuel for each material : "<<endl;
+ for( it_s_D = fStreamListEqMMassFractionMin.begin(); it_s_D != fStreamListEqMMassFractionMin.end(); it_s_D++)
{
- INFO << (* it_s_D).first << " " << fStreamListEqMMassFractionMin[(* it_s_D).first] << endl;
- }
+ INFO <<(* it_s_D).first<<" "<<fStreamListEqMMassFractionMin[(* it_s_D).first]<<endl;
+ }
- INFO << "Maximum fraction in the fuel for each material : " << endl;
- for ( it_s_D = fStreamListEqMMassFractionMax.begin(); it_s_D != fStreamListEqMMassFractionMax.end(); it_s_D++)
+ INFO<<"Maximum fraction in the fuel for each material : "<<endl;
+ for( it_s_D = fStreamListEqMMassFractionMax.begin(); it_s_D != fStreamListEqMMassFractionMax.end(); it_s_D++)
{
- INFO << (* it_s_D).first << " " << fStreamListEqMMassFractionMax[(* it_s_D).first] << endl;
- }
+ INFO <<(* it_s_D).first<<" "<<fStreamListEqMMassFractionMax[(* it_s_D).first]<<endl;
+ }
- INFO << "ZAI limits (validity domain)[prop in fresh fuel] (Z A I min max) :" << endl;
- for (map< ZAI, pair<double, double> >::iterator Domain_it = fZAILimits.begin(); Domain_it != fZAILimits.end(); Domain_it++)
- {
+ INFO<<"ZAI limits (validity domain)[prop in fresh fuel] (Z A I min max) :"<<endl;
+ for (map< ZAI,pair<double,double> >::iterator Domain_it = fZAILimits.begin(); Domain_it != fZAILimits.end(); Domain_it++)
+ {
double ThatZAIMin = Domain_it->second.first;
double ThatZAIMax = Domain_it->second.second;
int Z = Domain_it->first.Z();
int A = Domain_it->first.A();
int I = Domain_it->first.I();
- INFO << ThatZAIMin << " < ZAI (" << Z << " " << A << " " << I << ")" << " < " << ThatZAIMax << endl;
+ INFO <<ThatZAIMin<<" < ZAI ("<< Z<< " " << A << " " << I<<")"<<" < "<<ThatZAIMax<< endl;
}
}
diff --git a/source/Model/Equivalence/EQ_OneParameter.hxx b/source/Model/Equivalence/EQ_OneParameter.hxx
index 6922e56bb8baff759a96715afce03fef4261960c..392eccd44d7f0140aaa8090aa97b9075e54e528a 100644
--- a/source/Model/Equivalence/EQ_OneParameter.hxx
+++ b/source/Model/Equivalence/EQ_OneParameter.hxx
@@ -12,6 +12,7 @@
@version 3.0
*/
+#include "EquivalenceModel.hxx"
#include "IsotopicVector.hxx"
#include <math.h>
#include "TTree.h"
@@ -111,8 +112,11 @@ class EQ_OneParameter : public EquivalenceModel
double GetStreamListEqMMassFractionMin(string keyword){return fStreamListEqMMassFractionMin[keyword] ;}
double GetPCMPrecision(){return fPCMprecision/1e5;}//!< Get the precision on @f$\langle k \rangle@f$ prediction []. Neural network predictor constructors
- void SetModelParameter(string sMP, double dMP) { fModelParameter[sMP] = dMP; } //!< Set Model Parameters precised in NFO file
- map<string, double> GetModelParameter() { return fModelParameter; } //!< Get Model Parameters precised in NFO file
+ void SetModelParameter(string sMP, double dMP) { fModelParameter[sMP] = dMP; } //!< Set Model Parameters precised in NFO file
+ map<string, double> GetModelParameter() { return fModelParameter; } //!< Get Model Parameters precised in NFO file
+
+ void SetNonZaiTMVAVariable(string snZP, double dnZP); //!< Set NonZaiTMVAVariables
+ vector< pair<double, string> > GetNonZaiTMVAVariables() { return fListOfNonZaiTMVAVariables; } //!< Get NonZaiTMVAVariables
void SetMaxIterration(int val) { fMaxIterration = val; } //!< Max iteration in build fuel algorithm
void SetTargetParameterStDev(double TPSD){fTargetParameterStDev = TPSD;} //!< Set the precision on Target Parameter
@@ -176,7 +180,13 @@ class EQ_OneParameter : public EquivalenceModel
*/
void ReadTargetParameter(const string &line);
//}
-
+ //{
+ /// ReadNonZaiTMVAVariables : read the NonZai variables for the predictor (ex : Nbatch, Specific power)
+ /*!
+ \param line : line suppossed to contain the NonZai variables for TMVA starts with "k_nonZAIforTMVA" keyword
+ */
+ void ReadNonZaiTMVAVariables(const string &line);
+ //}
//{
/// ReadOutput : read the output type of the predictor (ex : kinf)
/*!
@@ -244,6 +254,7 @@ class EQ_OneParameter : public EquivalenceModel
string fBuffer ; //!< Name of material used as buffer in fuel
map<string, double> fModelParameter ; ///< Map of equivalence model parameter
+ vector< pair<double, string> > fListOfNonZaiTMVAVariables ; ///!< List of TMVA input variable names that are not ZAIs
map<ZAI, string> fMapOfTMVAVariableNames; //!< List of TMVA input variable names (read from fMLPInformationFile ) , name depends on the training step
diff --git a/source/Model/XS/XSM_SFR.cxx b/source/Model/XS/XSM_SFR.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..61c778cfe1d3d030ab510fd23cceed9ec852f35d
--- /dev/null
+++ b/source/Model/XS/XSM_SFR.cxx
@@ -0,0 +1,568 @@
+#include <dirent.h>
+#include <errno.h>
+#include <sstream>
+#include <string>
+#include <iostream>
+#include <fstream>
+#include <algorithm>
+#include <cmath>
+
+#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 "XSModel.hxx"
+#include "XSM_SFR.hxx"
+#include "CLASSLogger.hxx"
+#include "CLASSMethod.hxx"
+#include "CLASSReader.hxx"
+#include "external/StringLine.hxx"
+
+//________________________________________________________________________
+//
+// XSM_SFR
+//
+//________________________________________________________________________
+XSM_SFR::XSM_SFR(string TMVA_Weight_Directory,map<string,double> FixedParameters,string InformationFile, bool IsTimeStep):XSModel(new CLASSLogger("XSM_SFR.log"))
+{
+
+ fIsStepTime = IsTimeStep;
+ fTMVAWeightFolder = TMVA_Weight_Directory;
+
+ fInformationFile = fTMVAWeightFolder+InformationFile;
+
+ GetMLPWeightFiles();
+
+ INFO << "__A cross section interpolator using" << endl;
+ INFO << "Neural Network has been define__" << endl;
+ INFO << " \t His TMVA folder is : \" " << fTMVAWeightFolder << "\"" << endl;
+
+ LoadKeyword();
+ ReadNFO();
+
+ int number_of_elements = fTMVAVariableNames.size();
+ double default_value = -1;
+ vector<double> v(number_of_elements, default_value);
+ fTMVAFixedVariableValues = v;
+ for(unsigned int i=0;i<fTMVAVariableNames.size();i++){
+ if(FixedParameters.find(fTMVAVariableNames[i]) != FixedParameters.end() ){
+ fTMVAFixedVariableValues[i]=FixedParameters[fTMVAVariableNames[i]];
+ fTMVAFixedVariable[i]=true;
+ }
+ }
+}
+//________________________________________________________________________
+XSM_SFR::XSM_SFR(CLASSLogger* Log,string TMVA_Weight_Directory,map<string,double> FixedParameters,string InformationFile, bool IsTimeStep):XSModel(Log)
+{
+
+ fIsStepTime = IsTimeStep;
+ fTMVAWeightFolder = TMVA_Weight_Directory;
+
+ fInformationFile = fTMVAWeightFolder+InformationFile;
+
+ GetMLPWeightFiles();
+
+ INFO << "__A cross section interpolator using" << endl;
+ INFO << "Neural Network has been define__" << endl;
+ INFO << " \t His TMVA folder is : \" " << fTMVAWeightFolder << "\"" << endl;
+
+ LoadKeyword();
+ ReadNFO();
+
+ int number_of_elements = fTMVAVariableNames.size();
+ double default_value = -1;
+ vector<double> v(number_of_elements, default_value);
+ fTMVAFixedVariableValues = v;
+ for(unsigned int i=0;i<fTMVAVariableNames.size();i++){
+ if(FixedParameters.find(fTMVAVariableNames[i]) != FixedParameters.end() ){
+ fTMVAFixedVariableValues[i]=FixedParameters[fTMVAVariableNames[i]];
+ fTMVAFixedVariable[i]=true;
+ }
+ }
+}
+//________________________________________________________________________
+XSM_SFR::~XSM_SFR()
+{
+ DBGL
+ fTMVAVariableNames.clear();
+ fDKeyword.clear();
+ DBGL
+}
+//________________________________________________________________________
+void XSM_SFR::SetFixedVariablesValues(map<string,double> FixedParameters)
+{
+ int number_of_elements = fTMVAVariableNames.size();
+ double default_value = -1;
+ vector<double> v(number_of_elements, default_value);
+ fTMVAFixedVariableValues=v;
+ for(unsigned int i=0;i<fTMVAVariableNames.size();i++){
+ if(FixedParameters.find(fTMVAVariableNames[i]) != FixedParameters.end() ){
+ fTMVAFixedVariableValues[i]=FixedParameters[fTMVAVariableNames[i]];
+ fTMVAFixedVariable[i]=true;
+ }
+ }
+}
+//________________________________________________________________________
+void XSM_SFR::LoadKeyword()
+{
+ DBGL
+ fDKeyword.insert( pair<string, XS_SFR_DMthPtr>( "k_timestep", & XSM_SFR::ReadTimeSteps));
+ fDKeyword.insert( pair<string, XS_SFR_DMthPtr>( "k_inputparameter", & XSM_SFR::ReadInputParameter) );
+ DBGL
+}
+//________________________________________________________________________
+void XSM_SFR::ReadInputParameter(const string &line)
+{
+ DBGL
+ int pos = 0;
+ string keyword = tlc(StringLine::NextWord(line, pos, ' '));
+ if( keyword != "k_inputparameter" ) // Check the keyword
+ {
+ ERROR << " Bad keyword : \"k_inputparameter\" not found !" << endl;
+ exit(1);
+ }
+
+ string name = StringLine::NextWord(line, pos, ' ');
+
+ fTMVAVariableNames.push_back( name );
+
+ string IsFixed = StringLine::NextWord(line, pos, ' ');
+
+ if(IsFixed=="F"){
+ fTMVAFixedVariable.push_back(true);
+ }
+ else{
+ fTMVAFixedVariable.push_back(false);
+ }
+
+ DBGL
+}
+//________________________________________________________________________
+void XSM_SFR::ReadTimeSteps(const string &line)
+{
+ DBGL
+ int pos = 0;
+ string keyword = tlc(StringLine::NextWord(line, pos, ' '));
+ if( keyword != "k_timestep" ) // Check the keyword
+ {
+ ERROR << " Bad keyword : \"k_timestep\" not found !" << endl;
+ exit(1);
+ }
+
+ while( pos < (int)line.size() )
+ fMLP_Time.push_back( atof( (StringLine::NextWord(line,pos,' ')).c_str() ));
+ DBGL
+}
+//________________________________________________________________________
+void XSM_SFR::ReadLine(string line)
+{
+ DBGL
+
+ int pos = 0;
+ string keyword = tlc(StringLine::NextWord(line, pos, ' '));
+
+ map<string, XS_SFR_DMthPtr>::iterator it = fDKeyword.find(keyword);
+ if(it != fDKeyword.end())
+ (this->*(it->second))( line );
+
+ DBGL
+}
+//________________________________________________________________________
+void XSM_SFR::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)
+ {
+ ERROR << " Reading error for TMVA weight folder " << fTMVAWeightFolder.c_str() << " : " << strerror(errno) << endl;
+ 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' && FileName[0] != '.' )
+ fWeightFiles.push_back(FileName);
+
+ }
+ }
+ DBGL
+}
+//________________________________________________________________________
+//________________________________________________________________________
+//
+// Time (MLP take time as parameter)
+//________________________________________________________________________
+//________________________________________________________________________
+void XSM_SFR::ReadWeightFile(string Filename, int &Z, int &A, int &I, int &Reaction)
+{
+ int tempZ = -1;
+ int tempA = -1;
+ int tempI = -1;
+ int tempReaction = -1;
+
+ size_t found = Filename.find("XS");
+ string NameJOB;
+ NameJOB = Filename.substr(found);
+ int pos = 0;
+
+ StringLine::NextWord(NameJOB, pos, '_');
+ tempZ = atof( (StringLine::NextWord(NameJOB,pos,'_') ).c_str() );
+ tempA = atof( (StringLine::NextWord(NameJOB,pos,'_') ).c_str() );
+ tempI = 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")
+ tempReaction = 0;
+ if(sReaction == "cap")
+ tempReaction = 1;
+ if(sReaction == "n2n")
+ tempReaction = 2;
+
+ if(tempZ<= 0 || tempA<= 0 || tempI<0 || tempReaction == -1)
+ {
+ ERROR << " wrong TMVA weight format " << endl;
+ exit(0);
+ }
+ else{
+ Z = tempZ;
+ A = tempA;
+ I = tempI;
+ Reaction = tempReaction;
+ }
+
+}
+//________________________________________________________________________
+TTree* XSM_SFR::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(unsigned int i = 0 ; i< (int)fTMVAVariableNames.size() ; i++){
+ InputTMVA.push_back(0);
+ }
+ float Time = 0;
+ IsotopicVector IVInputTMVA;
+ vector<ZAI*> myZAIs;
+
+ for(unsigned int j = 0 ; j< fTMVAVariableNames.size() ; j++)
+ {
+ InputTree->Branch( fTMVAVariableNames[j].c_str() ,&InputTMVA[j], (fTMVAVariableNames[j] + "/F").c_str());
+
+ string ZAIname=fTMVAVariableNames[j];
+ // For this function to work, ZAI name in the xml file have to be AAAName*
+ // AAA is A, the number of nucleon
+ // Name is the symbol of the element (with letters : U, Pu, Am, Cm...) and give Z
+ // * => I=1, no * => I=0
+ int Z = -1;
+ int A = -1;
+ int I = 0;
+
+ if(ZAIname.back()=='*'){//Check last charcter to see if it's *
+ I = 1;
+ ZAIname.pop_back();
+ }
+ if (isdigit(ZAIname.c_str()[0])){//read all numbers in the name to find A
+ A = atoi (ZAIname.c_str());
+ }
+ //remove all numbers in the ZAIname
+ ZAIname.erase(remove_if(ZAIname.begin(), ZAIname.end(),[](unsigned char c){ return isdigit(c); }),ZAIname.end());
+ //remaing (i.e. the text in ZAIname) is the element symbol
+ if(ZAIname=="Am"){
+ Z = 95;
+ }
+ else if(ZAIname=="Pu"){
+ Z = 94;
+ }
+ else if(ZAIname=="U"){
+ Z = 92;
+ }
+ if((Z<= 0 || A<= 0 )&& fTMVAFixedVariable[j]==false)
+ {
+ ERROR << " wrong TMVA weight format " << endl;
+ exit(0);
+ }
+ ZAI* myZAI= new ZAI(Z,A,I);
+ myZAIs.push_back(myZAI);
+
+ if(fTMVAFixedVariable[j]==false){
+ IVInputTMVA.Add(*myZAI,1);
+ }
+ }
+
+ if( !fIsStepTime){
+ InputTree->Branch( "Time" ,&Time ,"Time/F" );
+ }
+ IsotopicVector IVAccordingToUserInfoFile = isotopicvector.GetThisComposition(IVInputTMVA);
+ double Ntot = IVAccordingToUserInfoFile.GetSumOfAll();
+ IVAccordingToUserInfoFile = IVAccordingToUserInfoFile/Ntot;
+
+ DBGV("INPUT TMVA");
+ for(unsigned int j = 0 ; j< fTMVAVariableNames.size() ; j++)
+ {
+ if(fTMVAFixedVariable[j]){
+ InputTMVA[j] = fTMVAFixedVariableValues[j];
+ }
+ else{
+ InputTMVA[j] = IVAccordingToUserInfoFile.GetZAIIsotopicQuantity( *myZAIs[j] ) ;
+ }
+ }
+ Time = fMLP_Time[TimeStep];
+
+ InputTree->Fill();
+
+ return InputTree;
+}
+//________________________________________________________________________
+double XSM_SFR::ExecuteTMVA(string WeightFile,TTree* InputTree)
+{
+ DBGV( "File :" << WeightFile);
+ // --- 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(unsigned int i = 0 ; i< fTMVAVariableNames.size() ; i++){
+ InputTMVA.push_back(0);
+ }
+ Float_t Time;
+
+ for(unsigned int j = 0 ; j< fTMVAVariableNames.size() ; j++){
+ reader->AddVariable( fTMVAVariableNames[j].c_str(),&InputTMVA[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;
+ for(unsigned int j = 0 ; j< fTMVAVariableNames.size() ; j++)
+ {
+ InputTree->SetBranchAddress(fTMVAVariableNames[j].c_str(),&InputTMVA[j]);
+ }
+ if(!fIsStepTime){
+ InputTree->SetBranchAddress( "Time" ,&Time );
+ }
+ InputTree->GetEntry(0);
+ Float_t val = (reader->EvaluateRegression( methodName ))[0];
+
+ delete reader;
+ DBGL
+
+ return (double)val;
+}
+//________________________________________________________________________
+EvolutionData XSM_SFR::GetCrossSectionsTime(IsotopicVector IV)
+{
+ DBGL
+
+ string dir = fTMVAWeightFolder;
+ if(dir[dir.size()-1] != '/'){
+ dir+= "/";
+ }
+
+ EvolutionData EvolutionDataFromMLP = EvolutionData();
+
+ map<ZAI,TGraph*> ExtrapolatedXS[3];
+ /*************DATA BASE INFO****************/
+ EvolutionDataFromMLP.SetReactorType(fDBRType);
+ EvolutionDataFromMLP.SetFuelType(fDBFType);
+ EvolutionDataFromMLP.SetPower(fDBPower);
+ EvolutionDataFromMLP.SetHeavyMetalMass(fDBHMMass);
+ /************* The Cross sections***********/
+ for(unsigned int i = 0;i<fWeightFiles.size();i++)
+ {
+ int Z = -2;
+ int A = -2;
+ int I = -2;
+ int Reaction = -2;
+ ReadWeightFile( fWeightFiles[i], Z, A, I, Reaction);
+ if( Z >= GetZAIThreshold() )
+ {
+ CLASSReader * reader = new CLASSReader( fTMVAVariableNames );
+ if(!fIsStepTime) { reader->AddVariable( "Time" ); }
+
+ reader->BookMVA( "MLP method" , dir + fWeightFiles[i] );
+
+ for(unsigned int TimeStep = 0;TimeStep<fMLP_Time.size();TimeStep++)
+ {
+ TTree* InputTree = CreateTMVAInputTree(IV,TimeStep);
+ reader->SetInputData( InputTree );
+ double XSValue = reader->EvaluateRegression( "MLP method" )[0];
+
+ pair< map<ZAI, TGraph*>::iterator, bool> IResult = ExtrapolatedXS[Reaction].insert( pair<ZAI ,TGraph* >(ZAI(Z,A,I), new TGraph()) );
+
+ if(IResult.second )
+ {
+ (IResult.first)->second->SetPoint(0, (double)fMLP_Time[TimeStep], XSValue );
+
+ }
+ else
+ {
+ (IResult.first)->second->SetPoint( (IResult.first)->second->GetN(), (double)fMLP_Time[TimeStep], XSValue );
+ }
+
+ delete InputTree;
+ }
+
+ delete reader;
+ }
+ }
+
+ /**********Sorting TGraph*********/
+ for(unsigned 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 NN per time step)
+//________________________________________________________________________
+//________________________________________________________________________
+void XSM_SFR::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_SFR::GetCrossSectionsStep(IsotopicVector IV)
+{
+ DBGL
+ TTree* InputTree = CreateTMVAInputTree(IV);
+ EvolutionData EvolutionDataFromMLP = EvolutionData();
+
+ map<ZAI,TGraph*> ExtrapolatedXS[3];
+ /*************DATA BASE INFO****************/
+ EvolutionDataFromMLP.SetReactorType(fDBRType);
+ EvolutionDataFromMLP.SetFuelType(fDBFType);
+ EvolutionDataFromMLP.SetPower(fDBPower);
+ EvolutionDataFromMLP.SetHeavyMetalMass(fDBHMMass);
+
+ /************* The Cross sections***********/
+ for(unsigned int i = 0;i<fWeightFiles.size();i++){ // JM2016 : need to booker the reader at each step because the file is different each time so we don't use CLASSReader
+ int Z = -2;
+ int A = -2;
+ int I = -2;
+ int Reaction = -2;
+ int TimeStep = -2;
+ ReadWeightFileStep( fWeightFiles[i], Z, A, I, Reaction, TimeStep);
+ if( Z >= GetZAIThreshold() ){
+ 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)fMLP_Time[TimeStep], ExecuteTMVA(fWeightFiles[i],InputTree) );
+ }
+ else{
+ (IResult.first)->second->SetPoint( (IResult.first)->second->GetN(), (double)fMLP_Time[TimeStep], ExecuteTMVA(fWeightFiles[i],InputTree) );
+ }
+ }
+ }
+
+ /**********Sorting TGraph*********/
+ for(unsigned 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_SFR::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;
+}
+//________________________________________________________________________
diff --git a/source/Model/XS/XSM_SFR.hxx b/source/Model/XS/XSM_SFR.hxx
new file mode 100644
index 0000000000000000000000000000000000000000..b4b05b78efb48501dd1ead0a89fbe3691ed73eca
--- /dev/null
+++ b/source/Model/XS/XSM_SFR.hxx
@@ -0,0 +1,143 @@
+#ifndef _XSM_SFR_HXX
+#define _XSM_SFR_HXX
+
+/*!
+ \file
+ \brief Header file for XSM_SFR class.
+
+
+ @authors Marc
+ @version 1.0
+ */
+#include <string>
+#include <fstream>
+#include <iostream>
+#include <map>
+#include <vector>
+
+#include "TTree.h"
+
+#include "XSModel.hxx"
+
+typedef long long int cSecond;
+using namespace std;
+
+
+class XSM_SFR;
+#ifndef __CINT__
+typedef void (XSM_SFR::*XS_SFR_DMthPtr)( const string & ) ;
+#endif
+//-----------------------------------------------------------------------------//
+//! Defines a XSModel getting mean cross sections from neural network execution
+
+/*!
+ Define a XSM_SFR.
+ This is the class to predict cross sections with a
+ set of Multi Layer Perceptrons (MLP)
+
+ @authors Marc
+ @version 1.0
+ */
+//________________________________________________________________________
+
+
+class XSM_SFR : public XSModel
+{
+ public :
+
+ /*!
+ \name Constructor/Desctructor
+ */
+ //@{
+
+ //{
+ /// Normal Constructor
+ /*!
+ \param TMVA_Weight_Directory : The directory where all the TMVA weight are located
+ \param InformationFile : Name of the information file located in TMVA_Weight_Directory (default : Data_Base_Info.nfo)
+ \param IsTimeStep : if true , one TMVA weihgt per step time is requiered otherwise it assumes time is part of the MLP inputs
+
+ */
+ XSM_SFR(string TMVA_Weight_Directory,map<string,double> FixedParameters,string InformationFile = "/Data_Base_Info.nfo",bool IsTimeStep = false);
+ //}
+
+ //{
+ /// CLASSLogger Constructor
+ /*!
+ \param log : The CLASSLogger
+ \param TMVA_Weight_Directory : The directory where all the TMVA weight are located
+ \param InformationFile : Name of the information file located in TMVA_Weight_Directory (default : Data_Base_Info.nfo)
+ \param IsTimeStep : if true , one TMVA weihgt per step time is requiered otherwise it assumes time is part of the MLP inputs
+
+ */
+ XSM_SFR(CLASSLogger* Log,string TMVA_Weight_Directory,map<string,double> FixedParameters,string InformationFile = "/Data_Base_Info.nfo",bool IsTimeStep = false);
+ //}
+
+ ~XSM_SFR();
+ //@}
+
+ /*!
+ \name Reading NFO related Method
+ */
+ //@{
+
+ //{
+ /// LoadKeyword() : make the correspondance between keyword and reading method
+ void LoadKeyword();
+ //}
+
+ //{
+ /// ReadTimeSteps : read the time step of the model
+ /*!
+ \param line : line suppossed to contain the time step information starts with "k_timestep" keyword
+ */
+ void ReadTimeSteps(const string &line);
+ //}
+
+ //{
+ /// ReadZAIName : read the zai name in the TMWA MLP model
+ /*!
+ \param line : line suppossed to contain the ZAI name starts with "k_zainame" keyword
+ */
+ void ReadInputParameter(const string &line);
+ //}
+ //{
+ /// ReadLine : read a line
+ /*!
+ \param line : line to read
+ */
+ void ReadLine(string line);
+ //}
+
+ //@}
+ void FixTMVAVariable(string VariableName,double VariableValue);
+
+ void SetFixedVariablesValues(map<string,double> FixedParameters);
+
+ EvolutionData GetCrossSections(IsotopicVector IV,double t = 0); //!< Return calculated cross section by the MLP regression
+
+ private :
+
+ void GetMLPWeightFiles(); //!< Find all .xml file in TMVA_Weight_Directory
+ EvolutionData GetCrossSectionsStep(IsotopicVector IV); //!< Return calculated cross section by the MLP regression when fIsTimeStep == true
+ EvolutionData GetCrossSectionsTime(IsotopicVector IV); //!< Return calculated cross section by the MLP regression when fIsTimeStep == false
+ void ReadWeightFile(string Filename, int &Z, int &A, int &I, int &Reaction) ; //!< Select the reaction according to the weight file name
+ void ReadWeightFileStep(string Filename, int &Z, int &A, int &I, int &Reaction, int &TimeStep);; //!< Select the reaction according to the weight file name
+ double ExecuteTMVA(string WeightFile, TTree* InputTree); //!<Execute the MLP according to the input tree created
+ TTree* CreateTMVAInputTree(IsotopicVector isotopicvector,int TimeStep = 0); //!<Create input tmva tree to be read by ExecuteTMVA
+
+ vector<double> fMLP_Time; //!< Time vector of the data base
+ vector<string> fWeightFiles; //!< All the weight file contains in fTMVAWeightFolder
+ string fTMVAWeightFolder; //!< folder containing all the weight file
+ bool fIsStepTime; //!< true if one TMVA weihgt per step time is requiered otherwise it assumes time is part of the MLP inputs
+ vector<string> fTMVAVariableNames;//!< List of TMVA input variable names (read from fMLPInformationFile ) , name depends on the training step
+ vector<bool> fTMVAFixedVariable;//!< List of value for TMVA that have to be used for all
+ vector<double> fTMVAFixedVariableValues;//!< List of value for TMVA that have to be used for all
+
+#ifndef __CINT__
+ map<string, XS_SFR_DMthPtr> fDKeyword;
+#endif
+};
+
+#endif
+