From c126c0ed982e284abf7ce0d03482c112d6d34b83 Mon Sep 17 00:00:00 2001
From: Baptiste Mouginot <mouginot.baptiste@gmail.com>
Date: Wed, 9 Jul 2014 13:03:25 +0000
Subject: [PATCH] git-svn-id: svn+ssh://svn.in2p3.fr/class@328
0e7d625b-0364-4367-a6be-d5be4a48d228
---
.../branches/CLASSV3/src/EQM_LIN_PWR_MOX.cxx | 262 ---------
.../branches/CLASSV3/src/EQM_MLP_PWR_MOX.cxx | 182 -------
.../branches/CLASSV3/src/EQM_QUAD_PWR_MOX.cxx | 133 -----
source/branches/CLASSV3/src/IM_Matrix.cxx | 315 -----------
source/branches/CLASSV3/src/IM_RK4.cxx | 408 --------------
source/branches/CLASSV3/src/XSM_CLOSEST.cxx | 366 -------------
source/branches/CLASSV3/src/XSM_MLP.cxx | 509 ------------------
source/branches/CLASSV3/src/XSModel.cxx | 11 +-
8 files changed, 8 insertions(+), 2178 deletions(-)
delete mode 100644 source/branches/CLASSV3/src/EQM_LIN_PWR_MOX.cxx
delete mode 100644 source/branches/CLASSV3/src/EQM_MLP_PWR_MOX.cxx
delete mode 100644 source/branches/CLASSV3/src/EQM_QUAD_PWR_MOX.cxx
delete mode 100644 source/branches/CLASSV3/src/IM_Matrix.cxx
delete mode 100644 source/branches/CLASSV3/src/IM_RK4.cxx
delete mode 100644 source/branches/CLASSV3/src/XSM_CLOSEST.cxx
delete mode 100644 source/branches/CLASSV3/src/XSM_MLP.cxx
diff --git a/source/branches/CLASSV3/src/EQM_LIN_PWR_MOX.cxx b/source/branches/CLASSV3/src/EQM_LIN_PWR_MOX.cxx
deleted file mode 100644
index 0d6dbaca0..000000000
--- a/source/branches/CLASSV3/src/EQM_LIN_PWR_MOX.cxx
+++ /dev/null
@@ -1,262 +0,0 @@
-#include "EQM_LIN_PWR_MOX.hxx"
-
-#include "CLASSConstante.hxx"
-
-#include <vector>
-
-#include "StringLine.hxx"
-#include "CLASSLogger.hxx"
-#include "IsotopicVector.hxx"
-
-
-
-EQM_LIN_PWR_MOX::EQM_LIN_PWR_MOX(string WeightPath):EquivalenceModel(new CLASSLogger("EQM_LIN_PWR_MOX.log"))
-{
- fWeightPath = WeightPath;
-
- ifstream DataDB(fWeightPath.c_str()); // Open the File
- if(!DataDB)
- WARNING << "Can't open \"" << fWeightPath << "\"\n" << endl;
-
- string line;
- int start = 0; // First Get Fuel Parameter
- getline(DataDB, line);
-
- if( StringLine::NextWord(line, start, ' ') != "PARAM")
- {
- ERROR << " Bad Database file : " << fWeightPath << " Can't find the Parameter of the DataBase " << endl;
- exit (1);
- }
- while(start < (int)line.size())
- fFuelParameter.push_back(atof(StringLine::NextWord(line, start, ' ').c_str()));
-
- INFO << 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(CLASSLogger* log, string WeightPath):EquivalenceModel(log)
-{
- fWeightPath = WeightPath;
-
- ifstream DataDB(fWeightPath.c_str()); // Open the File
- if(!DataDB)
- WARNING << " Can't open \"" << fWeightPath << "\"\n" << endl;
-
- string line;
- int start = 0; // First Get Fuel Parameter
- getline(DataDB, line);
-
- if( StringLine::NextWord(line, start, ' ') != "PARAM")
- {
- ERROR << " Bad Database file : " << fWeightPath << " Can't find the Parameter of the DataBase"<< endl;
- exit (1);
- }
- while(start < (int)line.size())
- fFuelParameter.push_back(atof(StringLine::NextWord(line, start, ' ').c_str()));
-
- INFO << fFuelParameter.size() << " have been read"<< endl;
-
-
-
- //-----------------------------------------------------------------------------//
- //-----------------------------------------------------------------------------//
- //-----------------------------------------------------------------------------//
- //-----------------------------------------------------------------------------//
- //-----------------------------------------------------------------------------//
- // ADD ENrichment of the U reading !!!!!!!!!!!!!!!!!!!!!!!!!!!! //
- //-----------------------------------------------------------------------------//
- //-----------------------------------------------------------------------------//
- //-----------------------------------------------------------------------------//
- //-----------------------------------------------------------------------------//
-
- ZAI U8(92,238,0);
- ZAI U5(92,235,0);
- double U5_enrich= 0.0025;
- fFertileList = U5*U5_enrich + U8*(1-U5_enrich);
-
-
- ZAI Pu8(94,238,0);
- ZAI Pu9(94,239,0);
- ZAI Pu0(94,240,0);
- ZAI Pu1(94,241,0);
- ZAI Pu2(94,242,0);
- fFissileList = Pu8*1+Pu9*1+Pu0*1+Pu1*1+Pu2*1;
-
-
-}
-
-EQM_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;
- if(FissilArray.size() == 0)
- {
- for(int i = 0; i < (int)lambda.size(); i++)
- lambda[i] = -1;
-
- FuelBuild = true;
- }
- int N_FissilStock_OnCheck = 0;
-
- while(!FuelBuild)
- {
-
- double nPu_0 = 0;
- double MPu_0 = 0;
- {
- map<ZAI ,double>::iterator it;
-
- map<ZAI ,double> isotopicquantity = FullUsedStock.GetSpeciesComposition(94).GetIsotopicQuantity();
- for( it = isotopicquantity.begin(); it != isotopicquantity.end(); it++ )
- nPu_0 += (*it).second;
-
- 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.GetMass( ZAI(92,238,0) ) * 0.997
- + cZAIMass.GetMass( ZAI(92,235,0) ) * 0.003 );
-
- double N1 = (BurnUp - fFuelParameter[6]) * NT;
- double N2 = -Sum_AlphaI_nPuI0;
- double N3 = -fFuelParameter[0] * Na / (cZAIMass.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)
- {
- WARNING << "!!!FabricationPlant!!! Oups Bug in calculating stock fraction to use "<< endl;
- lambda[N_FissilStock_OnCheck] = 0.;
- N_FissilStock_OnCheck++;
- FuelBuild = false;
- }
- else if( StockFactionToUse > 1 )
- {
-
- FullUsedStock += stock;
- lambda[N_FissilStock_OnCheck] = 1;
- N_FissilStock_OnCheck++;
- FuelBuild = false;
- }
- else
- {
- lambda[N_FissilStock_OnCheck] = StockFactionToUse;
-
- FuelBuild = true;
-
- double U8_Quantity = (HMMass - (MPu_0+StockFactionToUse*MPu_1 ))/(cZAIMass.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();
- }
-
-
- if( N_FissilStock_OnCheck == (int) FissilArray.size() ) // Check if the last Fissil stock has been tested... quit if so...
- {
- for(int i = 0; i < (int)lambda.size(); i++)
- lambda[i] = -1;
-
- FuelBuild = true;
- }
- }
-
-
-
- return lambda;
-}
diff --git a/source/branches/CLASSV3/src/EQM_MLP_PWR_MOX.cxx b/source/branches/CLASSV3/src/EQM_MLP_PWR_MOX.cxx
deleted file mode 100644
index b604ff50f..000000000
--- a/source/branches/CLASSV3/src/EQM_MLP_PWR_MOX.cxx
+++ /dev/null
@@ -1,182 +0,0 @@
-#include "EquivalenceModel.hxx"
-#include "EQM_MLP_PWR_MOX.hxx"
-#include "CLASSLogger.hxx"
-#include "StringLine.hxx"
-
-#include <string>
-#include <iostream>
-#include <fstream>
-#include <algorithm>
-#include <cmath>
-#include <cassert>
-
-#include "TSystem.h"
-#include "TMVA/Reader.h"
-#include "TMVA/Tools.h"
-#include "TMVA/MethodCuts.h"
-
-
-//________________________________________________________________________
-//
-// EQM_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):EquivalenceModel(new CLASSLogger("EQM_MLP_MOX.log"))
-{
- fTMVAWeightPath = TMVAWeightPath;
-
- ZAI U8(92,238,0);
- ZAI U5(92,235,0);
- double U5_enrich= 0.0025;
-
- ZAI Pu8(94,238,0);
- ZAI Pu9(94,239,0);
- ZAI Pu0(94,240,0);
- ZAI Pu1(94,241,0);
- ZAI Pu2(94,242,0);
-
- fFissileList = Pu8*1+Pu9*1+Pu0*1+Pu1*1+Pu2*1;
- fFertileList = U5*U5_enrich + U8*(1-U5_enrich);
-
- INFO<<"__An equivalence model of PWR MOX has been define__"<<endl;
- INFO<<"\tThis model is based on a multi layer perceptron"<<endl;
- INFO<<"\t\tThe TMVA weight file is :"<<endl;
- INFO<<"\t\t\t"<<fTMVAWeightPath<<endl;
-
-}
-
-//________________________________________________________________________
-EQM_MLP_MOX::EQM_MLP_MOX(CLASSLogger* log, string TMVAWeightPath):EquivalenceModel(log)
-{
- fTMVAWeightPath = TMVAWeightPath;
-
- ZAI U8(92,238,0);
- ZAI U5(92,235,0);
- double U5_enrich= 0.0025;
-
- ZAI Pu8(94,238,0);
- ZAI Pu9(94,239,0);
- ZAI Pu0(94,240,0);
- ZAI Pu1(94,241,0);
- ZAI Pu2(94,242,0);
-
- fFissileList = Pu8*1+Pu9*1+Pu0*1+Pu1*1+Pu2*1;
- fFertileList = U5*U5_enrich + U8*(1-U5_enrich);
-
- INFO<<"__An equivalence model of PWR MOX has been define__"<<endl;
- INFO<<"\tThis model is based on a multi layer perceptron"<<endl;
- INFO<<"\t\tThe TMVA weight file is :"<<endl;
- INFO<<"\t\t\t"<<fTMVAWeightPath<<endl;
-
-}
-
-//________________________________________________________________________
-TTree* EQM_MLP_MOX::CreateTMVAInputTree(IsotopicVector Fissil,IsotopicVector Fertil,double BurnUp)
-{
- TTree* InputTree = new TTree("EQTMP", "EQTMP");
- float Pu8 = 0;
- float Pu9 = 0;
- float Pu10 = 0;
- float Pu11 = 0;
- float Pu12 = 0;
- float Am1 = 0;
- float U5_enrichment = 0;
- float BU = 0;
-
- InputTree->Branch( "Pu8" ,&Pu8 ,"Pu8/F" );
- InputTree->Branch( "Pu9" ,&Pu9 ,"Pu9/F" );
- InputTree->Branch( "Pu10" ,&Pu10 ,"Pu10/F" );
- InputTree->Branch( "Pu11" ,&Pu11 ,"Pu11/F" );
- InputTree->Branch( "Pu12" ,&Pu12 ,"Pu12/F" );
- InputTree->Branch( "Am1" ,&Am1 ,"Am1/F" );
- InputTree->Branch( "U5_enrichment" ,&U5_enrichment ,"U5_enrichment/F" );
- InputTree->Branch( "BU" ,&BU ,"BU/F" );
-
-
- float U8 = 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;
-//cout<<"Pu8 "<<Pu8 <<" Pu9 "<< Pu9 <<" Pu10 "<< Pu10 << " Pu11 "<< Pu11 <<" Pu12 "<<Pu12 <<" Am1 "<<Am1<<endl;
-//cout<<"BU "<<BU<<" U5_enrichment "<<U5_enrichment<<endl;
- if(Pu8 + Pu9 + Pu10 + Pu11 + Pu12 + Am1 > 1.00001 )//?????1.00001??? I don't know it! goes in condition if =1 !! may be float/double issue ...
- {
- ERROR << Pu8 << " " << Pu9 << " " << Pu10 << " " << Pu11 << " " << Pu12 << " " << Am1 << endl;
- exit(0);
- }
- // All value are molar (!weight)
-
- InputTree->Fill();
- return InputTree;
-}
-//________________________________________________________________________
-double EQM_MLP_MOX::ExecuteTMVA(TTree* theTree)
-{
- // --- Create the Reader object
- TMVA::Reader *reader = new TMVA::Reader( "Silent" );
- // Create a set of variables and declare them to the reader
- // - the variable names MUST corresponds in name and type to those given in the weight file(s) used
- Float_t Pu8,Pu9,Pu10,Pu11,Pu12,Am1,BU,U5_enrichment;
-
- reader->AddVariable( "BU" ,&BU );
- reader->AddVariable( "U5_enrichment",&U5_enrichment );
- reader->AddVariable( "Pu8" ,&Pu8 );
- reader->AddVariable( "Pu9" ,&Pu9 );
- reader->AddVariable( "Pu10" ,&Pu10);
- reader->AddVariable( "Pu11" ,&Pu11);
- reader->AddVariable( "Pu12" ,&Pu12);
- reader->AddVariable( "Am1" ,&Am1 );
-
- // --- Book the MVA methods
-
- // Book method MLP
- TString methodName = "MLP method";
- reader->BookMVA( methodName, fTMVAWeightPath );
- //theTree->SetBranchAddress("teneur",&teneur);
- theTree->SetBranchAddress( "BU" ,&BU );
- theTree->SetBranchAddress( "U5_enrichment" ,&U5_enrichment ) ;
- theTree->SetBranchAddress( "Pu8" ,&Pu8 );
- theTree->SetBranchAddress( "Pu9" ,&Pu9 );
- theTree->SetBranchAddress( "Pu10" ,&Pu10 );
- theTree->SetBranchAddress( "Pu11" ,&Pu11 );
- theTree->SetBranchAddress( "Pu12" ,&Pu12 );
- theTree->SetBranchAddress( "Am1" ,&Am1 );
- theTree->GetEntry(0);
-
- Float_t val = (reader->EvaluateRegression( methodName ))[0];
-
- delete reader;
- delete theTree;
-
- return (double)val; //retourne teneur
-}
-//________________________________________________________________________
-double EQM_MLP_MOX::GetFissileMolarFraction(IsotopicVector Fissil,IsotopicVector Fertil,double BurnUp)
-{DBGL
- return ExecuteTMVA(CreateTMVAInputTree(Fissil,Fertil,BurnUp));
-}
diff --git a/source/branches/CLASSV3/src/EQM_QUAD_PWR_MOX.cxx b/source/branches/CLASSV3/src/EQM_QUAD_PWR_MOX.cxx
deleted file mode 100644
index 79ce73946..000000000
--- a/source/branches/CLASSV3/src/EQM_QUAD_PWR_MOX.cxx
+++ /dev/null
@@ -1,133 +0,0 @@
-#include "EQM_QUAD_PWR_MOX.hxx"
-
-#include <vector>
-
-#include "StringLine.hxx"
-#include "CLASSLogger.hxx"
-
-
-
-
-EQM_QUAD_PWR_MOX::EQM_QUAD_PWR_MOX(string WeightPath):EquivalenceModel(new CLASSLogger("EQM_QUAD_PWR_MOX.log"))
-{
- fWeightPath = WeightPath;
-
- ifstream DataDB(fWeightPath.c_str()); // Open the File
- if(!DataDB)
- WARNING << " Can't open \"" << fWeightPath << "\"" << endl;
-
- string line;
- int start = 0; // First Get Fuel Parameter
- getline(DataDB, line);
-
- if( StringLine::NextWord(line, start, ' ') != "PARAM")
- {
- ERROR << "Bad Database file : " << fWeightPath << " Can't find the Parameter of the DataBase " << endl;
- exit (1);
- }
- while(start < (int)line.size())
- fFuelParameter.push_back(atof(StringLine::NextWord(line, start, ' ').c_str()));
-
- INFO << fFuelParameter.size() << " have been read " << endl;
-
-
- ZAI U8(92,238,0);
- ZAI U5(92,235,0);
- double U5_enrich= 0.0025;
- fFertileList = U5*U5_enrich + U8*(1-U5_enrich);
-
-
- ZAI Pu8(94,238,0);
- ZAI Pu9(94,239,0);
- ZAI Pu0(94,240,0);
- ZAI Pu1(94,241,0);
- ZAI Pu2(94,242,0);
- fFissileList = Pu8*1+Pu9*1+Pu0*1+Pu1*1+Pu2*1;
-
-
-}
-
-EQM_QUAD_PWR_MOX::EQM_QUAD_PWR_MOX(CLASSLogger* log, string WeightPath):EquivalenceModel(log)
-{
- fWeightPath = WeightPath;
-
- ifstream DataDB(fWeightPath.c_str()); // Open the File
- if(!DataDB)
- WARNING << " Can't open \"" << fWeightPath << "\"" << endl;
-
- string line;
- int start = 0; // First Get Fuel Parameter
- getline(DataDB, line);
-
- if( StringLine::NextWord(line, start, ' ') != "PARAM")
- {
- ERROR << "Bad Database file : " << fWeightPath << " Can't find the Parameter of the DataBase " << endl;
- exit (1);
- }
- while(start < (int)line.size())
- fFuelParameter.push_back(atof(StringLine::NextWord(line, start, ' ').c_str()));
-
- INFO << fFuelParameter.size() << " have been read " << endl;
-
-
- ZAI U8(92,238,0);
- ZAI U5(92,235,0);
- double U5_enrich= 0.0025;
- fFertileList = U5*U5_enrich + U8*(1-U5_enrich);
-
-
- ZAI Pu8(94,238,0);
- ZAI Pu9(94,239,0);
- ZAI Pu0(94,240,0);
- ZAI Pu1(94,241,0);
- ZAI Pu2(94,242,0);
- fFissileList = Pu8*1+Pu9*1+Pu0*1+Pu1*1+Pu2*1;
-
-
-}
-
-
-
-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/IM_Matrix.cxx b/source/branches/CLASSV3/src/IM_Matrix.cxx
deleted file mode 100644
index 177e9b146..000000000
--- a/source/branches/CLASSV3/src/IM_Matrix.cxx
+++ /dev/null
@@ -1,315 +0,0 @@
-//
-// 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 "CLASSLogger.hxx"
-#include "StringLine.hxx"
-
-#include <TGraph.h>
-#include <TString.h>
-
-
-#include <sstream>
-#include <string>
-#include <iostream>
-#include <fstream>
-#include <algorithm>
-#include <cmath>
-
-
-
-using namespace std;
-
-
-//________________________________________________________________________
-IM_Matrix::IM_Matrix():IrradiationModel(new CLASSLogger("IM_Matrix.log")), DynamicalSystem()
-{
- fShorstestHalflife = 3600.*24*160.; //cut by default all nuclei with a shorter liftime than the Cm242 -> remain 33 actinides
-}
-
-
-IM_Matrix::IM_Matrix(CLASSLogger* log):IrradiationModel(log), 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)
-{
-DBGL
- 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 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/AVOGADRO*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
-DBGL
- for(int i = 0; i < NStep-1; i++)
- {
- double TStepMax = ( (DBTimeStep[i+1]-DBTimeStep[i] ) ) * Power_ref/M_ref / Power*M ; // Get the next Time step
-
-
- TMatrixT<double> BatemanMatrix = TMatrixT<double>(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();
-
-
- }
-DBGL
- FissionXSMatrix.push_back(GetFissionXsMatrix(XSSet, DBTimeStep[NStep-1])); //Feel the reaction Matrix
- CaptureXSMatrix.push_back(GetCaptureXsMatrix(XSSet, DBTimeStep[NStep-1])); //Feel the reaction Matrix
- n2nXSMatrix.push_back(Getn2nXsMatrix(XSSet, DBTimeStep[NStep-1])); //Feel the reaction Matrix
-
-
- EvolutionData GeneratedDB = EvolutionData(GetLog());
-
- double ESigmaN = 0;
- for (int j = 0; j < (int)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.SetReactorType(ReactorType );
- GeneratedDB.SetCycleTime(cycletime);
-
- for (int i = 0; i < (int) FissionXSMatrix.size(); i++)
- {
- FissionXSMatrix[i].Clear();
- CaptureXSMatrix[i].Clear();
- n2nXSMatrix[i].Clear();
- }
- FissionXSMatrix.clear();
- CaptureXSMatrix.clear();
- n2nXSMatrix.clear();
-DBGL
- return GeneratedDB;
-
-}
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
diff --git a/source/branches/CLASSV3/src/IM_RK4.cxx b/source/branches/CLASSV3/src/IM_RK4.cxx
deleted file mode 100644
index 95ac16189..000000000
--- a/source/branches/CLASSV3/src/IM_RK4.cxx
+++ /dev/null
@@ -1,408 +0,0 @@
-//
-// 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 "CLASSConstante.hxx"
-#include "CLASSLogger.hxx"
-
-#include "StringLine.hxx"
-
-#include <TGraph.h>
-#include <TString.h>
-
-
-#include <sstream>
-#include <string>
-#include <iostream>
-#include <fstream>
-#include <algorithm>
-#include <cmath>
-
-
-
-using namespace std;
-
-
-//________________________________________________________________________
-IM_RK4::IM_RK4():IrradiationModel(new CLASSLogger("IM_RK4.log")), DynamicalSystem()
-{
- fTheNucleiVector = 0;
- fTheMatrix = 0;
-
- SetForbidNegativeValue();
-
-}
-
-
-IM_RK4::IM_RK4(CLASSLogger* log):IrradiationModel(log), DynamicalSystem()
-{
- fTheNucleiVector = 0;
- fTheMatrix = 0;
-
- SetForbidNegativeValue();
-
-}
-
-
-
-
-
-
-
-
-
-//________________________________________________________________________
-/* Evolution Calculation */
-//________________________________________________________________________
-EvolutionData IM_RK4::GenerateEvolutionData(IsotopicVector isotopicvector, EvolutionData XSSet, double Power, double cycletime)
-{
-DBGL
- 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 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/AVOGADRO*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
-DBGL
- for(int i = 0; i < NStep-1; i++)
- {
- double TStepMax = ( (DBTimeStep[i+1]-DBTimeStep[i] ) ) * Power_ref/M_ref / Power*M ; // Get the next Time step
-
-
- TMatrixT<double> BatemanMatrix = TMatrixT<double>(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)));
- }
-DBGL
- GeneratedDB.SetPower(Power );
- 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();
-DBGL
- 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/src/XSM_CLOSEST.cxx b/source/branches/CLASSV3/src/XSM_CLOSEST.cxx
deleted file mode 100644
index 6cf7cf311..000000000
--- a/source/branches/CLASSV3/src/XSM_CLOSEST.cxx
+++ /dev/null
@@ -1,366 +0,0 @@
-#include "XSModel.hxx"
-#include "XSM_CLOSEST.hxx"
-#include "CLASSLogger.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(string DB_index_file, bool oldreadmethod ): XSModel(new CLASSLogger("XSM_CLOSEST.log"))
-{
- fOldReadMethod = oldreadmethod;
- fDataBaseIndex = DB_index_file;
- fDistanceType = 0;
- fWeightedDistance = false;
- fEvolutionDataInterpolation = false;
- ReadDataBase();
-
- // Warning
- INFO << " A EvolutionData has been define :" << endl;
- INFO << "\t His index is : \"" << DB_index_file << "\" " << endl;
- INFO << "\t " << fFuelDataBank.size() << " EvolutionData have been read."<< endl << endl;
-
-}
-
-//________________________________________________________________________
-XSM_CLOSEST::XSM_CLOSEST(CLASSLogger* Log,string DB_index_file, bool oldreadmethod ): XSModel(Log)
-{
- fOldReadMethod = oldreadmethod;
- fDataBaseIndex = DB_index_file;
- fDistanceType = 0;
- fWeightedDistance = false;
- fEvolutionDataInterpolation = false;
- ReadDataBase();
-
- // Warning
- INFO << " A EvolutionData has been define :" << endl;
- INFO << "\t His index is : \"" << DB_index_file << "\" " << endl;
- INFO << "\t " << fFuelDataBank.size() << " EvolutionData have been read."<< endl << endl;
-
-}
-
-//________________________________________________________________________
-XSM_CLOSEST::~XSM_CLOSEST()
-{
- for( int i = 0; i < (int)fFuelDataBank.size(); i++)
- fFuelDataBank[i].DeleteEvolutionData();
- fFuelDataBank.clear();
-}
-
-//________________________________________________________________________
-void XSM_CLOSEST::ReadDataBase()
-{
-
- if(fFuelDataBank.size() != 0)
- {
- for( int i = 0; i < (int)fFuelDataBank.size(); i++)
- fFuelDataBank[i].DeleteEvolutionData();
- fFuelDataBank.clear();
- }
-
-
- ifstream DataDB(fDataBaseIndex.c_str()); // Open the File
- if(!DataDB)
- WARNING << " 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")
- {
- ERROR << " 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(GetLog(), line, fOldReadMethod);
- fFuelDataBank.push_back(evolutionproduct);
- }
- }
-
-}
-//________________________________________________________________________
-//________________________________________________________________________
-//________________________________________________________________________
-/* Distance Calculation */
-//________________________________________________________________________
-//________________________________________________________________________
-//________________________________________________________________________
-map<double, int> XSM_CLOSEST::GetDistancesTo(IsotopicVector isotopicvector, double t)
-{
-
- map<double, int> distances;
-
- for( int i = 0; i < (int)fFuelDataBank.size(); i++)
- {
- pair<map<double, int>::iterator, bool> IResult;
-
- IsotopicVector ActinidesCompositionAtT = fFuelDataBank[i].GetIsotopicVectorAt(t).GetActinidesComposition();
- IsotopicVector IV_ActinidesComposition = isotopicvector.GetActinidesComposition();
-
- double NormalisationFactor = Norme(IV_ActinidesComposition) / Norme( ActinidesCompositionAtT );
-
-
- double distance = Distance( IV_ActinidesComposition,
- ActinidesCompositionAtT / NormalisationFactor,
- fDistanceType,
- fDistanceParameter);
-
- IResult = distances.insert( pair< double, int >(distance, i) );
- }
-
- return distances;
-
-}
-//________________________________________________________________________
-EvolutionData XSM_CLOSEST::GetCrossSections(IsotopicVector isotopicvector, double t)
-{
-DBGL
- double distance = 0;
- int N_BestEvolutionData = 0;
-
-
- if(fWeightedDistance)
- {
-DBGL
- IsotopicVector ActinidesCompositionAtT = fFuelDataBank[0].GetIsotopicVectorAt(t).GetActinidesComposition();
- IsotopicVector IV_ActinidesComposition = isotopicvector.GetActinidesComposition();
-
- double NormalisationFactor = Norme( ActinidesCompositionAtT ) / Norme(IV_ActinidesComposition);
-
-
- distance = Distance( IV_ActinidesComposition / NormalisationFactor,
- fFuelDataBank[0]);
-
-
- for( int i = 1; i < (int)fFuelDataBank.size(); i++)
- {
- double D = 0;
- ActinidesCompositionAtT = fFuelDataBank[i].GetIsotopicVectorAt(t).GetActinidesComposition();
- IV_ActinidesComposition = isotopicvector.GetActinidesComposition();
-
- D = Distance( IV_ActinidesComposition / NormalisationFactor,
- fFuelDataBank[i]);
-
-
- if (D< distance)
- {
- distance = D;
- N_BestEvolutionData = i;
- }
- }
-DBGL
- return fFuelDataBank[N_BestEvolutionData];
- }
- else if (fEvolutionDataInterpolation)
- {
-DBGL
- map<double, int> distance_map = GetDistancesTo(isotopicvector, t);
- map<double, int>::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++)
- {
-
- double distance = (*it_distance).first;
- int ED_Indice = (*it_distance).second;
-
- if(distance == 0 )
- {
- EvolInterpolate = Multiply(1,fFuelDataBank[ED_Indice]);
- return EvolInterpolate;
- }
-
- if(Nstep == 0)
- EvolInterpolate = Multiply(1./distance, fFuelDataBank[ED_Indice]);
- else
- {
-
- EvolutionData Evoltmp = EvolInterpolate;
- EvolutionData Evoltmp2 = Multiply(1./distance, fFuelDataBank[ED_Indice]);
-
- EvolInterpolate = Sum(Evoltmp, Evoltmp2);
- Evoltmp.DeleteEvolutionData();
- Evoltmp2.DeleteEvolutionData();
-
-
- }
-
- SumOfDistance += 1./distance;
- Nstep++;
- if(Nstep == NClose) break;
-
- }
-
- EvolutionData Evoltmp = EvolInterpolate;
- EvolInterpolate.Clear();
-
- EvolInterpolate = Multiply(1/SumOfDistance, Evoltmp);
-
- Evoltmp.DeleteEvolutionData();
-DBGL
- return EvolInterpolate;
-
- }
- else
- {
-DBGL
- IsotopicVector ActinidesCompositionAtT = fFuelDataBank[0].GetIsotopicVectorAt(t).GetActinidesComposition();
- IsotopicVector IV_ActinidesComposition = isotopicvector.GetActinidesComposition();
-
- double NormalisationFactor = Norme(IV_ActinidesComposition) / Norme( ActinidesCompositionAtT );
-
-
- distance = Distance( IV_ActinidesComposition,
- ActinidesCompositionAtT / NormalisationFactor,
- fDistanceType,
- fDistanceParameter);
-
- for( int i = 1; i < (int)fFuelDataBank.size(); i++)
- {
- double D = 0;
- ActinidesCompositionAtT = fFuelDataBank[i].GetIsotopicVectorAt(t).GetActinidesComposition();
- IV_ActinidesComposition = isotopicvector.GetActinidesComposition();
-
- D = Distance( IV_ActinidesComposition,
- ActinidesCompositionAtT / NormalisationFactor,
- fDistanceType,
- fDistanceParameter);
-
- if (D< distance)
- {
- distance = D;
- N_BestEvolutionData = i;
- }
- }
-
-DBGL
- return fFuelDataBank[N_BestEvolutionData];
- }
-
-}
-//________________________________________________________________________
-void XSM_CLOSEST::CalculateDistanceParameter()
-{
-DBGL
- if(fDistanceType!=1){
- WARNING << " 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.
- int NevolutionDatainFuelDataBank = 0;
-
- for( int i = 0; i < (int)fFuelDataBank.size(); i++)
- {
- NevolutionDatainFuelDataBank++;
- map<ZAI ,double>::iterator itit;
- map<ZAI ,double> isovector = fFuelDataBank[i].GetIsotopicVectorAt(0).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+= fFuelDataBank[i].GetXSForAt(0, (*itit).first, i);
-
- fDistanceParameter.Add((*itit).first,TmpXS);
- }
-
-
- }
- fDistanceParameter.Multiply( (double)1.0/NevolutionDatainFuelDataBank );
-
- if(GetLog())
- {
- INFO <<"!!INFO!! Distance Parameters "<<endl;
- map<ZAI ,double >::iterator it2;
- for(it2 = fDistanceParameter.GetIsotopicQuantity().begin();it2 != fDistanceParameter.GetIsotopicQuantity().end(); it2++)
- {
- INFO << (*it2).first.Z() << " ";
- INFO << (*it2).first.A() << " ";
- INFO << (*it2).first.I() << " ";
- INFO << ": " << (*it2).second;
- INFO << endl;
- }
- INFO << endl;
- }
-
-DBGL
-}
-//________________________________________________________________________
-void XSM_CLOSEST::SetDistanceParameter(IsotopicVector DistanceParameter)
-{
-
- fDistanceParameter = DistanceParameter;
-
- INFO <<"!!INFO!! Distance Parameters "<<endl;
- map<ZAI ,double >::iterator it2;
- for(it2 = fDistanceParameter.GetIsotopicQuantity().begin();it2 != fDistanceParameter.GetIsotopicQuantity().end(); it2++)
- {
- INFO << (*it2).first.Z() << " ";
- INFO << (*it2).first.A() << " ";
- INFO << (*it2).first.I() << " ";
- INFO << ": " << (*it2).second;
- INFO << endl;
- }
- INFO << 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....
- WARNING << " 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 ){
- ERROR << " 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/src/XSM_MLP.cxx b/source/branches/CLASSV3/src/XSM_MLP.cxx
deleted file mode 100644
index 5dbf23d7f..000000000
--- a/source/branches/CLASSV3/src/XSM_MLP.cxx
+++ /dev/null
@@ -1,509 +0,0 @@
-
-#include "XSModel.hxx"
-#include "XSM_MLP.hxx"
-#include "CLASSLogger.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 "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
-//
-//
-//
-//
-//________________________________________________________________________
-XSM_MLP::XSM_MLP(string TMVA_Weight_Directory,string InformationFile, bool IsTimeStep):XSModel(new CLASSLogger("XSM_MLP.log"))
-{
-
- fIsStepTime=IsTimeStep;
- fTMVAWeightFolder = TMVA_Weight_Directory;
- if(InformationFile=="")
- fMLPInformationFile = TMVA_Weight_Directory+"/Data_Base_Info.nfo";
- else
- fMLPInformationFile=fTMVAWeightFolder+InformationFile;
-
- GetMLPWeightFiles();
- GetDataBaseInformation();
-
- INFO<<"__A cross section interpolator using" <<endl;
- INFO<<"Multi Layer Perceptron has been define__"<<endl;
- INFO << " \t His TMVA folder is : \" " << fTMVAWeightFolder << "\"" << endl;
-
-}
-//________________________________________________________________________
-XSM_MLP::XSM_MLP(CLASSLogger* Log,string TMVA_Weight_Directory,string InformationFile, bool IsTimeStep):XSModel(Log)
-{
-
- fIsStepTime=IsTimeStep;
- fTMVAWeightFolder = TMVA_Weight_Directory;
- if(InformationFile=="")
- fMLPInformationFile = TMVA_Weight_Directory+"/Data_Base_Info.nfo";
- else
- fMLPInformationFile=fTMVAWeightFolder+InformationFile;
-
- GetMLPWeightFiles();
- GetDataBaseInformation();
-
- INFO<<"__A cross section interpolator using" <<endl;
- INFO<<"Multi Layer Perceptron has been define__"<<endl;
- INFO << " \t His TMVA folder is : \" " << fTMVAWeightFolder << "\"" << endl;
-
-}
-//________________________________________________________________________
-XSM_MLP::~XSM_MLP()
-{
- fMapOfTMVAVariableNames.clear();
-}
-//________________________________________________________________________
-void XSM_MLP::GetDataBaseInformation()
-{
- ifstream FILE(fMLPInformationFile.c_str());
-
- if(FILE.good())
- {
- while(!FILE.eof())
- {
- string line;
- getline(FILE, line);
- size_t foundRType = line.find("ReactorType :");
- size_t foundFType = line.find("FuelType :");
- size_t foundHM = line.find("Heavy Metal (t) :");
- size_t foundPower = line.find("Thermal Power (W) :");
- size_t foundTime = line.find("Time (s) :");
- size_t foundZAI = line.find("Z A I Name (input MLP) :");
- int pos=0;
- if(foundRType != std::string::npos)
- { StringLine::NextWord(line,pos,':');
- fDataBaseRType = atof( (StringLine::NextWord(line,pos,':')).c_str() );
- }
- if(foundFType != std::string::npos)
- { StringLine::NextWord(line,pos,':');
- fDataBaseFType = atof( (StringLine::NextWord(line,pos,':')).c_str() );
- }
- if(foundHM != std::string::npos)
- { StringLine::NextWord(line,pos,':');
- fDataBaseHMMass = atof( (StringLine::NextWord(line,pos,':')).c_str() );
- }
- if(foundPower !=std::string::npos)
- { StringLine::NextWord(line,pos,':');
- fDataBasePower = atof( (StringLine::NextWord(line,pos,':') ).c_str() );
- }
- pos=0;
- if(foundTime!=std::string::npos)
- {
- StringLine::NextWord(line,pos,':');
- while( pos< (int)line.size() )
- fMLP_Time.push_back( atof( (StringLine::NextWord(line,pos,' ')).c_str() ));
- }
-
- if(foundZAI != std::string::npos)
- {
- while(!FILE.eof())
- {
- int Z=-4;
- int A=-4;
- int I=-4;
- string Name;
- getline(FILE, line);
- stringstream ssline;
- ssline<<line;
- ssline>>Z>>A>>I>>Name;
- //cout<<Z<<" "<<A<<" "<<I<<" "<<Name<<endl;
- fMapOfTMVAVariableNames.insert( pair<ZAI,string>(ZAI(Z,A,I),Name) );
-
- }
- }
-
- }
- }
- else
- {
- ERROR << "Can't find/open file " << fMLPInformationFile << endl;
- exit(0);
- }
-
- /********DEBUG*************************************/
- INFO<<"\tMLP XS Data Base Information : "<<endl;
- INFO<<"\t\tHeavy Metal (t) :"<<fDataBaseHMMass<<endl;
- INFO<<"\t\tThermal Power (W) :"<<fDataBasePower<<endl;
- INFO<<"\t\tTime (s) :"<<endl;
- for (int i = 0; i < (int)fMLP_Time.size(); ++i)
- INFO<<"\t\t\t"<<fMLP_Time[i]<<endl;
- INFO<<"\t\tZ A I Name (input MLP) :"<<endl;
-
- map<ZAI ,string >::iterator it;
-
- for (it= fMapOfTMVAVariableNames.begin();it!=fMapOfTMVAVariableNames.end();it++)
- INFO<<"\t\t\t"<< it->first.Z()<<" "<<it->first.A()<<" "<<it->second<<endl;
-
-
-}
-//________________________________________________________________________
-void XSM_MLP::GetMLPWeightFiles()
-{DBGL
- /**********Get All TMVA weight files*******************/
- //check if the folder containing weights exists
- DIR* rep = NULL;
- struct dirent* fichierLu = NULL;
- rep = opendir(fTMVAWeightFolder.c_str());
- if (rep == NULL)
- {
- 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);
-
- }
- }
-DBGL
-}
-//________________________________________________________________________
-//________________________________________________________________________
-//
-// Time (MLP take time as parameter)
-//________________________________________________________________________
-//________________________________________________________________________
-void XSM_MLP::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() ;
- size_t foundext = sReaction.find(".weights.xml");
- sReaction = sReaction.substr(0,foundext);
-
- if(sReaction=="fis")
- Reaction=0;
- if(sReaction=="cap")
- Reaction=1;
- if(sReaction=="n2n")
- Reaction=2;
-
- if(Z<=0 || A<=0 || I<0 || Reaction==-1)
- {
- ERROR << " wrong TMVA weight format " << endl;
- exit(0);
- }
-
-}
-//________________________________________________________________________
-TTree* XSM_MLP::CreateTMVAInputTree(IsotopicVector isotopicvector,int TimeStep)
-{
- /******Create Input data tree to be interpreted by TMVA::Reader***/
- TTree* InputTree = new TTree("InTMP", "InTMP");
-
- vector<float> InputTMVA;
- for(int i = 0 ; i< (int)fMapOfTMVAVariableNames.size() ; i++)
- InputTMVA.push_back(0);
-
- float Time=0;
-
- IsotopicVector IVInputTMVA;
- map<ZAI ,string >::iterator it;
- int j=0;
-
- for( it = fMapOfTMVAVariableNames.begin() ; it!=fMapOfTMVAVariableNames.end() ; it++)
- {
- InputTree->Branch( ((*it).second).c_str() ,&InputTMVA[j], ((*it).second + "/F").c_str());
- IVInputTMVA+= ((*it).first)*1;
- j++;
- }
-
- if( !fIsStepTime)
- InputTree->Branch( "Time" ,&Time ,"Time/F" );
-
- IsotopicVector IVAccordingToUserInfoFile = isotopicvector.GetThisComposition(IVInputTMVA);
-
- double Ntot = IVAccordingToUserInfoFile.GetTotalMass()*1e6/IVAccordingToUserInfoFile.MeanMolar()*6.02214129e23;
-
- IVAccordingToUserInfoFile = IVAccordingToUserInfoFile/Ntot;
-
- j=0;
- map<ZAI ,string >::iterator it2;
- for( it2 = fMapOfTMVAVariableNames.begin() ; it2!=fMapOfTMVAVariableNames.end() ; it2++)
- {
- InputTMVA[j] = IVAccordingToUserInfoFile.GetZAIIsotopicQuantity( (*it2).first ) ;
- INFO<< (*it2).first.Z()<<" "<<(*it2).first.A()<<" "<<InputTMVA[j]<<endl;
- j++;
- }
-
- Time=fMLP_Time[TimeStep];
-
- InputTree->Fill();
-
- return InputTree;
-}
-//________________________________________________________________________
-double XSM_MLP::ExecuteTMVA(string WeightFile,TTree* InputTree)
-{
- // --- Create the Reader object
- TMVA::Reader *reader = new TMVA::Reader( "Silent" );
-
- // Create a set of variables and declare them to the reader
- // - the variable names MUST corresponds in name and type to those given in the weight file(s) used
- vector<float> InputTMVA;
- for(int i = 0 ; i< (int)fMapOfTMVAVariableNames.size() ; i++)
- InputTMVA.push_back(0);
- Float_t Time;
-
- map<ZAI ,string >::iterator it;
- int j=0;
- for( it = fMapOfTMVAVariableNames.begin() ; it!=fMapOfTMVAVariableNames.end() ; it++)
- { reader->AddVariable( ( (*it).second ).c_str(),&InputTMVA[j]);
- j++;
- }
- if(!fIsStepTime)
- reader->AddVariable( "Time" ,&Time);
-
- // --- Book the MVA methods
-
- string dir = fTMVAWeightFolder;
- if(dir[dir.size()-1]!='/')
- dir+="/";
-
- // Book method MLP
- TString methodName = "MLP method";
- TString weightpath = dir + WeightFile ;
- reader->BookMVA( methodName, weightpath );
-
- map<ZAI ,string >::iterator it2;
- j=0;
- for( it2 = fMapOfTMVAVariableNames.begin() ; it2!=fMapOfTMVAVariableNames.end() ; it2++)
- {
- InputTree->SetBranchAddress(( (*it2).second ).c_str(),&InputTMVA[j]);
- j++;
- }
-
- if(!fIsStepTime)
- InputTree->SetBranchAddress( "Time" ,&Time );
-
- InputTree->GetEntry(0);
- Float_t val = (reader->EvaluateRegression( methodName ))[0];
-
- delete reader;
-
- return (double)val;
-}
-//________________________________________________________________________
-EvolutionData XSM_MLP::GetCrossSectionsTime(IsotopicVector IV)
-{DBGL
-
- EvolutionData EvolutionDataFromMLP = EvolutionData();
-
- map<ZAI,TGraph*> ExtrapolatedXS[3];
- /*************DATA BASE INFO****************/
- EvolutionDataFromMLP.SetReactorType(fDataBaseRType);
- EvolutionDataFromMLP.SetFuelType(fDataBaseFType);
- 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);
-
- for(int TimeStep=0;TimeStep<int(fMLP_Time.size());TimeStep++)
- {
- TTree* InputTree = CreateTMVAInputTree(IV,TimeStep);
-
- pair< map<ZAI, TGraph*>::iterator, bool> IResult;
-
- IResult = ExtrapolatedXS[Reaction].insert(pair<ZAI ,TGraph* >(ZAI(Z,A,I), new TGraph() ) );
-
- double XSValue = ExecuteTMVA(fWeightFiles[i],InputTree );
- if(IResult.second )
- {
- (IResult.first)->second->SetPoint(0, (double)GetMLPTime()[TimeStep], XSValue );
-
- }
- else
- {
- (IResult.first)->second->SetPoint( (IResult.first)->second->GetN(), (double)GetMLPTime()[TimeStep], XSValue );
- }
-
- delete InputTree;
- }
-
- }
-
- /**********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]);
-
-
-DBGL
- return EvolutionDataFromMLP;
-}
-//________________________________________________________________________
-//________________________________________________________________________
-//
-// Time step (1 MLP per time step)
-//________________________________________________________________________
-//________________________________________________________________________
-void XSM_MLP::ReadWeightFileStep(string Filename, int &Z, int &A, int &I, int &Reaction, int &TimeStep)
-{
- Z=-1;
- A=-1;
- I=-1;
- Reaction=-1;
-
- size_t found = Filename.find("XS");
-
- string NameJOB;
- NameJOB=Filename.substr(found);
- int pos=0;
-
- StringLine::NextWord(NameJOB, pos, '_');
-
- Z = atof( (StringLine::NextWord(NameJOB,pos,'_') ).c_str() );
- A = atof( (StringLine::NextWord(NameJOB,pos,'_') ).c_str() );
- I = atof( (StringLine::NextWord(NameJOB,pos,'_') ).c_str() );
-
- string sReaction = (StringLine::NextWord(NameJOB,pos,'_') ).c_str() ;
- if(sReaction=="fis")
- Reaction=0;
- if(sReaction=="cap")
- Reaction=1;
- if(sReaction=="n2n")
- Reaction=2;
-
- TimeStep = atof( (StringLine::NextWord(NameJOB,pos,'_') ).c_str() );
-
- if(Z==-1 || A==-1 || I==-1 || Reaction==-1 || TimeStep==-1)
- {
- ERROR << " wrong TMVA weight format " << endl;
- exit(0);
- }
-}
-
-//________________________________________________________________________
-EvolutionData XSM_MLP::GetCrossSectionsStep(IsotopicVector IV)
-{DBGL
- TTree* InputTree=CreateTMVAInputTree(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], ExecuteTMVA(fWeightFiles[i],InputTree) );
- }
- else
- {
- (IResult.first)->second->SetPoint( (IResult.first)->second->GetN(), (double)GetMLPTime()[TimeStep], ExecuteTMVA(fWeightFiles[i],InputTree) );
- }
-
- }
-
- /**********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]);
-
- delete InputTree;
-DBGL
- return EvolutionDataFromMLP;
-}
-//________________________________________________________________________
-EvolutionData XSM_MLP::GetCrossSections(IsotopicVector IV ,double t)
-{DBGL
- if(t!=0)
- WARNING << " Argument t has non effect here " << endl;
-
- EvolutionData EV;
- if(fIsStepTime)
- EV=GetCrossSectionsStep(IV);
-
- else
- EV=GetCrossSectionsTime(IV);
-
-DBGL
- return EV;
-}
-//________________________________________________________________________
diff --git a/source/branches/CLASSV3/src/XSModel.cxx b/source/branches/CLASSV3/src/XSModel.cxx
index 606d7a125..723261331 100644
--- a/source/branches/CLASSV3/src/XSModel.cxx
+++ b/source/branches/CLASSV3/src/XSModel.cxx
@@ -7,12 +7,17 @@
//
#include "XSModel.hxx"
-#include "CLASSHeaders.hxx"
-
using namespace std;
-XSModel::XSModel(): CLASSObject()
+XSModel::XSModel(): CLASSObject(new CLASSLogger("XSModel.log"))
{
}
+
+
+
+XSModel::XSModel(CLASSLogger* log): CLASSObject(log)
+{
+
+}
\ No newline at end of file
--
GitLab