From 286f39d165486ee503f527f01ce0b9aa283023a8 Mon Sep 17 00:00:00 2001 From: Baptiste Mouginot <mouginot.baptiste@gmail.com> Date: Wed, 9 Jul 2014 12:13:33 +0000 Subject: [PATCH] git-svn-id: svn+ssh://svn.in2p3.fr/class@320 0e7d625b-0364-4367-a6be-d5be4a48d228 --- .../branches/CLASSV3/include/CLASSObject.hxx | 6 +- .../CLASSV3/include/EvolutionData.hxx | 4 +- .../branches/CLASSV3/include/PhysicModels.hxx | 1 + source/branches/CLASSV3/src/CLASSFacility.cxx | 2 +- source/branches/CLASSV3/src/CLASSObject.cxx | 5 +- source/branches/CLASSV3/src/DecayDataBank.cxx | 11 +- .../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 ----- .../branches/CLASSV3/src/EquivalenceModel.cxx | 2 +- source/branches/CLASSV3/src/EvolutionData.cxx | 8 +- .../branches/CLASSV3/src/FabricationPlant.cxx | 2 +- source/branches/CLASSV3/src/IM_Matrix.cxx | 315 ----------- source/branches/CLASSV3/src/IM_RK4.cxx | 408 -------------- .../branches/CLASSV3/src/IrradiationModel.cxx | 2 +- source/branches/CLASSV3/src/PhysicModels.cxx | 13 +- source/branches/CLASSV3/src/Pool.cxx | 2 +- source/branches/CLASSV3/src/Reactor.cxx | 4 +- source/branches/CLASSV3/src/Scenario.cxx | 2 +- source/branches/CLASSV3/src/Storage.cxx | 2 +- source/branches/CLASSV3/src/XSM_CLOSEST.cxx | 372 ------------- .../branches/CLASSV3/src/XSM_MLP_PWR_MOX.cxx | 509 ------------------ source/branches/CLASSV3/src/XSModel.cxx | 4 +- 23 files changed, 35 insertions(+), 2216 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_PWR_MOX.cxx diff --git a/source/branches/CLASSV3/include/CLASSObject.hxx b/source/branches/CLASSV3/include/CLASSObject.hxx index 9f8217c57..9009bc686 100644 --- a/source/branches/CLASSV3/include/CLASSObject.hxx +++ b/source/branches/CLASSV3/include/CLASSObject.hxx @@ -45,10 +45,9 @@ public : virtual CLASSObject* Clone() { return new CLASSObject(*this); } //!< Correct way to copy a CLASSObject in case of derivation - void SetLog(CLASSLogger* log) { fLog = log; fIsLog = true; } //!< Set the CLASSLogger + void SetLog(CLASSLogger* log) { fLog = log;} //!< Set the CLASSLogger CLASSLogger* GetLog() { return fLog; } //!< Return the Pointer to the Log - bool IsLog() { return fIsLog; } //!< reutrn true if a CLASSLogger is defined using TNamed::SetName; protected : @@ -56,8 +55,7 @@ protected : private : - bool fIsLog; //!< Set at true if a CLASSLogger are define - + ClassDef(CLASSObject,0); }; diff --git a/source/branches/CLASSV3/include/EvolutionData.hxx b/source/branches/CLASSV3/include/EvolutionData.hxx index 97fb8b8be..5b4dc8c15 100755 --- a/source/branches/CLASSV3/include/EvolutionData.hxx +++ b/source/branches/CLASSV3/include/EvolutionData.hxx @@ -70,7 +70,7 @@ public : Use create an empty EvolutionData loading a CLASSLogger \param CLASSLogger CLASSLogger used for the log... */ - EvolutionData(CLASSLogger* Log); ///< Make a new Evolutive Product evolution + EvolutionData(CLASSLogger* log); ///< Make a new Evolutive Product evolution //} //{ @@ -82,7 +82,7 @@ public : \param oldread true if the oldmethod should be use to read the DatBase File \param zai set the ZAI if you want to add a stable nuclei. */ - EvolutionData(CLASSLogger* Log, string DB_file, bool oldread = true, ZAI zai = ZAI(0,0,0) ); + EvolutionData(CLASSLogger* log, string DB_file, bool oldread = true, ZAI zai = ZAI(0,0,0) ); //} diff --git a/source/branches/CLASSV3/include/PhysicModels.hxx b/source/branches/CLASSV3/include/PhysicModels.hxx index 6df8a9b2f..19def5543 100644 --- a/source/branches/CLASSV3/include/PhysicModels.hxx +++ b/source/branches/CLASSV3/include/PhysicModels.hxx @@ -53,6 +53,7 @@ class PhysicModels : public CLASSObject //@{ PhysicModels(XSModel* XS, EquivalenceModel* EM, IrradiationModel* IM ); + PhysicModels(CLASSLogger* log, XSModel* XS, EquivalenceModel* EM, IrradiationModel* IM ); ~PhysicModels() {;} //{ diff --git a/source/branches/CLASSV3/src/CLASSFacility.cxx b/source/branches/CLASSV3/src/CLASSFacility.cxx index bf110feac..66ed554be 100644 --- a/source/branches/CLASSV3/src/CLASSFacility.cxx +++ b/source/branches/CLASSV3/src/CLASSFacility.cxx @@ -17,7 +17,7 @@ ClassImp(CLASSFacility) -CLASSFacility::CLASSFacility(int type):CLASSObject() +CLASSFacility::CLASSFacility(int type):CLASSObject(new CLASSLogger("CLASSFacility.log")) { fParc = 0; diff --git a/source/branches/CLASSV3/src/CLASSObject.cxx b/source/branches/CLASSV3/src/CLASSObject.cxx index 2fbface2c..b7b6ef196 100644 --- a/source/branches/CLASSV3/src/CLASSObject.cxx +++ b/source/branches/CLASSV3/src/CLASSObject.cxx @@ -17,11 +17,10 @@ ClassImp(CLASSObject) CLASSObject::CLASSObject() { - fLog = 0; - fIsLog = false; + fLog = new CLASSLogger("CLASSObject.log"); } CLASSObject::CLASSObject(CLASSLogger* log) { - SetLog(log); + fLog = log; } \ No newline at end of file diff --git a/source/branches/CLASSV3/src/DecayDataBank.cxx b/source/branches/CLASSV3/src/DecayDataBank.cxx index 64172f423..359eed080 100644 --- a/source/branches/CLASSV3/src/DecayDataBank.cxx +++ b/source/branches/CLASSV3/src/DecayDataBank.cxx @@ -25,7 +25,7 @@ // //________________________________________________________________________ -DecayDataBank::DecayDataBank():CLASSObject() +DecayDataBank::DecayDataBank():CLASSObject(new CLASSLogger("DecayDataBank.log")) { } @@ -35,21 +35,16 @@ DecayDataBank::DecayDataBank():CLASSObject() //________________________________________________________________________ //________________________________________________________________________ -DecayDataBank::DecayDataBank(CLASSLogger* Log, string DB_index_file, bool setlog, bool olfreadmethod) +DecayDataBank::DecayDataBank(CLASSLogger* log, string DB_index_file, bool setlog, bool olfreadmethod):CLASSObject(log) { - SetLog(Log); fDataBaseIndex = DB_index_file; - fOldReadMethod = olfreadmethod; // Warning - if(IsLog()) - { INFO << " A EvolutionData has been define :" << endl; INFO << "\t His index is : \"" << DB_index_file << "\"" << endl << endl; - } - + } //________________________________________________________________________ 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/EquivalenceModel.cxx b/source/branches/CLASSV3/src/EquivalenceModel.cxx index d4d77f63d..fb646da27 100644 --- a/source/branches/CLASSV3/src/EquivalenceModel.cxx +++ b/source/branches/CLASSV3/src/EquivalenceModel.cxx @@ -5,7 +5,7 @@ -EquivalenceModel::EquivalenceModel():CLASSObject() +EquivalenceModel::EquivalenceModel():CLASSObject(new CLASSLogger("EquivalenceModel.log")) { } diff --git a/source/branches/CLASSV3/src/EvolutionData.cxx b/source/branches/CLASSV3/src/EvolutionData.cxx index c9faec54d..3f45095cc 100755 --- a/source/branches/CLASSV3/src/EvolutionData.cxx +++ b/source/branches/CLASSV3/src/EvolutionData.cxx @@ -343,7 +343,7 @@ EvolutionData Sum(EvolutionData const& evol1, EvolutionData const& evol2) ClassImp(EvolutionData) -EvolutionData::EvolutionData():CLASSObject() +EvolutionData::EvolutionData():CLASSObject(new CLASSLogger("EvolutionData.log")) { fIsCrossSection = false; fPower = 0; @@ -353,10 +353,9 @@ EvolutionData::EvolutionData():CLASSObject() } //________________________________________________________________________ -EvolutionData::EvolutionData(CLASSLogger* Log) +EvolutionData::EvolutionData(CLASSLogger* log):CLASSObject(log) { - SetLog(Log); fIsCrossSection = false; fPower = 0; fCycleTime = 0; @@ -367,10 +366,9 @@ EvolutionData::EvolutionData(CLASSLogger* Log) } //________________________________________________________________________ -EvolutionData::EvolutionData(CLASSLogger* Log, string DB_file, bool oldread, ZAI zai) +EvolutionData::EvolutionData(CLASSLogger* log, string DB_file, bool oldread, ZAI zai):CLASSObject(log) { - SetLog(Log); fIsCrossSection = false; fDB_file = DB_file; fPower = 0; diff --git a/source/branches/CLASSV3/src/FabricationPlant.cxx b/source/branches/CLASSV3/src/FabricationPlant.cxx index 2ed215842..5cf7aaf4a 100644 --- a/source/branches/CLASSV3/src/FabricationPlant.cxx +++ b/source/branches/CLASSV3/src/FabricationPlant.cxx @@ -36,7 +36,7 @@ ClassImp(FabricationPlant) -FabricationPlant::FabricationPlant():CLASSFacility(16) +FabricationPlant::FabricationPlant():CLASSFacility(new CLASSLogger("FabricationPlant.log"), 16) { SetName("F_FabricationPLant."); diff --git a/source/branches/CLASSV3/src/IM_Matrix.cxx b/source/branches/CLASSV3/src/IM_Matrix.cxx deleted file mode 100644 index 4e6c13fd9..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():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 79fe063ff..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(), 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/IrradiationModel.cxx b/source/branches/CLASSV3/src/IrradiationModel.cxx index a727720c1..c7f829a6a 100644 --- a/source/branches/CLASSV3/src/IrradiationModel.cxx +++ b/source/branches/CLASSV3/src/IrradiationModel.cxx @@ -26,7 +26,7 @@ using namespace std; -IrradiationModel::IrradiationModel():CLASSObject() +IrradiationModel::IrradiationModel():CLASSObject(new CLASSLogger("IrradiationModel.log")) { fShorstestHalflife = 3600.*24*2.; fZAIThreshold = 90; diff --git a/source/branches/CLASSV3/src/PhysicModels.cxx b/source/branches/CLASSV3/src/PhysicModels.cxx index 2823ca5ea..96968fd22 100644 --- a/source/branches/CLASSV3/src/PhysicModels.cxx +++ b/source/branches/CLASSV3/src/PhysicModels.cxx @@ -8,7 +8,7 @@ // // //________________________________________________________________________ -PhysicModels::PhysicModels(XSModel* XS, EquivalenceModel* EM, IrradiationModel* IM ) +PhysicModels::PhysicModels(XSModel* XS, EquivalenceModel* EM, IrradiationModel* IM ):CLASSObject(new CLASSLogger("PhysicsModel.log")) { fXSModel = XS; @@ -17,6 +17,17 @@ PhysicModels::PhysicModels(XSModel* XS, EquivalenceModel* EM, IrradiationModel* } +//________________________________________________________________________ +PhysicModels::PhysicModels(CLASSLogger* log, XSModel* XS, EquivalenceModel* EM, IrradiationModel* IM ):CLASSObject(log) +{ + + fXSModel = XS; + fEquivalenceModel = EM; + fIrradiationModel = IM; + + +} + //________________________________________________________________________ EvolutionData PhysicModels::GenerateEvolutionData(IsotopicVector IV, double cycletime, double Power) { diff --git a/source/branches/CLASSV3/src/Pool.cxx b/source/branches/CLASSV3/src/Pool.cxx index 1299771e1..458591251 100755 --- a/source/branches/CLASSV3/src/Pool.cxx +++ b/source/branches/CLASSV3/src/Pool.cxx @@ -22,7 +22,7 @@ ClassImp(Pool) -Pool::Pool():CLASSBackEnd(8) +Pool::Pool():CLASSBackEnd(new CLASSLogger("Pool.log"), 8) { fOutBackEndFacility = 0; SetName("P_Pool."); diff --git a/source/branches/CLASSV3/src/Reactor.cxx b/source/branches/CLASSV3/src/Reactor.cxx index 647234ce7..696d8d980 100755 --- a/source/branches/CLASSV3/src/Reactor.cxx +++ b/source/branches/CLASSV3/src/Reactor.cxx @@ -27,13 +27,11 @@ ClassImp(Reactor) -Reactor::Reactor():CLASSFacility() +Reactor::Reactor():CLASSFacility(new CLASSLogger("Reactor.log"), 4) { - SetFacilityType(4); SetName("R_Reactor."); - fOutBackEndFacility = 0; fStorage = 0; fFuelTypeDB = 0; diff --git a/source/branches/CLASSV3/src/Scenario.cxx b/source/branches/CLASSV3/src/Scenario.cxx index e8da21322..715c947c6 100755 --- a/source/branches/CLASSV3/src/Scenario.cxx +++ b/source/branches/CLASSV3/src/Scenario.cxx @@ -44,7 +44,7 @@ string dtoa(double num) } //________________________________________________________________________ -Scenario::Scenario():CLASSObject(new CLASSLogger()) +Scenario::Scenario():CLASSObject(new CLASSLogger("CLASS_OUTPUT.log")) { fNewTtree = true; diff --git a/source/branches/CLASSV3/src/Storage.cxx b/source/branches/CLASSV3/src/Storage.cxx index 93ce0a5df..dcfda34f8 100644 --- a/source/branches/CLASSV3/src/Storage.cxx +++ b/source/branches/CLASSV3/src/Storage.cxx @@ -20,7 +20,7 @@ //________________________________________________________________________ ClassImp(Storage) -Storage::Storage():CLASSBackEnd(-1) +Storage::Storage():CLASSBackEnd(new CLASSLogger("Storage.log"), -1) { SetIsStorageType(); SetName("S_Storage."); diff --git a/source/branches/CLASSV3/src/XSM_CLOSEST.cxx b/source/branches/CLASSV3/src/XSM_CLOSEST.cxx deleted file mode 100644 index 6b65d0fae..000000000 --- a/source/branches/CLASSV3/src/XSM_CLOSEST.cxx +++ /dev/null @@ -1,372 +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(); - - if(IsLog()) - { - // 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(); - - if(IsLog()) - { - // 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_PWR_MOX.cxx b/source/branches/CLASSV3/src/XSM_MLP_PWR_MOX.cxx deleted file mode 100644 index b61a17f6c..000000000 --- a/source/branches/CLASSV3/src/XSM_MLP_PWR_MOX.cxx +++ /dev/null @@ -1,509 +0,0 @@ - -#include "XSModel.hxx" -#include "XSM_MLP_PWR_MOX.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_PWR_MOX -// -// -// -// -//________________________________________________________________________ -XSM_MLP_PWR_MOX::XSM_MLP_PWR_MOX(string TMVA_Weight_Directory,string InformationFile, bool IsTimeStep):XSModel(new CLASSLogger("XSM_MLP_PWR_MOX.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_PWR_MOX::XSM_MLP_PWR_MOX(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_PWR_MOX::~XSM_MLP_PWR_MOX() -{ - fMapOfTMVAVariableNames.clear(); -} -//________________________________________________________________________ -void XSM_MLP_PWR_MOX::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_PWR_MOX::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_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() ; - 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_PWR_MOX::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_PWR_MOX::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_PWR_MOX::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_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() ); - - if(Z==-1 || A==-1 || I==-1 || Reaction==-1 || TimeStep==-1) - { - ERROR << " wrong TMVA weight format " << endl; - exit(0); - } -} - -//________________________________________________________________________ -EvolutionData XSM_MLP_PWR_MOX::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_PWR_MOX::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 19b7c12ce..1299790c0 100644 --- a/source/branches/CLASSV3/src/XSModel.cxx +++ b/source/branches/CLASSV3/src/XSModel.cxx @@ -10,14 +10,14 @@ using namespace std; -XSModel::XSModel(): CLASSObject() +XSModel::XSModel():CLASSObject(new CLASSLogger("XSModel.log")) { } -XSModel::XSModel(CLASSLogger* log): CLASSObject(log) +XSModel::XSModel(CLASSLogger* log):CLASSObject(log) { } \ No newline at end of file -- GitLab