diff --git a/source/branches/CLASSV3/include/EQM_MLP_PWR_MOX.hxx b/source/branches/CLASSV3/include/EQM_MLP_PWR_MOX.hxx
index cc7cfeb05c6f7a531948e11d444724b4c90bdd3f..aa5eb0e50236e9b3d95138a41e4dad7717705f97 100644
--- a/source/branches/CLASSV3/include/EQM_MLP_PWR_MOX.hxx
+++ b/source/branches/CLASSV3/include/EQM_MLP_PWR_MOX.hxx
@@ -1,7 +1,7 @@
#ifndef _EQM_MLP_MOX_HXX
#define _EQM_MLP_MOX_HXX
-#include "CLASSHeaders.hxx"
+#include "EquivalenceModel.hxx"
/*!
\file
@@ -26,10 +26,8 @@ using namespace std;
*/
//________________________________________________________________________
-class IsotopicVector;
-
-class EQM_MLP_MOX : public EquivalenceModel
+class EQM_MLP_MOX : public EquivalenceModel
{
public :
@@ -47,11 +45,11 @@ class EQM_MLP_MOX : public EquivalenceModel
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();
- void GuessLambda(double& lambda, int& StockID,int FirstStockID, int LastStockID, double DeltaM,double StockHM);
string fTMVAWeightPath;
diff --git a/source/branches/CLASSV3/include/EQM_QUAD_PWR_MOX.hxx b/source/branches/CLASSV3/include/EQM_QUAD_PWR_MOX.hxx
new file mode 100644
index 0000000000000000000000000000000000000000..9a13fca25583a74334fd5a8da93ac21533520864
--- /dev/null
+++ b/source/branches/CLASSV3/include/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/include/EquivalenceModel.hxx b/source/branches/CLASSV3/include/EquivalenceModel.hxx
index 74044a80b8127fe8198e0298f1c8e0bce3c7af23..07a6be8d853e57efbe3c7525099638eb4cac0d9a 100644
--- a/source/branches/CLASSV3/include/EquivalenceModel.hxx
+++ b/source/branches/CLASSV3/include/EquivalenceModel.hxx
@@ -12,6 +12,7 @@
*/
#include "IsotopicVector.hxx"
+#include "CLASSObject.hxx"
using namespace std;
@@ -28,11 +29,15 @@ using namespace std;
//________________________________________________________________________
-class EquivalenceModel : public TObject
+class EquivalenceModel : public CLASSObject
{
public :
- EquivalenceModel() {}
+ EquivalenceModel();
+ EquivalenceModel(LogFile* log);
+
+ virtual ~EquivalenceModel();
+
/// virtueal 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....
@@ -43,9 +48,10 @@ class EquivalenceModel : public TObject
\param vector<IsotopicVector> FertilArray, isotopicvectors to use to get the fertil part of the fuel (if empty take it from the god)
*/
- virtual vector<double> BuildFuel(double BurnUp, double HMMass, vector<IsotopicVector> FissilArray, vector<IsotopicVector> FertilArray )
- { vector<double> empty; return empty;}
+ virtual vector<double> BuildFuel(double BurnUp, double HMMass, vector<IsotopicVector> FissilArray, vector<IsotopicVector> FertilArray );
//}
+ virtual void GuessLambda(double& lambda, int& StockID,int FirstStockID, int LastStockID, double DeltaM,double StockHM);
+ virtual double GetFissileMolarFraction(IsotopicVector Fissil,IsotopicVector Fertil,double BurnUp) {return 0;}
@@ -58,11 +64,13 @@ class EquivalenceModel : public TObject
protected :
+ IsotopicVector fFertileList; //!< contain the list of zai, needed as fertile, taken in a stock before fabrication
+ //!< if no stock are provided, will take the isotopic vector in the Park income
+
+ IsotopicVector fFissileList; //!< contain the list of zai, needed as fissile, taken in a stock before fabrication
+ //!< if no stock are provided, will take the isotopic vector in the Park income
- IsotopicVector fFertileList;
- IsotopicVector fFissileList;
-
};
#endif
diff --git a/source/branches/CLASSV3/src/EQM_MLP_PWR_MOX.cxx b/source/branches/CLASSV3/src/EQM_MLP_PWR_MOX.cxx
index 6357e560f6ad25911ec3977c90683cfde4313116..a2c57dde9d9b1eba3500a2a403a89e871da9887e 100644
--- a/source/branches/CLASSV3/src/EQM_MLP_PWR_MOX.cxx
+++ b/source/branches/CLASSV3/src/EQM_MLP_PWR_MOX.cxx
@@ -1,4 +1,3 @@
-#include "EquivalenceModel.hxx"
#include "EQM_MLP_PWR_MOX.hxx"
#include "LogFile.hxx"
#include "StringLine.hxx"
@@ -14,8 +13,9 @@
#include "TMVA/Reader.h"
#include "TMVA/Tools.h"
#include "TMVA/MethodCuts.h"
+#include "TFile.h"
-#include "CLASSHeaders.hxx"
+#include "IsotopicVector.hxx"
//________________________________________________________________________
//
@@ -47,6 +47,7 @@ EQM_MLP_MOX::EQM_MLP_MOX(string TMVAWeightPath)
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)
{
@@ -191,7 +192,7 @@ void EQM_MLP_MOX::GuessLambda(double& lambda, int& StockID,int FirstStockID, int
if(lambda >= 1) //this stock is not enought go to next one
{
- lambda=1;
+ lambda = 1;
if( !( StockID+1 == LastStockID) ) //if its the last stock don't try to look the next one (segfault otherwize)
StockID++;
@@ -219,6 +220,7 @@ void EQM_MLP_MOX::GuessLambda(double& lambda, int& StockID,int FirstStockID, int
}
}
+
//________________________________________________________________________
vector<double> EQM_MLP_MOX::BuildFuel(double BurnUp, double HMMass,vector<IsotopicVector> FissilArray, vector<IsotopicVector> FertilArray)
{
diff --git a/source/branches/CLASSV3/src/EQM_QUAD_PWR_MOX.cxx b/source/branches/CLASSV3/src/EQM_QUAD_PWR_MOX.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..aba54ecba9607ee6b90cc9ed7de1fd6476606ded
--- /dev/null
+++ b/source/branches/CLASSV3/src/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/src/EquivalenceModel.cxx b/source/branches/CLASSV3/src/EquivalenceModel.cxx
index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..fd698a155342695b04cac5384d5a7576e526fc36 100644
--- a/source/branches/CLASSV3/src/EquivalenceModel.cxx
+++ b/source/branches/CLASSV3/src/EquivalenceModel.cxx
@@ -0,0 +1,164 @@
+#include "EquivalenceModel.hxx"
+
+
+
+
+
+
+EquivalenceModel::EquivalenceModel():CLASSObject()
+{
+
+}
+
+EquivalenceModel::EquivalenceModel(LogFile* log):CLASSObject(log)
+{
+
+}
+
+EquivalenceModel::~EquivalenceModel()
+{
+
+}
+
+
+
+
+//________________________________________________________________________
+vector<double> EquivalenceModel::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)
+ Fertile.Add(U5, U5_enrich *HMMass/MeanMolarDepletedU*AVOGADRO);
+ Fertile.Add(U8, (1-U5_enrich)*HMMass/MeanMolarDepletedU*AVOGADRO);
+ */
+ IsotopicVector Fertile;
+ IsotopicVector Fissile;
+ 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
+ Fissile.Clear();
+ for(int i=0;i<=StockID;i++)
+ Fissile += lambda[i] * FissilArray[i];
+
+ AvailablePuMass = Fissile.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;
+ Fertile.Clear();
+ for(int i=int(FissilArray.size());i<=FertileStockID;i++)
+ { Fertile += lambda[i] * FertilArray[j];
+ j++;
+ }
+ AvailableFertilMass=Fertile.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]=0;
+ break;
+ }
+
+ }
+ /*Calcul the quantity of this composition needed to reach the burnup*/
+ double MolarPuContent = GetFissileMolarFraction(Fissile, Fertile, BurnUp);
+
+ double MeanMolarPu = Fissile.MeanMolar();
+ double MeanMolarDepletedU = Fertile.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]=0;
+ break;
+ }
+
+ }
+
+ return lambda;
+}
+
+
+
+//________________________________________________________________________
+void EquivalenceModel::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);
+ }
+
+ }
+
+ }
+
+}
\ No newline at end of file
diff --git a/source/branches/CLASSV3/src/Makefile b/source/branches/CLASSV3/src/Makefile
index 3108e59e4ce1c8c1d18dd5cd083f9ab04f20d9da..67706e2ce7ac5291f34e9524f9b77ae3425c06ad 100755
--- a/source/branches/CLASSV3/src/Makefile
+++ b/source/branches/CLASSV3/src/Makefile
@@ -28,7 +28,7 @@ OBJS = CLASS.o \
DynamicalSystem.o\
EvolutionData.o EvolutionDataDict.o \
IrradiationModel.o IM_RK4.o \
- EquivalenceModel.o EQM_MLP_PWR_MOX.o \
+ EquivalenceModel.o EQM_MLP_PWR_MOX.o EQM_QUAD_PWR_MOX.o EQM_LIN_PWR_MOX.o \
XSModel.o XSM_MLP_PWR_MOX.o XSM_CLOSEST.o \
PhysicModels.o \
LogFile.o