This is the next CLASS release (4.2). Its the Fanny's(FaC) version of CLASS...

This is the next CLASS release (4.2). Its the Fanny's(FaC) version of CLASS with the multi stream FP. This version should replace the Trunk version soon (once its completed, debbuged and tested

git-svn-id: svn+ssh://svn.in2p3.fr/class@726 0e7d625b-0364-4367-a6be-d5be4a48d228

source/branches/CLASSV4.2/Model/Equivalence/EQM_FBR_BakerRoss_MOX.cxx

+#include "EQM_FBR_BakerRoss_MOX.hxx"
+#include "CLASSLogger.hxx"
+
+#include <vector>
+
+//________________________________________________________________________
+//
+//		EQM_FBR_BakerRoss_MOX
+//
+//	Equivalenve Model based on Plutonium 239 equivalent  
+//	using the formula of Baker&Ross
+//
+//________________________________________________________________________
+
+
+EQM_FBR_BakerRoss_MOX::EQM_FBR_BakerRoss_MOX(double Weight_U_235, double Weight_Pu_238, double Weight_Pu_240, double Weight_Pu_241, double Weight_Pu_242, double Weight_Am_241, double EquivalentFissile):EquivalenceModel(new CLASSLogger("EQM_FBR_BakerRoss_MOX.log"))
+{
+	ZAI U8(92,238,0);
+	ZAI U5(92,235,0);
+	double U5_enrich = 0.0025;
+
+	ZAI Pu8(94,238,0);
+	ZAI Pu9(94,239,0);
+	ZAI Pu0(94,240,0);
+	ZAI Pu1(94,241,0);
+	ZAI Pu2(94,242,0);
+
+	FertileList = U5*U5_enrich + U8*(1-U5_enrich); //Default Fertile composition (if no Fertile Storage is set for the FabricationPlant  )	
+	FissileList = Pu8*1+Pu9*1+Pu0*1+Pu1*1+Pu2*1; //ZAI to be extracted by the FabricationPlant in its affiliated Fissile Storage 
+	fStreamList["Fissile"] = FissileList;
+	fStreamList["Fertile"] = FertileList;
+
+	fReferenceFissilContent = EquivalentFissile;//Equivalent fissile content
+
+	fFirstGuessContent["Fissile"] = fReferenceFissilContent;
+	fFirstGuessContent["Fertile"] = 1 - fFirstGuessContent["Fissile"];
+
+	fWeight_U_235  = Weight_U_235 ;
+	fWeight_Pu_238 = Weight_Pu_238;
+	fWeight_Pu_240 = Weight_Pu_240;
+	fWeight_Pu_241 = Weight_Pu_241;
+	fWeight_Pu_242 = Weight_Pu_242;
+	fWeight_Am_241 = Weight_Am_241;
+
+
+	INFO << "__An equivalence model of FBR MOX has been defined__" << endl;
+	INFO << "\tThis model is based on Plutonium 239 equivalent" << endl;
+	INFO << "\t\t Weigt values : " << endl;
+	INFO << "\t\t Fertile :  " << endl;
+	INFO << "\t\t\tU_235 " <<   Weight_U_235 << endl;
+	INFO << "\t\t\tU_238 0 (by definition)" << endl;
+	INFO << "\t\t Fissile :  " << endl;
+	INFO << "\t\t\tPu_238 " << Weight_Pu_238 << endl;
+	INFO << "\t\t\tPu_239 1 (by definition)" << endl;
+	INFO << "\t\t\tPu_240 " << Weight_Pu_240 << endl;
+	INFO << "\t\t\tPu_241 " << Weight_Pu_241 << endl;
+	INFO << "\t\t\tPu_242 " << Weight_Pu_242 << endl;
+	INFO << "\t\t\tAm_241 " << Weight_Am_241 << endl;
+
+}
+//________________________________________________________________________
+EQM_FBR_BakerRoss_MOX::EQM_FBR_BakerRoss_MOX(CLASSLogger* log, double Weight_U_235, double Weight_Pu_238, double Weight_Pu_240, double Weight_Pu_241, double Weight_Pu_242, double Weight_Am_241, double EquivalentFissile ):EquivalenceModel(log)
+{
+	ZAI U8(92,238,0);
+	ZAI U5(92,235,0);
+	double U5_enrich = 0.0025;
+	
+	ZAI Pu8(94,238,0);
+	ZAI Pu9(94,239,0);
+	ZAI Pu0(94,240,0);
+	ZAI Pu1(94,241,0);
+	ZAI Pu2(94,242,0);
+
+	FertileList = U5*U5_enrich + U8*(1-U5_enrich); 	//Default Fertile composition (if no Fertile Storage is set for the FabricationPlant  )	
+	FissileList = Pu8*1+Pu9*1+Pu0*1+Pu1*1+Pu2*1; 	//ZAI to be extracted by the FabricationPlant in its affiliated Fissile Storage 
+	fStreamList["Fissile"] = FissileList;
+	fStreamList["Fertile"] = FertileList;
+
+	fReferenceFissilContent = EquivalentFissile;		//Equivalent fissile content
+
+	fFirstGuessContent["Fissile"] = fReferenceFissilContent;
+	fFirstGuessContent["Fertile"] = 1 - fFirstGuessContent["Fissile"];
+
+	fWeight_U_235  = Weight_U_235 ;
+	fWeight_Pu_238 = Weight_Pu_238;
+	fWeight_Pu_240 = Weight_Pu_240;
+	fWeight_Pu_241 = Weight_Pu_241;
+	fWeight_Pu_242 = Weight_Pu_242;
+	fWeight_Am_241 = Weight_Am_241;
+
+
+	INFO << "__An equivalence model of FBR MOX has been defined__" << endl;
+	INFO << "\tThis model is based on Plutonium 239 equivalent" << endl;
+	INFO << "\t\t Weigt values : " << endl;
+	INFO << "\t\t Fertile : " << endl;
+	INFO << "\t\t\tU_235 " <<   Weight_U_235 << endl;
+	INFO << "\t\t\tU_238 0 (by definition)" << endl;
+	INFO << "\t\t Fissile : " << endl;
+	INFO << "\t\t\tPu_238 " << Weight_Pu_238 << endl;
+	INFO << "\t\t\tPu_239 1 (by definition)" << endl;
+	INFO << "\t\t\tPu_240 " << Weight_Pu_240 << endl;
+	INFO << "\t\t\tPu_241 " << Weight_Pu_241 << endl;
+	INFO << "\t\t\tPu_242 " << Weight_Pu_242 << endl;
+	INFO << "\t\t\tAm_241 " << Weight_Am_241 << endl;
+
+
+}
+//________________________________________________________________________
+map < string , double> EQM_FBR_BakerRoss_MOX::GetMolarFraction(map < string , IsotopicVector> IVStream, double BurnUp)
+{
+
+	if( BurnUp != 0 )
+		WARNING << "Burn up (third argument) has no effect " << endl;
+	
+	IsotopicVector Fissile = IVStream["Fissile"];
+	IsotopicVector Fertile = IVStream["Fertile"];
+
+	double FissileContent = 0.;       
+
+	IsotopicVector FissileListPlusDecay;
+	FissileListPlusDecay.Add(94,238,0,1);
+	FissileListPlusDecay.Add(94,239,0,1);
+	FissileListPlusDecay.Add(94,240,0,1);
+	FissileListPlusDecay.Add(94,241,0,1);
+	FissileListPlusDecay.Add(94,242,0,1);
+	FissileListPlusDecay.Add(95,241,0,1);
+
+	IsotopicVector FertileList;
+	FertileList.Add(92,238,0,1);
+	FertileList.Add(92,239,0,1);
+
+	//Getting the fissile from the Fissile input & normalize it
+	IsotopicVector FissileFromInput = Fissile.GetThisComposition(FissileListPlusDecay);
+	FissileFromInput = FissileFromInput/FissileFromInput.GetSumOfAll();
+
+	//Getting the fissile from the Fissile input & normalize it
+	IsotopicVector FertileFromInput = Fertile.GetThisComposition(FertileList);
+	FertileFromInput = FertileFromInput/FertileFromInput.GetSumOfAll();
+
+	double  SumWeightNFissile = fWeight_Pu_238 * FissileFromInput.GetZAIIsotopicQuantity(94,238,0) 
+								+ 1		 * FissileFromInput.GetZAIIsotopicQuantity(94,239,0)
+								+ fWeight_Pu_240 * FissileFromInput.GetZAIIsotopicQuantity(94,240,0) 
+								+ fWeight_Pu_241 * FissileFromInput.GetZAIIsotopicQuantity(94,241,0) 
+								+ fWeight_Pu_242 * FissileFromInput.GetZAIIsotopicQuantity(94,242,0) 
+								+ fWeight_Am_241 * FissileFromInput.GetZAIIsotopicQuantity(95,241,0) ;
+
+	double  SumWeightNFertile = fWeight_U_235 * FertileFromInput.GetZAIIsotopicQuantity(92,235,0);
+
+	FissileContent = (fReferenceFissilContent - SumWeightNFertile )/(SumWeightNFissile-SumWeightNFertile); //Baker & Ross formula
+
+	map < string , double> MolarFraction;
+	
+	MolarFraction["Fissile"] = FissileContent;
+	MolarFraction["Fertile"] = 1.- FissileContent;
+
+	return MolarFraction; //return Molar content of each component in the fuel
+}

source/branches/CLASSV4.2/Model/Equivalence/EQM_FBR_BakerRoss_MOX.hxx

+#ifndef _EQM_FBR_BakerRoss_MOX_HXX
+#define _EQM_FBR_BakerRoss_MOX_HXX
+
+#include "EquivalenceModel.hxx"
+
+#include <string>
+
+using namespace std;
+
+//-----------------------------------------------------------------------------//
+//! Defines an EquivalenceModel based on Baker and Ross formula
+
+/*!
+ The aim of these class is to constuct a fuel from an equivalence model
+ based on a @f$^{239}Pu@f$ equivalent from Baker and Ross formula :
+ 
+ It returns Pu content (E) needed for the FBR-Na loaded with a given Pu vector according
+  :
+ 
+ @f$E = \frac{E_{ref} - \sum_{fertile}N_{i}W_{i} }{\sum_{fissile}N_{i}W_{i}-\sum_{fertile}N_{i}W_{i}}@f$
+ 
+ with :
+ 
+ @f$W_i = \frac{\alpha_{i} - \alpha_{^{238}U} }{\alpha_{^{239}Pu}-\alpha_{^{238}U}}@f$
+ and
+ @f$\alpha_{i} = \bar{\nu_{i}}\cdot\sigma_{i}^{fis} - \sigma_{i}^{abs}@f$
+
+ @author BLG
+ @version 3.0
+ */
+//________________________________________________________________________
+
+class EQM_FBR_BakerRoss_MOX : public EquivalenceModel
+{
+	public :
+	
+	/*!
+	 \name Constructor
+	 */
+	//@{
+	//{
+	/// Logger constructor
+	
+	/*!
+	 Make a new EQM_FBR_BakerRoss_MOX : the default values have been calculated for a FBR-Na of type : ESFR like (without blanket)
+	 \param log : use for the log
+	 \param Weight_U_235  : reactivity weight @f$W_{^{235}U} = \frac{\alpha_{{^{235}U}} - \alpha_{^{238}U} }{\alpha_{^{239}Pu}-\alpha_{^{238}U}}@f$
+	 \param Weight_Pu_238 : reactivity weight @f$W_{^{238}Pu} = \frac{\alpha_{{^{238}Pu}} - \alpha_{^{238}U} }{\alpha_{^{239}Pu}-\alpha_{^{238}U}}@f$
+	 \param Weight_Pu_240 : reactivity weight @f$W_{^{240}Pu} = \frac{\alpha_{{^{240}Pu}} - \alpha_{^{238}U} }{\alpha_{^{239}Pu}-\alpha_{^{238}U}}@f$
+	 \param Weight_Pu_241 : reactivity weight @f$W_{^{241}Pu} = \frac{\alpha_{{^{241}Pu}} - \alpha_{^{238}U} }{\alpha_{^{239}Pu}-\alpha_{^{238}U}}@f$
+	 \param Weight_Pu_242 : reactivity weight @f$W_{^{242}Pu} = \frac{\alpha_{{^{242}Pu}} - \alpha_{^{238}U} }{\alpha_{^{239}Pu}-\alpha_{^{238}U}}@f$
+	 \param Weight_Am_241 : reactivity weight @f$W_{^{241}Pu} = \frac{\alpha_{{^{241}Am}} - \alpha_{^{238}U} }{\alpha_{^{239}Pu}-\alpha_{^{238}U}}@f$
+	 \param EquivalentFissile : reference fresh fuel @f$^{239}Pu@f$ content neeed in a @f$^{238}U@f$ + @f$^{239}Pu@f$ fuel to satisfy criticality at t = 0
+	 */
+	EQM_FBR_BakerRoss_MOX(CLASSLogger* log, double Weight_U_235 = 0.791135, double Weight_Pu_238 = 0.686385, double Weight_Pu_240 = 0.13553, double Weight_Pu_241 = 1.54572 , double Weight_Pu_242 = 0.0829001, double Weight_Am_241 = -0.336945, double EquivalentFissile = 0.103213);
+	//}
+	
+	//{
+	/// normal constructor
+	
+	/*!
+	 Make a new EQM_FBR_BakerRoss_MOX : the default values have been calculated for a FBR-Na of type : ESFR like (without blanket)
+	 \param Weight_U_235  : reactivity weight @f$W_{^{235}U} = \frac{\alpha_{{^{235}U}} - \alpha_{^{238}U} }{\alpha_{^{239}Pu}-\alpha_{^{238}U}}@f$
+	 \param Weight_Pu_238 : reactivity weight @f$W_{^{238}Pu} = \frac{\alpha_{{^{238}Pu}} - \alpha_{^{238}U} }{\alpha_{^{239}Pu}-\alpha_{^{238}U}}@f$
+	 \param Weight_Pu_240 : reactivity weight @f$W_{^{240}Pu} = \frac{\alpha_{{^{240}Pu}} - \alpha_{^{238}U} }{\alpha_{^{239}Pu}-\alpha_{^{238}U}}@f$
+	 \param Weight_Pu_241 : reactivity weight @f$W_{^{241}Pu} = \frac{\alpha_{{^{241}Pu}} - \alpha_{^{238}U} }{\alpha_{^{239}Pu}-\alpha_{^{238}U}}@f$
+	 \param Weight_Pu_242 : reactivity weight @f$W_{^{242}Pu} = \frac{\alpha_{{^{242}Pu}} - \alpha_{^{238}U} }{\alpha_{^{239}Pu}-\alpha_{^{238}U}}@f$
+	 \param Weight_Am_241 : reactivity weight @f$W_{^{241}Pu} = \frac{\alpha_{{^{241}Am}} - \alpha_{^{238}U} }{\alpha_{^{239}Pu}-\alpha_{^{238}U}}@f$
+	 \param EquivalentFissile : reference fresh fuel @f$^{239}Pu@f$ content neeed in a @f$^{238}U@f$ + @f$^{239}Pu@f$ fuel to satisfy criticality at t = 0
+	 */
+	EQM_FBR_BakerRoss_MOX(double Weight_U_235 = 0.791135, double Weight_Pu_238 = 0.686385, double Weight_Pu_240 = 0.13553, double Weight_Pu_241 = 1.54572 , double Weight_Pu_242 = 0.0829001, double Weight_Am_241 = -0.336945, double EquivalentFissile = 0.103213);
+	//}
+	//@}
+	
+	virtual map < string , double> GetMolarFraction(map < string , IsotopicVector> IVStream, double BurnUp=0)
+
+	private :
+	double fReferenceFissilContent;   //!<The reference fraction of Pu for BurnUp = 100Gwj/t with a ideal model(only @f$^{239}Pu@f$ and &@f$^{238}U@f$ are taken into account)
+
+	/*!
+	 \name Reactivity coefficients :
+	 */
+	//@{
+	double fWeight_U_235;	//!< weight for @f$^{235}U@f$
+	double fWeight_Pu_238;	//!< weight for @f$^{238}Pu@f$
+	double fWeight_Pu_240;	//!< weight for @f$^{240}Pu@f$
+	double fWeight_Pu_241;	//!< weight for @f$^{241}Pu@f$
+	double fWeight_Pu_242;	//!< weight for @f$^{242}Pu@f$
+	double fWeight_Am_241;	//!< weight for @f$^{241}Am@f$
+	//@}
+};
+
+#endif
+
+

source/branches/CLASSV4.2/Model/Equivalence/EQM_FBR_MLP_Keff.cxx

+#include "EQM_FBR_MLP_Keff.hxx"
+#include "CLASSLogger.hxx"
+#include "CLASSMethod.hxx"
+#include "StringLine.hxx"
+
+#include <string>
+#include <iostream>
+#include <fstream>
+#include <algorithm>
+#include <cmath>
+#include <cassert>
+#include <map>
+
+#include "TSystem.h"
+#include "TMVA/Reader.h"
+#include "TMVA/Tools.h"
+#include "TMVA/MethodCuts.h"
+
+//________________________________________________________________________
+//
+//		EQM_FBR_MLP_Keff
+//
+//	Equivalenve Model based on multi layer perceptron from TMVA (root cern)
+//	For FBR 
+//
+//________________________________________________________________________
+
+
+//________________________________________________________________________
+//
+//	Objects & Methods for content prediction with calculation of Keff at
+//	one given time (often BOC or EOC)
+//
+//________________________________________________________________________
+
+//________________________________________________________________________
+EQM_FBR_MLP_Keff::EQM_FBR_MLP_Keff(string TMVAWeightPath, double keff_target, string InformationFile):EquivalenceModel(new CLASSLogger("EQM_FBR_MLP_Keff.log"))
+{
+	DBGL
+
+	/**The tmva weight*/
+	fTMVAWeightPath = TMVAWeightPath;
+
+
+	/* INFORMATION FILE HANDLING */
+	
+	if(InformationFile == "")
+		InformationFile = StringLine::ReplaceAll(TMVAWeightPath,".xml",".nfo");
+
+	fMaximalContent = 0;
+	fInformationFile = InformationFile;
+	LoadKeyword();
+	ReadNFO();//Getting information from fMLPInformationFile
+	
+	if(fMaximalContent == 0 )
+	{
+		ERROR<<"Can't find the k_maxfiscontent keyword in .nfo file\n this is mandatory"<<endl;
+		exit(0);
+	}
+	
+	fTargetKeff = keff_target;
+	
+	SetPCMprecision(10);
+	SetBuildFuelFirstGuess(0.15);	//First fissile content guess for the EquivalenceModel::BuildFuel algorithm
+	fActualFissileContent = fFirstGuessFissilContent ;
+
+	INFO << "__An equivalence model has been define__" << endl;
+	INFO << "\tThis model is based on the prediction of keff at a specific time" << endl;
+	INFO << "\tThe TMVA (weight | information) files are :" << endl;
+	INFO << "\t" << "( " << fTMVAWeightPath[0] << " | " << fInformationFile << " )" << endl;
+
+	DBGL
+}
+//________________________________________________________________________
+EQM_FBR_MLP_Keff::EQM_FBR_MLP_Keff(CLASSLogger* log, string TMVAWeightPath, double keff_target, string InformationFile):EquivalenceModel(log)
+{
+	DBGL
+
+	/**The tmva weight*/
+	fTMVAWeightPath = TMVAWeightPath;
+
+
+	/* INFORMATION FILE HANDLING */
+	
+	if(InformationFile == "")
+		InformationFile = StringLine::ReplaceAll(TMVAWeightPath,".xml",".nfo");
+
+	fMaximalContent = 0;
+	fInformationFile = InformationFile;
+	LoadKeyword();
+	ReadNFO();//Getting information from fMLPInformationFile
+	
+	if(fMaximalContent == 0 )
+	{
+		ERROR<<"Can't find the k_maxfiscontent keyword in .nfo file\n this is mandatory"<<endl;
+		exit(0);
+	}
+
+	fTargetKeff = keff_target;
+	
+	SetPCMprecision(10);
+	SetBuildFuelFirstGuess(0.15);	//First fissile content guess for the EquivalenceModel::BuildFuel algorithm
+	fActualFissileContent = fFirstGuessFissilContent ;
+
+	INFO << "__An equivalence model has been define__" << endl;
+	INFO << "\tThis model is based on the prediction of keff at a specific time" << endl;
+	INFO << "\tThe TMVA (weight | information) files are :" << endl;
+	INFO << "\t" << "( " << fTMVAWeightPath[0] << " | "  << fInformationFile << " )" << endl;
+
+	DBGL
+}
+
+
+//________________________________________________________________________
+TTree* EQM_FBR_MLP_Keff::CreateTMVAInputTree(IsotopicVector TheFreshfuel, double ThisTime)
+{
+	DBGL
+
+	/******Create Input data tree to be interpreted by TMVA::Reader***/
+	TTree*   InputTree = new TTree("InTMPKef", "InTMPKef");
+
+	vector<float> 	InputTMVA;
+	for(int i = 0 ; i< (int)fMapOfTMVAVariableNames.size() ; i++)
+		InputTMVA.push_back(0);
+
+	float Time = 0;
+
+	IsotopicVector IVInputTMVA;
+	map<ZAI ,string >::iterator it;
+	int j = 0;
+
+	for( it = fMapOfTMVAVariableNames.begin()  ; it != fMapOfTMVAVariableNames.end() ; it++)
+	{
+		InputTree->Branch( ((*it).second).c_str()	,&InputTMVA[j], ((*it).second + "/F").c_str());
+		IVInputTMVA+=  ((*it).first)*1;
+		j++;
+	}
+
+	if(ThisTime != -1)
+		InputTree->Branch(	"Time"	,&Time	,"Time/F"	);
+
+	IsotopicVector IVAccordingToUserInfoFile = TheFreshfuel.GetThisComposition(IVInputTMVA);
+
+	double Ntot = IVAccordingToUserInfoFile.GetSumOfAll();
+
+	IVAccordingToUserInfoFile = IVAccordingToUserInfoFile/Ntot;
+
+	j = 0;
+	map<ZAI ,string >::iterator it2;
+
+	for( it2 = fMapOfTMVAVariableNames.begin() ; it2 != fMapOfTMVAVariableNames.end() ; it2++)
+	{
+		InputTMVA[j] = IVAccordingToUserInfoFile.GetZAIIsotopicQuantity( (*it2).first ) ;
+		j++;
+	}
+
+	Time = ThisTime;
+
+	InputTree->Fill();
+
+	DBGL
+	return InputTree;
+
+}
+//________________________________________________________________________
+double EQM_FBR_MLP_Keff::ExecuteTMVA(TTree* InputTree, bool IsTimeDependent)
+{
+	DBGL
+
+	// --- Create the Reader object
+	TMVA::Reader *reader = new TMVA::Reader( "Silent" );
+
+	// Create a set of variables and declare them to the reader
+	// - the variable names MUST corresponds in name and type to those given in the weight file(s) used
+	vector<float> 	InputTMVA;
+	for(int i = 0 ; i< (int)fMapOfTMVAVariableNames.size() ; i++)
+		InputTMVA.push_back(0);
+	
+	Float_t Time = 0;
+
+	
+	map<ZAI ,string >::iterator it;
+	int j = 0;
+	
+	
+	for( it = fMapOfTMVAVariableNames.begin()  ; it != fMapOfTMVAVariableNames.end() ; it++)
+	{	reader->AddVariable( ( (*it).second ).c_str(),&InputTMVA[j]);
+		j++;
+	}
+
+	if(IsTimeDependent)
+		reader->AddVariable( "Time" ,&Time);
+
+	
+	
+	// --- Book the MVA methods
+
+	// Book method MLP
+	TString methodName = "MLP method";
+	reader->BookMVA( methodName, fTMVAWeightPath );
+
+	map<ZAI ,string >::iterator it2;
+	j = 0;
+	for( it2 = fMapOfTMVAVariableNames.begin()  ; it2 != fMapOfTMVAVariableNames.end() ; it2++)
+	{
+		InputTree->SetBranchAddress(( (*it2).second ).c_str(),&InputTMVA[j]);
+		j++;
+	}
+
+	if(IsTimeDependent)
+		InputTree->SetBranchAddress( "Time" ,&Time );
+
+	InputTree->GetEntry(0);
+	
+	Float_t val = (reader->EvaluateRegression( methodName ))[0];
+
+	delete reader;
+
+	DBGL
+	return (double)val;	//return k_{eff}(t = Time)
+}
+
+
+//________________________________________________________________________
+void EQM_FBR_MLP_Keff::LoadKeyword()
+{
+	DBGL
+
+	fDKeyword.insert( pair<string, FBR_MLP_Keff_DMthPtr>( "k_zainame",	& EQM_FBR_MLP_Keff::ReadZAIName)	 );
+
+	fDKeyword.insert( pair<string, FBR_MLP_Keff_DMthPtr>( "k_maxfiscontent", &EQM_FBR_MLP_Keff::ReadMaxFisContent) );
+
+	DBGL
+}
+
+
+//________________________________________________________________________
+void EQM_FBR_MLP_Keff::ReadZAIName(const string &line)
+{
+	DBGL
+	int pos = 0;
+	string keyword = tlc(StringLine::NextWord(line, pos, ' '));
+	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());
+	
+	fFissileList.Add(Z, A, I, 1.0);
+	fStreamList.push_back(fFissileList);
+	
+	DBGL
+}
+
+//________________________________________________________________________
+void EQM_FBR_MLP_Keff::ReadMaxFisContent(const string &line)
+{
+	DBGL
+	int pos = 0;
+	string keyword = tlc(StringLine::NextWord(line, pos, ' '));
+	if( keyword != "k_maxfiscontent" )	// Check the keyword
+	{
+		ERROR << " Bad keyword : \"k_maxfiscontent\" not found !" << endl;
+		exit(1);
+	}
+	
+	fMaximalContent = atof(StringLine::NextWord(line, pos, ' ').c_str());
+	
+	DBGL
+}
+
+//________________________________________________________________________
+void EQM_FBR_MLP_Keff::ReadLine(string line)
+{
+	DBGL
+	
+	int pos = 0;
+	string keyword = tlc(StringLine::NextWord(line, pos, ' '));
+	
+	map<string, FBR_MLP_Keff_DMthPtr>::iterator it = fDKeyword.find(keyword);
+	
+	if(it != fDKeyword.end())
+		(this->*(it->second))( line );
+	
+	DBGL
+}
+
+
+
+
+//________________________________________________________________________
+map < string , double> EQM_FBR_MLP_Keff::GetMolarFraction(vector <IsotopicVector> IVStream,double TargetBU)
+{
+	DBGL
+	
+	if(TargetBU != 0)
+		WARNING << "The third arguement : Burnup has no effect here.";
+
+	IsotopicVector Fissile = IVStream["Fissile"];
+	IsotopicVector Fertile = IVStream["Fertile"];
+
+	//initialization
+	double FissileContent = GetActualFissileContent();
+	double OldFissileContentMinus = 0;
+	double OldFissileContentPlus = fMaximalContent;
+	double PredictedKeff = 0 ;
+	IsotopicVector FreshFuel = (1-FissileContent)*(Fertile/Fertile.GetSumOfAll()) + FissileContent*(Fissile/Fissile.GetSumOfAll());
+	double OldPredictedKeff = GetKeffAtFixedTime(FreshFuel);
+	
+	double Precision = fPCMprecision/1e5*fTargetKeff; //pcm to 1
+	
+	int count = 0;
+	int MaximumLoopCount = 100;
+	do
+	{
+		if(count > MaximumLoopCount )
+		{
+			ERROR << "CRITICAL ! Can't manage to predict fissile content\nHint : Try to decrease the precision on keff using :\nYourEQM_FBR_MLP_Keff->SetPCMPrecision(prop); with prop the precision  (default 0.5percent :  0.005) INCREASE IT \n If this message still appear mail to leniau@subatech.in2p3.fr\nor nicolas.thiolliere@subatech.in2p3.fr " << endl;
+			exit(1);
+		}
+		
+		if( (OldPredictedKeff - fTargetKeff) < 0 ) //The Content can be increased
+		{
+			OldFissileContentMinus = FissileContent;
+			FissileContent = FissileContent + fabs(OldFissileContentPlus-FissileContent)/2.;
+		}
+		else if( (OldPredictedKeff - fTargetKeff) > 0) //The Content is too high
+		{
+			OldFissileContentPlus = FissileContent;
+			FissileContent = FissileContent - fabs(OldFissileContentMinus-FissileContent)/2.;
+		}
+		
+		IsotopicVector FreshFuel = (1-FissileContent)*(Fertile/Fertile.GetSumOfAll()) + FissileContent*(Fissile/Fissile.GetSumOfAll());
+		
+		PredictedKeff = GetKeffAtFixedTime(FreshFuel);
+		
+		OldPredictedKeff = PredictedKeff;
+		count ++;
+		
+	}while(fabs(fTargetKeff-PredictedKeff)>Precision);
+	
+	DBGV( "Predicted keff " << PredictedKeff << " FissileContent " << FissileContent << endl);
+	
+	map < string , double> MolarFraction;
+	MolarFraction["Fissile"] = FissileContent;
+	MolarFraction["Fertile"] = 1.- FissileContent;
+
+	return MolarFraction; //return Molar content of each component in the fuel
+
+
+}
+//________________________________________________________________________

source/branches/CLASSV4.2/Model/Equivalence/EQM_FBR_MLP_Keff.hxx

+#ifndef _EQM_FBR_MLP_Keff_HXX
+#define _EQM_FBR_MLP_Keff_HXX
+
+#include "EquivalenceModel.hxx"
+#include "TTree.h"
+#include "TGraph.h"
+#include <string>
+#include <fstream>
+#include <iostream>
+#include <map>
+#include <vector>
+
+using namespace std;
+
+class EQM_FBR_MLP_Keff;
+#ifndef __CINT__
+typedef void (EQM_FBR_MLP_Keff::*FBR_MLP_Keff_DMthPtr)( const string & ) ;
+#endif
+
+//-----------------------------------------------------------------------------//
+//! Defines an EquivalenceModel based on neural network to predict @f$k_{eff}@f$.
+
+/*!
+The aim of these class is to constuct a fuel from an equivalence model
+based on a  Multi layer perceptron (MLP).
+This MLP aims to predict the Pu content such as it has to verify
+ @f$ k(WantedTime) = k_{target}@f$  with @f$k_{target}@f$ is often close to 1.0 
+ but can be set by user. The wanted time is often either 
+ the begining of cycle or end of cycle. WantedTime can't be set by user since it is
+ contain in the .xml file. Indeed this method suppose you have trained your MLP to predict 
+ the @f$k_{eff}@f$ either at BOC or EOC (or any other time)
+
+ @author BLG
+ @author BaM
+ @version 1.0
+ */
+//________________________________________________________________________
+
+class EQM_FBR_MLP_Keff : public EquivalenceModel
+{
+	public:
+	/*!
+	 \name Constructor
+	 */
+	//@{
+
+
+
+	//{
+	///  Create a EQM_FBR_MLP_Keff using Keffective at a given time
+ 	/// (see class desctiption)
+	 /*!
+	 Create a EQM_FBR_MLP_Keff 
+	 \param  TMVAWeightPath :  Path to the .xml file containing neural network informations for prediction of keff(t = tfixed)
+	 \param  keff_target : Wanted  @f$k_{eff}@f$ (see detailed description)
+	 \param  InformationFile : Total path to the file containing time steps, fissile and ferile list (ante and post fabrication time cooling). Default is the same total path as TMVAWeightPath but extension is replaced by .nfo (see manual for format)
+	 */
+	EQM_FBR_MLP_Keff(string TMVAWeightPath, double keff_target = 1.00, string InformationFile = "");
+	//}
+
+	//{
+	///  Create a EQM_FBR_MLP_Keff using Keffective at a given time
+ 	/// (see class desctiption)
+	 /*!
+	 Create a EQM_FBR_MLP_Keff 
+	 \param  log:  CLASSLogger object to handle log messages
+	 \param  TMVAWeightPath :  Path to the .xml file containing neural network informations for prediction of keff(t = tfixed)
+	 \param  keff_target : Wanted  @f$k_{eff}@f$ (see detailed description)
+	 \param  InformationFile : Total path to the file containing time steps, fissile and ferile list (ante and post fabrication time cooling). Default is the same total path as TMVAWeightPath but extension is replaced by .nfo (see manual for format)
+	 */
+	EQM_FBR_MLP_Keff(CLASSLogger* log, string TMVAWeightPath, double keff_target = 1.00, string InformationFile = "");
+	//}
+
+	//@}
+	
+	//{
+	/// Return the molar fissile fraction according fissile & ferile content using @f$<k_{\infty}>(t)@f$ prediction
+	/*!
+	 \param Fissil : The composition of the fissile matter
+	 \param Fertil : The composition of the Fertil matter
+	 \param BurnUp : Maximum achievable burn up envisaged
+	 */
+	virtual double GetFissileMolarFraction(vector <IsotopicVector> IVStream, double BurnUp = 0);
+	//}
+	
+	/*!
+	 \name Get/Set methods
+	 */
+	//@{
+	void 	SetPCMprecision(double pcm){fPCMprecision = pcm;}	//!< Set the precision on @f$\langle k \rangle@f$ prediction [pcm]. Neural network predictor constructors
+	double 	GetPCMprecision(){return fPCMprecision/1e5;}		//!< Get the precision on @f$\langle k \rangle@f$ prediction []. Neural network predictor constructors
+
+	//@}
+	
+	
+	
+	
+	
+	/*!
+	 \name TMVA related  methods
+	 */
+	//@{
+	TTree* 	CreateTMVAInputTree(IsotopicVector FreshFuel, double ThisTime);//!<Create input tmva tree to be read by ExecuteTMVA
+	
+	double 	ExecuteTMVA(TTree* theTree, bool IsTimeDependant);//!<Execute the MLP according to the input tree created by CreateTMVAInputTree
+	//@}
+	
+	
+	
+	/*!
+	 \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 ReadZAIName(const string &line);
+	//}
+	
+	//{
+	/// ReadMaxFisContent  : read a guessed (very overestimated) maximum fissile content (purpose : algorithm initialization)
+	/*!
+	 \param line : line suppossed to contain the ZAI name  starts with "k_zainame" keyword
+	 */
+	void ReadMaxFisContent(const string &line);
+	//}
+	
+	//{
+	/// ReadLine : read a line
+	/*!
+	 \param line : line to read
+	 */
+	void ReadLine(string line);
+	//}
+	
+	//@}
+	
+	
+	
+	private :
+
+	string fTMVAWeightPath;					//!<The weight needed by TMVA to construct and execute the multilayer perceptron
+
+#ifndef __CINT__
+	map<string, FBR_MLP_Keff_DMthPtr> fDKeyword;
+#endif
+
+
+	map<ZAI,string> fMapOfTMVAVariableNames;//!<  List of TMVA input variable names (read from fMLPInformationFile ) , name depends on the training step
+
+	vector<double> fMLP_Time;	//!< Time (in seconds) when the MLP(t) = keff(t) has been trained.
+	double 	fSpecificPower; 	//!< The specific power in W/gHM (HM: heavy Metal)
+	double  fMaximalContent;	//!< The approx. maximum fissile content reachable by the MLP model
+	
+
+	
+	int 	fNumberOfBatch;		//!< The number of batches for the loading plan
+	
+	double 	fKThreshold;		//!< The @f$k_{Threshold}@f$
+	double 	fPCMprecision;		//!< precision on @f$\langle k \rangle@f$ prediction [pcm]
+	
+	double 	fTargetKeff;		//!< Use for Varying Fissile content to reach fTargetKeff at time used in the MLP Training
+	
+
+	
+	double 	GetKeffAtFixedTime(IsotopicVector FreshFuel){return ExecuteTMVA( CreateTMVAInputTree(FreshFuel,-1), false );} //!<time independant since the MLP is trained for 1 time
+
+	TGraph* BuildKeffGraph(IsotopicVector FreshFuel);
+	TGraph* BuildAverageKeffGraph(TGraph* GRAPH_KEFF);
+
+	double 	GetKeffAt(TGraph* GRAPH_KEFF, int Step);
+
+
+
+	
+};
+
+#endif
+

source/branches/CLASSV4.2/Model/Equivalence/EQM_FBR_MLP_Keff_BOUND.cxx

+#include "EQM_FBR_MLP_Keff_BOUND.hxx"
+#include "CLASSMethod.hxx"
+#include "CLASSLogger.hxx"
+#include "StringLine.hxx"
+
+#include <string>
+#include <iostream>
+#include <fstream>
+#include <algorithm>
+#include <cmath>
+#include <cassert>
+
+#include "TSystem.h"
+#include "TMVA/Reader.h"
+#include "TMVA/Tools.h"
+#include "TMVA/MethodCuts.h"
+
+//________________________________________________________________________
+//
+//		EQM_FBR_MLP_Keff_BOUND
+//
+//	Equivalenve Model based on multi layer perceptron from TMVA (root cern)
+//	For FBR
+//
+//________________________________________________________________________
+
+
+
+//________________________________________________________________________
+//
+//	Objects & Methods for content prediction with calculation of Keff bounded
+//	over batches.
+//
+//________________________________________________________________________
+
+//________________________________________________________________________
+EQM_FBR_MLP_Keff_BOUND::EQM_FBR_MLP_Keff_BOUND(string TMVAWeightPath,  int NumOfBatch, double LowerKeffective, double UpperKeffective, string InformationFile):EquivalenceModel(new CLASSLogger("EQM_FBR_MLP_Keff_BOUND.log"))
+{
+	DBGL
+
+	/** The tmva weight **/
+	
+	fTMVAWeightPath = TMVAWeightPath;
+	
+	
+	/* INFORMATION FILE HANDLING */
+	
+	if(InformationFile == "")
+		InformationFile = StringLine::ReplaceAll(TMVAWeightPath,".xml",".nfo");
+	
+	fInformationFile = InformationFile;
+	
+	LoadKeyword();
+	ReadNFO();//Getting information from fMLPInformationFile
+	
+	
+	/* OTHER MODEL PARAMETERS */
+	
+	fNumberOfBatch = NumOfBatch;
+	fKmin = LowerKeffective ;
+	fKmax = UpperKeffective ;
+	
+
+	/* MODEL PARAMETERS INITIALIZATION */
+	
+	SetPCMprecision(10);
+	SetBuildFuelFirstGuess(0.15);//First fissile content guess for the EquivalenceModel::BuildFuel algorithm
+	fActualFissileContent = fFirstGuessFissilContent ;
+	
+	
+	/* INFO */
+
+	INFO << "__An equivalence model has been define__" << endl;
+	INFO << "\tThis model is based on the prediction of keff averaged over the number of batch" << endl;
+	INFO << "\tThe TMVA (weight | information) files are :" << endl;
+	INFO << "\t" << "( " << fTMVAWeightPath[0] << " | " << fMLPInformationFile << " )" << endl;
+
+	
+	DBGL
+}
+//________________________________________________________________________
+EQM_FBR_MLP_Keff_BOUND::EQM_FBR_MLP_Keff_BOUND(CLASSLogger* log, string TMVAWeightPath,  int NumOfBatch, double LowerKeffective, double UpperKeffective, string InformationFile):EquivalenceModel(log)
+{
+	DBGL
+	
+	/** The tmva weight **/
+	
+	fTMVAWeightPath = TMVAWeightPath;
+	
+	
+	/* INFORMATION FILE HANDLING */
+	
+	if(InformationFile == "")
+		InformationFile = StringLine::ReplaceAll(TMVAWeightPath,".xml",".nfo");
+	
+	fInformationFile = InformationFile;
+	LoadKeyword();
+	ReadNFO();//Getting information from fMLPInformationFile
+	
+	
+	
+	/* OTHER MODEL PARAMETERS */
+	
+	fNumberOfBatch = NumOfBatch;
+	fKmin = LowerKeffective ;
+	fKmax = UpperKeffective ;
+	
+	
+	/* INFO */
+	
+	SetPCMprecision(10);
+	SetBuildFuelFirstGuess(0.15);//First fissile content guess for the EquivalenceModel::BuildFuel algorithm
+	fActualFissileContent = fFirstGuessFissilContent ;
+	
+	INFO << "__An equivalence model has been define__" << endl;
+	INFO << "\tThis model is based on the prediction of keff averaged over the number of batch" << endl;
+	INFO << "\tThe TMVA (weight | information) files are :" << endl;
+	INFO << "\t" << "( " << fTMVAWeightPath[0] << " | " << fMLPInformationFile << " )" << endl;
+	
+	
+	DBGL
+}
+
+
+//________________________________________________________________________
+TGraph* EQM_FBR_MLP_Keff_BOUND::BuildKeffGraph(IsotopicVector FreshFuel)
+{
+	DBGL
+	
+	TGraph * keffGraph = new TGraph();
+	for(int i = 0 ; i < (int) fMLP_Time.size() ; i++)
+	{
+		double keff_t = ExecuteTMVA( CreateTMVAInputTree(FreshFuel,(float) fMLP_Time[i]), true );
+		keffGraph->SetPoint(i, (double)fMLP_Time[i], keff_t );
+	}
+	
+	DBGL
+	return keffGraph;
+}
+//________________________________________________________________________
+TGraph* EQM_FBR_MLP_Keff_BOUND::BuildAverageKeffGraph(TGraph* GRAPH_KEFF)
+{
+	DBGL
+	
+	TGraph * AveragekeffGraph = new TGraph();
+	int NumberOfPoint = 50;
+	
+	int NumOfInputGraphPoint = GRAPH_KEFF->GetN()-1 ;
+	double TimeFinal = 0;
+	double KFinal = 0;
+	
+	GRAPH_KEFF->GetPoint(NumOfInputGraphPoint, TimeFinal, KFinal);
+	
+	double step = TimeFinal/NumberOfPoint;
+	int p = 0;
+	for(int n = 0 ;n<NumberOfPoint;n++)
+	{
+		double k_av = 0;
+		for(int b = 0;b<fNumberOfBatch;b++)
+		{
+			if((step*p)> TimeFinal/(double)fNumberOfBatch)
+				p = 0;
+			k_av += GRAPH_KEFF->Eval( (step*p + b*TimeFinal/(double)fNumberOfBatch) , 0 , "S" );
+			
+		}
+		p++;
+		k_av/= (double)fNumberOfBatch;
+		
+		AveragekeffGraph->SetPoint(n, step*n, k_av);
+	}
+	
+	DBGL
+	return AveragekeffGraph;
+}
+//________________________________________________________________________
+double EQM_FBR_MLP_Keff_BOUND::GetKeffAt(TGraph* GRAPH_KEFF, int Step)
+{
+	DBGL
+	
+	double Time = 0;
+	double Keff = 0;
+	GRAPH_KEFF->GetPoint(Step, Time,  Keff);
+	
+	DBGL
+	return Keff;
+}
+
+
+//________________________________________________________________________
+TTree* EQM_FBR_MLP_Keff_BOUND::CreateTMVAInputTree(IsotopicVector TheFreshfuel, double ThisTime)
+{
+	DBGL
+
+	/******Create Input data tree to be interpreted by TMVA::Reader***/
+	TTree*   InputTree = new TTree("InTMPKef", "InTMPKef");
+	
+	vector<float> 	InputTMVA;
+	for(int i = 0 ; i< (int)fMapOfTMVAVariableNames.size() ; i++)
+		InputTMVA.push_back(0);
+	
+	float Time = 0;
+	
+	IsotopicVector IVInputTMVA;
+	map<ZAI ,string >::iterator it;
+	int j = 0;
+	
+	for( it = fMapOfTMVAVariableNames.begin()  ; it != fMapOfTMVAVariableNames.end() ; it++)
+	{
+		InputTree->Branch( ((*it).second).c_str()	,&InputTMVA[j], ((*it).second + "/F").c_str());
+		IVInputTMVA+=  ((*it).first)*1;
+		j++;
+	}
+	
+	if(ThisTime != -1)
+		InputTree->Branch(	"Time"	,&Time	,"Time/F"	);
+	
+	IsotopicVector IVAccordingToUserInfoFile = TheFreshfuel.GetThisComposition(IVInputTMVA);
+	
+	double Ntot = IVAccordingToUserInfoFile.GetSumOfAll();
+	
+	IVAccordingToUserInfoFile = IVAccordingToUserInfoFile/Ntot;
+	
+	j = 0;
+	map<ZAI ,string >::iterator it2;
+	
+	for( it2 = fMapOfTMVAVariableNames.begin() ; it2 != fMapOfTMVAVariableNames.end() ; it2++)
+	{
+		InputTMVA[j] = IVAccordingToUserInfoFile.GetZAIIsotopicQuantity( (*it2).first ) ;
+		j++;
+	}
+	
+	Time = ThisTime;
+	
+	InputTree->Fill();
+	
+	DBGL
+	return InputTree;
+	
+}
+//________________________________________________________________________
+double EQM_FBR_MLP_Keff_BOUND::ExecuteTMVA(TTree* InputTree, bool IsTimeDependent)
+{
+	DBGL
+
+	// --- Create the Reader object
+	TMVA::Reader *reader = new TMVA::Reader( "Silent" );
+	
+	// Create a set of variables and declare them to the reader
+	// - the variable names MUST corresponds in name and type to those given in the weight file(s) used
+	vector<float> 	InputTMVA;
+	for(int i = 0 ; i< (int)fMapOfTMVAVariableNames.size() ; i++)
+		InputTMVA.push_back(0);
+	Float_t Time;
+	
+	map<ZAI ,string >::iterator it;
+	int j = 0;
+	for( it = fMapOfTMVAVariableNames.begin()  ; it != fMapOfTMVAVariableNames.end() ; it++)
+	{	reader->AddVariable( ( (*it).second ).c_str(),&InputTMVA[j]);
+		j++;
+	}
+	
+	if(IsTimeDependent)
+		reader->AddVariable( "Time" ,&Time);
+	
+	// --- Book the MVA methods
+	
+	// Book method MLP
+	TString methodName = "MLP method";
+	reader->BookMVA( methodName, fTMVAWeightPath );
+	
+	map<ZAI ,string >::iterator it2;
+	j = 0;
+	for( it2 = fMapOfTMVAVariableNames.begin()  ; it2 != fMapOfTMVAVariableNames.end() ; it2++)
+	{
+		InputTree->SetBranchAddress(( (*it2).second ).c_str(),&InputTMVA[j]);
+		j++;
+	}
+	
+	if(IsTimeDependent)
+		InputTree->SetBranchAddress( "Time" ,&Time );
+	
+	InputTree->GetEntry(0);
+	Float_t val = (reader->EvaluateRegression( methodName ))[0];
+	
+	delete reader;
+	
+	DBGL
+	return (double)val;	//return k_{eff}(t = Time)
+}
+
+//________________________________________________________________________
+void EQM_FBR_MLP_Keff_BOUND::LoadKeyword()
+{
+	DBGL
+
+	fDKeyword.insert( pair<string, FBR_MLP_Keff_BOUND_DMthPtr>( "k_timestep",	& EQM_FBR_MLP_Keff_BOUND::ReadTimeSteps));
+	fDKeyword.insert( pair<string, FBR_MLP_Keff_BOUND_DMthPtr>( "k_zainame",	& EQM_FBR_MLP_Keff_BOUND::ReadZAIName)	 );
+	
+	DBGL
+}
+
+
+//________________________________________________________________________
+void EQM_FBR_MLP_Keff_BOUND::ReadZAIName(const string &line)
+{
+	DBGL
+	int pos = 0;
+	string keyword = tlc(StringLine::NextWord(line, pos, ' '));
+	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 ) );
+	DBGL
+}
+
+
+//________________________________________________________________________
+void EQM_FBR_MLP_Keff_BOUND::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 EQM_FBR_MLP_Keff_BOUND::ReadLine(string line)
+{
+	DBGL
+	
+	int pos = 0;
+	string keyword = tlc(StringLine::NextWord(line, pos, ' '));
+	
+	map<string, FBR_MLP_Keff_BOUND_DMthPtr>::iterator it = fDKeyword.find(keyword);
+	
+	if(it != fDKeyword.end())
+		(this->*(it->second))( line );
+	
+	DBGL
+}
+
+
+//________________________________________________________________________
+double EQM_FBR_MLP_Keff_BOUND::GetFissileMolarFraction(vector <IsotopicVector> IVStream, double TargetBU)
+{
+	DBGL
+
+	IsotopicVector Fissile = IVStream[0];
+	IsotopicVector Fertile = IVStream[1];
+
+	if(TargetBU != 0)
+		WARNING << "The third arguement : Burnup has no effect here.";
+	
+	/**Algorithm  not so clever ...**/
+	/**need improvements to make it faster*/
+	
+	double FissileContent = 0.01;
+	double test_Keff_beg = fKmin - 0.01;
+	double test_keff_end = fKmax + 0.01;
+	
+	double speedstep = 1;
+	
+	while(fKmin >= test_Keff_beg || fKmax <=  test_keff_end)
+	{
+		IsotopicVector FreshFuel = (1-FissileContent)*(Fertile/Fertile.GetSumOfAll()) + FissileContent*(Fissile/Fissile.GetSumOfAll());
+		
+		TGraph* KEFF  = BuildKeffGraph(FreshFuel);
+		TGraph* KEFF_avg = BuildAverageKeffGraph(KEFF);
+		
+		test_Keff_beg = GetKeffAt(KEFF_avg , 0);
+		test_keff_end = GetKeffAt(KEFF_avg , fMLP_Time.size()-1);
+		delete KEFF;
+		delete KEFF_avg;
+		
+		if(test_Keff_beg < 0.9) //why 0.9 ? exactly
+			speedstep = 10;
+		else
+			speedstep = 1;
+		
+		FissileContent+= 0.001*speedstep;
+		
+		if( test_Keff_beg > 1.30 )
+		{
+			ERROR << "This plutonium can not satisfy the criticality condition imposed" << endl;
+			FissileContent = -1;
+			break;
+		}
+		
+	}
+	
+	DBGL
+	return FissileContent;
+}
+//________________________________________________________________________

source/branches/CLASSV4.2/Model/Equivalence/EQM_FBR_MLP_Keff_BOUND.hxx

+#ifndef _EQM_FBR_MLP_Keff_BOUND_HXX_
+#define _EQM_FBR_MLP_Keff_BOUND_HXX_
+
+#include "EquivalenceModel.hxx"
+#include "TTree.h"
+#include "TGraph.h"
+
+using namespace std;
+
+class EQM_FBR_MLP_Keff_BOUND;
+#ifndef __CINT__
+typedef void (EQM_FBR_MLP_Keff_BOUND::*FBR_MLP_Keff_BOUND_DMthPtr)( const string & ) ;
+#endif
+
+//-----------------------------------------------------------------------------//
+//! Defines an EquivalenceModel based on neural network to predict @f$k_{\infty}@f$
+
+/*!
+ The aim of these class is to constuct a fuel from an equivalence model
+ based on a  Multi layer perceptron (MLP).
+ This MLP aims to predict :
+ 
+ The @f$k_{\infty}(t)@f$ of a FBR-MOX from a given fresh fuel composition.
+ With this MLP prediction and a given number of batch (for the loading plan) an
+ average @f$<k_{\infty}(t)>@f$ is calculated according :
+ @f$<k_{\infty}>^{batch}(t) = \frac{1}{N}\sum_{i}^{N} k_{\infty}(t + \frac{iT}{N} )@f$
+ The Fissile content has to verify this condition :
+ @f$ k_{\infty Max} \geq <k_{\infty}>^{batch}(T/N) \geq k_{\infty Min} @f$
+ Where @f$ k_{\infty Max}@f$ and @f$k_{\infty Min}@f$ are arguments of the constructor.
+ 
+ 
+ \warning
+ it is not guaranted that there is a solution for Pu content verifying :
+ @f$<k_{\infty}>^{batch}(t) = \frac{1}{N}\sum_{i}^{N}k_{\infty}(t+\frac{iT}{N} )@f$
+ 
+ @author BLG
+ @author BaM
+ @version 1.0
+ */
+//________________________________________________________________________
+
+class EQM_FBR_MLP_Keff_BOUND : public EquivalenceModel
+{
+	public:
+	/*!
+	 \name Constructor
+	 */
+	//@{
+	
+	//{
+	///  Create a EQM_FBR_MLP_Keff_BOUND using @f$k_{\infty}@f$ average over batch
+	/// (see class desctiption)
+	/*!
+	 Create a EQM_FBR_MLP_Keff_BOUND using @f$k_{\infty}@f$ average over batch
+	 \param  TMVAWeightPath0:  Path to the .xml file containing neural network informations for prediction of keff(t)
+	 \param  NumOfBatch : Number of batch for the loading plan (often 5 for FBR-Na)
+	 \param  KeffMin : Lower average @f$<k_{\infty}>@f$ value
+	 \param  KeffMax : Upper average @f$<k_{\infty}>@f$ value
+	 \param  InformationFile : Total path to the file containing time steps, fissile and ferile list (ante and post fabrication time cooling). Default is the same total path as TMVAWeightPath but extension is replaced by .nfo (see manual for format)
+	 */
+	EQM_FBR_MLP_Keff_BOUND(string TMVAWeightPath,  int NumOfBatch , double LowerKeffective, double UpperKeffective, string InformationFile = "");
+	//}
+	
+	//{
+	///  Create a EQM_FBR_MLP_Keff_BOUND using @f$k_{\infty}@f$ average over batch
+	/// (see class desctiption)
+	/*!
+	 Create a EQM_FBR_MLP_Keff_BOUND using @f$k_{\infty}@f$ average over batch
+	 \param  log:  CLASSLogger object to handle log messages
+	 \param  TMVAWeightPath0:  Path to the .xml file containing neural network informations for prediction of keff(t)
+	 \param  NumOfBatch : Number of batch for the loading plan (often 5 for FBR-Na)
+	 \param  KeffMin : Lower average @f$<k_{\infty}>@f$ value
+	 \param  KeffMax : Upper average @f$<k_{\infty}>@f$ value
+	 \param  InformationFile : Total path to the file containing time steps, fissile and ferile list (ante and post fabrication time cooling). Default is the same total path as TMVAWeightPath but extension is replaced by .nfo (see manual for format)
+	 */
+	EQM_FBR_MLP_Keff_BOUND(CLASSLogger* log, string TMVAWeightPath,  int NumOfBatch , double LowerKeffective, double UpperKeffective, string InformationFile = "");
+	//}
+	
+	
+	//@}
+	
+	//{
+	/// Return the molar fissile fraction according fissile & ferile content using @f$<k_{\infty}>(t)@f$ prediction
+	/*!
+	 \param Fissil : The composition of the fissile matter
+	 \param Fertil : The composition of the Fertil matter
+	 \param BurnUp : Maximum achievable burn up envisaged
+	 */
+	virtual double GetFissileMolarFraction(vector <IsotopicVector> IVStream,double BurnUp = 0);
+	//}
+	
+	/*!
+	 \name Get/Set methods
+	 */
+	//@{
+	void 	SetPCMprecision(double pcm){fPCMprecision = pcm;}	//!< Set the precision on @f$\langle k \rangle@f$ prediction [pcm]. Neural network predictor constructors
+	double 	GetPCMprecision(){return fPCMprecision/1e5;}		//!< Get the precision on @f$\langle k \rangle@f$ prediction []. Neural network predictor constructors
+	
+	//@}
+	
+	
+	
+	
+	
+	/*!
+	 \name TMVA related  methods
+	 */
+	//@{
+	TTree* 	CreateTMVAInputTree(IsotopicVector FreshFuel, double ThisTime);//!<Create input tmva tree to be read by ExecuteTMVA
+	
+	double 	ExecuteTMVA(TTree* theTree, bool IsTimeDependant);//!<Execute the MLP according to the input tree created by CreateTMVAInputTree
+	//@}
+	
+	
+	
+	/*!
+	 \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 ReadZAIName(const string &line);
+	//}
+	
+
+	
+	//{
+	/// ReadLine : read a line
+	/*!
+	 \param line : line to read
+	 */
+	void ReadLine(string line);
+	//}
+	
+	//@}
+	
+	
+	
+	private :
+	
+	string fTMVAWeightPath;					//!<The weight needed by TMVA to construct and execute the multilayer perceptron
+	string fMLPInformationFile;				//!<The path to the informations necessary to execute the MLP
+	
+#ifndef __CINT__
+	map<string, FBR_MLP_Keff_BOUND_DMthPtr> fDKeyword;
+#endif
+	
+	map<ZAI,string> fMapOfTMVAVariableNames;//!<  List of TMVA input variable names (read from fMLPInformationFile ) , name depends on the training step
+	
+	vector<double> fMLP_Time;	//!< Time (in seconds) when the MLP(t) = keff(t) has been trained.
+	
+	
+	
+	
+	
+	int 	fNumberOfBatch;		//!< The number of batches for the loading plan
+	
+	double 	fKThreshold;		//!< The @f$k_{Threshold}@f$
+	double 	fPCMprecision;		//!< precision on @f$\langle k \rangle@f$ prediction [pcm]
+	double  fKmin;				//!< Lower edge of kedd Used by second constructor (fissile content prediction using keff at BOC (or other time)
+	double  fKmax;				//!< Upper edge of kedd Used by second constructor (fissile content prediction using keff at BOC (or other time)
+	double 	fTargetKeff;		//!< Use for Varying Fissile content to reach fTargetKeff at time used in the MLP Training
+	
+	
+	
+	
+	
+	
+	double 	GetKeffAtFixedTime(IsotopicVector FreshFuel){return ExecuteTMVA( CreateTMVAInputTree(FreshFuel,-1), false );} //!<time independant since the MLP is trained for 1 time
+	
+	TGraph* BuildKeffGraph(IsotopicVector FreshFuel);
+	TGraph* BuildAverageKeffGraph(TGraph* GRAPH_KEFF);
+	
+	double 	GetKeffAt(TGraph* GRAPH_KEFF, int Step);
+	
+	
+	
+	
+};
+
+#endif
+
+
+
+
+
+
+
+

source/branches/CLASSV4.2/Model/Equivalence/EQM_MLP_Kinf.cxx

+#include "EQM_MLP_Kinf.hxx"
+#include "CLASSLogger.hxx"
+#include "CLASSMethod.hxx"
+#include "StringLine.hxx"
+
+#include <string>
+#include <iostream>
+#include <fstream>
+#include <algorithm>
+#include <cmath>
+#include <cassert>
+
+#include "TSystem.h"
+#include "TMVA/Reader.h"
+#include "TMVA/Tools.h"
+#include "TMVA/MethodCuts.h"
+
+
+//________________________________________________________________________
+//
+//		EQM_MLP_Kinf
+//
+//	Equivalenve Model based on multi layer perceptron from TMVA (root cern)
+//	For REP MOX use
+//
+//________________________________________________________________________
+EQM_MLP_Kinf::EQM_MLP_Kinf(string WeightPathAlpha0, string WeightPathAlpha1, string WeightPathAlpha2, string InformationFile,  int NumOfBatch, double CriticalityThreshold):EquivalenceModel(new CLASSLogger("EQM_MLP_Kinf.log"))
+{
+	/**The information file and tmva weight*/
+	fTMVAWeightPath.push_back(WeightPathAlpha0);
+	fTMVAWeightPath.push_back(WeightPathAlpha1);
+	fTMVAWeightPath.push_back(WeightPathAlpha2);
+	
+	fMaximalBU = 0;
+	fMaximalContent = 0;
+	fInformationFile = InformationFile;
+	LoadKeyword();
+	ReadNFO();//Getting information from fMLPInformationFile
+	
+	if( fMaximalBU == 0 || fMaximalContent == 0 )
+	{
+		ERROR<<"Can't find the k_maxfiscontent and/or k_maxburnup keyword(s) in .nfo file\n this is mandatory"<<endl;
+		exit(0);
+	}
+	
+	fNumberOfBatch = NumOfBatch;
+	fKThreshold = CriticalityThreshold ;
+	SetBurnUpPrecision(0.005);//1 % of the targeted burnup
+	SetBuildFuelFirstGuess(0.04);//First fissile content guess for the EquivalenceModel::BuildFuel algorithm
+	
+	INFO << "__An equivalence model has been define__" << endl;
+	INFO << "\tThis model is based on the prediction of kinf" << endl;
+	INFO << "\tThe TMVA weight  files are :" << endl;
+	INFO << "\t" << fTMVAWeightPath[0] << endl;
+	INFO << "\t" << fTMVAWeightPath[1] << endl;
+	INFO << "\t" << fTMVAWeightPath[2] << endl;
+	INFO << "\tThe Information file is :" << endl;
+	INFO << "\t" << fInformationFile << endl;
+	
+}
+//________________________________________________________________________
+EQM_MLP_Kinf::EQM_MLP_Kinf(CLASSLogger* log, string WeightPathAlpha0, string WeightPathAlpha1, string WeightPathAlpha2, string InformationFile,  int NumOfBatch, double CriticalityThreshold):EquivalenceModel(log)
+{
+	/**The information file and tmva weight*/
+	fTMVAWeightPath.push_back(WeightPathAlpha0);
+	fTMVAWeightPath.push_back(WeightPathAlpha1);
+	fTMVAWeightPath.push_back(WeightPathAlpha2);
+	
+	fMaximalBU = 0;
+	fMaximalContent = 0;
+	fInformationFile = InformationFile;
+	LoadKeyword();
+	ReadNFO();//Getting information from fMLPInformationFile
+	
+	if( fMaximalBU == 0 || fMaximalContent == 0 )
+	{
+		ERROR<<"Can't find the k_maxfiscontent and/or k_maxburnup keyword(s) in .nfo file\n this is mandatory"<<endl;
+		exit(0);
+	}
+	
+	fNumberOfBatch = NumOfBatch;
+	fKThreshold = CriticalityThreshold ;
+	SetBurnUpPrecision(0.005);//1 % of the targeted burnup
+	SetBuildFuelFirstGuess(0.04);//First fissile content guess for the EquivalenceModel::BuildFuel algorithm
+	
+	INFO << "__An equivalence model has been define__" << endl;
+	INFO << "\tThis model is based on the prediction of kinf" << endl;
+	INFO << "\tThe TMVA weight  files are :" << endl;
+	INFO << "\t" << fTMVAWeightPath[0] << endl;
+	INFO << "\t" << fTMVAWeightPath[1] << endl;
+	INFO << "\t" << fTMVAWeightPath[2] << endl;
+	INFO << "\tThe Information file is :" << endl;
+	INFO << "\t" << fInformationFile << endl;
+	
+}
+//________________________________________________________________________
+EQM_MLP_Kinf::EQM_MLP_Kinf(string TMVAWeightPath,  int NumOfBatch, string InformationFile, double CriticalityThreshold):EquivalenceModel(new CLASSLogger("EQM_MLP_Kinf.log"))
+{
+	/**The information file and tmva weight*/
+	fTMVAWeightPath.push_back(TMVAWeightPath);
+	
+	/* INFORMATION FILE HANDLING */
+	
+	if(InformationFile == "")
+		InformationFile = StringLine::ReplaceAll(TMVAWeightPath,".xml",".nfo");
+	
+	fMaximalBU = 0;
+	fMaximalContent = 0;
+	fInformationFile = InformationFile;
+	LoadKeyword();
+	ReadNFO();//Getting information from fMLPInformationFile
+	
+	if( fMaximalBU == 0 || fMaximalContent == 0 )
+	{
+		ERROR<<"Can't find the k_maxfiscontent and/or k_maxburnup keyword(s) in .nfo file\n this is mandatory"<<endl;
+		exit(0);
+	}
+	
+	fNumberOfBatch = NumOfBatch;
+	fKThreshold = CriticalityThreshold ;
+	SetBurnUpPrecision(0.005);//1 % of the targeted burnup
+	SetPCMprecision(10);
+	SetBuildFuelFirstGuess(0.04);//First fissile content guess for the EquivalenceModel::BuildFuel algorithm
+	
+	INFO << "__An equivalence model has been define__" << endl;
+	INFO << "\tThis model is based on the prediction of kinf" << endl;
+	INFO << "\tThe TMVA (weight | information) files are :" << endl;
+	INFO << "\t" << "( " << fTMVAWeightPath[0] << " | " << fInformationFile << " )" << endl;
+	
+}
+//________________________________________________________________________
+EQM_MLP_Kinf::EQM_MLP_Kinf(CLASSLogger* log, string TMVAWeightPath,  int NumOfBatch, string InformationFile, double CriticalityThreshold):EquivalenceModel(log)
+{
+	
+	/**The information file and tmva weight*/
+	fTMVAWeightPath.push_back(TMVAWeightPath);
+	
+	if(InformationFile == "")
+		InformationFile = StringLine::ReplaceAll(TMVAWeightPath,".xml",".nfo");
+	
+	fMaximalBU = 0;
+	fMaximalContent = 0;
+	fInformationFile = InformationFile;
+	LoadKeyword();
+	ReadNFO();//Getting information from fMLPInformationFile
+	
+	if( fMaximalBU == 0 || fMaximalContent == 0 )
+	{
+		ERROR<<"Can't find the k_maxfiscontent and/or k_maxburnup keyword(s) in .nfo file\n this is mandatory"<<endl;
+		exit(0);
+	}
+	
+	fNumberOfBatch = NumOfBatch;
+	fKThreshold = CriticalityThreshold ;
+	SetBurnUpPrecision(0.005);//1 % of the targeted burnup
+	SetPCMprecision(10);
+	SetBuildFuelFirstGuess(0.04);//First fissile content guess for the EquivalenceModel::BuildFuel algorithm
+	
+	INFO << "__An equivalence model has been define__" << endl;
+	INFO << "\tThis model is based on the prediction of kinf" << endl;
+	INFO << "\tThe TMVA (weight | information) files are :" << endl;
+	INFO << "\t" << "( " << fTMVAWeightPath[0] << " | " << fInformationFile << " )" << endl;
+	
+	
+}
+
+//________________________________________________________________________
+void EQM_MLP_Kinf::LoadKeyword()
+{
+	DBGL
+	
+	fDKeyword.insert( pair<string, PWR_MLP_KINF_DMthPtr>( "k_zainame", &EQM_MLP_Kinf::ReadZAIName) );
+	fDKeyword.insert( pair<string, PWR_MLP_KINF_DMthPtr>( "k_maxburnup", &EQM_MLP_Kinf::ReadMaxBurnUp) );
+	fDKeyword.insert( pair<string, PWR_MLP_KINF_DMthPtr>( "k_maxfiscontent", &EQM_MLP_Kinf::ReadMaxFisContent) );
+	DBGL
+}
+
+
+//________________________________________________________________________
+void EQM_MLP_Kinf::ReadZAIName(const string &line)
+{
+	DBGL
+	int pos = 0;
+	string keyword = tlc(StringLine::NextWord(line, pos, ' '));
+	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());
+	
+	fFissileList.Add(Z, A, I, 1.0);
+	
+	DBGL
+}
+//________________________________________________________________________
+void EQM_MLP_Kinf::ReadMaxBurnUp(const string &line)
+{
+	DBGL
+	int pos = 0;
+	string keyword = tlc(StringLine::NextWord(line, pos, ' '));
+	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
+}
+//________________________________________________________________________
+void EQM_MLP_Kinf::ReadMaxFisContent(const string &line)
+{
+	DBGL
+	int pos = 0;
+	string keyword = tlc(StringLine::NextWord(line, pos, ' '));
+	if( keyword != "k_maxfiscontent" )	// Check the keyword
+	{
+		ERROR << " Bad keyword : \"k_maxfiscontent\" not found !" << endl;
+		exit(1);
+	}
+	
+	fMaximalContent = atof(StringLine::NextWord(line, pos, ' ').c_str());
+	
+	DBGL
+}
+//________________________________________________________________________
+void EQM_MLP_Kinf::ReadLine(string line)
+{
+	DBGL
+	
+	int pos = 0;
+	string keyword = tlc(StringLine::NextWord(line, pos, ' '));
+	
+	map<string, PWR_MLP_KINF_DMthPtr>::iterator it = fDKeyword.find(keyword);
+	
+	if(it != fDKeyword.end())
+		(this->*(it->second))( line );
+	
+	DBGL
+}
+
+
+
+//________________________________________________________________________
+TTree* EQM_MLP_Kinf::CreateTMVAInputTree(IsotopicVector TheFreshfuel, double ThisTime)
+{
+	/******Create Input data tree to be interpreted by TMVA::Reader***/
+	TTree*   InputTree = new TTree("InTMPKinf", "InTMPKinf");
+	
+	vector<float> 	InputTMVA;
+	for(int i = 0 ; i< (int)fMapOfTMVAVariableNames.size() ; i++)
+		InputTMVA.push_back(0);
+	
+	float Time = 0;
+	
+	IsotopicVector IVInputTMVA;
+	map<ZAI ,string >::iterator it;
+	int j = 0;
+	
+	for( it = fMapOfTMVAVariableNames.begin()  ; it != fMapOfTMVAVariableNames.end() ; it++)
+	{
+		InputTree->Branch( ((*it).second).c_str()	,&InputTMVA[j], ((*it).second + "/F").c_str());
+		IVInputTMVA+=  ((*it).first)*1;
+		j++;
+	}
+	
+	if(ThisTime != -1)
+		InputTree->Branch(	"Time"	,&Time	,"Time/F"	);
+	
+	IsotopicVector IVAccordingToUserInfoFile = TheFreshfuel.GetThisComposition(IVInputTMVA);
+	
+	double Ntot = IVAccordingToUserInfoFile.GetSumOfAll();
+	
+	IVAccordingToUserInfoFile = IVAccordingToUserInfoFile/Ntot;
+	
+	j = 0;
+	map<ZAI ,string >::iterator it2;
+	
+	for( it2 = fMapOfTMVAVariableNames.begin() ; it2 != fMapOfTMVAVariableNames.end() ; it2++)
+	{
+		InputTMVA[j] = IVAccordingToUserInfoFile.GetZAIIsotopicQuantity( (*it2).first ) ;
+		j++;
+	}
+	
+	Time = ThisTime;
+	
+	InputTree->Fill();
+	
+	return InputTree;
+	
+}
+//________________________________________________________________________
+double EQM_MLP_Kinf::ExecuteTMVA(TTree* InputTree,string WeightPath, bool IsTimeDependent)
+{
+	
+	// --- Create the Reader object
+	TMVA::Reader *reader = new TMVA::Reader( "Silent" );
+	
+	// Create a set of variables and declare them to the reader
+	// - the variable names MUST corresponds in name and type to those given in the weight file(s) used
+	vector<float> 	InputTMVA;
+	for(int i = 0 ; i< (int)fMapOfTMVAVariableNames.size() ; i++)
+		InputTMVA.push_back(0);
+	Float_t Time;
+	
+	map<ZAI ,string >::iterator it;
+	int j = 0;
+	for( it = fMapOfTMVAVariableNames.begin()  ; it != fMapOfTMVAVariableNames.end() ; it++)
+	{	reader->AddVariable( ( (*it).second ).c_str(),&InputTMVA[j]);
+		j++;
+	}
+	
+	if(IsTimeDependent)
+		reader->AddVariable( "Time" ,&Time);
+	
+	// --- Book the MVA methods
+	
+	// Book method MLP
+	TString methodName = "MLP method";
+	reader->BookMVA( methodName, WeightPath );
+	
+	map<ZAI ,string >::iterator it2;
+	j = 0;
+	for( it2 = fMapOfTMVAVariableNames.begin()  ; it2 != fMapOfTMVAVariableNames.end() ; it2++)
+	{
+		InputTree->SetBranchAddress(( (*it2).second ).c_str(),&InputTMVA[j]);
+		j++;
+	}
+	
+	if(IsTimeDependent)
+		InputTree->SetBranchAddress( "Time" ,&Time );
+	
+	InputTree->GetEntry(0);
+	Float_t val = (reader->EvaluateRegression( methodName ))[0];
+	
+	delete reader;
+	
+	return (double)val;//retourn k_{inf}(t = Time)
+}
+
+
+//________________________________________________________________________
+double EQM_MLP_Kinf::GetMaximumBurnUp_MLP(IsotopicVector TheFuel, double TargetBU)
+{
+	/**************************************************************************/
+	//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 TheFinalTime = BurnupToSecond(TargetBU);
+	double OldFinalTimeMinus = 0;
+	double MaximumBU = fMaximalBU;
+	double OldFinalTimePlus = BurnupToSecond(MaximumBU);
+	double k_av = 0; //average kinf
+	double OldPredictedk_av = 0;
+	for(int b = 0;b<fNumberOfBatch;b++)
+	{
+		float TheTime = (b+1)*TheFinalTime/fNumberOfBatch;
+		TTree* InputTree = CreateTMVAInputTree(TheFuel,TheTime);
+		OldPredictedk_av += ExecuteTMVA(InputTree,fTMVAWeightPath[0],true);
+		delete InputTree;
+	}
+	OldPredictedk_av/= fNumberOfBatch;
+	
+	//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-fKThreshold)  > 0) //The burnup can be increased
+		{
+			OldFinalTimeMinus = TheFinalTime;
+			TheFinalTime = TheFinalTime + fabs(OldFinalTimePlus - TheFinalTime)/2.;
+			
+			if(SecondToBurnup(TheFinalTime) >= (MaximumBU-MaximumBU*GetBurnUpPrecision() ) )
+				return MaximumBU;
+		}
+		
+		else if( (OldPredictedk_av-fKThreshold)  < 0)//The burnup is too high
+		{
+			OldFinalTimePlus = TheFinalTime;
+			TheFinalTime = TheFinalTime - fabs(OldFinalTimeMinus-TheFinalTime)/2.;
+			if( SecondToBurnup(TheFinalTime) < TargetBU*GetBurnUpPrecision() )
+				return 0;
+		}
+		
+		k_av = 0;
+		for(int b = 0;b<fNumberOfBatch;b++)
+		{
+			float TheTime = (b+1)*TheFinalTime/fNumberOfBatch;
+			TTree* InputTree = CreateTMVAInputTree(TheFuel,TheTime);
+			k_av += ExecuteTMVA(InputTree,fTMVAWeightPath[0],true);
+			delete InputTree;
+		}
+		k_av/= fNumberOfBatch;
+		
+		OldPredictedk_av = k_av;
+		count++;
+	}	while( fabs(OldPredictedk_av-fKThreshold) > GetPCMprecision() )  ;
+	
+	
+	return SecondToBurnup(TheFinalTime);
+}
+//________________________________________________________________________
+double EQM_MLP_Kinf::GetMaximumBurnUp_Pol2(IsotopicVector TheFuel,double TargetBU)
+{
+	
+	TTree* InputTree = CreateTMVAInputTree(TheFuel,-1);
+	double Alpha_0 = ExecuteTMVA(InputTree,fTMVAWeightPath[0],false);
+	double Alpha_1 = ExecuteTMVA(InputTree,fTMVAWeightPath[1],false);
+	double Alpha_2 = ExecuteTMVA(InputTree,fTMVAWeightPath[2],false);
+	delete InputTree;
+	
+	if(Alpha_0 < fKThreshold) //not enought fissile for sure !!
+		return 0;
+	
+	double Sum = 0;
+	double SumSquare = 0;
+	
+	for(int i = 1 ; i<=   fNumberOfBatch; i++)
+	{	Sum+= double(i)/double(fNumberOfBatch);
+		SumSquare+= double(i*i)/double(fNumberOfBatch*fNumberOfBatch);
+	}
+	
+	double C = (Alpha_0-fKThreshold)*fNumberOfBatch;
+	double B = Alpha_1*Sum;
+	double A = Alpha_2*SumSquare;
+	
+	double Delta = B*B-4*A*C;
+	
+	double T = 0;
+	if( Delta > 0)
+	{
+		double T_1 = (-B + sqrt(Delta))/(2*A);
+		double T_2 = (-B - sqrt(Delta))/(2*A);
+		
+		if(T_1 < 0 && T_2 > 0)
+			T = T_2;
+		else if(T_1 > 0 && T_2 < 0)
+			T = T_1;
+		
+		else if( T_1 > 0 && T_2 > 0 )
+		{
+			if(T_2 < T_1)
+				T = T_2;
+			else
+				T = T_1;
+		}
+		else
+		{
+			ERROR << "No positive solution" << endl;
+			exit(1);
+		}
+	}
+	else if(Delta == 0)
+	{	T = -B/(2*A);
+		if(T<0)
+		{	ERROR << "No positive solution" << endl;
+			exit(1);
+		}
+	}
+	else
+	{
+		WARNING << "No real solution" << endl;
+		double K_LongTime = Alpha_0+BurnupToSecond(10*TargetBU)*Alpha_1+Alpha_2*BurnupToSecond(10*TargetBU)*BurnupToSecond(10*TargetBU);
+		DBGV("K_LongTime " << K_LongTime)
+		
+		if(K_LongTime > fKThreshold)
+			return 10000;
+		else
+		{
+			ERROR << " CRITICAL ! Can't find a physical solution ! \n Should not happening please contact BLG :" << endl;
+			ERROR << "mail to baptiste.leniau@subatech.in2p2.fr\nor nicolas.thiolliere@subatech.in2p3.fr " << endl;
+			exit(1);
+		}
+		
+	}
+	
+	return SecondToBurnup(T);
+}
+//________________________________________________________________________
+double EQM_MLP_Kinf::GetMaximumBurnUp(IsotopicVector TheFuel, double TargetBU)
+{
+	double TheBurnUp = -1;
+	if(fTMVAWeightPath.size() == 1)
+		TheBurnUp  = 	GetMaximumBurnUp_MLP(TheFuel,TargetBU);
+	else if(fTMVAWeightPath.size() == 3)
+		TheBurnUp = GetMaximumBurnUp_Pol2(TheFuel,TargetBU);
+	
+	else
+	{
+		ERROR << "This method is not yet set up" << endl;
+		exit(0);
+	}
+	
+	return TheBurnUp;
+}
+//________________________________________________________________________
+double EQM_MLP_Kinf::GetFissileMolarFraction(vector <IsotopicVector> IVStream, double TargetBU)
+{
+	
+	//initialization
+	IsotopicVector Fissil = IVStream[0];
+	IsotopicVector Fertil = IVStream[1];
+
+	double FissileContent = GetActualFissileContent();
+	double OldFissileContentMinus = 0;
+	double OldFissileContentPlus = fMaximalContent;
+	double PredictedBU = 0 ;
+	IsotopicVector FreshFuel = (1-FissileContent)*(Fertil/Fertil.GetSumOfAll()) + FissileContent*(Fissil/Fissil.GetSumOfAll());
+	double OldPredictedBU = GetMaximumBurnUp(FreshFuel,TargetBU);
+	
+	double Precision = GetBurnUpPrecision()*TargetBU; //1 % of the targeted burnup
+	int count = 0;
+	int MaximumLoopCount = 500;
+	do
+	{
+		if(count > MaximumLoopCount )
+		{
+			ERROR << "CRITICAL ! Can't manage to predict fissile content\nHint : Try to decrease the precision on burnup using :\nYourEQM_MLP_Kinf->SetBurnUpPrecision(prop); with prop the precision  (default 0.5percent :  0.005) INCREASE IT\nIf this message still appear mail to leniau@subatech.in2p3.fr\nor nicolas.thiolliere@subatech.in2p3.fr " << endl;
+			exit(1);
+		}
+		
+		if( (OldPredictedBU - TargetBU) < 0 ) //The Content can be increased
+		{
+			OldFissileContentMinus = FissileContent;
+			FissileContent = FissileContent + fabs(OldFissileContentPlus-FissileContent)/2.;
+		}
+		else if( (OldPredictedBU - TargetBU) > 0) //The Content is too high
+		{
+			OldFissileContentPlus = FissileContent;
+			FissileContent = FissileContent - fabs(OldFissileContentMinus-FissileContent)/2.;
+		}
+		
+		IsotopicVector FreshFuel = (1-FissileContent)*(Fertil/Fertil.GetSumOfAll()) + FissileContent*(Fissil/Fissil.GetSumOfAll());
+		PredictedBU = GetMaximumBurnUp(FreshFuel,TargetBU);
+		
+		OldPredictedBU = PredictedBU;
+		count ++;
+		
+	}while(fabs(TargetBU-PredictedBU)>Precision);
+	
+	
+	
+	DBGV("Predicted BU " << PredictedBU << " FissileContent " << FissileContent);
+	return FissileContent;
+}
+//________________________________________________________________________

source/branches/CLASSV4.2/Model/Equivalence/EQM_MLP_Kinf.hxx

+#ifndef _EQM_MLP_Kinf_HXX
+#define _EQM_MLP_Kinf_HXX
+
+#include "EquivalenceModel.hxx"
+#include "TTree.h"
+#include <map>
+
+using namespace std;
+
+class EQM_MLP_Kinf;
+#ifndef __CINT__
+typedef void (EQM_MLP_Kinf::*PWR_MLP_KINF_DMthPtr)( const string & ) ;
+#endif
+
+//-----------------------------------------------------------------------------//
+//! Defines an EquivalenceModel based on neural network to predict @f$k_{\infty}@f$.
+
+/*!
+The aim of these class is to constuct a fuel from an equivalence model
+based on a  Multi layer perceptron (MLP)
+This MLP aims to predict the @f$k_{\infty}(t)@f$ of a PWR-MOX from a given fresh fuel 
+composition
+With this MLP prediction and a given number of batch (for the loading plan) an 
+average @f$\langle k_{\infty}\rangle (t)@f$ is calculated according :
+
+@f$\langle k_{\infty}\rangle ^{batch}(t) = \frac{1}{N}\sum_{i}^{N}k_{\infty}(t+\frac{iT}{N})@f$
+The maximal reachable burnup has to verify the following conditions :
+@f$\langle k_{\infty}\rangle ^{batch}(T/N) = <\langle k_{\infty}\rangle ^{batch}(2T/N)  = ... = k_{Threshold}@f$
+Where @f$k_{Threshold}@f$ is the criticality threshold which take into account leakage and capture
+in non simulated devices such as control rods and mixing grid.
+
+ @author BLG
+ @version 1.0
+ */
+//________________________________________________________________________
+
+
+class EQM_MLP_Kinf : public EquivalenceModel
+{
+	public:
+	/*!
+	 \name Constructor
+	 */
+	//@{
+
+	//{
+	/// Polynnomial 2nd order constructor @f$k_{\infty} = \alpha_{0} + \alpha_{1}t + \alpha_{2}t^{2}@f$
+ 	/// @f$\alpha_{0}@f$, @f$\alpha_{1}@f$, @f$\alpha_{2}@f$ are predict by 3 MLP (one for each)
+ 	/*!
+	 Create a EQM_MLP_Kinf 
+	 \param  TMVAWeightPath0 :  PAth to the .xml file containing neural network informations for @f$\alpha_{0}@f$ prediction : PATH/TMVAWeight.xml (total path to tmva weight)
+	 \param  TMVAWeightPath1 :  PAth to the .xml file containing neural network informations for @f$\alpha_{1}@f$ prediction: PATH/TMVAWeight.xml (total path to tmva weight)
+	 \param  TMVAWeightPath2 :  PAth to the .xml file containing neural network informations for @f$\alpha_{2}@f$ prediction: PATH/TMVAWeight.xml (total path to tmva weight)
+	 \param  InformationFile : Total path to the file containing time steps, fissile and ferile list (ante and post fabrication time cooling). Default is the same total path as TMVAWeightPath but extension is replaced by .nfo
+	 \param  NumOfBatch : Number of batch for the loading plan (often 3 or 4 for PWR)
+	 \param  CriticalityThreshold : Threshold for the @f$k_{\infty}@f$ (see detailed description)
+	 */
+	EQM_MLP_Kinf(string WeightPathAlpha0, string WeightPathAlpha1, string WeightPathAlpha2, string InformationFile,  int NumOfBatch, double CriticalityThreshold = 1.01);
+	//}
+	//{
+	/// Polynnomial 2nd order constructor @f$k_{\infty} = \alpha_{0} + \alpha_{1}t + \alpha_{2}t^{2}@f$
+ 	/// @f$\alpha_{0}@f$, @f$\alpha_{1}@f$, @f$\alpha_{2}@f$ are predict by 3 MLP (one for each)
+	 /*!
+	 Create a EQM_MLP_Kinf 
+	 \param log : use for log
+	 \param  TMVAWeightPath0 :  PAth to the .xml file containing neural network informations for @f$\alpha_{0}@f$ prediction : PATH/TMVAWeight.xml (total path to tmva weight)
+	 \param  TMVAWeightPath1 :  PAth to the .xml file containing neural network informations for @f$\alpha_{1}@f$ prediction: PATH/TMVAWeight.xml (total path to tmva weight)
+	 \param  TMVAWeightPath2 :  PAth to the .xml file containing neural network informations for @f$\alpha_{2}@f$ prediction: PATH/TMVAWeight.xml (total path to tmva weight)
+	 \param  InformationFile : Total path to the file containing time steps, fissile and ferile list (ante and post fabrication time cooling). Default is the same total path as TMVAWeightPath but extension is replaced by .nfo
+	 \param  NumOfBatch : Number of batch for the loading plan (often 3 or 4 for PWR)
+	 \param  CriticalityThreshold : Threshold for the @f$k_{\infty}@f$ (see detailed description)
+	 */
+	EQM_MLP_Kinf(CLASSLogger* log, string WeightPathAlpha0, string WeightPathAlpha1, string WeightPathAlpha2, string InformationFile,  int NumOfBatch, double CriticalityThreshold = 1.01 );
+	//}
+	//{
+	/// Neural network predictor. The kinf(t) is predicted with a MLP 
+	/*!
+	 Create a EQM_MLP_Kinf 
+	 \param  TMVAWeightPath :  PAth to the .xml file containing neural network informations : PATH/TMVAWeight.xml (total path to tmva weight)
+	 \param  NumOfBatch : Number of batch for the loading plan (often 3 or 4 for PWR)
+	 \param  InformationFile : Total path to the file containing time steps, fissile and ferile list (ante and post fabrication time cooling). Default is the same total path as TMVAWeightPath but extension is replaced by .nfo
+	 \param  CriticalityThreshold : Threshold for the @f$k_{\infty}@f$ (see detailed description)
+	 */
+	EQM_MLP_Kinf(string TMVAWeightPath,int NumOfBatch, string InformationFile = "", double CriticalityThreshold = 1.01);
+	//}
+	
+	//{
+	/// Neural network predictor. The kinf(t) is predicted with a MLP 
+	/*!
+	 Create a EQM_MLP_Kinf
+	 \param log : use for log
+	 \param  TMVAWeightPath :  PAth to the .xml file containing neural network informations : PATH/TMVAWeight.xml (total path to tmva weight)
+	 \param  NumOfBatch : Number of batch for the loading plan (often 3 or 4 for PWR)
+	 \param  InformationFile : Total path to the file containing time steps, fissile and ferile list (ante and post fabrication time cooling). Default is the same total path as TMVAWeightPath but extension is replaced by .nfo
+	 \param  CriticalityThreshold : Threshold for the @f$k_{\infty}@f$ (see detailed description)
+	 */
+	EQM_MLP_Kinf(CLASSLogger* log, string TMVAWeightPath,int NumOfBatch, string InformationFile = "", double CriticalityThreshold = 1.01);
+	//}
+	//@}
+	
+	//{
+	/// Return the molar fissile fraction according fissile & ferile content 
+	/*!
+	 \param Fissil : The composition of the fissile matter
+	 \param Fertil : The composition of the Fertil matter
+	 \param BurnUp : Maximum achievable burn up envisaged
+	 */
+	virtual double GetFissileMolarFraction(vector <IsotopicVector> IVStream, double BurnUp);
+	//}
+	
+	/*!
+	 \name Get/Set methods
+	 */
+	//@{
+
+	void SetBurnUpPrecision(double prop){fBurnUpPrecision = prop;} //!< Set the precision on Burnup : proportion of the targeted burnup
+	void SetPCMprecision(double pcm){fPCMprecision = pcm;}		  //!< Set the precision on @f$\langle k \rangle@f$ prediction [pcm]. Neural network predictor constructors
+
+	double GetBurnUpPrecision(){return fBurnUpPrecision;}//!< Get the precision on Burnup : proportion of the targeted burnup
+	double GetPCMprecision(){return fPCMprecision/1e5;}//!< Get the precision on @f$\langle k \rangle@f$ prediction []. Neural network predictor constructors
+	double GetMaximumBurnUp_MLP(IsotopicVector TheFuel, double TargetBU);
+
+	//@}
+	TTree* CreateTMVAInputTree(IsotopicVector FreshFuel, double ThisTime);//!<Create input tmva tree to be read by ExecuteTMVA
+	double ExecuteTMVA(TTree* theTree, string WeightPath, bool IsTimeDependant);//!<Execute the MLP according to the input tree created by CreateTMVAInputTree
+	double SecondToBurnup(double Second){return Second*fSpecificPower/(24*3.6e6);}
+
+	
+	/*!
+	 \name Reading NFO related Method
+	 */
+	//@{
+	
+	//{
+	/// LoadKeyword() : make the correspondance between keyword and reading method
+	void LoadKeyword();
+	//}
+	
+	//{
+	/// 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 ReadZAIName(const string &line);
+	//}
+	
+	//{
+	/// ReadMaxBurnUp : read a guessed (very overestimated) maximum burnup a fuel can reach (purpose : algorithm initialization)
+	/*!
+	 \param line : line suppossed to contain the ZAI name  starts with "k_zainame" keyword
+	 */
+	void ReadMaxBurnUp(const string &line);
+	//}
+	
+	//{
+	/// ReadMaxFisContent  : read a guessed (very overestimated) maximum fissile content (purpose : algorithm initialization)
+	/*!
+	 \param line : line suppossed to contain the ZAI name  starts with "k_zainame" keyword
+	 */
+	void ReadMaxFisContent(const string &line);
+	//}
+	
+	//{
+	/// ReadLine : read a line
+	/*!
+	 \param line : line to read
+	 */
+	void ReadLine(string line);
+	//}
+	
+	//@}
+	
+	private :
+	vector <string> fTMVAWeightPath;		//!<The weight needed by TMVA to construct and execute the multilayer perceptron
+
+#ifndef __CINT__
+	map<string, PWR_MLP_KINF_DMthPtr> fDKeyword;
+#endif
+	
+	map<ZAI,string> fMapOfTMVAVariableNames;//!<  List of TMVA input variable names (read from fMLPInformationFile ) , name depends on the training step
+
+	
+	
+	double GetMaximumBurnUp_Pol2(IsotopicVector TheFuel, double TargetBU);
+	double GetMaximumBurnUp(IsotopicVector TheFuel, double TargetBU);//!<Get the maximum reachable burnup according the freshfuel composition
+	void   GetModelInformation();//!<Read the fMLPInformationFile and fill containers and variables
+	
+	double BurnupToSecond(double BurnUp){return BurnUp/fSpecificPower*(24*3.6e6);}
+
+	int 	fNumberOfBatch;		//!< The number of batches for the loading plan
+	double 	fKThreshold;		//!< The @f$k_{Threshold}@f$
+	double 	fSpecificPower; 	//!< The specific power in W/gHM (HM: heavy Metal)
+	double 	fMaximalBU;			//!< The approx. maximum burnup reachable by the MLP model		
+	double  fMaximalContent;	//!< The approx. maximum fissile content reachable by the MLP model
+	double 	fBurnUpPrecision;	//!< precision on Burnup 
+	double 	fPCMprecision;		//!< precision on @f$\langle k \rangle@f$ prediction [pcm]
+			
+
+};
+
+#endif
+

source/branches/CLASSV4.2/Model/Equivalence/EQM_PWR_LIN_MOX.cxx

+#include "EQM_PWR_LIN_MOX.hxx"
+
+#include "CLASSConstante.hxx"
+
+#include <vector>
+
+#include "StringLine.hxx"
+#include "CLASSLogger.hxx"
+#include "IsotopicVector.hxx"
+
+
+
+EQM_PWR_LIN_MOX::EQM_PWR_LIN_MOX(string WeightPath):EquivalenceModel(new CLASSLogger("EQM_PWR_LIN_MOX.log"))
+{
+	fWeightPath = WeightPath;
+
+	ifstream DataDB(fWeightPath.c_str());							// Open the File
+	if(!DataDB)
+		WARNING << "Can't open \"" << fWeightPath << "\"\n" << endl;
+
+	string line;
+	int start = 0;	// First Get Fuel Parameter
+	getline(DataDB, line);
+
+	if( StringLine::NextWord(line, start, ' ') != "PARAM")
+	{
+		ERROR << " Bad Database file : " <<  fWeightPath << " Can't find the Parameter of the DataBase " << endl;
+		exit (1);
+	}
+	while(start < (int)line.size())
+		fFuelParameter.push_back(atof(StringLine::NextWord(line, start, ' ').c_str()));
+
+	INFO << " " << (int)fFuelParameter.size() << " parameters have been read " << endl;
+
+
+
+	//-----------------------------------------------------------------------------//
+	//-----------------------------------------------------------------------------//
+	//-----------------------------------------------------------------------------//
+	//-----------------------------------------------------------------------------//
+	//-----------------------------------------------------------------------------//
+	// ADD ENrichment of the U reading !!!!!!!!!!!!!!!!!!!!!!!!!!!!		       //
+	//-----------------------------------------------------------------------------//
+	//-----------------------------------------------------------------------------//
+	//-----------------------------------------------------------------------------//
+	//-----------------------------------------------------------------------------//
+
+	ZAI U8(92,238,0);
+	ZAI U5(92,235,0);
+	double U5_enrich = 0.0025;
+	fFertileList = U5*U5_enrich + U8*(1-U5_enrich);
+
+
+	ZAI Pu8(94,238,0);
+	ZAI Pu9(94,239,0);
+	ZAI Pu0(94,240,0);
+	ZAI Pu1(94,241,0);
+	ZAI Pu2(94,242,0);
+	fFissileList = Pu8*1+Pu9*1+Pu0*1+Pu1*1+Pu2*1;
+
+
+}
+
+
+EQM_PWR_LIN_MOX::EQM_PWR_LIN_MOX(CLASSLogger* log, string WeightPath):EquivalenceModel(log)
+{
+	fWeightPath = WeightPath;
+
+	ifstream DataDB(fWeightPath.c_str());							// Open the File
+	if(!DataDB)
+		WARNING << " Can't open \"" << fWeightPath << "\"\n" << endl;
+
+	string line;
+	int start = 0;	// First Get Fuel Parameter
+	getline(DataDB, line);
+
+	if( StringLine::NextWord(line, start, ' ') != "PARAM")
+	{
+		ERROR << " Bad Database file : " <<  fWeightPath << " Can't find the Parameter of the DataBase" <<  endl;
+		exit (1);
+	}
+	while(start < (int)line.size())
+		fFuelParameter.push_back(atof(StringLine::NextWord(line, start, ' ').c_str()));
+
+	INFO << fFuelParameter.size() << " have been read" <<  endl;
+
+
+
+	//-----------------------------------------------------------------------------//
+	//-----------------------------------------------------------------------------//
+	//-----------------------------------------------------------------------------//
+	//-----------------------------------------------------------------------------//
+	//-----------------------------------------------------------------------------//
+	// ADD ENrichment of the U reading !!!!!!!!!!!!!!!!!!!!!!!!!!!!		       //
+	//-----------------------------------------------------------------------------//
+	//-----------------------------------------------------------------------------//
+	//-----------------------------------------------------------------------------//
+	//-----------------------------------------------------------------------------//
+
+	ZAI U8(92,238,0);
+	ZAI U5(92,235,0);
+	double U5_enrich = 0.0025;
+	fFertileList = U5*U5_enrich + U8*(1-U5_enrich);
+
+
+	ZAI Pu8(94,238,0);
+	ZAI Pu9(94,239,0);
+	ZAI Pu0(94,240,0);
+	ZAI Pu1(94,241,0);
+	ZAI Pu2(94,242,0);
+	fFissileList = Pu8*1+Pu9*1+Pu0*1+Pu1*1+Pu2*1;
+
+
+}
+
+EQM_PWR_LIN_MOX::~EQM_PWR_LIN_MOX()
+{
+
+}
+
+//________________________________________________________________________
+vector<double> EQM_PWR_LIN_MOX::BuildFuel(double BurnUp, double HMMass,vector <IsotopicVector> IVStream)
+{
+
+	//-----------------------------------------------------------------------------//
+	//-----------------------------------------------------------------------------//
+	//-----------------------------------------------------------------------------//
+	//-----------------------------------------------------------------------------//
+	//-----------------------------------------------------------------------------//
+	// ADD ENrichment of the U check ++ Check Un seul fertile !!!!		       //
+	//-----------------------------------------------------------------------------//
+	//-----------------------------------------------------------------------------//
+	//-----------------------------------------------------------------------------//
+	//-----------------------------------------------------------------------------//
+
+	vector<IsotopicVector> FissilArray = StreamArray[0];
+	vector<IsotopicVector> FertilArray = StreamArray[1];
+
+	vector<double> lambda;
+	for(int i = 0; i < (int) (FissilArray.size() + FertilArray.size()); i++)
+		lambda.push_back(0);
+
+
+	double Na = 6.02214129e23;	//N Avogadro
+
+	IsotopicVector FullUsedStock;
+	IsotopicVector stock;
+
+	bool FuelBuild = false;
+	if(FissilArray.size() ==  0)
+	{
+		for(int i = 0; i < (int)lambda.size(); i++)
+			lambda[i] = -1;
+
+		FuelBuild = true;
+	}
+	int N_FissilStock_OnCheck = 0;
+
+	while(!FuelBuild)
+	{
+
+		double nPu_0 = 0;
+		double MPu_0 = 0;
+		{
+			map<ZAI ,double>::iterator it;
+
+			map<ZAI ,double> isotopicquantity = FullUsedStock.GetSpeciesComposition(94).GetIsotopicQuantity();
+			for( it = isotopicquantity.begin(); it != isotopicquantity.end(); it++ )
+				nPu_0 += (*it).second;
+
+			IsotopicVector IV_Pm_Am = (FullUsedStock.GetSpeciesComposition(94)
+						   + ZAI(94,241,0)*FullUsedStock.GetZAIIsotopicQuantity(95,241,0));
+			//Add the 241Am as 241Pu... the Pu is not old in the Eq Model but is in the FissileArray....;
+			MPu_0 += cZAIMass.GetMass(IV_Pm_Am);
+		}
+		stock = FissilArray[N_FissilStock_OnCheck];
+		double nPu_1 = 0;
+		double MPu_1 = 0;
+		double Sum_AlphaI_nPuI = 0;
+		double Sum_AlphaI_nPuI0 = 0;
+		{
+			map<ZAI ,double>::iterator it;
+			map<ZAI ,double> isotopicquantity = stock.GetSpeciesComposition(94).GetIsotopicQuantity();
+
+			for( it = isotopicquantity.begin(); it != isotopicquantity.end(); it++ )
+			{
+				if ((*it).first.A() >= 238 && (*it).first.A() <= 242)
+				{
+					nPu_1 += (*it).second;
+					Sum_AlphaI_nPuI += fFuelParameter[(*it).first.A() -237]*(*it).second;
+				}
+			}
+
+			isotopicquantity = (stock.GetSpeciesComposition(94) + ZAI(94,241,0)*stock.GetZAIIsotopicQuantity(95,241,0)).GetIsotopicQuantity();
+			IsotopicVector IV_Pm_Am = (stock.GetSpeciesComposition(94)
+						   + ZAI(94,241,0)*stock.GetZAIIsotopicQuantity(95,241,0));
+			//Add the 241Am as 241Pu... the Pu is not old in the Eq Model but is in the FissileArray....
+			MPu_1 += cZAIMass.GetMass(IV_Pm_Am);
+
+			isotopicquantity = FullUsedStock.GetSpeciesComposition(94).GetIsotopicQuantity();
+			for( it = isotopicquantity.begin(); it != isotopicquantity.end(); it++ )
+				if ((*it).first.A() >= 238 && (*it).first.A() <= 242)
+				{
+					Sum_AlphaI_nPuI0 += fFuelParameter[(*it).first.A() -237]*(*it).second;
+				}
+		}
+
+		double StockFactionToUse = 0;
+
+		double NT = HMMass*1e6 * Na / (cZAIMass.GetMass( ZAI(92,238,0) ) * 0.997
+					       + cZAIMass.GetMass( ZAI(92,235,0) ) * 0.003 );
+
+		double N1 = (BurnUp - fFuelParameter[6]) * NT;
+		double N2 = -Sum_AlphaI_nPuI0;
+		double N3 = -fFuelParameter[0] * Na / (cZAIMass.GetMass( ZAI(92,238,0) )*0.997
+						       + cZAIMass.GetMass( ZAI(92,235,0) )*0.003 )
+							* (HMMass*1e6 - MPu_0*1e6);
+
+		double D1 = Sum_AlphaI_nPuI;
+		double D2 = -fFuelParameter[0] * MPu_1*1e6 * Na / (cZAIMass.GetMass( ZAI(92,238,0) )*0.997
+								   + cZAIMass.GetMass( ZAI(92,235,0) )*0.003 );
+
+		StockFactionToUse = (N1 + N2 + N3) / (D1 + D2);
+
+		if(StockFactionToUse < 0)
+		{
+			WARNING << "!!!FabricationPlant!!! Oups Bug in calculating stock fraction to use " <<  endl;
+			lambda[N_FissilStock_OnCheck] = 0.;
+			N_FissilStock_OnCheck++;
+			FuelBuild = false;
+		}
+		else if( StockFactionToUse > 1 )
+		{
+
+			FullUsedStock += stock;
+			lambda[N_FissilStock_OnCheck] = 1;
+			N_FissilStock_OnCheck++;
+			FuelBuild = false;
+		}
+		else
+		{
+			lambda[N_FissilStock_OnCheck] = StockFactionToUse;
+
+			FuelBuild = true;
+
+			double U8_Quantity = (HMMass - (MPu_0+StockFactionToUse*MPu_1 ))/(cZAIMass.GetMass( ZAI(92,238,0))*0.997 + cZAIMass.GetMass( ZAI(92,235,0))*0.003 )*Na/1e-6;
+
+			lambda.back() = U8_Quantity / FertilArray[0].GetSumOfAll();
+		}
+
+
+		if( N_FissilStock_OnCheck ==  (int) FissilArray.size() )	// Check if the last Fissil stock has been tested... quit if so...
+		{
+			for(int i = 0; i < (int)lambda.size(); i++)
+				lambda[i] = -1;
+			
+			FuelBuild = true;
+		}
+	}
+	
+	
+	
+	return lambda;
+}

source/branches/CLASSV4.2/Model/Equivalence/EQM_PWR_LIN_MOX.hxx

+#ifndef _EQM_PWR_LIN_MOX_HXX
+#define _EQM_PWR_LIN_MOX_HXX
+
+#include "EquivalenceModel.hxx"
+
+#include <string>
+
+using namespace std;
+
+//-----------------------------------------------------------------------------//
+//! Defines an EquivalenceModel based on a linear fit
+
+/*!
+ The aim of these class is to constuct a fuel from an equivalence model
+ based on a Linear Eq Model @f$BU = \alpha_{0} + \sum_{i\in fissile}\alpha_{i}\cdot n_{i} @f$
+ For one set of  @f$\alpha @f$ values the fabrication time is fixed in order to decrease 
+ the dimentionality by 1 (@f$^{241}Am@f$ is thus fixed by this fixed decay time) .
+ @author BaM
+ @version 3.0
+ */
+//________________________________________________________________________
+
+class EQM_PWR_LIN_MOX : public EquivalenceModel
+{
+	public :
+	/*!
+	 \name Constructor
+	 */
+	//@{
+	
+	//{
+	/// Simple constructor
+	
+	/*!
+	 Make a new EQM_PWR_LIN_MOX
+	 \param WeightPath : Path to the file containing the @f$\alpha_{i}@f$. The file is format as :
+	 
+	 @f$PARAM@f$   @f$\alpha_{0}@f$  @f$\alpha_{^{238}Pu}@f$   @f$\alpha_{^{239}Pu}@f$  @f$\alpha_{^{240}Pu}@f$   @f$\alpha_{^{241}Pu}@f$   @f$\alpha_{^{242}Pu}@f$
+	 
+	 */
+	EQM_PWR_LIN_MOX(string WeightPath);
+	//}
+	
+	//{
+	/// Logger constructor
+	
+	/*!
+	 Make a new EQM_PWR_LIN_MOX
+	 \param log : use for the log
+	 \param WeightPath : Path to the file containing the @f$\alpha_{i}@f$. The file is format as :
+	 
+	  @f$PARAM@f$   @f$\alpha_{0}@f$  @f$\alpha_{^{238}Pu}@f$   @f$\alpha_{^{239}Pu}@f$  @f$\alpha_{^{240}Pu}@f$   @f$\alpha_{^{241}Pu}@f$   @f$\alpha_{^{242}Pu}@f$
+	 
+	 */
+	EQM_PWR_LIN_MOX(CLASSLogger* log, string WeightPath);
+	//}
+	
+	~EQM_PWR_LIN_MOX();
+	//@}
+
+	virtual vector<double> BuildFuel(double BurnUp, double HMMass, vector < vector <IsotopicVector> > StreamArray);
+
+	private :
+
+	string fWeightPath;		//!< The full path to the file containing the @f$\alpha_{i}@f$
+	vector<double> fFuelParameter;	//!< The vector of @f$\alpha_{i}@f$
+
+};
+
+#endif
+

source/branches/CLASSV4.2/Model/Equivalence/EQM_PWR_MLP_MOX.cxx

+#include "EquivalenceModel.hxx"
+#include "EQM_PWR_MLP_MOX.hxx"
+#include "CLASSLogger.hxx"
+#include "StringLine.hxx"
+
+#include <string>
+#include <iostream>
+#include <fstream>
+#include <algorithm>
+#include <cmath>
+#include <cassert>
+
+#include "TSystem.h"
+#include "TMVA/Reader.h"
+#include "TMVA/Tools.h"
+#include "TMVA/MethodCuts.h"
+
+
+//________________________________________________________________________
+//
+//		EQM_PWR_MLP_MOX
+//
+//	Equivalenve Model based on multi layer perceptron from TMVA (root cern)
+//	For REP MOX use
+//
+//________________________________________________________________________
+
+EQM_PWR_MLP_MOX::EQM_PWR_MLP_MOX(string TMVAWeightPath):EquivalenceModel(new CLASSLogger("EQM_PWR_MLP_MOX.log"))
+{
+	fTMVAWeightPath = TMVAWeightPath;
+
+	ZAI U8(92,238,0);
+	ZAI U5(92,235,0);
+	double U5_enrich = 0.0025;
+
+	ZAI Pu8(94,238,0);
+	ZAI Pu9(94,239,0);
+	ZAI Pu0(94,240,0);
+	ZAI Pu1(94,241,0);
+	ZAI Pu2(94,242,0);
+
+	fStreamList["Fissile"] = Pu8*1+Pu9*1+Pu0*1+Pu1*1+Pu2*1;
+	fStreamList["Fertile"] = U5*U5_enrich + U8*(1-U5_enrich);
+	
+	fFirstGuessContent["Fissile"] = 0.02;
+	fFirstGuessContent["Fertile"] = 1 - fFirstGuessContent["Fissile"];
+
+	INFO << "__An equivalence model of PWR MOX has been define__" << endl;
+	INFO << "\tThis model is based on a multi layer perceptron" << endl;
+	INFO << "\t\tThe TMVA weight file is :" << endl;
+	INFO << "\t\t\t" << fTMVAWeightPath << endl;
+
+}
+
+//________________________________________________________________________
+EQM_PWR_MLP_MOX::EQM_PWR_MLP_MOX(CLASSLogger* log, string TMVAWeightPath):EquivalenceModel(log)
+{
+	fTMVAWeightPath = TMVAWeightPath;
+
+	ZAI U8(92,238,0);
+	ZAI U5(92,235,0);
+	double U5_enrich = 0.0025;
+
+	ZAI Pu8(94,238,0);
+	ZAI Pu9(94,239,0);
+	ZAI Pu0(94,240,0);
+	ZAI Pu1(94,241,0);
+	ZAI Pu2(94,242,0);
+
+	fStreamList["Fissile"] = Pu8*1+Pu9*1+Pu0*1+Pu1*1+Pu2*1;
+	fStreamList["Fertile"] = U5*U5_enrich + U8*(1-U5_enrich);
+	
+	fFirstGuessContent["Fissile"] = 0.02;
+	fFirstGuessContent["Fertile"] = 1 - fFirstGuessContent["Fissile"];
+
+	INFO << "__An equivalence model of PWR MOX has been define__" << endl;
+	INFO << "\tThis model is based on a multi layer perceptron" << endl;
+	INFO << "\t\tThe TMVA weight file is :" << endl;
+	INFO << "\t\t\t" << fTMVAWeightPath << endl;
+
+}
+
+//________________________________________________________________________
+TTree* EQM_PWR_MLP_MOX::CreateTMVAInputTree(map < string , IsotopicVector> IVStream, double BurnUp)
+{
+	TTree*   InputTree = new TTree("EQTMP", "EQTMP");
+	float Pu8   		 = 0;
+	float Pu9   		 = 0;
+	float Pu10  		 = 0;
+	float Pu11  		 = 0;
+	float Pu12  		 = 0;
+	float Am1   		 = 0;
+	float U5_enrichment  	 = 0;
+	float BU  		 = 0;
+
+	InputTree->Branch(	"Pu8"			,&Pu8			,"Pu8/F"		);
+	InputTree->Branch(	"Pu9"			,&Pu9			,"Pu9/F"		);
+	InputTree->Branch(	"Pu10"			,&Pu10			,"Pu10/F"		);
+	InputTree->Branch(	"Pu11"			,&Pu11			,"Pu11/F"		);
+	InputTree->Branch(	"Pu12"			,&Pu12			,"Pu12/F"		);
+	InputTree->Branch(	"Am1"			,&Am1			,"Am1/F"		);
+	InputTree->Branch(	"U5_enrichment"	,&U5_enrichment	,"U5_enrichment/F"	);
+	InputTree->Branch(	"BU"			,&BU			,"BU/F"			);
+
+
+	float U8     = IVStream["Fertile"].GetZAIIsotopicQuantity(92,238,0);
+	float U5     = IVStream["Fertile"].GetZAIIsotopicQuantity(92,235,0);
+	float U4     = IVStream["Fertile"].GetZAIIsotopicQuantity(92,234,0);
+
+	float UTOT = U8 + U5 + U4;
+
+	Pu8    	= IVStream["Fissile"].GetZAIIsotopicQuantity(94,238,0);
+	Pu9    	= IVStream["Fissile"].GetZAIIsotopicQuantity(94,239,0);
+	Pu10   	= IVStream["Fissile"].GetZAIIsotopicQuantity(94,240,0);
+	Pu11   	= IVStream["Fissile"].GetZAIIsotopicQuantity(94,241,0);
+	Pu12   	= IVStream["Fissile"].GetZAIIsotopicQuantity(94,242,0);
+	Am1   	= IVStream["Fissile"].GetZAIIsotopicQuantity(95,241,0);
+
+	double TOTPU = (Pu8+Pu9+Pu10+Pu11+Pu12+Am1);
+
+	Pu8 	= Pu8	 / TOTPU ;
+	Pu9 	= Pu9  	 / TOTPU ;
+	Pu10 	= Pu10	 / TOTPU ;
+	Pu11 	= Pu11	 / TOTPU ;
+	Pu12	= Pu12	 / TOTPU ;
+	Am1 	= Am1  / TOTPU ;
+
+	U5_enrichment = U5 / UTOT;
+
+	BU = BurnUp;
+	if(Pu8 + Pu9 + Pu10 + Pu11 + Pu12 + Am1 > 1.00001 )//?????1.00001??? I don't know it! goes in condition if  = 1 !! may be float/double issue ...
+	{
+		ERROR << Pu8 << " " << Pu9 << " " << Pu10 << " " << Pu11 << " " << Pu12 << " " << Am1 << endl;
+		exit(0);
+	}
+	// All value are molar (!weight)
+
+	InputTree->Fill();
+	return InputTree;
+}
+//________________________________________________________________________
+map < string , double> EQM_PWR_MLP_MOX::ExecuteTMVA(TTree* theTree)
+{
+	map < string , double> MolarFraction;
+	// --- Create the Reader object
+	TMVA::Reader *reader = new TMVA::Reader( "Silent" );
+	// Create a set of variables and declare them to the reader
+	// - the variable names MUST corresponds in name and type to those given in the weight file(s) used
+	Float_t Pu8,Pu9,Pu10,Pu11,Pu12,Am1,BU,U5_enrichment;
+
+	reader->AddVariable( "BU"   			,&BU 			);
+	reader->AddVariable( "U5_enrichment"	,&U5_enrichment 	);
+	reader->AddVariable( "Pu8"  			,&Pu8 			);
+	reader->AddVariable( "Pu9"  			,&Pu9 			);
+	reader->AddVariable( "Pu10" 			,&Pu10			);
+	reader->AddVariable( "Pu11" 			,&Pu11			);
+	reader->AddVariable( "Pu12" 			,&Pu12			);
+	reader->AddVariable( "Am1"  			,&Am1 			);
+
+	// --- Book the MVA methods
+
+	// Book method MLP
+	TString methodName = "MLP method";
+	reader->BookMVA( methodName, fTMVAWeightPath );
+	theTree->SetBranchAddress( "BU"   			,&BU 			);
+	theTree->SetBranchAddress( "U5_enrichment"  	,&U5_enrichment  	);
+	theTree->SetBranchAddress( "Pu8"  			,&Pu8  			);
+	theTree->SetBranchAddress( "Pu9"  			,&Pu9  			);
+	theTree->SetBranchAddress( "Pu10" 			,&Pu10 		);
+	theTree->SetBranchAddress( "Pu11" 			,&Pu11 		);
+	theTree->SetBranchAddress( "Pu12" 			,&Pu12 		);
+	theTree->SetBranchAddress( "Am1"  			,&Am1 			);
+	theTree->GetEntry(0);
+
+	Float_t val = (reader->EvaluateRegression( methodName ))[0];
+
+	delete reader;
+	delete theTree;
+	
+	MolarFraction["Fissile"] = (double)val;
+	MolarFraction["Fertile"] = 1.- (double)val;
+
+	return MolarFraction; //return Molar content of each component in the fuel
+}
+//________________________________________________________________________
+map < string , double> EQM_PWR_MLP_MOX::GetMolarFraction(map < string , IsotopicVector> IVStream, double BurnUp)
+{DBGL
+	return	ExecuteTMVA(CreateTMVAInputTree(IVStream,BurnUp));
+}

source/branches/CLASSV4.2/Model/Equivalence/EQM_PWR_MLP_MOX.hxx

+#ifndef _EQM_PWR_MLP_MOX_HXX
+#define _EQM_PWR_MLP_MOX_HXX
+
+#include "EquivalenceModel.hxx"
+#include "TTree.h"
+
+using namespace std;
+
+//-----------------------------------------------------------------------------//
+//! Defines an EquivalenceModel based on neural network
+
+/*!
+ The aim of these class is to constuct a fuel from an equivalence model
+ based on a  Multi layer perceptron
+
+ @author BLG
+ @version 3.0
+ */
+//________________________________________________________________________
+
+class EQM_PWR_MLP_MOX : public EquivalenceModel
+{
+	public :
+	/*!
+	 \name Constructor
+	 */
+	//@{
+	
+	//{
+	/// normal constructor
+	/*!
+	 Create a EQM_PWR_MLP_MOX 
+	 \param  TMVAWeightPath :  PAth to the .xml file containing neural network informations : PATH/TMVAWeight.xml (total path to tmva weight)
+	 */
+	EQM_PWR_MLP_MOX(string TMVAWeightPath);
+	//}
+	
+	//{
+	/// Logger constructor
+	/*!
+	 Create a EQM_PWR_MLP_MOX
+	 \param log : use for log
+	 \param  TMVAWeightPath :  PAth to the .xml file containing neural network informations : PATH/TMVAWeight.xml (total path to tmva weight)
+	 */
+	EQM_PWR_MLP_MOX(CLASSLogger* log, string TMVAWeightPath);
+	//}
+	//@}
+	
+	//{
+	/// Return the molar fissile fraction according fissile & ferile content using a Multi Layer Peceptron (MLP)
+	/*!
+	 \param Fissil : The composition of the fissile matter
+	 \param Fertil : The composition of the Fertil matter
+	 \param BurnUp : Maximum achievable burn up envisaged
+	 */
+	map < string , double> GetMolarFraction(map < string , IsotopicVector> IVStream, double BurnUp);
+	//}
+	
+	private :
+	
+	TTree* CreateTMVAInputTree(map < string , IsotopicVector> IVStream, double BurnUp);//!<Create input tmva tree to be read by ExecuteTMVA
+	map < string , double> ExecuteTMVA(TTree* theTree);//!<Execute the MLP according to the input tree created by CreateTMVAInputTree
+
+
+	string fTMVAWeightPath;;//!<The weight needed by TMVA to construct and execute the multilayer perceptron
+
+};
+
+#endif
+

source/branches/CLASSV4.2/Model/Equivalence/EQM_PWR_MLP_MOX_Am.cxx

+#include "EquivalenceModel.hxx"
+#include "EQM_PWR_MLP_MOX_Am.hxx"
+#include "CLASSLogger.hxx"
+#include "StringLine.hxx"
+
+#include <string>
+#include <iostream>
+#include <fstream>
+#include <algorithm>
+#include <cmath>
+#include <cassert>
+
+#include "TSystem.h"
+#include "TMVA/Reader.h"
+#include "TMVA/Tools.h"
+#include "TMVA/MethodCuts.h"
+
+
+//________________________________________________________________________
+//
+//		EQM_PWR_MLP_MOX_AM
+//
+//	Equivalenve Model based on multi layer perceptron from TMVA (root cern)
+//	For REP MOX use
+//
+//________________________________________________________________________
+
+EQM_PWR_MLP_MOX_AM::EQM_PWR_MLP_MOX_AM(string TMVAWeightPath):EquivalenceModel(new CLASSLogger("EQM_PWR_MLP_MOX_AM.log"))
+{
+	fTMVAWeightPath = TMVAWeightPath;
+
+	ZAI U8(92,238,0);
+	ZAI U5(92,235,0);
+	double U5_enrich = 0.0025;
+
+	ZAI Pu8(94,238,0);
+	ZAI Pu9(94,239,0);
+	ZAI Pu0(94,240,0);
+	ZAI Pu1(94,241,0);
+	ZAI Pu2(94,242,0);
+	ZAI Am1(95,241,0);
+	ZAI Am2(95,242,1);
+	ZAI Am3(95,243,0);
+
+	fFissileList = Pu8*1 + Pu9*1 + Pu0*1 + Pu1*1 + Pu2*1 + Am1*1 + Am2*1 + Am3*1;
+	fFertileList = U5*U5_enrich + U8*(1-U5_enrich);
+
+	SetBuildFuelFirstGuess(0.04);
+
+	fStreamList.push_back(FissileList);
+	fStreamList.push_back(FertileList);
+	
+	INFO << "__An equivalence model of PWR MOX has been define__" << endl;
+	INFO << "\tThis model is based on a multi layer perceptron" << endl;
+	INFO << "\t\tThe TMVA weight file is :" << endl;
+	INFO << "\t\t\t"<<fTMVAWeightPath << endl;
+
+}
+
+//________________________________________________________________________
+EQM_PWR_MLP_MOX_AM::EQM_PWR_MLP_MOX_AM(CLASSLogger* log, string TMVAWeightPath):EquivalenceModel(log)
+{
+	fTMVAWeightPath = TMVAWeightPath;
+
+	ZAI U8(92,238,0);
+	ZAI U5(92,235,0);
+	double U5_enrich = 0.0025;
+
+	ZAI Pu8(94,238,0);
+	ZAI Pu9(94,239,0);
+	ZAI Pu0(94,240,0);
+	ZAI Pu1(94,241,0);
+	ZAI Pu2(94,242,0);
+	ZAI Am1(95,241,0);
+	ZAI Am2(95,242,1);
+	ZAI Am3(95,243,0);
+	
+	fFissileList = Pu8*1 + Pu9*1 + Pu0*1 + Pu1*1 + Pu2*1 + Am1*1 + Am2*1 + Am3*1;
+	fFertileList = U5*U5_enrich + U8*(1-U5_enrich);
+
+	SetBuildFuelFirstGuess(0.04);
+
+	fStreamList.push_back(FissileList);
+	fStreamList.push_back(FertileList);
+
+	INFO << "__An equivalence model of PWR MOX has been define__" << endl;
+	INFO << "\tThis model is based on a multi layer perceptron" << endl;
+	INFO << "\t\tThe TMVA weight file is :" << endl;
+	INFO << "\t\t\t"<<fTMVAWeightPath << endl;
+
+}
+
+//________________________________________________________________________
+TTree* EQM_PWR_MLP_MOX_AM::CreateTMVAInputTree(vector <IsotopicVector> IVStream,double BurnUp)
+{
+	IsotopicVector Fissil = IVStream[0];
+	IsotopicVector Fertil = IVStream[1];
+
+	TTree*   InputTree = new TTree("EQTMP", "EQTMP");
+	float Pu8   		 = 0;
+	float Pu9   		 = 0;
+	float Pu10  		 = 0;
+	float Pu11  		 = 0;
+	float Pu12  		 = 0;
+	float Am1   		 = 0;
+	float Am2   		 = 0;
+	float Am3   		 = 0;
+	float U5_enrichment  = 0;
+	float BU  			  = 0;
+
+	InputTree->Branch(	"Pu8"	,&Pu8	,"Pu8/F"	);
+	InputTree->Branch(	"Pu9"	,&Pu9	,"Pu9/F"	);
+	InputTree->Branch(	"Pu10"	,&Pu10	,"Pu10/F"	);
+	InputTree->Branch(	"Pu11"	,&Pu11	,"Pu11/F"	);
+	InputTree->Branch(	"Pu12"	,&Pu12	,"Pu12/F"	);
+	InputTree->Branch(	"Am1"	,&Am1	,"Am1/F"	);
+	InputTree->Branch(	"Am2"	,&Am2	,"Am2/F"	);
+	InputTree->Branch(	"Am3"	,&Am3	,"Am3/F"	);
+	InputTree->Branch(	"U5_enrichment"	,&U5_enrichment	,"U5_enrichment/F"	);
+	InputTree->Branch(	"BU"	,&BU	,"BU/F"	);
+
+
+	float U8     = Fertil.GetZAIIsotopicQuantity(92,238,0);
+	float U5     = Fertil.GetZAIIsotopicQuantity(92,235,0);
+	float U4     = Fertil.GetZAIIsotopicQuantity(92,234,0);
+
+	float UTOT = U8 + U5 + U4;
+
+	Pu8    	   = Fissil.GetZAIIsotopicQuantity(94,238,0);
+	Pu9    	   = Fissil.GetZAIIsotopicQuantity(94,239,0);
+	Pu10   	   = Fissil.GetZAIIsotopicQuantity(94,240,0);
+	Pu11   	   = Fissil.GetZAIIsotopicQuantity(94,241,0);
+	Pu12   	   = Fissil.GetZAIIsotopicQuantity(94,242,0);
+	Am1        = Fissil.GetZAIIsotopicQuantity(95,241,0);
+	Am2        = Fissil.GetZAIIsotopicQuantity(95,242,1);
+	Am3        = Fissil.GetZAIIsotopicQuantity(95,243,0);
+
+	double TOTPU = (Pu8+Pu9+Pu10+Pu11+Pu12+Am1+Am2+Am3);
+
+	Pu8 = Pu8  / TOTPU;
+	Pu9 = Pu9  / TOTPU;
+	Pu10 = Pu10 / TOTPU;
+	Pu11 = Pu11 / TOTPU;
+	Pu12 = Pu12 / TOTPU;
+	Am1 = Am1  / TOTPU;
+	Am2 = Am2  / TOTPU;
+	Am3 = Am3  / TOTPU;
+	
+	U5_enrichment = U5 / UTOT;
+
+	BU = BurnUp;
+	if(Pu8 + Pu9 + Pu10 + Pu11 + Pu12 + Am1 + Am2 + Am3 > 1.00001 )//?????1.00001??? I don't know it! goes in condition if  = 1 !! may be float/double issue ...
+	{
+		ERROR << Pu8 << " " << Pu9 << " " << Pu10 << " " << Pu11 << " " << Pu12 << " " << Am1 << " " << Am2 << " " << Am3 << endl;
+		exit(0);
+	}
+	// All value are molar (!weight)
+
+	InputTree->Fill();
+	return InputTree;
+}
+//________________________________________________________________________
+double EQM_PWR_MLP_MOX_AM::ExecuteTMVA(TTree* theTree)
+{
+	// --- Create the Reader object
+	TMVA::Reader *reader = new TMVA::Reader( "Silent" );
+	// Create a set of variables and declare them to the reader
+	// - the variable names MUST corresponds in name and type to those given in the weight file(s) used
+	Float_t Pu8,Pu9,Pu10,Pu11,Pu12,Am1, Am2,Am3,BU,U5_enrichment;
+
+	reader->AddVariable( "BU"   		,&BU );
+	reader->AddVariable( "U5_enrichment",&U5_enrichment );
+	reader->AddVariable( "Pu8"  		,&Pu8 );
+	reader->AddVariable( "Pu9"  		,&Pu9 );
+	reader->AddVariable( "Pu10" 		,&Pu10);
+	reader->AddVariable( "Pu11" 		,&Pu11);
+	reader->AddVariable( "Pu12" 		,&Pu12);
+	reader->AddVariable( "Am1"  		,&Am1 );
+	reader->AddVariable( "Am2"  		,&Am2 );
+	reader->AddVariable( "Am3"  		,&Am3 );
+
+	// --- Book the MVA methods
+
+	// Book method MLP
+	TString methodName = "MLP method";
+	reader->BookMVA( methodName, fTMVAWeightPath );
+	theTree->SetBranchAddress( "BU"   			,&BU 	);
+	theTree->SetBranchAddress( "U5_enrichment"		,&U5_enrichment  )	;
+	theTree->SetBranchAddress( "Pu8"  			,&Pu8  );
+	theTree->SetBranchAddress( "Pu9"  			,&Pu9  );
+	theTree->SetBranchAddress( "Pu10" 			,&Pu10 );
+	theTree->SetBranchAddress( "Pu11" 			,&Pu11 );
+	theTree->SetBranchAddress( "Pu12" 			,&Pu12 );
+	theTree->SetBranchAddress( "Am1"  			,&Am1  );
+	theTree->SetBranchAddress( "Am2"  			,&Am2  );
+	theTree->SetBranchAddress( "Am3"  			,&Am3  );
+	theTree->GetEntry(0);
+
+	Float_t val = (reader->EvaluateRegression( methodName ))[0];
+
+	delete reader;
+	delete theTree;
+
+	return (double)val; //retourne teneur
+}
+//________________________________________________________________________
+double EQM_PWR_MLP_MOX_AM::GetFissileMolarFraction(vector <IsotopicVector> IVStream,double BurnUp)
+{DBGL
+	return	ExecuteTMVA(CreateTMVAInputTree(IVStream,BurnUp));
+}

source/branches/CLASSV4.2/Model/Equivalence/EQM_PWR_MLP_MOX_Am.hxx

+#ifndef _EQM_PWR_MLP_MOX_AM_HXX
+#define _EQM_PWR_MLP_MOX_AM_HXX
+
+#include "EquivalenceModel.hxx"
+#include "TTree.h"
+
+using namespace std;
+
+//-----------------------------------------------------------------------------//
+//! Defines an EquivalenceModel based on neural network
+
+/*!
+ The aim of these class is to constuct a fuel from an equivalence model
+ based on a  Multi layer perceptron
+
+ @author BLG
+ @version 3.0
+ */
+//________________________________________________________________________
+
+class EQM_PWR_MLP_MOX_AM : public EquivalenceModel
+{
+	public :
+	/*!
+	 \name Constructor
+	 */
+	//@{
+	
+	//{
+	/// normal constructor
+	/*!
+	 Create a EQM_PWR_MLP_MOX_AM 
+	 \param  TMVAWeightPath :  PAth to the .xml file containing neural network informations : PATH/TMVAWeight.xml (total path to tmva weight)
+	 */
+	EQM_PWR_MLP_MOX_AM(string TMVAWeightPath);
+	//}
+	
+	//{
+	/// Logger constructor
+	/*!
+	 Create a EQM_PWR_MLP_MOX_AM
+	 \param log : use for log
+	 \param  TMVAWeightPath :  PAth to the .xml file containing neural network informations : PATH/TMVAWeight.xml (total path to tmva weight)
+	 */
+	EQM_PWR_MLP_MOX_AM(CLASSLogger* log, string TMVAWeightPath);
+	//}
+	//@}
+	
+	//{
+	/// Return the molar fissile fraction according fissile & ferile content using a Multi Layer Peceptron (MLP)
+	/*!
+	 \param Fissil : The composition of the fissile matter
+	 \param Fertil : The composition of the Fertil matter
+	 \param BurnUp : Maximum achievable burn up envisaged
+	 */
+	virtual double GetFissileMolarFraction(vector <IsotopicVector> IVStream,double BurnUp);
+	//}
+	
+	private :
+	
+	TTree* CreateTMVAInputTree(vector <IsotopicVector> IVStream,double BurnUp);//!<Create input tmva tree to be read by ExecuteTMVA
+	double ExecuteTMVA(TTree* theTree);//!<Execute the MLP according to the input tree created by CreateTMVAInputTree
+
+
+	string fTMVAWeightPath;;//!<The weight needed by TMVA to construct and execute the multilayer perceptron
+
+};
+
+#endif
+

source/branches/CLASSV4.2/Model/Equivalence/EQM_PWR_POL_UO2.cxx

+#include "EquivalenceModel.hxx"
+#include "EQM_PWR_POL_UO2.hxx"
+#include "CLASSLogger.hxx"
+#include "StringLine.hxx"
+
+
+// ________________________________________________________________________
+// EQM_PWR_POL_UO2
+//
+// ________________________________________________________________________
+
+
+
+
+//Constructor(s)
+EQM_PWR_POL_UO2::EQM_PWR_POL_UO2(string PathToWeightFile):EquivalenceModel(new CLASSLogger("EQM_PWR_POL_UO2.log"))
+{
+	
+	// Fertile
+	ZAI U8(92 ,238 ,0) ;
+	fFertileList = U8*1;
+	// Fissile
+	ZAI U5(92 ,235 ,0) ;
+	// ...
+	fFissileList = U5*1;
+
+	ReadWeightFile(PathToWeightFile);
+
+}
+// _______________________________________________________________________
+EQM_PWR_POL_UO2::EQM_PWR_POL_UO2(CLASSLogger* log,string PathToWeightFile):EquivalenceModel(log)
+{
+	
+	// Fertile
+	ZAI U8(92 ,238 ,0) ;
+	fFertileList = U8*1;
+	// Fissile
+	ZAI U5(92 ,235 ,0) ;
+	// ...
+	fFissileList = U5*1;
+
+	ReadWeightFile(PathToWeightFile);
+
+}
+// _______________________________________________________________________
+void EQM_PWR_POL_UO2::ReadWeightFile(string PathToWeightFile)
+{
+	ifstream DataDB(PathToWeightFile.c_str());							// Open the File
+	if(!DataDB)
+		WARNING << " Can't open \"" << PathToWeightFile << "\"" << endl;
+	
+	string line;
+	int start = 0;	// First Get Fuel Parameter
+	getline(DataDB, line);
+	
+	if( StringLine::NextWord(line, start, ' ') != "PARAM")
+	{
+		ERROR << "Bad Database file : " <<  PathToWeightFile << " Can't find the Parameter of the DataBase " << endl;
+		exit (1);
+	}
+	fParam_Bu_0 = atof(StringLine::NextWord(line, start, ' ').c_str()) ;
+	fParam_Bu = atof(StringLine::NextWord(line, start, ' ').c_str());
+	fParam_BuSquare = atof(StringLine::NextWord(line, start, ' ').c_str());
+	
+	INFO << "Weight parameters has been read " << endl;
+	INFO  << "\t U enrichment = " << fParam_Bu_0 << " + " << fParam_Bu << "*Burnup + " <<  fParam_BuSquare << "*Burnup*Burnup" << endl;
+}
+// _______________________________________________________________________
+double EQM_PWR_POL_UO2::GetFissileMolarFraction ( vector <IsotopicVector> IVStream , double BurnUp )
+{
+	return fParam_Bu_0 + fParam_Bu*BurnUp + fParam_BuSquare*BurnUp*BurnUp ;
+	
+}
\ No newline at end of file

source/branches/CLASSV4.2/Model/Equivalence/EQM_PWR_POL_UO2.hxx

+#ifndef _EQM_PWR_POL_UO2_HXX
+#define _EQM_PWR_POL_UO2_HXX
+
+#include "EquivalenceModel.hxx"
+
+using namespace std;
+
+//−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−--------------−−−−−−−−−−−−−−−//
+//! Define an EquivalenceModel based on a polynomial fit
+
+/*!
+ Defines a EQM_PWR_POL_UO2
+ It returns the @f$^{235}U@f$ enrichment e according to this polynom :
+ 
+ @f$e = \alpha_{0} + \alpha_{1}\cdot Burnup + \alpha_{2}\cdot Burnup \cdot Burnup @f$
+ 
+ BU : Maximum achievable burnup
+ 
+ @author BaM
+ @version 3.0
+ */
+//________________________________________________________________________
+
+class EQM_PWR_POL_UO2 : public EquivalenceModel
+{
+	
+public:
+	/*!
+	 \name Constructor
+	 */
+	//@{
+	
+	//{
+	/// normal constructor
+	/*!
+	 Create a EQM_PWR_POL_UO2
+	 \param  PathToWeightFile :  Path to the file containing the  @f$\alpha_{i}@f$
+	 Format :  @f$PARAM@f$  @f$\alpha_{0}@f$  @f$\alpha_{1}@f$  @f$\alpha_{2}@f$
+	 */
+	EQM_PWR_POL_UO2(string PathToWeightFile);
+	//}
+	
+	//{
+	/// logger constructor
+	/*!
+	 Create a EQM_PWR_POL_UO2
+	 \param log : Use for the log
+	 \param  PathToWeightFile :  Path to the file containing the  @f$\alpha_{i}@f$
+	 Format :  @f$PARAM@f$  @f$\alpha_{0}@f$  @f$\alpha_{1}@f$  @f$\alpha_{2}@f$
+	 */
+	EQM_PWR_POL_UO2(CLASSLogger* log, string PathToWeightFile);
+	//}
+	
+	//@}
+	/**This function IS the equivalence model**/
+	double GetFissileMolarFraction(vector <IsotopicVector> IVStream,double BurnUp) ; // !<Return the molar fraction of fissile element
+	
+	
+private:
+	
+	void ReadWeightFile(string PathToWeightFile); //!< Function to read the weight file & file the parameters
+	
+	double fParam_Bu_0 ;		//!<  @f$\alpha_{0}@f$
+	double fParam_Bu ;		//!<  @f$\alpha_{1}@f$
+	double fParam_BuSquare ;	//!<  @f$\alpha_{2}@f$
+
+};
+
+#endif
\ No newline at end of file

source/branches/CLASSV4.2/Model/Equivalence/EQM_PWR_QUAD_MOX.cxx

+#include "EQM_PWR_QUAD_MOX.hxx"
+
+#include <vector>
+
+#include "StringLine.hxx"
+#include "CLASSLogger.hxx"
+
+
+
+
+EQM_PWR_QUAD_MOX::EQM_PWR_QUAD_MOX(string WeightPath):EquivalenceModel(new CLASSLogger("EQM_PWR_QUAD_MOX.log"))
+{
+	fWeightPath = WeightPath;
+
+	ifstream DataDB(fWeightPath.c_str());							// Open the File
+	if(!DataDB)
+		WARNING << " Can't open \"" << fWeightPath << "\"" << endl;
+
+	string line;
+	int start = 0;	// First Get Fuel Parameter
+	getline(DataDB, line);
+
+	if( StringLine::NextWord(line, start, ' ') != "PARAM")
+	{
+		ERROR << "Bad Database file : " <<  fWeightPath << " Can't find the Parameter of the DataBase " << endl;
+		exit (1);
+	}
+	while(start < (int)line.size())
+		fFuelParameter.push_back(atof(StringLine::NextWord(line, start, ' ').c_str()));
+
+	INFO << fFuelParameter.size() << " have been read " << endl;
+
+
+	ZAI U8(92,238,0);
+	ZAI U5(92,235,0);
+	double U5_enrich = 0.0025;
+	fFertileList = U5*U5_enrich + U8*(1-U5_enrich);
+
+
+	ZAI Pu8(94,238,0);
+	ZAI Pu9(94,239,0);
+	ZAI Pu0(94,240,0);
+	ZAI Pu1(94,241,0);
+	ZAI Pu2(94,242,0);
+	fFissileList = Pu8*1+Pu9*1+Pu0*1+Pu1*1+Pu2*1;
+
+	SetBuildFuelFirstGuess(0.04);
+	
+	fStreamList.push_back(FissileList);
+	fStreamList.push_back(FertileList);
+
+}
+
+EQM_PWR_QUAD_MOX::EQM_PWR_QUAD_MOX(CLASSLogger* log, string WeightPath):EquivalenceModel(log)
+{
+	fWeightPath = WeightPath;
+
+	ifstream DataDB(fWeightPath.c_str());							// Open the File
+	if(!DataDB)
+		WARNING << " Can't open \"" << fWeightPath << "\"" << endl;
+
+	string line;
+	int start = 0;	// First Get Fuel Parameter
+	getline(DataDB, line);
+
+	if( StringLine::NextWord(line, start, ' ') != "PARAM")
+	{
+		ERROR << "Bad Database file : " <<  fWeightPath << " Can't find the Parameter of the DataBase " << endl;
+		exit (1);
+	}
+	while(start < (int)line.size())
+		fFuelParameter.push_back(atof(StringLine::NextWord(line, start, ' ').c_str()));
+
+	INFO << fFuelParameter.size() << " have been read " << endl;
+
+
+	ZAI U8(92,238,0);
+	ZAI U5(92,235,0);
+	double U5_enrich = 0.0025;
+	fFertileList = U5*U5_enrich + U8*(1-U5_enrich);
+
+
+	ZAI Pu8(94,238,0);
+	ZAI Pu9(94,239,0);
+	ZAI Pu0(94,240,0);
+	ZAI Pu1(94,241,0);
+	ZAI Pu2(94,242,0);
+	fFissileList = Pu8*1+Pu9*1+Pu0*1+Pu1*1+Pu2*1;
+
+	SetBuildFuelFirstGuess(0.04);
+
+	fStreamList.push_back(FissileList);
+	fStreamList.push_back(FertileList);
+
+}
+
+
+
+EQM_PWR_QUAD_MOX::~EQM_PWR_QUAD_MOX()
+{
+
+}
+
+
+
+
+double EQM_PWR_QUAD_MOX::GetFissileMolarFraction(vector <IsotopicVector> IVStream,double BurnUp)
+{
+
+	IsotopicVector Fissile = IVStream[0];
+	IsotopicVector Fertile = IVStream[1];
+
+	ZAI ZAIList[6] = {ZAI(94,238,0), ZAI(94,239,0), ZAI(94,240,0), ZAI(94,241,0), ZAI(94,242,0), ZAI(95,241,0)  };
+
+	vector<double> PuCompo;
+	double Sum = Fissile.GetSumOfAll();
+
+
+	for(int i = 0; i< 5; i++)
+		PuCompo.push_back( Fissile.GetZAIIsotopicQuantity(ZAIList[i])/Sum);
+
+	PuCompo[2] += Fissile.GetZAIIsotopicQuantity(ZAIList[5])/Sum;
+	double A = 0;
+
+	if(PuCompo[0] <= PuCompo[2] && PuCompo[0] <= PuCompo[4] && PuCompo[1] + PuCompo[3] >= 0.40 && PuCompo[1] >0 )
+	{
+		int par = 0;
+		for(int j = 0 ; j < 5 ; j++)
+		{
+			A += fFuelParameter[par]  * PuCompo[j] ;
+			par++;
+			for(int i = j ; i < 5 ; i++)
+			{
+				A += fFuelParameter[par] *PuCompo[i] *PuCompo[j];
+				par++;
+			}
+		}
+		A += fFuelParameter[par];
+	}
+	else
+	{
+		cout << "the composition is not in the range of the Model" << endl;
+	}
+	return A;
+}
+

source/branches/CLASSV4.2/Model/Equivalence/EQM_PWR_QUAD_MOX.hxx

+#ifndef _EQM_PWR_QUAD_MOX_HXX
+#define _EQM_PWR_QUAD_MOX_HXX
+
+#include "EquivalenceModel.hxx"
+
+#include <string>
+
+using namespace std;
+
+//-----------------------------------------------------------------------------//
+//! Defines an EquivalenceModel based on a quadratic fit
+
+/*!
+ The aim of these class is to constuct a fuel from an equivalence model
+ based on a Quadratic Pu equivalent Model
+ The Plutonium content e is calculated using :
+ 
+ @f$ e = \alpha_{0} + \sum_{i\in Pu}^{N} \left(\alpha_{i} \cdot n_{i}\ + \sum_{j\leq i} \alpha_{ij} \cdot n_{i}\cdot n_{j}\right)@f$
+ 
+ @author BaM
+ @version 3.0
+ */
+//________________________________________________________________________
+
+class EQM_PWR_QUAD_MOX : public EquivalenceModel
+{
+	public :
+
+	/*!
+	 \name Constructor
+	 */
+	//@{
+	
+	//{
+	/// normal constructor
+	/*!
+	 Create a EQM_PWR_POL_UO2
+	 \param  WeightPath :  Path to the file containing the  @f$\alpha_{i}@f$
+	 Format :  @f$PARAM@f$   @f$\alpha_{^{238}Pu}@f$  @f$\alpha_{^{238}Pu ^{238}Pu}@f$  @f$\alpha_{^{238}Pu ^{239}Pu}@f$ .... @f$\alpha_{^{238}Pu ^{242}Pu}@f$ @f$\alpha_{^{239}Pu}@f$ ...
+	 @f$\alpha_{^{242}Pu}@f$  @f$\alpha_{^{242}Pu ^{242}Pu}@f$ @f$\alpha_{0}@f$
+	 */
+	EQM_PWR_QUAD_MOX(string WeightPath);
+	//}
+	
+	//{
+	/// Logger constructor
+	/*!
+	 Create a EQM_PWR_POL_UO2
+	 \param log : Use for the log
+	 \param  WeightPath :  Path to the file containing the  @f$\alpha_{i}@f$
+	 Format :  @f$PARAM@f$   @f$\alpha_{^{238}Pu}@f$  @f$\alpha_{^{238}Pu ^{238}Pu}@f$  @f$\alpha_{^{238}Pu ^{239}Pu}@f$ .... @f$\alpha_{^{238}Pu ^{242}Pu}@f$ @f$\alpha_{^{239}Pu}@f$ ...
+	 @f$\alpha_{^{242}Pu}@f$  @f$\alpha_{^{242}Pu ^{242}Pu}@f$ @f$\alpha_{0}@f$
+	 */
+	EQM_PWR_QUAD_MOX(CLASSLogger* log, string WeightPath);
+	//}
+	
+	~EQM_PWR_QUAD_MOX();
+	//@}
+	
+	virtual double GetFissileMolarFraction(vector <IsotopicVector> IVStream,double BurnUp);
+
+	private :
+
+	string fWeightPath;		//!< Path to the file containing the @f$\alpha_{ij}@f$
+	vector<double> fFuelParameter;	//!< vector containing the @f$\alpha_{ij}@f$
+
+};
+
+#endif
+

source/branches/CLASSV4.2/Model/Irradiation/IM_Matrix.cxx

+//
+//  IM_Matrix.cxx
+//  CLASSSource
+//
+//  Created by BaM on 04/05/2014.
+//  Copyright (c) 2014 BaM. All rights reserved.
+//
+
+#include "IM_Matrix.hxx"
+
+#include "IsotopicVector.hxx"
+#include "CLASSConstante.hxx"
+#include "CLASSLogger.hxx"
+#include "StringLine.hxx"
+
+#include <TGraph.h>
+#include <TString.h>
+
+
+#include <sstream>
+#include <string>
+#include <iostream>
+#include <fstream>
+#include <algorithm>
+#include <cmath>
+
+
+
+using namespace std;
+
+
+//________________________________________________________________________
+IM_Matrix::IM_Matrix():IrradiationModel(new CLASSLogger("IM_Matrix.log"))
+{
+	fShorstestHalflife = 3600.*24*160.; //cut by default all nuclei with a shorter liftime than the Cm242 -> remain 33 actinides
+}
+
+
+IM_Matrix::IM_Matrix(CLASSLogger* log):IrradiationModel(log)
+{
+	fShorstestHalflife = 3600.*24*160.; //cut by default all nuclei with a shorter liftime than the Cm242 -> remain 33 actinides
+}
+
+
+
+
+//________________________________________________________________________
+/*			Evolution Calculation			*/
+//________________________________________________________________________
+EvolutionData IM_Matrix::GenerateEvolutionData(IsotopicVector isotopicvector, EvolutionData XSSet, double Power, double cycletime)
+{
+	DBGL
+	if(fFastDecay.size() ==  0)
+		NuclearDataInitialization();
+
+
+	string ReactorType;
+
+
+	vector< TMatrixT<double> > NMatrix ;//  TMatrixT<double>(decayindex.size(),1))
+	{	// Filling the t = 0 State;
+		map<ZAI, double > isotopicquantity = isotopicvector.GetIsotopicQuantity();
+		TMatrixT<double>  N_0Matrix = TMatrixT<double>( fReverseMatrixIndex.size(),1) ;
+		for(int i = 0; i < (int)fReverseMatrixIndex.size(); i++)
+			N_0Matrix[i] = 0;
+
+		map<ZAI, double >::iterator it ;
+		for(int i = 0; i < (int)fReverseMatrixIndex.size(); i++)
+			N_0Matrix[i] = 0;
+
+		for(it = isotopicquantity.begin(); it != isotopicquantity.end(); it++)
+		{
+			/// Need TO change with FP managment
+			map<ZAI, int >::iterator it2;
+
+			if( (*it).first.Z() < fZAIThreshold )
+				it2 = fMatrixIndex.find( ZAI(-2,-2,-2) );
+			else it2 = fMatrixIndex.find( (*it).first );
+
+			if(it2 ==  fMatrixIndex.end() )		//If not in index should be TMP, can't be fast decay for new Fuel !!!
+				it2 = fMatrixIndex.find( ZAI(-3,-3,-3) );
+			N_0Matrix[ (*it2).second ][0] = (*it).second ;
+
+
+		}
+
+		isotopicquantity.clear();
+
+		NMatrix.push_back(N_0Matrix);
+		N_0Matrix.Clear();
+
+	}
+
+
+	//-------------------------//
+	//--- Perform Evolution ---//
+	//-------------------------//
+	ReactorType = XSSet.GetReactorType();
+
+	double M_ref = XSSet.GetHeavyMetalMass();
+	double M = cZAIMass.GetMass(isotopicvector.GetActinidesComposition());
+	double Power_ref = XSSet.GetPower();
+
+	int NStep = XSSet.GetFissionXS().begin()->second->GetN();
+	double* DBTimeStep = XSSet.GetFissionXS().begin()->second->GetX();
+
+	int InsideStep = 10;
+
+	double timevector[NStep];
+	timevector[0] = 0;
+
+	double  Flux[NStep];
+
+	TMatrixT<double> FissionEnergy = TMatrixT<double>(fReverseMatrixIndex.size(),1);
+	for(int i = 0; i < (int)fReverseMatrixIndex.size(); i++)
+		FissionEnergy[i] = 0;
+
+	{
+		map< ZAI, int >::iterator it;
+		for(it = fMatrixIndex.begin(); it != fMatrixIndex.end(); it++)
+		{
+			map< ZAI, double >::iterator it2 = fFissionEnergy.find(it->first);
+			if(it2 ==  fFissionEnergy.end())
+			{
+				if(it->first.Z() > fZAIThreshold)
+					FissionEnergy[it->second][0] = 1.9679e6*it->first.A()-2.601e8; // //simple linear fit to known values ;extrapolation to unknown isotopes
+				else FissionEnergy[it->second][0] = 0;
+			}
+			else
+				FissionEnergy[it->second][0] = it2->second;
+
+		}
+	}
+
+	vector< TMatrixT<double> > FissionXSMatrix;	// Store The Fisison XS Matrix
+	vector< TMatrixT<double> > CaptureXSMatrix;	// Store The Capture XS Matrix
+	vector< TMatrixT<double> > n2nXSMatrix;		// Store The n2N XS Matrix
+	DBGL
+	for(int i = 0; i < NStep-1; i++)
+	{
+		double TStepMax = ( (DBTimeStep[i+1]-DBTimeStep[i] ) ) * Power_ref/M_ref / Power*M ;	// Get the next Time step
+
+
+		TMatrixT<double> BatemanMatrix = TMatrixT<double>(fReverseMatrixIndex.size(),fReverseMatrixIndex.size());
+		TMatrixT<double> BatemanReactionMatrix = TMatrixT<double>(fReverseMatrixIndex.size(),fReverseMatrixIndex.size());
+
+		TMatrixT<double> NEvolutionMatrix = TMatrixT<double>(fReverseMatrixIndex.size(),1);
+		NEvolutionMatrix = NMatrix.back();
+
+
+
+		FissionXSMatrix.push_back(GetFissionXsMatrix(XSSet, DBTimeStep[i]));	//Feel the fission reaction Matrix
+		CaptureXSMatrix.push_back(GetCaptureXsMatrix(XSSet, DBTimeStep[i]));	//Feel the capture reaction Matrix
+		n2nXSMatrix.push_back(Getn2nXsMatrix(XSSet, DBTimeStep[i]));		//Feel the (n,2n)  reaction Matrix
+
+		// ----------------   Evolution
+
+		BatemanReactionMatrix = FissionXSMatrix[i];
+		BatemanReactionMatrix += CaptureXSMatrix[i];
+		BatemanReactionMatrix += n2nXSMatrix[i];
+
+		for(int k = 0; k < InsideStep; k++)
+		{
+			double ESigmaN = 0;
+			for (int j = 0; j < (int)fReverseMatrixIndex.size() ; j++)
+				ESigmaN -=  FissionXSMatrix[i][j][j]*NEvolutionMatrix[j][0]*1.6e-19*FissionEnergy[j][0];
+			// Update Flux
+			double Flux_k = Power/ESigmaN;
+
+			if(k == 0)
+				Flux[i] = Flux_k;
+
+			BatemanMatrix = BatemanReactionMatrix;
+			BatemanMatrix *=  Flux_k;
+			BatemanMatrix += fDecayMatrix ;
+			BatemanMatrix *=  TStepMax/InsideStep ;
+
+
+			TMatrixT<double> IdMatrix = TMatrixT<double>(fReverseMatrixIndex.size(),fReverseMatrixIndex.size());
+			for(int j = 0; j < (int)fReverseMatrixIndex.size(); j++)
+				for(int k = 0; k < (int)fReverseMatrixIndex.size(); k++)
+				{
+					if(k ==  j)	IdMatrix[j][k] = 1;
+					else 		IdMatrix[j][k] = 0;
+				}
+
+
+			TMatrixT<double> BatemanMatrixDL = TMatrixT<double>(fReverseMatrixIndex.size(),fReverseMatrixIndex.size());   // Order 0 Term from the DL : Id
+			TMatrixT<double> BatemanMatrixDLTermN = TMatrixT<double>(fReverseMatrixIndex.size(),fReverseMatrixIndex.size());  // Addind it;
+
+			{
+				BatemanMatrix *=  TStepMax ;
+				BatemanMatrixDLTermN = IdMatrix;
+				BatemanMatrixDL = BatemanMatrixDLTermN;
+				int j = 1;
+				double NormN;
+
+				do
+				{
+					TMatrixT<double> BatemanMatrixDLTermtmp = TMatrixT<double>(fReverseMatrixIndex.size(),fReverseMatrixIndex.size());  // Adding it;
+					BatemanMatrixDLTermtmp = BatemanMatrixDLTermN;
+
+					BatemanMatrixDLTermN.Mult(BatemanMatrixDLTermtmp, BatemanMatrix );
+
+					BatemanMatrixDLTermN *=  1./j;
+					BatemanMatrixDL += BatemanMatrixDLTermN;
+
+					NormN = 0;
+					for(int m = 0; m < (int)fReverseMatrixIndex.size(); m++)
+						for(int n = 0; n < (int)fReverseMatrixIndex.size(); n++)
+							NormN += BatemanMatrixDLTermN[m][n]*BatemanMatrixDLTermN[m][n];
+					j++;
+
+				} while ( NormN != 0 );
+			}
+			NEvolutionMatrix = BatemanMatrixDL * NEvolutionMatrix ;
+		}
+		NMatrix.push_back(NEvolutionMatrix);
+
+
+		timevector[i+1] = timevector[i] + TStepMax;
+
+		BatemanMatrix.Clear();
+		BatemanReactionMatrix.Clear();
+		NEvolutionMatrix.Clear();
+
+
+	}
+	DBGL
+	FissionXSMatrix.push_back(GetFissionXsMatrix(XSSet, DBTimeStep[NStep-1])); //Feel the reaction Matrix
+	CaptureXSMatrix.push_back(GetCaptureXsMatrix(XSSet, DBTimeStep[NStep-1])); //Feel the reaction Matrix
+	n2nXSMatrix.push_back(Getn2nXsMatrix(XSSet, DBTimeStep[NStep-1])); //Feel the reaction Matrix
+
+
+	EvolutionData GeneratedDB = EvolutionData(GetLog());
+
+	double ESigmaN = 0;
+	for (int j = 0; j < (int)fReverseMatrixIndex.size() ; j++)
+		ESigmaN -=  FissionXSMatrix.back()[j][j]*NMatrix.back()[j][0]*1.6e-19*FissionEnergy[j][0];
+
+	Flux[NStep-1] = Power/ESigmaN;
+
+	GeneratedDB.SetFlux( new TGraph(NStep, timevector, Flux) );
+
+	for(int i = 0; i < (int)fReverseMatrixIndex.size(); i++)
+	{
+		double ZAIQuantity[NMatrix.size()];
+		double FissionXS[NStep];
+		double CaptureXS[NStep];
+		double n2nXS[NStep];
+		for(int j = 0; j < (int)NMatrix.size(); j++)
+			ZAIQuantity[j] = (NMatrix[j])[i][0];
+
+		for(int j = 0; j < NStep; j++)
+		{
+			FissionXS[j] = FissionXSMatrix[j][i][i];
+			CaptureXS[j] = CaptureXSMatrix[j][i][i];
+			n2nXS[j] = n2nXSMatrix[j][i][i];
+		}
+
+		GeneratedDB.NucleiInsert(pair<ZAI, TGraph*> (fReverseMatrixIndex[i], new TGraph(NMatrix.size(), timevector, ZAIQuantity)));
+/*		GeneratedDB.FissionXSInsert(pair<ZAI, TGraph*> (fReverseMatrixIndex[i], new TGraph(NStep, timevector, FissionXS)));
+		GeneratedDB.CaptureXSInsert(pair<ZAI, TGraph*> (fReverseMatrixIndex[i], new TGraph(NStep, timevector, CaptureXS)));
+		GeneratedDB.n2nXSInsert(pair<ZAI, TGraph*> (fReverseMatrixIndex[i], new TGraph(NStep, timevector, n2nXS)));
+*/	}
+
+	GeneratedDB.SetPower(Power );
+	GeneratedDB.SetHeavyMetalMass(M);
+	GeneratedDB.SetReactorType(ReactorType );
+	GeneratedDB.SetCycleTime(cycletime);
+
+	for (int i = 0; i < (int) FissionXSMatrix.size(); i++)
+	{
+		FissionXSMatrix[i].Clear();
+		CaptureXSMatrix[i].Clear();
+		n2nXSMatrix[i].Clear();
+	}
+	FissionXSMatrix.clear();
+	CaptureXSMatrix.clear();
+	n2nXSMatrix.clear();
+	DBGL
+	return GeneratedDB;
+	
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+

source/branches/CLASSV4.2/Model/Irradiation/IM_Matrix.hxx

+#ifndef _IMMATRIX_
+#define _IMMATRIX_
+
+
+/*!
+ \file
+ \brief Header file for IrradiationModel class.
+
+
+ @author BaM
+ @version 2.0
+ */
+#include "IrradiationModel.hxx"
+
+
+using namespace std;
+
+
+class CLASSLogger;
+
+//-----------------------------------------------------------------------------//
+//! Defines an IrradiationModel based on power series of the exponential of the Bateman matrix
+
+/*!
+ Define a IM_Matrix.
+ The aim of these class is to solve numericaly the Bateman equations using the
+ development in a power series of the exponential of the Bateman matrix.
+
+ @author BaM
+ @version 3.0
+ */
+//________________________________________________________________________
+
+class EvolutionData;
+
+class IM_Matrix : public IrradiationModel
+{
+
+	public :
+	/*!
+	 \name Constructor
+	 */
+	//@{
+	
+	//{
+	/// Default constructor
+	
+	/*!
+	 Make a new IM_Matrix : */
+	IM_Matrix();
+	//}
+	
+	/// Logger constructor
+	
+	/*!
+	 Make a new IM_Matrix : */
+	//param log : Use for the log
+	IM_Matrix(CLASSLogger* log);
+	//}
+	
+	//@}
+
+
+	/// virtual method called to perform the irradiation calculation using a set of cross section.
+	/*!
+	 Perform the Irradiation Calcultion using the XSSet data
+	 \param IsotopicVector IV isotopic vector to irradiate
+	 \param EvolutionData XSSet set of corss section to use to perform the evolution calculation
+	 */
+	virtual EvolutionData GenerateEvolutionData(IsotopicVector IV, EvolutionData XSSet, double Power, double cycletime);
+	//}
+
+
+
+
+	private :
+
+	
+	
+};
+
+#endif
+