EQ_OneParameter.cxx

#include "EquivalenceModel.hxx"
#include "EQ_OneParameter.hxx"
#include "StringLine.hxx"
#include "CLASSMethod.hxx"

#include <string>
#include <iostream>
#include <fstream>
#include <algorithm>
#include <cmath>
#include <cassert>

#include "TSystem.h"
#include "TMVA/Reader.h"
#include "TMVA/Tools.h"
#include "TMVA/MethodCuts.h"

#include "CLASSReader.hxx"

//________________________________________________________________________
//________________________________________________________________________
EQ_OneParameter::EQ_OneParameter(string TMVAXMLFilePath, string TMVANFOFilePath):EquivalenceModel(new CLASSLogger("XSM_MLP.log"))
{
    fUseTMVAPredictor = true;

    fMaxIterration  = 500;
    fPCMprecision  = 10;

    fTMVAXMLFilePath = TMVAXMLFilePath;
    fTMVANFOFilePath = TMVANFOFilePath;

    fDBRType = "";
    fDBFType = "";
    fSpecificPower = 0;
    fMaximalBU = 0;
    fTargetParameter = "";
    fTargetParameterStDev = 0;
    fBuffer = "";
    fPredictorType = "";
    fOutput = "";

    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);}

    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();

}
//________________________________________________________________________
EQ_OneParameter::EQ_OneParameter(CLASSLogger* log, string TMVAXMLFilePath, string TMVANFOFilePath):EquivalenceModel(log)
{
    fUseTMVAPredictor = true;

    fMaxIterration  = 500;
    freaded                = false;     
    fPCMprecision         = 10;

    fTMVAXMLFilePath = TMVAXMLFilePath;
    fTMVANFOFilePath = TMVANFOFilePath;

    fDBRType = "";
    fDBFType = "";
    fSpecificPower = 0;
    fMaximalBU = 0;
    fTargetParameter = "";
    fTargetParameterStDev = 0;
    fBuffer = "";
    fPredictorType = "";
    fOutput = "";

    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);}

    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();}
//________________________________________________________________________
EQ_OneParameter::EQ_OneParameter(string TMVANFOFilePath):EquivalenceModel(new CLASSLogger("XSM_MLP.log"))
{
    fUseTMVAPredictor = false;

    fTMVANFOFilePath = TMVANFOFilePath;

    fDBRType = "";
    fDBFType = "";
    fSpecificPower = 0;
    fMaximalBU = 0;
    fTargetParameter = "";
    fTargetParameterStDev = 0;
    fBuffer = "";

    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);}

    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)
{
    fUseTMVAPredictor = false;

    fTMVANFOFilePath = TMVANFOFilePath;

    fDBRType = "";
    fDBFType = "";
    fSpecificPower = 0;
    fMaximalBU = 0;
    fBuffer = "";

    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);}

    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();}
//________________________________________________________________________
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)
{
    DBGL
    //Iterators declaration
    map < string , vector  <IsotopicVector> >::const_iterator it_s_vIV;
    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++)
    {
        if((*it_s_B ).second==true){ BufferDenomination = (*it_s_B).first; }    
    }

    for(int i = 0; i< lambda[BufferDenomination].size(); i++)
    {
        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++)
        {
            IV  +=  lambda[(*it_s_vIV).first][i] * StreamArray.at( it_s_vIV->first )[i];    
        }
    }

    //Calculate MassBuffer
    double MassBuffer = HMMass - IV.GetTotalMass()*1e06;

    //Set buffer lambda according to MassBuffer
    ConvertMassToLambdaVector(BufferDenomination, lambda[BufferDenomination], MassBuffer, StreamArray.at(BufferDenomination));

    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++)
        {
            IV  +=  lambda[(*it_s_vIV).first][i] * StreamArray.at( it_s_vIV->first )[i];    
        }
    }
    DBGL
    return IV; 

}

//________________________________________________________________________
map <string , vector<double> > EQ_OneParameter::BuildFuel(double BurnUp, double HMMass, map < string , vector <IsotopicVector> > StreamArray,  map < string , double> StreamListFPMassFractionMin, map < string , double> StreamListFPMassFractionMax, map < int , string > StreamListPriority, map < string , bool> StreamListIsBuffer)
{
    DBGL

    HMMass *=  1e6; //Unit conversion : tons to gram

    map <string , vector<double> > lambda ; // map containing name of the list and associated vector of proportions taken from stocks
    //Iterators declaration
    map < string , vector  <IsotopicVector> >::iterator it_s_vIV;
    map < string , vector  <double> >::iterator it_s_vD;
    map < string , IsotopicVector >::iterator it_s_IV;
    map < string , double >::iterator it_s_D;
    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++)
        {
            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 = false; 
    for( it_s_vIV = StreamArray.begin();  it_s_vIV != StreamArray.end(); it_s_vIV++)
    {
        if(StreamArray[(*it_s_vIV).first].size() == 0)
        {
            WARNING << " No stock available for stream : "<< (*it_s_vIV).first <<".  Fuel not built." << endl;
            SetLambdaToErrorCode(lambda[(*it_s_vIV).first]);
            BreakReturnLambda = true;   
        }
    }
    if(BreakReturnLambda) { return lambda;} 

    // Check if the targeted burn-up is lower than maximum burn-up of model //
    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);    
    }

// Fissile fraction calculation is needed.
    if (fUseTMVAPredictor)
    {
        // Check if EquivalenceModel->SetTMVAXMLFilePath() and/or EquivalenceModel->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 EquivalenceModel->SetTMVAXMLFilePath() and/or EquivalenceModel->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)
            {
                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"];}
        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);            
        }
            
        /// 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
            {
                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
            {
                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++)
        {
            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;
                SetLambdaToErrorCode(lambda[(*it_s_D).first]);
                BreakReturnLambda = true;   
            }
        }
        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++)
        {
            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; 
    
        IsotopicVector FuelToTest;
    
        bool TargetParameterIncluded = false;
        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();
            TargetParameterMin[(*it_i_s ).second] =  CalculateTargetParameter(FuelToTest, fTargetParameter);
    
            //Check is TargetParameterMin < TargetParameter
            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
                }
                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;
                    it_i_s --;
                    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;
                }
            }
            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];    
            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();
            TargetParameterMax[(*it_i_s ).second]   =  CalculateTargetParameter(FuelToTest, fTargetParameter);
        
            if(TargetParameterMax[(*it_i_s ).second]>=TargetParameterValue)
            {
                TargetParameterIncluded = true ; 
                break;
            }
        }
    
        //Check if target parameter increases monotously with the material mass
        CheckTargetParameterConsistency(StreamListPriority, TargetParameterMin, TargetParameterMax);
    
        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;
        }
    
        //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    
    
        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;
                }
                */
    
        do
        {
            if(count > fMaxIterration)
            {
                ERROR << "CRITICAL ! Can't manage to predict fissile content\nHint : Try to decrease the precision on the target parameter using :\nYourEquivalenceModel->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; 
                exit(1);
            }
    
            if( (CalculatedTargetParameter - TargetParameterValue) < 0 ) //Need to add more fissile material in fuel
            {
                LastMassMinus = MassToAdd;
                MassToAdd   = MassToAdd + fabs(LastMassPlus - MassToAdd)/2.;
            }
            else if( (CalculatedTargetParameter - TargetParameterValue) > 0) //Need to add less fissile material in fuel
            {
                LastMassPlus    = MassToAdd;
                MassToAdd   = MassToAdd - fabs(LastMassMinus - MassToAdd)/2.;
            }
            ConvertMassToLambdaVector(MaterialToSearch, lambda[MaterialToSearch], MassToAdd, StreamArray[MaterialToSearch]);    
            FuelToTest          = BuildFuelToTest(lambda, StreamArray, HMMass, StreamListIsBuffer);
            FuelToTest          = FuelToTest/FuelToTest.GetSumOfAll();
            CalculatedTargetParameter   = CalculateTargetParameter(FuelToTest, fTargetParameter);
            
            count ++;
    
        }while(fabs(TargetParameterValue - CalculatedTargetParameter) > GetTargetParameterStDev()*TargetParameterValue);    
    }
// Fissile fraction is imposed by the FP
// No need to use algo
// Simplified fuel building
    else
    {
        // Check if EquivalenceModel->SetTMVANFOFilePath() have been defined
        if (fTMVANFOFilePath.empty())
        {
            ERROR << " TMVA NFO File path is not defined..."<< endl;
            ERROR << " You have to use EquivalenceModel->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
            {
                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
            {
                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++)
        {
            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;
                SetLambdaToErrorCode(lambda[(*it_s_D).first]);
                BreakReturnLambda = true;   
            }
        }
        if(BreakReturnLambda) { return lambda;}

        IsotopicVector FuelToTest;
        // Build Fuel
        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]);
            FuelToTest = BuildFuelToTest(lambda, StreamArray, HMMass, StreamListIsBuffer);
        }
    }

    //Final builded fuel 
    IsotopicVector IVStream;
    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++) 
        {
            IVStream +=lambda[(*it_s_vD).first][i] * StreamArray[(*it_s_vD).first][i];
        }
    }   
    
    //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(int i=0; i < (int)lambda[(*it_s_vD).first].size(); i++)
        {
            DBGV(lambda[(*it_s_vD).first][i]); 
        }
    }
    
    DBGL
    return lambda;
}

//________________________________________________________________________
TTree* EQ_OneParameter::CreateTMVAInputTree(IsotopicVector TheFreshfuel, double ThisTime)
{
    /******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++)
        InputTMVA.push_back(0);
    
    float Time = 0;
    
    IsotopicVector IVInputTMVA;
    map<ZAI ,string >::iterator it_ZAI_s;
    int j = 0;
    
    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;
        j++;
    }
    
    if(ThisTime != -1)
        InputTree->Branch("Time" ,&Time ,"Time/F");
    
    IsotopicVector IVAccordingToUserInfoFile    = TheFreshfuel.GetThisComposition(IVInputTMVA);
    double Ntot                     = IVAccordingToUserInfoFile.GetSumOfAll();
    IVAccordingToUserInfoFile           = IVAccordingToUserInfoFile/Ntot;
    
    j = 0;
    
    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)
{
    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++)
    {   
        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())
        {
             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;
            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;
            exit(1);
        }
    }
}
//________________________________________________________________________
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);
    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);            
    }

    return ParameterToCalculate ;
}
//________________________________________________________________________
double EQ_OneParameter::CalculateBurnUpMax(IsotopicVector TheFuel, map<string, double> ModelParameter)
{
    /**************************************************************************/
    //With a dichotomy, the maximal irradiation time (TheFinalTime) is calculated
    //When average Kinf is very close (according "Precision") to the threshold
    //then the corresponding irradiation time is convert in burnup and returned
    /**************************************************************************/
    //Algorithm initialization
    double KThreshold       = fModelParameter["kThreshold"];
    int NumberOfBatch       = (int)fModelParameter["NumberOfBatch"];
    double OldFinalTimeMinus    = 0;
    double MaximumBU        = fMaximalBU;
    double MinimumBU        = 0 ;
    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++)
    {
        float TheTime = (b+1)*TheFinalTime/NumberOfBatch;

        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 )
        {
            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
        {
            OldFinalTimeMinus = TheFinalTime;
            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
        {
            OldFinalTimePlus = TheFinalTime;
            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++)
        {
            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;
        //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() )  ;
    
    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;
} 
//________________________________________________________________________
void EQ_OneParameter::ReadNFO()
{
    DBGL
    ifstream NFO(fTMVANFOFilePath.c_str());
    
    if(!NFO)
    {
        ERROR << "Can't find/open file " << fTMVANFOFilePath << endl;
        exit(0);
    }
    
    do
    {
        string line;
        getline(NFO,line);
        
        EquivalenceModel::ReadLine(line);
        
    } 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, EQM_MthPtr>::iterator it = fKeyword.find(keyword);
        
        if(it != fKeyword.end())
            (this->*(it->second))( line );
        
        freaded = true;
        ReadLine(line);
        
    }
    
    freaded = false;
    
    DBGL
}
//________________________________________________________________________
void EQ_OneParameter::LoadKeyword() 
{
    DBGL
    fKeyword.insert( pair<string, EQM_MthPtr>( "k_zail",                & EquivalenceModel::ReadZAIlimits)           );
    fKeyword.insert( pair<string, EQM_MthPtr>( "k_reactor",         & EquivalenceModel::ReadType)            );
    fKeyword.insert( pair<string, EQM_MthPtr>( "k_fuel",                & EquivalenceModel::ReadType)            );
    fKeyword.insert( pair<string, EQM_MthPtr>( "k_massfractionmin",     & EquivalenceModel::ReadEqMinFraction)       );
    fKeyword.insert( pair<string, EQM_MthPtr>( "k_massfractionmax",     & EquivalenceModel::ReadEqMaxFraction)       );
    fKeyword.insert( pair<string, EQM_MthPtr>( "k_list",                & EquivalenceModel::ReadList)            );
    fKeyword.insert( pair<string, EQM_MthPtr>( "k_specpower",           & EquivalenceModel::ReadSpecificPower)       );
    if (fUseTMVAPredictor) fKeyword.insert( pair<string, EQM_MthPtr>( "k_zainame",          & EquivalenceModel::ReadZAIName)             );
    fKeyword.insert( pair<string, EQM_MthPtr>( "k_maxburnup",           & EquivalenceModel::ReadMaxBurnUp)       ); 
    if (fUseTMVAPredictor) fKeyword.insert( pair<string, EQM_MthPtr>( "k_targetparameter",      & EquivalenceModel::ReadTargetParameter)         );
    if (fUseTMVAPredictor) fKeyword.insert( pair<string, EQM_MthPtr>( "k_predictortype",        & EquivalenceModel::ReadPredictorType)       );
    if (fUseTMVAPredictor) fKeyword.insert( pair<string, EQM_MthPtr>( "k_output",           & EquivalenceModel::ReadOutput)              );
    fKeyword.insert( pair<string, EQM_MthPtr>( "k_buffer",          & EquivalenceModel::ReadBuffer)          ); 
    if (fUseTMVAPredictor) fKeyword.insert( pair<string, EQM_MthPtr>( "k_modelparameter",       & EquivalenceModel::ReadModelParameter)          ); 
    if (fUseTMVAPredictor) fKeyword.insert( pair<string, EQM_MthPtr>( "k_targetparameterstdev",     & EquivalenceModel::ReadTargetParameterStDev)    ); 

    DBGL
}
//________________________________________________________________________
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
    {
        ERROR << " Bad keyword : " << keyword << " Not found !" << endl;
        exit(1);
    }
    if( keyword ==  "k_fuel" )
        fDBFType = StringLine::NextWord(line, pos, ' ');
    else if( keyword ==  "k_reactor" )
        fDBRType = StringLine::NextWord(line, pos, ' ');
    
    DBGL
}
//________________________________________________________________________
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
    {
        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)));
    DBGL
}
//________________________________________________________________________
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
    {
        ERROR << " Bad keyword : \"k_list\" not found !" << endl;
        exit(1);
    }
    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
}
//________________________________________________________________________
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
    {
        ERROR << " Bad keyword : \"k_massfractionmin\" not found !" << endl;
        exit(1);
    }
    string ListName= StringLine::NextWord(line, pos, ' ');
    double Q     = atof(StringLine::NextWord(line, pos, ' ').c_str());
    fStreamListEqMMassFractionMin[ListName] = Q;

    DBGL
}

//________________________________________________________________________
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
    {
        ERROR << " Bad keyword : \"k_massfractionmax\" not found !" << endl;
        exit(1);
    }
    string ListName= StringLine::NextWord(line, pos, ' ');
    double Q     = atof(StringLine::NextWord(line, pos, ' ').c_str());
    fStreamListEqMMassFractionMax[ListName] = Q;

    DBGL
}

//________________________________________________________________________
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
    {
        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
    {