diff --git a/source/branches/CLASSV3/include/CLASSObject.hxx b/source/branches/CLASSV3/include/CLASSObject.hxx index 9f8217c57def5d57cb1584e433e58b440131bb0c..9009bc686d79a4dfc2e85445308474e68d937475 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 97fb8b8be58ea449ff7cc7443caaac4c90bd5106..5b4dc8c15c8114cbab20fb9f5f0ad72e149e0014 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 6df8a9b2fcc574dfeb286aa4f2f1553131946c0f..19def55433f0ad7ef67e39c57cb23c35e233c227 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 bf110feaca4a9256f5e6a920e7b229816bb8a8e9..66ed554be86c3df43bf5c43127e1c0a501f09f81 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 2fbface2cd38878ce1e199ea4f89e231928e9d5d..b7b6ef1968dc39a3f81ab651c011fd251132872b 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 64172f423ff2d46799ba7ba656eb5e17341aa794..359eed08044b5496158b96907fae1134af521298 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 0d6dbaca0e8d9c9db144b4c95003ed942fe30be4..0000000000000000000000000000000000000000 --- 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 b604ff50fab2a5052b1eab9a837ac760d2164aa1..0000000000000000000000000000000000000000 --- 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 79ce73946849e91d319b45a7353bea0f2ab10e72..0000000000000000000000000000000000000000 --- 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 d4d77f63d0c2f3fdaf13f07d7059c6c448e1c111..fb646da2768a6fabaa06d21668746f1dc8e7c6bf 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 c9faec54d390afa835449ee421d0a62b596431af..3f45095cc5365f3a5c5b385a51969f230459539c 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 2ed21584203c8a9740cb3278e8bae9d3bb1865ba..5cf7aaf4acdac6c82988eb3ef9725317cc6c8c14 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 4e6c13fd9ef2800f951c364834c5fb60c8a89ebe..0000000000000000000000000000000000000000 --- 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 79fe063ff4290ac738f220192acc5f55e0ccea49..0000000000000000000000000000000000000000 --- 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 a727720c13080c5c598a4a753520aaa99dc693e1..c7f829a6a9b3de3a1feea639b12b64a25084af6d 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 2823ca5ea407c9616a044075c92d5bb5c19f9457..96968fd227e0a6de94a7959ae3af02512b9358da 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 1299771e18f038800be21826446ba0f414f444a4..45859125183920c12901d33fdce9bcffd4201760 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 647234ce72e2ab2d8c30d4f83848cea397f095cd..696d8d9808bffae252ae9614e607f6f58254ba14 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 e8da2132244e4f1b50d6e3d8e2d5724f20495737..715c947c687280c3f7a58ad94241ccb3b3d731ed 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 93ce0a5df9dd7bdc02272ea4250425a2cf799731..dcfda34f87ddaa315ead90d308c4a8036d1066da 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 6b65d0fae1e78d67319e102366fcbfe225a02638..0000000000000000000000000000000000000000 --- 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 b61a17f6cf5fc2a1728ca07c425bb69bcd0d9507..0000000000000000000000000000000000000000 --- 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 19b7c12ce271f49a9c763e54cc523002e4a30bff..1299790c022e071470cf4284bd04e29d6f9b14f2 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