diff --git a/source/branches/CLASSV3/Model/Equivalence/EQM_LIN_PWR_MOX.cxx b/source/branches/CLASSV3/Model/Equivalence/EQM_LIN_PWR_MOX.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..5fc357c4ff3b0fdf5d4776c6b8027710e44dac27
--- /dev/null
+++ b/source/branches/CLASSV3/Model/Equivalence/EQM_LIN_PWR_MOX.cxx
@@ -0,0 +1,200 @@
+#include "EQM_LIN_PWR_MOX.hxx"
+
+#include <vector>
+
+#include "StringLine.hxx"
+#include "LogFile.hxx"
+#include "IsotopicVector.hxx"
+#include "CLASSHeaders.hxx"
+
+
+
+EQM_LIN_PWR_MOX::EQM_LIN_PWR_MOX(string WeightPath):EquivalenceModel()
+{
+ fWeightPath = WeightPath;
+
+ ifstream DataDB(fWeightPath.c_str()); // Open the File
+ if(!DataDB)
+ cout << "!!Warning!! !!!EQM QUAD PWR MOX!!! \n Can't open \"" << fWeightPath << "\"\n" << endl;
+
+ string line;
+ int start = 0; // First Get Fuel Parameter
+ getline(DataDB, line);
+
+ if( StringLine::NextWord(line, start, ' ') != "PARAM")
+ {
+ cout << "!!Bad Trouble!! !!!EQM QUAD PWR MOX!!! 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()));
+
+ cout << "!!INFO!! !!!EQM QUAD PWR MOX!!! " << fFuelParameter.size() << " have been read"<< endl;
+
+
+
+ //-----------------------------------------------------------------------------//
+ //-----------------------------------------------------------------------------//
+ //-----------------------------------------------------------------------------//
+ //-----------------------------------------------------------------------------//
+ //-----------------------------------------------------------------------------//
+ // ADD ENrichment of the U reading !!!!!!!!!!!!!!!!!!!!!!!!!!!! //
+ //-----------------------------------------------------------------------------//
+ //-----------------------------------------------------------------------------//
+ //-----------------------------------------------------------------------------//
+ //-----------------------------------------------------------------------------//
+
+ ZAI U8(92,238,0);
+ ZAI U5(92,235,0);
+ double U5_enrich= 0.0025;
+ fFertileList = U5*U5_enrich + U8*(1-U5_enrich);
+
+
+ ZAI Pu8(94,238,0);
+ ZAI Pu9(94,239,0);
+ ZAI Pu0(94,240,0);
+ ZAI Pu1(94,241,0);
+ ZAI Pu2(94,242,0);
+ fFissileList = Pu8*1+Pu9*1+Pu0*1+Pu1*1+Pu2*1;
+
+
+}
+
+EQM_LIN_PWR_MOX::~EQM_LIN_PWR_MOX()
+{
+
+}
+
+//________________________________________________________________________
+vector<double> EQM_LIN_PWR_MOX::BuildFuel(double BurnUp, double HMMass,vector<IsotopicVector> FissilArray, vector<IsotopicVector> FertilArray)
+{
+
+ //-----------------------------------------------------------------------------//
+ //-----------------------------------------------------------------------------//
+ //-----------------------------------------------------------------------------//
+ //-----------------------------------------------------------------------------//
+ //-----------------------------------------------------------------------------//
+ // ADD ENrichment of the U check ++ Check Un seul fertile !!!! //
+ //-----------------------------------------------------------------------------//
+ //-----------------------------------------------------------------------------//
+ //-----------------------------------------------------------------------------//
+ //-----------------------------------------------------------------------------//
+
+
+ vector<double> lambda;
+ for(int i = 0; i < (int) (FissilArray.size() + FertilArray.size()); i++)
+ lambda.push_back(0);
+
+
+ double Na = 6.02214129e23; //N Avogadro
+
+ IsotopicVector FullUsedStock;
+ IsotopicVector stock;
+
+ bool FuelBuild = false;
+
+ int N_FissilStock_OnCheck = 0;
+
+ while(!FuelBuild)
+ {
+ double nPu_0 = 0;
+ double MPu_0 = 0;
+ {
+ map<ZAI ,double>::iterator it;
+
+ map<ZAI ,double> isotopicquantity = FullUsedStock.GetSpeciesComposition(94).GetIsotopicQuantity();
+ for( it = isotopicquantity.begin(); it != isotopicquantity.end(); it++ )
+ nPu_0 += (*it).second;
+
+ isotopicquantity = (FullUsedStock.GetSpeciesComposition(94) + ZAI(94,241,0)*FullUsedStock.GetZAIIsotopicQuantity(95,241,0)).GetIsotopicQuantity(); //Add the 241Am as 241Pu... the Pu is not old in the Eq Model but is in the FissileArray....;
+ for( it = isotopicquantity.begin(); it != isotopicquantity.end(); it++ )
+ MPu_0 += (*it).second*cZAIMass.fZAIMass.find( (*it).first )->second/Na*1e-6;
+ }
+
+ stock = FissilArray[N_FissilStock_OnCheck];
+ double nPu_1 = 0;
+ double MPu_1 = 0;
+ double Sum_AlphaI_nPuI = 0;
+ double Sum_AlphaI_nPuI0 = 0;
+ {
+ map<ZAI ,double>::iterator it;
+ map<ZAI ,double> isotopicquantity = stock.GetSpeciesComposition(94).GetIsotopicQuantity();
+
+ for( it = isotopicquantity.begin(); it != isotopicquantity.end(); it++ )
+ {
+ if ((*it).first.A() >= 238 && (*it).first.A() <= 242)
+ {
+ nPu_1 += (*it).second;
+ Sum_AlphaI_nPuI += fFuelParameter[(*it).first.A() -237]*(*it).second;
+ }
+ }
+
+ isotopicquantity = (stock.GetSpeciesComposition(94) + ZAI(94,241,0)*stock.GetZAIIsotopicQuantity(95,241,0)).GetIsotopicQuantity(); //Add the 241Am as 241Pu... the Pu is not old in the Eq Model but is in the FissileArray....
+ for( it = isotopicquantity.begin(); it != isotopicquantity.end(); it++ )
+ if ((*it).first.A() >= 238 && (*it).first.A() <= 242)
+ {
+ MPu_1 += (*it).second * (cZAIMass.fZAIMass.find( (*it).first )->second)/Na*1e-6;
+ }
+
+ isotopicquantity = FullUsedStock.GetSpeciesComposition(94).GetIsotopicQuantity();
+ for( it = isotopicquantity.begin(); it != isotopicquantity.end(); it++ )
+ if ((*it).first.A() >= 238 && (*it).first.A() <= 242)
+ {
+ Sum_AlphaI_nPuI0 += fFuelParameter[(*it).first.A() -237]*(*it).second;
+ }
+ }
+
+ double StockFactionToUse = 0;
+
+ double NT = HMMass*1e6 * Na / (cZAIMass.fZAIMass.find( ZAI(92,238,0) )->second*0.997 + cZAIMass.fZAIMass.find( ZAI(92,235,0) )->second*0.003 );
+
+ double N1 = (BurnUp - fFuelParameter[6]) * NT;
+ double N2 = -Sum_AlphaI_nPuI0;
+ double N3 = -fFuelParameter[0] * Na / (cZAIMass.fZAIMass.find( ZAI(92,238,0) )->second*0.997 + cZAIMass.fZAIMass.find( ZAI(92,235,0) )->second*0.003 ) * (HMMass*1e6 - MPu_0*1e6);
+
+ double D1 = Sum_AlphaI_nPuI;
+ double D2 = -fFuelParameter[0] * MPu_1*1e6 * Na / (cZAIMass.fZAIMass.find( ZAI(92,238,0) )->second*0.997 + cZAIMass.fZAIMass.find( ZAI(92,235,0) )->second*0.003 ) ;
+
+ StockFactionToUse = (N1 + N2 + N3) / (D1 + D2);
+
+ if(StockFactionToUse < 0)
+ {
+ cout << "!!Bad Trouble!! !!!FabricationPlant!!! Oups Bug in calculating stock fraction to use "<< endl;
+ GetLog()->fLog << "!!Bad Trouble!! !!!FabricationPlant!!! Oups Bug in calculating stock fraction to use" << endl;
+ lambda[N_FissilStock_OnCheck] = 0.;
+ FuelBuild = false;
+ }
+ else if( StockFactionToUse > 1 )
+ {
+
+ FullUsedStock += stock;
+ lambda[N_FissilStock_OnCheck] = 1;
+ N_FissilStock_OnCheck++;
+ FuelBuild = false;
+
+
+ if( N_FissilStock_OnCheck > (int) lambda.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;
+ }
+ }
+ else
+ {
+ lambda[N_FissilStock_OnCheck] = StockFactionToUse;
+
+ FuelBuild = true;
+
+ double U8_Quantity = (HMMass - (MPu_0+StockFactionToUse*MPu_1 ))/(cZAIMass.fZAIMass.find( ZAI(92,238,0) )->second*0.997 + cZAIMass.fZAIMass.find( ZAI(92,235,0) )->second*0.003 )*Na/1e-6;
+
+ lambda.back() = U8_Quantity / FertilArray[0].GetSumOfAll();
+ }
+
+ }
+
+
+
+ return lambda;
+}
diff --git a/source/branches/CLASSV3/Model/Equivalence/EQM_LIN_PWR_MOX.hxx b/source/branches/CLASSV3/Model/Equivalence/EQM_LIN_PWR_MOX.hxx
new file mode 100644
index 0000000000000000000000000000000000000000..f61dce3d4b8e3cbc4d38bad843af91603dc9c90e
--- /dev/null
+++ b/source/branches/CLASSV3/Model/Equivalence/EQM_LIN_PWR_MOX.hxx
@@ -0,0 +1,47 @@
+#ifndef _EQM_LIN_PWR_MOX_HXX
+#define _EQM_LIN_PWR_MOX_HXX
+
+#include "EquivalenceModel.hxx"
+
+#include <string>
+
+/*!
+ \file
+ \brief Header file for EQM_MLP_MOX class.
+
+
+ @author BaM
+ @version 1.0
+ */
+
+using namespace std;
+
+//-----------------------------------------------------------------------------//
+/*!
+ Define a EQM_MLP_MOX.
+ The aim of these class is to constuct a fuel from an equivalence model
+ based on a Linear Eq Model BU = SUM (alpha_i*n_i}
+
+ @author BaM
+ @version 3.0
+ */
+//________________________________________________________________________
+
+class EQM_LIN_PWR_MOX : public EquivalenceModel
+{
+ public :
+
+ EQM_LIN_PWR_MOX(string WeightPath);
+ ~EQM_LIN_PWR_MOX();
+
+ vector<double> BuildFuel(double BurnUp, double HMMass, vector<IsotopicVector> FissilArray, vector<IsotopicVector> FertilArray );
+
+ private :
+
+ string fWeightPath;
+ vector<double> fFuelParameter;
+
+};
+
+#endif
+
diff --git a/source/branches/CLASSV3/Model/Equivalence/EQM_MLP_PWR_MOX.cxx b/source/branches/CLASSV3/Model/Equivalence/EQM_MLP_PWR_MOX.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..a2c57dde9d9b1eba3500a2a403a89e871da9887e
--- /dev/null
+++ b/source/branches/CLASSV3/Model/Equivalence/EQM_MLP_PWR_MOX.cxx
@@ -0,0 +1,316 @@
+#include "EQM_MLP_PWR_MOX.hxx"
+#include "LogFile.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"
+#include "TFile.h"
+
+#include "IsotopicVector.hxx"
+
+//________________________________________________________________________
+//
+// EQM_MLP_MOX
+//
+// Equivalenve Model based on multi layer perceptron from TMVA (root cern)
+// For REP MOX use
+//
+//________________________________________________________________________
+
+EQM_MLP_MOX::EQM_MLP_MOX(string TMVAWeightPath)
+{
+ fTMVAWeightPath = TMVAWeightPath;
+
+
+
+ ZAI U8(92,238,0);
+ ZAI U5(92,235,0);
+ double U5_enrich= 0.0025;
+
+ ZAI Pu6(94,236,0);
+ 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 Pu4(94,244,0);
+
+ fFissileList = Pu6*1+ Pu8*1+Pu9*1+Pu0*1+Pu1*1+Pu2*1+Pu4*1;
+ fFertileList = U5*U5_enrich + U8*(1-U5_enrich);
+}
+
+//________________________________________________________________________
+void EQM_MLP_MOX::CreateTMVAInputTree(IsotopicVector Fissil,IsotopicVector Fertil,double BurnUp)
+{
+ /******Create Input data tree to be interpreted by TMVA::Reader***/
+
+ TFile* InputFile = new TFile("./EQMTMP.root","RECREATE");
+ TTree* InputTree = new TTree("EQTMP", "EQTMP");
+ float Pu8 = 0;
+ float Pu9 = 0;
+ float Pu10 = 0;
+ float Pu11 = 0;
+ float Pu12 = 0;
+ float Am1 = 0;
+ float U5_enrichment = 0;
+ float BU = 0;
+
+ InputTree->Branch( "Pu8" ,&Pu8 ,"Pu8/F" );
+ InputTree->Branch( "Pu9" ,&Pu9 ,"Pu9/F" );
+ InputTree->Branch( "Pu10" ,&Pu10 ,"Pu10/F" );
+ InputTree->Branch( "Pu11" ,&Pu11 ,"Pu11/F" );
+ InputTree->Branch( "Pu12" ,&Pu12 ,"Pu12/F" );
+ InputTree->Branch( "Am1" ,&Am1 ,"Am1/F" );
+ InputTree->Branch( "U5_enrichment" ,&U5_enrichment ,"U5_enrichment/F" );
+ InputTree->Branch( "BU" ,&BU ,"BU/F" );
+
+
+
+ float U8 = Fertil.GetZAIIsotopicQuantity(92,238,0);
+ float U5 = Fertil.GetZAIIsotopicQuantity(92,235,0);
+ float U4 = Fertil.GetZAIIsotopicQuantity(92,234,0);
+
+ float UTOT = U8 + U5 + U4;
+
+ Pu8 = Fissil.GetZAIIsotopicQuantity(94,238,0);
+ Pu9 = Fissil.GetZAIIsotopicQuantity(94,239,0);
+ Pu10 = Fissil.GetZAIIsotopicQuantity(94,240,0);
+ Pu11 = Fissil.GetZAIIsotopicQuantity(94,241,0);
+ Pu12 = Fissil.GetZAIIsotopicQuantity(94,242,0);
+ Am1 = Fissil.GetZAIIsotopicQuantity(95,241,0);
+
+ double TOTPU=(Pu8+Pu9+Pu10+Pu11+Pu12+Am1);
+
+ Pu8 = Pu8 / TOTPU;
+ Pu9 = Pu9 / TOTPU;
+ Pu10= Pu10 / TOTPU;
+ Pu11= Pu11 / TOTPU;
+ Pu12= Pu12 / TOTPU;
+ Am1 = Am1 / TOTPU;
+
+ U5_enrichment = U5 / UTOT;
+
+ BU=BurnUp;
+
+ 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 ...
+ {
+ cout<<"!!!!!!!!!!!ERRORR!!!!!!!!!!!!"<<endl;
+ cout<<Pu8<<" "<<Pu9<<" "<<Pu10<<" "<<Pu11<<" "<<Pu12<<" "<<Am1<<endl;
+ exit(0);
+ }
+ // All value are molar (!weight)
+
+ InputTree->Fill();
+
+ InputFile->Write();
+ delete InputTree;
+ InputFile-> Close();
+ delete InputFile;
+}
+//________________________________________________________________________
+double EQM_MLP_MOX::ExecuteTMVA()
+{
+ // --- 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,teneur;
+ 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( "BU" ,&BU );
+
+
+ // --- Book the MVA methods
+
+ // Book method MLP
+ TString methodName = "MLP method";
+ reader->BookMVA( methodName, fTMVAWeightPath );
+
+ // Prepare input tree
+ TFile *input(0);
+ if (!gSystem->AccessPathName( "./EQMTMP.root" )) {
+ input = TFile::Open( "./EQMTMP.root" ); // check if file in local directory exists
+ }
+ if (!input) {
+ std::cout << "ERROR: could not open data file" << std::endl;
+ exit(1);
+ }
+
+ TTree* theTree = (TTree*)input->Get("EQTMP");
+
+ theTree->SetBranchAddress("teneur",&teneur);
+ 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( "BU" ,&BU );
+
+ theTree->GetEntry(0);
+ Float_t val = (reader->EvaluateRegression( methodName ))[0];
+
+ delete reader;
+ delete theTree;
+ input->Close();
+ delete input;
+ //cout<<"....done"<<endl;
+
+ //cout<<val<<endl;
+ return (double)val; //retourne teneur
+}
+//________________________________________________________________________
+void EQM_MLP_MOX::GuessLambda(double& lambda, int& StockID,int FirstStockID, int LastStockID, double DeltaM,double StockHM)
+{
+
+ double Threshold = 50 ; //50 grams
+
+
+ /****Initialization***/
+ if(lambda==0)
+ lambda=0.5;
+
+ /********dichotomie**************/
+ if(DeltaM>0) //MassNeeded - AvailableMass
+ {
+ lambda+=lambda/2.;
+
+ if(lambda >= 1) //this stock is not enought go to next one
+ {
+ lambda = 1;
+
+ if( !( StockID+1 == LastStockID) ) //if its the last stock don't try to look the next one (segfault otherwize)
+ StockID++;
+ }
+
+ }
+
+ else if(DeltaM<0)
+ {
+ lambda-=lambda/2.;
+
+ if(lambda*StockHM < Threshold) //if only 50 grams left in actual stock go to previous one
+ {
+ lambda=0;
+ StockID--;
+ if(StockID<FirstStockID)
+ {
+ cout<<"Critical error EQM_MLP_MOX::GuessLambda"<<endl;
+ cout<<"Contact BLG"<<endl;
+ exit(1);
+ }
+
+ }
+
+ }
+
+}
+
+//________________________________________________________________________
+vector<double> EQM_MLP_MOX::BuildFuel(double BurnUp, double HMMass,vector<IsotopicVector> FissilArray, vector<IsotopicVector> FertilArray)
+{
+
+ HMMass*=1e6;//tons to gram
+
+ vector<double> lambda ; //vector of portion of stocks taken (fissile & fertil)
+ for(int i=0;FissilArray.size() + FertilArray.size();i++ );
+ lambda.push_back(0);
+
+ /*******************Depleted Uranium Vector**************************/
+ /* ZAI U5(92,235,0);
+ ZAI U8(92,238,0);
+ double U5_enrich= 0.0025;
+ double MeanMolarDepletedU = U5.GetMass()*U5_enrich + (1-U5_enrich)*U8.GetMass();
+ double AVOGADRO = 6.02214129e23;
+ //Building a isotopic vector with a total mass of HMMass (we are going to take just a share of it but want to be sure we have enought)
+ DepletedUranium.Add(U5, U5_enrich *HMMass/MeanMolarDepletedU*AVOGADRO);
+ DepletedUranium.Add(U8, (1-U5_enrich)*HMMass/MeanMolarDepletedU*AVOGADRO);
+ */
+ IsotopicVector DepletedUranium;
+ IsotopicVector PlutoniumVector;
+ int StockID = 0 ;
+ double AvailablePuMass=0;
+ double PuMassNeeded=1000; //At present time I have no clue what is the requiered Pu mass, I assume at least 1kg is needed
+ double WeightPuContent=0;
+ double MassPrecision=100; //Mass precision is 100 grams
+ while( abs(PuMassNeeded - AvailablePuMass) > MassPrecision )
+ {
+ //Increase the portion of the stock stokID taken, according to the followings variables
+ double DeltaM=PuMassNeeded-AvailablePuMass;
+ GuessLambda( lambda[StockID], StockID,0,FissilArray.size()-1, DeltaM, FissilArray[StockID].GetTotalMass() * 1e6 );
+
+ //Build the Plutonium vector from stocks
+ PlutoniumVector.Clear();
+ for(int i=0;i<=StockID;i++)
+ PlutoniumVector += lambda[i] * FissilArray[i];
+
+ AvailablePuMass = PlutoniumVector.GetTotalMass() * 1e6; //in grams
+
+ //Build uranium vector from stocks
+ int FertileStockID = FissilArray.size();
+ double FertilMassNeeded = HMMass - AvailablePuMass;
+ double AvailableFertilMass = 0;
+
+ while( abs(FertilMassNeeded - AvailableFertilMass) > MassPrecision )
+ {
+ double DeltaM=FertilMassNeeded-AvailableFertilMass;
+ GuessLambda( lambda[FertileStockID], FertileStockID,FissilArray.size(),FertilArray.size()-1, DeltaM, FertilArray[FertileStockID-FissilArray.size()].GetTotalMass() * 1e6 ) ;
+
+ int j=-1;
+ DepletedUranium.Clear();
+ for(int i=int(FissilArray.size());i<=FertileStockID;i++)
+ { DepletedUranium += lambda[i] * FertilArray[j];
+ j++;
+ }
+ AvailableFertilMass=DepletedUranium.GetTotalMass() * 1e6; //in grams
+
+ if( j+1 == int( FertilArray.size() ) && FertilMassNeeded > AvailableFertilMass ) //if this is the last stock and it's not enought
+ {
+ cout<<"You requiere more depleted uranium "<<"("<<DeltaM/1e6<<" t needed) ! Reactor not fill"<<endl;
+ for(int i=0 ; i<int(lambda.size()) ; i++)
+ lambda[i]=-1;
+ break;
+ }
+
+ }
+ /*Calcul the quantity of this composition needed to reach the burnup*/
+ CreateTMVAInputTree(PlutoniumVector,DepletedUranium,BurnUp);
+ double MolarPuContent = ExecuteTMVA();
+ system("rm EQMTMP.root ");
+
+ double MeanMolarPu = PlutoniumVector.MeanMolar();
+ double MeanMolarDepletedU = DepletedUranium.MeanMolar();
+ double MeanMolar = MeanMolarPu*MolarPuContent + (1-MolarPuContent)*MeanMolarDepletedU;
+
+ WeightPuContent = MolarPuContent * MeanMolarPu/MeanMolar ;
+
+ PuMassNeeded = WeightPuContent * HMMass ;
+
+ if( StockID+1 == int( FissilArray.size() ) && PuMassNeeded > AvailablePuMass ) //if this is the last stock and it's not enought
+ {
+ cout<<"You requiere more (or better) plutonium !! Reactor not fill"<<endl;
+ for(int i=0 ; i<int(lambda.size()) ; i++)
+ lambda[i]=-1;
+ break;
+ }
+
+ }
+
+return lambda;
+}
+//________________________________________________________________________
diff --git a/source/branches/CLASSV3/Model/Equivalence/EQM_MLP_PWR_MOX.hxx b/source/branches/CLASSV3/Model/Equivalence/EQM_MLP_PWR_MOX.hxx
new file mode 100644
index 0000000000000000000000000000000000000000..aa5eb0e50236e9b3d95138a41e4dad7717705f97
--- /dev/null
+++ b/source/branches/CLASSV3/Model/Equivalence/EQM_MLP_PWR_MOX.hxx
@@ -0,0 +1,60 @@
+#ifndef _EQM_MLP_MOX_HXX
+#define _EQM_MLP_MOX_HXX
+
+#include "EquivalenceModel.hxx"
+
+/*!
+ \file
+ \brief Header file for EQM_MLP_MOX class.
+
+
+ @author BaM
+ @version 1.0
+ */
+
+
+using namespace std;
+
+//-----------------------------------------------------------------------------//
+/*!
+ Define a EQM_MLP_MOX.
+ The aim of these class is to constuct a fuel from an equivalence model
+ based on a Multi layer perceptron
+
+ @author BaM
+ @version 3.0
+ */
+//________________________________________________________________________
+
+
+class EQM_MLP_MOX : public EquivalenceModel
+{
+ public :
+
+ EQM_MLP_MOX(string TMVAWeightPath);
+
+ /// method called to build a reprocessed fuel as a function of the burnup requierement the stock, mass....
+ /*!
+ Build the fuel following the equivalance model with the proper requierment in term of mass burnup....
+ \param double BurnUp desireted burnup reached by the fuel at the end of irradiation
+ \param double HMMass, needed Heavy metal mass needed
+ \param vector<double> &lambda, fraction of the stock to take (initialy should be 0)
+ \param vector<IsotopicVector> FissilArray, isotopicvectors to use to get the fissil part of the fuel
+ \param vector<IsotopicVector> FertilArray, isotopicvectors to use to get the fertil part of the fuel (if empty take it from the god)
+ */
+
+ vector<double> BuildFuel(double BurnUp, double HMMass,vector<IsotopicVector> FissilArray, vector<IsotopicVector> FertilArray = vector<IsotopicVector>());
+ //}
+ void GuessLambda(double& lambda, int& StockID,int FirstStockID, int LastStockID, double DeltaM,double StockHM);
+
+ private :
+ void CreateTMVAInputTree(IsotopicVector Fissil,IsotopicVector Fertil,double BurnUp);
+ double ExecuteTMVA();
+
+
+ string fTMVAWeightPath;
+
+};
+
+#endif
+
diff --git a/source/branches/CLASSV3/Model/Equivalence/EQM_QUAD_PWR_MOX.cxx b/source/branches/CLASSV3/Model/Equivalence/EQM_QUAD_PWR_MOX.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..aba54ecba9607ee6b90cc9ed7de1fd6476606ded
--- /dev/null
+++ b/source/branches/CLASSV3/Model/Equivalence/EQM_QUAD_PWR_MOX.cxx
@@ -0,0 +1,93 @@
+#include "EQM_QUAD_PWR_MOX.hxx"
+
+#include <vector>
+
+#include "StringLine.hxx"
+#include "LogFile.hxx"
+
+
+
+
+EQM_QUAD_PWR_MOX::EQM_QUAD_PWR_MOX(string WeightPath):EquivalenceModel()
+{
+ fWeightPath = WeightPath;
+
+ ifstream DataDB(fWeightPath.c_str()); // Open the File
+ if(!DataDB)
+ cout << "!!Warning!! !!!EQM QUAD PWR MOX!!! \n Can't open \"" << fWeightPath << "\"\n" << endl;
+
+ string line;
+ int start = 0; // First Get Fuel Parameter
+ getline(DataDB, line);
+
+ if( StringLine::NextWord(line, start, ' ') != "PARAM")
+ {
+ cout << "!!Bad Trouble!! !!!EQM QUAD PWR MOX!!! 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()));
+
+ cout << "!!INFO!! !!!EQM QUAD PWR MOX!!! " << fFuelParameter.size() << " have been read"<< endl;
+
+
+ ZAI U8(92,238,0);
+ ZAI U5(92,235,0);
+ double U5_enrich= 0.0025;
+ fFertileList = U5*U5_enrich + U8*(1-U5_enrich);
+
+
+ ZAI Pu8(94,238,0);
+ ZAI Pu9(94,239,0);
+ ZAI Pu0(94,240,0);
+ ZAI Pu1(94,241,0);
+ ZAI Pu2(94,242,0);
+ fFissileList = Pu8*1+Pu9*1+Pu0*1+Pu1*1+Pu2*1;
+
+
+}
+
+
+EQM_QUAD_PWR_MOX::~EQM_QUAD_PWR_MOX()
+{
+
+}
+
+
+
+
+double EQM_QUAD_PWR_MOX::GetFissileMolarFraction(IsotopicVector Fissile,IsotopicVector Fertile,double BurnUp)
+{
+
+
+ ZAI ZAIList[6] = {ZAI(94,238,0), ZAI(94,239,0), ZAI(94,240,0), ZAI(94,241,0), ZAI(94,242,0), ZAI(95,241,0) };
+
+ vector<double> PuCompo;
+ double Sum = Fissile.GetSumOfAll();
+
+
+ for(int i = 0; i< 5; i++)
+ PuCompo.push_back( Fissile.GetZAIIsotopicQuantity(ZAIList[i])/Sum *100 );
+
+ PuCompo[2] += Fissile.GetZAIIsotopicQuantity(ZAIList[5])/Sum *100;
+ double A = 0;
+
+ if(PuCompo[0] <= PuCompo[2] && PuCompo[0] <= PuCompo[4] && PuCompo[1] + PuCompo[3] >= 40 && PuCompo[1] >0 )
+ {
+ int par = 0;
+ for(int j = 0 ; j < 5 ; j++)
+ {
+ A += fFuelParameter[par] * PuCompo[j]/100 ;
+ par++;
+ for(int i = j ; i < 5 ; i++)
+ {
+ A += fFuelParameter[par] *PuCompo[i]/100 *PuCompo[j]/100;
+ par++;
+ }
+ }
+ A += fFuelParameter[par];
+ }
+
+ return A;
+}
+
diff --git a/source/branches/CLASSV3/Model/Equivalence/EQM_QUAD_PWR_MOX.hxx b/source/branches/CLASSV3/Model/Equivalence/EQM_QUAD_PWR_MOX.hxx
new file mode 100644
index 0000000000000000000000000000000000000000..9a13fca25583a74334fd5a8da93ac21533520864
--- /dev/null
+++ b/source/branches/CLASSV3/Model/Equivalence/EQM_QUAD_PWR_MOX.hxx
@@ -0,0 +1,47 @@
+#ifndef _EQM_QUAD_PWR_MOX_HXX
+#define _EQM_QUAD_PWR_MOX_HXX
+
+#include "EquivalenceModel.hxx"
+
+#include <string>
+
+/*!
+ \file
+ \brief Header file for EQM_MLP_MOX class.
+
+
+ @author BaM
+ @version 1.0
+ */
+
+using namespace std;
+
+//-----------------------------------------------------------------------------//
+/*!
+ Define a EQM_MLP_MOX.
+ The aim of these class is to constuct a fuel from an equivalence model
+ based on a Quadratic Pu equivalent Model
+
+ @author BaM
+ @version 3.0
+ */
+//________________________________________________________________________
+
+class EQM_QUAD_PWR_MOX : public EquivalenceModel
+{
+ public :
+
+ EQM_QUAD_PWR_MOX(string WeightPath);
+ ~EQM_QUAD_PWR_MOX();
+
+ double GetFissileMolarFraction(IsotopicVector Fissil,IsotopicVector Fertil,double BurnUp);
+
+ private :
+
+ string fWeightPath;
+ vector<double> fFuelParameter;
+
+};
+
+#endif
+
diff --git a/source/branches/CLASSV3/Model/Irradiation/IM_Matrix.cxx b/source/branches/CLASSV3/Model/Irradiation/IM_Matrix.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..f10b124b9f50d448ce732da5669875c8b43fa92e
--- /dev/null
+++ b/source/branches/CLASSV3/Model/Irradiation/IM_Matrix.cxx
@@ -0,0 +1,317 @@
+//
+// 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 "CLASSHeaders.hxx"
+#include "LogFile.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():DynamicalSystem()
+{
+ fShorstestHalflife = 3600.*24*160.; //cut by default all nuclei with a shorter liftime than the Cm242 -> remain 33 actinides
+}
+
+
+IM_Matrix::IM_Matrix(LogFile* log):IrradiationModel(log), DynamicalSystem()
+{
+ fShorstestHalflife = 3600.*24*160.; //cut by default all nuclei with a shorter liftime than the Cm242 -> remain 33 actinides
+}
+
+
+
+
+//________________________________________________________________________
+/* Evolution Calculation */
+//________________________________________________________________________
+EvolutionData IM_Matrix::GenerateEvolutionData(IsotopicVector isotopicvector, EvolutionData XSSet, double Power, double cycletime)
+{
+
+ if(fFastDecay.size() == 0)
+ {
+ BuildDecayMatrix();
+ fNVar = findex_inver.size();
+ }
+
+ string ReactorType;
+
+
+ vector< TMatrixT<double> > NMatrix ;// TMatrixT<double>(decayindex.size(),1))
+ { // Filling the t=0 State;
+ map<ZAI, double > isotopicquantity = isotopicvector.GetIsotopicQuantity();
+ TMatrixT<double> N_0Matrix = TMatrixT<double>( findex.size(),1) ;
+ for(int i = 0; i < (int)findex.size(); i++)
+ N_0Matrix[i] = 0;
+
+ map<ZAI, double >::iterator it ;
+ for(int i = 0; i < (int)findex.size(); i++)
+ N_0Matrix[i] = 0;
+
+ for(it = isotopicquantity.begin(); it != isotopicquantity.end(); it++)
+ {
+ /// Need TO change with FP managment
+ map<ZAI, int >::iterator it2;
+
+ if( (*it).first.Z() < fZAIThreshold )
+ it2 = findex_inver.find( ZAI(-2,-2,-2) );
+ else it2 = findex_inver.find( (*it).first );
+
+ if(it2 == findex_inver.end() ) //If not in index should be TMP, can't be fast decay for new Fuel !!!
+ it2 = findex_inver.find( ZAI(-3,-3,-3) );
+ N_0Matrix[ (*it2).second ][0] = (*it).second ;
+
+
+ }
+
+ isotopicquantity.clear();
+
+ NMatrix.push_back(N_0Matrix);
+ N_0Matrix.Clear();
+
+ }
+
+
+ //-------------------------//
+ //--- Perform Evolution ---//
+ //-------------------------//
+ ReactorType = XSSet.GetReactorType();
+
+ double Na = 6.02214129e23; //N Avogadro
+ double M_ref = XSSet.GetHeavyMetalMass();
+ double M = 0;
+ double Power_ref = XSSet.GetPower();
+
+ // Get the mass of the fuel to irradiate in order to perform the evolution at fixed burnup (the burnup step of the calculation will match the burnup step of the XSSet
+ {
+ map<ZAI, double >::iterator it ;
+
+
+ IsotopicVector IVtmp = isotopicvector.GetActinidesComposition();
+ map<ZAI, double >isotopicquantity = IVtmp.GetIsotopicQuantity();
+
+ for( it = isotopicquantity.begin(); it != isotopicquantity.end(); it++ )
+ M += isotopicvector.GetActinidesComposition().GetZAIIsotopicQuantity( (*it).first )*cZAIMass.fZAIMass.find( (*it).first )->second/Na*1e-6;
+ isotopicquantity.clear();
+
+ }
+
+ int NStep = XSSet.GetFissionXS().begin()->second->GetN();
+ double* DBTimeStep = XSSet.GetFissionXS().begin()->second->GetX();
+
+ int InsideStep = 10;
+
+ double timevector[NStep];
+ timevector[0] = 0;
+
+ double Flux[NStep];
+
+ TMatrixT<double> FissionEnergy = TMatrixT<double>(findex.size(),1);
+ for(int i = 0; i < (int)findex.size(); i++)
+ FissionEnergy[i] = 0;
+
+ {
+ map< ZAI, int >::iterator it;
+ for(it = findex_inver.begin(); it != findex_inver.end(); it++)
+ {
+ map< ZAI, double >::iterator it2 = fFissionEnergy.find(it->first);
+ if(it2 == fFissionEnergy.end())
+ {
+ if(it->first.Z() > fZAIThreshold)
+ FissionEnergy[it->second][0] = 1.9679e6*it->first.A()-2.601e8; // //simple linear fit to known values ;extrapolation to unknown isotopes
+ else FissionEnergy[it->second][0] = 0;
+ }
+ else
+ FissionEnergy[it->second][0] = it2->second;
+
+ }
+ }
+
+ vector< TMatrixT<double> > FissionXSMatrix; // Store The Fisison XS Matrix
+ vector< TMatrixT<double> > CaptureXSMatrix; // Store The Capture XS Matrix
+ vector< TMatrixT<double> > n2nXSMatrix; // Store The n2N XS Matrix
+
+ for(int i = 0; i < NStep-1; i++)
+ {
+ double TStepMax = ( (DBTimeStep[i+1]-DBTimeStep[i] ) ) * Power_ref/M_ref / Power*M ; // Get the next Time step
+
+
+ TMatrixT<double> BatemanMatrix = TMatrixT<double>(findex.size(),findex.size());
+ TMatrixT<double> BatemanReactionMatrix = TMatrixT<double>(findex.size(),findex.size());
+
+ TMatrixT<double> NEvolutionMatrix = TMatrixT<double>(findex.size(),1);
+ NEvolutionMatrix = NMatrix.back();
+
+
+
+ FissionXSMatrix.push_back(GetFissionXsMatrix(XSSet, DBTimeStep[i])); //Feel the fission reaction Matrix
+ CaptureXSMatrix.push_back(GetCaptureXsMatrix(XSSet, DBTimeStep[i])); //Feel the capture reaction Matrix
+ n2nXSMatrix.push_back(Getn2nXsMatrix(XSSet, DBTimeStep[i])); //Feel the (n,2n) reaction Matrix
+
+ // ---------------- Evolution
+
+ BatemanReactionMatrix = FissionXSMatrix[i];
+ BatemanReactionMatrix += CaptureXSMatrix[i];
+ BatemanReactionMatrix += n2nXSMatrix[i];
+
+ for(int k=0; k < InsideStep; k++)
+ {
+ double ESigmaN = 0;
+ for (int j = 0; j < (int)findex.size() ; j++)
+ ESigmaN -= FissionXSMatrix[i][j][j]*NEvolutionMatrix[j][0]*1.6e-19*FissionEnergy[j][0];
+ // Update Flux
+ double Flux_k = Power/ESigmaN;
+
+ if(k==0)
+ Flux[i]=Flux_k;
+
+ BatemanMatrix = BatemanReactionMatrix;
+ BatemanMatrix *= Flux_k;
+ BatemanMatrix += fDecayMatrix ;
+ BatemanMatrix *= TStepMax/InsideStep ;
+
+
+ TMatrixT<double> IdMatrix = TMatrixT<double>(findex.size(),findex.size());
+ for(int j = 0; j < (int)findex.size(); j++)
+ for(int k = 0; k < (int)findex.size(); k++)
+ {
+ if(k == j) IdMatrix[j][k] = 1;
+ else IdMatrix[j][k] = 0;
+ }
+
+
+ TMatrixT<double> BatemanMatrixDL = TMatrixT<double>(findex.size(),findex.size()); // Order 0 Term from the DL : Id
+ TMatrixT<double> BatemanMatrixDLTermN = TMatrixT<double>(findex.size(),findex.size()); // Addind it;
+
+ {
+ BatemanMatrix *= TStepMax ;
+ BatemanMatrixDLTermN = IdMatrix;
+ BatemanMatrixDL = BatemanMatrixDLTermN;
+ int j = 1;
+ double NormN;
+
+ do
+ {
+ TMatrixT<double> BatemanMatrixDLTermtmp = TMatrixT<double>(findex.size(),findex.size()); // Adding it;
+ BatemanMatrixDLTermtmp = BatemanMatrixDLTermN;
+
+ BatemanMatrixDLTermN.Mult(BatemanMatrixDLTermtmp, BatemanMatrix );
+
+ BatemanMatrixDLTermN *= 1./j;
+ BatemanMatrixDL += BatemanMatrixDLTermN;
+
+ NormN = 0;
+ for(int m = 0; m < (int)findex.size(); m++)
+ for(int n = 0; n < (int)findex.size(); n++)
+ NormN += BatemanMatrixDLTermN[m][n]*BatemanMatrixDLTermN[m][n];
+ j++;
+
+ } while ( NormN != 0 );
+ }
+ NEvolutionMatrix = BatemanMatrixDL * NEvolutionMatrix ;
+ }
+ NMatrix.push_back(NEvolutionMatrix);
+
+
+ timevector[i+1] = timevector[i] + TStepMax;
+
+ BatemanMatrix.Clear();
+ BatemanReactionMatrix.Clear();
+ NEvolutionMatrix.Clear();
+
+
+ }
+ FissionXSMatrix.push_back(GetFissionXsMatrix(XSSet, DBTimeStep[NStep-1])); //Feel the reaction Matrix
+ CaptureXSMatrix.push_back(GetCaptureXsMatrix(XSSet, DBTimeStep[NStep-1])); //Feel the reaction Matrix
+ n2nXSMatrix.push_back(Getn2nXsMatrix(XSSet, DBTimeStep[NStep-1])); //Feel the reaction Matrix
+
+
+ EvolutionData GeneratedDB = EvolutionData(GetLog());
+
+ double ESigmaN = 0;
+ for (int j = 0; j < (int)findex.size() ; j++)
+ ESigmaN -= FissionXSMatrix.back()[j][j]*NMatrix.back()[j][0]*1.6e-19*FissionEnergy[j][0];
+
+ Flux[NStep-1] = Power/ESigmaN;
+
+ GeneratedDB.SetFlux( new TGraph(NStep, timevector, Flux) );
+
+ for(int i = 0; i < (int)findex.size(); i++)
+ {
+ double ZAIQuantity[NMatrix.size()];
+ double FissionXS[NStep];
+ double CaptureXS[NStep];
+ double n2nXS[NStep];
+ for(int j = 0; j < (int)NMatrix.size(); j++)
+ ZAIQuantity[j] = (NMatrix[j])[i][0];
+
+ for(int j = 0; j < NStep; j++)
+ {
+ FissionXS[j] = FissionXSMatrix[j][i][i];
+ CaptureXS[j] = CaptureXSMatrix[j][i][i];
+ n2nXS[j] = n2nXSMatrix[j][i][i];
+ }
+
+ GeneratedDB.NucleiInsert(pair<ZAI, TGraph*> (findex.find(i)->second, new TGraph(NMatrix.size(), timevector, ZAIQuantity)));
+ GeneratedDB.FissionXSInsert(pair<ZAI, TGraph*> (findex.find(i)->second, new TGraph(NStep, timevector, FissionXS)));
+ GeneratedDB.CaptureXSInsert(pair<ZAI, TGraph*> (findex.find(i)->second, new TGraph(NStep, timevector, CaptureXS)));
+ GeneratedDB.n2nXSInsert(pair<ZAI, TGraph*> (findex.find(i)->second, new TGraph(NStep, timevector, n2nXS)));
+ }
+
+ GeneratedDB.SetPower(Power );
+ // GeneratedDB.SetFuelType(fFuelType );
+ 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();
+
+ return GeneratedDB;
+
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/source/branches/CLASSV3/Model/Irradiation/IM_Matrix.hxx b/source/branches/CLASSV3/Model/Irradiation/IM_Matrix.hxx
new file mode 100644
index 0000000000000000000000000000000000000000..3b14f64aebdd3ac1e959361330624127043f6629
--- /dev/null
+++ b/source/branches/CLASSV3/Model/Irradiation/IM_Matrix.hxx
@@ -0,0 +1,64 @@
+#ifndef _IMMATRIX_HXX
+#define _IMMATRIX_HXX
+
+
+/*!
+ \file
+ \brief Header file for IrradiationModel class.
+
+
+ @author BaM
+ @version 2.0
+ */
+#include "DynamicalSystem.hxx"
+#include "IrradiationModel.hxx"
+
+
+using namespace std;
+
+
+class LogFile;
+
+//-----------------------------------------------------------------------------//
+/*!
+ Define a IM_Matrix.
+ The aim of these class is perform the calculation of the evolution of a fuel trough irradiation solving numericaly the Bateman using Matrix.
+
+
+ @author BaM
+ @version 3.0
+ */
+//________________________________________________________________________
+
+class EvolutionData;
+
+class IM_Matrix : public IrradiationModel, DynamicalSystem
+{
+
+ public :
+
+ IM_Matrix();
+ IM_Matrix(LogFile* log);
+
+
+
+ /// virtueal 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
+ */
+ EvolutionData GenerateEvolutionData(IsotopicVector IV, EvolutionData XSSet, double Power, double cycletime);
+ //}
+
+
+
+
+ private :
+
+
+
+};
+
+#endif
+
diff --git a/source/branches/CLASSV3/Model/Irradiation/IM_RK4.cxx b/source/branches/CLASSV3/Model/Irradiation/IM_RK4.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..46000361cd19167bfbe07a1d1dbcb1ae1d54910a
--- /dev/null
+++ b/source/branches/CLASSV3/Model/Irradiation/IM_RK4.cxx
@@ -0,0 +1,410 @@
+//
+// IM_RK4.cxx
+// CLASSSource
+//
+// Created by BaM on 04/05/2014.
+// Copyright (c) 2014 BaM. All rights reserved.
+//
+
+#include "IM_RK4.hxx"
+
+#include "IsotopicVector.hxx"
+#include "CLASSHeaders.hxx"
+#include "LogFile.hxx"
+#include "StringLine.hxx"
+
+#include <TGraph.h>
+#include <TString.h>
+
+
+#include <sstream>
+#include <string>
+#include <iostream>
+#include <fstream>
+#include <algorithm>
+#include <cmath>
+
+
+
+using namespace std;
+
+
+//________________________________________________________________________
+IM_RK4::IM_RK4():DynamicalSystem(), IrradiationModel()
+{
+ fTheNucleiVector = 0;
+ fTheMatrix = 0;
+
+ SetForbidNegativeValue();
+
+}
+
+
+IM_RK4::IM_RK4(LogFile* log):DynamicalSystem(), IrradiationModel(log)
+{
+
+ fTheNucleiVector = 0;
+ fTheMatrix = 0;
+
+ SetForbidNegativeValue();
+
+}
+
+
+
+
+
+
+
+
+
+//________________________________________________________________________
+/* Evolution Calculation */
+//________________________________________________________________________
+EvolutionData IM_RK4::GenerateEvolutionData(IsotopicVector isotopicvector, EvolutionData XSSet, double Power, double cycletime)
+{
+
+ if(fFastDecay.size() == 0)
+ {
+ BuildDecayMatrix();
+ fNVar = findex_inver.size();
+ }
+
+ SetTheMatrixToZero();
+ SetTheNucleiVectorToZero();
+
+ string ReactorType;
+
+
+ vector< TMatrixT<double> > NMatrix ;// TMatrixT<double>(decayindex.size(),1))
+ { // Filling the t=0 State;
+ map<ZAI, double > isotopicquantity = isotopicvector.GetIsotopicQuantity();
+ TMatrixT<double> N_0Matrix = TMatrixT<double>( findex.size(),1) ;
+ for(int i = 0; i < (int)findex.size(); i++)
+ N_0Matrix[i] = 0;
+
+ map<ZAI, double >::iterator it ;
+ for(int i = 0; i < (int)findex.size(); i++)
+ N_0Matrix[i] = 0;
+
+ for(it = isotopicquantity.begin(); it != isotopicquantity.end(); it++)
+ {
+ /// Need TO change with FP managment
+ map<ZAI, int >::iterator it2;
+
+ if( (*it).first.Z() < fZAIThreshold )
+ it2 = findex_inver.find( ZAI(-2,-2,-2) );
+ else it2 = findex_inver.find( (*it).first );
+
+ if(it2 == findex_inver.end() ) //If not in index should be TMP, can't be fast decay for new Fuel !!!
+ it2 = findex_inver.find( ZAI(-3,-3,-3) );
+ N_0Matrix[ (*it2).second ][0] = (*it).second ;
+
+
+ }
+
+ isotopicquantity.clear();
+
+ NMatrix.push_back(N_0Matrix);
+ N_0Matrix.Clear();
+
+ }
+
+
+ //-------------------------//
+ //--- Perform Evolution ---//
+ //-------------------------//
+ ReactorType = XSSet.GetReactorType();
+
+ double Na = 6.02214129e23; //N Avogadro
+ double M_ref = XSSet.GetHeavyMetalMass();
+ double M = 0;
+ double Power_ref = XSSet.GetPower();
+
+ // Get the mass of the fuel to irradiate in order to perform the evolution at fixed burnup (the burnup step of the calculation will match the burnup step of the XSSet
+ {
+ map<ZAI, double >::iterator it ;
+
+
+ IsotopicVector IVtmp = isotopicvector.GetActinidesComposition();
+ map<ZAI, double >isotopicquantity = IVtmp.GetIsotopicQuantity();
+
+ for( it = isotopicquantity.begin(); it != isotopicquantity.end(); it++ )
+ M += isotopicvector.GetActinidesComposition().GetZAIIsotopicQuantity( (*it).first )*cZAIMass.fZAIMass.find( (*it).first )->second/Na*1e-6;
+ isotopicquantity.clear();
+
+ }
+
+ int NStep = XSSet.GetFissionXS().begin()->second->GetN();
+ double* DBTimeStep = XSSet.GetFissionXS().begin()->second->GetX();
+
+ int InsideStep = 10;
+
+ double timevector[NStep];
+ timevector[0] = 0;
+
+ double Flux[NStep];
+
+ TMatrixT<double> FissionEnergy = TMatrixT<double>(findex.size(),1);
+ for(int i = 0; i < (int)findex.size(); i++)
+ FissionEnergy[i] = 0;
+
+ {
+ map< ZAI, int >::iterator it;
+ for(it = findex_inver.begin(); it != findex_inver.end(); it++)
+ {
+ map< ZAI, double >::iterator it2 = fFissionEnergy.find(it->first);
+ if(it2 == fFissionEnergy.end())
+ {
+ if(it->first.Z() > fZAIThreshold)
+ FissionEnergy[it->second][0] = 1.9679e6*it->first.A()-2.601e8; // //simple linear fit to known values ;extrapolation to unknown isotopes
+ else FissionEnergy[it->second][0] = 0;
+ }
+ else
+ FissionEnergy[it->second][0] = it2->second;
+
+ }
+ }
+
+ vector< TMatrixT<double> > FissionXSMatrix; // Store The Fisison XS Matrix
+ vector< TMatrixT<double> > CaptureXSMatrix; // Store The Capture XS Matrix
+ vector< TMatrixT<double> > n2nXSMatrix; // Store The n2N XS Matrix
+
+ for(int i = 0; i < NStep-1; i++)
+ {
+ double TStepMax = ( (DBTimeStep[i+1]-DBTimeStep[i] ) ) * Power_ref/M_ref / Power*M ; // Get the next Time step
+
+
+ TMatrixT<double> BatemanMatrix = TMatrixT<double>(findex.size(),findex.size());
+ TMatrixT<double> BatemanReactionMatrix = TMatrixT<double>(findex.size(),findex.size());
+
+ TMatrixT<double> NEvolutionMatrix = TMatrixT<double>(findex.size(),1);
+ NEvolutionMatrix = NMatrix.back();
+
+
+
+ FissionXSMatrix.push_back(GetFissionXsMatrix(XSSet, DBTimeStep[i])); //Feel the fission reaction Matrix
+ CaptureXSMatrix.push_back(GetCaptureXsMatrix(XSSet, DBTimeStep[i])); //Feel the capture reaction Matrix
+ n2nXSMatrix.push_back(Getn2nXsMatrix(XSSet, DBTimeStep[i])); //Feel the (n,2n) reaction Matrix
+
+ // ---------------- Evolution
+
+ BatemanReactionMatrix = FissionXSMatrix[i];
+ BatemanReactionMatrix += CaptureXSMatrix[i];
+ BatemanReactionMatrix += n2nXSMatrix[i];
+
+ for(int k=0; k < InsideStep; k++)
+ {
+ double ESigmaN = 0;
+ for (int j = 0; j < (int)findex.size() ; j++)
+ ESigmaN -= FissionXSMatrix[i][j][j]*NEvolutionMatrix[j][0]*1.6e-19*FissionEnergy[j][0];
+ // Update Flux
+ double Flux_k = Power/ESigmaN;
+
+ if(k==0)
+ Flux[i]=Flux_k;
+
+ BatemanMatrix = BatemanReactionMatrix;
+ BatemanMatrix *= Flux_k;
+ BatemanMatrix += fDecayMatrix ;
+
+ SetTheMatrixToZero();
+ SetTheNucleiVectorToZero();
+
+ SetTheMatrix(BatemanMatrix);
+ SetTheNucleiVector(NEvolutionMatrix);
+
+
+ RungeKutta(fTheNucleiVector, timevector[i]+TStepMax/InsideStep*k, timevector[i]+TStepMax/InsideStep*(k+1), fNVar);
+ NEvolutionMatrix = GetTheNucleiVector();
+
+ }
+ NEvolutionMatrix = GetTheNucleiVector();
+ NMatrix.push_back(NEvolutionMatrix);
+
+ timevector[i+1] = timevector[i] + TStepMax;
+
+ BatemanMatrix.Clear();
+ BatemanReactionMatrix.Clear();
+ NEvolutionMatrix.Clear();
+
+
+ }
+ FissionXSMatrix.push_back(GetFissionXsMatrix(XSSet, DBTimeStep[NStep-1])); //Feel the reaction Matrix
+ CaptureXSMatrix.push_back(GetCaptureXsMatrix(XSSet, DBTimeStep[NStep-1])); //Feel the reaction Matrix
+ n2nXSMatrix.push_back(Getn2nXsMatrix(XSSet, DBTimeStep[NStep-1])); //Feel the reaction Matrix
+
+
+ EvolutionData GeneratedDB = EvolutionData(GetLog());
+
+ double ESigmaN = 0;
+ for (int j = 0; j < (int)findex.size() ; j++)
+ ESigmaN -= FissionXSMatrix.back()[j][j]*NMatrix.back()[j][0]*1.6e-19*FissionEnergy[j][0];
+
+ Flux[NStep-1] = Power/ESigmaN;
+
+ GeneratedDB.SetFlux( new TGraph(NStep, timevector, Flux) );
+
+ for(int i = 0; i < (int)findex.size(); i++)
+ {
+ double ZAIQuantity[NMatrix.size()];
+ double FissionXS[NStep];
+ double CaptureXS[NStep];
+ double n2nXS[NStep];
+ for(int j = 0; j < (int)NMatrix.size(); j++)
+ ZAIQuantity[j] = (NMatrix[j])[i][0];
+
+ for(int j = 0; j < NStep; j++)
+ {
+ FissionXS[j] = FissionXSMatrix[j][i][i];
+ CaptureXS[j] = CaptureXSMatrix[j][i][i];
+ n2nXS[j] = n2nXSMatrix[j][i][i];
+ }
+
+ GeneratedDB.NucleiInsert(pair<ZAI, TGraph*> (findex.find(i)->second, new TGraph(NMatrix.size(), timevector, ZAIQuantity)));
+ GeneratedDB.FissionXSInsert(pair<ZAI, TGraph*> (findex.find(i)->second, new TGraph(NStep, timevector, FissionXS)));
+ GeneratedDB.CaptureXSInsert(pair<ZAI, TGraph*> (findex.find(i)->second, new TGraph(NStep, timevector, CaptureXS)));
+ GeneratedDB.n2nXSInsert(pair<ZAI, TGraph*> (findex.find(i)->second, new TGraph(NStep, timevector, n2nXS)));
+ }
+
+ GeneratedDB.SetPower(Power );
+// GeneratedDB.SetFuelType(fFuelType );
+ GeneratedDB.SetReactorType(ReactorType );
+ GeneratedDB.SetCycleTime(cycletime);
+
+ ResetTheMatrix();
+ ResetTheNucleiVector();
+
+ for (int i = 0; i < (int) FissionXSMatrix.size(); i++)
+ {
+ FissionXSMatrix[i].Clear();
+ CaptureXSMatrix[i].Clear();
+ n2nXSMatrix[i].Clear();
+ }
+ FissionXSMatrix.clear();
+ CaptureXSMatrix.clear();
+ n2nXSMatrix.clear();
+
+ return GeneratedDB;
+
+}
+
+
+void IM_RK4::ResetTheMatrix()
+{
+
+ if(fTheMatrix)
+ {
+ for(int i= 0; i<fNVar; i++)
+ delete [] fTheMatrix[i];
+ delete [] fTheMatrix;
+ }
+ fTheMatrix = 0;
+}
+
+void IM_RK4::SetTheMatrixToZero()
+{
+ ResetTheMatrix();
+
+ fNVar = findex.size();
+ fTheMatrix = new double*[fNVar];
+
+#pragma omp parallel for
+ for(int i= 0; i < fNVar; i++)
+ fTheMatrix[i] = new double[fNVar];
+
+ for(int i = 0; i < fNVar; i++)
+ for(int k = 0; k < fNVar; k++)
+ {
+ fTheMatrix[i][k]=0.0;
+ }
+
+}
+
+//________________________________________________________________________
+void IM_RK4::ResetTheNucleiVector()
+{
+ if(fTheNucleiVector)
+ delete [] fTheNucleiVector;
+ fTheNucleiVector = 0;
+}
+
+//________________________________________________________________________
+void IM_RK4::SetTheNucleiVectorToZero()
+{
+ ResetTheNucleiVector();
+ fTheNucleiVector = new double[fNVar];
+
+#pragma omp parallel for
+ for(int i = 0; i < fNVar; i++)
+ fTheNucleiVector[i]=0.0;
+
+}
+
+//________________________________________________________________________
+void IM_RK4::BuildEqns(double t, double *N, double *dNdt)
+{
+ double sum=0;
+ // pragma omp parallel for reduction(+:sum)
+ for(int i = 0; i < fNVar; i++)
+ {
+ sum=0;
+ for(int k = 0; k < fNVar; k++)
+ {
+ sum += fTheMatrix[i][k]*N[k];
+ }
+ dNdt[i] = sum;
+ }
+}
+
+//________________________________________________________________________
+void IM_RK4::SetTheMatrix(TMatrixT<double> BatemanMatrix)
+{
+ for (int k = 0; k < (int)fNVar; k++)
+ for (int l = 0; l < (int)findex_inver.size(); l++)
+ fTheMatrix[l][k] = BatemanMatrix[l][k];
+}
+
+//________________________________________________________________________
+TMatrixT<double> IM_RK4::GetTheMatrix()
+{
+ TMatrixT<double> BatemanMatrix = TMatrixT<double>(findex.size(),findex.size());
+ for (int k = 0; k < (int)fNVar; k++)
+ for (int l = 0; l < (int)findex_inver.size(); l++)
+ BatemanMatrix[l][k] = fTheMatrix[l][k];
+
+ return BatemanMatrix;
+}
+
+//________________________________________________________________________
+void IM_RK4::SetTheNucleiVector(TMatrixT<double> NEvolutionMatrix)
+{
+ for (int k = 0; k < (int)fNVar; k++)
+ fTheNucleiVector[k] = NEvolutionMatrix[k][0];
+}
+
+//________________________________________________________________________
+TMatrixT<double> IM_RK4::GetTheNucleiVector()
+{
+ TMatrixT<double> NEvolutionMatrix = TMatrixT<double>(findex.size(),1);
+ for (int k = 0; k < (int)fNVar; k++)
+ NEvolutionMatrix[k][0] = fTheNucleiVector[k];
+
+ return NEvolutionMatrix;
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/source/branches/CLASSV3/Model/Irradiation/IM_RK4.hxx b/source/branches/CLASSV3/Model/Irradiation/IM_RK4.hxx
new file mode 100644
index 0000000000000000000000000000000000000000..d829105d3c4e179b5bce5232e34b960252c914c4
--- /dev/null
+++ b/source/branches/CLASSV3/Model/Irradiation/IM_RK4.hxx
@@ -0,0 +1,102 @@
+#ifndef _IMRK4_HXX
+#define _IMRK4_HXX
+
+
+/*!
+ \file
+ \brief Header file for IrradiationModel class.
+
+
+ @author BaM
+ @version 2.0
+ */
+#include "DynamicalSystem.hxx"
+#include "IrradiationModel.hxx"
+
+
+using namespace std;
+
+
+class LogFile;
+
+//-----------------------------------------------------------------------------//
+/*!
+ Define a IM_RK4.
+ The aim of these class is perform the calculation of the evolution of a fuel trough irradiation solving numericaly the Bateman using RK4.
+
+
+ @author BaM
+ @version 3.0
+ */
+//________________________________________________________________________
+
+class EvolutionData;
+
+class IM_RK4 : public DynamicalSystem, IrradiationModel
+{
+
+ public :
+
+ IM_RK4();
+ IM_RK4(LogFile* Log);
+
+
+
+ /// virtueal 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
+ */
+ EvolutionData GenerateEvolutionData(IsotopicVector IV, EvolutionData XSSet, double Power, double cycletime);
+ //}
+
+ //********* RK4 Method *********//
+
+ //@}
+ /*!
+ \name RK4 Method
+ */
+ //@{
+
+
+ using DynamicalSystem::RungeKutta;
+ //! Pre-treatment Runge-Kutta method.
+ /*!
+ // This method does initialisation and then call DynamicalSystem::RungeKutta
+ // \param t1: initial time
+ // \param t2: final time
+ */
+
+
+ void BuildEqns(double t, double *N, double *dNdt);
+ void SetTheMatrixToZero(); //!< Initialize the evolution Matrix
+ void ResetTheMatrix();
+ void SetTheMatrix(TMatrixT<double> BatemanMatrix); //!< Set the Evolution Matrix (Bateman equations)
+ TMatrixT<double> GetTheMatrix(); //!< return the Evolution Matrix (Bateman equations)
+
+ void SetTheNucleiVectorToZero(); //!< Initialize the evolution Matrix
+ void ResetTheNucleiVector();
+ void SetTheNucleiVector(TMatrixT<double> NEvolutionMatrix); //!< Set the Evolution Matrix (Bateman equations)
+ TMatrixT<double> GetTheNucleiVector(); //!< return the Evolution Matrix (Bateman equations)
+ //@}
+
+
+
+ private :
+
+ double *fTheNucleiVector; //!< The evolving atoms copied from Material proportions.
+ double **fTheMatrix; //!< The evolution Matrix
+
+ double fPrecision; //!< Precision of the RungeKutta
+ double fHestimate; //!< RK Step estimation.
+ double fHmin; //!< RK minimum Step.
+ double fMaxHdid; //!< store the effective RK max step
+ double fMinHdid; //!< store the effective RK min step
+ bool fIsNegativeValueAllowed; //!< whether or not negative value are physical.
+
+
+};
+
+#endif
+
diff --git a/source/branches/CLASSV3/Model/XS/XSM_CLOSEST.cxx b/source/branches/CLASSV3/Model/XS/XSM_CLOSEST.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..756c448cae9881e9b2d25e70648e0a9db9b8ae4b
--- /dev/null
+++ b/source/branches/CLASSV3/Model/XS/XSM_CLOSEST.cxx
@@ -0,0 +1,366 @@
+#include "XSModel.hxx"
+#include "XSM_CLOSEST.hxx"
+#include "LogFile.hxx"
+#include "StringLine.hxx"
+
+
+#include <TGraph.h>
+#include <TString.h>
+#include "TSystem.h"
+#include "TROOT.h"
+
+#include <sstream>
+#include <string>
+#include <iostream>
+#include <fstream>
+#include <algorithm>
+#include <cmath>
+
+
+
+
+//________________________________________________________________________
+//
+// XSM_CLOSEST
+//
+//
+//
+//
+//________________________________________________________________________
+XSM_CLOSEST::XSM_CLOSEST(LogFile* Log,string DB_index_file, bool oldreadmethod )
+{
+ SetLog(Log);
+ fOldReadMethod = oldreadmethod;
+ fDataBaseIndex = DB_index_file;
+ fDistanceType = 0;
+ fWeightedDistance = false;
+ fEvolutionDataInterpolation = false;
+ ReadDataBase();
+ if(IsLog())
+ {
+ // Warning
+ cout << "!!INFO!! !!!XSM_CLOSEST!!! A EvolutionData has been define :" << endl;
+ cout << "\t His index is : \"" << DB_index_file << "\"" << endl;
+ cout << "\t " << fFuelDataBank.size() << " EvolutionData have been read."<< endl << endl;
+
+ GetLog()->fLog << "!!INFO!! !!!XSM_CLOSEST!!! A EvolutionData has been define :" << endl;
+ GetLog()->fLog << "\t His index is : \"" << DB_index_file << "\"" << endl;
+ GetLog()->fLog << "\t " << fFuelDataBank.size() << " EvolutionData have been read."<< endl << endl;
+ }
+
+}
+//________________________________________________________________________
+XSM_CLOSEST::~XSM_CLOSEST()
+{
+ map<IsotopicVector ,EvolutionData >::iterator it_del;
+ for( it_del = fFuelDataBank.begin(); it_del != fFuelDataBank.end(); it_del++)
+ (*it_del).second.DeleteEvolutionData();
+ fFuelDataBank.clear();
+
+
+ for( it_del = fFuelDataBankCalculated.begin(); it_del != fFuelDataBankCalculated.end(); it_del++)
+ (*it_del).second.DeleteEvolutionData();
+ fFuelDataBankCalculated.clear();
+}
+//________________________________________________________________________
+void XSM_CLOSEST::ReadDataBase()
+{
+
+ if(fFuelDataBank.size() != 0)
+ {
+ map<IsotopicVector ,EvolutionData >::iterator it_del;
+ for( it_del = fFuelDataBank.begin(); it_del != fFuelDataBank.end(); it_del++)
+ (*it_del).second.DeleteEvolutionData();
+ fFuelDataBank.clear();
+ }
+
+ if(fFuelDataBankCalculated.size() != 0)
+ {
+ map<IsotopicVector ,EvolutionData >::iterator it_del;
+ for( it_del = fFuelDataBankCalculated.begin(); it_del != fFuelDataBankCalculated.end(); it_del++)
+ (*it_del).second.DeleteEvolutionData();
+ fFuelDataBankCalculated.clear();
+ }
+
+
+ ifstream DataDB(fDataBaseIndex.c_str()); // Open the File
+ if(!DataDB)
+ {
+ cout << "!!Warning!! !!!FuelDataBank!!! \n Can't open \"" << fDataBaseIndex << "\"\n" << endl;
+ GetLog()->fLog << "!!Warning!! !!!FuelDataBank!!! \n Can't open \"" << fDataBaseIndex << "\"\n" << endl;
+ }
+ vector<double> vTime;
+ vector<double> vTimeErr;
+
+ string line;
+ int start = 0;
+
+
+ // First Get Fuel Type
+ getline(DataDB, line);
+ if( StringLine::NextWord(line, start, ' ') != "TYPE")
+ {
+ cout << "!!Bad Trouble!! !!!FuelDataBank!!! Bad Database file : " << fDataBaseIndex << " Can't find the type of the DataBase"<< endl;
+ GetLog()->fLog << "!!Bad Trouble!! !!!FuelDataBank!!! Bad Database file : " << fDataBaseIndex << " Can't find the type of the DataBase"<< endl;
+ exit (1);
+ }
+ fFuelType = StringLine::NextWord(line, start, ' ');
+
+ //Then Get All the Database
+
+ while (!DataDB.eof())
+ {
+ getline(DataDB, line);
+ if(line != "")
+ {
+ EvolutionData* evolutionproduct = new EvolutionData(GetLog(), line, fOldReadMethod);
+ IsotopicVector ivtmp = evolutionproduct->GetIsotopicVectorAt(0.).GetActinidesComposition();
+ fFuelDataBank.insert( pair<IsotopicVector, EvolutionData >(ivtmp , (*evolutionproduct) ));
+ }
+ }
+
+}
+//________________________________________________________________________
+//________________________________________________________________________
+//________________________________________________________________________
+/* Distance Calculation */
+//________________________________________________________________________
+//________________________________________________________________________
+//________________________________________________________________________
+map<double, EvolutionData> XSM_CLOSEST::GetDistancesTo(IsotopicVector isotopicvector, double t) const
+{
+
+ map<double, EvolutionData> distances;
+
+ map<IsotopicVector, EvolutionData > evolutiondb = fFuelDataBank;
+
+ map<IsotopicVector, EvolutionData >::iterator it;
+ for( it = evolutiondb.begin(); it != evolutiondb.end(); it++ )
+ {
+ pair<map<double, EvolutionData>::iterator, bool> IResult;
+
+ double D = Distance(isotopicvector.GetActinidesComposition(), (*it).second.GetIsotopicVectorAt(t).GetActinidesComposition()/ Norme( (*it).second.GetIsotopicVectorAt(t).GetActinidesComposition() )*Norme(isotopicvector.GetActinidesComposition())
+ ,fDistanceType, fDistanceParameter);
+
+ IResult = distances.insert( pair<double, EvolutionData>( D , (*it).second ) );
+ }
+
+ return distances;
+
+}
+//________________________________________________________________________
+EvolutionData XSM_CLOSEST::GetCrossSections(IsotopicVector isotopicvector, double t)
+{
+
+ map<IsotopicVector, EvolutionData > evolutiondb = fFuelDataBank;
+ double distance = 0;
+
+ map<IsotopicVector, EvolutionData >::iterator it_close = evolutiondb.begin();
+
+
+ map<IsotopicVector, EvolutionData >::iterator it;
+
+
+ if(fWeightedDistance)
+ {
+ distance = Distance(isotopicvector.GetActinidesComposition()
+ * evolutiondb.begin()->second.GetIsotopicVectorAt(t).GetActinidesComposition().GetSumOfAll()
+ / isotopicvector.GetActinidesComposition().GetSumOfAll(),
+ evolutiondb.begin()->second);
+
+
+ for( it = evolutiondb.begin(); it != evolutiondb.end(); it++ )
+ {
+ double D = 0;
+ D = Distance(isotopicvector.GetActinidesComposition()
+ * (*it).second.GetIsotopicVectorAt(t).GetActinidesComposition().GetSumOfAll()
+ / isotopicvector.GetActinidesComposition().GetSumOfAll(),
+ (*it).second);
+
+ if (D< distance)
+ {
+ distance = D;
+ it_close = it;
+ }
+ }
+
+ return (*it_close).second;
+ }
+ else if (fEvolutionDataInterpolation)
+ {
+
+ map<double, EvolutionData> distance_map = GetDistancesTo(isotopicvector, t);
+ map<double, EvolutionData>::iterator it_distance;
+ int NClose = 64;
+ int Nstep = 0;
+ EvolutionData EvolInterpolate;
+ double SumOfDistance = 0;
+ for( it_distance = distance_map.begin(); it_distance != distance_map.end(); it_distance++)
+ {
+
+ if((*it_distance).first == 0 )
+ {
+ EvolInterpolate = Multiply(1,(*it_distance).second);
+ return EvolInterpolate;
+ }
+ if(Nstep == 0)
+ EvolInterpolate = Multiply(1./(*it_distance).first, (*it_distance).second);
+ else
+ {
+
+ EvolutionData Evoltmp = EvolInterpolate;
+ EvolutionData Evoltmp2 = Multiply(1./(*it_distance).first, (*it_distance).second);
+
+ EvolInterpolate = Sum(Evoltmp, Evoltmp2);
+ Evoltmp.DeleteEvolutionData();
+ Evoltmp2.DeleteEvolutionData();
+
+
+ }
+
+ SumOfDistance += 1./(*it_distance).first;
+ Nstep++;
+ if(Nstep == NClose) break;
+
+ }
+
+ EvolutionData Evoltmp = EvolInterpolate;
+ EvolInterpolate.Clear();
+
+ EvolInterpolate = Multiply(1/SumOfDistance, Evoltmp);
+
+ Evoltmp.DeleteEvolutionData();
+ return EvolInterpolate;
+
+ }
+ else
+ {
+ distance = Distance(isotopicvector.GetActinidesComposition(),
+ evolutiondb.begin()->second.GetIsotopicVectorAt(t).GetActinidesComposition()
+ / evolutiondb.begin()->second.GetIsotopicVectorAt(t).GetActinidesComposition().GetSumOfAll()
+ * isotopicvector.GetActinidesComposition().GetSumOfAll(),
+ fDistanceType, fDistanceParameter);
+ for( it = evolutiondb.begin(); it != evolutiondb.end(); it++ )
+ {
+
+ double D = 0;
+
+
+ D = Distance(isotopicvector.GetActinidesComposition(),
+ (*it).second.GetIsotopicVectorAt(t).GetActinidesComposition()
+ / (*it).second.GetIsotopicVectorAt(t).GetActinidesComposition().GetSumOfAll()
+ * isotopicvector.GetActinidesComposition().GetSumOfAll(),
+ fDistanceType, fDistanceParameter);
+
+ if (D< distance)
+ {
+ distance = D;
+ it_close = it;
+ }
+ }
+ return (*it_close).second;
+
+ }
+
+}
+//________________________________________________________________________
+void XSM_CLOSEST::CalculateDistanceParameter()
+{
+
+ if(fDistanceType!=1){
+ cout << "!!Warning!! !!!CalculateDistanceParameter!!!"
+ << " Distance Parameter will be calculate even if the distance type is not the good one. Any Distance Parameters given by the user will be overwriten"<<endl;
+
+ GetLog()->fLog << "!!Warning!! !!!CalculateDistanceParameter!!!"
+ << " Distance Parameter will be calculate even if the distance type is not the good one. Any Distance Parameters given by the user will be overwriten"<<endl;
+ }
+
+ fDistanceParameter.Clear();
+
+ //We calculate the weight for the distance calculation.
+ map<IsotopicVector ,EvolutionData >::iterator it;
+ map<IsotopicVector ,EvolutionData > FuelDataBank = (*this).GetFuelDataBank();
+ int NevolutionDatainFuelDataBank = 0;
+
+ for( it = FuelDataBank.begin(); it != FuelDataBank.end(); it++ )
+ {
+ NevolutionDatainFuelDataBank++;
+ map<ZAI ,double>::iterator itit;
+ map<ZAI ,double> isovector=(*it).first.GetIsotopicQuantity();
+ for(itit=isovector.begin(); itit != isovector.end(); itit++) //Boucle sur ZAI
+ {
+ double TmpXS=0;
+
+ for( int i=1; i<4; i++ ) //Loop on Reactions 1==fission, 2==capture, 3==n2n
+ TmpXS+= (*it).second.GetXSForAt(0, (*itit).first, i);
+
+ fDistanceParameter.Add((*itit).first,TmpXS);
+ }
+
+
+ }
+ fDistanceParameter.Multiply( (double)1.0/NevolutionDatainFuelDataBank );
+
+ if(GetLog()){
+ GetLog()->fLog <<"!!INFO!! Distance Parameters "<<endl;
+ map<ZAI ,double >::iterator it2;
+ for(it2 = fDistanceParameter.GetIsotopicQuantity().begin();it2 != fDistanceParameter.GetIsotopicQuantity().end(); it2++)
+ {
+ GetLog()->fLog << (*it2).first.Z() << " ";
+ GetLog()->fLog << (*it2).first.A() << " ";
+ GetLog()->fLog << (*it2).first.I() << " ";
+ GetLog()->fLog << ": " << (*it2).second;
+ GetLog()->fLog << endl;
+ }
+ GetLog()->fLog << endl;
+ }
+
+
+}
+//________________________________________________________________________
+void XSM_CLOSEST::SetDistanceParameter(IsotopicVector DistanceParameter)
+{
+
+ fDistanceParameter = DistanceParameter;
+
+ GetLog()->fLog <<"!!INFO!! Distance Parameters "<<endl;
+ map<ZAI ,double >::iterator it2;
+ for(it2 = fDistanceParameter.GetIsotopicQuantity().begin();it2 != fDistanceParameter.GetIsotopicQuantity().end(); it2++)
+ {
+ GetLog()->fLog << (*it2).first.Z() << " ";
+ GetLog()->fLog << (*it2).first.A() << " ";
+ GetLog()->fLog << (*it2).first.I() << " ";
+ GetLog()->fLog << ": " << (*it2).second;
+ GetLog()->fLog << endl;
+ }
+ GetLog()->fLog << endl;
+
+}
+
+//________________________________________________________________________
+void XSM_CLOSEST::SetDistanceType(int DistanceType)
+{
+
+ fDistanceType=DistanceType;
+ if(fDistanceType==1){
+ CalculateDistanceParameter();
+ }
+ else if(fDistanceType == 2 && Norme(fDistanceParameter)==0){
+ // This is so bad!! You will probably unsynchronize all the reactor....
+ cout << "!!Warning!! !!!DistanceType!!!"
+ << " Distance use weight defined by user for each isotope, but no weight have been given" << endl<<"Use SetDistanceParameter()"<<endl;
+
+ GetLog()->fLog << "!!Warning!! !!!DistanceType!!!"
+ << " Distance use weight defined by user for each isotope, but no weight have been given" << endl<<"Use SetDistanceParameter()"<<endl;
+ exit(1);
+ }
+ else if (fDistanceType != 0 && fDistanceType != 1 && fDistanceType != 2 ){
+ cout << "!!ERROR!! !!!DistanceType!!!"
+ << " Distancetype defined by the user isn't recognized by the code"<<endl;
+
+ GetLog()->fLog << "!!ERROR!! !!!DistanceType!!!"
+ << " Distancetype defined by the user isn't recognized by the code"<<endl;
+ exit(1);
+ }
+
+}
\ No newline at end of file
diff --git a/source/branches/CLASSV3/Model/XS/XSM_CLOSEST.hxx b/source/branches/CLASSV3/Model/XS/XSM_CLOSEST.hxx
new file mode 100644
index 0000000000000000000000000000000000000000..2f861248d83532cc0d9bc212a5042093fcf0291e
--- /dev/null
+++ b/source/branches/CLASSV3/Model/XS/XSM_CLOSEST.hxx
@@ -0,0 +1,136 @@
+
+#ifndef _XSM_CLOSEST_HXX
+#define _XSM_CLOSEST_HXX
+
+
+/*!
+ \file
+ \brief Header file for XSM_MLP_PWR_MOX class.
+
+
+ @authors BaM,BLG
+ @version 1.0
+ */
+#include "XSModel.hxx"
+#include <string>
+#include <fstream>
+#include <iostream>
+#include <map>
+#include <vector>
+typedef long long int cSecond;
+using namespace std;
+
+//-----------------------------------------------------------------------------//
+/*!
+ Define a XSM_CLOSEST.
+CLASS to get cross sections from a set of pre-calculation
+(with MURE,or other depletion code)
+get cross sections of the closest (in composition) calculation
+
+ @authors BaM,BLG
+ @version 1.0
+ */
+//________________________________________________________________________
+
+
+class XSM_CLOSEST : public XSModel
+{
+
+public :
+
+ /*!
+ \name Constructor/Desctructor
+ */
+ //@{
+
+ XSM_CLOSEST(LogFile* Log, string DB_index_file, bool oldreadmethod = true );
+ ~XSM_CLOSEST();
+ //{
+
+
+ //********* Get Method *********//
+ /*!
+ \name Get Method
+ */
+ //@{
+ EvolutionData GetCrossSections(IsotopicVector isotopicvector,double t=0) ; //!< Reason to live of this CLASS Return the closest Evolutiondata
+ map<IsotopicVector ,EvolutionData > GetFuelDataBank() const { return fFuelDataBank; } //!< Return the FuelDataBank
+ string GetDataBaseIndex() const { return fDataBaseIndex; } //!< Return the index Name
+ string GetFuelType() const { return fFuelType; } //!< Return the fuel type of the DB
+ vector<double> GetFuelParameter() const { return fFuelParameter; } //!< Return the Fuel parameter of the DB
+ pair<double,double> GetBurnUpRange() const { return fBurnUpRange;} //!< Return the BurnUp range of the DB
+ bool IsDefine(IsotopicVector IV) const; //!< True the key is define, false unstead
+
+ map<double, EvolutionData> GetDistancesTo(IsotopicVector isotopicvector, double t = 0) const; //! Return a map containing the distance of each EvolutionData in the DataBase to the set IV at the t time
+ //@}
+
+ //********* Set Method *********//
+
+ /*!
+ \name Set Method
+ */
+ //@{
+
+ void SetFuelDataBank(map< IsotopicVector ,EvolutionData > mymap) { fFuelDataBank = mymap; } //!< Set the FuelDataBank map
+
+ void SetDataBaseIndex(string database) { fDataBaseIndex = database;; ReadDataBase(); } //!< Set the Name of the database index
+
+ void SetOldReadMethod(bool val) { fOldReadMethod = val; ReadDataBase();} ///< use the old reading method
+
+
+
+ void SetWeightedDistanceCalculation(bool val = true) { fWeightedDistance = val;} ///< Set weighted Distance calculation
+ void SetEvolutionDataInterpolation(bool val = true) { fEvolutionDataInterpolation = val;} ///< Set weighted Distance calculation
+ void SetDistanceParameter(IsotopicVector DistanceParameter); ///< Define mannually the weight for each ZAI in the distance calculation
+
+
+ //{
+ /// Define the way to decide if two isotopic vectors are close.
+ /*!
+ // The different algorythm are:
+ // \li 0 is for the standard norme,
+ // \li 1 for each ZAI weighted with its XS,
+ // \li 2 for each ZAI weighted with coefficient given by the user.
+ */
+ void SetDistanceType(int DistanceType);
+ //}
+
+ //********* Other Method *********//
+ /*!
+ \name Other Method
+ */
+ //@{
+ void ReadDataBase(); ///< read the index file and fill the evolutionData map
+ void CalculateDistanceParameter(); ///< Calcul of the weight for each ZAI in the distance calculation from the mean XS of the FuelDataBank
+
+ //@}
+
+ private :
+
+ map<IsotopicVector, EvolutionData> fFuelDataBank; ///< DataBanck map
+ map<IsotopicVector, EvolutionData> fFuelDataBankCalculated; ///< Map of the already calculated EvolutionData (to avoid recalculation...)
+
+ string fDataBaseIndex; ///< Name of the index
+
+ bool fOldReadMethod; ///< use old DB format
+ bool fWeightedDistance; ///< USe XS weighted distance calculation
+ bool fEvolutionDataInterpolation; ///< USe XS weighted distance calculation
+
+
+ string fFuelType; ///< Type of fuel of the FuelDataBank
+ pair<double,double> fBurnUpRange; ///< Range of the Burn-up range of the FuelDataBank
+ vector<double> fFuelParameter; ///< Parameter needed by the equivalence model
+
+
+
+ int fDistanceType; ///< Set the distance calculation algorytm
+ /// \li 0 is for the standard norm (Default = 0),
+ /// \li 1 for each ZAI weighted with its XS,
+ /// \li 2 for each ZAI weighted with coefficient given by the user.
+
+ IsotopicVector fDistanceParameter; ///< weight for each ZAI in the distance calculation
+
+};
+
+#endif
+
diff --git a/source/branches/CLASSV3/Model/XS/XSM_MLP_PWR_MOX.cxx b/source/branches/CLASSV3/Model/XS/XSM_MLP_PWR_MOX.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..a2e67b4fb54deba29360132fd7bdf293bc7171b2
--- /dev/null
+++ b/source/branches/CLASSV3/Model/XS/XSM_MLP_PWR_MOX.cxx
@@ -0,0 +1,604 @@
+
+#include "XSModel.hxx"
+#include "XSM_MLP_PWR_MOX.hxx"
+#include "LogFile.hxx"
+#include "StringLine.hxx"
+
+#include "TMVA/Reader.h"
+#include "TMVA/Tools.h"
+#include "TMVA/MethodCuts.h"
+
+#include <TGraph.h>
+#include <TString.h>
+#include "TFile.h"
+#include "TTree.h"
+#include "TSystem.h"
+#include "TROOT.h"
+#include "TStopwatch.h"
+
+#include <dirent.h>
+#include <errno.h>
+#include <sstream>
+#include <string>
+#include <iostream>
+#include <fstream>
+#include <algorithm>
+#include <cmath>
+
+//________________________________________________________________________
+//
+// XSM_MLP_PWR_MOX
+//
+//
+//
+//
+//________________________________________________________________________
+XSM_MLP_PWR_MOX::XSM_MLP_PWR_MOX(LogFile* Log,string TMVA_Weight_Directory,string InformationFile, bool IsTimeStep)
+{
+
+ SetLog(Log);
+
+ fIsStepTime=IsTimeStep;
+ fTMVAWeightFolder = TMVA_Weight_Directory;
+ if(InformationFile=="")
+ fMLPInformationFile=TMVA_Weight_Directory+"/Data_Base_Info.nfo";
+ else
+ fMLPInformationFile=InformationFile;
+
+ GetMLPWeightFiles();
+ SetDataBaseInformation();
+
+ if(IsLog())
+ {
+ // Warning
+ cout << "!!INFO!! !!!XSM_MLP_PWR_MOX!!! A EvolutionData has been define :" << endl;
+ cout << "\t His TMVA folder is : \"" << fTMVAWeightFolder << "\"" << endl;
+
+
+ GetLog()->fLog << "!!INFO!! !!!XSM_MLP_PWR_MOX!!! A EvolutionData has been define :" << endl;
+ GetLog()->fLog <<"\t His TMVA folder is : \"" << fTMVAWeightFolder << "\"" << endl;
+ }
+}
+
+//________________________________________________________________________
+void XSM_MLP_PWR_MOX::SetDataBaseInformation()
+{
+ ifstream FILE(fMLPInformationFile.c_str());
+
+ if(FILE.good())
+ {
+ double HM_Mass_tonne=0;
+ double Power_watt=0;
+
+ FILE>>HM_Mass_tonne;
+ FILE>>Power_watt;
+
+ fDataBaseHMMass=HM_Mass_tonne;
+ fDataBasePower=Power_watt;
+
+ while(FILE.eof())
+ {
+ double TIME=-1;
+ FILE>>TIME;
+
+ if(TIME!=-1)
+ fMLP_Time.push_back(TIME);
+ }
+ }
+ else
+ {
+ cout<<"Can't find/open file "<<fMLPInformationFile<<endl;
+ exit(0);
+ }
+
+}
+//________________________________________________________________________
+void XSM_MLP_PWR_MOX::GetMLPWeightFiles()
+{
+
+ /**********Get All TMVA weight files*******************/
+
+ //check if the folder containing weights exists
+ DIR* rep = NULL;
+ struct dirent* fichierLu = NULL;
+ rep = opendir(fTMVAWeightFolder.c_str());
+ if (rep == NULL)
+ {
+ perror("");
+ exit(1);
+ }
+
+ /**Save file names of TMVA weights*/
+ fWeightFiles.resize(0);
+ while ((fichierLu = readdir(rep)) != NULL)
+ {
+ string FileName= fichierLu->d_name ;
+ if(FileName != "." && FileName != ".." )
+ {
+ if(FileName[FileName.size()-3]=='x' && FileName[FileName.size()-2]=='m' && FileName[FileName.size()-1]=='l' )
+ fWeightFiles.push_back(FileName);
+
+ }
+ }
+
+}
+//________________________________________________________________________
+//________________________________________________________________________
+//
+// Time (MLP take time as parameter)
+//________________________________________________________________________
+//________________________________________________________________________
+void XSM_MLP_PWR_MOX::ReadWeightFile(string Filename, int &Z, int &A, int &I, int &Reaction)
+{
+ Z=-1;
+ A=-1;
+ I=-1;
+ Reaction=-1;
+
+ size_t found = Filename.find("XS");
+
+ string NameJOB;
+ NameJOB=Filename.substr(found);
+ int pos=0;
+
+ StringLine::NextWord(NameJOB, pos, '_');
+
+ Z = atof( (StringLine::NextWord(NameJOB,pos,'_') ).c_str() );
+
+ A = atof( (StringLine::NextWord(NameJOB,pos,'_') ).c_str() );
+
+ I = atof( (StringLine::NextWord(NameJOB,pos,'_') ).c_str() );
+
+ string sReaction = (StringLine::NextWord(NameJOB,pos,'_') ).c_str() ;
+ if(sReaction=="fis")
+ Reaction=0;
+ if(sReaction=="cap")
+ Reaction=1;
+ if(sReaction=="n2n")
+ Reaction=2;
+
+}
+//________________________________________________________________________
+void XSM_MLP_PWR_MOX::CreateTMVAInputTree(IsotopicVector isotopicvector,int TimeStep)
+{
+ /******Create Input data tree to be interpreted by TMVA::Reader***/
+
+ TFile* InputFile = new TFile("./InputTMP.root","RECREATE");
+ TTree* InputTree = new TTree("InTMP", "InTMP");
+ float teneur = 0;
+ float Pu8 = 0;
+ float Pu9 = 0;
+ float Pu10 = 0;
+ float Pu11 = 0;
+ float Pu12 = 0;
+ float Am1 = 0;
+ float U5_enrichment = 0;
+ float Time = 0;
+
+ InputTree->Branch( "teneur",&teneur,"teneur/F" );
+ 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( "Time" ,&Time ,"Time/F" );
+
+
+
+ float U8 = isotopicvector.GetZAIIsotopicQuantity(92,238,0);
+ float U5 = isotopicvector.GetZAIIsotopicQuantity(92,235,0);
+ float U4 = isotopicvector.GetZAIIsotopicQuantity(92,234,0);
+
+ float UTOT = U8 + U5 + U4;
+
+ Pu8 = isotopicvector.GetZAIIsotopicQuantity(94,238,0);
+ Pu9 = isotopicvector.GetZAIIsotopicQuantity(94,239,0);
+ Pu10 = isotopicvector.GetZAIIsotopicQuantity(94,240,0);
+ Pu11 = isotopicvector.GetZAIIsotopicQuantity(94,241,0);
+ Pu12 = isotopicvector.GetZAIIsotopicQuantity(94,242,0);
+ Am1 = isotopicvector.GetZAIIsotopicQuantity(95,241,0);
+
+ teneur = (Pu8+Pu9+Pu10+Pu11+Pu12+Am1)/(Pu8+Pu9+Pu10+Pu11+Pu12+Am1+U8+U5+U4); //prop mol. Pu
+
+ 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;
+
+ Time=fMLP_Time[TimeStep];
+
+ 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 ...
+ {
+ cout<<"!!!!!!!!!!!ERRORR!!!!!!!!!!!!"<<endl;
+ cout<<Pu8<<" "<<Pu9<<" "<<Pu10<<" "<<Pu11<<" "<<Pu12<<" "<<Am1<<endl;
+ exit(0);
+ }
+ // All value are molar (!weight)
+
+ InputTree->Fill();
+
+ InputFile->Write();
+ delete InputTree;
+ InputFile-> Close();
+ delete InputFile;
+}
+///________________________________________________________________________
+double XSM_MLP_PWR_MOX::ExecuteTMVA(string WeightFile)
+{
+ // --- Create the Reader object
+
+ TMVA::Reader *reader = new TMVA::Reader( "Silent" );
+
+ // Create a set of variables and declare them to the reader
+ // - the variable names MUST corresponds in name and type to those given in the weight file(s) used
+ Float_t Pu8,Pu9,Pu10,Pu11,Pu12,Am1,Time,teneur;
+ reader->AddVariable("teneur",&teneur);
+ 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( "Time" ,&Time);
+
+ // --- Book the MVA methods
+
+ string dir = fTMVAWeightFolder;
+ if(dir[dir.size()-1]!='/')
+ dir+="/";
+
+ // Book method MLP
+ TString methodName = "MLP method";
+ TString weightpath = dir + WeightFile ;
+ reader->BookMVA( methodName, weightpath );
+
+ // Prepare input tree
+ TFile *input(0);
+ if (!gSystem->AccessPathName( "./InputTMP.root" )) {
+ input = TFile::Open( "./InputTMP.root" ); // check if file in local directory exists
+ }
+ if (!input) {
+ std::cout << "ERROR: could not open data file" << std::endl;
+ exit(1);
+ }
+
+ TTree* theTree = (TTree*)input->Get("InTMP");
+
+ theTree->SetBranchAddress("teneur",&teneur);
+ 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( "Time" ,&Time );
+
+ theTree->GetEntry(0);
+ Float_t val = (reader->EvaluateRegression( methodName ))[0];
+
+ delete reader;
+ delete theTree;
+ input->Close();
+ delete input;
+ //cout<<"....done"<<endl;
+
+ //cout<<val<<endl;
+ return (double)val;
+}
+//________________________________________________________________________
+EvolutionData XSM_MLP_PWR_MOX::GetCrossSectionsTime(IsotopicVector IV)
+{
+
+ //cout<<"=====Building Evolution Data From TMVA MLP====="<<endl;
+
+ EvolutionData EvolutionDataFromMLP = EvolutionData();
+
+ map<ZAI,TGraph*> ExtrapolatedXS[3];
+ /*************DATA BASE INFO****************/
+ EvolutionDataFromMLP.SetReactorType("PWR");
+ EvolutionDataFromMLP.SetFuelType("MOX");
+ EvolutionDataFromMLP.SetPower(fDataBasePower);
+ EvolutionDataFromMLP.SetHeavyMetalMass(fDataBaseHMMass);
+ /************* The Cross sections***********/
+ for(int i=0;i<int(fWeightFiles.size());i++)
+ {
+ int Z=-2;
+ int A=-2;
+ int I=-2;
+ int Reaction=-2;
+ ReadWeightFile( fWeightFiles[i], Z, A, I, Reaction);
+
+ ZAI zaitmp = ZAI(Z,A,I);
+
+ for(int TimeStep=0;TimeStep<int(fMLP_Time.size());TimeStep++)
+ {
+ std::ifstream ifs ("InputTMP.root");
+ if (ifs.is_open())
+ { ifs.close();
+ system( "rm InputTMP.root" );
+ }
+ CreateTMVAInputTree(IV,TimeStep);
+
+ pair< map<ZAI, TGraph*>::iterator, bool> IResult;
+
+ IResult = ExtrapolatedXS[Reaction].insert(pair<ZAI ,TGraph* >(ZAI(Z,A,I), new TGraph() ) );
+
+
+ if(IResult.second )
+ {
+ (IResult.first)->second->SetPoint(0, (double)GetMLPTime()[TimeStep], ExecuteTMVA(fWeightFiles[i]) );
+ }
+ else
+ {
+ (IResult.first)->second->SetPoint( (IResult.first)->second->GetN(), (double)GetMLPTime()[TimeStep], ExecuteTMVA(fWeightFiles[i]) );
+ }
+
+ system( "rm InputTMP.root" );
+ }
+
+ }
+
+ /**********Sorting TGraph*********/
+ for(int x=0;x<3;x++)
+ { map<ZAI,TGraph*>::iterator it;
+ for(it = ExtrapolatedXS[x].begin(); it != ExtrapolatedXS[x].end(); it++)
+ it->second->Sort();
+ }
+ /**********Filling Matrices*/
+ EvolutionDataFromMLP.SetFissionXS(ExtrapolatedXS[0]);
+ EvolutionDataFromMLP.SetCaptureXS(ExtrapolatedXS[1]);
+ EvolutionDataFromMLP.Setn2nXS(ExtrapolatedXS[2]);
+
+ //cout<<"=====Evolution Data Built====="<<endl;
+ return EvolutionDataFromMLP;
+}
+//________________________________________________________________________
+//________________________________________________________________________
+//
+// Time step (1 MLP per time step)
+//________________________________________________________________________
+//________________________________________________________________________
+void XSM_MLP_PWR_MOX::ReadWeightFileStep(string Filename, int &Z, int &A, int &I, int &Reaction, int &TimeStep)
+{
+ Z=-1;
+ A=-1;
+ I=-1;
+ Reaction=-1;
+
+ size_t found = Filename.find("XS");
+
+ string NameJOB;
+ NameJOB=Filename.substr(found);
+ int pos=0;
+
+ StringLine::NextWord(NameJOB, pos, '_');
+
+ Z = atof( (StringLine::NextWord(NameJOB,pos,'_') ).c_str() );
+
+ A = atof( (StringLine::NextWord(NameJOB,pos,'_') ).c_str() );
+
+ I = atof( (StringLine::NextWord(NameJOB,pos,'_') ).c_str() );
+
+ string sReaction = (StringLine::NextWord(NameJOB,pos,'_') ).c_str() ;
+ if(sReaction=="fis")
+ Reaction=0;
+ if(sReaction=="cap")
+ Reaction=1;
+ if(sReaction=="n2n")
+ Reaction=2;
+
+ TimeStep = atof( (StringLine::NextWord(NameJOB,pos,'_') ).c_str() );
+
+}
+//________________________________________________________________________
+void XSM_MLP_PWR_MOX::CreateTMVAInputTreeStep(IsotopicVector isotopicvector)
+{
+ /******Create Input data tree to be interpreted by TMVA::Reader***/
+
+ TFile* InputFile = new TFile("./InputTMP.root","RECREATE");
+ TTree* InputTree = new TTree("InTMP", "InTMP");
+ float teneur = 0;
+ float Pu8 = 0;
+ float Pu9 = 0;
+ float Pu10 = 0;
+ float Pu11 = 0;
+ float Pu12 = 0;
+ float Am1 = 0;
+ float U5_enrichment = 0;
+
+ InputTree->Branch( "teneur",&teneur,"teneur/F" );
+ 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" );
+
+
+ float U8 = isotopicvector.GetZAIIsotopicQuantity(92,238,0);
+ float U5 = isotopicvector.GetZAIIsotopicQuantity(92,235,0);
+ float U4 = isotopicvector.GetZAIIsotopicQuantity(92,234,0);
+
+ float UTOT = U8 + U5 + U4;
+
+ Pu8 = isotopicvector.GetZAIIsotopicQuantity(94,238,0);
+ Pu9 = isotopicvector.GetZAIIsotopicQuantity(94,239,0);
+ Pu10 = isotopicvector.GetZAIIsotopicQuantity(94,240,0);
+ Pu11 = isotopicvector.GetZAIIsotopicQuantity(94,241,0);
+ Pu12 = isotopicvector.GetZAIIsotopicQuantity(94,242,0);
+ Am1 = isotopicvector.GetZAIIsotopicQuantity(95,241,0);
+
+ teneur = (Pu8+Pu9+Pu10+Pu11+Pu12+Am1)/(Pu8+Pu9+Pu10+Pu11+Pu12+Am1+U8+U5+U4); //prop mol. Pu
+
+ 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;
+
+ 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 ...
+ {
+ cout<<"!!!!!!!!!!!ERRORR!!!!!!!!!!!!"<<endl;
+ cout<<Pu8<<" "<<Pu9<<" "<<Pu10<<" "<<Pu11<<" "<<Pu12<<" "<<Am1<<endl;
+ exit(0);
+ }
+ // All value are molar (!weight)
+
+ InputTree->Fill();
+
+ InputFile->Write();
+ delete InputTree;
+ InputFile-> Close();
+ delete InputFile;
+}
+///________________________________________________________________________
+double XSM_MLP_PWR_MOX::ExecuteTMVAStep(string WeightFile)
+{
+ // --- Create the Reader object
+
+ TMVA::Reader *reader = new TMVA::Reader( "Silent" );
+
+ // Create a set of variables and declare them to the reader
+ // - the variable names MUST corresponds in name and type to those given in the weight file(s) used
+ Float_t Pu8,Pu9,Pu10,Pu11,Pu12,Am1,teneur;
+ reader->AddVariable("teneur",&teneur);
+ 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
+
+ string dir = fTMVAWeightFolder;
+ if(dir[dir.size()-1]!='/')
+ dir+="/";
+
+ // Book method MLP
+ TString methodName = "MLP method";
+ TString weightpath = dir + WeightFile ;
+ reader->BookMVA( methodName, weightpath );
+
+ // Prepare input tree
+ TFile *input(0);
+ if (!gSystem->AccessPathName( "./InputTMP.root" )) {
+ input = TFile::Open( "./InputTMP.root" ); // check if file in local directory exists
+ }
+ if (!input) {
+ std::cout << "ERROR: could not open data file" << std::endl;
+ exit(1);
+ }
+
+ TTree* theTree = (TTree*)input->Get("InTMP");
+
+ theTree->SetBranchAddress("teneur",&teneur);
+ 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;
+ input->Close();
+ delete input;
+ //cout<<"....done"<<endl;
+
+ //cout<<val<<endl;
+ return (double)val;
+}
+//________________________________________________________________________
+EvolutionData XSM_MLP_PWR_MOX::GetCrossSectionsStep(IsotopicVector IV)
+{
+
+ std::ifstream ifs ("InputTMP.root");
+ if (ifs.is_open())
+ { ifs.close();
+ system( "rm InputTMP.root" );
+ }
+ CreateTMVAInputTreeStep(IV);
+ //cout<<"=====Building Evolution Data From TMVA MLP====="<<endl;
+
+ EvolutionData EvolutionDataFromMLP = EvolutionData();
+
+ map<ZAI,TGraph*> ExtrapolatedXS[3];
+ /*************DATA BASE INFO****************/
+ EvolutionDataFromMLP.SetReactorType("PWR");
+ EvolutionDataFromMLP.SetFuelType("MOX");
+ EvolutionDataFromMLP.SetPower(fDataBasePower);
+ EvolutionDataFromMLP.SetHeavyMetalMass(fDataBaseHMMass);
+
+ /************* The Cross sections***********/
+
+ for(int i=0;i<int(fWeightFiles.size());i++)
+ {
+ int Z=-2;
+ int A=-2;
+ int I=-2;
+ int Reaction=-2;
+ int TimeStep=-2;
+ ReadWeightFileStep( fWeightFiles[i], Z, A, I, Reaction, TimeStep);
+
+ ZAI zaitmp = ZAI(Z,A,I);
+
+ pair< map<ZAI, TGraph*>::iterator, bool> IResult;
+
+ IResult = ExtrapolatedXS[Reaction].insert(pair<ZAI ,TGraph* >(ZAI(Z,A,I), new TGraph() ) );
+
+ if( IResult.second )
+ {
+ (IResult.first)->second->SetPoint(0, (double)GetMLPTime()[TimeStep], ExecuteTMVAStep(fWeightFiles[i]) );
+ }
+ else
+ {
+ (IResult.first)->second->SetPoint( (IResult.first)->second->GetN(), (double)GetMLPTime()[TimeStep], ExecuteTMVAStep(fWeightFiles[i]) );
+ }
+
+ }
+
+ system( "rm InputTMP.root" );
+ /**********Sorting TGraph*********/
+ for(int x=0;x<3;x++)
+ { map<ZAI,TGraph*>::iterator it;
+ for(it = ExtrapolatedXS[x].begin(); it != ExtrapolatedXS[x].end(); it++)
+ it->second->Sort();
+ }
+ /**********Filling Matrices*/
+ EvolutionDataFromMLP.SetFissionXS(ExtrapolatedXS[0]);
+ EvolutionDataFromMLP.SetCaptureXS(ExtrapolatedXS[1]);
+ EvolutionDataFromMLP.Setn2nXS(ExtrapolatedXS[2]);
+
+ //cout<<"=====Evolution Data Built====="<<endl;
+ return EvolutionDataFromMLP;
+}
+//________________________________________________________________________
+EvolutionData XSM_MLP_PWR_MOX::GetCrossSections(IsotopicVector IV)
+{
+ EvolutionData EV;
+ if(fIsStepTime)
+ EV=GetCrossSectionsStep(IV);
+
+ else
+ EV=GetCrossSectionsTime(IV);
+
+return EV;
+}
+//________________________________________________________________________
diff --git a/source/branches/CLASSV3/Model/XS/XSM_MLP_PWR_MOX.hxx b/source/branches/CLASSV3/Model/XS/XSM_MLP_PWR_MOX.hxx
new file mode 100644
index 0000000000000000000000000000000000000000..fe5495bd159a47fb40d11b609a0410ed43c80cec
--- /dev/null
+++ b/source/branches/CLASSV3/Model/XS/XSM_MLP_PWR_MOX.hxx
@@ -0,0 +1,80 @@
+
+#ifndef _XSM_MLP_PWR_MOX_HXX
+#define _XSM_MLP_PWR_MOX_HXX
+
+
+/*!
+ \file
+ \brief Header file for XSM_MLP_PWR_MOX class.
+
+
+ @authors BLG
+ @version 1.0
+ */
+#include "XSModel.hxx"
+#include <string>
+#include <fstream>
+#include <iostream>
+#include <map>
+#include <vector>
+typedef long long int cSecond;
+using namespace std;
+
+//-----------------------------------------------------------------------------//
+/*!
+ Define a XSM_MLP_PWR_MOX.
+This is the class to predict cross sections with a set of MultiLayerPerceptrons
+Design for a PWR MOX reactor
+
+ @authors BLG
+ @version 1.0
+ */
+//________________________________________________________________________
+
+
+class XSM_MLP_PWR_MOX : public XSModel
+{
+public :
+
+ /*!
+ \name Constructor/Desctructor
+ */
+ //@{
+
+ XSM_MLP_PWR_MOX(LogFile* Log,string TMVA_Weight_Directory,string InformationFile="",bool IsTimeStep=true);
+ ~XSM_MLP_PWR_MOX() {;}
+ //{
+
+
+ EvolutionData GetCrossSections(IsotopicVector IV) ;
+
+
+ private :
+ void SetDataBaseInformation(); //<! Read information file and fill HM mass, Power, time vector
+ vector<double> GetMLPTime() {return fMLP_Time; }
+
+ void GetMLPWeightFiles();
+ void ReadWeightFile(string Filename, int &Z, int &A, int &I, int &Reaction) ;
+ double ExecuteTMVA(string WeightFile);
+ void CreateTMVAInputTree(IsotopicVector isotopicvector,int TimeStep);
+ EvolutionData GetCrossSectionsTime(IsotopicVector IV) ;
+
+ void ReadWeightFileStep(string Filename, int &Z, int &A, int &I, int &Reaction, int &TimeStep) ;
+ double ExecuteTMVAStep(string WeightFile);
+ void CreateTMVAInputTreeStep(IsotopicVector isotopicvector);
+ EvolutionData GetCrossSectionsStep(IsotopicVector IV) ;
+
+
+
+ vector<double> fMLP_Time;
+ vector<string> fWeightFiles;
+ string fTMVAWeightFolder;
+ string fMLPInformationFile;
+ double fDataBasePower;
+ double fDataBaseHMMass;
+ bool fIsStepTime;
+
+};
+
+#endif
+
diff --git a/source/branches/CLASSV3/Model/XS/XSModel.cxx b/source/branches/CLASSV3/Model/XS/XSModel.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..606d7a125e2fd2134c09c6d0c80141ce6f62713a
--- /dev/null
+++ b/source/branches/CLASSV3/Model/XS/XSModel.cxx
@@ -0,0 +1,18 @@
+//
+// XSModel.cxx
+// CLASSSource
+//
+// Created by BaM on 04/05/2014.
+// Copyright (c) 2014 BLG. All rights reserved.
+//
+
+#include "XSModel.hxx"
+#include "CLASSHeaders.hxx"
+
+
+using namespace std;
+
+XSModel::XSModel(): CLASSObject()
+{
+
+}
diff --git a/source/branches/CLASSV3/Model/XS/XSModel.hxx b/source/branches/CLASSV3/Model/XS/XSModel.hxx
new file mode 100644
index 0000000000000000000000000000000000000000..e5f425ee5b5946f308cdd9e66c94789c6380b3ab
--- /dev/null
+++ b/source/branches/CLASSV3/Model/XS/XSModel.hxx
@@ -0,0 +1,52 @@
+
+#ifndef _XSMODEL_HXX
+#define _XSMODEL_HXX
+
+
+/*!
+ \file
+ \brief Header file for XSMODEL class.
+
+
+ @author BLG
+ @version 1.0
+ */
+#include "CLASSHeaders.hxx"
+
+
+using namespace std;
+
+//-----------------------------------------------------------------------------//
+/*!
+ Define a XS Interpolator
+This is the class for method related to XS prediction
+
+see @/Models/XS/include/MLPXSModel.hxx
+see @/Models/XS/include/ClosestXSModel.hxx
+...
+
+!!!!!!!!!!!!!!!!CAUTION!!!!!!!!!!!!!
+Never instantiate XSModel in your CLASS input but it's derivated class
+
+
+ @author BLG
+ @version 1.0
+ */
+//________________________________________________________________________
+
+
+class XSModel : public CLASSObject
+{
+
+ public :
+
+ XSModel();
+
+ virtual EvolutionData GetCrossSections(IsotopicVector IV,double t=0) {return 0;}
+
+
+
+};
+
+#endif
+