From 05226e12450925cf030d2fa4492f34061d56b7cc Mon Sep 17 00:00:00 2001
From: Baptiste Mouginot <mouginot.baptiste@gmail.com>
Date: Thu, 22 Jan 2015 09:04:13 +0000
Subject: [PATCH] updating source: step 2a copying the new source : src
git-svn-id: svn+ssh://svn.in2p3.fr/class@478 0e7d625b-0364-4367-a6be-d5be4a48d228
---
source/trunk/src/CLASSBackEnd.cxx | 126 +++
source/trunk/src/CLASSFacility.cxx | 76 ++
source/trunk/src/CLASSFuel.cxx | 26 +
source/trunk/src/CLASSFuelPlan.cxx | 69 ++
source/trunk/src/CLASSLogger.cxx | 201 ++++
source/trunk/src/CLASSNucleiFiliation.cxx | 328 ++++++
source/trunk/src/CLASSObject.cxx | 24 +
source/trunk/src/DecayDataBank.cxx | 172 +++
source/trunk/src/DynamicalSystem.cxx | 214 ++++
source/trunk/src/EquivalenceModel.cxx | 232 +++++
source/trunk/src/EvolutionData.cxx | 1160 +++++++++++++++++++++
source/trunk/src/FabricationPlant.cxx | 626 +++++++++++
source/trunk/src/IrradiationModel.cxx | 638 ++++++++++++
source/trunk/src/IsotopicVector.cxx | 681 ++++++++++++
source/trunk/src/Makefile | 125 +++
source/trunk/src/PhysicsModels.cxx | 57 +
source/trunk/src/Pool.cxx | 226 ++++
source/trunk/src/Reactor.cxx | 607 +++++++++++
source/trunk/src/Scenario.cxx | 888 ++++++++++++++++
source/trunk/src/SeparationPlant.cxx | 160 +++
source/trunk/src/Storage.cxx | 215 ++++
source/trunk/src/XSModel.cxx | 61 ++
source/trunk/src/ZAI.cxx | 62 ++
source/trunk/src/ZAIMass.cxx | 82 ++
24 files changed, 7056 insertions(+)
create mode 100644 source/trunk/src/CLASSBackEnd.cxx
create mode 100644 source/trunk/src/CLASSFacility.cxx
create mode 100644 source/trunk/src/CLASSFuel.cxx
create mode 100644 source/trunk/src/CLASSFuelPlan.cxx
create mode 100755 source/trunk/src/CLASSLogger.cxx
create mode 100644 source/trunk/src/CLASSNucleiFiliation.cxx
create mode 100644 source/trunk/src/CLASSObject.cxx
create mode 100644 source/trunk/src/DecayDataBank.cxx
create mode 100644 source/trunk/src/DynamicalSystem.cxx
create mode 100644 source/trunk/src/EquivalenceModel.cxx
create mode 100755 source/trunk/src/EvolutionData.cxx
create mode 100644 source/trunk/src/FabricationPlant.cxx
create mode 100644 source/trunk/src/IrradiationModel.cxx
create mode 100755 source/trunk/src/IsotopicVector.cxx
create mode 100755 source/trunk/src/Makefile
create mode 100644 source/trunk/src/PhysicsModels.cxx
create mode 100755 source/trunk/src/Pool.cxx
create mode 100755 source/trunk/src/Reactor.cxx
create mode 100755 source/trunk/src/Scenario.cxx
create mode 100644 source/trunk/src/SeparationPlant.cxx
create mode 100644 source/trunk/src/Storage.cxx
create mode 100644 source/trunk/src/XSModel.cxx
create mode 100755 source/trunk/src/ZAI.cxx
create mode 100644 source/trunk/src/ZAIMass.cxx
diff --git a/source/trunk/src/CLASSBackEnd.cxx b/source/trunk/src/CLASSBackEnd.cxx
new file mode 100644
index 000000000..faf588e85
--- /dev/null
+++ b/source/trunk/src/CLASSBackEnd.cxx
@@ -0,0 +1,126 @@
+#include "CLASSBackEnd.hxx"
+
+#include "DecayDataBank.hxx"
+#include "Scenario.hxx"
+#include "CLASSLogger.hxx"
+
+
+#include <sstream>
+#include <string>
+#include <iostream>
+#include <cmath>
+#include <algorithm>
+
+//________________________________________________________________________
+//
+// CLASSBackEnd
+//
+//
+//
+//
+//________________________________________________________________________
+ClassImp(CLASSBackEnd)
+
+CLASSBackEnd::CLASSBackEnd(int type):CLASSFacility(type)
+{
+ fDecayDataBase = 0;
+}
+
+CLASSBackEnd::CLASSBackEnd(CLASSLogger* log, int type):CLASSFacility(log, type)
+{
+ fDecayDataBase = 0;
+}
+
+
+CLASSBackEnd::CLASSBackEnd(CLASSLogger* log, cSecond cycletime, int type):CLASSFacility(log, cycletime, type)
+{
+ fDecayDataBase = 0;
+}
+
+//________________________________________________________________________
+void CLASSBackEnd::ClearIVArray()
+{
+
+ IsotopicVector EmptyIV;
+ fInsideIV = EmptyIV;
+ fIVArray.clear();
+ fIVArrayArrivalTime.clear();
+}
+
+//________________________________________________________________________
+void CLASSBackEnd::AddIV(IsotopicVector isotopicvector)
+{
+
+ AddCumulativeIVIn(isotopicvector);
+
+ fIVArray.push_back(isotopicvector);
+ fIVArrayArrivalTime.push_back(fInternalTime);
+
+
+}
+//________________________________________________________________________
+void CLASSBackEnd::UpdateInsideIV()
+{
+ DBGL
+ fInsideIV = IsotopicVector();
+ for(int i = 0; i < (int)fIVArray.size(); i++)
+ fInsideIV += fIVArray[i];
+ DBGL
+}
+
+
+//________________________________________________________________________
+// Get Decay
+//________________________________________________________________________
+IsotopicVector CLASSBackEnd::GetDecay(IsotopicVector isotopicvector, cSecond t)
+{
+ DBGL
+
+ IsotopicVector IV;
+
+ map<ZAI ,double> isotopicquantity = isotopicvector.GetIsotopicQuantity();
+ map<ZAI ,double >::iterator it;
+ for( it = isotopicquantity.begin(); it != isotopicquantity.end(); it++)
+ {
+ if((*it).second > 0)
+ {
+ IsotopicVector ivtmp = fDecayDataBase->Evolution(it->first, t) * (*it).second ;
+ IV += ivtmp;
+ }
+ }
+
+ DBGL
+ return IV;
+}
+
+
+map<cSecond,int> CLASSBackEnd::GetTheBackEndTimePath()
+{
+ DBGL
+ map<cSecond, int> TheBackEndTimePath;
+
+ if(!fIsStorageType)
+ {
+ int FacilityType = GetFacilityType();
+ cSecond step = GetCycleTime();
+
+ pair< map<cSecond, int>::iterator, bool > IResult = TheBackEndTimePath.insert(pair<cSecond,int> (step, FacilityType));
+ if( !IResult.second ) IResult.first->second |= FacilityType;
+
+ map<cSecond, int> TheBackEndTimePath_tmp = GetOutBackEndFacility()->GetTheBackEndTimePath();
+
+ map<cSecond, int>::iterator it;
+ for (it = TheBackEndTimePath_tmp.begin(); it != TheBackEndTimePath_tmp.end(); it++)
+ {
+ pair< map<cSecond, int>::iterator, bool > IResult;
+
+ IResult = TheBackEndTimePath.insert( pair<cSecond ,int>(step + (*it).first, (*it).second) );
+ if( !IResult.second )
+ IResult.first->second |= (*it).second;
+
+ }
+ }
+
+ DBGL
+ return TheBackEndTimePath;
+}
\ No newline at end of file
diff --git a/source/trunk/src/CLASSFacility.cxx b/source/trunk/src/CLASSFacility.cxx
new file mode 100644
index 000000000..9759b0d3f
--- /dev/null
+++ b/source/trunk/src/CLASSFacility.cxx
@@ -0,0 +1,76 @@
+#include "CLASSFacility.hxx"
+#include "Scenario.hxx"
+
+using namespace std;
+
+ //________________________________________________________________________
+ //
+ // CLASSFacility
+ //
+ //
+ //
+ //
+ //________________________________________________________________________
+
+ClassImp(CLASSFacility)
+
+
+
+
+CLASSFacility::CLASSFacility(int type):CLASSObject()
+{
+ fParc = 0;
+ fFacilityType = type;
+ fInternalTime = 0;
+ fInCycleTime = 0;
+ fCycleTime = -1;
+}
+
+CLASSFacility::CLASSFacility(CLASSLogger* log, int type):CLASSObject(log)
+{
+ fParc = 0;
+
+ fFacilityType = type;
+
+ fInternalTime = 0;
+ fInCycleTime = 0;
+ fCycleTime = -1;
+}
+
+
+CLASSFacility::CLASSFacility(CLASSLogger* log, cSecond cycletime, int type):CLASSObject(log)
+{
+ fParc = 0;
+ fFacilityType = type;
+ fInternalTime = 0;
+ fInCycleTime = 0;
+ fCycleTime = cycletime;
+
+}
+
+CLASSFacility::CLASSFacility(CLASSLogger* log, cSecond creationtime, cSecond lifetime, int type):CLASSObject(log)
+{
+ fParc = 0;
+
+ fFacilityType = type;
+
+ fInternalTime = 0;
+ fInCycleTime = 0;
+ fCycleTime = -1;
+ fCreationTime = creationtime;
+ fLifeTime = lifetime;
+}
+
+CLASSFacility::CLASSFacility(CLASSLogger* log, cSecond creationtime, cSecond lifetime, cSecond cycletime, int type):CLASSObject(log)
+{
+ fParc = 0;
+
+ fFacilityType = type;
+
+ fInternalTime = 0;
+ fInCycleTime = 0;
+ fCycleTime = cycletime;
+ fCreationTime = creationtime;
+ fInternalTime = fCreationTime;
+ fLifeTime = lifetime;
+}
\ No newline at end of file
diff --git a/source/trunk/src/CLASSFuel.cxx b/source/trunk/src/CLASSFuel.cxx
new file mode 100644
index 000000000..de1d23c3a
--- /dev/null
+++ b/source/trunk/src/CLASSFuel.cxx
@@ -0,0 +1,26 @@
+#include "CLASSFuel.hxx"
+
+#include "CLASSLogger.hxx"
+
+using namespace std;
+
+//________________________________________________________________________
+//
+// CLASSFuel
+//
+//
+//
+//
+//________________________________________________________________________
+
+CLASSFuel::CLASSFuel(EvolutionData* evo)
+{
+ fEvolutionData = evo;
+ fPhysicsModels = 0;
+}
+
+CLASSFuel::CLASSFuel(PhysicsModels* evo)
+{
+ fEvolutionData = 0;
+ fPhysicsModels = evo;
+}
\ No newline at end of file
diff --git a/source/trunk/src/CLASSFuelPlan.cxx b/source/trunk/src/CLASSFuelPlan.cxx
new file mode 100644
index 000000000..a57a24d78
--- /dev/null
+++ b/source/trunk/src/CLASSFuelPlan.cxx
@@ -0,0 +1,69 @@
+#include "CLASSFuelPlan.hxx"
+
+#include "CLASSLogger.hxx"
+
+using namespace std;
+
+//________________________________________________________________________
+//
+// CLASSFuelPlan
+//
+//
+//
+//
+//________________________________________________________________________
+
+CLASSFuelPlan::CLASSFuelPlan():CLASSObject()
+{
+ DBGL
+}
+
+CLASSFuelPlan::CLASSFuelPlan(CLASSLogger* log):CLASSObject(log)
+{
+ DBGL
+ DBGL
+}
+
+void CLASSFuelPlan::AddFuel(cSecond time, CLASSFuel fuel, double BurnUp)
+{
+ DBGL
+ fLoadingPlan.insert(pair< cSecond, pair< CLASSFuel, double > >(time, pair< CLASSFuel, double > (fuel, BurnUp)) );
+ DBGL
+}
+
+pair< CLASSFuel, double > CLASSFuelPlan::GetFuelAt(cSecond t)
+{
+ DBGL
+
+ pair< CLASSFuel, double > FuelAtT = fLoadingPlan.begin()->second;
+
+ map< cSecond, pair< CLASSFuel, double > >::iterator it;
+
+ cSecond ThisFuelTime;
+
+ for (it = fLoadingPlan.begin(); it != fLoadingPlan.end(); it++ )
+ {
+ if( it == fLoadingPlan.begin())
+ {
+ FuelAtT = (*it).second;
+ }
+ else
+ {
+ ThisFuelTime = (*it).first;
+
+ if( t < ThisFuelTime )
+ {
+ DBGL
+ return FuelAtT;
+ }
+ else
+ {
+ FuelAtT = (*it).second;
+ }
+
+ }
+ }
+
+ DBGL
+ return FuelAtT;
+}
diff --git a/source/trunk/src/CLASSLogger.cxx b/source/trunk/src/CLASSLogger.cxx
new file mode 100755
index 000000000..d5b95630f
--- /dev/null
+++ b/source/trunk/src/CLASSLogger.cxx
@@ -0,0 +1,201 @@
+
+#include "CLASSLogger.hxx"
+
+#include <iostream>
+#include <ostream>
+#include <algorithm>
+using namespace std;
+
+//________________________________________________________________________
+//
+// CLASSLogger
+//
+//
+//
+//
+//________________________________________________________________________
+CLASSLogger::CLASSLogger()
+{
+ string CLASSLoggerName = "CLASS_OUTPUT.log";
+ int VerboseLvl = 0;
+ int OutputLvl = 1;
+
+
+ fCLASSLoggerName = CLASSLoggerName;
+ fOutPutFile.open(CLASSLoggerName.c_str());
+ fVerboseLVL = VerboseLvl;
+ if(!fOutPutFile)
+ {
+ cout << "Could not open the CLASSLogger: " << CLASSLoggerName << " !" << endl;
+ exit(-1);
+ }
+
+ fError = 0;
+ fWarning = 0;
+ fDebug = 0;
+ fInfo = 0;
+
+ if (VerboseLvl <= OutputLvl)
+ fMaxOutPutLVL = OutputLvl;
+ else
+ fMaxOutPutLVL = VerboseLvl;
+
+ if(VerboseLvl >= 0)
+ {
+ if (!fError)
+ fError = new LogType(std::cout);
+ else
+ fError->SetSecondOutput(std::cout);
+ }
+ if(VerboseLvl >= 1)
+ {
+ if (!fWarning)
+ fWarning = new LogType(std::cout);
+ else
+ fWarning->SetSecondOutput(std::cout);
+ }
+ if(VerboseLvl >= 2)
+ {
+ if (!fInfo)
+ fInfo = new LogType(std::cout);
+ else
+ fInfo->SetSecondOutput(std::cout);
+ }
+ if(VerboseLvl >= 3)
+ {
+ if (!fDebug)
+ fDebug = new LogType(std::cout);
+ else
+ fDebug->SetSecondOutput(std::cout);
+ }
+ if(OutputLvl >= 0)
+ {
+ if (!fError)
+ fError = new LogType(fOutPutFile);
+ else
+ fError->SetSecondOutput(fOutPutFile);
+ }
+ if(OutputLvl >= 1)
+ {
+ if (!fWarning)
+ fWarning = new LogType(fOutPutFile);
+ else
+ fWarning->SetSecondOutput(fOutPutFile);
+ }
+ if(OutputLvl >= 2)
+ {
+ if (!fInfo)
+ fInfo = new LogType(fOutPutFile);
+ else
+ fInfo->SetSecondOutput(fOutPutFile);
+ }
+ if(OutputLvl >= 3)
+ {
+ if (!fDebug)
+ fDebug = new LogType(fOutPutFile);
+ else
+ fDebug->SetSecondOutput(fOutPutFile);
+ }
+
+
+
+}
+
+
+
+CLASSLogger::CLASSLogger(string CLASSLoggerName, int VerboseLvl, int OutputLvl )
+{
+ fCLASSLoggerName = CLASSLoggerName;
+ fOutPutFile.open(CLASSLoggerName.c_str());
+ fVerboseLVL = VerboseLvl;
+ if(!fOutPutFile)
+ {
+ cout << "Could not open the CLASSLogger: " << CLASSLoggerName << " !" << endl;
+ exit(-1);
+ }
+
+ fError = 0;
+ fWarning = 0;
+ fDebug = 0;
+ fInfo = 0;
+
+ if (VerboseLvl <= OutputLvl)
+ fMaxOutPutLVL = OutputLvl;
+ else
+ fMaxOutPutLVL = VerboseLvl;
+
+ if(VerboseLvl >= 0)
+ {
+ if (!fError)
+ fError = new LogType(std::cout);
+ else
+ fError->SetSecondOutput(std::cout);
+ }
+
+ if(VerboseLvl >= 1)
+ {
+ if (!fWarning)
+ fWarning = new LogType(std::cout);
+ else
+ fWarning->SetSecondOutput(std::cout);
+ }
+ if(VerboseLvl >= 2)
+ {
+ if (!fInfo)
+ fInfo = new LogType(std::cout);
+ else
+ fInfo->SetSecondOutput(std::cout);
+ }
+ if(VerboseLvl >= 3)
+ {
+ if (!fDebug)
+ fDebug = new LogType(std::cout);
+ else
+ fDebug->SetSecondOutput(std::cout);
+ }
+
+
+ if(OutputLvl >= 0)
+ {
+ if (!fError)
+ fError = new LogType(fOutPutFile);
+ else
+ fError->SetSecondOutput(fOutPutFile);
+ }
+ if(OutputLvl >= 1)
+ {
+ if (!fWarning)
+ fWarning = new LogType(fOutPutFile);
+ else
+ fWarning->SetSecondOutput(fOutPutFile);
+ }
+ if(OutputLvl >= 2)
+ {
+ if (!fInfo)
+ fInfo = new LogType(fOutPutFile);
+ else
+ fInfo->SetSecondOutput(fOutPutFile);
+ }
+ if(OutputLvl >= 3)
+ {
+ if (!fDebug)
+ fDebug = new LogType(fOutPutFile);
+ else
+ fDebug->SetSecondOutput(fOutPutFile);
+ }
+
+
+
+}
+
+//________________________________________________________________________
+CLASSLogger::~CLASSLogger()
+{
+
+ fOutPutFile.close();
+
+}
+
+
+
+
diff --git a/source/trunk/src/CLASSNucleiFiliation.cxx b/source/trunk/src/CLASSNucleiFiliation.cxx
new file mode 100644
index 000000000..d267b37a9
--- /dev/null
+++ b/source/trunk/src/CLASSNucleiFiliation.cxx
@@ -0,0 +1,328 @@
+#include "CLASSNucleiFiliation.hxx"
+#include "ZAI.hxx"
+#include "IsotopicVector.hxx"
+
+#include <map>
+#include <vector>
+#include <cmath>
+
+using namespace std;
+
+//const string DEFAULTDATABASE = "DecayBase.dat";
+//________________________________________________________________________
+//
+// CLASSNucleiFiliation
+//
+//
+//
+//
+//________________________________________________________________________
+//____________________________InClass Operator____________________________
+//________________________________________________________________________
+
+CLASSNucleiFiliation::CLASSNucleiFiliation():CLASSObject()
+{
+}
+CLASSNucleiFiliation::CLASSNucleiFiliation(CLASSLogger* log):CLASSObject(log)
+{
+}
+
+CLASSNucleiFiliation::CLASSNucleiFiliation(const CLASSNucleiFiliation& CNF):CLASSObject()
+{
+ fNucleiFiliation = CNF.GetNucleiFIliation();
+}
+
+
+
+//________________________________________________________________________
+CLASSNucleiFiliation::~CLASSNucleiFiliation()
+{
+
+}
+
+
+//________________________________________________________________________
+void CLASSNucleiFiliation::Add( ZAI Mother, IsotopicVector Daughter )
+{
+ DBGL
+
+ pair< map<ZAI, IsotopicVector>::iterator, bool> IResult;
+
+
+ IResult = fNucleiFiliation.insert( pair<ZAI, IsotopicVector> ( Mother, Daughter ) );
+
+ if( !IResult.second)
+ (*IResult.first).second += Daughter ;
+
+
+ DBGL
+}
+
+
+//________________________________________________________________________
+IsotopicVector CLASSNucleiFiliation::GetFiliation(ZAI Mother)
+{
+ DBGV(Mother.Z()<<" "<<Mother.A()<<" "<<Mother.I());
+ map<ZAI, IsotopicVector>::iterator it_Filiation;
+
+ it_Filiation = fNucleiFiliation.find(Mother); // search for the ZAI in the map
+
+ DBGL
+ if(it_Filiation != fNucleiFiliation.end()) // test if it is present in the map
+ return it_Filiation->second;
+ else
+ return ZAI(-1,-1,-1)*1; // return -1 -1 _1 ZAI if unknown nuclei....
+
+}
+
+//________________________________________________________________________
+void CLASSNucleiFiliation::FiliationCleanUp(map<ZAI, int> GoodNuclei, CLASSNucleiFiliation CuttedNuclei)
+{
+ DBGL
+ map<ZAI, IsotopicVector>::iterator it_Filiation;
+ for(it_Filiation = fNucleiFiliation.begin(); it_Filiation != fNucleiFiliation.end(); it_Filiation++) // loop on the filiation map (on the Mother ZAI)
+ {
+ vector<ZAI> DautherList = it_Filiation->second.GetZAIList(); // recover for each mother ZAI, the list of its daughter ZAI
+
+ for (int i = 0; i < (int)DautherList.size(); i++) // Loop on the Daughter ZAI
+ {
+ if(GoodNuclei.find(DautherList[i]) == GoodNuclei.end() ) // if the ZAI is not in a dealed nuclei (cutted or unknown)
+ {
+ double Daughter_BR = it_Filiation->second.GetQuantity(DautherList[i]); // Get the quantity of the ZAI
+ it_Filiation->second -= Daughter_BR * DautherList[i]; // Remove it from the daughter list
+
+
+ IsotopicVector FastDecayChain = CuttedNuclei.GetFiliation(DautherList[i]); // Get the fast decay chain of it
+
+ if(FastDecayChain.GetQuantity(-1, -1, -1) == 0) // Check if the FastDecayChain is known
+ it_Filiation->second += Daughter_BR * FastDecayChain; // Add the FastDecayCHain & apply the BR for the cutted Daughter
+ else
+ {
+
+ ZAI Mother = DautherList[i];
+ while (FastDecayChain.GetQuantity(-1, -1, -1) != 0 && GoodNuclei.find(Mother) == GoodNuclei.end())
+ {
+ Mother = GetArtificialDecay(Mother); // Do an Artifial decay on the nuclei
+ FastDecayChain = CuttedNuclei.GetFiliation(Mother); // Get the fast decay chain of it
+ }
+
+ if(GoodNuclei.find(Mother) != GoodNuclei.end())
+ it_Filiation->second += Mother * Daughter_BR;
+
+ else if ( FastDecayChain.GetQuantity(-1, -1, -1) == 0)
+ it_Filiation->second += FastDecayChain * Daughter_BR;
+
+ else
+ {
+ ERROR << "Problem in Articial Decay!! check it!!" << endl;
+ exit(1);
+ }
+
+ }
+ }
+ }
+
+ }
+ DBGL
+}
+
+ZAI CLASSNucleiFiliation::GetArtificialDecay(ZAI Mother)
+{
+ DBGL
+
+ int A = Mother.A();
+ int Z = Mother.Z();
+ int I = Mother.I();
+
+ if(I!=0)
+ return ZAI(Z,A,I-1);
+ else
+ {
+ //Coef Ac & As of Bette & Weisacker are approximativ but enough precise for this application....
+ double Ac = 0.695;
+ double As = 23.2;
+
+ double ZTh = A/2 * ( 1 )/ ( 1 + Ac / (4*As) * pow(A,2/3) ); // Stable Z from isobarn calculation using Bette & Weisacker formula.
+
+
+ if( Z > ZTh ) // Then Beta+
+ return ZAI(Z-1,A,I);
+ else // Then Beta-
+ return ZAI(Z+1,A,I);
+
+ }
+
+
+ DBGL
+}
+
+
+//________________________________________________________________________
+void CLASSNucleiFiliation::SelfFiliationCleanUp(map<ZAI, int> GoodNuclei)
+{
+ DBGL
+
+ bool Cleaned = false;
+
+ while (!Cleaned) // loop until all the map is cleaned (all cut have been done)
+ {
+ int count = 0;
+ map<ZAI, IsotopicVector> CopyfNucleiFiliation = fNucleiFiliation;
+ map<ZAI, IsotopicVector>::iterator it_Filiation;
+ for(it_Filiation = fNucleiFiliation.begin(); it_Filiation != fNucleiFiliation.end(); it_Filiation++) // Loop on the mother ZAI
+ {
+ vector<ZAI> DautherList = it_Filiation->second.GetZAIList(); // Get the list of daughter ZAI
+
+ for (int i = 0; i < (int)DautherList.size(); i++) //Loop on daughter
+ {
+ if(GoodNuclei.find(DautherList[i]) == GoodNuclei.end() ) // if the ZAI is not in a dealed nuclei (cutted or unknown)
+ { count++;
+ map<ZAI, IsotopicVector>::iterator it_FiliationCopy = CopyfNucleiFiliation.find(it_Filiation->first) ;
+
+ double Daughter_BR = it_FiliationCopy->second.GetQuantity(DautherList[i]); // Get the quantity of the ZAI
+ it_FiliationCopy->second -= Daughter_BR * DautherList[i]; // Remove it from the daughter list
+ IsotopicVector FastDecayChain = (*this).GetFiliation(DautherList[i]); // Get the fast decay chain of it
+
+ if(FastDecayChain.GetQuantity(-1, -1, -1) == 0) // Check if the FastDecayChain is known
+ it_FiliationCopy->second += Daughter_BR * FastDecayChain; // Add the FastDecayCHain & apply the BR for the cutted Daughter
+ else
+ {
+ ZAI Mother = DautherList[i];
+ while (FastDecayChain.GetQuantity(-1, -1, -1) != 0 && GoodNuclei.find(Mother) == GoodNuclei.end())
+ {
+ Mother = GetArtificialDecay(Mother); // Do an Artifial decay on the nuclei
+ FastDecayChain = (*this).GetFiliation(Mother); // Get the fast decay chain of it
+ }
+
+ if(GoodNuclei.find(Mother) != GoodNuclei.end())
+ it_FiliationCopy->second += Mother * Daughter_BR;
+
+ else if ( FastDecayChain.GetQuantity(-1, -1, -1) == 0)
+ it_FiliationCopy->second += FastDecayChain * Daughter_BR;
+
+ else
+ {
+ ERROR << "Problem in Articial Decay!! check it!!" << endl;
+ exit(1);
+ }
+
+ }
+ }
+ }
+
+ }
+ fNucleiFiliation = CopyfNucleiFiliation;
+ Cleaned = true;
+ for(it_Filiation = fNucleiFiliation.begin(); it_Filiation != fNucleiFiliation.end(); it_Filiation++) // Loop on the mother ZAI
+ {
+ vector<ZAI> DautherList = it_Filiation->second.GetZAIList(); // Get the list of daughter ZAI
+
+ for (int i = 0; i < (int)DautherList.size(); i++) //Loop on daughter
+ if(GoodNuclei.find(DautherList[i]) == GoodNuclei.end() ) // if the ZAI is not in a dealed nuclei (cutted or unknown)
+ Cleaned = false;
+ }
+
+ }
+
+ NormalizeBranchingRatio();
+ DBGL
+}
+
+
+
+//________________________________________________________________________
+void CLASSNucleiFiliation::NormalizeBranchingRatio(double Value)
+{
+ DBGL
+ map<ZAI, IsotopicVector>::iterator it_Filiation;
+ for(it_Filiation = fNucleiFiliation.begin(); it_Filiation != fNucleiFiliation.end(); it_Filiation++)
+ {
+ it_Filiation->second *= Value/it_Filiation->second.GetSumOfAll();
+ }
+
+ DBGL
+}
+
+
+//________________________________________________________________________
+void CLASSNucleiFiliation::NormalizeBranchingRatio(ZAI Mother, double Value)
+{
+ DBGL
+ map<ZAI, IsotopicVector>::iterator it_Filiation = fNucleiFiliation.find(Mother);
+
+ if( it_Filiation != fNucleiFiliation.end())
+ it_Filiation->second *= Value/it_Filiation->second.GetSumOfAll();
+ else
+ WARNING << "Trying to normaliza a Branching Ratio for a Mother wich are not present in the Filiatiuon List...." << endl;
+
+ DBGL
+}
+
+
+
+
+//________________________________________________________________________
+
+vector<ZAI> CLASSNucleiFiliation::GetZAIList() const
+{
+
+ map<ZAI ,IsotopicVector > IsotopicQuantity = GetNucleiFIliation();
+ map<ZAI ,IsotopicVector >::iterator it;
+ vector<ZAI> zailist;
+ for( it = IsotopicQuantity.begin(); it != IsotopicQuantity.end(); it++)
+ zailist.push_back( (*it).first );
+
+ return zailist;
+
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/source/trunk/src/CLASSObject.cxx b/source/trunk/src/CLASSObject.cxx
new file mode 100644
index 000000000..85413d5b2
--- /dev/null
+++ b/source/trunk/src/CLASSObject.cxx
@@ -0,0 +1,24 @@
+#include "CLASSObject.hxx"
+
+using namespace std;
+
+ //________________________________________________________________________
+ //
+ // CLASSObject
+ //
+ //
+ //
+ //
+ //________________________________________________________________________
+
+ClassImp(CLASSObject)
+
+CLASSObject::CLASSObject()
+{
+ fLog = 0;
+}
+
+CLASSObject::CLASSObject(CLASSLogger* log)
+{
+ fLog = log;
+}
\ No newline at end of file
diff --git a/source/trunk/src/DecayDataBank.cxx b/source/trunk/src/DecayDataBank.cxx
new file mode 100644
index 000000000..0a184b104
--- /dev/null
+++ b/source/trunk/src/DecayDataBank.cxx
@@ -0,0 +1,172 @@
+#include "DecayDataBank.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>
+
+
+//________________________________________________________________________
+//
+// DecayDataBank
+//
+//
+//
+//
+//________________________________________________________________________
+
+DecayDataBank::DecayDataBank():CLASSObject(new CLASSLogger("DecayDataBank.log"))
+{
+}
+
+//________________________________________________________________________
+//________________________________________________________________________
+//________________________________________________________________________
+//________________________________________________________________________
+
+//________________________________________________________________________
+
+DecayDataBank::DecayDataBank(string DB_index_file, bool olfreadmethod):CLASSObject(new CLASSLogger("DecayDataBank.log"))
+{
+
+ fDataBaseIndex = DB_index_file;
+ fOldReadMethod = olfreadmethod;
+
+ // Warning
+ INFO << " A EvolutionData has been define :" << endl;
+ INFO << "\t His index is : \"" << DB_index_file << "\"" << endl << endl;
+
+}
+//________________________________________________________________________
+
+DecayDataBank::DecayDataBank(CLASSLogger* log, string DB_index_file, bool olfreadmethod):CLASSObject(log)
+{
+
+ fDataBaseIndex = DB_index_file;
+ fOldReadMethod = olfreadmethod;
+
+ // Warning
+ INFO << " A EvolutionData has been define :" << endl;
+ INFO << "\t His index is : \"" << DB_index_file << "\"" << endl << endl;
+
+}
+
+//________________________________________________________________________
+DecayDataBank::~DecayDataBank()
+{
+
+}
+
+//________________________________________________________________________
+IsotopicVector DecayDataBank::Evolution(const ZAI& zai, double dt)
+{
+DBGL
+ IsotopicVector returnIV;
+
+ map<ZAI ,EvolutionData >::iterator it = fDecayDataBank.find(zai);
+
+ if (it == fDecayDataBank.end() )
+ {
+ ifstream DB_index(fDataBaseIndex.c_str());
+ if( !DB_index)
+ {
+ ERROR << " Can't open \"" << fDataBaseIndex << "\"" << endl;
+ exit (1);
+ }
+ bool zaifind = false;
+ string tmp;
+ getline(DB_index,tmp); // Read first line
+ while (!DB_index.eof())
+ {
+ string line;
+ int start=0;
+ getline(DB_index,line); // Read first line
+ string first=StringLine::NextWord(line,start); // Read first word
+
+ if(first.size()==0) break; // If First word is null.... quit
+
+ int rZ=atoi(first.c_str()); // Get Z
+ int rA=atoi(StringLine::NextWord(line,start).c_str()); // Get A
+ int rI=atoi(StringLine::NextWord(line,start).c_str()); // Get Isomeric State
+
+ if(rZ == zai.Z() && rA == zai.A() && rI == zai.I() )
+ {
+ string file_name = StringLine::NextWord(line,start);
+ EvolutionData evolutionproduct = EvolutionData(GetLog(), file_name, fOldReadMethod);
+#pragma omp critical(DBupdate)
+ {fDecayDataBank.insert( pair<ZAI ,EvolutionData >(zai, evolutionproduct) );}
+ returnIV = evolutionproduct.GetIsotopicVectorAt(dt);
+ zaifind = true;
+ }
+ }
+
+ if(!zaifind)
+ {
+ WARNING << " Oups... Can't Find the ZAI : "
+ << zai.Z() << " " << zai.A() << " " << zai.I() << "!!! It will be considered as stable !!" << endl;
+
+ EvolutionData evolutionproduct = EvolutionData(GetLog()," " , false, zai);
+ {fDecayDataBank.insert( pair<ZAI, EvolutionData >(zai, evolutionproduct) );}
+ returnIV = evolutionproduct.GetIsotopicVectorAt(dt);
+
+
+ }
+
+
+ }
+ else returnIV = (*it).second.GetIsotopicVectorAt(dt);
+
+DBGL
+ return returnIV;
+}
+
+bool DecayDataBank::IsDefine(const ZAI& zai) const
+{
+
+ map<ZAI ,EvolutionData > evolutiondb = (*this).GetDecayDataBank();
+ if (evolutiondb.find(zai) != evolutiondb.end())
+ return true;
+ else
+ return false;
+
+}
+
+
+
+//________________________________________________________________________
+// Get Decay
+//________________________________________________________________________
+IsotopicVector DecayDataBank::GetDecay(IsotopicVector isotopicvector, cSecond t)
+{
+DBGL
+ IsotopicVector IV;
+
+ map<ZAI ,double> isotopicquantity = isotopicvector.GetIsotopicQuantity();
+ map<ZAI ,double >::iterator it;
+ for( it = isotopicquantity.begin(); it != isotopicquantity.end(); it++)
+ {
+ if((*it).second > 0)
+ {
+ IsotopicVector ivtmp = Evolution(it->first, t) * (*it).second ;
+ IV += ivtmp;
+ }
+ }
+
+DBGL
+ return IV;
+}
+//________________________________________________________________________
+//________________________________________________________________________
+//________________________________________________________________________
+//________________________________________________________________________
+//________________________________________________________________________
diff --git a/source/trunk/src/DynamicalSystem.cxx b/source/trunk/src/DynamicalSystem.cxx
new file mode 100644
index 000000000..5ec154052
--- /dev/null
+++ b/source/trunk/src/DynamicalSystem.cxx
@@ -0,0 +1,214 @@
+#include <cstdio>
+#include <iostream>
+#include <fstream>
+#include <cmath>
+#include <vector>
+#include <algorithm>
+
+#include "DynamicalSystem.hxx"
+
+using namespace std;
+
+//________________________________________________________________________
+DynamicalSystem::DynamicalSystem()
+{
+
+ SetPrecision();
+ fHestimate=1.; //this value is change in DynamicalSystem::RungeKutta
+ fHmin=0.;
+ fMaxHdid=-1e30;
+ fMinHdid=1e30;
+ fIsNegativeValueAllowed=true;
+}
+
+//________________________________________________________________________
+DynamicalSystem::DynamicalSystem(const DynamicalSystem & DS)
+{
+ fNVar=DS.fNVar;
+ fPrecision=DS.fPrecision;
+ fHestimate=DS.fHestimate;
+ fHmin=DS.fHmin;
+ fMaxHdid=DS.fMaxHdid;
+ fMinHdid=DS.fMinHdid;
+ fIsNegativeValueAllowed=DS.fIsNegativeValueAllowed;
+}
+
+//________________________________________________________________________
+DynamicalSystem::~DynamicalSystem()
+{
+
+}
+
+//________________________________________________________________________
+void DynamicalSystem::RungeKutta(double *YStart,double t1, double t2, int EquationNumber)
+{
+ //double shortestHalfLife=gMURE->GetShortestHalfLife();
+ fNVar = EquationNumber;
+ int nstp;
+ double t,hnext,hdid,h;
+ double *yscal=new double[fNVar];
+ double *y=new double[fNVar];
+ double *dydt=new double[fNVar];
+
+ const double MAXSTP = 10000;
+ const double TINY = 1.0e-30;
+
+
+
+ if(fabs(t1-t2)<=TINY)
+ cout << "Integration time is 0." << endl;
+
+ t=t1;
+ fHestimate=(t2-t1)/100;
+ //h=(t2 > t1) ? fabs(fHestimate) : -fabs(fHestimate);
+ h=fHestimate;
+
+ // pragma omp parallel for
+ for (int i = 0; i < fNVar; i++) y[i] = YStart[i];
+
+ for ( nstp = 1; nstp <= MAXSTP; nstp++)
+ {
+ BuildEqns(t,y,dydt);
+
+ // pragma omp parallel for
+ for ( int i = 0; i < fNVar; i++)
+ yscal[i] = fabs(y[i]) + fabs(dydt[i]*h)+TINY;
+
+ if ( (t+h-t2) * (t+h-t1) > 0.0 ) h=t2-t;
+
+ AdaptStepSize(y,dydt,&t,h,fPrecision,yscal,&hdid,&hnext);
+
+ if(fMaxHdid<hdid) fMaxHdid=hdid;
+ if(fMinHdid>hdid) fMinHdid=hdid;
+
+ if ((t-t2)*(t2-t1) >= 0.0)
+ {
+ // pragma omp parallel for
+ for (int i=0;i<fNVar;i++) YStart[i]=y[i];
+
+ delete [] dydt;
+ delete [] y;
+ delete [] yscal;
+ //cout << "The maximum step used in RK was "<<fMaxHdid<<" Step NUM in RK was "<<nstp << endl;
+ return;
+ }
+ if (fabs(hnext) <= fHmin)
+ cout << "Step size too small in RungeKutta" << endl;
+
+ h=hnext;
+ }
+ cout << "Too many steps in routine RungeKutta" << endl;
+}
+
+//________________________________________________________________________
+void DynamicalSystem::RK4(double *y, double *dydx, double x, double h, double *yout)
+{
+ //cout<<"Calling Function RK4"<<endl;
+ double xh,hh,h6;
+
+ double *dym=new double[fNVar];
+ double *dyt=new double[fNVar];
+ double *yt=new double[fNVar];
+ hh=h*0.5;
+ h6=h/6.0;
+ xh=x+hh;
+
+ // pragma omp parallel for
+ for (int i=0;i<fNVar;i++) yt[i]=y[i]+hh*dydx[i];
+
+ BuildEqns(xh,yt,dyt);
+
+ // pragma omp parallel for
+ for (int i=0;i<fNVar;i++) yt[i]=y[i]+hh*dyt[i];
+
+ BuildEqns(xh,yt,dym);
+
+ for (int i=0;i<fNVar;i++)
+ {
+ yt[i]=y[i]+h*dym[i];
+ dym[i] += dyt[i];
+ }
+
+ BuildEqns(x+h,yt,dyt);
+ // pragma omp parallel for
+ for (int i=0;i<fNVar;i++)
+ {
+ yout[i]=y[i]+h6*(dydx[i]+dyt[i]+2.0*dym[i]);
+ if(!fIsNegativeValueAllowed && yout[i]<0)
+ {
+ //cout << "Material composition is negative "
+ //<<"i="<<i<<" ("/*<<fEvolvingMaterial->GetComposition()[i]->GetZAI()->PrintName()
+ // */<<") old proportion="<<y[i]<<" new="<<yout[i]
+ //<<". Setting to 0." << endl;
+ yout[i]=0.;
+ }
+ }
+ delete [] yt;
+ delete [] dyt;
+ delete [] dym;
+}
+
+//________________________________________________________________________
+void DynamicalSystem::AdaptStepSize(double *y, double *dydx, double *x, double htry,
+ double eps, double *yscal, double *hdid, double *hnext)
+{
+ //cout<<"Calling Function AdaptStepSize()"<<endl;
+ double xsav,hh,h,temp,errmax;
+ double *dysav=new double[fNVar];
+ double *ysav=new double[fNVar];
+ double *ytemp=new double[fNVar];
+
+ const double PGROW =-0.20;
+ const double PSHRNK =-0.25;
+ const double FCOR =0.06666666 ;
+ const double SAFETY =0.9;
+ const double ERRCON =6.0e-4;
+
+ xsav=(*x);
+
+ // pragma omp parallel for
+ for (int i=0;i<fNVar;i++)
+ {
+ ysav[i]=y[i];
+ dysav[i]=dydx[i];
+ }
+ h=htry;
+ for (;;)
+ {
+ hh=0.5*h;
+ RK4(ysav,dysav,xsav,hh,ytemp);
+ *x=xsav+hh;
+
+ BuildEqns(*x,ytemp,dydx);
+ RK4(ytemp,dydx,*x,hh,y);
+ *x=xsav+h;
+ if (*x == xsav )
+ {
+ //cout << "Step size ("<<h<<") too small in routine AdaptStepSize" << endl;
+ }
+ RK4(ysav,dysav,xsav,h,ytemp);
+ errmax=0.0;
+
+ for (int i=0;i<fNVar;i++)
+ {
+ ytemp[i]=y[i]-ytemp[i];
+ temp=fabs(ytemp[i]/yscal[i]);
+ if (errmax < temp) errmax=temp;
+ }
+ errmax /= eps;
+ if (errmax <= 1.0)
+ {
+ *hdid=h;
+ *hnext=(errmax > ERRCON ? SAFETY*h*exp(PGROW*log(errmax)) : 4.0*h);
+ break;
+ }
+ h=SAFETY*h*exp(PSHRNK*log(errmax));
+ }
+
+ for (int i=0;i<fNVar;i++) y[i] += ytemp[i]*FCOR;
+
+ delete [] ytemp;
+ delete [] dysav;
+ delete [] ysav;
+}
+
diff --git a/source/trunk/src/EquivalenceModel.cxx b/source/trunk/src/EquivalenceModel.cxx
new file mode 100644
index 000000000..3c85ee56f
--- /dev/null
+++ b/source/trunk/src/EquivalenceModel.cxx
@@ -0,0 +1,232 @@
+#include "EquivalenceModel.hxx"
+
+
+
+
+
+
+EquivalenceModel::EquivalenceModel():CLASSObject()
+{
+
+}
+
+EquivalenceModel::EquivalenceModel(CLASSLogger* log):CLASSObject(log)
+{
+
+}
+
+EquivalenceModel::~EquivalenceModel()
+{
+
+}
+//________________________________________________________________________
+vector<double> EquivalenceModel::BuildFuel(double BurnUp, double HMMass,vector<IsotopicVector> FissilArray, vector<IsotopicVector> FertilArray)
+{
+DBGL
+ HMMass*=1e6;//tons to gram
+
+ /****Some initializations**/
+ vector<double> lambda ; //vector of portion of stocks taken (fissile & fertil)
+ for(int i = 0 ; i < (int)FissilArray.size() + (int)FertilArray.size() ; i++ )
+ lambda.push_back(0);
+
+
+ if( (int)FissilArray.size()==0 )
+ { WARNING<<" No fissile stocks available ! Fuel not build"<<endl;
+ lambda[0] = -1;
+ return lambda;
+ }
+
+
+ fOld_Lambda_Tot = 0;
+ fLambda_max = FindLambdaMax( FissilArray, HMMass );
+ DBGV("fLambda_max "<<fLambda_max);
+
+ IsotopicVector Fertile;
+ IsotopicVector Fissile;
+ double AvailablePuMass=0;
+ double PuMassNeeded=1000; //At present time I have no clue what is the requiered Pu mass, I assume at least 1kg is needed
+ double WeightPuContent=0;
+ double MassPrecision=100; //Mass precision is 100 grams
+
+ while( fabs(PuMassNeeded - AvailablePuMass) > MassPrecision )
+ { //Increase the portion of the stock stokID taken, according to the followings variables
+ double DeltaM=PuMassNeeded-AvailablePuMass;
+ GuessLambda( lambda,0,FissilArray.size()-1, DeltaM, FissilArray ,HMMass);
+ if(lambda[0]==-1)
+ {
+ for(int i=0 ; i < (int)FissilArray.size() + (int)FertilArray.size();i++ )
+ lambda[i ]= -1;
+ WARNING<<"Not enought fissile material to build fuel"<<endl;
+ return lambda;
+ }
+ //Build the Plutonium vector from stocks
+ Fissile.Clear();
+
+ for( int i = 0 ; i < (int)FissilArray.size() ; i++ )
+ Fissile +=lambda[i] * FissilArray[i];
+
+ AvailablePuMass = Fissile.GetTotalMass() * 1e6; //in grams
+ //Building complementary Fertile from stocks
+ double FertilMassNeeded = HMMass - AvailablePuMass;
+ double LAMBDA_FERTILE = FindLambdaMax(FertilArray, FertilMassNeeded);
+ SetLambda(lambda,FissilArray.size(), lambda.size()-1,LAMBDA_FERTILE);
+ int j=-1;
+ Fertile.Clear();
+ for(int i = (int)FissilArray.size() ; i < (int)FissilArray.size()+(int)FertilArray.size() ; i++)
+ { j++;
+ Fertile += lambda[i] * FertilArray[j];
+ }
+
+ if( fabs(Fertile.GetTotalMass()*1e6 - FertilMassNeeded) > FertilMassNeeded * 1e-6)//Not enought fertile in stocks
+ {
+ for( int i = 0 ; i < (int)FissilArray.size() + (int)FertilArray.size() ; i++ )
+ lambda[i] = -1; //error code (read by FabricationPlant)
+ WARNING<<"Not enought fertile material to build fuel"<<endl;
+ return lambda;
+ }
+ /*Calcul the quantity of this composition needed to reach the burnup*/
+ double MolarPuContent = GetFissileMolarFraction(Fissile, Fertile, BurnUp);
+ double MeanMolarPu = Fissile.MeanMolar();
+ double MeanMolarDepletedU = Fertile.MeanMolar();
+ double MeanMolar = MeanMolarPu * MolarPuContent + (1-MolarPuContent)*MeanMolarDepletedU;
+ DBGV("MolarPuContent "<<MolarPuContent);
+ WeightPuContent = MolarPuContent * MeanMolarPu / MeanMolar ;
+ PuMassNeeded = WeightPuContent * HMMass ;
+ }
+
+
+ DBGV("Weight percent fissil : "<<PuMassNeeded/HMMass );
+ DBGV("Lambda vector: ");
+ for(int i = 0 ; i < (int)FissilArray.size() + (int)FertilArray.size() ; i++ )
+ DBGV(lambda[i]);
+
+DBGL
+ return lambda;
+}
+//________________________________________________________________________
+void EquivalenceModel::SetLambda(vector<double>& lambda ,int FirstStockID, int LastStockID, double LAMBDA_TOT)
+{
+ if(LAMBDA_TOT > LastStockID - FirstStockID + 1 )
+ {
+ ERROR << " FATAL ERROR " <<endl;
+ exit(0);
+ }
+
+ for(int i = FirstStockID ; i <= LastStockID ; i++) //set to 0 all non touched value (to be sure)
+ lambda[i] = 0 ;
+
+ int PartieEntier = floor( LAMBDA_TOT );
+ double PartieDecimal = LAMBDA_TOT - PartieEntier;
+
+ for(int i=FirstStockID ; i < FirstStockID +PartieEntier ;i++ )
+ lambda[i]=1;
+
+ lambda[FirstStockID + PartieEntier] = PartieDecimal ;
+
+}
+//________________________________________________________________________
+void EquivalenceModel::GuessLambda(vector<double>& lambda,int FirstStockID, int LastStockID, double DeltaM, vector<IsotopicVector> Stocks, double HMMass)
+{
+ double LAMBDA_TOT=0;
+ for (int i = FirstStockID; i <=LastStockID ; i++)
+ LAMBDA_TOT+=lambda[i] ;
+
+
+ if( LAMBDA_TOT == 0 ) //Initialization je cherche les lambda telle que j'ai 5% de HM !!!
+ { //cout<<"INITIALIZATION"<<endl;
+ double MASS = 0 ;
+ int ID_max = 0;
+ while( MASS < 0.05*HMMass )
+ {
+ double StockMass = Stocks[ID_max].GetTotalMass() * 1e6;
+
+ if( StockMass > HMMass )
+ LAMBDA_TOT += 0.01*HMMass/StockMass;
+ else
+ LAMBDA_TOT += 0.01;
+
+ if(LAMBDA_TOT > ID_max+1)
+ {
+ LAMBDA_TOT=ID_max+1;
+ ID_max++ ;
+ if( ID_max >=(int) Stocks.size())
+ { lambda[0]=-1;//error code;
+ break;
+ }
+ }
+
+ SetLambda(lambda,FirstStockID,LastStockID,LAMBDA_TOT );
+
+ IsotopicVector test;
+ int j=0;
+ for(int i=FirstStockID;i<=LastStockID;i++)
+ {
+ IsotopicVector IVTMP = lambda[i] * Stocks[j];
+ test += IVTMP;
+ j++;
+ }
+ MASS = test.GetTotalMass() * 1e6; //in grams
+ //cout<<"LAMBDA_TOT "<<LAMBDA_TOT<<" MASS "<<MASS<<endl;
+ }
+ }
+
+ else if( DeltaM > 0)
+ {
+ fOld_Lambda_Tot = LAMBDA_TOT;
+ LAMBDA_TOT += (fLambda_max - LAMBDA_TOT)/2.;
+
+ if(LAMBDA_TOT > 0.99*fLambda_max) //if we get close to the total of the stocks
+ {
+ double MasseTot=0;
+ for(int i=0;i<(int)Stocks.size();i++)
+ MasseTot+=Stocks[i].GetTotalMass()*1e6;
+
+ if( (fLambda_max -LAMBDA_TOT )*MasseTot < 5 ) //If it remains less than 5gram in stocks => the stocks are not enought
+ lambda[0]=-1;//error code;
+ }
+ //cout<<"DM>0 LAMBDA_TOT "<<LAMBDA_TOT<<endl;
+ else
+ SetLambda(lambda,FirstStockID,LastStockID,LAMBDA_TOT );
+ }
+
+ else if( DeltaM < 0)
+ {
+ //cout <<"fOld_Lambda_Tot[f] "<<fOld_Lambda_Tot[f]<<" LAMBDA_TOT "<<LAMBDA_TOT<< endl;
+ LAMBDA_TOT = LAMBDA_TOT - fabs(fOld_Lambda_Tot - LAMBDA_TOT) / 2. ;
+ //cout<<"DM<0 LAMBDA_TOT "<<LAMBDA_TOT<<endl;
+ SetLambda(lambda,FirstStockID,LastStockID,LAMBDA_TOT );
+ }
+}
+//________________________________________________________________________
+double EquivalenceModel::FindLambdaMax(vector<IsotopicVector> Stocks, double HMMass)
+{
+ //je cherche les lambda telle que j'ai 100% de HMass ou bien la mass total des stocks
+ double LAMBDA = 0;
+ double TotStockMass = 0;
+ for (int i = 0 ; i < (int)Stocks.size() ; i++)
+ TotStockMass += Stocks[i].GetTotalMass() * 1e6;
+
+
+ if(TotStockMass <= HMMass)
+ return Stocks.size() ;
+
+ else
+ {
+ double RemainHM=HMMass;
+ for(int i=0 ;i<(int)Stocks.size();i++ )
+ {
+ if( RemainHM - Stocks[i].GetTotalMass() * 1e6 > 0)
+ {
+ RemainHM -= Stocks[i].GetTotalMass() * 1e6;
+ LAMBDA++;
+ }
+ else
+ { LAMBDA+= RemainHM / (Stocks[i].GetTotalMass() * 1e6);
+ return LAMBDA;
+ }
+ }
+ }
+
+return -1;
+}
\ No newline at end of file
diff --git a/source/trunk/src/EvolutionData.cxx b/source/trunk/src/EvolutionData.cxx
new file mode 100755
index 000000000..371b63b30
--- /dev/null
+++ b/source/trunk/src/EvolutionData.cxx
@@ -0,0 +1,1160 @@
+#include "EvolutionData.hxx"
+
+#include "CLASSLogger.hxx"
+#include "CLASSConstante.hxx"
+#include "StringLine.hxx"
+
+#include <TGraph.h>
+#include <TString.h>
+
+#include <cmath>
+#include <iostream>
+#include <fstream>
+#include <algorithm>
+
+ //________________________________________________________________________
+ //
+ // EvolutionData
+ //
+ //
+ //
+ //
+ //________________________________________________________________________
+
+
+
+struct my_tolower
+{
+ char operator()(char c) const
+ {
+ return std::tolower(static_cast<unsigned char>(c));
+ }
+};
+
+ //To Lower Case, convert any string in lower case
+string tlc(string data)
+{
+ transform(data.begin(), data.end(), data.begin(), my_tolower());
+ return data;
+}
+
+
+
+double Distance(IsotopicVector IV1, EvolutionData Evd1 )
+{
+
+ double d2=0;
+ IsotopicVector IV2 = Evd1.GetIsotopicVectorAt(0.);
+
+ IsotopicVector IVtmp = IV1;
+ map<ZAI ,double> IVtmpIsotopicQuantity = IVtmp.GetIsotopicQuantity();
+ map<ZAI ,double >::iterator it;
+
+ double SumOfXs = 0;
+ for( it = IVtmpIsotopicQuantity.begin(); it != IVtmpIsotopicQuantity.end(); it++)
+ {
+
+ SumOfXs += Evd1.GetXSForAt(0., (*it).first, 1);
+ SumOfXs += Evd1.GetXSForAt(0., (*it).first, 2);
+ SumOfXs += Evd1.GetXSForAt(0., (*it).first, 3);
+
+ }
+
+
+ for( it = IVtmpIsotopicQuantity.begin(); it != IVtmpIsotopicQuantity.end(); it++)
+ {
+ double Z1 = 0.0;
+ double Z2 = 0.0;
+ double XS = 0.0;
+
+ Z1 = IV1.GetZAIIsotopicQuantity( (*it).first );
+ Z2 = IV2.GetZAIIsotopicQuantity( (*it).first );
+ XS = Evd1.GetXSForAt(0., (*it).first, 1)
+ + Evd1.GetXSForAt(0., (*it).first, 2)
+ + Evd1.GetXSForAt(0., (*it).first, 3);
+
+ d2 += pow(Z1-Z2 , 2 ) * pow(XS,2);
+
+ }
+
+ return sqrt(d2)/SumOfXs;
+
+}
+
+double Distance(EvolutionData Evd1, IsotopicVector IV1 )
+{
+ return Distance(IV1,Evd1);
+}
+ //________________________________________________________________________
+ //________________________________________________________________________
+ //________________________________________________________________________
+EvolutionData operator*(EvolutionData const& evol, double F)
+{
+
+ EvolutionData evoltmp;
+
+ map<ZAI ,TGraph* > EvolutionData = evol.GetInventoryEvolution();
+ map<ZAI ,TGraph* >::iterator it;
+ for(it = EvolutionData.begin(); it != EvolutionData.end(); it++)
+ {
+ double X[(*it).second->GetN()];
+ double Y[(*it).second->GetN()];
+
+ for(int i = 0; i < (*it).second->GetN(); i++)
+ {
+ double y;
+ (*it).second->GetPoint( i, X[i], y );
+ Y[i] = y*F;
+ }
+ evoltmp.NucleiInsert( pair<ZAI, TGraph*> ( (*it).first,new TGraph((*it).second->GetN(), X, Y) ) );
+
+ }
+ evoltmp.SetPower(evol.GetPower()*F);
+ evoltmp.SetFissionXS(evol.GetFissionXS());
+ evoltmp.SetCaptureXS(evol.GetCaptureXS());
+ evoltmp.Setn2nXS(evol.Getn2nXS());
+
+
+ return evoltmp;
+
+}
+//________________________________________________________________________
+EvolutionData operator*(double F, EvolutionData const& evol)
+{
+
+ return evol*F;
+
+}
+//________________________________________________________________________
+EvolutionData operator/(EvolutionData const& evol, double F)
+{
+
+ return evol*(1./F);
+
+}
+
+EvolutionData Multiply(EvolutionData const& evol, double F)
+{
+
+ EvolutionData evoltmp;
+ map<ZAI ,TGraph* > EvolutionData = evol.GetInventoryEvolution();
+ map<ZAI ,TGraph* >::iterator it;
+ for(it = EvolutionData.begin(); it != EvolutionData.end(); it++)
+ {
+ double X[(*it).second->GetN()];
+ double Y[(*it).second->GetN()];
+
+ for(int i = 0; i < (*it).second->GetN(); i++)
+ {
+ double y;
+ (*it).second->GetPoint( i, X[i], y );
+ Y[i] = y*F;
+ }
+ evoltmp.NucleiInsert( pair<ZAI, TGraph*> ( (*it).first,new TGraph((*it).second->GetN(), X, Y) ) );
+
+ }
+
+
+ EvolutionData = evol.GetFissionXS();
+ for(it = EvolutionData.begin(); it != EvolutionData.end(); it++)
+ {
+ double X[(*it).second->GetN()];
+ double Y[(*it).second->GetN()];
+
+ for(int i = 0; i < (*it).second->GetN(); i++)
+ {
+ double y;
+ (*it).second->GetPoint( i, X[i], y );
+ Y[i] = y*F;
+ }
+ evoltmp.FissionXSInsert( pair<ZAI, TGraph*> ( (*it).first,new TGraph((*it).second->GetN(), X, Y) ) );
+
+ }
+
+ EvolutionData = evol.GetCaptureXS();
+ for(it = EvolutionData.begin(); it != EvolutionData.end(); it++)
+ {
+ double X[(*it).second->GetN()];
+ double Y[(*it).second->GetN()];
+
+ for(int i = 0; i < (*it).second->GetN(); i++)
+ {
+ double y;
+ (*it).second->GetPoint( i, X[i], y );
+ Y[i] = y*F;
+ }
+ evoltmp.CaptureXSInsert( pair<ZAI, TGraph*> ( (*it).first,new TGraph((*it).second->GetN(), X, Y) ) );
+
+ }
+ EvolutionData = evol.Getn2nXS();
+ for(it = EvolutionData.begin(); it != EvolutionData.end(); it++)
+ {
+ double X[(*it).second->GetN()];
+ double Y[(*it).second->GetN()];
+
+ for(int i = 0; i < (*it).second->GetN(); i++)
+ {
+ double y;
+ (*it).second->GetPoint( i, X[i], y );
+ Y[i] = y*F;
+ }
+ evoltmp.n2nXSInsert( pair<ZAI, TGraph*> ( (*it).first,new TGraph((*it).second->GetN(), X, Y) ) );
+
+ }
+
+ evoltmp.SetPower(evol.GetPower()*F);
+
+
+
+ return evoltmp;
+
+}
+
+//________________________________________________________________________
+EvolutionData Multiply(double F, EvolutionData const& evol)
+{
+
+ return Multiply(evol,F);
+
+}
+
+EvolutionData Sum(EvolutionData const& evol1, EvolutionData const& evol2)
+{
+ EvolutionData EvolSum = evol1;
+ map<ZAI ,TGraph* > EvolutionData1 = EvolSum.GetInventoryEvolution();
+ map<ZAI ,TGraph* > EvolutionData2 = evol2.GetInventoryEvolution();
+ map<ZAI ,TGraph* >::iterator it;
+
+ for(it = EvolutionData2.begin(); it != EvolutionData2.end(); it++)
+ {
+ pair<map<ZAI, TGraph*>::iterator, bool> IResult;
+
+ IResult = EvolutionData1.insert( pair<ZAI, TGraph*> ( *it ) );
+ if(!(IResult.second) )
+ {
+ double X[(*it).second->GetN()];
+ double Y[(*it).second->GetN()];
+ map<ZAI ,TGraph* >::iterator it2 = EvolutionData1.find( (*it).first );
+
+
+ for(int i = 0; i < (*it).second->GetN(); i++)
+ {
+ double y;
+ (*it).second->GetPoint( i, X[i], y );
+ Y[i] = y + (*it2).second->Eval(X[i]);
+ }
+ IResult.first->second = new TGraph((*it).second->GetN(), X, Y);
+
+ }
+
+ }
+ EvolSum.SetInventoryEvolution(EvolutionData1);
+
+
+ EvolutionData1 = evol1.GetFissionXS();
+ EvolutionData2 = evol2.GetFissionXS();
+
+ for(it = EvolutionData2.begin(); it != EvolutionData2.end(); it++)
+ {
+ pair<map<ZAI, TGraph*>::iterator, bool> IResult;
+
+ IResult = EvolutionData1.insert( pair<ZAI, TGraph*> ( *it ) );
+
+ if(!(IResult.second) )
+ {
+ double X[(*it).second->GetN()];
+ double Y[(*it).second->GetN()];
+ map<ZAI ,TGraph* >::iterator it2 = EvolutionData1.find( (*it).first );
+
+
+ for(int i = 0; i < (*it).second->GetN(); i++)
+ {
+ double y;
+ (*it).second->GetPoint( i, X[i], y );
+ Y[i] = y + (*it2).second->Eval(X[i]);
+ }
+ IResult.first->second = new TGraph((*it).second->GetN(), X, Y);
+ }
+ }
+ EvolSum.SetFissionXS(EvolutionData1);
+
+
+ EvolutionData1 = EvolSum.GetCaptureXS();
+ EvolutionData2 = evol2.GetCaptureXS();
+
+ for(it = EvolutionData2.begin(); it != EvolutionData2.end(); it++)
+ {
+ pair<map<ZAI, TGraph*>::iterator, bool> IResult;
+
+ IResult = EvolutionData1.insert( pair<ZAI, TGraph*> ( *it ) );
+
+ if(!(IResult.second) )
+ {
+ double X[(*it).second->GetN()];
+ double Y[(*it).second->GetN()];
+ map<ZAI ,TGraph* >::iterator it2 = EvolutionData1.find( (*it).first );
+
+
+ for(int i = 0; i < (*it).second->GetN(); i++)
+ {
+ double y;
+ (*it).second->GetPoint( i, X[i], y );
+ Y[i] = y + (*it2).second->Eval(X[i]);
+ }
+ IResult.first->second = new TGraph((*it).second->GetN(), X, Y);
+ }
+ }
+ EvolSum.SetCaptureXS(EvolutionData1);
+
+
+ EvolutionData1 = EvolSum.Getn2nXS();
+ EvolutionData2 = evol2.Getn2nXS();
+
+ for(it = EvolutionData2.begin(); it != EvolutionData2.end(); it++)
+ {
+ pair<map<ZAI, TGraph*>::iterator, bool> IResult;
+
+ IResult = EvolutionData1.insert( pair<ZAI, TGraph*> ( *it ) );
+
+ if(!(IResult.second) )
+ {
+ double X[(*it).second->GetN()];
+ double Y[(*it).second->GetN()];
+ map<ZAI ,TGraph* >::iterator it2 = EvolutionData1.find( (*it).first );
+
+
+ for(int i = 0; i < (*it).second->GetN(); i++)
+ {
+ double y;
+ (*it).second->GetPoint( i, X[i], y );
+ Y[i] = y + (*it2).second->Eval(X[i]);
+ }
+ IResult.first->second = new TGraph((*it).second->GetN(), X, Y);
+ }
+ }
+ EvolSum.Setn2nXS(EvolutionData1);
+
+ return EvolSum;
+}
+
+
+
+
+ //________________________________________________________________________
+ //________________________________________________________________________
+
+
+ClassImp(EvolutionData)
+
+
+
+EvolutionData::EvolutionData():CLASSObject()
+{
+ fIsCrossSection = false;
+ fPower = 0;
+ fCycleTime = 0;
+ fKeff = 0;
+ fFlux = 0;
+}
+
+ //________________________________________________________________________
+EvolutionData::EvolutionData(CLASSLogger* log):CLASSObject(log)
+{
+ fIsCrossSection = false;
+ fPower = 0;
+ fCycleTime = 0;
+ fKeff = 0;
+ fFlux = 0;
+}
+
+ //________________________________________________________________________
+EvolutionData::EvolutionData(CLASSLogger* log, string DB_file, bool oldread, ZAI zai):CLASSObject(log)
+{
+
+ fIsCrossSection = false;
+ fDB_file = DB_file;
+ fPower = 0;
+ fCycleTime = 0;
+ fKeff = 0;
+ fFlux = 0;
+
+ if(zai != ZAI(0,0,0))
+ AddAsStable(zai);
+ else
+ ReadDB( fDB_file, oldread); // Read Evolution Produc DB file name
+
+
+
+
+}
+
+ //________________________________________________________________________
+EvolutionData::~EvolutionData()
+{
+
+}
+//________________________________________________________________________
+void EvolutionData::DeleteEvolutionData()
+{
+
+ map<ZAI ,TGraph* >::iterator it_del;
+
+ for( it_del = fInventoryEvolution.begin(); it_del != fInventoryEvolution.end(); it_del++)
+ {
+ delete (*it_del).second;
+ (*it_del).second = 0;
+ }
+ for( it_del = fFissionXS.begin(); it_del != fFissionXS.end(); it_del++)
+ {
+ delete (*it_del).second;
+ (*it_del).second = 0;
+ }
+ for( it_del = fCaptureXS.begin(); it_del != fCaptureXS.end(); it_del++)
+ {
+ delete (*it_del).second;
+ (*it_del).second = 0;
+ }
+ for( it_del = fn2nXS.begin(); it_del != fn2nXS.end(); it_del++)
+ {
+ delete (*it_del).second;
+ (*it_del).second = 0;
+ }
+
+
+ delete fKeff;
+ delete fFlux;
+
+ fInventoryEvolution.clear();
+ fFissionXS.clear();
+ fCaptureXS.clear();
+ fn2nXS.clear();
+
+ fFlux = 0;
+ fKeff = 0;
+
+}
+
+
+bool EvolutionData::NucleiInsert(pair<ZAI, TGraph*> zaitoinsert)
+{
+
+ pair<map<ZAI, TGraph*>::iterator, bool> IResult;
+ IResult = fInventoryEvolution.insert( zaitoinsert);
+ return IResult.second;
+
+}
+
+bool EvolutionData::FissionXSInsert(pair<ZAI, TGraph*> zaitoinsert)
+{
+
+ pair<map<ZAI, TGraph*>::iterator, bool> IResult;
+ IResult = fFissionXS.insert( zaitoinsert);
+ return IResult.second;
+
+}
+
+bool EvolutionData::CaptureXSInsert(pair<ZAI, TGraph*> zaitoinsert)
+{
+
+ pair<map<ZAI, TGraph*>::iterator, bool> IResult;
+ IResult = fCaptureXS.insert( zaitoinsert);
+ return IResult.second;
+
+}
+
+bool EvolutionData::n2nXSInsert(pair<ZAI, TGraph*> zaitoinsert)
+{
+
+ pair<map<ZAI, TGraph*>::iterator, bool> IResult;
+ IResult = fn2nXS.insert( zaitoinsert);
+ return IResult.second;
+
+}
+
+ //________________________________________________________________________
+void EvolutionData::AddAsStable(ZAI zai)
+{
+
+ double time[2] = {0, (500*365.25*3600*24)};
+ double quantity[2] = {1., 1.};
+
+ fInventoryEvolution.insert(pair<ZAI ,TGraph* >(zai, new TGraph(2, time, quantity) ) );
+
+}
+
+ //________________________________________________________________________
+Double_t EvolutionData::Interpolate(double t, TGraph& EvolutionGraph)
+{
+
+ TString fOption;
+ return (double)EvolutionGraph.Eval(t,0x0,fOption);
+
+}
+
+ //________________________________________________________________________
+TGraph* EvolutionData::GetEvolutionTGraph(const ZAI& zai)
+{
+
+ map<ZAI ,TGraph *>::iterator it = GetInventoryEvolution().find(zai) ;
+
+ if ( it != GetInventoryEvolution().end() )
+ return it->second;
+ else
+ return new TGraph();
+
+
+}
+
+ //________________________________________________________________________
+IsotopicVector EvolutionData::GetIsotopicVectorAt(double t)
+{
+
+ IsotopicVector IsotopicVectorTmp;
+ map<ZAI ,TGraph* >::iterator it;
+ for( it = fInventoryEvolution.begin(); it != fInventoryEvolution.end(); it++ )
+ {
+ IsotopicVectorTmp.Add( (*it).first, Interpolate(t, *((*it).second)) );
+ }
+
+
+ return IsotopicVectorTmp;
+}
+
+ //________________________________________________________________________
+double EvolutionData::GetXSForAt(double t, ZAI zai, int ReactionId)
+{
+
+ map<ZAI ,TGraph* > XSEvol;
+ switch(ReactionId)
+ {
+ case 1: XSEvol = GetFissionXS();
+ break;
+ case 2: XSEvol = GetCaptureXS();
+ break;
+ case 3: XSEvol = Getn2nXS();
+ break;
+ default:ERROR << " Wrong ReactionId !!" << endl;
+ exit(1);
+ }
+
+ map<ZAI ,TGraph* >::iterator it = XSEvol.find(zai);
+
+
+ if (it == XSEvol.end())
+ return 0.;
+ else
+ return Interpolate(t, *((*it).second) );
+
+}
+
+
+
+
+
+
+
+
+
+
+
+ //________________________________________________________________________
+
+ //________________________________________________________________________
+
+
+
+
+void EvolutionData::ReadDB(string DBfile, bool oldread)
+{
+
+ if(oldread)
+ {
+
+ OldReadDB(DBfile);
+ return;
+ }
+
+ ReadInfo(); // Read the .info associeted
+
+ ifstream DecayDB(DBfile.c_str()); // Open the File
+ if(!DecayDB) //check if file is correctly open
+ {
+ INFO << " \n Can't open \"" << DBfile << "\"\n" << endl;
+ }
+ vector<double> vTime;
+
+ string line;
+ int start = 0;
+
+ getline(DecayDB, line);
+ if( tlc(StringLine::NextWord(line, start, ' ')) != "time")
+ {
+ ERROR << " Bad Database file : " << DBfile << endl;
+ ERROR << " The first Line MUST be the time line !!!" << endl;
+ exit (1);
+ }
+
+ while(start < (int)line.size())
+ vTime.push_back(atof(StringLine::NextWord(line, start, ' ').c_str()));
+
+ fFinalTime = vTime.back();
+ double Time[vTime.size()];
+ for(int i=0; i < (int)vTime.size();i++)
+ Time[i] = vTime[i];
+ const int NTimeStep = sizeof(Time)/sizeof(double);
+
+
+ enum { Keff=1, Flux, Inv, XSFis, XSCap, XSn2n };
+
+ map<string, int> keyword_map;
+ keyword_map["keff"] = Keff;
+ keyword_map["flux"] = Flux;
+ keyword_map["xsfis"] = XSFis;
+ keyword_map["xscap"] = XSCap;
+ keyword_map["xsn2n"] = XSn2n;
+ keyword_map["inv"] = Inv;
+
+ getline(DecayDB, line);
+ do
+ {
+ start = 0;
+ switch (keyword_map[tlc(StringLine::NextWord(line, start, ' '))])
+ {
+ case Keff:
+ ReadKeff(line, Time, NTimeStep);
+ break;
+
+ case Flux:
+ ReadFlux(line, Time, NTimeStep);
+ break;
+
+ case Inv:
+ ReadInv(line, Time, NTimeStep);
+ break;
+
+ case XSFis:
+ ReadXSFis(line, Time, NTimeStep);
+ break;
+
+ case XSCap:
+ ReadXSCap(line, Time, NTimeStep);
+ break;
+
+ case XSn2n:
+ ReadXSn2n(line, Time, NTimeStep);
+ break;
+
+
+ default:
+ break;
+ }
+
+ getline(DecayDB, line);
+
+
+ }while ( !DecayDB.eof() );
+
+ fHeavyMetalMass = cZAIMass.GetMass( GetIsotopicVectorAt(0.).GetActinidesComposition() );
+
+ DecayDB.close();
+
+
+}
+
+void EvolutionData::ReadKeff(string line, double* time, int NTimeStep)
+{
+ if(fKeff != 0)
+ delete fKeff;
+
+ int start = 0;
+ if( tlc(StringLine::NextWord(line, start, ' ')) != "keff" ) // Check the keyword
+ {
+ ERROR << " Bad keyword : \"keff\" not found !" << endl;
+ exit(1);
+ }
+
+
+ double Keff[NTimeStep];
+
+ int i = 0;
+ while(start < (int)line.size() && i < NTimeStep ) // Read the Data
+ {
+ Keff[i] = atof(StringLine::NextWord(line, start, ' ').c_str()) ;
+ i++;
+ }
+
+
+ fKeff = new TGraph(NTimeStep, time, Keff); // Add the TGraph
+
+
+
+
+}
+
+void EvolutionData::ReadFlux(string line, double* time, int NTimeStep)
+{
+
+ if(fFlux != 0)
+ delete fFlux;
+
+ int start = 0;
+
+ if( tlc(StringLine::NextWord(line, start, ' ')) != "flux" ) // Check the keyword
+ {
+ ERROR << " Bad keyword : \"flux\" not found !" << endl;
+ exit(1);
+ }
+
+ double Flux[NTimeStep];
+
+ int i = 0;
+ while(start < (int)line.size() && i < NTimeStep ) // Read the Data
+ {
+ Flux[i] = atof(StringLine::NextWord(line, start, ' ').c_str()) ;
+ i++;
+ }
+
+
+ fFlux = new TGraph(NTimeStep, time, Flux); // Add the TGraph
+
+
+
+}
+
+
+void EvolutionData::ReadInv(string line, double* time, int NTimeStep)
+{
+
+
+ int start = 0;
+ if( tlc(StringLine::NextWord(line, start, ' ')) != "inv" ) // Check the keyword
+ {
+ ERROR << " Bad keyword : \"inv\" not found !" << endl;
+ exit(1);
+ }
+ // Read the Z A I
+ int Z = atoi(StringLine::NextWord(line, start, ' ').c_str());
+ int A = atoi(StringLine::NextWord(line, start, ' ').c_str());
+ int I = atoi(StringLine::NextWord(line, start, ' ').c_str());
+
+ if(A!=0 && Z!=0)
+ {
+ double Inv[NTimeStep];
+
+ int i = 0;
+ while(start < (int)line.size() && i < NTimeStep ) // Read the Data
+ {
+ Inv[i] = atof(StringLine::NextWord(line, start, ' ').c_str()) ;
+ i++;
+ }
+ // Add the TGraph
+ fInventoryEvolution.insert(pair<ZAI ,TGraph* >(ZAI(Z,A,I), new TGraph(NTimeStep, time, Inv) ) );
+ }
+
+
+}
+
+
+void EvolutionData::ReadXSFis(string line, double* time, int NTimeStep)
+{
+
+
+ int start = 0;
+ if( tlc(StringLine::NextWord(line, start, ' ')) != "xsfis" ) // Check the keyword
+ {
+ ERROR << " Bad keyword : \"xsfis\" not found !" << endl;
+ exit(1);
+ }
+ // Read the Z A I
+ int Z = atoi(StringLine::NextWord(line, start, ' ').c_str());
+ int A = atoi(StringLine::NextWord(line, start, ' ').c_str());
+ int I = atoi(StringLine::NextWord(line, start, ' ').c_str());
+
+ if(A!=0 && Z!=0)
+ {
+ double XSFis[NTimeStep];
+
+ int i = 0;
+ while(start < (int)line.size() && i < NTimeStep ) // Read the Data
+ {
+ XSFis[i] = atof(StringLine::NextWord(line, start, ' ').c_str()) ;
+ i++;
+ }
+
+ // Add the TGraph
+ fFissionXS.insert(pair<ZAI ,TGraph* >(ZAI(Z,A,I), new TGraph(NTimeStep, time, XSFis) ) );
+ }
+
+
+}
+
+void EvolutionData::ReadXSCap(string line, double* time, int NTimeStep)
+{
+
+
+ int start = 0;
+ if( tlc(StringLine::NextWord(line, start, ' ')) != "xscap" ) // Check the keyword
+ {
+ ERROR << " Bad keyword : \"xscap\" not found !" << endl;
+ exit(1);
+ }
+ // Read the Z A I
+ int Z = atoi(StringLine::NextWord(line, start, ' ').c_str());
+ int A = atoi(StringLine::NextWord(line, start, ' ').c_str());
+ int I = atoi(StringLine::NextWord(line, start, ' ').c_str());
+
+ if(A!=0 && Z!=0)
+ {
+ double XSCap[NTimeStep];
+
+ int i = 0;
+ while(start < (int)line.size() && i < NTimeStep ) // Read the Data
+ {
+ XSCap[i] = atof(StringLine::NextWord(line, start, ' ').c_str()) ;
+ i++;
+ }
+
+ // Add the TGraph
+ fCaptureXS.insert(pair<ZAI ,TGraph* >(ZAI(Z,A,I), new TGraph(NTimeStep, time, XSCap) ) );
+ }
+
+
+}
+
+void EvolutionData::ReadXSn2n(string line, double* time, int NTimeStep)
+{
+
+
+ int start = 0;
+ if( tlc(StringLine::NextWord(line, start, ' ')) != "xsn2n" ) // Check the keyword
+ {
+ ERROR << " Bad keyword : \"xsn2n\" not found !" << endl;
+ exit(1);
+ }
+ // Read the Z A I
+ int Z = atoi(StringLine::NextWord(line, start, ' ').c_str());
+ int A = atoi(StringLine::NextWord(line, start, ' ').c_str());
+ int I = atoi(StringLine::NextWord(line, start, ' ').c_str());
+
+ if(A!=0 && Z!=0)
+ {
+ double XSn2n[NTimeStep];
+
+ int i = 0;
+ while(start < (int)line.size() && i < NTimeStep ) // Read the Data
+ {
+ XSn2n[i] = atof(StringLine::NextWord(line, start, ' ').c_str()) ;
+ i++;
+ }
+ // Add the TGraph
+ fn2nXS.insert(pair<ZAI ,TGraph* >(ZAI(Z,A,I), new TGraph(NTimeStep, time, XSn2n) ) );
+ }
+
+
+}
+
+
+
+void EvolutionData::ReadInfo()
+{
+ string InfoDBFile = fDB_file.erase(fDB_file.size()-3,fDB_file.size());
+ InfoDBFile += "Info";
+ ifstream InfoDB_tmp(InfoDBFile.c_str()); // Open the File
+
+ if(!InfoDB_tmp)
+ {
+ InfoDBFile = InfoDBFile.erase(InfoDBFile.size()-4,InfoDBFile.size());
+ InfoDBFile += "info";
+
+ }
+ InfoDB_tmp.close();
+
+ ifstream InfoDB(InfoDBFile.c_str()); // Open the File
+ if(!InfoDB)
+ {
+ WARNING << "!!ERROR!! !!!EvolutionData!!! \n Can't open \"" << InfoDBFile << "\"\n" << endl;
+ }
+
+ int start = 0;
+ string line;
+ getline(InfoDB, line);
+ if ( tlc(StringLine::NextWord(line, start, ' ')) == "reactor")
+ fReactorType = StringLine::NextWord(line, start, ' ');
+
+ start = 0;
+ getline(InfoDB, line);
+ if (tlc(StringLine::NextWord(line, start, ' ')) == "fueltype")
+ fFuelType = StringLine::NextWord(line, start, ' ');
+ start = 0;
+ getline(InfoDB, line);
+ if ( tlc(StringLine::NextWord(line, start, ' ')) == "cycletime")
+ fCycleTime = atof(StringLine::NextWord(line, start, ' ').c_str());;
+ getline(InfoDB, line); // Assembly HM Mass
+ start = 0;
+ getline(InfoDB, line);
+ if ( tlc(StringLine::NextWord(line, start, ' ')) == "constantpower")
+ fPower = atof(StringLine::NextWord(line, start, ' ').c_str());
+ InfoDB.close();
+}
+
+void EvolutionData::OldReadDB(string DBfile)
+{
+
+
+ ifstream DecayDB(DBfile.c_str()); // Open the File
+ if(!DecayDB)
+ {
+ WARNING << " Can't open \"" << DBfile << "\"\n" << endl;
+ }
+ vector<double> vTime;
+ vector<double> vTimeErr;
+
+ string line;
+ int start = 0;
+
+ getline(DecayDB, line);
+ if( StringLine::NextWord(line, start, ' ') != "time")
+ {
+ ERROR << " Bad Database file : " << DBfile << endl;
+ exit (1);
+ }
+
+ while(start < (int)line.size())
+ vTime.push_back(atof(StringLine::NextWord(line, start, ' ').c_str()));
+
+ fFinalTime = vTime.back();
+ double Time[vTime.size()];
+ for(int i=0; i < (int)vTime.size();i++)
+ Time[i] = vTime[i];
+ vector<double> vFlux;
+ start = 0;
+ getline(DecayDB, line);
+ string tmp = StringLine::NextWord(line, start, ' ');
+ if ( tmp == "keff" || tmp == "Keff" )
+ {
+ vector<double> vKeff;
+ while(start < (int)line.size())
+ vKeff.push_back(atof(StringLine::NextWord(line, start, ' ').c_str()));
+
+ double Keff[vKeff.size()];
+ for(int i=0; i < (int)vKeff.size();i++)
+ Keff[i] = vKeff[i];
+
+ fKeff = new TGraph(vTime.size(), Time, Keff);
+
+ start = 0;
+ getline(DecayDB, line);
+ if (StringLine::NextWord(line, start, ' ') == "flux")
+ {
+
+
+ while(start < (int)line.size())
+ vFlux.push_back(atof(StringLine::NextWord(line, start, ' ').c_str()));
+
+ double Flux[vFlux.size()];
+ for(int i=0; i < (int)vFlux.size();i++)
+ Flux[i] = vFlux[i];
+
+ fFlux = new TGraph(vTime.size(), Time, Flux);
+
+ }
+ }
+
+
+ do
+ {
+
+
+ start = 0;
+ int Z = atoi(StringLine::NextWord(line, start, ' ').c_str());
+ int A = atoi(StringLine::NextWord(line, start, ' ').c_str());
+ int I = atoi(StringLine::NextWord(line, start, ' ').c_str());
+
+ if(A!=0 && Z!=0)
+ {
+ double DPQuantity[vTime.size()];
+ for(int k = 0; k < (int)vTime.size(); k++ )
+ DPQuantity[k] = 0;
+
+
+ ZAI zaitmp(Z, A, I);
+ int i=0;
+ while(start < (int)line.size())
+ {
+ double DPQuantityTmp = atof(StringLine::NextWord(line, start, ' ').c_str());
+ DPQuantity[i] = (double)DPQuantityTmp;
+ i++;
+
+ }
+ TGraph* tgraphtmp = new TGraph((int)vTime.size()-1, Time, DPQuantity);
+ fInventoryEvolution.insert(pair<ZAI ,TGraph* >(zaitmp, tgraphtmp) );
+ }
+
+ getline(DecayDB, line);
+ if(line == "" || line == "CrossSection" ) break;
+ }while (!DecayDB.eof() );
+
+ if(line == "CrossSection")
+ {
+ fIsCrossSection = true;
+ getline(DecayDB, line);
+
+ if (line == "Fission")
+ {
+ getline(DecayDB, line);
+
+ do
+ {
+ double DPQuantity[vTime.size()];
+ for(int k = 0; k < (int)vTime.size(); k++ )
+ DPQuantity[k] = 0;
+
+ start = 0;
+ int Z = atoi(StringLine::NextWord(line, start, ' ').c_str());
+ int A = atoi(StringLine::NextWord(line, start, ' ').c_str());
+ int I = atoi(StringLine::NextWord(line, start, ' ').c_str());
+ if(A!=0 && Z!=0)
+ {
+
+
+ ZAI zaitmp(Z, A, I);
+ int i=0;
+ while(start < (int)line.size())
+ {
+ long double DPQuantityTmp = atof(StringLine::NextWord(line, start, ' ').c_str());
+ DPQuantity[i] = (double)DPQuantityTmp;
+ i++;
+
+ }
+ fFissionXS.insert(pair<ZAI ,TGraph* >(zaitmp, new TGraph(vTime.size()-1, Time, DPQuantity) ) );
+ }
+ getline(DecayDB, line);
+ if(line == "" || line == "Capture" ) break;
+ }while ( !DecayDB.eof() );
+ }
+
+ if (line == "Capture")
+ {
+ getline(DecayDB, line); // Nuclei is given with "A Z"
+
+ do
+ {
+ double DPQuantity[vTime.size()];
+ for(int k = 0; k < (int)vTime.size(); k++ )
+ DPQuantity[k] = 0;
+
+
+ start = 0;
+ int Z = atoi(StringLine::NextWord(line, start, ' ').c_str());
+ int A = atoi(StringLine::NextWord(line, start, ' ').c_str());
+ int I = atoi(StringLine::NextWord(line, start, ' ').c_str());
+
+ if(A!=0 && Z!=0)
+ {
+
+
+ ZAI zaitmp(Z, A, I);
+ int i=0;
+ while(start < (int)line.size())
+ {
+ long double DPQuantityTmp = atof(StringLine::NextWord(line, start, ' ').c_str());
+ DPQuantity[i] = (double)DPQuantityTmp;
+ i++;
+
+ }
+ fCaptureXS.insert(pair<ZAI ,TGraph* >(zaitmp, new TGraph(vTime.size()-1, Time, DPQuantity) ) );
+ }
+ getline(DecayDB, line); // Nuclei is given with "A Z"
+ if(line == "" || line == "n2n" ) break;
+ }while ( !DecayDB.eof() );
+
+ }
+
+ if (line == "n2n")
+ {
+
+ getline(DecayDB, line); // Nuclei is given with "A Z"
+
+ do
+ {
+ double DPQuantity[vTime.size()];
+ for(int k = 0; k < (int)vTime.size(); k++ )
+ DPQuantity[k] = 0;
+
+ start = 0;
+ int Z = atoi(StringLine::NextWord(line, start, ' ').c_str());
+ int A = atoi(StringLine::NextWord(line, start, ' ').c_str());
+ int I = atoi(StringLine::NextWord(line, start, ' ').c_str());
+
+ if(A!=0 && Z!=0)
+ {
+
+
+ ZAI zaitmp(Z, A, I);
+ int i=0;
+ while(start < (int)line.size())
+ {
+ long double DPQuantityTmp = atof(StringLine::NextWord(line, start, ' ').c_str());
+ DPQuantity[i] = (double)DPQuantityTmp;
+ i++;
+
+ }
+ fn2nXS.insert(pair<ZAI ,TGraph* >(zaitmp, new TGraph(vTime.size()-1, Time, DPQuantity) ) );
+ }
+ getline(DecayDB, line); // Nuclei is given with "A Z"
+ if(line == "" ) break;
+
+ }while ( !DecayDB.eof() );
+ }
+
+ }
+ DecayDB.close();
+ start = 0;
+
+ string InfoDBFile = DBfile.erase(DBfile.size()-3,DBfile.size());
+ InfoDBFile += "info";
+ ifstream InfoDB(InfoDBFile.c_str()); // Open the File
+ if(!InfoDB)
+ {
+ INFO << " Can't open \"" << InfoDBFile << "\"\n" << endl;
+ return;
+ }
+
+ start = 0;
+ getline(InfoDB, line);
+ if (StringLine::NextWord(line, start, ' ') == "Reactor")
+ fReactorType = StringLine::NextWord(line, start, ' ');
+ start = 0;
+ getline(InfoDB, line);
+ if (StringLine::NextWord(line, start, ' ') == "Fueltype")
+ fFuelType = StringLine::NextWord(line, start, ' ');
+ start = 0;
+ getline(InfoDB, line);
+ if (StringLine::NextWord(line, start, ' ') == "CycleTime")
+ fCycleTime = atof(StringLine::NextWord(line, start, ' ').c_str());;
+ getline(InfoDB, line); // Assembly HM Mass
+ start = 0;
+ getline(InfoDB, line);
+ if (StringLine::NextWord(line, start, ' ') == "ConstantPower")
+ fPower = atof(StringLine::NextWord(line, start, ' ').c_str());
+ getline(InfoDB, line); // cutoff
+ getline(InfoDB, line); // NUmber of Nuclei
+ start = 0;
+ getline(InfoDB, line);
+ if (StringLine::NextWord(line, start, ' ') == "NormalizationFactor")
+ {
+ double NormFactor = atof(StringLine::NextWord(line, start, ' ').c_str());
+ fPower = fPower * NormFactor;
+ double Flux[vFlux.size()];
+ for(int i=0; i < (int)vFlux.size();i++)
+ Flux[i] = vFlux[i];
+
+ fFlux = new TGraph(vTime.size()-1, Time, Flux);
+ }
+ InfoDB.close();
+}
+
+
+
+
+
diff --git a/source/trunk/src/FabricationPlant.cxx b/source/trunk/src/FabricationPlant.cxx
new file mode 100644
index 000000000..804ba2e5f
--- /dev/null
+++ b/source/trunk/src/FabricationPlant.cxx
@@ -0,0 +1,626 @@
+#include "FabricationPlant.hxx"
+
+#include "Storage.hxx"
+#include "Reactor.hxx"
+#include "EvolutionData.hxx"
+#include "DecayDataBank.hxx"
+#include "PhysicsModels.hxx"
+#include "IsotopicVector.hxx"
+#include "Scenario.hxx"
+#include "CLASSLogger.hxx"
+#include "CLASSConstante.hxx"
+
+
+
+
+#include "TMatrixT.h"
+
+#include <sstream>
+#include <string>
+#include <iostream>
+#include <cmath>
+#include <algorithm>
+
+ClassImp(FabricationPlant)
+
+ //________________________________________________________________________
+ //________________________________________________________________________
+ //________________________________________________________________________
+ //
+ // FabricationPlant
+ //
+ //
+ //
+ //
+ //________________________________________________________________________
+ //________________________________________________________________________
+
+
+FabricationPlant::FabricationPlant():CLASSFacility(16)
+{
+ SetName("F_FabricationPLant.");
+
+ fReUsable = 0;
+ fIsReusable = false;
+}
+
+
+FabricationPlant::FabricationPlant(CLASSLogger* log, double fabricationtime):CLASSFacility(log, fabricationtime, 16)
+{
+DBGL
+ SetName("F_FabricationPLant.");
+
+ fFiFo = false;
+ fSubstitutionFuel = false;
+
+ fReUsable = 0;
+ fIsReusable = false;
+
+ INFO << " A FabricationPlant has been define :" << endl;
+ INFO << "\t Chronological Stock Priority has been set! "<< endl;
+ INFO << "\t Fabrication time set to \t " << (double)(GetCycleTime()/3600/24/365.25) << " year" << endl << endl;
+DBGL
+}
+
+
+
+ //________________________________________________________________________
+FabricationPlant::~FabricationPlant()
+{
+
+
+}
+
+
+ //________________________________________________________________________
+void FabricationPlant::SetSeparartionEfficiencyIV(ZAI zai, double factor)
+{
+
+ pair<map<ZAI, double>::iterator, bool> IResult;
+ if(factor > 1) factor = 1;
+
+ if(factor > 0)
+ {
+ IResult = fSeparationLostFraction.GetIsotopicQuantity().insert( pair<ZAI ,double>(zai, 1 - factor));
+ if(!IResult.second)
+ IResult.first->second = 1 - factor;
+ }
+
+}
+
+
+
+ //________________________________________________________________________
+ //_______________________________ Evolution ______________________________
+ //________________________________________________________________________
+void FabricationPlant::Evolution(cSecond t)
+{
+
+ // Check if the FabricationPlant has been created ...
+ if(t == fInternalTime && t != 0) return;
+ // Make the evolution for the FabricationPlant ...
+ FabricationPlantEvolution(t);
+ //Update Inside IsotopicVector
+ UpdateInsideIV();
+ // ... And Finaly update the AbsoluteInternalTime
+ fInternalTime = t;
+
+}
+
+ //________________________________________________________________________
+void FabricationPlant::FabricationPlantEvolution(cSecond t)
+{
+DBGL
+ map<int ,cSecond >::iterator it;
+ for( it = fReactorNextStep.begin(); it!= fReactorNextStep.end(); it++ )
+ {
+ double R_CreactionTime = GetParc()->GetReactor()[ (*it).first ]->GetCreationTime();
+ double R_LifeTime = GetParc()->GetReactor()[ (*it).first ]->GetLifeTime();
+
+ int ReactorId = (*it).first;
+ pair<CLASSFuel, double> R_Fuel = GetParc()->GetReactor()[ReactorId]->GetFuelPlan()->GetFuelAt( t + GetCycleTime() );
+ double R_BU = R_Fuel.second;
+ double R_Power = GetParc()->GetReactor()[ReactorId]->GetPower();
+ double R_HMMass = GetParc()->GetReactor()[ReactorId]->GetHeavyMetalMass();
+ cSecond R_CycleTime = cSecond (R_BU / R_Power * R_HMMass * 1e9 *3600*24);
+ if( R_CycleTime < GetCycleTime())
+ {
+ ERROR << "Reactor Cycle Time is shorter than Fabrication Time of the fuel, we cannot deal it upto now!!!"<< endl;
+ exit(1);
+ }
+
+ if( t + GetCycleTime() >= R_CreactionTime
+ && t + GetCycleTime() < R_CreactionTime + R_LifeTime)
+ {
+ if( (*it).second == t )
+ {
+#pragma omp critical(FuelBuild)
+ {
+ if( R_Fuel.first.GetPhysicsModels() )
+ {
+ BuildFuelForReactor( (*it).first, t );
+ }
+ (*it).second += R_CycleTime;
+ }
+
+ }
+ else if ( (*it).second - R_CycleTime + GetCycleTime() >= t && (*it).second - R_CycleTime < t )
+ {
+ map<int ,IsotopicVector >::iterator it2 = fReactorFuturIV.find( (*it).first );
+ if (it2 != fReactorFuturIV.end())
+ (*it2).second = GetDecay((*it2).second, t - fInternalTime );
+ }
+ }
+ }
+
+
+DBGL
+}
+
+void FabricationPlant::UpdateInsideIV()
+{
+ DBGL
+ fInsideIV = IsotopicVector();
+
+ map< int,IsotopicVector >::iterator it;
+ for( it = fReactorFuturIV.begin(); it != fReactorFuturIV.end(); it++ )
+ fInsideIV += (*it).second;
+
+ DBGL
+}
+
+
+ //________________________________________________________________________
+void FabricationPlant::BuildFuelForReactor(int ReactorId, cSecond t)
+{
+ DBGL
+ if(fFissileStorage.size() == 0)
+ {
+ ERROR << " One need at least one Fissile storage to build fuel " << endl;
+ ERROR << " Use AddFissileStorage to add a stock to provide fissil material... " << endl;
+ ERROR << " One need at least one Fissile storage to build fuel " << endl;
+ exit(1);
+ }
+
+
+
+ double R_HM_Mass = GetParc()->GetReactor()[ ReactorId ]->GetHeavyMetalMass();
+ double R_CycleTime = GetParc()->GetReactor()[ ReactorId ]->GetCycleTime();
+ double R_Power = GetParc()->GetReactor()[ ReactorId ]->GetPower();
+
+ pair<CLASSFuel, double > FuelBU = GetParc()->GetReactor()[ReactorId]->GetFuelPlan()->GetFuelAt(t+GetCycleTime()) ;
+ PhysicsModels FuelType = *FuelBU.first.GetPhysicsModels();
+ double R_BU = FuelBU.second;
+
+ fFissileList = FuelType.GetEquivalenceModel()->GetFissileList();
+ BuildFissileArray();
+
+
+ fFertileList = FuelType.GetEquivalenceModel()->GetFertileList();
+
+
+ if(fFertileStorage.size() != 0) // If the fertile need to be taken in stock
+ BuildFertileArray();
+ else // if their is not stock and the fertile come from outside of the park
+ {
+ fFertileArray.push_back( fFertileList / fFertileList.GetTotalMass() * R_HM_Mass );
+ }
+
+ vector<double> LambdaArray = FuelType.GetEquivalenceModel()->BuildFuel(R_BU, R_HM_Mass, fFissileArray, fFertileArray);
+
+ double LambdaSum = 0;
+ for(int i = 0; i < (int)fFissileArray.size();i++)
+ LambdaSum += LambdaArray[i];
+
+ if(LambdaArray[0] != -1 && LambdaSum > 0 )
+ {
+ IsotopicVector IV = BuildFuelFromEqModel(LambdaArray);
+ EvolutionData EvolDB = FuelType.GenerateEvolutionData( GetDecay(IV,fCycleTime), R_CycleTime, R_Power);
+
+ {
+ pair<map<int, IsotopicVector>::iterator, bool> IResult;
+ IResult = fReactorFuturIV.insert( pair<int, IsotopicVector>(ReactorId, IV) );
+ if(!IResult.second)
+ IResult.first->second = IV;
+ }
+ {
+ pair<map<int, EvolutionData>::iterator, bool> IResult;
+ IResult = fReactorFuturDB.insert( pair<int, EvolutionData>(ReactorId,EvolDB) );
+ if(!IResult.second)
+ IResult.first->second = EvolDB;
+ }
+ fInsideIV += IV;
+ AddCumulativeIVIn(IV);
+ DBGL
+ return;
+ }
+ else
+ {
+
+ if (!fSubstitutionFuel)
+ {
+ {
+ EvolutionData EmptyDB;
+ pair<map<int, EvolutionData>::iterator, bool> IResult;
+ IResult = fReactorFuturDB.insert( pair<int, EvolutionData>(ReactorId,EmptyDB) );
+ if(!IResult.second)
+ IResult.first->second = EmptyDB;
+ }
+ {
+ IsotopicVector EmptyIV;
+ pair<map<int, IsotopicVector>::iterator, bool> IResult;
+ IResult = fReactorFuturIV.insert( pair<int, IsotopicVector>(ReactorId,EmptyIV) );
+ if(!IResult.second)
+ IResult.first->second = EmptyIV;
+ }
+ ResetArrays();
+ }
+ else
+ {
+
+ IsotopicVector IV = fSubstitutionEvolutionData.GetIsotopicVectorAt(0);
+ EvolutionData evolutiondb = fSubstitutionEvolutionData * R_HM_Mass / IV.GetTotalMass();
+
+ IV = IV* R_HM_Mass / IV.GetTotalMass();
+ {
+ pair<map<int, IsotopicVector>::iterator, bool> IResult;
+ IResult = fReactorFuturIV.insert( pair<int, IsotopicVector>(ReactorId, IV) );
+ if(!IResult.second)
+ IResult.first->second = IV;
+ }
+ {
+ pair<map<int, EvolutionData>::iterator, bool> IResult;
+ IResult = fReactorFuturDB.insert( pair<int, EvolutionData>(ReactorId,evolutiondb) );
+ if(!IResult.second)
+ IResult.first->second = evolutiondb;
+ }
+ GetParc()->AddOutIncome( IV );
+ fInsideIV += IV;
+ AddCumulativeIVIn(IV);
+
+
+
+ }DBGL
+ return;
+ }
+DBGL
+}
+
+
+
+void FabricationPlant::BuildFissileArray()
+{
+
+ for(int i = 0; i < (int)fFissileStorage.size(); i++)
+ {
+ vector<IsotopicVector> IVArray = fFissileStorage[i]->GetIVArray();
+
+ for(int j = 0; j < (int)IVArray.size(); j++)
+ {
+
+ IsotopicVector SeparatedIV = Separation(IVArray[j], fFissileList).first;
+
+ if(Norme(SeparatedIV) != 0)
+ {
+ IsotopicVector CooledSeparatedIV = GetDecay( SeparatedIV , GetCycleTime());
+
+ fFissileArray.push_back( CooledSeparatedIV );
+ fFissileArrayAdress.push_back( pair<int,int>(i,j) );
+ fFissileArrayTime.push_back(fFissileStorage[i]->GetIVArrayArrivalTime()[j]);
+ }
+
+ }
+
+ }
+
+ SortArray(0);
+}
+
+
+void FabricationPlant::BuildFertileArray()
+{
+
+
+ for(int i = 0; i < (int)fFertileStorage.size(); i++)
+ {
+ vector<IsotopicVector> IVArray = fFertileStorage[i]->GetIVArray();
+ for(int j = 0; j < (int)IVArray.size(); j++)
+ {
+
+ IsotopicVector SeparatedIV = Separation(IVArray[j], fFertileList).first;
+ if(Norme(SeparatedIV) != 0)
+ {
+ IsotopicVector CooledSeparatedIV = GetDecay( SeparatedIV , GetCycleTime());
+
+ fFertileArray.push_back( CooledSeparatedIV );
+ fFertileArrayAdress.push_back( pair<int,int>(i,j) );
+ fFertileArrayTime.push_back(fFertileStorage[i]->GetIVArrayArrivalTime()[j]);
+ }
+ }
+
+ }
+
+ SortArray(1);
+
+}
+
+void FabricationPlant::SortArray(int i)
+{
+
+
+ vector<IsotopicVector> IVArray;
+ vector<cSecond> TimeArray;
+ vector< pair<int,int> > AdressArray;
+
+ if(i==0) //Fissile
+ {
+ IVArray = fFissileArray;
+ TimeArray = fFissileArrayTime;
+ AdressArray = fFissileArrayAdress;
+ }
+ else if (i==1) //Fertile
+ {
+ IVArray = fFertileArray;
+ TimeArray = fFertileArrayTime;
+ AdressArray = fFertileArrayAdress;
+
+ }
+
+ if(fFiFo)
+ {
+ for(int j = 0; j < (int)TimeArray.size(); j++)
+ {
+ for (int k = j+1; k < (int)TimeArray.size(); k++)
+ {
+ cSecond time_tmp = TimeArray[j];
+ pair<int,int> Adress_tmp = AdressArray[j];
+ IsotopicVector IV_tmp = IVArray[j];
+
+ if(time_tmp > TimeArray[k])
+ {
+ TimeArray[j] = TimeArray[k];
+ TimeArray[k] = time_tmp;
+
+ AdressArray[j] = AdressArray[k];
+ AdressArray[k] = Adress_tmp;
+
+ IVArray[j] = IVArray[k];
+ IVArray[k] = IV_tmp;
+ }
+
+ }
+ }
+ }
+ else
+ {
+ for(int j = 0; j < (int)fFissileArrayTime.size(); j++)
+ {
+ for (int k = j+1; k < (int)TimeArray.size(); k++)
+ {
+ cSecond time_tmp = TimeArray[j];
+ pair<int,int> Adress_tmp = AdressArray[j];
+ IsotopicVector IV_tmp = IVArray[j];
+
+ if(time_tmp < TimeArray[k])
+ {
+ TimeArray[j] = TimeArray[k];
+ TimeArray[k] = time_tmp;
+
+ AdressArray[j] = AdressArray[k];
+ AdressArray[k] = Adress_tmp;
+
+ IVArray[j] = IVArray[k];
+ IVArray[k] = IV_tmp;
+ }
+
+ }
+ }
+ }
+
+
+ if(i==0) //Fissile
+ {
+ fFissileArray = IVArray;
+ fFissileArrayTime = TimeArray;
+ fFissileArrayAdress = AdressArray;
+ }
+ else if (i==1) //Fertile
+ {
+ fFertileArray = IVArray;
+ fFertileArrayTime = TimeArray;
+ fFertileArrayAdress = AdressArray;
+
+ }
+
+
+}
+
+
+
+
+//________________________________________________________________________
+void FabricationPlant::SetSubstitutionFuel(EvolutionData fuel)
+{
+
+ fSubstitutionFuel = true;
+
+ double M0 = cZAIMass.GetMass( fuel.GetIsotopicVectorAt(0.).GetActinidesComposition() );
+ fSubstitutionEvolutionData = fuel / M0;
+
+}
+
+
+ //________________________________________________________________________
+ //_____________________________ Reactor & DB _____________________________
+ //________________________________________________________________________
+ //________________________________________________________________________
+void FabricationPlant::TakeReactorFuel(int Id)
+{
+DBGL
+ IsotopicVector IV;
+ map<int ,IsotopicVector >::iterator it2 = fReactorFuturIV.find( Id );
+
+ AddCumulativeIVOut(it2->second);
+
+ if (it2 != fReactorFuturIV.end())
+ (*it2).second = IV;
+
+ map< int,EvolutionData >::iterator it = fReactorFuturDB.find(Id);
+ (*it).second = EvolutionData();
+
+ UpdateInsideIV();
+DBGL
+}
+
+//________________________________________________________________________
+EvolutionData FabricationPlant::GetReactorEvolutionDB(int ReactorId)
+{
+
+ map< int,EvolutionData >::iterator it = fReactorFuturDB.find(ReactorId);
+ return (*it).second;
+}
+ //________________________________________________________________________
+ //_______________________________ Storage ________________________________
+ //________________________________________________________________________
+IsotopicVector FabricationPlant::BuildFuelFromEqModel(vector<double> LambdaArray)
+{
+DBGL
+ IsotopicVector BuildedFuel;
+ IsotopicVector Lost;
+
+ for(int i = 0; i < (int)fFissileArray.size(); i++)
+ {
+ if(LambdaArray[i] != 0)
+ {
+ int Stor_N = fFissileArrayAdress[i].first;
+ int IV_N = fFissileArrayAdress[i].second;
+
+ pair<IsotopicVector, IsotopicVector> Separated_Lost;
+ Separated_Lost = Separation( fFissileStorage[Stor_N]->GetIVArray()[IV_N]*LambdaArray[i], fFissileList );
+ BuildedFuel += Separated_Lost.first;
+ Lost += Separated_Lost.second;
+
+ }
+ }
+
+ if(fFertileStorage.size() != 0)
+ {
+ for(int i = fFissileArray.size(); i < (int)(fFertileArray.size()+fFissileArray.size()); i++)
+ {
+ if(LambdaArray[i] != 0)
+ {
+ int Stor_N = fFertileArrayAdress[i].first;
+ int IV_N = fFertileArrayAdress[i].second;
+
+ pair<IsotopicVector, IsotopicVector> Separated_Lost;
+ Separated_Lost = Separation( fFertileStorage[Stor_N]->GetIVArray()[IV_N]*LambdaArray[i], fFertileList);
+ BuildedFuel += Separated_Lost.first;
+ Lost += Separated_Lost.second;
+ }
+ }
+ }
+ else
+ BuildedFuel += fFertileArray[0]*LambdaArray.back();
+
+ if(fIsReusable)
+ fReUsable->AddIV(Lost);
+ else
+ GetParc()->AddWaste(Lost);
+
+ DumpStock(LambdaArray);
+
+DBGL
+ return BuildedFuel;
+}
+
+
+ //________________________________________________________________________
+void FabricationPlant::DumpStock(vector<double> LambdaArray)
+{
+DBGL
+
+ for(int i = 0; i < (int)fFissileArray.size(); i++)
+ {
+ if(LambdaArray[i] != 0)
+ {
+ int Stor_N = fFissileArrayAdress[i].first;
+ int IV_N = fFissileArrayAdress[i].second;
+ fFissileStorage[Stor_N]->TakeFractionFromStock( IV_N, LambdaArray[i] );
+ }
+ }
+ if(fFertileStorage.size() != 0)
+ {
+ for(int i = fFissileArray.size(); i < (int)(fFertileArray.size()+fFissileArray.size()); i++)
+ {
+ if(LambdaArray[i] != 0)
+ {
+ int Stor_N = fFertileArrayAdress[i].first;
+ int IV_N = fFertileArrayAdress[i].second;
+
+ fFertileStorage[Stor_N]->TakeFractionFromStock( IV_N, LambdaArray[i] );
+ }
+ }
+ }
+ else
+ GetParc()->AddOutIncome( fFertileArray[0]*LambdaArray.back() );
+
+ ResetArrays();
+
+
+DBGL
+}
+//________________________________________________________________________
+void FabricationPlant::ResetArrays()
+{
+ //Clear the Building Array (Fissile and Fertile)
+ fFissileArray.clear();
+ fFissileArrayTime.clear();
+ fFissileArrayAdress.clear();
+ fFertileArray.clear();
+ fFertileArrayTime.clear();
+ fFertileArrayAdress.clear();
+
+ fFertileList = fFissileList = IsotopicVector();
+}
+//________________________________________________________________________
+pair<IsotopicVector, IsotopicVector> FabricationPlant::Separation(IsotopicVector isotopicvector, IsotopicVector ExtractedList)
+{
+DBGL
+ //[0] = re-use ; [1] = waste
+ IsotopicVector LostInReprocessing = isotopicvector.GetThisComposition(ExtractedList) * fSeparationLostFraction;
+ IsotopicVector SeparatedPart = isotopicvector.GetThisComposition(ExtractedList) - LostInReprocessing;
+ IsotopicVector LostPart = isotopicvector - SeparatedPart;
+DBGL
+ return pair<IsotopicVector, IsotopicVector> (SeparatedPart, LostPart);
+}
+
+
+
+//________________________________________________________________________
+// Get Decay
+//________________________________________________________________________
+IsotopicVector FabricationPlant::GetDecay(IsotopicVector isotopicvector, cSecond t)
+{
+
+ IsotopicVector IV;
+
+ map<ZAI ,double> isotopicquantity = isotopicvector.GetIsotopicQuantity();
+ map<ZAI ,double >::iterator it;
+ for( it = isotopicquantity.begin(); it != isotopicquantity.end(); it++)
+ {
+ if((*it).second > 0)
+ {
+ IsotopicVector ivtmp = fDecayDataBase->Evolution(it->first, t) * (*it).second ;
+ IV += ivtmp;
+ }
+ }
+
+ return IV;
+
+}
+
+
+
+
+
+
diff --git a/source/trunk/src/IrradiationModel.cxx b/source/trunk/src/IrradiationModel.cxx
new file mode 100644
index 000000000..677c43594
--- /dev/null
+++ b/source/trunk/src/IrradiationModel.cxx
@@ -0,0 +1,638 @@
+//
+// IrradiationModel.cxx
+// CLASSSource
+//
+// Created by BaM on 04/05/2014.
+// Copyright (c) 2014 BaM. All rights reserved.
+//
+
+#include "IrradiationModel.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;
+
+IrradiationModel::IrradiationModel():CLASSObject()
+{
+ fShorstestHalflife = 3600.*24*2.;
+ fZAIThreshold = 90;
+
+
+ fDataDirectoryName = getenv("CLASS_PATH");
+ fDataDirectoryName += "/data/";
+ fDataFileName = "chart.JEF3T";
+
+}
+
+IrradiationModel::IrradiationModel(CLASSLogger* log):CLASSObject(log)
+{
+ fShorstestHalflife = 3600.*24*2.;
+ fZAIThreshold = 90;
+
+
+ fDataDirectoryName = getenv("CLASS_PATH");
+ fDataDirectoryName += "/data/";
+ fDataFileName = "chart.JEF3T";
+
+ fNormalDecay = CLASSNucleiFiliation( log );
+ fFastDecay = CLASSNucleiFiliation( log );
+
+ fCaptureReaction= CLASSNucleiFiliation( log );
+ fn2nReaction = CLASSNucleiFiliation( log );
+}
+//________________________________________________________________________
+//________________________________________________________________________
+/* Physics */
+//________________________________________________________________________
+//________________________________________________________________________
+
+
+
+//________________________________________________________________________
+/* Decay Stuff */
+//________________________________________________________________________
+string IrradiationModel::GetDecay(string DecayModes, double &BR,int &Iso, int &StartPos)
+{
+ string header;
+
+ BR = 0;
+ string DecayBR = StringLine::NextWord(DecayModes,StartPos,','); //extraction of the decay mode and the BR
+
+ int ss = 0;
+ string Decay = StringLine::NextWord(DecayBR,ss,':'); //extraction of the decay
+
+
+ if( ss < (int)DecayBR.size() ) //extraction of the BR if exist (i.e. for non stable isotop)
+ BR = atof(DecayBR.substr(ss+1).c_str());
+
+ BR /= 100.; //BR in % -> BR
+
+
+ Iso = 0;
+ if( Decay.find("/",0) < string::npos ) //find the Isomeric state of Daughter
+ {
+ Iso = atoi( Decay.substr( Decay.find("/") + 1 ).c_str() );
+ Decay = Decay.substr( 0, Decay.find("/") );
+ }
+ return Decay;
+}
+
+//________________________________________________________________________
+void IrradiationModel::LoadDecay()
+{
+ DBGL
+
+ // Add TMP and PF as a stable nuclei in the normal decay lsit
+ fNormalDecay.Add(ZAI(-3,-3,-3), IsotopicVector() ); // TMP
+ fNormalDecay.Add(ZAI(-2,-2,-2), IsotopicVector() ); // PF
+
+
+
+ string DataFullPathName = GetDataDirectoryName()+ GetDataFileName();
+ ifstream infile(DataFullPathName.c_str());
+
+ if(!infile)
+ WARNING << " Can't open \"" << DataFullPathName<< "\"" << endl;
+
+
+ do
+ {
+ int A = -1;
+ int Z = -1;
+ int I = -1;
+ string zainame;
+ string Iname;
+ int unkown;
+ double HalfLife;
+ string DecayModes;
+
+ infile >> A >> Z >> zainame >> Iname >> unkown >> HalfLife >> DecayModes;
+ if(Z >= fZAIThreshold )
+ {
+ // Get Isomeric State;
+
+ if(Iname=="gs")
+ I = 0;
+ else
+ {
+ if(Iname[0]=='m')
+ {
+ if( atoi( Iname.substr(1).c_str() )==0 )
+ I = 1;
+ else
+ I = atoi( Iname.substr(1).c_str() );
+ }
+ }
+
+ int start = 0;
+ double branch_test = 0;
+ double branch_test_f = 0;
+
+ ZAI ParentZAI = ZAI(Z,A,I);
+ IsotopicVector DaughtersIV;
+ bool stable = true;
+
+ while(start<int(DecayModes.size()))
+ {
+ double BR;
+ int daughter_A=0;
+ int daughter_N=0;
+ int daughter_Z=0;
+ int Iso=0;
+ int DM=-1;
+ // FPDistribution *FP=0;
+ string decay_name = GetDecay(DecayModes, BR, Iso, start);
+
+ if (decay_name == "s") {DM=0;daughter_N=(A-Z); daughter_Z=Z;BR=1;}
+ if (decay_name == "b-") {DM=1;stable=false; daughter_N=(A-Z)-1; daughter_Z=Z+1;}
+ if (decay_name == "n") {DM=2;stable=false; daughter_N=(A-Z)-1; daughter_Z=Z;}
+ if (decay_name == "nn") {DM=3;stable=false; daughter_N=(A-Z)-2; daughter_Z=Z;}
+ if (decay_name == "b-n"){DM=4;stable=false; daughter_N=(A-Z)-2; daughter_Z=Z+1;}
+ if (decay_name == "p") {DM=5;stable=false; daughter_N=(A-Z); daughter_Z=Z-1;}
+ if (decay_name == "b-a"){DM=6;stable=false; daughter_N=(A-Z)-3; daughter_Z=Z-1;}
+ if (decay_name == "pp") {DM=7;stable=false; daughter_N=(A-Z); daughter_Z=Z-2;}
+ if (decay_name == "ce") {DM=8;stable=false; daughter_N=(A-Z)+1; daughter_Z=Z-1;}
+ if (decay_name == "a") {DM=9;stable=false; daughter_N=(A-Z)-2; daughter_Z=Z-2;}
+ if (decay_name == "cen"){DM=10;stable=false; daughter_N=(A-Z); daughter_Z=Z-1;}
+ if (decay_name == "cep"){DM=11;stable=false; daughter_N=(A-Z)+1; daughter_Z=Z-2;}
+ if (decay_name == "it") {DM=12;stable=false; daughter_N=(A-Z); daughter_Z=Z;}
+ if (decay_name == "db-"){DM=13;stable=false; daughter_N=(A-Z)-2; daughter_Z=Z+2;}
+ if (decay_name == "db+"){DM=14;stable=false; daughter_N=(A-Z)+2; daughter_Z=Z-2;}
+ if (decay_name == "ita"){DM=15;stable=false; daughter_N=(A-Z)-2; daughter_Z=Z-2;}
+ if (decay_name == "sf") {DM=16;stable=false; daughter_N=0; daughter_Z=-2; Iso = -2;}
+ if (decay_name == "cesf"){DM=17;stable=false; daughter_N=0; daughter_Z=-2; Iso = -2;}
+ if (decay_name == "b-sf"){DM=18;stable=false; daughter_N=0; daughter_Z=-2; Iso = -2;}
+
+ daughter_A = daughter_Z + daughter_N;
+ {
+ if( daughter_Z < fZAIThreshold && daughter_Z!=-2 )
+ daughter_A = daughter_Z = Iso = -3;
+ // not spontaneous fission
+
+ if((BR>1e-10) && (!stable))
+ {
+ if(DM <= 15)
+ {
+ ZAI DaughterZAI = ZAI(daughter_Z,daughter_A,Iso);
+ DaughtersIV += BR * DaughterZAI;
+ branch_test+=BR;
+ }
+ else if( DM <= 18)
+ {
+ if(fSpontaneusYield.size() == 0 )
+ {
+ DaughtersIV += 2*BR * ZAI(-2,-2,-2);
+
+ branch_test_f += 2*BR;
+ }
+ else
+ {
+
+ IsotopicVector SpontanuesFissionProduct = fSpontaneusYield.GetFiliation( ParentZAI );
+ if(SpontanuesFissionProduct.GetQuantity(-1,-1,-1) == 0)
+ {
+ DaughtersIV += BR* SpontanuesFissionProduct ;
+ branch_test_f += BR* SpontanuesFissionProduct.GetSumOfAll();
+
+ }
+ else
+ {
+ WARNING << " Unknwon Spontanues yield for ZAI : " << ParentZAI.Z() << " " << ParentZAI.A() << " " << ParentZAI.I() << endl;
+ DaughtersIV += 2*BR * ZAI(-2,-2,-2);
+ branch_test_f += 2*BR;
+ }
+ }
+ }
+
+ }
+
+ }
+ if (DM !=0)
+ stable = false;
+ // End of While loop
+ }
+
+ double btest = fabs(branch_test + branch_test_f/2.-1.0);
+ if ( btest > 1e-8 && !stable )
+ DaughtersIV = DaughtersIV / (branch_test+branch_test_f/2.);
+
+
+
+ if (HalfLife < fShorstestHalflife && !stable)
+ fFastDecay.Add( ParentZAI, DaughtersIV ); // Fill the FastDecay by the Daughter IV
+ else if (stable)
+ {
+ fNormalDecay.Add(ParentZAI, IsotopicVector() ); // Fill the NormalDecay with a empty IV (mother is stable)
+ }
+ else
+ {
+ fNormalDecay.Add( ParentZAI, DaughtersIV ); // FIll the NormalDecay with the daughter IV scaled by the decay constante.
+ fDecayConstante += ParentZAI*log(2)/HalfLife;
+ }
+ }
+
+ } while (!infile.eof());
+ DBGL
+
+ //Build the Matrix index :
+ fReverseMatrixIndex = fNormalDecay.GetZAIList();
+ for(int i = 0; i< (int)fReverseMatrixIndex.size(); i++)
+ fMatrixIndex.insert(pair<ZAI, int> (fReverseMatrixIndex[i], i) );
+
+ DBGL
+ fFastDecay.SelfFiliationCleanUp(fMatrixIndex);
+ DBGL
+ fNormalDecay.FiliationCleanUp(fMatrixIndex, fFastDecay);
+
+
+ DBGL
+}
+
+
+void IrradiationModel::BuildDecayMatrix()
+{
+ DBGL
+
+ fDecayMatrix.Clear();
+
+ fDecayMatrix.ResizeTo( fMatrixIndex.size(), fMatrixIndex.size() );
+ for(int i = 0; i < (int)fMatrixIndex.size(); i++)
+ for(int j = 0; j < (int)fMatrixIndex.size(); j++)
+ fDecayMatrix[i][j] = 0;
+
+ for(int i = 0; i < (int)fReverseMatrixIndex.size(); i++)
+ {
+
+ IsotopicVector DaughterIV = fNormalDecay.GetFiliation(fReverseMatrixIndex[i]);
+ (*this).GetNuclearProcessMatrix(fDecayMatrix, fReverseMatrixIndex[i], DaughterIV, fDecayConstante.GetQuantity(fReverseMatrixIndex[i]) );
+
+ }
+
+ DBGL
+}
+
+
+void IrradiationModel::GetNuclearProcessMatrix(TMatrixT<double> &NuclearProcessMatrix, ZAI Mother, IsotopicVector ProductedIV, double XSValue)
+{
+ DBGL
+
+ vector<ZAI> ProductedZAIList = ProductedIV.GetZAIList();
+
+ if(fMatrixIndex.find(Mother) != fMatrixIndex.end())
+ {
+
+
+ int i = fMatrixIndex[Mother];
+
+ NuclearProcessMatrix[i][i] -= XSValue;
+
+ for(int j = 0; j < (int)ProductedZAIList.size(); j++)
+ {
+ if(fMatrixIndex.find(ProductedZAIList[j]) != fMatrixIndex.end() )
+ { NuclearProcessMatrix[fMatrixIndex[ ProductedZAIList[j] ]][i] += ProductedIV.GetQuantity(ProductedZAIList[j])*XSValue;
+ }
+ else
+ NuclearProcessMatrix[0][i] += ProductedIV.GetQuantity(ProductedZAIList[j])*XSValue;
+ }
+ }
+ else
+ WARNING << " Can't have nuclear process on this nucleus, ZAI : " << Mother.Z() << " " << Mother.A() << " " << Mother.I() << " its halflife seems to be below the threshold!" << endl;
+
+ DBGL
+}
+
+
+//________________________________________________________________________
+/* Fission Stuff */
+//________________________________________________________________________
+void IrradiationModel::SetFissionEnergy(ZAI zai, double E)
+{
+ pair<map<ZAI, double>::iterator, bool> IResult;
+ IResult = fFissionEnergy.insert( pair<ZAI ,double>(zai, E));
+ if(!IResult.second)
+ IResult.first->second = E;
+
+}
+
+//________________________________________________________________________
+void IrradiationModel::SetFissionEnergy(string FissionEnergyFile)
+{
+ ifstream FissionFile( FissionEnergyFile.c_str() ); // Open the File
+ if(!FissionFile) //check if file is correctly open
+ WARNING << " Can't open \"" << FissionEnergyFile << "\" File" << endl;
+
+ do
+ {
+ int Z = 0;
+ int A = 0;
+ int I = 0;
+ double E = 0;
+ FissionFile >> Z >> A >> I >> E;
+ SetFissionEnergy(Z, A, I, E);
+ } while (!FissionFile.eof());
+}
+
+//________________________________________________________________________
+CLASSNucleiFiliation IrradiationModel::ReadFPYield(string Yield)
+{
+ DBGL
+
+ CLASSNucleiFiliation MyYield = CLASSNucleiFiliation( fLog );
+
+ ifstream infile(Yield.c_str());
+ if(!infile)
+ WARNING << " Can't open \"" << Yield<< "\" File !!!" << endl;
+
+
+ string line;
+ int start = 0;
+
+ getline(infile, line);
+ vector<ZAI> Fissile;
+ vector<IsotopicVector> FPYields;
+ do
+ {
+ int Z = atof(StringLine::NextWord(line, start, ' ').c_str());
+ int A = atof(StringLine::NextWord(line, start, ' ').c_str());
+ int I = atof(StringLine::NextWord(line, start, ' ').c_str());
+ Fissile.push_back(ZAI(Z,A,I));
+ FPYields.push_back( IsotopicVector() );
+
+ }while(start < (int)line.size()-1);
+
+ getline(infile, line);
+ do
+ {
+ start = 0;
+
+ int Z = atof(StringLine::NextWord(line, start, ' ').c_str());
+ int A = atof(StringLine::NextWord(line, start, ' ').c_str());
+ int I = atof(StringLine::NextWord(line, start, ' ').c_str());
+ int i=0;
+
+ do
+ {
+ double Yield_values = atof(StringLine::NextWord(line, start, ' ').c_str());
+
+ FPYields[i] += Yield_values * ZAI(Z,A,I);
+
+ i++;
+
+ }while(start < (int)line.size()-1);
+
+ getline(infile, line);
+
+ } while (!infile.eof());
+
+
+ for(int i=0 ; i<(int) Fissile.size() ; i++) // Fill the CLASSNucleiFiliation
+ MyYield.Add( Fissile[i], FPYields[i] );
+
+
+ DBGL
+ return MyYield;
+}
+
+//________________________________________________________________________
+void IrradiationModel::LoadFPYield(string SpontaneusYield, string ReactionYield)
+{
+ fSpontaneusYieldFile = SpontaneusYield;
+ fReactionYieldFile = ReactionYield;
+ fZAIThreshold = 0;
+}
+//________________________________________________________________________
+void IrradiationModel::NuclearDataInitialization()
+{
+ DBGL
+
+ if(fSpontaneusYieldFile!="")
+ fSpontaneusYield = ReadFPYield(fSpontaneusYieldFile);
+ if(fReactionYieldFile!="")
+ fReactionYield = ReadFPYield(fReactionYieldFile);
+
+ LoadDecay();
+
+
+ if(fSpontaneusYieldFile!="")
+ fSpontaneusYield.FiliationCleanUp(fMatrixIndex, fFastDecay); // remove the cutted nuclei....
+ if(fReactionYieldFile!="")
+ fReactionYield.FiliationCleanUp(fMatrixIndex, fFastDecay); // remove the cutted nuclei....
+
+
+ BuildDecayMatrix();
+ BuildReactionFiliation();
+
+ DBGL
+}
+
+
+void IrradiationModel::BuildReactionFiliation()
+{
+ DBGL
+
+ // (n,Gamma) Special Reaction.....
+ {
+ // 241Am(n,Gamma)
+ {
+ fCaptureReaction.Add( ZAI(95,241,0), ZAI(96,242,0) * 0.8733*0.827 ); //directly cut the Am242 as in MURE
+ fCaptureReaction.Add( ZAI(95,241,0), ZAI(94,242,0) * 0.8733*0.173 ); //directly cut the Am242 as in MURE
+ fCaptureReaction.Add( ZAI(95,241,0), ZAI(95,242,1) * 0.1267 );
+ }
+
+ if(fReactionYieldFile!="")
+ {// 165Ho(n,Gamma)
+ {
+ fCaptureReaction.Add( ZAI(67,165,0), ZAI(68,166,0) * 0.9490 ); //
+ fCaptureReaction.Add( ZAI(67,165,0), ZAI(67,166,1) * 0.0510 ); //
+ }
+
+ // 147Pm(n,Gamma)
+ {
+ fCaptureReaction.Add( ZAI(61,147,0), ZAI(61,148,0) * 0.5330 );
+ fCaptureReaction.Add( ZAI(61,147,0), ZAI(61,148,1) * 0.4670 );
+
+ }
+
+ // 109Ag(n, Gamma)
+ {
+ fCaptureReaction.Add( ZAI(47,109,0), ZAI(48,110,0) * 0.9970*0.9508);
+ fCaptureReaction.Add( ZAI(47,109,0), ZAI(46,110,0) * 0.0030*0.9508);
+ fCaptureReaction.Add( ZAI(47,109,0), ZAI(47,110,1) *0.0492);
+ }
+
+
+ // 107Ag(n, Gamma)
+ {
+ fCaptureReaction.Add( ZAI(47,107,0), ZAI(48,108,0) * 0.9715*0.9895 );
+ fCaptureReaction.Add( ZAI(47,107,0), ZAI(46,108,0) * 0.0285*0.9895 );
+ fCaptureReaction.Add( ZAI(47,107,0), ZAI(47,108,1) *0.0105 );
+ }
+ }
+ }
+
+
+
+ // (n,2n) Special Reaction.....
+ {
+ // 237Np(n,2n)
+ {
+ fn2nReaction.Add(ZAI(93,237,0), ZAI(93,236,0) * 0.2 );
+ fn2nReaction.Add(ZAI(93,237,0), ZAI(93,236,1) * 0.8 );
+ }
+
+ }
+
+
+
+ for(int i = 2; i < (int)fReverseMatrixIndex.size(); i++) // Start at 2 to skeep "TMP" ZAI and "PF" ZAI
+ {
+
+ int Z = fReverseMatrixIndex[i].Z();
+ int A = fReverseMatrixIndex[i].A();
+
+ if(fCaptureReaction.GetFiliation(fReverseMatrixIndex[i]).GetQuantity(ZAI(-1,-1,-1)) == 1 )
+ {
+ fCaptureReaction.Add(fReverseMatrixIndex[i], ZAI(Z,A+1)*1);
+ }
+ if(fn2nReaction.GetFiliation(fReverseMatrixIndex[i]).GetQuantity(ZAI(-1,-1,-1)) == 1 )
+ {
+ if(A>1)
+ fn2nReaction.Add(fReverseMatrixIndex[i], ZAI(Z,A-1)*1);
+ }
+ }
+
+
+ fCaptureReaction.FiliationCleanUp(fMatrixIndex, fFastDecay); // clean the filiation link
+ fCaptureReaction.NormalizeBranchingRatio(); // normalize it
+
+
+ fn2nReaction.FiliationCleanUp(fMatrixIndex, fFastDecay); // clean the filiation link
+ fn2nReaction.NormalizeBranchingRatio(); // normalize it
+
+
+ DBGL
+}
+
+
+//________________________________________________________________________
+/* Reaction Stuff */
+//________________________________________________________________________
+TMatrixT<double> IrradiationModel::GetFissionXsMatrix(EvolutionData EvolutionDataStep,double TStep)
+{
+ DBGL
+ TMatrixT<double> FissionMatrix = TMatrixT<double>( fReverseMatrixIndex.size(),fReverseMatrixIndex.size() );
+ for(int i = 0; i < (int)fReverseMatrixIndex.size(); i++)
+ for(int j = 0; j < (int)fReverseMatrixIndex.size(); j++)
+ FissionMatrix[i][j] = 0;
+
+ // ---------------- A(n,.) X+Y
+
+ map<ZAI ,TGraph* > FissionXS = EvolutionDataStep.GetFissionXS();
+ map<ZAI ,TGraph* >::iterator it_XS;
+
+ for(it_XS = FissionXS.begin() ; it_XS != FissionXS.end(); it_XS++) //loop on fissionable nuclei
+ {
+ ZAI Mother = (*it_XS).first; // Note the Mother ZAI (not necessary but help for reading the code)
+ double XS_Value = (*it_XS).second->Eval(TStep) * 1e-24; // Get Cross section values
+
+ IsotopicVector FissionProductIV = fReactionYield.GetFiliation(Mother); // Get the Isotopicvector produced by the reaction
+
+ if(FissionProductIV.GetQuantity(ZAI(-1,-1,-1)) != 1) // Check if ZAI is dealed
+ GetNuclearProcessMatrix( FissionMatrix, Mother, FissionProductIV, XS_Value ); // add the Nuclear process in the Reaction Matrix
+ else
+ {
+ WARNING << "Don't have fission Yield for this nuclei, ZAI : " << Mother.Z() << " " << Mother.A() << " " << Mother.I() << endl;
+ GetNuclearProcessMatrix(FissionMatrix, Mother, ZAI(-2, -2, -2) * 2 , XS_Value ); // add the Nuclear process in the Reaction Matrix
+ }
+
+
+ }
+
+ DBGL
+ return FissionMatrix;
+}
+
+//________________________________________________________________________
+TMatrixT<double> IrradiationModel::GetCaptureXsMatrix(EvolutionData EvolutionDataStep,double TStep)
+{
+ DBGL
+ TMatrixT<double> CaptureMatrix = TMatrixT<double>( fReverseMatrixIndex.size(),fReverseMatrixIndex.size() );
+ for(int i = 0; i < (int)fReverseMatrixIndex.size(); i++)
+ for(int j = 0; j < (int)fReverseMatrixIndex.size(); j++)
+ CaptureMatrix[i][j] = 0;
+
+ // ---------------- A(n,Gamma) A+1
+
+ map<ZAI ,TGraph* > CaptureXS = EvolutionDataStep.GetCaptureXS();
+ map<ZAI ,TGraph* >::iterator it_XS;
+
+ for(it_XS = CaptureXS.begin() ; it_XS != CaptureXS.end(); it_XS++) //loop on nuclei
+ {
+ ZAI Mother = (*it_XS).first; // Note the Mother ZAI (not necessary but help for reading the code)
+ double XS_Value = (*it_XS).second->Eval(TStep) * 1e-24; // Get Cross section values
+
+ IsotopicVector CaptureProductIV = fCaptureReaction.GetFiliation(Mother); // Get the Isotopicvector produced by the reaction
+
+ if(CaptureProductIV.GetQuantity(ZAI(-1,-1,-1)) != 1) // Check if ZAI is dealed
+ GetNuclearProcessMatrix(CaptureMatrix, Mother, CaptureProductIV, XS_Value ); // add the Nuclear process in the Reaction Matrix
+ else
+ WARNING << "Can't have capture reaction on this nuclei, ZAI : " << Mother.Z() << " " << Mother.A() << " " << Mother.I() << endl;
+
+
+
+
+
+ }
+
+ DBGL
+ return CaptureMatrix;
+}
+
+
+//________________________________________________________________________
+TMatrixT<double> IrradiationModel::Getn2nXsMatrix(EvolutionData EvolutionDataStep,double TStep)
+{
+ DBGL
+
+ TMatrixT<double> n2nMatrix = TMatrixT<double>( fReverseMatrixIndex.size(),fReverseMatrixIndex.size() );
+ for(int i = 0; i < (int)fReverseMatrixIndex.size(); i++)
+ for(int j = 0; j < (int)fReverseMatrixIndex.size(); j++)
+ n2nMatrix[i][j] = 0;
+
+ // ---------------- A(n,2n) A-1
+
+ map<ZAI ,TGraph* > CaptureXS = EvolutionDataStep.Getn2nXS();
+ map<ZAI ,TGraph* >::iterator it_XS;
+
+ for(it_XS = CaptureXS.begin() ; it_XS != CaptureXS.end(); it_XS++) //loop on nuclei
+ {
+ ZAI Mother = (*it_XS).first; // Note the Mother ZAI (not necessary but help for reading the code)
+ double XS_Value = (*it_XS).second->Eval(TStep) * 1e-24; // Get Cross section values
+
+ IsotopicVector n2nProductIV = fn2nReaction.GetFiliation(Mother); // Get the Isotopicvector produced by the reaction
+
+ if(n2nProductIV.GetQuantity(ZAI(-1,-1,-1)) != 1) // Check if ZAI is dealed
+ GetNuclearProcessMatrix(n2nMatrix, Mother, n2nProductIV, XS_Value ); // add the Nuclear process in the Reaction Matrix
+ else
+ WARNING << "Can't have n,2n reaction on this nuclei, ZAI : " << Mother.Z() << " " << Mother.A() << " " << Mother.I() << endl;
+
+
+ }
+
+ DBGL
+ return n2nMatrix;
+}
+
diff --git a/source/trunk/src/IsotopicVector.cxx b/source/trunk/src/IsotopicVector.cxx
new file mode 100755
index 000000000..960db6952
--- /dev/null
+++ b/source/trunk/src/IsotopicVector.cxx
@@ -0,0 +1,681 @@
+#include "IsotopicVector.hxx"
+
+#include "CLASSLogger.hxx"
+#include "CLASSConstante.hxx"
+
+
+#include <cmath>
+#include <iostream>
+#include <fstream>
+#include <string>
+#include <algorithm>
+//________________________________________________________________________
+//________________________________________________________________________
+//
+//
+//
+// IsotopicVector
+//
+//
+//________________________________________________________________________
+//________________________________________________________________________
+
+
+
+
+//________________________________________________________________________
+//__________________________Operator Overlaoding__________________________
+//________________________________________________________________________
+
+
+
+//____________________________General Operator____________________________
+//________________________________________________________________________
+
+ClassImp(IsotopicVector)
+
+double Norme(IsotopicVector IV1, int DistanceType, IsotopicVector DistanceParameter)
+{
+
+ IsotopicVector IV;
+
+ return Distance(IV1, IV, DistanceType,DistanceParameter);
+
+}
+double DistanceStandard(IsotopicVector IV1, IsotopicVector IV2)
+{
+
+ double d2=0;
+ IsotopicVector IVtmp = IV1 + IV2;
+ map<ZAI ,double> IVtmpIsotopicQuantity = IVtmp.GetIsotopicQuantity();
+ map<ZAI ,double >::iterator it;
+ for( it = IVtmpIsotopicQuantity.begin(); it != IVtmpIsotopicQuantity.end(); it++)
+ {
+ double Z1 = IV1.GetZAIIsotopicQuantity( (*it).first );
+ double Z2 = IV2.GetZAIIsotopicQuantity( (*it).first );
+ d2 += pow(Z1-Z2 , 2 );
+ }
+ return sqrt(d2);
+
+}
+double DistanceAdjusted(IsotopicVector IV1, IsotopicVector IV2, IsotopicVector DistanceParameter)
+{
+
+ double d2=0;
+ IsotopicVector IVtmp = IV1 + IV2;
+ map<ZAI ,double> IVtmpIsotopicQuantity = IVtmp.GetIsotopicQuantity();
+ map<ZAI ,double >::iterator it;
+ for( it = IVtmpIsotopicQuantity.begin(); it != IVtmpIsotopicQuantity.end(); it++)
+ {
+ double Z1 = IV1.GetZAIIsotopicQuantity( (*it).first );
+ double Z2 = IV2.GetZAIIsotopicQuantity( (*it).first );
+ double lambda = DistanceParameter.GetZAIIsotopicQuantity( (*it).first );
+ d2 += lambda*abs(Z1-Z2);
+ }
+ return d2;
+
+}
+
+
+
+double Distance(IsotopicVector IV1, IsotopicVector IV2, int DistanceType, IsotopicVector DistanceParameter)
+{
+
+ if(DistanceType==0)
+ {
+ return DistanceStandard(IV1,IV2);
+ }
+ else if(DistanceType==1||DistanceType==2){
+ return DistanceAdjusted(IV1,IV2,DistanceParameter);
+ }
+ else
+ {
+ cout << " DistanceType defined by the user isn't recognized by the code" << endl;
+
+ exit(1);
+ }
+
+}
+
+
+double RelativDistance(IsotopicVector IV1, IsotopicVector IV2 )
+{
+
+ double d2 = 0;
+
+ IsotopicVector IVtmp = IV1 + IV2;
+ map<ZAI ,double> IVtmpIsotopicQuantity = IVtmp.GetIsotopicQuantity();
+
+ double Z1total = 0;
+ double Z2total = 0;
+ map<ZAI ,double >::iterator it;
+ for( it = IVtmpIsotopicQuantity.begin(); it != IVtmpIsotopicQuantity.end(); it++)
+ {
+ Z1total += IV1.GetZAIIsotopicQuantity( (*it).first );
+ Z2total += IV2.GetZAIIsotopicQuantity( (*it).first );
+ }
+ for( it = IVtmpIsotopicQuantity.begin(); it != IVtmpIsotopicQuantity.end(); it++)
+ {
+ double Z1 = IV1.GetZAIIsotopicQuantity( (*it).first );
+ double Z2 = IV2.GetZAIIsotopicQuantity( (*it).first );
+ d2 += pow( (Z1/Z1total - Z2/Z2total) , 2 );
+ }
+
+
+ return sqrt(d2);
+}
+
+
+IsotopicVector operator+(IsotopicVector const& IVa, IsotopicVector const& IVb)
+{
+
+ IsotopicVector IVtmp;
+ IVtmp = IVa;
+ return IVtmp += IVb;
+}
+
+//________________________________________________________________________
+IsotopicVector operator-(IsotopicVector const& IVa, IsotopicVector const& IVb)
+{
+
+ IsotopicVector IVtmp;
+ IVtmp = IVa;
+ return IVtmp -= IVb;
+}
+
+
+//________________________________________________________________________
+IsotopicVector operator*(ZAI const& zai, double F)
+{
+
+ IsotopicVector IVtmp;
+
+ IVtmp.Add( zai, F);
+ return IVtmp;
+}
+
+
+//________________________________________________________________________
+IsotopicVector operator/(ZAI const& zai, double F)
+{
+ IsotopicVector IVtmp;
+
+ IVtmp.Add( zai, 1./F);
+
+ return IVtmp;
+}
+
+
+//________________________________________________________________________
+IsotopicVector operator*(double F, IsotopicVector const& IVA) {return IVA*F;}
+
+//________________________________________________________________________
+IsotopicVector operator*(IsotopicVector const& IVA, double F)
+{
+
+ IsotopicVector IV = IVA;
+ IV.Multiply(F);
+ return IV;
+}
+
+//________________________________________________________________________
+IsotopicVector operator*(double F, ZAI const& zai)
+{
+ return zai*F;
+}
+
+//________________________________________________________________________
+IsotopicVector operator*(IsotopicVector const& IVa, IsotopicVector const& IVb)
+{
+
+ IsotopicVector IVtmp;
+ IVtmp = IVa;
+ IVtmp *= IVb;
+ return IVtmp;
+}
+
+//________________________________________________________________________
+IsotopicVector operator/(IsotopicVector const& IVA, double F)
+{
+
+ IsotopicVector IV = IVA;
+ IV.Multiply(1./F);
+ return IV;
+}
+
+
+//____________________________InClass Operator____________________________
+
+//________________________________________________________________________
+IsotopicVector& IsotopicVector::operator+=(const IsotopicVector& IVa)
+{
+
+ Add(IVa);
+ return *this;
+
+}
+
+//________________________________________________________________________
+IsotopicVector& IsotopicVector::operator-=(const IsotopicVector& IVa)
+{
+
+ Remove(IVa);
+ return *this;
+
+}
+//________________________________________________________________________
+IsotopicVector& IsotopicVector::operator*=(const IsotopicVector& IVa)
+{
+ map<ZAI, double> IVA_isotopicquantity = IVa.GetIsotopicQuantity();
+
+ map<ZAI, double> isotopicquantity = (*this).GetIsotopicQuantity(); // get the isotopic quantity to loop on it
+ map<ZAI, double>::iterator isotopicIT; // iterator on a isotopic quantity map
+
+ for(isotopicIT = isotopicquantity.begin(); isotopicIT != isotopicquantity.end(); isotopicIT++) // loop on the isotopicquantity...
+ {
+ map<ZAI, double>::iterator IVa_isotopicIT = IVA_isotopicquantity.find( (*isotopicIT).first );
+
+ if( IVa_isotopicIT != IVA_isotopicquantity.end() )
+ (*isotopicIT).second *= (*IVa_isotopicIT).second;
+ else
+ (*isotopicIT).second = 0;
+
+ }
+
+ fIsotopicQuantity = isotopicquantity;
+
+ return *this;
+
+}
+
+//________________________________________________________________________
+IsotopicVector& IsotopicVector::operator*=(const double& factor)
+{
+
+ Multiply(factor);
+
+ return *this;
+
+}
+
+
+//________________________________________________________________________
+bool IsotopicVector::operator<(const IsotopicVector& isotopicvector) const
+{
+
+ if( Norme(*this) != Norme(isotopicvector) )
+ return Norme(*this) < Norme(isotopicvector);
+ else if( (*this).GetIsotopicQuantity().size() != isotopicvector.GetIsotopicQuantity().size() )
+ return (*this).GetIsotopicQuantity().size() < isotopicvector.GetIsotopicQuantity().size();
+ else
+ {
+ map<ZAI ,double>::iterator it;
+ map<ZAI ,double>::iterator it2 = isotopicvector.GetIsotopicQuantity().begin();
+ map<ZAI ,double> IsotopicQuantity = (*this).GetIsotopicQuantity();
+ for( it = IsotopicQuantity.begin(); it != IsotopicQuantity.end(); it++ )
+ {
+ if( (*it).first != (*it2).first )
+ return (*it).first < (*it2).first;
+ else it2++;
+ }
+ return false;
+ }
+
+}
+
+
+//________________________________________________________________________
+//________________________Constructor & Destructor________________________
+//________________________________________________________________________
+IsotopicVector::IsotopicVector()
+{
+}
+
+
+//_____________________________________________________GetSpeciesComposition___________________
+IsotopicVector::~IsotopicVector()
+{
+ fIsotopicQuantity.clear();
+ fIsotopicQuantityNeeded.clear();
+}
+
+
+
+//________________________________________________________________________
+//_____________________________General Method_____________________________
+//________________________________________________________________________
+void IsotopicVector::Clear()
+{
+
+ fIsotopicQuantityNeeded.clear();
+ fIsotopicQuantity.clear();
+
+}
+//________________________________________________________________________
+void IsotopicVector::ClearNeed()
+{
+
+ fIsotopicQuantityNeeded.clear();
+
+}
+
+//________________________________________________________________________
+void IsotopicVector::Multiply(double factor)
+{
+
+ map<ZAI ,double >::iterator it;
+ for( it = fIsotopicQuantity.begin(); it != fIsotopicQuantity.end(); it++)
+ (*it).second = (*it).second * factor;
+ for( it = fIsotopicQuantityNeeded.begin(); it != fIsotopicQuantityNeeded.end(); it++)
+ (*it).second = (*it).second * factor;
+
+
+}
+
+//________________________________________________________________________
+
+double IsotopicVector::GetSumOfAll() const
+{
+ double Sum = 0;
+ map<ZAI ,double >::iterator it;
+ map<ZAI ,double > isotopicquantity = GetIsotopicQuantity();
+ for( it = isotopicquantity.begin(); it != isotopicquantity.end(); it++)
+ Sum += (*it).second;
+
+ return Sum;
+
+}
+//________________________________________________________________________
+void IsotopicVector::Add(const ZAI& zai, double quantity)
+{
+
+ if( ceil(quantity*1e50) - quantity*1e50 > quantity*1e50 - floor(quantity*1e50) )
+ quantity = floor(quantity*1e50)*1/1e50;
+ else quantity = ceil(quantity*1e50)*1/1e50;
+
+
+ if(quantity > 0)
+ {
+ pair<map<ZAI, double>::iterator, bool> IResult;
+ IResult = fIsotopicQuantity.insert( pair<ZAI ,double>(zai, quantity));
+ if(!IResult.second)
+ IResult.first->second += quantity;
+ }
+
+
+}
+//________________________________________________________________________
+
+void IsotopicVector::Add(const IsotopicVector& isotopicvector)
+{
+
+ map<ZAI ,double> isotopicquantity = isotopicvector.GetIsotopicQuantity();
+ map<ZAI ,double >::iterator it;
+ for( it = isotopicquantity.begin(); it != isotopicquantity.end(); it++)
+ Add( (*it).first, (*it).second);
+
+
+}
+//________________________________________________________________________
+
+void IsotopicVector::Add(const map<ZAI ,double>& quantity)
+{
+
+ map<ZAI ,double> isotopicquantity = quantity;
+ map<ZAI ,double >::iterator it;
+ for( it = isotopicquantity.begin(); it != isotopicquantity.end(); it++)
+ Add( (*it).first, (*it).second);
+
+
+}
+
+
+//________________________________________________________________________
+void IsotopicVector::Remove(const ZAI& zai, double quantity)
+{
+
+
+ map<ZAI ,double>::iterator it;
+ it = fIsotopicQuantity.find(zai);
+
+ if(quantity > 0)
+ {
+ if ( it != fIsotopicQuantity.end() )
+ {
+ if (it->second > quantity)
+ it->second = it->second - quantity;
+ else
+ {
+ if( (it->second - quantity)/it->second > 1e-16 ) // to fit with double precision : 16 digits
+ Need(zai, quantity - it->second );
+ it->second = 0;
+ }
+ }
+ else
+ {
+ Need(zai, quantity);
+ }
+ }
+
+ if(it->second == 0)
+ fIsotopicQuantity.erase(it);
+}
+
+//________________________________________________________________________
+void IsotopicVector::Remove(const IsotopicVector& isotopicvector)
+{
+
+ map<ZAI ,double> isotopicquantity = isotopicvector.GetIsotopicQuantity();
+ map<ZAI ,double >::iterator it;
+ for( it = isotopicquantity.begin(); it != isotopicquantity.end(); it++)
+ Remove( (*it).first, (*it).second);
+
+}
+
+//________________________________________________________________________
+void IsotopicVector::Need(const ZAI& zai, double quantity)
+{
+ if(quantity < 0.5) quantity = 0;
+
+ pair<map<ZAI, double>::iterator, bool> IResult;
+ if(quantity > 0)
+ { cout << "Negative quantity : "<<quantity <<" for ZAI " << zai.Z() << " " << zai.A() << " " << zai.I() << " in this IsotopicVector" << endl;
+ exit(0);
+ IResult = fIsotopicQuantityNeeded.insert( pair<ZAI ,double>(zai, quantity));
+ if(!IResult.second)
+ IResult.first->second += quantity;
+ }
+
+
+}
+
+//________________________________________________________________________
+void IsotopicVector::Need(const IsotopicVector& isotopicvector)
+{
+
+ map<ZAI ,double> isotopicquantity = isotopicvector.GetIsotopicQuantity();
+ map<ZAI ,double >::iterator it;
+ for( it = isotopicquantity.begin(); it != isotopicquantity.end(); it++)
+ Need( (*it).first, (*it).second);
+
+}
+
+
+//________________________________________________________________________
+double IsotopicVector::GetZAIIsotopicQuantity(const ZAI& zai) const
+{
+
+ map<ZAI ,double> IsotopicQuantity = fIsotopicQuantity;
+
+ map<ZAI ,double>::iterator it;
+ it = IsotopicQuantity.find(zai);
+
+
+ if ( it != IsotopicQuantity.end() )
+ {
+ return it->second;
+ }
+ else
+ {
+ return 0;
+ }
+}
+
+//________________________________________________________________________
+double IsotopicVector::GetZAIIsotopicQuantity(const int z, const int a, const int i) const
+{
+
+ ZAI zai(z, a, i);
+ return GetZAIIsotopicQuantity(zai);
+}
+
+IsotopicVector IsotopicVector::GetSpeciesComposition(int z) const
+{
+
+ IsotopicVector IV;
+ map<ZAI ,double > IsotopicQuantity = GetIsotopicQuantity();
+ map<ZAI ,double >::iterator it;
+ for( it = IsotopicQuantity.begin(); it != IsotopicQuantity.end(); it++)
+ if( (*it).first.Z() == z )
+ IV += (*it).first * (*it).second;
+
+ return IV;
+
+}
+
+
+IsotopicVector IsotopicVector::GetThisComposition(IsotopicVector IV) const
+{
+ IsotopicVector IVtmp;
+ map<ZAI ,double > IsotopicQuantity = IV.GetIsotopicQuantity();
+ map<ZAI ,double > THISIsotopicQuantity = (*this).GetIsotopicQuantity();
+
+ map<ZAI ,double >::iterator it;
+ for( it = IsotopicQuantity.begin(); it != IsotopicQuantity.end(); it++)
+ {
+ map<ZAI ,double >::iterator it2 = THISIsotopicQuantity.find( (*it).first );
+ if(it2 != THISIsotopicQuantity.end())
+ IVtmp += (*it2).first * (*it2).second ;
+
+ }
+
+ return IVtmp;
+
+}
+//________________________________________________________________________
+double IsotopicVector::GetTotalMass() const
+{
+ return cZAIMass.GetMass(*this);//in tons
+}
+//________________________________________________________________________
+
+double IsotopicVector::MeanMolar() const
+{
+ return GetTotalMass() * 1e6 * AVOGADRO / GetActinidesComposition().GetSumOfAll();
+}
+//________________________________________________________________________
+
+vector<ZAI> IsotopicVector::GetZAIList() const
+{
+
+ map<ZAI ,double > IsotopicQuantity = GetIsotopicQuantity();
+ map<ZAI ,double >::iterator it;
+ vector<ZAI> zailist;
+ for( it = IsotopicQuantity.begin(); it != IsotopicQuantity.end(); it++)
+ zailist.push_back( (*it).first );
+
+ return zailist;
+
+}
+
+IsotopicVector IsotopicVector::GetActinidesComposition() const
+{
+
+ IsotopicVector IV;
+ for (int i = 0; i <12; i++)
+ IV += GetSpeciesComposition(89+i);
+ return IV;
+
+}
+
+vector<int> IsotopicVector::GetChemicalSpecies() const
+{
+
+ vector<int> ChemicalSpecies;
+
+ map<ZAI ,double > IsotopicQuantity = GetIsotopicQuantity();
+ map<ZAI ,double >::iterator it;
+ for( it = IsotopicQuantity.begin(); it != IsotopicQuantity.end(); it++)
+ if( (int)ChemicalSpecies.size() ==0 || (*it).first.Z() != ChemicalSpecies.back() )
+ ChemicalSpecies.push_back((*it).first.Z());
+
+
+ return ChemicalSpecies;
+}
+
+
+//________________________________________________________________________
+void IsotopicVector::Write(string filename, cSecond time) const
+{
+ ofstream IVfile(filename.c_str(), ios_base::app); // Open the File
+ if(!IVfile)
+ cout << "!!Warning!! !!!IsotopicVector!!! \n Can't open \"" << filename << "\"\n" << endl;
+
+ if(time != -1)
+ IVfile << "Time "<< time/365.25/3600./24. << endl;
+
+ map<ZAI ,double> IsotopicQuantity = GetIsotopicQuantity();
+ map<ZAI ,double >::iterator it;
+ for(it = IsotopicQuantity.begin(); it != IsotopicQuantity.end(); it++)
+ {
+ IVfile << (*it).first.Z() << " ";
+ IVfile << (*it).first.A() << " ";
+ IVfile << (*it).first.I() << " ";
+ IVfile << (*it).second << " " << endl;
+ }
+ IVfile << endl;
+}
+//________________________________________________________________________
+void IsotopicVector::Print(string option) const
+{
+
+ cout<<sPrint();
+
+}
+//________________________________________________________________________
+string IsotopicVector::sPrint() const
+{
+ stringstream ss;
+ ss << "**************************" << endl;
+ ss << "*Isotopic Vector Property*" << endl;
+ ss << "**************************" << endl << endl;
+
+ bool QuantityPrint = false;
+ bool DBPrint = false;
+
+ QuantityPrint = true;
+
+ if(QuantityPrint)
+ {
+ ss << "*Isotopic Vector Quantity*" << endl;
+ map<ZAI ,double> IsotopicQuantity = GetIsotopicQuantity();
+ map<ZAI ,double >::iterator it;
+ for(it = IsotopicQuantity.begin();it != IsotopicQuantity.end(); it++)
+ {
+ ss << (*it).first.Z() << " ";
+ ss << (*it).first.A() << " ";
+ ss << (*it).first.I() << " ";
+ ss << ": " << (*it).second;
+ ss << endl;
+ }
+ ss << endl;
+ ss << "*Isotopic Vector Quantity Needed*" << endl;
+ map<ZAI ,double> IsotopicQuantityNeeded = GetIsotopicQuantityNeeded();
+ for(it = IsotopicQuantityNeeded.begin(); it != IsotopicQuantityNeeded.end(); it++)
+ {
+ ss << (*it).first.Z() << " ";
+ ss << (*it).first.A() << " ";
+ ss << (*it).first.I() << " ";
+ ss << ": " << (*it).second;
+ ss << endl;
+ }
+ ss << endl;
+ }
+ if(DBPrint)
+ {
+ ss << "****Isotopic Vector DB****" << endl;
+ }
+return ss.str();
+}
+//________________________________________________________________________
+void IsotopicVector::PrintList(string option) const
+{
+ bool QuantityPrint = false;
+ bool DBPrint = false;
+
+ QuantityPrint = true;
+
+ if(QuantityPrint)
+ {
+ map<ZAI ,double> IsotopicQuantity = GetIsotopicQuantity();
+ map<ZAI ,double >::iterator it;
+ for(it = IsotopicQuantity.begin();it != IsotopicQuantity.end(); it++)
+ {
+ cout << (*it).first.Z() << " ";
+ cout << (*it).first.A() << " ";
+ cout << (*it).first.I() << " ";
+ cout << endl;
+ }
+ cout << endl;
+
+ }
+ if(DBPrint)
+ {
+ cout << "****Isotopic Vector DB****" << endl;
+ }
+
+}
+
+
+
+
diff --git a/source/trunk/src/Makefile b/source/trunk/src/Makefile
new file mode 100755
index 000000000..8c2be94fe
--- /dev/null
+++ b/source/trunk/src/Makefile
@@ -0,0 +1,125 @@
+# Directory containing libraries
+LIBDIR = $(CLASS_lib)
+# Directory containing includes for CLASS
+LOCALINC = $(CLASS_include)
+
+
+EQM = ../Model/Equivalence
+IM = ../Model/Irradiation
+XSM = ../Model/XS
+
+
+ROOTCFLAGS = `root-config --cflags`
+ROOTLIBS = `root-config --libs`
+ROOTGLIBS = `root-config --glibs`
+
+######### nothing to change from here #########
+INCLUDES = $(LOCALINC)/*.hxx
+LIBNAME = CLASSpkg
+OBJS = CLASSLogger.o \
+ ZAI.o ZAIDict.o \
+ IsotopicVector.o IsotopicVectorDict.o \
+ ZAIMass.o \
+ CLASSNucleiFiliation.o \
+ CLASSObject.o CLASSObjectDict.o\
+ CLASSFacility.o CLASSFacilityDict.o\
+ FabricationPlant.o FabricationPlantDict.o \
+ Reactor.o ReactorDict.o \
+ CLASSBackEnd.o CLASSBackEndDict.o\
+ SeparationPlant.o SeparationPlantDict.o\
+ Storage.o StorageDict.o\
+ Pool.o PoolDict.o\
+ DecayDataBank.o \
+ DynamicalSystem.o\
+ IrradiationModel.o \
+ EquivalenceModel.o \
+ XSModel.o \
+ CLASSFuel.o\
+ PhysicsModels.o \
+ EvolutionData.o EvolutionDataDict.o \
+ CLASSFuelPlan.o\
+ Scenario.o
+
+OBJMODEL = $(EQM)/EQM_MLP_PWR_MOX.o $(EQM)/EQM_QUAD_PWR_MOX.o $(EQM)/EQM_LIN_PWR_MOX.o \
+ $(XSM)/XSM_MLP.o $(XSM)/XSM_CLOSEST.o \
+ $(IM)/IM_RK4.o $(IM)/IM_Matrix.o
+
+
+ROOTOBJS = CLASSLogger.o \
+ ZAI.o ZAIDict.o \
+ IsotopicVector.o IsotopicVectorDict.o \
+ ZAIMass.o \
+ CLASSObject.o CLASSObjectDict.o\
+ CLASSFacility.o CLASSFacilityDict.o\
+ Reactor.o ReactorDict.o \
+ FabricationPlant.o FabricationPlantDict.o \
+ CLASSBackEnd.o CLASSBackEndDict.o\
+ Storage.o StorageDict.o\
+ Pool.o PoolDict.o\
+ SeparationPlant.o SeparationPlantDict.o\
+ DecayDataBank.o \
+ DynamicalSystem.o\
+ CLASSFuel.o\
+ EvolutionData.o EvolutionDataDict.o \
+ CLASSFuelPlan.o\
+ PhysicsModels.o
+
+
+CXX = g++
+CXXFLAGS = -O2 -g -Wall -fopenmp -fPIC -I$(LOCALINC) $(ROOTCFLAGS)
+LD = g++
+LDFLAGS = -g -Wall -fPIC $(ROOTLIBS) -lTMVA -shared -lgomp
+
+all: $(OBJS) $(OBJMODEL)
+ $(LD) $(LDFLAGS) $(OBJS) $(OBJMODEL) -o $(LIBDIR)/lib$(LIBNAME).so
+ @echo "lib$(LIBNAME).so done"
+ $(LD) $(LDFLAGS) $(ROOTOBJS) -o $(LIBDIR)/lib$(LIBNAME)_root.so
+ @echo "lib$(LIBNAME)_root.so done"
+
+clean:
+ @rm -vf $(OBJS) $(OBJMODEL) *~ core *Dict.cxx *Dict.h
+
+
+install:
+ @ln -sf ../Model/* ../include/
+
+
+
+CLASSObjectDict.cxx: $(LOCALINC)/CLASSObject.hxx
+ rootcint -f $@ -c $^
+
+CLASSFacilityDict.cxx: $(LOCALINC)/CLASSFacility.hxx
+ rootcint -f $@ -c $^
+
+CLASSBackEndDict.cxx: $(LOCALINC)/CLASSBackEnd.hxx
+ rootcint -f $@ -c $^
+
+StorageDict.cxx: $(LOCALINC)/Storage.hxx
+ rootcint -f $@ -c $^
+
+ReactorDict.cxx: $(LOCALINC)/Reactor.hxx
+ rootcint -f $@ -c $^
+
+FabricationPlantDict.cxx: $(LOCALINC)/FabricationPlant.hxx
+ rootcint -f $@ -c $^
+
+PoolDict.cxx: $(LOCALINC)/Pool.hxx
+ rootcint -f $@ -c $^
+
+SeparationPlantDict.cxx: $(LOCALINC)/SeparationPlant.hxx
+ rootcint -f $@ -c $^
+
+IsotopicVectorDict.cxx: $(LOCALINC)/IsotopicVector.hxx
+ rootcint -f $@ -c $^
+
+ZAIDict.cxx: $(LOCALINC)/ZAI.hxx
+ rootcint -f $@ -c $^
+
+EvolutionDataDict.cxx: $(LOCALINC)/EvolutionData.hxx
+ rootcint -f $@ -c $^
+
+
+.SUFFIXES: .cxx
+
+%.o: %.cxx $(INCLUDES)
+ $(CXX) $(CXXFLAGS) -c $*.cxx -o $*.o
diff --git a/source/trunk/src/PhysicsModels.cxx b/source/trunk/src/PhysicsModels.cxx
new file mode 100644
index 000000000..ab9f23061
--- /dev/null
+++ b/source/trunk/src/PhysicsModels.cxx
@@ -0,0 +1,57 @@
+#include "PhysicsModels.hxx"
+
+//________________________________________________________________________
+//
+// PhysicsModels
+//
+//
+//
+//
+//________________________________________________________________________
+
+
+
+PhysicsModels::PhysicsModels():CLASSObject()
+{
+
+ fXSModel = 0;
+ fEquivalenceModel = 0;
+ fIrradiationModel = 0;
+
+
+}
+//________________________________________________________________________
+PhysicsModels::PhysicsModels(XSModel* XS, EquivalenceModel* EM, IrradiationModel* IM ):CLASSObject()
+{
+
+ fXSModel = XS;
+ fEquivalenceModel = EM;
+ fIrradiationModel = IM;
+
+ int Z_ZAIThreshold = fIrradiationModel->GetZAIThreshold();
+ fXSModel->SetZAIThreshold(Z_ZAIThreshold);
+
+
+}
+//________________________________________________________________________
+PhysicsModels::PhysicsModels(CLASSLogger* log, XSModel* XS, EquivalenceModel* EM, IrradiationModel* IM ):CLASSObject(log)
+{
+
+ fXSModel = XS;
+ fEquivalenceModel = EM;
+ fIrradiationModel = IM;
+
+ int Z_ZAIThreshold = fIrradiationModel->GetZAIThreshold();
+ fXSModel->SetZAIThreshold(Z_ZAIThreshold);
+
+
+}
+
+//________________________________________________________________________
+EvolutionData PhysicsModels::GenerateEvolutionData(IsotopicVector IV, double cycletime, double Power)
+{
+ fXSModel->isIVInDomain(IV);
+
+ return fIrradiationModel->GenerateEvolutionData(IV, fXSModel->GetCrossSections(IV), Power, cycletime);
+}
+//________________________________________________________________________
diff --git a/source/trunk/src/Pool.cxx b/source/trunk/src/Pool.cxx
new file mode 100755
index 000000000..cbbb3cc0b
--- /dev/null
+++ b/source/trunk/src/Pool.cxx
@@ -0,0 +1,226 @@
+#include "Pool.hxx"
+
+#include "IsotopicVector.hxx"
+#include "Storage.hxx"
+#include "Scenario.hxx"
+#include "CLASSLogger.hxx"
+
+#include <sstream>
+#include <string>
+#include <iostream>
+#include <cmath>
+#include <algorithm>
+
+//________________________________________________________________________
+//
+// Pool
+//
+//
+//
+//
+//________________________________________________________________________
+ClassImp(Pool)
+
+
+Pool::Pool():CLASSBackEnd(8)
+{
+ fOutBackEndFacility = 0;
+ SetName("P_Pool.");
+}
+
+ //________________________________________________________________________
+Pool::Pool(CLASSLogger* log, cSecond coolingtime):CLASSBackEnd(log, coolingtime, 8)
+{
+ DBGL
+
+
+ fCycleTime = (cSecond)coolingtime;
+ fPutToWaste = true;
+ fCoolingLastIndex = 0;
+
+ fOutBackEndFacility = 0;
+ SetName("P_Pool.");
+
+
+ INFO << " A new Pool has been define :" << endl;
+ INFO << "\t The Cooling Time set at\t " << (double)(fCycleTime/3600/24/365.25) << " year" << endl;
+ WARNING << " All Cooled Fuel goes directly to WASTE after cooling !! " << endl;
+
+ DBGL
+
+}
+
+//________________________________________________________________________
+Pool::Pool(CLASSLogger* log, CLASSBackEnd* storage, cSecond coolingtime):CLASSBackEnd(log, coolingtime, 8)
+{
+ DBGL
+
+ fOutBackEndFacility = storage;
+ SetIsStorageType(false);
+
+ fPutToWaste = false;
+ fCoolingLastIndex = 0;
+ SetName("P_Pool.");
+
+
+
+ INFO << " A new Pool has been define :" << endl;
+ INFO << "\t The Cooling Time set at\t " << (double)(fCycleTime/3600/24/365.25) << " year" << endl;
+
+ DBGL
+
+}
+
+//________________________________________________________________________
+Pool::~Pool()
+{
+
+
+}
+
+//________________________________________________________________________
+//________________________________________________________________________
+void Pool::SetIVArray(vector<IsotopicVector> ivarray)
+{
+ INFO << "This method as no effect !!!" << endl;
+ INFO << "Use SetIVArray(vector<IsotopicVector> ivarray, vector<cSecond>n timearray) unstead!!!!"<<endl;
+}
+
+
+//________________________________________________________________________
+void Pool::SetIVArray(vector<IsotopicVector> ivarray, vector<cSecond> timearray)
+{
+ fIVArray = ivarray;
+ fIVArrayArrivalTime = timearray;
+
+}
+//________________________________________________________________________
+// Add Temporary IV :
+// Cooling
+//
+//________________________________________________________________________
+void Pool::AddIV(IsotopicVector IV)
+{
+
+ fIVArray.push_back(IV);
+ fInsideIV += IV;
+ fIVArrayArrivalTime.push_back(fInternalTime);
+ fCoolingLastIndex++;
+ fCoolingIndex.push_back(fCoolingLastIndex);
+
+ AddCumulativeIVIn(IV);
+
+}
+
+
+//________________________________________________________________________
+void Pool::RemoveIVCooling(int i) //!< Remove a Cooling IsotopicVector
+{
+ AddCumulativeIVOut(fIVArray[i]);
+
+ fIVArray.erase(fIVArray.begin()+i);
+ fIVArrayArrivalTime.erase( fIVArrayArrivalTime.begin()+i);
+ fCoolingIndex.erase(fCoolingIndex.begin()+i);
+ UpdateInsideIV();
+}
+
+
+
+
+//________________________________________________________________________
+// Time Action with the reste of the Universe :
+// Out Storage
+// Evolution :
+// Cooling
+//________________________________________________________________________
+
+
+
+
+//________________________________________________________________________
+void Pool::CoolingEvolution(cSecond t)
+{
+DBGL
+
+ if(t == fInternalTime && t!=0) return;
+ int RemainingCoolingTime;
+ cSecond EvolutionTime = t - fInternalTime;
+
+#pragma omp parallel for
+ for ( int i = 0 ; i < (int)fIVArray.size() ; i++)
+ {
+ if ( abs(t - fIVArrayArrivalTime[i] - fCycleTime) < 3600 ) // ">" should not append, only "=" is normal...
+ {
+ if (t - fIVArrayArrivalTime[i] > fCycleTime) // Warning & Quit
+ {
+ ERROR << " Cooling Step : " << t/3600./24/365.25<< " :" << " an evolution Step is probably missing ! " << " " << endl;
+ exit (1);
+ }
+
+ RemainingCoolingTime = fCycleTime - (fInternalTime - fIVArrayArrivalTime[i]);
+ //Cooling Decay
+ fIVArray[i] = GetDecay( fIVArray[i], RemainingCoolingTime);
+
+
+#pragma omp critical(DeleteCoolingIVPB)
+ {fCoolingEndOfCycle.push_back(i);}
+
+ }
+ else if ( fIVArrayArrivalTime[i] != t )
+ {
+ fIVArray[i] = GetDecay( fIVArray[i] , EvolutionTime);
+ }
+ }
+#pragma omp critical(DeleteCoolingIVPB)
+ {sort (fCoolingEndOfCycle.begin(), fCoolingEndOfCycle.end());}
+
+DBGL
+}
+
+
+//________________________________________________________________________
+void Pool::Evolution(cSecond t)
+{
+
+ // Check if the Pool has been created ...
+ if(t == fInternalTime && t!=0) return;
+ // Make the evolution for the Cooling IV ...
+ CoolingEvolution(t);
+ // Update Inside IV
+ UpdateInsideIV();
+ // ... And Finaly update the AbsoluteInternalTime
+ fInternalTime = t;
+
+
+
+}
+
+
+//________________________________________________________________________
+void Pool::Dump()
+{
+DBGL
+//------ Cooling ------//
+ for(int i = (int)fCoolingEndOfCycle.size()-1; i >=0 ; i--) // IV End Of Cooling
+ {
+
+ int idx = fCoolingEndOfCycle[i]; // Get Index number
+
+ if(!fPutToWaste)
+ fOutBackEndFacility->AddIV(fIVArray[idx]);
+ else
+ GetParc()->AddWaste(fIVArray[idx]);
+
+ fCoolingEndOfCycle.erase(fCoolingEndOfCycle.begin()+i); // Remove index entry
+ RemoveIVCooling(idx); // Remove IVcooling
+ }
+
+ if((int)fCoolingEndOfCycle.size() != 0 )// Control
+ {
+ ERROR << "Problem while Dumping Cooling"<< endl;
+ exit (1);
+ }
+DBGL
+}
+
+
diff --git a/source/trunk/src/Reactor.cxx b/source/trunk/src/Reactor.cxx
new file mode 100755
index 000000000..adfc85a24
--- /dev/null
+++ b/source/trunk/src/Reactor.cxx
@@ -0,0 +1,607 @@
+#include "Reactor.hxx"
+
+#include "EvolutionData.hxx"
+#include "Pool.hxx"
+#include "FabricationPlant.hxx"
+#include "Storage.hxx"
+#include "Scenario.hxx"
+#include "CLASSConstante.hxx"
+
+#include <iostream>
+#include <cmath>
+#include <omp.h>
+#include <typeinfo>
+
+//________________________________________________________________________
+//
+// Reactor
+//
+//
+//
+//
+//________________________________________________________________________
+
+
+
+ClassImp(Reactor)
+
+
+Reactor::Reactor():CLASSFacility(4)
+{
+
+ SetName("R_Reactor.");
+
+ fOutBackEndFacility = 0;
+ fStorage = 0;
+ fFabricationPlant = 0;
+
+}
+
+Reactor::Reactor(CLASSLogger* log):CLASSFacility(log, 4)
+{
+ DBGL
+
+ fOutBackEndFacility = 0;
+ fStorage = 0;
+ fFabricationPlant = 0;
+ SetName("R_Reactor.");
+
+ DBGL
+}
+
+Reactor::Reactor(CLASSLogger* log,
+ CLASSBackEnd* Pool,
+ cSecond creationtime,
+ cSecond lifetime,
+ double power, double HMMass, double CapacityFactor ):CLASSFacility(log, creationtime, lifetime, 4)
+{
+ DBGL
+ (*this).SetName("R_Reactor.");
+
+
+ fIsStarted = false;
+ fIsShutDown = false;
+ fIsAtEndOfCycle = false;
+
+ fStorage = 0;
+ fFabricationPlant = 0;
+
+ fFixedFuel = true;
+ fIsStorage = false;
+
+ fOutBackEndFacility = Pool;
+
+ fPower = power * CapacityFactor;
+
+ fHeavyMetalMass = HMMass;
+
+ fBurnUp = -1;
+ fCycleTime = (-1);
+
+ fIVBeginCycle = fEvolutionDB.GetIsotopicVectorAt(0);
+ fIVInCycle = fEvolutionDB.GetIsotopicVectorAt(0);
+ fIVOutCycle = fEvolutionDB.GetIsotopicVectorAt( (cSecond)(fCycleTime/fEvolutionDB.GetPower()*fPower) );
+
+
+ fFuelPlan = 0;
+
+ INFO << " A Reactor has been define :" << endl;
+ INFO << "\t Fuel Composition is fixed (for now)! "<< endl;
+ INFO << "\t Creation time set at \t " << (double)(GetCreationTime()/3600/24/365.25) << " year" << endl;
+ INFO << "\t Life time (Operating's Duration) set at \t " << (double)(GetLifeTime()/3600/24/365.25) << " year" << endl;
+ INFO << "\t The Effective Thermal Power is \t " << (double)(fPower *1e-6) << " MW (with Full Power " << power << " and " << CapacityFactor << " capacity factor)"<< endl;
+ INFO << "\t The Heavy Metal Mass in the Core set at " << (double)(fHeavyMetalMass) << " tons" << endl << endl;
+
+ DBGL
+
+}
+
+Reactor::Reactor(CLASSLogger* log,
+ FabricationPlant* fabricationplant, CLASSBackEnd* Pool,
+ cSecond creationtime, cSecond lifetime,
+ double Power, double HMMass, double CapacityFactor):CLASSFacility(log, creationtime, lifetime, 4)
+{
+ DBGL
+ (*this).SetName("R_Reactor.");
+
+
+ fStorage = 0;
+ fIsStarted = false;
+ fIsShutDown = false;
+ fIsAtEndOfCycle = false;
+
+ fFabricationPlant = fabricationplant;
+ fFixedFuel = false;
+
+ fOutBackEndFacility = Pool;
+
+ fBurnUp = -1;
+ fHeavyMetalMass = HMMass;
+ fPower = Power*CapacityFactor;
+ fCycleTime = -1; //BU in GWd/t
+
+ fFuelPlan = 0;
+
+
+
+ INFO << " A Reactor has been define :" << endl;
+ INFO << "\t Fuel Composition is not fixed (for now)! "<< endl;
+ INFO << "\t Creation time set at \t " << (double)(GetCreationTime()/3600/24/365.25) << " year" << endl;
+ INFO << "\t Life time (Operating's Duration) set at \t " << (double)(GetLifeTime()/3600/24/365.25) << " year" << endl;
+ INFO << "\t The Effective Thermal Power is \t " << (double)(fPower *1e-6) << " MW (with Full Power " << Power << " and " << CapacityFactor << " capacity factor)"<< endl;
+ INFO << "\t The Heavy Metal Mass in the Core set at " << (double)(fHeavyMetalMass) << " tons" << endl << endl;
+
+
+ DBGL
+
+}
+
+
+Reactor::Reactor(CLASSLogger* log, PhysicsModels* fueltypeDB, FabricationPlant* fabricationplant, CLASSBackEnd* Pool,
+ cSecond creationtime, cSecond lifetime,
+ double Power, double HMMass, double BurnUp, double CapacityFactor):CLASSFacility(log, creationtime, lifetime, 4)
+{
+ DBGL
+ (*this).SetName("R_Reactor.");
+
+
+ fStorage = 0;
+ fIsStarted = false;
+ fIsShutDown = false;
+ fIsAtEndOfCycle = false;
+
+ fFabricationPlant = fabricationplant;
+
+ fFixedFuel = false;
+
+ fOutBackEndFacility = Pool;
+
+ fBurnUp = BurnUp;
+ fHeavyMetalMass = HMMass;
+ fPower = Power*CapacityFactor;
+ fCycleTime = (cSecond) (fBurnUp*1e9 / (fPower) * fHeavyMetalMass *3600*24); //BU in GWd/t
+
+ fFuelPlan = new CLASSFuelPlan(log);
+ fFuelPlan->AddFuel(creationtime, CLASSFuel(fueltypeDB), fBurnUp);
+
+
+
+ INFO << " A Reactor has been define :" << endl;
+ INFO << "\t Fuel Composition is not fixed ! "<< endl;
+ INFO << "\t Creation time set at \t " << (double)(GetCreationTime()/3600/24/365.25) << " year" << endl;
+ INFO << "\t Life time (Operating's Duration) set at \t " << (double)(GetLifeTime()/3600/24/365.25) << " year" << endl;
+ INFO << "\t The Effective Thermal Power is \t " << (double)(fPower *1e-6) << " MW (with Full Power " << Power << " and " << CapacityFactor << " capacity factor)"<< endl;
+ INFO << "\t Burn-Up at end of Cycle set at \t " << (double)(fBurnUp) << " GWj/t" << endl;
+ INFO << "\t The corresponding Cycle Time is\t " << (double)(fCycleTime/3600/24/365.25) << " year" << endl;
+ INFO << "\t The Heavy Metal Mass in the Core set at " << (double)(fHeavyMetalMass) << " tons" << endl << endl;
+
+
+ DBGL
+
+}
+
+Reactor::Reactor(CLASSLogger* log, PhysicsModels* fueltypeDB,
+ FabricationPlant* fabricationplant,
+ CLASSBackEnd* Pool,
+ cSecond creationtime, cSecond lifetime, cSecond cycletime,
+ double HMMass, double BurnUp):CLASSFacility(log, creationtime, lifetime, cycletime, 4)
+{
+ DBGL
+ (*this).SetName("R_Reactor.");
+
+
+ fIsStarted = false;
+ fIsShutDown = false;
+ fIsAtEndOfCycle = false;
+
+ fStorage = 0;
+
+ fFabricationPlant = fabricationplant;
+ fFixedFuel = false;
+ fBurnUp = BurnUp;
+ fHeavyMetalMass = HMMass;
+
+ fOutBackEndFacility = Pool;
+ fPower = BurnUp*3600.*24. / (fCycleTime) * HMMass *1e9; //BU in GWd/t
+
+ fFuelPlan = new CLASSFuelPlan(log);
+ fFuelPlan->AddFuel(creationtime, CLASSFuel(fueltypeDB), fBurnUp);
+
+
+ INFO << " A Reactor has been define :" << endl;
+ INFO << "\t Fuel Composition is not fixed ! "<< endl;
+ INFO << "\t Creation time set at \t " << (double)(GetCreationTime()/3600/24/365.25) << " year" << endl;
+ INFO << "\t Life time (Operating's Duration) set at \t " << (double)(GetCreationTime()/3600/24/365.25) << " year" << endl;
+ INFO << "\t The Cycle Time set at\t " << (double)(fCycleTime/3600/24/365.25) << " year" << endl;
+ INFO << "\t Burn-Up at end of Cycle set at \t " << (double)(fBurnUp) << " GWj/t" << endl;
+ INFO << "\t The corresponding Effective Thermal Power is \t " << (double)(fPower *1e-6) << " MW" << endl;
+ INFO << "\t The Heavy Metal Mass in the Core set at " << (double)(fHeavyMetalMass) << " tons" << endl << endl;
+
+
+ DBGL
+
+}
+
+
+Reactor::Reactor(CLASSLogger* log, EvolutionData* evolutivedb,
+ CLASSBackEnd* Pool,
+ cSecond creationtime,
+ cSecond lifetime,
+ double power, double HMMass, double BurnUp, double CapacityFactor ):CLASSFacility(log, creationtime, lifetime, 4)
+{
+ DBGL
+ (*this).SetName("R_Reactor.");
+
+
+ fIsStarted = false;
+ fIsShutDown = false;
+ fIsAtEndOfCycle = false;
+
+ fStorage = 0;
+ fFabricationPlant = 0;
+
+ fFixedFuel = true;
+ fIsStorage = false;
+
+ fOutBackEndFacility = Pool;
+
+ fPower = power * CapacityFactor;
+
+ fHeavyMetalMass = HMMass;
+
+ double M0 = cZAIMass.GetMass( evolutivedb->GetIsotopicVectorAt(0.).GetActinidesComposition() );
+
+ fEvolutionDB = (*evolutivedb) * (fHeavyMetalMass/M0);
+
+ fBurnUp = BurnUp;
+ fCycleTime = (cSecond) (fBurnUp*1e9 / (fPower) * fHeavyMetalMass *3600*24);
+
+ fIVBeginCycle = fEvolutionDB.GetIsotopicVectorAt(0);
+ fIVInCycle = fEvolutionDB.GetIsotopicVectorAt(0);
+ fIVOutCycle = fEvolutionDB.GetIsotopicVectorAt( (cSecond)(fCycleTime/fEvolutionDB.GetPower()*fPower) );
+
+ fFuelPlan = new CLASSFuelPlan(log);
+ fFuelPlan->AddFuel(creationtime, CLASSFuel(evolutivedb), fBurnUp);
+
+ INFO << " A Reactor has been define :" << endl;
+ INFO << "\t Fuel Composition is fixed ! "<< endl;
+ INFO << "\t Creation time set at \t " << (double)(GetCreationTime()/3600/24/365.25) << " year" << endl;
+ INFO << "\t Life time (Operating's Duration) set at \t " << (double)(GetLifeTime()/3600/24/365.25) << " year" << endl;
+ INFO << "\t The Cycle Time set at\t " << (double)(fCycleTime/3600/24/365.25) << " year" << endl;
+ INFO << "\t The Effective Thermal Power is \t " << (double)(fPower *1e-6) << " MW (with Full Power " << power << " and " << CapacityFactor << " capacity factor)"<< endl;
+ INFO << "\t The Heavy Metal Mass in the Core set at " << (double)(fHeavyMetalMass) << " tons" << endl << endl;
+
+
+ DBGL
+}
+
+Reactor::Reactor(CLASSLogger* log, EvolutionData* evolutivedb,
+ CLASSBackEnd* Pool,
+ cSecond creationtime, cSecond lifetime, cSecond cycletime,
+ double HMMass, double BurnUp ):CLASSFacility(log, creationtime, lifetime, cycletime, 4)
+{
+ DBGL
+ (*this).SetName("R_Reactor.");
+
+
+ fIsStarted = false;
+ fIsShutDown = false;
+ fIsAtEndOfCycle = false;
+
+ fStorage = 0;
+ fFabricationPlant = 0;
+
+ fFixedFuel = true;
+ fIsStorage = false;
+
+ fOutBackEndFacility = Pool;
+
+ fPower = BurnUp*3600.*24. / (fCycleTime) * HMMass *1e9; //BU in GWd/t
+
+ fHeavyMetalMass = HMMass;
+
+ double M0 = cZAIMass.GetMass( evolutivedb->GetIsotopicVectorAt(0.).GetActinidesComposition() );
+
+ fEvolutionDB = (*evolutivedb) * (fHeavyMetalMass/M0);
+
+ fBurnUp = BurnUp;
+
+ fIVBeginCycle = fEvolutionDB.GetIsotopicVectorAt(0);
+ fIVInCycle = fEvolutionDB.GetIsotopicVectorAt(0);
+ fIVOutCycle = fEvolutionDB.GetIsotopicVectorAt( (cSecond)(fCycleTime/fEvolutionDB.GetPower()*fPower) );
+
+
+ fFuelPlan->AddFuel(creationtime, CLASSFuel(evolutivedb), fBurnUp);
+
+ INFO << " A Reactor has been define :" << endl;
+ INFO << "\t Fuel Composition is fixed ! "<< endl;
+ INFO << "\t Creation time set at \t " << (double)(GetCreationTime()/3600/24/365.25) << " year" << endl;
+ INFO << "\t Life time (Operating's Duration) set at \t " << (double)(GetLifeTime()/3600/24/365.25) << " year" << endl;
+ INFO << "\t The Cycle Time set at\t " << (double)(fCycleTime/3600/24/365.25) << " year" << endl;
+ INFO << "\t The Effective Thermal Power is \t " << (double)(fPower *1e-6) << " MW (with Full Power " << fPower << endl;
+ INFO << "\t The Heavy Metal Mass in the Core set at " << (double)(fHeavyMetalMass) << " tons" << endl << endl;
+
+ DBGL
+}
+
+
+//________________________________________________________________________
+Reactor::~Reactor()
+{
+
+
+}
+
+//________________________________________________________________________
+void Reactor::SetCycleTime(double cycletime)
+{
+ fCycleTime = (cSecond)cycletime;
+
+ if(fFixedFuel==true)
+ {
+ fIVOutCycle = fEvolutionDB.GetIsotopicVectorAt(fCycleTime/fEvolutionDB.GetPower()*fPower);
+ fBurnUp = fPower*fCycleTime/3600./24./fHeavyMetalMass;
+ }
+ else
+ {
+ fBurnUp = fPower*fCycleTime/3600./24./fHeavyMetalMass;
+ }
+}
+//________________________________________________________________________
+void Reactor::SetPower(double Power)
+{
+ fPower = Power;
+
+ if(fFixedFuel==true)
+ {
+ fCycleTime = (cSecond) (fBurnUp*1e9 / (fPower) * fHeavyMetalMass *3600*24);
+ fIVOutCycle = fEvolutionDB.GetIsotopicVectorAt( (cSecond)(fCycleTime/fEvolutionDB.GetPower()*fPower) );
+ }
+ else
+ fCycleTime = (cSecond)(fBurnUp*1e9 / (fPower) * fHeavyMetalMass *3600*24); //BU in GWd/t
+
+
+}
+//________________________________________________________________________
+void Reactor::SetBurnUp(double BU)
+{
+
+ fBurnUp = BU;
+
+ if(fFixedFuel==true)
+ {
+ fCycleTime = (cSecond) (fBurnUp*1e9 / (fPower) * fHeavyMetalMass *3600*24);
+ fIVOutCycle = fEvolutionDB.GetIsotopicVectorAt( (cSecond)(fCycleTime/fEvolutionDB.GetPower()*fPower) );
+ }
+ else
+ fCycleTime = (cSecond) (fBurnUp*1e9 / (fPower) * fHeavyMetalMass *3600*24);
+
+}
+
+//________________________________________________________________________
+void Reactor::SetEvolutionDB(EvolutionData evolutionDB)
+{
+ double M0 = cZAIMass.GetMass( evolutionDB.GetIsotopicVectorAt(0.).GetActinidesComposition() );
+ fEvolutionDB = evolutionDB * (fHeavyMetalMass/M0);
+
+ fIVOutCycle = fEvolutionDB.GetIsotopicVectorAt( (cSecond)(fCycleTime/fEvolutionDB.GetPower()*fPower) );
+ fIVBeginCycle = fEvolutionDB.GetIsotopicVectorAt(0);
+
+
+
+}
+
+//________________________________________________________________________
+void Reactor::SetNewFuel(EvolutionData ivdb)
+{
+
+ SetEvolutionDB(ivdb);
+
+}
+
+//________________________________________________________________________
+void Reactor::Evolution(cSecond t)
+{
+DBGL
+
+
+ if( fIsShutDown || t < GetCreationTime() ) return; // Reactor stop or not started...
+
+ if(Norme(fInsideIV)!=0)
+ {
+#pragma omp critical(ParcPowerUpdate)
+ {GetParc()->AddToPower(fPower);}
+ }
+ else if(fIsStarted==true)
+ {
+ WARNING << " Reactor should be working but have no Heavy Nucleus Inside. It's not working so have a zero power..."
+ << " Time : "<< t/365.25/3600/24 << " years" << endl;
+ }
+
+
+ if( t < fInternalTime ) return;
+ if( t == fInternalTime && t!=GetCreationTime() ) return;
+
+
+
+ if( t == GetCreationTime() && !fIsStarted) // Start of the Reactor
+ {
+ fIsAtEndOfCycle = true;
+ fInsideIV = fIVBeginCycle;
+ fInternalTime = t;
+ fInCycleTime = 0;
+
+ }
+
+ // Check if the Reactor if started ...
+ if(!fIsStarted) return; // If the reactor just start don't need to make Fuel evolution
+
+
+ cSecond EvolutionTime = t - fInternalTime; // Calculation of the evolution time (relativ)
+
+
+ if( EvolutionTime + fInCycleTime == fCycleTime ) //End of Cycle
+ {
+ fIsAtEndOfCycle = true;
+ fInternalTime += EvolutionTime; // Update Internal Time
+ fInCycleTime += EvolutionTime; // Update InCycleTime
+
+ if(t >= GetCreationTime() + GetLifeTime()) // if the Next Cycle don't 'Exist...
+ fIsShutDown = true;
+
+
+ }
+ else if(EvolutionTime + fInCycleTime < fCycleTime ) // During Cycle
+ {
+
+ fInternalTime += EvolutionTime; // Update Internal Time
+ fInCycleTime += EvolutionTime; // Update InCycleTime
+
+ fInsideIV = fEvolutionDB.GetIsotopicVectorAt( (cSecond)(fInCycleTime/fEvolutionDB.GetPower()*fPower) ); // update the fuel composition
+ if(t>=GetCreationTime() + GetLifeTime()) fIsShutDown = true;
+ }
+ else
+ {
+ // Evolution goes after the end of cycle.... check it
+ ERROR << " " << (*this).GetName() << endl;
+ ERROR << " Evolution Time is " << t << endl;
+ ERROR << " Evolution is too long! There is a problem in Reactor evolution at " << t/365.25/3600/24 << endl;
+ ERROR << " This is too long of : " << EvolutionTime + fInCycleTime - fCycleTime << endl;
+ ERROR << " I have spend " << fInCycleTime +EvolutionTime << " and should have been " << fCycleTime << endl;
+ exit(1);
+ }
+
+DBGL
+}
+
+//________________________________________________________________________
+void Reactor::Dump()
+{
+DBGL
+
+
+ if(fInternalTime < GetCreationTime()) return;
+ if(fIsShutDown && !fIsStarted) return; // Reactor stopped...
+ if(!fIsAtEndOfCycle && !fIsShutDown) return;
+
+// First trash the irradiated fuel
+ if(fIsAtEndOfCycle && !fIsShutDown )
+ {
+ if(fIsStarted ) // A Cycle has already been done
+ {
+ fOutBackEndFacility->AddIV(fInsideIV);
+ AddCumulativeIVOut(fInsideIV);
+ }
+ else fIsStarted = true; // Just start the first cycle
+
+ }
+ else if (fIsAtEndOfCycle && fIsShutDown ) //shutdown at end of Cycle
+ {
+
+ fOutBackEndFacility->AddIV(fIVOutCycle);
+ AddCumulativeIVOut(fIVOutCycle);
+ fInsideIV.Clear();
+ fInCycleTime = 0;
+ fIsStarted = false; // shut down the Reactor
+ }
+ else if (!fIsAtEndOfCycle && fIsShutDown ) //shutdown during Cycle
+ {
+ fOutBackEndFacility->AddIV(fInsideIV);
+ AddCumulativeIVOut(fInsideIV);
+ fInsideIV.Clear();
+ fInCycleTime = 0;
+ fIsStarted = false; // shut down the Reactor
+ }
+
+
+ DBGL
+
+
+// Get the new Fuel !
+ pair<CLASSFuel, double> NextFuel = fFuelPlan->GetFuelAt(fInternalTime);
+ SetBurnUp((NextFuel).second);
+
+ if( NextFuel.first.GetPhysicsModels() )
+ fFixedFuel = false;
+ else if( NextFuel.first.GetEvolutionData() )
+ fFixedFuel = true;
+ else
+ {
+ ERROR << typeid(NextFuel.first).name() << endl;
+ ERROR << "WRONG Fuel Format Correct it !! " << endl;
+ exit(1);
+ }
+ DBGL
+
+ if(fFixedFuel )
+ {
+ DBGL
+ if(fIsAtEndOfCycle && !fIsShutDown )
+ {
+ SetEvolutionDB( *(NextFuel.first.GetEvolutionData()) );
+ fIsAtEndOfCycle = false;
+
+
+ if(!GetParc()->GetStockManagement() && fIsStorage )
+ {
+ IsotopicVector BuildIVtmp ;
+ IsotopicVector OutIncomePart;
+
+ //Get The Storage Compostion
+ BuildIVtmp.Add(fStorage->GetInsideIV().GetIsotopicQuantity());
+ //Get the rest after IVIn creation
+ BuildIVtmp -= fIVInCycle;
+ //Get the OutIncome part form this rest
+ OutIncomePart.Add(BuildIVtmp.GetIsotopicQuantityNeeded()) ;
+ //Take what you can from Storage...
+ fStorage->TakeFromStock( fIVInCycle - OutIncomePart);
+ //And Get the rest from OutIncome
+ GetParc()->AddOutIncome(OutIncomePart);
+
+ }
+ else GetParc()->AddOutIncome(fIVInCycle);
+
+
+ fInsideIV = fIVBeginCycle;
+ AddCumulativeIVIn(fIVBeginCycle);
+
+ fInCycleTime = 0;
+ }
+ DBGL
+
+ }
+ else
+ {
+ DBGL
+ if(!GetParc()->GetStockManagement())
+ {
+ ERROR << " Can't have unfixedFuel without stock management" << endl;
+ exit(1);
+ }
+
+ if(fIsAtEndOfCycle && !fIsShutDown )
+ {
+ fIsAtEndOfCycle = false;
+
+ SetNewFuel(fFabricationPlant->GetReactorEvolutionDB(GetId()));
+ fFabricationPlant->TakeReactorFuel(GetId());
+
+
+ fInsideIV = fIVBeginCycle;
+ AddCumulativeIVIn(fIVBeginCycle);
+
+ fInCycleTime = 0;
+ }
+
+ DBGL
+
+
+
+ }
+
+DBGL
+}
+
+
+
+
diff --git a/source/trunk/src/Scenario.cxx b/source/trunk/src/Scenario.cxx
new file mode 100755
index 000000000..487c45ec0
--- /dev/null
+++ b/source/trunk/src/Scenario.cxx
@@ -0,0 +1,888 @@
+#include "Scenario.hxx"
+
+#include <ctime>
+#include "time.h"
+#include <cmath>
+#include <iomanip>
+#include <fstream>
+#include <sstream>
+#include <algorithm>
+#include <omp.h>
+#include "stdlib.h"
+
+#include "Storage.hxx"
+#include "Reactor.hxx"
+#include "CLASSBackEnd.hxx"
+#include "Pool.hxx"
+#include "FabricationPlant.hxx"
+#include "SeparationPlant.hxx"
+#include "CLASSLogger.hxx"
+
+
+//________________________________________________________________________
+//
+// CLASS
+//
+//
+//
+//
+//________________________________________________________________________
+
+
+
+float random(float a, float b) //peak random numebr between a and b
+{
+ float range = pow(2., 31);
+ srand(time(NULL)); //initialize the srand
+ return (float)a + (float)(b-a)*rand()/range;
+}
+
+string dtoa(double num)
+{
+ ostringstream os(ostringstream::out);
+ os<<setprecision(3)<<num;
+ return os.str();
+}
+
+//________________________________________________________________________
+Scenario::Scenario():CLASSObject(new CLASSLogger("CLASS_OUTPUT.log"))
+{
+
+ fNewTtree = true;
+ fPrintStep = (cSecond)(3600*24*365.25); // One Step per Year
+ fAbsoluteTime = 0;
+ fStartingTime = fAbsoluteTime;
+
+ fStockManagement = true;
+ fLogTimeStep = false;
+
+ fOutputFileName = "CLASS_Default.root";
+ fOutputTreeName = "Data";
+ fOutFile = 0;
+ fOutT = 0;
+ fParcPower = 0;
+
+ // Warning
+
+ INFO << "!!INFO!! !!!Scenario!!! Parc has been define :" << endl;
+ INFO << "\t Print set at : " << (double)(fPrintStep/3600/24/365.25) << " year" << endl;
+ INFO << "\t StockManagement set at : true" << endl;
+ INFO << "\t OutPut will be in \"" << fOutputFileName << "\" File and \"" << fOutputTreeName << "\" TTree" << endl;
+ INFO << "\t Log will be in " << GetLog()->GetCLASSLoggerName() << endl;
+
+}
+
+
+//________________________________________________________________________
+Scenario::Scenario(cSecond abstime):CLASSObject(new CLASSLogger())
+{
+
+ fNewTtree = true;
+ fPrintStep = (cSecond)(3600*24*365.25); // One Step per Year
+ fAbsoluteTime = abstime;
+ fStartingTime = fAbsoluteTime;
+
+ fStockManagement = true;
+ fLogTimeStep = false;
+
+ fOutputFileName = "CLASS_Default.root";
+ fOutputTreeName = "Data";
+ fOutFile = 0;
+ fOutT = 0;
+ fParcPower = 0;
+
+ // Warning
+
+ INFO << "!!INFO!! !!!Scenario!!! Parc has been define :" << endl;
+ INFO << "\t Print set at : " << (double)(fPrintStep/3600/24/365.25) << " year" << endl;
+ INFO << "\t StockManagement set at : true" << endl;
+ INFO << "\t OutPut will be in \"" << fOutputFileName << "\" File and \"" << fOutputTreeName << "\" TTree" << endl;
+ INFO << "\t Log will be in " << GetLog()->GetCLASSLoggerName() << endl;
+
+}
+
+
+//________________________________________________________________________
+Scenario::Scenario(CLASSLogger* log, cSecond abstime):CLASSObject(log)
+{
+
+
+ fNewTtree = true;
+ fPrintStep = (cSecond)(3600*24*365.25); // One Step per Year
+ fAbsoluteTime = abstime;
+ fStartingTime = fAbsoluteTime;
+
+ fStockManagement = true;
+ fLogTimeStep = false;
+
+ fOutputFileName = "CLASS_Default.root";
+ fOutputTreeName = "Data";
+ fOutFile = 0;
+ fOutT = 0;
+
+ fParcPower = 0;
+
+
+ // Warning
+
+
+ INFO << " Parc has been define :" << endl;
+ INFO << " Print set at : " << (double)(fPrintStep/3600/24/365.25) << " year" << endl;
+ INFO << " StockManagement set at : true" << endl;
+ INFO << " OutPut will be in \"" << fOutputFileName << "\" File and \"" << fOutputTreeName << "\" TTree" << endl;
+ INFO << " Log will be in " << GetLog()->GetCLASSLoggerName() << endl;
+
+
+
+}
+//________________________________________________________________________
+Scenario::Scenario(cSecond abstime, CLASSLogger* log):CLASSObject(log)
+{
+
+ fNewTtree = true;
+ fPrintStep = (cSecond)(3600*24*365.25); // One Step per Year
+ fAbsoluteTime = abstime;
+ fStartingTime = fAbsoluteTime;
+
+ fStockManagement = true;
+ fLogTimeStep = false;
+
+ fOutputFileName = "CLASS_Default.root";
+ fOutputTreeName = "Data";
+ fOutFile = 0;
+ fOutT = 0;
+ fParcPower = 0;
+
+ // Warning
+
+ INFO << "!!INFO!! !!!Scenario!!! Parc has been define :" << endl;
+ INFO << "\t Print set at : " << (double)(fPrintStep/3600/24/365.25) << " year" << endl;
+ INFO << "\t StockManagement set at : true" << endl;
+ INFO << "\t OutPut will be in \"" << fOutputFileName << "\" File and \"" << fOutputTreeName << "\" TTree" << endl;
+ INFO << "\t Log will be in " << GetLog()->GetCLASSLoggerName() << endl;
+
+}
+
+
+//________________________________________________________________________
+Scenario::~Scenario()
+{
+
+#pragma omp single
+ {CloseOutputTree();}
+
+
+}
+
+//________________________________________________________________________
+void Scenario::AddPool(Pool* Pool)
+{
+
+
+ fPool.push_back(Pool);
+ fPool.back()->SetParc(this);
+ fPool.back()->SetDecayDataBank( (*this).GetDecayDataBase() );
+ fPool.back()->SetLog(GetLog());
+ fPool.back()->SetId((int)fPool.size()-1);
+ fPool.back()->SetInternalTime(fAbsoluteTime);
+
+
+ string Pool_name = fPool.back()->GetName();
+ if(Pool_name == "P_Pool.")
+ {
+ Pool_name = "P_Pool";
+ Pool_name += dtoa(fPool.back()->GetId());
+ Pool_name += ".";
+ fPool.back()->SetName(Pool_name.c_str());
+ }
+ else
+ {
+ string name_tmp = Pool_name;
+ Pool_name = "P_";
+ Pool_name += name_tmp;
+ Pool_name += ".";
+ fPool.back()->SetName(Pool_name.c_str());
+ }
+
+ if(!fNewTtree)
+ fOutT->Branch(fPool.back()->GetName(), "Pool", &fPool.back());
+
+
+}
+
+//________________________________________________________________________
+void Scenario::AddReactor(Reactor* reactor)
+{
+
+ fReactor.push_back(reactor);
+ fReactor.back()->SetParc(this);
+ fReactor.back()->SetLog(GetLog());
+ fReactor.back()->SetId((int)fReactor.size()-1);
+ fReactor.back()->SetInternalTime(fAbsoluteTime);
+
+
+ string Reactor_name = fReactor.back()->GetName();
+ if(Reactor_name == "R_Reactor.")
+ {
+ Reactor_name = "R_Reactor";
+ Reactor_name += dtoa(fReactor.back()->GetId());
+ Reactor_name += ".";
+ fReactor.back()->SetName(Reactor_name.c_str());
+ }
+ else
+ {
+ string name_tmp = Reactor_name;
+ Reactor_name = "R_";
+ Reactor_name += name_tmp;
+ Reactor_name += ".";
+ fReactor.back()->SetName(Reactor_name.c_str());
+ }
+
+ if(!fNewTtree)
+ fOutT->Branch(fReactor.back()->GetName(), "Reactor", &fReactor.back());
+
+
+}
+
+//________________________________________________________________________
+void Scenario::AddStorage(Storage* storage)
+{
+
+ fStorage.push_back(storage);
+ fStorage.back()->SetParc(this);
+ fStorage.back()->SetDecayDataBank( (*this).GetDecayDataBase() );
+ fStorage.back()->SetLog(GetLog());
+ fStorage.back()->SetId((int)fStorage.size()-1);
+ fStorage.back()->SetInternalTime(fAbsoluteTime);
+
+ string Storage_name = fStorage.back()->GetName();
+
+ if(Storage_name == "S_Storage.")
+ {
+ Storage_name = "S_Storage";
+ Storage_name += dtoa(fStorage.back()->GetId());
+ Storage_name += ".";
+ fStorage.back()->SetName(Storage_name.c_str());
+ }
+ else
+ {
+ string name_tmp = Storage_name;
+ Storage_name = "S_";
+ Storage_name += name_tmp;
+ Storage_name += ".";
+ fStorage.back()->SetName(Storage_name.c_str());
+ }
+
+ if(!fNewTtree)
+ fOutT->Branch(fStorage.back()->GetName(), "Storage", &fStorage.back());
+
+
+}
+//________________________________________________________________________
+void Scenario::AddFabricationPlant(FabricationPlant* fabricationplant)
+{
+
+ fFabricationPlant.push_back(fabricationplant);
+ fFabricationPlant.back()->SetParc(this);
+ fFabricationPlant.back()->SetDecayDataBank( (*this).GetDecayDataBase() );
+ fFabricationPlant.back()->SetLog(GetLog());
+ fFabricationPlant.back()->SetId((int)fStorage.size()-1);
+ fFabricationPlant.back()->SetInternalTime(fAbsoluteTime);
+
+
+ string FP_name = fFabricationPlant.back()->GetName();
+ if(FP_name == "F_FabricationPlant.")
+ {
+ FP_name = "F_FabricationPlant";
+ FP_name += dtoa(fFabricationPlant.back()->GetId());
+ FP_name += ".";
+ fFabricationPlant.back()->SetName(FP_name.c_str());
+ }
+ else
+ {
+ string name_tmp = FP_name;
+ FP_name = "F_";
+ FP_name += name_tmp;
+ FP_name += ".";
+ fFabricationPlant.back()->SetName(FP_name.c_str());
+ }
+
+ if(!fNewTtree)
+ fOutT->Branch(fFabricationPlant.back()->GetName(), "FabricationPlant", &fFabricationPlant.back());
+}
+//________________________________________________________________________
+void Scenario::AddSeparationPlant(SeparationPlant* SeparationPlant)
+{
+
+
+ fSeparationPlant.push_back(SeparationPlant);
+ fSeparationPlant.back()->SetParc(this);
+ fSeparationPlant.back()->SetDecayDataBank( (*this).GetDecayDataBase() );
+ fSeparationPlant.back()->SetLog(GetLog());
+ fSeparationPlant.back()->SetId((int)fSeparationPlant.size()-1);
+ fSeparationPlant.back()->SetInternalTime(fAbsoluteTime);
+
+
+ string SeparationPlant_name = fSeparationPlant.back()->GetName();
+ if(SeparationPlant_name == "C_SepPlant.")
+ {
+ SeparationPlant_name = "C_SepPlant";
+ SeparationPlant_name += dtoa(fSeparationPlant.back()->GetId());
+ SeparationPlant_name += ".";
+ fSeparationPlant.back()->SetName(SeparationPlant_name.c_str());
+ }
+ else
+ {
+ string name_tmp = SeparationPlant_name;
+ SeparationPlant_name = "C_";
+ SeparationPlant_name += name_tmp;
+ SeparationPlant_name += ".";
+ fSeparationPlant.back()->SetName(SeparationPlant_name.c_str());
+ }
+
+ if(!fNewTtree)
+ fOutT->Branch(fSeparationPlant.back()->GetName(), "SeparationPlant", &fSeparationPlant.back());
+
+
+}
+
+
+//________________________________________________________________________
+
+
+//________________________________________________________________________
+void Scenario::BuildTimeVector(cSecond t)
+{
+ DBGL
+ fTimeStep.clear();
+ fTimeStep.insert( pair<cSecond ,int>(t,1) );
+ //********* Printing Step *********//
+ {
+ DBGL
+ cSecond step = fStartingTime;
+
+ if(step >= fAbsoluteTime )
+ fTimeStep.insert( pair<cSecond ,int>(step,1) );
+ step += fPrintStep;
+ cSecond timescale = 1;
+
+ do
+ {
+
+ if(step >= fAbsoluteTime )
+ fTimeStep.insert( pair<cSecond ,int>(step,1) );
+
+ if(fLogTimeStep)
+ {
+ timescale *= 10;
+ step = fPrintStep*timescale;
+ }
+ else
+ step += fPrintStep;
+
+ }
+ while( step < t );
+ DBGL
+ }
+
+
+ for(int i = 0; i < (int)fReactor.size();i++)
+ {
+ DBGL
+ cSecond R_StartingTime = fReactor[i]->GetCreationTime();
+ cSecond R_ShutDownTime = fReactor[i]->GetCreationTime() + fReactor[i]->GetLifeTime();
+
+ double R_Power = fReactor[i]->GetPower();
+ double R_HMMass = fReactor[i]->GetHeavyMetalMass();
+ pair<CLASSFuel, double> R_Fuel = fReactor[i]->GetFuelPlan()->GetFuelAt(R_StartingTime);
+
+ double R_BU = R_Fuel.second;
+ cSecond R_CycleTime = (cSecond) (R_BU / R_Power * R_HMMass * 1e9 *3600*24);
+ if(R_CycleTime == 0)
+ {
+ ERROR << " Be carefull a reactor cycletime is set to 0 second....\"\n" << endl;
+ exit(1);
+ }
+
+ int R_FacilityType = fReactor[i]->GetFacilityType();
+
+
+ cSecond F_CycleTime = 0;
+
+
+ cSecond step = R_StartingTime;
+
+ map< cSecond, int > R_BackEndTimePath = fReactor[i]->GetOutBackEndFacility()->GetTheBackEndTimePath();
+
+ if( R_Fuel.first.GetPhysicsModels() )
+ F_CycleTime = fReactor[i]->GetFabricationPlant()->GetCycleTime();
+
+
+ //********* Reactor Evolution Step *********//
+ // ShutDown of a reactor
+
+ // Test if the sutdown of the reactor is after the actual time (AbsolutreTime) and before the end of the evolution (t)
+ if( R_ShutDownTime < t )
+ {
+ // Shutdown
+ if( R_ShutDownTime > fAbsoluteTime)
+ {
+ pair< map<cSecond, int>::iterator, bool > IResult;
+ IResult = fTimeStep.insert( pair<cSecond ,int>(R_ShutDownTime, 2) );
+ if( !IResult.second )
+ IResult.first->second |= 2;
+ }
+
+ // BackEnd fuel Cycle after reactor Shutdown
+
+ map< cSecond, int >::iterator TV_it; // the time vector iterator
+
+ // Loop on the BackEnd fuel Cycle Time path
+ for(TV_it = R_BackEndTimePath.begin(); TV_it != R_BackEndTimePath.end(); TV_it++)
+ {
+ // Test if each step of the Fuel Cycle BackEnd is after the actual time (AbsolutreTime) and before the end of the evolution (t)
+ if( R_ShutDownTime + (*TV_it).first >= fAbsoluteTime && R_ShutDownTime + (*TV_it).first <= t)
+ {
+ pair< map<cSecond, int>::iterator, bool > IResult;
+ IResult = fTimeStep.insert( pair<cSecond ,int>(R_ShutDownTime + (*TV_it).first, (*TV_it).second) );
+ if( !IResult.second )
+ IResult.first->second |= (*TV_it).second;
+ }
+ }
+
+
+ }
+
+ // Start the reactor and the Fuel Fabrication
+ if(step >= fAbsoluteTime && step <= t && step < R_ShutDownTime)
+ {
+ pair< map<cSecond, int>::iterator, bool > IResult;
+ IResult = fTimeStep.insert( pair<cSecond ,int>(step, R_FacilityType) );
+ if( !IResult.second )
+ IResult.first->second |= R_FacilityType;
+ }
+
+
+ //********* FabricationPlant Evolution Step *********//
+
+
+ if( R_Fuel.first.GetPhysicsModels() )
+ {
+
+ fReactor[i]->GetFabricationPlant()->AddReactor( i, step );
+
+
+ F_CycleTime = fReactor[i]->GetFabricationPlant()->GetCycleTime();
+
+ if(step - F_CycleTime >= fAbsoluteTime && step - F_CycleTime <= t && step < R_ShutDownTime)
+ { // Set End of reactor cycle
+ pair< map<cSecond, int>::iterator, bool > IResult;
+ IResult = fTimeStep.insert( pair<cSecond ,int>(step - F_CycleTime,16) );
+ if( !IResult.second ) IResult.first->second |= 16;
+ }
+ else if( step - F_CycleTime < fStartingTime )
+ {
+ ERROR << " Can't Build Fuel before Scenario's start\"\n" << endl;
+ exit(1);
+ }
+ }
+
+ step += R_CycleTime;
+ //Prepare the first Cycle
+ R_Fuel = fReactor[i]->GetFuelPlan()->GetFuelAt(step);
+
+ R_BU = fReactor[i]->GetFuelPlan()->GetFuelAt(step).second;
+ R_CycleTime = (cSecond) (R_BU / R_Power * R_HMMass * 1e9 *3600*24);
+
+ if(R_CycleTime == 0)
+ {
+ ERROR << " Be carefull a reactor cycletime is set to 0 second....\"\n" << endl;
+ exit(1);
+ }
+
+
+ while(step <= t && step <= R_ShutDownTime )
+ {
+ DBGL
+
+ // FabricationPlant Evolution Step
+ if( R_Fuel.first.GetPhysicsModels() )
+ {
+ fReactor[i]->GetFabricationPlant()->AddReactor( i, step );
+
+ F_CycleTime = fReactor[i]->GetFabricationPlant()->GetCycleTime();
+
+ if(step - F_CycleTime >= fAbsoluteTime && step - F_CycleTime <= t && step < R_ShutDownTime)
+ { // Set End of reactor cycle
+ pair< map<cSecond, int>::iterator, bool > IResult;
+ IResult = fTimeStep.insert( pair<cSecond ,int>(step - F_CycleTime,16) );
+ if( !IResult.second ) IResult.first->second |= 16;
+ }
+ }
+
+
+ if(step >= fAbsoluteTime && step <= t && step < R_ShutDownTime)
+ { // Set End of reactor cycle
+ pair< map<cSecond, int>::iterator, bool > IResult = fTimeStep.insert( pair<cSecond ,int>(step,4) );
+ if( !IResult.second ) IResult.first->second |= 4;
+ }
+
+ // End/Start Of Reactor Cycle Step //
+
+
+ map< cSecond, int >::iterator TV_it; // the time vector iterator
+ // BackEnd fuel Cycle
+ // Loop on the BackEnd fuel Cycle Time path
+ for(TV_it = R_BackEndTimePath.begin(); TV_it != R_BackEndTimePath.end(); TV_it++)
+ {
+
+ if(step + (*TV_it).first >= fAbsoluteTime &&
+ step + (*TV_it).first <= t)
+ { // Test if each step of the Fuel Cycle BackEnd is after the actual time (AbsolutreTime) and before the end of the evolution (t)
+ pair< map<cSecond, int>::iterator, bool > IResult;
+
+ IResult = fTimeStep.insert( pair<cSecond ,int>(step + (*TV_it).first, (*TV_it).second) );
+ if( !IResult.second )
+ IResult.first->second |= (*TV_it).second;
+ }
+ }
+
+
+ step += R_CycleTime;
+
+ // Update to the next fuel
+ R_Fuel = fReactor[i]->GetFuelPlan()->GetFuelAt(step);
+
+ R_BU = fReactor[i]->GetFuelPlan()->GetFuelAt(step).second;
+ R_CycleTime = (cSecond) (R_BU / R_Power * R_HMMass * 1e9 *3600*24);
+ if(R_CycleTime == 0)
+ {
+ ERROR << " Be carefull a reactor cycletime is set to 0 second....\"\n" << endl;
+ exit(1);
+ }
+
+ DBGL
+ }
+
+
+
+ DBGL
+ }
+ //****** Print the Time Index ******//
+ ofstream TimeStepfile("CLASS_TimeStep", ios_base::app); // Open the File
+
+ if(!TimeStepfile)
+ WARNING << " Can't open \" CLASS_TimeStep \"\n" << endl;
+
+ map<cSecond ,int >::iterator it;
+ for( it = fTimeStep.begin(); it != fTimeStep.end(); it++)
+ TimeStepfile << (*it).first << " " << (*it).second << endl;
+
+ DBGL
+}
+
+
+//________________________________________________________________________
+//___________________________ Evolution Method ___________________________
+//________________________________________________________________________
+
+void Scenario::BackEndEvolution()
+{
+ DBGL
+ StorageEvolution();
+ PoolEvolution();
+
+ PoolDump();
+ DBGL
+}
+
+void Scenario::PoolEvolution()
+{
+ DBGL
+#pragma omp parallel for
+ for(int i = 0; i < (int) fPool.size();i++)
+ fPool[i]->Evolution(fAbsoluteTime);
+
+ DBGL
+}
+
+void Scenario::StorageEvolution()
+{
+ DBGL
+#pragma omp parallel for
+ for(int i = 0; i < (int) fStorage.size();i++)
+ fStorage[i]->Evolution(fAbsoluteTime);
+
+ DBGL
+}
+
+void Scenario::FabricationPlantEvolution()
+{
+ DBGL
+ //#pragma omp parallel for
+ for(int i = 0; i < (int) fFabricationPlant.size();i++)
+ fFabricationPlant[i]->Evolution(fAbsoluteTime);
+
+ DBGL
+}
+
+
+void Scenario::PoolDump()
+{
+ DBGL
+ for(int i = 0; i < (int) fPool.size();i++)
+ fPool[i]->Dump();
+ DBGL
+}
+
+//________________________________________________________________________
+void Scenario::ReactorEvolution()
+{
+ DBGL
+ fParcPower = 0;
+#pragma omp parallel for
+ for(int i = 0; i < (int)fReactor.size(); i++)
+ fReactor[i]->Evolution(fAbsoluteTime);
+
+
+ for(int i = 0; i < (int)fReactor.size(); i++)
+ fReactor[i]->Dump();
+
+ DBGL
+}
+
+//________________________________________________________________________
+void Scenario::Evolution(cSecond t)
+{
+ DBGL
+
+ BuildTimeVector(t);
+
+ if(fNewTtree)
+ {
+ OpenOutputTree();
+ OutAttach();
+ UpdateParc();
+ fOutT->Fill();
+ }
+
+ map<cSecond ,int >::iterator it;
+ for(it = fTimeStep.begin(); it != fTimeStep.end(); it++)
+ {
+
+ fWaste = fDecayDataBase->GetDecay(fWaste, (*it).first - fAbsoluteTime);
+ fAbsoluteTime = (*it).first;
+
+ if( (*it).second & 1 || (*it).second & 2 || (*it).second & 4 || (*it).second & 8 || (*it).second & 16 )
+ BackEndEvolution();
+
+ if( (*it).second & 1 || (*it).second & 2 || (*it).second & 4 || (*it).second & 16 )
+ FabricationPlantEvolution();
+
+ if( (*it).second & 1 || (*it).second & 2 || (*it).second & 4 )
+ ReactorEvolution();
+
+ if( (*it).second & 1 || it == fTimeStep.begin() )
+ {
+#pragma omp single
+ {
+ UpdateParc();
+ fOutT->Fill();
+ ProgressPrintout( (cSecond)t);
+ }
+ }
+
+ }
+ cout << endl;
+
+ DBGL
+}
+
+void Scenario::ProgressPrintout(cSecond t)
+{
+
+ double Time = (fAbsoluteTime-fStartingTime)/3600/24/365.25 ;
+ double Total = (t-fStartingTime)/3600/24/365.25;
+
+ // Reset the line
+ for(int i = 0; i < 10; i++)
+ cout << " ";
+ cout << flush ;
+
+ cout << "\r[";
+ for(int i = 0; i < (int)(Time/Total*20.0); i++)
+ cout << "|";
+ for(int i = 20; i >= (int)(Time/Total*20.0); i--)
+ cout << "-";
+ cout << "] ";
+
+ cout << " Processed ";
+ if (Time < 10) cout << " ";
+ if (Time < 100) cout << " ";
+ cout << (int)Time << " / " << (int)Total << " Years \r";
+ if( fLog->GetVerboseLVL() < 2) cout << flush;
+ else cout << endl;
+
+
+ INFO << " Proccessed " << (int)Time << " / " << (int)Total << " Years \r" << endl;
+
+}
+
+//________________________________________________________________________
+//______________________________ Out Method ______________________________
+//________________________________________________________________________
+void Scenario::UpdateParc()
+{
+
+ ResetQuantity();
+
+ for(int i =0; i < (int)fFabricationPlant.size(); i++)
+ fFuelFabrication += fFabricationPlant[i]->GetInsideIV();
+
+ for(int i = 0; i < (int) fPool.size();i++)
+ fTotalCooling += fPool[i]->GetInsideIV();
+
+ for(int i = 0; i < (int)fStorage.size(); i++)
+ fTotalStorage += fStorage[i]->GetInsideIV();
+
+ for(int i = 0; i < (int)fReactor.size(); i++)
+ fTotalInReactor += fReactor[i]->GetIVReactor();
+
+ fIVTotal = fWaste + fTotalStorage + fTotalCooling + fFuelFabrication + fTotalInReactor;
+ fIVInCycleTotal = fTotalStorage + fTotalCooling + fFuelFabrication + fTotalInReactor;
+
+
+}
+
+
+
+void Scenario::ResetQuantity()
+{
+
+ fTotalInReactor.Clear();
+ fTotalStorage.Clear();
+ fTotalCooling.Clear();
+ fFuelFabrication.Clear();
+ fIVInCycleTotal.Clear();
+ fIVTotal.Clear();
+
+}
+
+//________________________________________________________________________
+void Scenario::OpenOutputTree()
+{
+
+
+ INFO << "Opening : " << fOutputFileName << " ...\t";
+ fOutFile = new TFile(fOutputFileName.c_str(),"UPDATE");
+
+ if(!fOutFile)
+ {
+ ERROR << "\nCould not open " << fOutputFileName <<endl;
+ exit(-1);
+ }
+ INFO << "\t ...OK!" << endl;
+
+
+ fOutT = new TTree(fOutputTreeName.c_str(), "Data Tree");
+ INFO << "Creating Data Tree ...\t";
+ if(!fOutT)
+ {
+ ERROR << "\nCould not create Data Tree in " << fOutputFileName << endl;
+ exit(-1);
+ }
+ fNewTtree = false;
+ INFO << "\t ...OK!" << endl;
+
+}
+void Scenario::CloseOutputTree()
+{
+
+
+ fOutFile->ls();
+ INFO << "Writing outTree " << fOutputFileName << endl;
+ fOutFile->Write();
+
+ if(fOutFile->IsOpen()) {
+ INFO << "Deleting outTree : " << endl;
+ delete fOutT;
+ INFO << "Closing file : " << fOutputFileName <<endl;
+ fOutFile-> Close();
+ delete fOutFile;
+ } else {
+ ERROR << "File was not opened " << fOutputFileName << endl;
+ exit(-1);
+ }
+}
+//________________________________________________________________________
+void Scenario::OutAttach()
+{
+
+ ResetQuantity();
+ //Branch Absolut Time
+ fOutT->Branch("AbsTime",&fAbsoluteTime,"AbsoluteTime/L");
+ //Branch The Power installed in the Parc
+ fOutT->Branch("ParcPower",&fParcPower,"ParcPower/D");
+
+
+
+ // Branch the Sum IV
+
+
+ fOutT->Branch("STOCK.", "IsotopicVector", &fTotalStorage);
+ fOutT->Branch("FUELFABRICATION.", "IsotopicVector", &fFuelFabrication);
+ fOutT->Branch("COOLING.", "IsotopicVector", &fTotalCooling);
+ fOutT->Branch("REACTOR.", "IsotopicVector", &fTotalInReactor);
+ fOutT->Branch("INCYCLE.", "IsotopicVector", &fIVInCycleTotal);
+ fOutT->Branch("TOTAL.", "IsotopicVector", &fIVTotal);
+
+ fOutT->Branch("OUTINCOME.", "IsotopicVector", &fOutIncome);
+ fOutT->Branch("WASTE.", "IsotopicVector", &fWaste);
+
+ // Branch the separate object
+
+ for(int i = 0; i < (int)fStorage.size(); i++)
+ fOutT->Branch(fStorage[i]->GetName(), "Storage", &fStorage[i]);
+
+ for(int i = 0; i < (int)fPool.size(); i++)
+
+ fOutT->Branch(fPool[i]->GetName(), "Pool", &fPool[i]);
+
+ for(int i = 0; i < (int)fReactor.size(); i++)
+
+ fOutT->Branch(fReactor[i]->GetName(), "Reactor", &fReactor[i]);
+
+ for(int i = 0; i < (int)fFabricationPlant.size(); i++)
+ fOutT->Branch(fFabricationPlant[i]->GetName(), "FabricationPlant", &fFabricationPlant[i]);
+
+
+}
+
+//________________________________________________________________________
+void Scenario::Write()
+{
+
+
+
+
+}
+
+//________________________________________________________________________
+void Scenario::Print()
+{
+
+ for(int i = 0; i < (int) fPool.size();i++)
+ {
+ INFO << "!!!!!!!!!STEP : " << fAbsoluteTime/(int)(3600*24*365.25) << endl;
+ INFO << "Pool : " << endl;
+ INFO << "Cooling ";
+ INFO << fPool[i]->GetIVArray().size()<< endl;
+ }
+
+ for(int i = 0; i < (int)fReactor.size(); i++)
+ {
+ INFO << "Reactor" << endl;
+ fReactor[i]->GetIVReactor().Print();
+ }
+
+}
diff --git a/source/trunk/src/SeparationPlant.cxx b/source/trunk/src/SeparationPlant.cxx
new file mode 100644
index 000000000..8e75cf5f7
--- /dev/null
+++ b/source/trunk/src/SeparationPlant.cxx
@@ -0,0 +1,160 @@
+#include "SeparationPlant.hxx"
+
+#include "IsotopicVector.hxx"
+#include "Storage.hxx"
+#include "Scenario.hxx"
+#include "CLASSLogger.hxx"
+
+#include <sstream>
+#include <string>
+#include <iostream>
+#include <cmath>
+#include <algorithm>
+
+//________________________________________________________________________
+//
+// SeparationPlant
+//
+//
+//
+//
+//________________________________________________________________________
+ClassImp(SeparationPlant)
+
+
+SeparationPlant::SeparationPlant():CLASSBackEnd(-2)
+{
+ fOutBackEndFacility = 0;
+ SetName("C_SeparationPlant.");
+ SetIsStorageType();
+}
+
+//________________________________________________________________________
+SeparationPlant::SeparationPlant(CLASSLogger* log):CLASSBackEnd(log, -2)
+{
+
+
+ fCycleTime = 0;
+ fIsStarted = false;
+ fPutToWaste = true;
+
+ fOutBackEndFacility = 0;
+ SetName("C_SeparationPlant.");
+
+ SetIsStorageType();
+
+ INFO << " A new SeparationPlant has been define :" << endl;
+ INFO << "\t The Separation Time set at\t " << (double)(fCycleTime/3600/24/365.25) << " year" << endl;
+ WARNING << " All Separated Fuel go directly to WASTE after cooling !! " << endl;
+
+
+}
+
+
+//________________________________________________________________________
+SeparationPlant::~SeparationPlant()
+{
+
+
+}
+
+//________________________________________________________________________
+
+//________________________________________________________________________
+void SeparationPlant::SetBackEndDestination(CLASSBackEnd* storagedestination, IsotopicVector isotopicvector, cSecond destinationstartingtime)
+{
+ DBGL
+
+ fDestinationStorageStartingTime.push_back(destinationstartingtime);
+ fDestinationStorage.push_back(storagedestination);
+ fDestinationStorageIV.push_back(isotopicvector);
+
+ if(fDestinationStorage.size()>=2)
+ {
+ //checker que pas 2 fois le même ZAI
+ //Dans VectorIsotopic, BaM a fait des trucs cools
+ //-> IV->GetThisComposition(IV2)
+ }
+ if (fDestinationStorage.size() != fDestinationStorageIV.size())
+ ERROR << " fDestinationStorage.size() != fDestinationStorageIV.size() !! " << endl;
+
+
+/*
+
+
+
+les pertes non gérées -> WASTE par défaut
+
+
+
+
+*/
+ DBGL
+
+}
+
+
+void SeparationPlant::AddIV(IsotopicVector IV)
+{
+DBGL
+ for(int fds=0; fds<(int)fDestinationStorage.size(); fds++)
+ {
+ cSecond CurrentTime = GetParc()->GetAbsoluteTime();
+
+ INFO << "Separation..." <<endl;
+ INFO << "Current Time : " << CurrentTime <<endl;
+ INFO << "IV Separation Time : " << fDestinationStorageStartingTime[fds] <<endl;
+
+ if(CurrentTime >= fDestinationStorageStartingTime[fds])
+ {
+ IsotopicVector IVtmp;
+ IVtmp = IV*fDestinationStorageIV[fds];
+ fDestinationStorage[fds]->AddIV(IVtmp);
+ IV -= IVtmp;
+ }
+
+ }
+
+ GetParc()->AddWaste(IV);
+
+ DBGL
+}
+
+
+
+map<cSecond,int> SeparationPlant::GetTheBackEndTimePath()
+{
+ DBGL
+
+ map<cSecond, int> TheBackEndTimePath;
+ for( int i = 0; i < (int)fDestinationStorage.size(); i++ )
+ {
+ map<cSecond, int> TheBackEndTimePath_tmp = fDestinationStorage[i]->GetTheBackEndTimePath();
+ map<cSecond, int>::iterator it;
+ for (it = TheBackEndTimePath_tmp.begin(); it != TheBackEndTimePath_tmp.end(); it++)
+ {
+ pair< map<cSecond, int>::iterator, bool > IResult;
+ IResult = TheBackEndTimePath.insert( pair<cSecond ,int>((*it).first, (*it).second) );
+ if( !IResult.second )
+ IResult.first->second |= (*it).second;
+
+ }
+ }
+
+ DBGL
+ return TheBackEndTimePath;
+}
+
+
+
+
+
+//________________________________________________________________________
+// Time Action with the reste of the Universe :
+// Out Storage
+// Evolution :
+// Cooling
+//________________________________________________________________________
+
+
+
diff --git a/source/trunk/src/Storage.cxx b/source/trunk/src/Storage.cxx
new file mode 100644
index 000000000..bd720c8cf
--- /dev/null
+++ b/source/trunk/src/Storage.cxx
@@ -0,0 +1,215 @@
+#include "Storage.hxx"
+
+#include "Scenario.hxx"
+#include "CLASSLogger.hxx"
+
+
+#include <sstream>
+#include <string>
+#include <iostream>
+#include <cmath>
+#include <algorithm>
+
+//________________________________________________________________________
+//
+// Storage
+//
+//
+//
+//
+//________________________________________________________________________
+ClassImp(Storage)
+
+
+
+
+Storage::Storage():CLASSBackEnd(-1)
+{
+ DBGL
+ SetIsStorageType();
+ SetName("S_Storage.");
+
+ DBGL
+}
+
+Storage::Storage(CLASSLogger* log):CLASSBackEnd(log, -1)
+{
+ DBGL
+ SetIsStorageType();
+
+ SetName("S_Storage.");
+
+ INFO << " A new Storage has been define." << endl;
+
+}
+//________________________________________________________________________
+Storage::Storage(CLASSLogger* log, DecayDataBank* evolutivedb):CLASSBackEnd(log, -1)
+{
+ DBGL
+ SetIsStorageType();
+
+ SetDecayDataBank(evolutivedb);
+
+ SetName("S_Storage.");
+
+ INFO << " A new Storage has been define." << endl;
+
+
+ DBGL
+}
+
+//________________________________________________________________________
+Storage::~Storage()
+{
+
+
+}
+
+//________________________________________________________________________
+void Storage::AddIV(IsotopicVector isotopicvector)
+{
+
+ AddCumulativeIVIn(isotopicvector);
+
+ if(GetParc())
+ {
+ if(GetParc()->GetStockManagement() )
+ {
+ fIVArray.push_back(isotopicvector);
+ fIVArrayArrivalTime.push_back(fInternalTime);
+ }
+ }
+ else
+ {
+ fIVArray.push_back(isotopicvector);
+ fIVArrayArrivalTime.push_back(fInternalTime);
+ }
+ UpdateInsideIV();
+}
+
+//________________________________________________________________________
+void Storage::TakeFractionFromStock(int IVId,double fraction)
+{
+DBGL
+ if(GetParc())
+ {
+ if(GetParc()->GetStockManagement() )
+ {
+ if(fraction > 1 || fraction < 0)
+ {
+ WARNING << " You try to remove fraction superior than 1 or a negative one..." << endl;
+ }
+ else
+ {
+ AddCumulativeIVOut(fIVArray[IVId]*fraction);
+ fIVArray[IVId] -= fIVArray[IVId] * fraction;
+ }
+
+ }
+ else
+ {
+ ERROR << " TakeFractionFromStock can't be DEFINE without REAL stock management" << endl;
+ exit(1);
+
+ }
+ }
+ else
+ {
+ if(fraction > 1 || fraction < 0)
+ {
+ WARNING << " You try to remove fraction superior than 1 or a negative one..." << endl;
+ }
+ else
+ {
+ AddCumulativeIVOut(fIVArray[IVId]*fraction);
+ fIVArray[IVId] -= fIVArray[IVId] * fraction;
+ }
+
+ }
+ UpdateInsideIV();
+ DBGL
+}
+
+void Storage::TakeFromStock(IsotopicVector isotopicvector)
+{
+
+ if(GetParc())
+ {
+ if(!GetParc()->GetStockManagement())
+ {
+
+ AddCumulativeIVOut(isotopicvector);
+ fInsideIV -= isotopicvector;
+ }
+ else
+ {
+ ERROR << " TakeFromStock can't be DEFINE WITH REAL stock management" << endl;
+ exit(1);
+ }
+ }
+ else
+ {
+ AddCumulativeIVOut(isotopicvector);
+ fInsideIV -= isotopicvector;
+
+ }
+
+}
+
+//________________________________________________________________________
+void Storage::StorageEvolution(cSecond t)
+{
+DBGL
+
+ if(t == fInternalTime && t !=0 ) return;
+
+ RemoveEmptyStocks();
+
+ cSecond EvolutionTime = t - fInternalTime;
+
+#pragma omp parallel for
+ for (int i=0; i <(int) fIVArray.size() ; i++)
+ {
+ fIVArray[i] = GetDecay(fIVArray[i] , EvolutionTime);
+ }
+DBGL
+}
+
+//________________________________________________________________________
+void Storage::Evolution(cSecond t)
+{
+
+ // Check if the Storage has been created ...
+ if(t == fInternalTime && t!=0) return;
+ // Make the evolution for the Storage ...
+
+ StorageEvolution(t);
+ // Update Inside IV;
+ UpdateInsideIV();
+
+ // ... And Finaly update the AbsoluteInternalTime
+ fInternalTime = t;
+
+}
+
+void Storage::Write(string filename, cSecond date)
+{
+
+ for(int i=0;i < (int)fIVArray.size(); i++)
+ {
+
+ fIVArray[i].Write(filename, date);
+ }
+
+}
+//________________________________________________________________________
+void Storage::RemoveEmptyStocks()
+{
+ for(int i = (int)fIVArray.size()-1 ; i >=0; i--) //Removing empty Stock
+ if(fIVArray[i].GetSumOfAll() == 0)
+ {
+ fIVArray.erase(fIVArray.begin()+i);
+ fIVArrayArrivalTime.erase(fIVArrayArrivalTime.begin()+i);
+
+ }
+}
diff --git a/source/trunk/src/XSModel.cxx b/source/trunk/src/XSModel.cxx
new file mode 100644
index 000000000..8f6aa2ee5
--- /dev/null
+++ b/source/trunk/src/XSModel.cxx
@@ -0,0 +1,61 @@
+//
+// XSModel.cxx
+// CLASSSource
+//
+// Created by BaM on 04/05/2014.
+// Copyright (c) 2014 BLG. All rights reserved.
+//
+
+#include "XSModel.hxx"
+
+using namespace std;
+
+XSModel::XSModel():CLASSObject()
+{
+
+}
+
+
+XSModel::XSModel(CLASSLogger* log):CLASSObject(log)
+{
+
+}
+
+bool XSModel::isIVInDomain(IsotopicVector IV)
+{DBGL
+ bool IsInDomain=true;
+
+ if(fFreshFuelDomain.empty())
+ {
+ WARNING << "Fresh Fuel variation domain is not set" << endl;
+ WARNING << "CLASS has no clue if the computed evolution for this fresh fuel is correct" << endl;
+ WARNING << "Proceed finger crossed !!" << endl;
+ return true;
+ }
+
+ else
+ {
+ IsotopicVector IVNorm = IV /IV.GetSumOfAll();
+ for (map< ZAI,pair<double,double> >::iterator Domain_it=fFreshFuelDomain.begin(); Domain_it!=fFreshFuelDomain.end(); Domain_it++)
+ {
+ double ThatZAIProp = IVNorm.GetIsotopicQuantity()[Domain_it->first] ;
+ double ThatZAIMin = Domain_it->second.first;
+ double ThatZAIMax = Domain_it->second.second;
+ if( (ThatZAIProp > ThatZAIMax) || (ThatZAIProp < ThatZAIMin) )
+ {
+ IsInDomain = false;
+
+ WARNING<<"Fresh fuel out of model range"<<endl;
+ WARNING<<"\t AT LEAST this ZAI is accused to be outrange :"<<endl;
+ WARNING<<"\t\t"<<Domain_it->first.Z()<<" "<<Domain_it->first.A()<<" "<<Domain_it->first.I()<<endl;
+ WARNING<<"\t\t min="<<ThatZAIMin<<" value="<<ThatZAIProp<<" max="<<ThatZAIMax<<endl;
+ WARNING<<"\t IV accused :"<<endl<<endl;
+ WARNING<<IVNorm.sPrint()<<endl;
+ break;
+ }
+ }
+ }
+DBGL
+return IsInDomain;
+
+}
\ No newline at end of file
diff --git a/source/trunk/src/ZAI.cxx b/source/trunk/src/ZAI.cxx
new file mode 100755
index 000000000..a75992d85
--- /dev/null
+++ b/source/trunk/src/ZAI.cxx
@@ -0,0 +1,62 @@
+#include "ZAI.hxx"
+#include "stdlib.h"
+
+
+//const string DEFAULTDATABASE = "DecayBase.dat";
+//________________________________________________________________________
+//
+// ZAI
+//
+//
+//
+//
+//________________________________________________________________________
+//____________________________InClass Operator____________________________
+//________________________________________________________________________
+ClassImp(ZAI)
+
+ZAI ZAI::operator=(ZAI IVa)
+{
+ fZ = IVa.Z();
+ fA = IVa.A();
+ fI = IVa.I();
+ return *this;
+}
+
+
+
+ZAI::ZAI()
+{
+
+ fZ=0;
+ fA=0;
+ fI=0;
+
+}
+
+
+//________________________________________________________________________
+ZAI::ZAI(int Z, int A, int I)
+{
+
+ if( Z > A )
+ {
+ cout << "!!!ERROR!!! " << "[" << __FILE__ << ":" << __FUNCTION__ << "]" << endl;
+ cout << "!!!ERROR!!! Z:" << Z << " is higher than A: " << A << endl;
+ cout << "!!!ERROR!!! CLASS did not manage yet anti-mater!!! Update comming soon !!!" << endl;
+ exit(1);
+ }
+
+ fZ=Z;
+ fA=A;
+ fI=I;
+
+}
+
+
+ZAI::~ZAI()
+{
+
+}
+
+
diff --git a/source/trunk/src/ZAIMass.cxx b/source/trunk/src/ZAIMass.cxx
new file mode 100644
index 000000000..ef2e009e7
--- /dev/null
+++ b/source/trunk/src/ZAIMass.cxx
@@ -0,0 +1,82 @@
+#include "ZAIMass.hxx"
+
+#include "CLASSConstante.hxx"
+
+#include "stdlib.h"
+#include <fstream>
+#include <string>
+
+#include "IsotopicVector.hxx"
+ //________________________________________________________________________
+ //
+ // ZAI
+ //
+ //
+ //
+ //
+ //________________________________________________________________________
+ //____________________________InClass Operator____________________________
+ //________________________________________________________________________
+
+
+ZAIMass::ZAIMass()
+{
+ string CLASSPATH = getenv("CLASS_PATH");
+ string MassDataFile = CLASSPATH + "/data/Mass.dat";
+
+ ifstream infile(MassDataFile.c_str());
+
+ if(!infile.good())
+ {
+ cout << " ZAIMass Error.\n can't find/open file " << MassDataFile << endl;
+ exit(1);
+ }
+
+ int Z,A;
+ string Name;
+ double MassUnity,MassDec,error;
+ while (infile>>Z>>A>>Name>>MassUnity>>MassDec>>error)
+ {
+ double Masse = MassUnity + MassDec * 1e-6;
+ fZAIMass.insert( pair< ZAI,double >( ZAI(Z,A,0), Masse ) );
+ }
+
+ infile.close();
+
+}
+
+
+ZAIMass::~ZAIMass()
+{
+ fZAIMass.clear();
+}
+
+
+double ZAIMass::GetMass(ZAI zai ) const
+{
+ map<ZAI,double> ZAIMasscpy = fZAIMass ;
+ map<ZAI,double>::iterator MassIT = ZAIMasscpy.find( ZAI(zai.Z(), zai.A(), 0) );
+
+ if(MassIT==fZAIMass.end())
+ return zai.A();
+ else
+ return MassIT->second;
+
+}
+
+
+double ZAIMass::GetMass(const IsotopicVector IV) const
+{
+
+ double TotalMass = 0;
+
+ map<ZAI ,double >::iterator it;
+ map<ZAI ,double > isotopicquantity = IV.GetIsotopicQuantity();
+ for( it = isotopicquantity.begin(); it != isotopicquantity.end(); it++)
+ {
+ TotalMass += (*it).second/AVOGADRO * GetMass( (*it).first ) ;
+ }
+
+ return TotalMass*1e-6;
+
+}
\ No newline at end of file
--
GitLab