From e7d09fd88466b1bf9e678d0a1b0b76b731edc69d Mon Sep 17 00:00:00 2001 From: Nicolas Thiolliere <nicolas.thiolliere@subatech.in2p3.fr> Date: Tue, 28 Mar 2017 15:18:02 +0200 Subject: [PATCH] First try for a generic ROOT2ROOT... --- Utils/ROOT2ROOT/GENERIC/CLASS_ROOT2ROOT.cxx | 521 +++++++++++++ Utils/ROOT2ROOT/GENERIC/Makefile | 14 + .../GENERIC/SSENAR/IsotopicVector.cxx | 691 ++++++++++++++++++ .../GENERIC/SSENAR/IsotopicVector.hxx | 202 +++++ Utils/ROOT2ROOT/GENERIC/SSENAR/Linkdef.hxx | 21 + Utils/ROOT2ROOT/GENERIC/SSENAR/Makefile | 15 + Utils/ROOT2ROOT/GENERIC/SSENAR/README.rst | 14 + Utils/ROOT2ROOT/GENERIC/SSENAR/Scenar_t.cxx | 509 +++++++++++++ Utils/ROOT2ROOT/GENERIC/SSENAR/Scenar_t.hxx | 127 ++++ Utils/ROOT2ROOT/GENERIC/SSENAR/ZAI.cxx | 62 ++ Utils/ROOT2ROOT/GENERIC/SSENAR/ZAI.hxx | 96 +++ 11 files changed, 2272 insertions(+) create mode 100644 Utils/ROOT2ROOT/GENERIC/CLASS_ROOT2ROOT.cxx create mode 100644 Utils/ROOT2ROOT/GENERIC/Makefile create mode 100644 Utils/ROOT2ROOT/GENERIC/SSENAR/IsotopicVector.cxx create mode 100644 Utils/ROOT2ROOT/GENERIC/SSENAR/IsotopicVector.hxx create mode 100644 Utils/ROOT2ROOT/GENERIC/SSENAR/Linkdef.hxx create mode 100644 Utils/ROOT2ROOT/GENERIC/SSENAR/Makefile create mode 100644 Utils/ROOT2ROOT/GENERIC/SSENAR/README.rst create mode 100644 Utils/ROOT2ROOT/GENERIC/SSENAR/Scenar_t.cxx create mode 100644 Utils/ROOT2ROOT/GENERIC/SSENAR/Scenar_t.hxx create mode 100644 Utils/ROOT2ROOT/GENERIC/SSENAR/ZAI.cxx create mode 100644 Utils/ROOT2ROOT/GENERIC/SSENAR/ZAI.hxx diff --git a/Utils/ROOT2ROOT/GENERIC/CLASS_ROOT2ROOT.cxx b/Utils/ROOT2ROOT/GENERIC/CLASS_ROOT2ROOT.cxx new file mode 100644 index 000000000..8ca793b0c --- /dev/null +++ b/Utils/ROOT2ROOT/GENERIC/CLASS_ROOT2ROOT.cxx @@ -0,0 +1,521 @@ +/* + + Simple code used to convert multi root output file in one root file + +Needs : + +Select some element of interest given a Z list + +If you are interested in a specific observable you will need to implement the method in Scenar_t.hxx +There is no need to recompile the BIG root file + +Authors: + +BaL +Nico. T. +ZaK + +*/ + +#include "CLASSHeaders.hxx" +#include <iostream> +#include <sstream> +#include <iomanip> +#include <math.h> +#include <vector> +#include <cstdio> +#include <memory> +#include <stdexcept> +#include <string> +#include "TROOT.h" +// +#include "SSENAR/Scenar_t.hxx" + +using namespace std; + +// Get the output of a linux command... +string exec(string s_cmd) +{ + char *cmd = new char[s_cmd.length()+1]; strcpy(cmd,s_cmd.c_str()); + char buffer[128]; + string result = ""; + shared_ptr<FILE> pipe(popen(cmd, "r"), pclose); + if (!pipe) throw std::runtime_error("popen() failed!"); + while (!feof(pipe.get())) { + if (fgets(buffer, 128, pipe.get()) != NULL) + result += buffer; + } + return result; +} + +// Show a progress bar ... + template <class T> +void exec_progress(T i, T Ni) +{ + cout<<"\r =====> "<< (T)((double)i/Ni*100. +1) <<" % "<<flush; +} + +// +// Retrieve quantities of interest +// + +// // Obtain inventories from facilities or from a isotopic vector pointer +// double GetInvOfTime (CLASSFacility *object, IsotopicVector isovec) +// { +// return object->GetInsideIV().GetThisComposition(isovec).GetTotalMass(); +// } +// // +// double GetInvOfTime (IsotopicVector *object, IsotopicVector isovec) +// { +// return object->GetThisComposition(isovec).GetTotalMass(); +// } +// // +double GetInvAtBOC (Reactor *object, IsotopicVector isovec) +{ + return double(object->GetIVBeginCycle().GetThisComposition(isovec).GetTotalMass()); +} +// // +// double GetInvAtEOC (Reactor *object, IsotopicVector isovec) +// { +// return object->GetIVOutCycle().GetThisComposition(isovec).GetTotalMass(); +// } +// +// // +// double GetDiffFlux (CLASSFacility *object, IsotopicVector isovec) +// { +// double cinflux, coutflux; +// +// cinflux = object->GetCumulativeIVIn().GetThisComposition(isovec).GetTotalMass(); +// coutflux = object->GetCumulativeIVOut().GetThisComposition(isovec).GetTotalMass(); +// +// return (cinflux - coutflux); +// } +// // +// + +int main(int argc, char** argv) +{ + + //--------------------------------------------------------------------------------------------------------------------------------------------- + //---------------------------------------------------------------INIT--------------------------------------------------------------------- + //--------------------------------------------------------------------------------------------------------------------------------------------- + + if (argc != 2){ + cout << " Usage : ./ex ROOTFILESDIR " << endl; + return EXIT_FAILURE; + } + // + string s_ROOTDIR = argv[1]; + string s_dirtag = s_ROOTDIR + "/OUT_"; + size_t dirtagsize = s_dirtag.length(); + size_t numformat = 9; // nber of bytes used to print the values in the full root files dir names + + //--------------------------------------------------------------------------------------------------------------------------------------------- + //---------------------------------------------------------------VARIABLES--------------------------------------------------------------------- + //--------------------------------------------------------------------------------------------------------------------------------------------- + + // Time and Power will Get dumped from original OUT.root files : respect the TYPE + Long64_t TimeSecond = 0; + double ParcPower = 0.0; + // + + + // Z list of interest : Only these elements will be dumped in the root file + vector<unsigned short> ZLIST;//{92,94,95,96}; + ZLIST.push_back(92); + ZLIST.push_back(94); + ZLIST.push_back(95); + ZLIST.push_back(96); + + // Isotopes of Interest + // + vector<IsotopicVector> ILIST;//{U,Pu,Am,Cm}; + // + IsotopicVector v_U,v_U5, v_U8; + ZAI U8 = ZAI(92,238,0); + ZAI U5 = ZAI(92,235,0); + v_U = 1*U8 + 1*U5; + ILIST.push_back(v_U); + // + v_U5 = 1*U5; + v_U8 = 1*U8; + + IsotopicVector v_Pu, v_Pu_Fis, v_Pu8, v_Pu9, v_Pu0, v_Pu1, v_Pu2; + ZAI Pu8 = ZAI(94,238,0); + ZAI Pu9 = ZAI(94,239,0); + ZAI Pu0 = ZAI(94,240,0); + ZAI Pu1 = ZAI(94,241,0); + ZAI Pu2 = ZAI(94,242,0); + v_Pu = 1*Pu8 + 1*Pu9 + 1*Pu0 + 1*Pu1 + 1*Pu2; + ILIST.push_back(v_Pu); + // + v_Pu_Fis = 1*Pu9 + 1*Pu1; + v_Pu8 = 1*Pu8; + v_Pu9 = 1*Pu9; + v_Pu0 = 1*Pu0; + v_Pu1 = 1*Pu1; + v_Pu2 = 1*Pu2; + + IsotopicVector v_Am, v_Am1, v_Am3; + ZAI Am1 = ZAI(95,241,0); + ZAI Am2x = ZAI(95,242,1); + ZAI Am3 = ZAI(95,243,0); + v_Am = 1*Am1 + 1*Am2x + 1*Am3; + ILIST.push_back(v_Am); + // + v_Am1 = 1*Am1; + v_Am3 = 1*Am3; + + IsotopicVector v_Cm, v_Cm2, v_Cm4; + ZAI Cm2 = ZAI(96,242,0); + ZAI Cm3 = ZAI(96,243,0); + ZAI Cm4 = ZAI(96,244,0); + v_Cm = 1*Cm2 +1*Cm3 + 1*Cm4; + ILIST.push_back(v_Cm); + // + v_Cm2 = 1*Cm2; + v_Cm4 = 1*Cm4; + + //--------------------------------------------------------------------------------------------------------------------------------------------- + //---------------------------------------------------------------OUTPUT ----------------------------------------------------------------------- + //--------------------------------------------------------------------------------------------------------------------------------------------- + + TFile *FileScenario = new TFile("SSenario.root","RECREATE"); + TTree *TreeScenario = new TTree("TreeScenario","TreeScenario"); + + // + + + //--------------------------------------------------------------------------------------------------------------------------------------------- + //---------------------------------------------------------------INPUT ------------------------------------------------------------------------ + //--------------------------------------------------------------------------------------------------------------------------------------------- + + // Names of ROOT Files + string s_Elements = "find " + s_ROOTDIR + " -type f -name \"*.root\"" ; + string s_Candidates = "find " + s_ROOTDIR + " -type f ! -size -20M -name \"*.root\""; + string s_Crashed = "find " + s_ROOTDIR + " -type f -size -20M -name \"*.root\""; + + // Number and Names of ROOT Files + string s_NumberOfElements = exec(s_Elements + " | wc -l"); + string s_NumberOfCandidates = exec(s_Candidates + " | wc -l"); + string s_NumberOfCrashed = exec(s_Crashed + " | wc -l"); + + + // Writing in corresponding files + s_Elements += " > ROOTFileList.txt"; + s_Candidates += " > ROOTFileCandidates.txt"; + s_Crashed += " > ROOTFileCrashed.txt"; + system(s_Elements.c_str() ); + system(s_Candidates.c_str()); + system(s_Crashed.c_str() ); + + unsigned short NumberOfElements = atoi(s_NumberOfElements.c_str()); + unsigned short NumberOfCandidates = atoi(s_NumberOfCandidates.c_str()); + unsigned short NumberOfCrashed = atoi(s_NumberOfCrashed.c_str()); + + // Only candidates root files will be analysed + string s_OneFileForBranches = exec("sed -n 1p ROOTFileCandidates.txt | tr -d '\040\011\012\015'"); + + // Load ROOT file number f and load TTree + TFile *TFileName = new TFile(s_OneFileForBranches.c_str()); + TTree *fData = new TTree(); + fData = (TTree*) gDirectory->Get(TFileName->GetListOfKeys()->At(TFileName->GetNkeys() - 1)->GetName()); + //Time Steps + Long64_t NTime = fData->GetEntries(); + TFileName->Close(); + + // Branching to Tree + //Scenar_t* Scen = new Scenar_t(NTime -1); + Scenar_t* Scen = new Scenar_t(); + TreeScenario->Branch("ssenar","Scenar_t", &Scen); + + //--------------------------------------------------------------------------------------------------------------------------------------------- + //---------------------------------------------------------------LOAD BRANCHES----------------------------------------------------------------- + //--------------------------------------------------------------------------------------------------------------------------------------------- + + cout<< endl << "#########################" << endl; + cout<<"Number Of ROOT Files : " << NumberOfElements <<endl; + cout<<"Number Of ROOT Candidates Files : " << NumberOfCandidates <<endl; + cout<<"Number Of ROOT Crashed Files : " << NumberOfCrashed <<endl; + cout<<"Number Of Time Step / Files : " << NTime <<endl<<endl; + cout<<"Progression : " << endl; + + string s_ROOTFileName; + + // Failed simulations are dealt by a null score in the observables + // the Tree is filled by default value which are set to zero + ifstream f_ROOTFileCrashed("ROOTFileCrashed.txt"); + for (int f=1; f<=NumberOfCrashed; f++) + { + // Get input arguments + getline(f_ROOTFileCrashed, s_ROOTFileName); + + Scen->BU_UOx = atof(s_ROOTFileName.substr(dirtagsize ,numformat).c_str()); + Scen->TC_UOx = atof(s_ROOTFileName.substr(dirtagsize + 1*(numformat + 1),numformat).c_str()); + Scen->BU_MOx = atof(s_ROOTFileName.substr(dirtagsize + 2*(numformat + 1),numformat).c_str()); + Scen->TC_MOx = atof(s_ROOTFileName.substr(dirtagsize + 3*(numformat + 1),numformat).c_str()); + Scen->TF_MOx = atof(s_ROOTFileName.substr(dirtagsize + 4*(numformat + 1),numformat).c_str()); + Scen->IsMOxAm = atoi(s_ROOTFileName.substr(dirtagsize + 5*(numformat + 1),numformat).c_str()); + Scen->Fr_MOx = atof(s_ROOTFileName.substr(dirtagsize + 6*(numformat + 1),numformat).c_str()); + Scen->NB_Fuel = atof(s_ROOTFileName.substr(dirtagsize + 7*(numformat + 1),numformat).c_str()); + Scen->Ks_Fuel = atof(s_ROOTFileName.substr(dirtagsize + 8*(numformat + 1),numformat).c_str()); + Scen->LF_Fuel = atof(s_ROOTFileName.substr(dirtagsize + 9*(numformat + 1),numformat).c_str()); + Scen->StMMOx = atoi(s_ROOTFileName.substr(dirtagsize + 10*(numformat + 1),numformat).c_str()); + + Scen->TimeStep = 1.0/12.0; // each month + + TreeScenario->Fill(); + Scen->Clear(); + + }f_ROOTFileCrashed.close(); + + // + ifstream f_ROOTFileCandidates("ROOTFileCandidates.txt"); + ofstream f_ROOTFileSelected("ROOTFileSelected.txt"); + // + Scen->NTimeStep = (unsigned short) (NTime - 1); // total nber of time step, just a caution: should be set by ctor + Scen->TimeStep = 1.0/12.0; // each month + + // + NumberOfCandidates = 1001; + for (unsigned short f=1; f<=NumberOfCandidates; f++) + { + unsigned short UOx_NLOAD_Theoric = 0, MOx_NLOAD_Theoric = 0; + Scen->UOx_NLOAD = 0; + Scen->MOx_NLOAD = 0; + Scen->MOx_MLOAD = 0; + + // Get input arguments + getline(f_ROOTFileCandidates, s_ROOTFileName); + + Scen->BU_UOx = atof(s_ROOTFileName.substr(dirtagsize ,numformat).c_str()); + Scen->TC_UOx = atof(s_ROOTFileName.substr(dirtagsize + 1*(numformat + 1),numformat).c_str()); + Scen->BU_MOx = atof(s_ROOTFileName.substr(dirtagsize + 2*(numformat + 1),numformat).c_str()); + Scen->TC_MOx = atof(s_ROOTFileName.substr(dirtagsize + 3*(numformat + 1),numformat).c_str()); + Scen->TF_MOx = atof(s_ROOTFileName.substr(dirtagsize + 4*(numformat + 1),numformat).c_str()); + Scen->IsMOxAm = atoi(s_ROOTFileName.substr(dirtagsize + 5*(numformat + 1),numformat).c_str()); + Scen->Fr_MOx = atof(s_ROOTFileName.substr(dirtagsize + 6*(numformat + 1),numformat).c_str()); + Scen->NB_Fuel = atof(s_ROOTFileName.substr(dirtagsize + 7*(numformat + 1),numformat).c_str()); + Scen->Ks_Fuel = atof(s_ROOTFileName.substr(dirtagsize + 8*(numformat + 1),numformat).c_str()); + Scen->LF_Fuel = atof(s_ROOTFileName.substr(dirtagsize + 9*(numformat + 1),numformat).c_str()); + Scen->StMMOx = atoi(s_ROOTFileName.substr(dirtagsize + 10*(numformat + 1),numformat).c_str()); + + // + Scen->S_IsOk = true; + + // Load ROOT file number f and load TTree + cout<<s_ROOTFileName.c_str()<<endl; + TFile *TFileName = new TFile(s_ROOTFileName.c_str()); + + // avoid corrupted files + if (TFileName->IsZombie()) {TFileName->Close(); Scen->S_IsOk = false; NumberOfCrashed++; continue;} + if (!TFileName->IsOpen()) {Scen->S_IsOk = false; NumberOfCrashed++; continue;} + + // Record Selected root files + f_ROOTFileSelected << s_ROOTFileName << endl; + + // + TTree *fData = new TTree(); + fData = (TTree*) gDirectory->Get(TFileName->GetListOfKeys()->At(TFileName->GetNkeys() - 1)->GetName()); + //fData->Print(); + fData->SetBranchStatus("*", 0); // All branches are unbranched + + // CONNECT BRANCHES + string Branchname, ActiveBranchName; + + fData->SetBranchStatus("AbsTime", 1); fData->SetBranchAddress("AbsTime", &TimeSecond); + fData->SetBranchStatus("ParcPower", 1); fData->SetBranchAddress("ParcPower", &ParcPower); + + // TOTAL IV + IsotopicVector *IV_TOTAL=0; Branchname = "TOTAL"; + fData->SetBranchStatus((Branchname+"*").c_str(), 1); + fData->SetBranchAddress((Branchname+".").c_str(), &IV_TOTAL); + + // WASTE IV + IsotopicVector *IV_WASTE=0; Branchname = "WASTE"; + fData->SetBranchStatus((Branchname+"*").c_str(), 1); + fData->SetBranchAddress((Branchname+".").c_str(), &IV_WASTE); + + // INCYCLE IV + IsotopicVector *IV_INCYCLE=0; Branchname = "INCYCLE"; + fData->SetBranchStatus((Branchname+"*").c_str(), 1); + fData->SetBranchAddress((Branchname+".").c_str(), &IV_INCYCLE); + + // FACILITIES + Reactor* PWR_UOx = new Reactor(); + Branchname = "R_PWR_UOx"; + ActiveBranchName = Branchname + "*"; + fData->SetBranchStatus(ActiveBranchName.c_str(), 1); + fData->SetBranchAddress((Branchname+".").c_str(), &PWR_UOx); + + Reactor* PWR_MOx = new Reactor(); + Branchname = "R_PWR_MOx"; + ActiveBranchName = Branchname + "*"; + fData->SetBranchStatus(ActiveBranchName.c_str(), 1); + fData->SetBranchAddress((Branchname+".").c_str(), &PWR_MOx); + + Storage* StockUOx = new Storage(); + Branchname = "S_StockUOx"; + ActiveBranchName = Branchname + "*"; + fData->SetBranchStatus(ActiveBranchName.c_str(), 1); + fData->SetBranchAddress((Branchname+".").c_str(), &StockUOx); + + Storage* StockMOx = new Storage(); + Branchname = "S_StockMOx"; + ActiveBranchName = Branchname + "*"; + fData->SetBranchStatus(ActiveBranchName.c_str(), 1); + fData->SetBranchAddress((Branchname+".").c_str(), &StockMOx); + + + + //--------------------------------------------------------------------------------------------------------------------------------------------- + //---------------------------------------------------------------LOOP ON EVENTS AND FILE WRITING----------------------------------------------- + //--------------------------------------------------------------------------------------------------------------------------------------------- + + exec_progress(f,NumberOfCandidates); + + // PROCEEDING + for (Long64_t t = 1; t < NTime; t++) //loop over scenario time + { + fData->GetEntry(t); //Update all branched object to the new CLASS time step t + + double tTime = double(TimeSecond)/double(cYear); + Scen->Time.push_back(tTime); //The time (in year) at this time step + Scen->Power.push_back((double) ParcPower); + + // + Scen->UOx_CT = double(PWR_UOx->GetCycleTime())/double(cYear); + Scen->MOx_CT = double(PWR_MOx->GetCycleTime())/double(cYear); + + + // + // Cumulated Number of Load + // + + // New Reactor load - UOx + if (tTime >= double(PWR_UOx->GetCreationTime())/double(cYear) + Scen->UOx_CT*(UOx_NLOAD_Theoric)) + { + UOx_NLOAD_Theoric++; + if (GetInvAtBOC(PWR_UOx, v_U5) > 0) + { + Scen->UOx_NLOAD += 1; + } + } + + // New Reactor load - MOx + if (tTime >= double(PWR_MOx->GetCreationTime())/double(cYear) + Scen->MOx_CT*(MOx_NLOAD_Theoric)) + { + MOx_NLOAD_Theoric++; + if (GetInvAtBOC(PWR_MOx, v_Pu) > 0) Scen->MOx_NLOAD += 1; + else Scen->MOx_MLOAD += 1; + } + // + + // + // + IsotopicVector* PWR_UOx_BOC = new IsotopicVector(); + IsotopicVector* PWR_MOx_BOC = new IsotopicVector(); + + IsotopicVector* PWR_UOx_EOC = new IsotopicVector(); + IsotopicVector* PWR_MOx_EOC = new IsotopicVector(); + + IsotopicVector* StockUOx_CIN = new IsotopicVector(); + IsotopicVector* StockUOx_COU = new IsotopicVector(); + IsotopicVector* StockUOx_INV = new IsotopicVector(); + + IsotopicVector* StockMOx_CIN = new IsotopicVector(); + IsotopicVector* StockMOx_COU = new IsotopicVector(); + IsotopicVector* StockMOx_INV = new IsotopicVector(); + + IsotopicVector* TOTAL_INV = new IsotopicVector(); + IsotopicVector* WASTE_INV = new IsotopicVector(); + IsotopicVector* INCYCLE_INV = new IsotopicVector(); + + //// By ZLIST + //for (unsigned short i = 0; i < ZLIST.size(); i++){ + + // PWR_UOx_BOC->Add(PWR_UOx->GetIVBeginCycle().GetSpeciesComposition(ZLIST[i])); + // PWR_MOx_BOC->Add(PWR_MOx->GetIVBeginCycle().GetSpeciesComposition(ZLIST[i])); + + // PWR_UOx_EOC->Add(PWR_UOx->GetIVOutCycle().GetSpeciesComposition(ZLIST[i])); + // PWR_MOx_EOC->Add(PWR_MOx->GetIVOutCycle().GetSpeciesComposition(ZLIST[i])); + + // StockUOx_CIN->Add(StockUOx->GetCumulativeIVIn().GetSpeciesComposition(ZLIST[i])); + // StockUOx_COU->Add(StockUOx->GetCumulativeIVOut().GetSpeciesComposition(ZLIST[i])); + // StockUOx_INV->Add(StockUOx->GetInsideIV().GetSpeciesComposition(ZLIST[i])); + + // StockMOx_CIN->Add(StockMOx->GetCumulativeIVIn().GetSpeciesComposition(ZLIST[i])); + // StockMOx_COU->Add(StockMOx->GetCumulativeIVOut().GetSpeciesComposition(ZLIST[i])); + // StockMOx_INV->Add(StockMOx->GetInsideIV().GetSpeciesComposition(ZLIST[i])); + + // TOTAL_INV->Add(IV_TOTAL->GetSpeciesComposition(ZLIST[i])); + // WASTE_INV->Add(IV_WASTE->GetSpeciesComposition(ZLIST[i])); + // INCYCLE_INV->Add(IV_INCYCLE->GetSpeciesComposition(ZLIST[i])); + + //} + + // By ILIST + for (unsigned short i = 0; i < ILIST.size(); i++){ + + PWR_UOx_BOC->Add(PWR_UOx->GetIVBeginCycle().GetThisComposition(ILIST[i])); + PWR_MOx_BOC->Add(PWR_MOx->GetIVBeginCycle().GetThisComposition(ILIST[i])); + + PWR_UOx_EOC->Add(PWR_UOx->GetIVOutCycle().GetThisComposition(ILIST[i])); + PWR_MOx_EOC->Add(PWR_MOx->GetIVOutCycle().GetThisComposition(ILIST[i])); + + StockUOx_CIN->Add(StockUOx->GetCumulativeIVIn().GetThisComposition(ILIST[i])); + StockUOx_COU->Add(StockUOx->GetCumulativeIVOut().GetThisComposition(ILIST[i])); + StockUOx_INV->Add(StockUOx->GetInsideIV().GetThisComposition(ILIST[i])); + + StockMOx_CIN->Add(StockMOx->GetCumulativeIVIn().GetThisComposition(ILIST[i])); + StockMOx_COU->Add(StockMOx->GetCumulativeIVOut().GetThisComposition(ILIST[i])); + StockMOx_INV->Add(StockMOx->GetInsideIV().GetThisComposition(ILIST[i])); + + TOTAL_INV->Add(IV_TOTAL->GetThisComposition(ILIST[i])); + WASTE_INV->Add(IV_WASTE->GetThisComposition(ILIST[i])); + INCYCLE_INV->Add(IV_INCYCLE->GetThisComposition(ILIST[i])); + + } + + // + Scen->v_PWR_UOx_BOC.push_back(PWR_UOx_BOC); + Scen->v_PWR_MOx_BOC.push_back(PWR_MOx_BOC); + + Scen->v_PWR_UOx_EOC.push_back(PWR_UOx_EOC); + Scen->v_PWR_MOx_EOC.push_back(PWR_MOx_EOC); + + Scen->v_StockUOx_CIN.push_back(StockUOx_CIN); + Scen->v_StockUOx_COU.push_back(StockUOx_COU); + Scen->v_StockUOx.push_back(StockUOx_INV); + + Scen->v_StockMOx_CIN.push_back(StockMOx_CIN); + Scen->v_StockMOx_COU.push_back(StockMOx_COU); + Scen->v_StockMOx.push_back(StockMOx_INV); + + Scen->v_IVTOTAL.push_back(TOTAL_INV); + Scen->v_IVINCYCLE.push_back(WASTE_INV); + Scen->v_IVWASTE.push_back(INCYCLE_INV); + // + + // + } + TreeScenario->Fill(); + Scen->Clear(); + TFileName->Close(); + } + f_ROOTFileCandidates.close(); + f_ROOTFileSelected.close(); + // + cout<<"ROOT Crashed Files Updated : " << NumberOfCrashed <<endl; + cout<<endl<<"END OF ..."<<endl<<"#########################"<<endl; + + // system("rm -f ROOTFileList.txt"); + // system("rm -f ROOTFileListCrashed.txt"); + + FileScenario->Write(); + FileScenario->Close(); +} + +/* + g++ -std=c++11 -o CLASS_ROOT2ROOT.ex CLASS_ROOT2ROOT.cxx -I $CLASS_include -L $CLASS_lib -L SCENAR -lCLASSpkg -lScenar `root-config --cflags` `root-config --libs` -fopenmp -lgomp -Wunused-result + */ diff --git a/Utils/ROOT2ROOT/GENERIC/Makefile b/Utils/ROOT2ROOT/GENERIC/Makefile new file mode 100644 index 000000000..b35d0bc1c --- /dev/null +++ b/Utils/ROOT2ROOT/GENERIC/Makefile @@ -0,0 +1,14 @@ +SRC=$(wildcard *.cxx) +OBJ=$(SRC:.cxx=.ex) +CC=g++ -std=c++11 +LIB=-I ${CLASS_include} -L ${CLASS_lib} -L SSENAR -lCLASSpkg -lSsenar `root-config --cflags` `root-config --libs` -fopenmp -lgomp -Wunused-result -lTMVA + +all: $(OBJ) + +%.ex: %.cxx + ${CC} -o $@ $< ${LIB} + +clean: + rm *.ex + + diff --git a/Utils/ROOT2ROOT/GENERIC/SSENAR/IsotopicVector.cxx b/Utils/ROOT2ROOT/GENERIC/SSENAR/IsotopicVector.cxx new file mode 100644 index 000000000..d3c2d0e4b --- /dev/null +++ b/Utils/ROOT2ROOT/GENERIC/SSENAR/IsotopicVector.cxx @@ -0,0 +1,691 @@ +#include "IsotopicVector.hxx" + +#include "CLASSLogger.hxx" +#include "CLASSConstante.hxx" + + +#include <cmath> +#include <iostream> +#include <fstream> +#include <string> +#include <algorithm> +//________________________________________________________________________ +//________________________________________________________________________ +// +// +// +// IsotopicVector +// +// +//________________________________________________________________________ +//________________________________________________________________________ + +ClassImp(IsotopicVector) + +//________________________________________________________________________ +//________________________Constructor & Destructor________________________ +//________________________________________________________________________ +IsotopicVector::IsotopicVector () +{ ; } + +IsotopicVector::IsotopicVector ( IsotopicVector const& IVa ) : + fIsotopicQuantity(IVa.fIsotopicQuantity) , fIsotopicQuantityNeeded(IVa.fIsotopicQuantityNeeded) +{ ; } + +//_____________________________________________________GetSpeciesComposition___________________ +IsotopicVector::~IsotopicVector () +{ + //fIsotopicQuantity.clear(); + //fIsotopicQuantityNeeded.clear(); +} + + +//____________________________InClass Operator____________________________ + +//________________________________________________________________________ +IsotopicVector& IsotopicVector::operator+= (const IsotopicVector& IVa) +{ + const_iterator end = IVa.end(); + for ( const_iterator it=IVa.begin() ; it!=end ; ++it ) + { Add( *it ); } + + return *this; +} + +//________________________________________________________________________ +IsotopicVector& IsotopicVector::operator-= (const IsotopicVector& IVa) +{ + const_iterator end = IVa.end(); + for ( const_iterator it=IVa.begin() ; it!=end ; ++it ) + { Remove( *it ); } + + return *this; +} +//________________________________________________________________________ +IsotopicVector& IsotopicVector::operator*= (const IsotopicVector& IVa) +{ + iterator f_end = fIsotopicQuantity.end(); + + const_iterator IVa_end = IVa.end(); + const_iterator jt; + + for( iterator it = fIsotopicQuantity.begin() ; it != f_end ; ++it ) + { + jt = IVa.find( it->first ); + + if ( jt != IVa_end ) + { it->second *= jt->second; } + else + { it->second = 0; } + } + + 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; + } +} + +//________________________________________________________________________ +//_____________________________General Method_____________________________ +//________________________________________________________________________ +void IsotopicVector::Clear() +{ + fIsotopicQuantityNeeded.clear(); + fIsotopicQuantity.clear(); +} +//________________________________________________________________________ +void IsotopicVector::ClearNeed() +{ + fIsotopicQuantityNeeded.clear(); +} + +//________________________________________________________________________ +void IsotopicVector::Multiply(double factor) +{ + 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 +{ + return std::accumulate( + fIsotopicQuantity.begin() , fIsotopicQuantity.end() , // input iterator + 0.0 , // first value for the sum + []( const double sum , const std::pair<ZAI,double> & zaiQ ) { return zaiQ.second + sum; } + ); +} + +//________________________________________________________________________ +void IsotopicVector::Add(const ZAI& zai, double q) +{ + if( ceil(q*1e50) - q*1e50 > q*1e50 - floor(q*1e50) ) + { q = floor(q*1e50)*1/1e50; } + else + { q = ceil(q*1e50)*1/1e50; } + + if ( q > 0 ) + { fIsotopicQuantity[zai] += q; } +} + +//________________________________________________________________________ +void IsotopicVector::Add(const IsotopicVector& isotopicvector) +{ + Add( isotopicvector.fIsotopicQuantity ); +} +//________________________________________________________________________ +void IsotopicVector::Add(const map<ZAI ,double>& quantity) +{ + const_iterator end = quantity.end(); + for ( const_iterator it = quantity.begin() ; it!=end ; ++it ) + { + Add( it->first , it->second ); + } +} + + +//________________________________________________________________________ +void IsotopicVector::Remove(const ZAI& zai, double quantity) +{ + iterator 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& IVa) +{ + const_iterator IVa_end = IVa.end(); + + iterator f_end = fIsotopicQuantity.end(); + iterator jt; + for ( const_iterator it=IVa.begin() ; it!=IVa_end ; ++it ) + { + jt = fIsotopicQuantity.find( it->first ); + if ( jt != f_end ) + { + jt->second -= it->second; + if ( jt->second <= 0 ) + { + fIsotopicQuantity.erase(jt); + f_end = fIsotopicQuantity.end(); + } + } + } +} + +//________________________________________________________________________ + +void IsotopicVector::ApplyZAIThreshold(int z) +{ + double quantity = 0.0; + iterator it = fIsotopicQuantity.begin(); + + while ( it != fIsotopicQuantity.end() ) + { + if ( it->first.Z() < z ) + { + quantity += it->second; + it = fIsotopicQuantity.erase(it++); //becarefull of postcrement incrementation + } + else + { + ++it; + } + } + fIsotopicQuantity.insert( std::pair<ZAI,double>( ZAI(-2,-2,-2) , quantity ) ); +} + + + + +//________________________________________________________________________ +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 +{ + const_iterator it = fIsotopicQuantity.find(zai); + if ( it != fIsotopicQuantity.end() ) + { return it->second; } + + return 0; +} + +//________________________________________________________________________ +double IsotopicVector::GetZAIIsotopicQuantity(const int z, const int a, const int i) const + { return GetZAIIsotopicQuantity(ZAI(z,a,i)); } + +//________________________________________________________________________ +IsotopicVector IsotopicVector::GetSpeciesComposition ( const int z ) const +{ + IsotopicVector tmp; + const_iterator end = fIsotopicQuantity.end(); + + for ( const_iterator it=fIsotopicQuantity.begin() ; it!=end ; ++it ) + { + if ( it->first.Z() == z ) + { + tmp.fIsotopicQuantity.insert( *it ); + } + } + + return tmp; +} + +//________________________________________________________________________ +IsotopicVector IsotopicVector::GetThisComposition(IsotopicVector IVa) const +{ + IsotopicVector tmp; + + const_iterator IVa_end = IVa.end(); + const_iterator f_end = fIsotopicQuantity.end(); + + const_iterator jt; + for ( const_iterator it = IVa.begin() ; it != IVa_end ; ++it ) + { + jt = fIsotopicQuantity.find( it->first ); + if ( jt != f_end ) + { tmp.fIsotopicQuantity.insert( *jt ); } + } + + return tmp; +} +//________________________________________________________________________ +double IsotopicVector::GetTotalMass() const +{ + return cZAIMass.GetMass(*this);//in tons +} +//________________________________________________________________________ + +double IsotopicVector::GetMeanMolarMass() const +{ + return GetTotalMass() * 1e6 * AVOGADRO / GetActinidesComposition().GetSumOfAll(); +} + +//________________________________________________________________________ +vector<ZAI> IsotopicVector::GetZAIList() const +{ + std::vector< ZAI > tmp( fIsotopicQuantity.size() ); + + std::transform ( + fIsotopicQuantity.begin() , fIsotopicQuantity.end() , // input iterator + tmp.begin() , // output iterator + []( const std::pair<ZAI,double> & zaiQ ) { return zaiQ.first; } + ); + + return tmp; +} +//________________________________________________________________________ +void IsotopicVector::Initiatlize(double val) +{ + if ( val == 0 ) { fIsotopicQuantity.clear(); } + else + { + iterator end = fIsotopicQuantity.end(); + for ( iterator it=fIsotopicQuantity.begin() ; it!=end ; ++it ) + { it->second = val; } + } +} + +//________________________________________________________________________ +IsotopicVector IsotopicVector::GetActinidesComposition() const +{ + IsotopicVector tmp; + const_iterator end = fIsotopicQuantity.end(); + + for ( const_iterator it=fIsotopicQuantity.begin() ; it!=end ; ++it ) + { + if ( it->first.Z() >= 89 && it->first.Z() <= 103 ) + { + tmp.fIsotopicQuantity.insert( *it ); + } + } + + return tmp; +} + +vector<int> IsotopicVector::GetChemicalSpecies() const +{ + std::vector< int > tmp( fIsotopicQuantity.size() ); + + std::transform ( + fIsotopicQuantity.begin() , fIsotopicQuantity.end() , // input iterator + tmp.begin() , // output iterator + []( const std::pair<ZAI,double> & zaiQ ) { return zaiQ.first.Z(); } + ); + + std::vector<int>::iterator end = std::unique( tmp.begin() , tmp.end() ); + tmp.resize( end - tmp.begin() ); + + return tmp; +} + + +//________________________________________________________________________ +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/cYear << 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; + } +} + +//________________________________________________________________________ +//__________________________Operator Overlaoding__________________________ +//________________________________________________________________________ + + +IsotopicVector operator+(IsotopicVector const& IVa, IsotopicVector const& IVb) +{ + IsotopicVector IVtmp = IVa; + return IVtmp += IVb; +} + +//________________________________________________________________________ +IsotopicVector operator-(IsotopicVector const& IVa, IsotopicVector const& IVb) +{ + + IsotopicVector 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/(IsotopicVector const& IVA, IsotopicVector const& IVB) +{ + IsotopicVector IVAFromB = IVA.GetThisComposition(IVB); + unsigned int IVAFromB_Size = IVAFromB.GetZAIList().size(); + unsigned int IVA_Size = IVA.GetZAIList().size(); + unsigned int IVB_Size = IVB.GetZAIList().size(); + + if( IVAFromB_Size < IVA_Size) + { + cout << "Something try to divide by zero" << endl; + cout << "IVA / IVB All ZAI in IVA have to be present in IVB " << endl; + exit(0); + } + else if ( IVB_Size > IVAFromB_Size ) + cout << " IVA / IVB : Size of IVB is bigger than IVA. Non commun nuclei will not be present in IVresult" << endl; + + IsotopicVector IVresult; + map<ZAI,double>::const_iterator end = IVAFromB.end(); + map<ZAI,double>::const_iterator begin = IVAFromB.begin(); + map<ZAI,double>::const_iterator it; + + for ( it=begin ; it!=end ; ++it ) + IVresult.Add(it->first,it->second/IVB.GetIsotopicQuantity()[it->first]); + + return IVresult; + +} + +//________________________________________________________________________ +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; +} + + +//____________________________General Operator____________________________ +//________________________________________________________________________ +double RelativDistance ( const IsotopicVector & a, const IsotopicVector & b ) +{ + double d2 = 0; + + IsotopicVector tmp = a + b; + + double Z1total = 0; + double Z2total = 0; + for( IsotopicVector::iterator it = tmp.begin(); it != tmp.end(); ++it ) + { + Z1total += a.GetZAIIsotopicQuantity( it->first ); + Z2total += b.GetZAIIsotopicQuantity( it->first ); + } + double Z1, Z2; + for( IsotopicVector::iterator it = tmp.begin(); it != tmp.end(); ++it ) + { + Z1 = a.GetZAIIsotopicQuantity( it->first ); + Z2 = b.GetZAIIsotopicQuantity( it->first ); + d2 += (Z1/Z1total - Z2/Z2total)*(Z1/Z1total - Z2/Z2total); + } + + return std::sqrt(d2); +} +//____________________________________________________________________________ +double Distance ( const IsotopicVector & a , const IsotopicVector & b , int DistanceType , const IsotopicVector & DistanceParameter ) +{ + if(DistanceType == 0) + { + return DistanceStandard(a,b); + } + else if(DistanceType == 1||DistanceType == 2){ + return DistanceAdjusted(a,b,DistanceParameter); + } + else + { + std::cout << " DistanceType defined by the user isn't recognized by the code" << std::endl; + + exit(1); + } +} +//____________________________________________________________________________ +double DistanceStandard( const IsotopicVector & a , const IsotopicVector & b ) +{ + double d2 = 0.0; + IsotopicVector tmp = a + b; + double Z1,Z2; + for( IsotopicVector::iterator it = tmp.begin(); it != tmp.end(); ++it ) + { + Z1 = a.GetZAIIsotopicQuantity( it->first ); + Z2 = b.GetZAIIsotopicQuantity( it->first ); + d2 += (Z1-Z2)*(Z1-Z2); + } + return std::sqrt(d2); +} +//____________________________________________________________________________ +double DistanceAdjusted( const IsotopicVector & a , const IsotopicVector & b , const IsotopicVector & DistanceParameter ) +{ + double d2 = 0; + IsotopicVector tmp = a + b; + double Z1,Z2; + double lambda; + for( IsotopicVector::iterator it = tmp.begin(); it != tmp.end(); ++it ) + { + Z1 = a.GetZAIIsotopicQuantity( it->first ); + Z2 = b.GetZAIIsotopicQuantity( it->first ); + + lambda = DistanceParameter.GetZAIIsotopicQuantity( it->first ); + d2 += lambda*std::abs(Z1-Z2); + } + return d2; +} +//____________________________________________________________________________ +double Norme( const IsotopicVector & a , int DistanceType , const IsotopicVector & DistanceParameter ) +{ + IsotopicVector zero; + + return Distance(a, zero, DistanceType,DistanceParameter); +} + + diff --git a/Utils/ROOT2ROOT/GENERIC/SSENAR/IsotopicVector.hxx b/Utils/ROOT2ROOT/GENERIC/SSENAR/IsotopicVector.hxx new file mode 100644 index 000000000..6b58094f8 --- /dev/null +++ b/Utils/ROOT2ROOT/GENERIC/SSENAR/IsotopicVector.hxx @@ -0,0 +1,202 @@ +#ifndef _ISOTOPICVECTOR_ +#define _ISOTOPICVECTOR_ + + +/*! + \file + \brief Header file for IsotopicVector class. + @version 2.0 + */ +#include "ZAI.hxx" + +#include "TObject.h" +#include <string> +#include <vector> +#include <map> + +using namespace std; +typedef long long int cSecond; + +//-----------------------------------------------------------------------------// +//! Allows to store & operate on radioactive sample + +/*! + Defines an Isotopicvector. + An isotopicVector is a map of ZAI and double (e.g number of atoms). + Its aim is to define a radioactive sample. + + @author BaM + @author BLG + @author Marc + @version 2.0 + */ +//________________________________________________________________________ + + + +class IsotopicVector : public TObject +{ + public : + typedef std::map<ZAI,double>::const_iterator const_iterator; + typedef std::map<ZAI,double>::iterator iterator; + + //********* Constructor/Destructor Method *********// + /*! + \name Constructor/Desctructor + */ + //@{ + + IsotopicVector(); ///< Normal Constructor. + IsotopicVector ( IsotopicVector const& IVa ); ///< Copy Constructor. + + + ~IsotopicVector(); ///< Normal Destructor. + + //@} + + //********* Operator Method *********// + /*! + \name Internal Operator + */ + //@{ + + IsotopicVector& operator += ( IsotopicVector const& IVb ); //!< Operator += definition + IsotopicVector& operator -= ( IsotopicVector const& IVb ); //!< Operator -= definition + IsotopicVector& operator *= ( IsotopicVector const& IVb ); //!< Operator *= definition + IsotopicVector& operator *= ( double const& factor ); //!< Operator *= definition (scalar) + bool operator < (const IsotopicVector& isotopicvector) const; //!< IsotopicVector Comparator + + //@} + + //********* Get Method *********// + /*! + \name Get Method + */ + //@{ + + map<ZAI,double> GetIsotopicQuantity () const + { return fIsotopicQuantity; } //!< Return the IsotopicVector as a map + map<ZAI,double> GetIsotopicQuantityNeeded () const + { return fIsotopicQuantityNeeded; } //!< Return the needed IsotopicVector as a map + IsotopicVector GetSpeciesComposition ( int z ) const; //!< Return the Species composition of the "z" atom + IsotopicVector GetThisComposition ( IsotopicVector IV ) const; //!< Return the composition according the IV list... + vector<ZAI> GetZAIList () const; //!< Return the list of ZAI present in the IV + IsotopicVector GetActinidesComposition() const; //!< Return the Actinides composition (Z >= 89) + + double GetZAIIsotopicQuantity ( const ZAI& zai ) const; ///< Return the quantity of the ZAI + double GetZAIIsotopicQuantity ( const int z, const int a, const int i) const; ///< Return the quantity of the ZAI + double GetQuantity ( const int z, const int a, const int i ) const { return GetZAIIsotopicQuantity(z,a,i); } + double GetQuantity ( const ZAI& zai ) const { return GetZAIIsotopicQuantity(zai); } + + void Initiatlize ( double val ); + + + double GetTotalMass () const; //!< Return the mass (in tons) of the isotopic vector + double GetMeanMolarMass () const; //<! Return the mean molar mass of the isotopic vector + + vector<int> GetChemicalSpecies () const; //!< Return the list of chemichal species contained + int GetZAIQuantity () const + {return fIsotopicQuantity.size(); } //!< Return the number of different ZAI in the IsotopicVector + + double GetSumOfAll () const; //!< Return the Sum of nuclei in the IsotopicVector + + iterator begin () + { return fIsotopicQuantity.begin(); } //!< Return an iterator on the begin of fIsotopicQuantity + const_iterator begin () const + { return fIsotopicQuantity.begin(); } //!< Return a constant iterator on the begin of fIsotopicQuantity + iterator end () + { return fIsotopicQuantity.end(); } //!< Return an iterator on the end of fIsotopicQuantity + const_iterator end () const + { return fIsotopicQuantity.end(); } //!< Return a constant iterator on the end of fIsotopicQuantity + + iterator find ( const ZAI & zai ) { return fIsotopicQuantity.find(zai); } + const_iterator find ( const ZAI & zai ) const { return fIsotopicQuantity.find(zai); } + //@} + + + + +//********* Internal Operation Method *********// + + /*! + \name Internal Operation Method + */ + //@{ + + + void Clear(); //!< Empty all the IV + void ClearNeed(); //!< Empty Need componant of the IV + + void Add(const ZAI& zai, double quantity); //!< Add Quantity gramme of the ZAI Element + void Add(const IsotopicVector& isotopicvector); //!< Add IsotopicVector to the existing IsotopicVector + void Add(const map<ZAI ,double>& quantity); //!< Add IsotopicVector to the existing IsotopicVector + void Add(int Z, int A, int I, double quantity) + { (*this).Add(ZAI(Z,A,I), quantity); } //!< Add Quantity gramme of the ZAI Element + void Add ( const pair<ZAI,double> & zaiQ ) { Add( zaiQ.first , zaiQ.second ); } + + + void Need(const ZAI& zai, double quantity); //!< Fill the fIsotopicQuantityNeeded + void Need(const IsotopicVector& isotopicvector); //!< Fill the fIsotopicQuantityNeeded + void Need(const map<ZAI ,double>& quantityneeded) { fIsotopicQuantityNeeded = quantityneeded; } + //!< Fill the fIsotopicQuantityNeeded + + void Remove(const ZAI& zai, double quantity); //!< Remove Quantity gramme of the ZAI Element + void Remove(const IsotopicVector& isotopicvector); //!< Remove IsotopicVector to the existing IsotopicVector + void Remove ( const pair<ZAI,double> & zaiQ ) { Remove( zaiQ.first , zaiQ.second ); } + + void Multiply(double factor); //!< Multiply the IV by a Factor + + void ApplyZAIThreshold(int z = 90); //!< Put all nuclei below the threshold in -2 -2 -2 ZAI... + //@} + + + +//********* In/Out related Method *********// + + /*! + \name In/Out Method + */ + //@{ + + void Write(string filename, cSecond time = -1 ) const; //!< Write the Content of the IV in the filename file + + void Print(string o = " ") const ; //!< Print the composition of the IV in terminal + string sPrint() const ; //!< Print the composition of the IV in a string + + void PrintList(string o = " ") const ; //!< Print the composition of the IV + + //@} + + +//***************************************************///< + + + protected : + + map<ZAI ,double> fIsotopicQuantity; ///< Isotopic vector composition in atomes Number + map<ZAI ,double> fIsotopicQuantityNeeded; ///< map where negative value are saved + + ClassDef(IsotopicVector,1); +}; + +IsotopicVector operator/(IsotopicVector const& IVA, double F); +IsotopicVector operator/(IsotopicVector const& IVA, IsotopicVector const& IVB); +IsotopicVector operator/(ZAI const& zai, double F); +IsotopicVector operator*(IsotopicVector const& IVA, double F); +IsotopicVector operator*(ZAI const& zai, double F); +IsotopicVector operator*(double F, IsotopicVector const& IVA); +IsotopicVector operator*(double F, ZAI const& zai); +IsotopicVector operator+(IsotopicVector const& IVa, IsotopicVector const& IVb); +IsotopicVector operator-(IsotopicVector const& IVa, IsotopicVector const& IVb); + +IsotopicVector operator*(IsotopicVector const& IVa, IsotopicVector const& IVb); + + +double RelativDistance ( const IsotopicVector & , const IsotopicVector & ); //!< return the euclidean distance between two IV. The two IV are normalize to unity +double Distance ( const IsotopicVector & , const IsotopicVector & , int DistanceType = 0 , const IsotopicVector & DistanceParameter = IsotopicVector() ); //!< return weighted euclidean distance between two IV +double DistanceStandard( const IsotopicVector & , const IsotopicVector & ); //!< return the euclidean distance between two IV +double DistanceAdjusted( const IsotopicVector & , const IsotopicVector & , const IsotopicVector & ); //!< return the weighted euclidean distance between two IV +double Norme( const IsotopicVector & , int DistanceType = 0 , const IsotopicVector & DistanceParameter = IsotopicVector() ); //!< return the norm of an IV + + +#endif diff --git a/Utils/ROOT2ROOT/GENERIC/SSENAR/Linkdef.hxx b/Utils/ROOT2ROOT/GENERIC/SSENAR/Linkdef.hxx new file mode 100644 index 000000000..51a6545d0 --- /dev/null +++ b/Utils/ROOT2ROOT/GENERIC/SSENAR/Linkdef.hxx @@ -0,0 +1,21 @@ +// ===================================================================================== + +#include "ZAI.hxx" +#include "IsotopicVector.hxx" +#include "Scenar_t.hxx" + +#ifdef __CINT__ + +#pragma link off all globals; +#pragma link off all classes; +#pragma link off all functions; +#pragma link C++ nestedclasses; + + +#pragma link C++ class pair<ZAI,double>+; +#pragma link C++ class map<ZAI,double>+; +#pragma link C++ class vector<IsotopicVector*>+; +#pragma link C++ class Scenar_t+; +#pragma link C++ all_functions Scenar_t; + +#endif diff --git a/Utils/ROOT2ROOT/GENERIC/SSENAR/Makefile b/Utils/ROOT2ROOT/GENERIC/SSENAR/Makefile new file mode 100644 index 000000000..cc32d372c --- /dev/null +++ b/Utils/ROOT2ROOT/GENERIC/SSENAR/Makefile @@ -0,0 +1,15 @@ +CC=g++ -std=c++11 -shared +LIBCINT=-I${CLASS_include} +LIB=`root-config --ldflags` -I${ROOTSYS}/include -I${CLASS_include} -L${CLASS_lib} -lCLASSpkg -fPIC -fopenmp -lgomp -Wunused-result + +all: libSsenar.so + +Scenar_t_dict.cxx: ZAI.hxx IsotopicVector.hxx Scenar_t.hxx Linkdef.hxx + rootcint -f $@ -c ${LIBCINT} -p $^ + +libSsenar.so: Scenar_t_dict.cxx ZAI.cxx IsotopicVector.cxx Scenar_t.cxx + ${CC} -o $@ ${LIB} $^ + +clean: + rm libSsenar.so *dict* + diff --git a/Utils/ROOT2ROOT/GENERIC/SSENAR/README.rst b/Utils/ROOT2ROOT/GENERIC/SSENAR/README.rst new file mode 100644 index 000000000..c0a6bf1db --- /dev/null +++ b/Utils/ROOT2ROOT/GENERIC/SSENAR/README.rst @@ -0,0 +1,14 @@ +class Scenar_t +============== + +Les objets de type **Scenar_t** sont utilisés pour décrire un *scénario* CLASS. En pratique les attributs de **Scenar_t** sont les paramètres qui identifient le scénario. Par exemple, dans une étude multiparamétrique ces attributs sont (au minimum) les paramètres avec lesquels on exécutes CLASS_Exec : + +* Burn-Up +* Nbre de rechargements +* Facteur de charge +* Temps de refroidissement +* Temps de fabrication +* Vecteurs Isotopiques (Dans les différentes étapes du cycle combustible) +* ... + +L'utilisation des objets **Scenar_t** implique que l'ensemble des observables d'intérêt, par exemple l'inventaire total des noyaux, ne soient accessibles que par des methodes internes de la class. L'objectif est de pouvoir calculer ce que l'on veut à la demande. diff --git a/Utils/ROOT2ROOT/GENERIC/SSENAR/Scenar_t.cxx b/Utils/ROOT2ROOT/GENERIC/SSENAR/Scenar_t.cxx new file mode 100644 index 000000000..b08404864 --- /dev/null +++ b/Utils/ROOT2ROOT/GENERIC/SSENAR/Scenar_t.cxx @@ -0,0 +1,509 @@ +// ===================================================================================== + +#include "Scenar_t.hxx" + +ClassImp(Scenar_t) // Needed by ROOT + +// +// Const +// ===================================================================================== +Scenar_t::Scenar_t():TObject() +{ + NTimeStep = 0; + S_IsOk = false ; + // + UOx_CT = 0.0; + MOx_CT = 0.0; + // Cumulated Number of Load + UOx_NLOAD = 0; + MOx_NLOAD = 0; + MOx_MLOAD = 0; +} +Scenar_t::Scenar_t(const unsigned short vlength):TObject() +{ + + NTimeStep = vlength; + Time.resize(vlength); // [Year] + Power.resize(vlength); + // + S_IsOk = false ; + // + UOx_CT = 0.0; + MOx_CT = 0.0; + // Cumulated Number of Load + UOx_NLOAD = 0; + MOx_NLOAD = 0; + MOx_MLOAD = 0; + + // Facilities IV of Interest + // + v_PWR_UOx_BOC.resize(vlength); + v_PWR_MOx_BOC.resize(vlength); + + v_PWR_UOx_EOC.resize(vlength); + v_PWR_MOx_EOC.resize(vlength); + + v_StockUOx_CIN.resize(vlength); + v_StockUOx_COU.resize(vlength); + v_StockUOx.resize(vlength); + + v_StockMOx_CIN.resize(vlength); + v_StockMOx_COU.resize(vlength); + v_StockMOx.resize(vlength); + + // Isotopic vectors, CYCLE + // + v_IVTOTAL.resize(vlength); + v_IVINCYCLE.resize(vlength); + v_IVWASTE.resize(vlength); +} +// +// Dest +// ===================================================================================== +Scenar_t::~Scenar_t() +{} + +// +// Methods +// ===================================================================================== + +void Scenar_t::Print() +{ + + cout << " =============================================================== " << endl; + cout << " Scenario Info " << endl; + cout << " BU_UOx : " << this->BU_UOx << endl; + cout << " TC_UOx : " << this->TC_UOx << endl; + cout << " BU_MOx : " << this->BU_MOx << endl; + cout << " TC_MOx : " << this->TC_MOx << endl; + cout << " TF_MOx : " << this->TF_MOx << endl; + cout << " IsMOxAm : " << this->IsMOxAm << endl; + cout << " Fr_MOx : " << this->Fr_MOx << endl; + cout << " NB_Fuel : " << this->NB_Fuel << endl; + cout << " Ks_Fuel : " << this->Ks_Fuel << endl; + cout << " LF_Fuel : " << this->LF_Fuel << endl; + cout << " StMMOx : " << this->StMMOx << endl; + cout << " NTimeStep: " << this->NTimeStep << endl; + cout << " TimeStep : " << this->TimeStep << endl; + cout << " S_IsOk : " << this->S_IsOk << endl; + cout << " =============================================================== " << endl; + +} + + +void Scenar_t::Clear() +{ + + Time.clear(); // [Year] + Power.clear(); + + // Facilities IV of Interest + // + v_PWR_UOx_BOC.clear(); + v_PWR_MOx_BOC.clear(); + + v_PWR_UOx_EOC.clear(); + v_PWR_MOx_EOC.clear(); + + v_StockUOx_CIN.clear(); + v_StockUOx_COU.clear(); + v_StockUOx.clear(); + + v_StockMOx_CIN.clear(); + v_StockMOx_COU.clear(); + v_StockMOx.clear(); + + // Isotopic vectors, CYCLE + // + v_IVTOTAL.clear(); + v_IVINCYCLE.clear(); + v_IVWASTE.clear(); +} + + +// ===================================================================================== + +vector<double> Scenar_t::GetPWR_UOx_BOC(unsigned short Z) +{// + vector<double> v_tBOC; + v_tBOC.reserve(this->NTimeStep); + for (unsigned short i = 0; i < this->NTimeStep; i++){ + v_tBOC.push_back(this->v_PWR_UOx_BOC[i]->GetSpeciesComposition(Z).GetTotalMass()); + } + return v_tBOC; +} + +vector<double> Scenar_t::GetPWR_UOx_BOC(unsigned short Z, unsigned short A, unsigned short I) +{// + IsotopicVector tIV; tIV.Add(Z,A,I,1); + vector<double> v_tBOC; + for (unsigned short i = 0; i < this->NTimeStep; i++){ + v_tBOC.push_back(this->v_PWR_UOx_BOC[i]->GetThisComposition(tIV).GetTotalMass()); + } + return v_tBOC; +} +/// + +vector<double> Scenar_t::GetPWR_MOx_BOC(unsigned short Z) +{// + vector<double> v_tBOC; + v_tBOC.reserve(this->NTimeStep); + for (unsigned short i = 0; i < this->NTimeStep; i++){ + v_tBOC.push_back(this->v_PWR_MOx_BOC[i]->GetSpeciesComposition(Z).GetTotalMass()); + } + return v_tBOC; +} + +vector<double> Scenar_t::GetPWR_MOx_BOC(unsigned short Z, unsigned short A, unsigned short I) +{// + IsotopicVector tIV; tIV.Add(Z,A,I,1); + vector<double> v_tBOC; + v_tBOC.reserve(this->NTimeStep); + for (unsigned short i = 0; i < this->NTimeStep; i++){ + v_tBOC.push_back(this->v_PWR_MOx_BOC[i]->GetThisComposition(tIV).GetTotalMass()); + } + return v_tBOC; +} + +//// + +vector<double> Scenar_t::GetPWR_UOx_EOC(unsigned short Z) +{// + vector<double> v_tEOC; + v_tEOC.reserve(this->NTimeStep); + for (unsigned short i = 0; i < this->NTimeStep; i++){ + v_tEOC.push_back(this->v_PWR_UOx_EOC[i]->GetSpeciesComposition(Z).GetTotalMass()); + } + return v_tEOC; +} + +vector<double> Scenar_t::GetPWR_UOx_EOC(unsigned short Z, unsigned short A, unsigned short I) +{// + IsotopicVector tIV; tIV.Add(Z,A,I,1); + vector<double> v_tEOC; + v_tEOC.reserve(this->NTimeStep); + for (unsigned short i = 0; i < this->NTimeStep; i++){ + v_tEOC.push_back(this->v_PWR_UOx_EOC[i]->GetThisComposition(tIV).GetTotalMass()); + } + return v_tEOC; +} +/// + +vector<double> Scenar_t::GetPWR_MOx_EOC(unsigned short Z) +{// + vector<double> v_tEOC; + v_tEOC.reserve(this->NTimeStep); + for (unsigned short i = 0; i < this->NTimeStep; i++){ + v_tEOC.push_back(this->v_PWR_MOx_EOC[i]->GetSpeciesComposition(Z).GetTotalMass()); + } + return v_tEOC; +} + +vector<double> Scenar_t::GetPWR_MOx_EOC(unsigned short Z, unsigned short A, unsigned short I) +{// + IsotopicVector tIV; tIV.Add(Z,A,I,1); + vector<double> v_tEOC; + v_tEOC.reserve(this->NTimeStep); + for (unsigned short i = 0; i < this->NTimeStep; i++){ + v_tEOC.push_back(this->v_PWR_MOx_EOC[i]->GetThisComposition(tIV).GetTotalMass()); + } + return v_tEOC; +} + + +// ===================================================================================== + +vector<double> Scenar_t::GetTotalInvOfIsotope(unsigned short Z) +{ + vector<double> v_tTotalInv; + v_tTotalInv.reserve(this->NTimeStep); + // + for (unsigned short i =0; i < this->NTimeStep; i++){// + v_tTotalInv.push_back((this->v_IVTOTAL)[i]->GetSpeciesComposition(Z).GetTotalMass()); + } + return v_tTotalInv; +} + +vector<double> Scenar_t::GetTotalInvOfIsotope(unsigned short Z, unsigned short A, unsigned short I) +{ + IsotopicVector tIV; tIV.Add(Z,A,I,1); + vector<double> v_tTotalInv; + v_tTotalInv.reserve(this->NTimeStep); + // + for (unsigned short i =0; i < this->NTimeStep; i++){// + v_tTotalInv.push_back((this->v_IVTOTAL)[i]->GetThisComposition(tIV).GetTotalMass()); + } + return v_tTotalInv; +} + +double Scenar_t::GetTotalInvOfIsotopeAtTime(unsigned short Z, double time) +{ + if (this->v_IVTOTAL.empty()) return 0.0; + else{ + unsigned short i = (unsigned short) (time/(this->TimeStep)); + return (this->v_IVTOTAL)[i]->GetSpeciesComposition(Z).GetTotalMass(); + } +} + +double Scenar_t::GetTotalInvOfIsotopeAtTime(unsigned short Z, unsigned short A, unsigned short I, double time) +{ + if (this->v_IVTOTAL.empty()) return 0.0; + else{ + unsigned short i = (unsigned short) (time/(this->TimeStep)); + IsotopicVector tIV; tIV.Add(Z,A,I,1); + // + return (this->v_IVTOTAL)[i]->GetThisComposition(tIV).GetTotalMass(); + } +} + +// ===================================================================================== + +vector<double> Scenar_t::GetIncycleInvOfIsotope(unsigned short Z) +{ + vector<double> v_tIncycleInv; + v_tIncycleInv.reserve(this->NTimeStep); + // + for (unsigned short i =0; i < this->NTimeStep; i++){// + v_tIncycleInv.push_back((this->v_IVINCYCLE)[i]->GetSpeciesComposition(Z).GetTotalMass()); + } + return v_tIncycleInv; +} + +vector<double> Scenar_t::GetIncycleInvOfIsotope(unsigned short Z, unsigned short A, unsigned short I) +{ + IsotopicVector tIV; tIV.Add(Z,A,I,1); + vector<double> v_tIncycleInv; + v_tIncycleInv.reserve(this->NTimeStep); + // + for (unsigned short i =0; i < this->NTimeStep; i++){// + v_tIncycleInv.push_back((this->v_IVINCYCLE)[i]->GetThisComposition(tIV).GetTotalMass()); + } + return v_tIncycleInv; +} + +double Scenar_t::GetIncycleInvOfIsotopeAtTime(unsigned short Z, double time) +{ + if (this->v_IVINCYCLE.empty()) return 0.0; + else{ + unsigned short i = (unsigned short) (time/(this->TimeStep)); + return (this->v_IVINCYCLE)[i]->GetSpeciesComposition(Z).GetTotalMass(); + } +} + +double Scenar_t::GetIncycleInvOfIsotopeAtTime(unsigned short Z, unsigned short A, unsigned short I, double time) +{ + if (this->v_IVINCYCLE.empty()) return 0.0; + else{ + unsigned short i = (unsigned short) (time/(this->TimeStep)); + IsotopicVector tIV; tIV.Add(Z,A,I,1); + // + return (this->v_IVINCYCLE)[i]->GetThisComposition(tIV).GetTotalMass(); + } +} + + +// ===================================================================================== + +vector<double> Scenar_t::GetWasteInvOfIsotope(unsigned short Z) +{ + vector<double> v_tWasteInv; + v_tWasteInv.reserve(this->NTimeStep); + // + for (unsigned short i =0; i < this->NTimeStep; i++){// + v_tWasteInv.push_back((this->v_IVWASTE)[i]->GetSpeciesComposition(Z).GetTotalMass()); + } + return v_tWasteInv; +} + +vector<double> Scenar_t::GetWasteInvOfIsotope(unsigned short Z, unsigned short A, unsigned short I) +{ + IsotopicVector tIV; tIV.Add(Z,A,I,1); + vector<double> v_tWasteInv; + v_tWasteInv.reserve(this->NTimeStep); + // + for (unsigned short i =0; i < this->NTimeStep; i++){// + v_tWasteInv.push_back((this->v_IVWASTE)[i]->GetThisComposition(tIV).GetTotalMass()); + } + return v_tWasteInv; +} + +double Scenar_t::GetWasteInvOfIsotopeAtTime(unsigned short Z, double time) +{ + if (this->v_IVWASTE.empty()) return 0.0; + else{ + unsigned short i = (unsigned short) (time/(this->TimeStep)); + return (this->v_IVWASTE)[i]->GetSpeciesComposition(Z).GetTotalMass(); + } +} + + +double Scenar_t::GetWasteInvOfIsotopeAtTime(unsigned short Z, unsigned short A, unsigned short I, double time) +{ + if (this->v_IVWASTE.empty()) return 0.0; + else{ + unsigned short i = (unsigned short) (time/(this->TimeStep)); + IsotopicVector tIV; tIV.Add(Z,A,I,1); + // + return (this->v_IVWASTE)[i]->GetThisComposition(tIV).GetTotalMass(); + } +} + + +// ===================================================================================== + +vector<double> Scenar_t::GetStockUOxFlux(unsigned short Z) +{ + vector<double> v_tStockUOxFlux; + v_tStockUOxFlux.reserve(this->NTimeStep); + // + for (unsigned short i =0; i < this->NTimeStep; i++){ + double cumin = (this->v_StockUOx_CIN[i])->GetSpeciesComposition(Z).GetTotalMass(); + double cumout = (this->v_StockUOx_COU[i])->GetSpeciesComposition(Z).GetTotalMass(); + + v_tStockUOxFlux.push_back(cumin - cumout); + } + return v_tStockUOxFlux; +} + +vector<double> Scenar_t::GetStockUOxFlux(unsigned short Z, unsigned short A, unsigned short I) +{ + IsotopicVector tIV; tIV.Add(Z,A,I,1); + vector<double> v_tStockUOxFlux; + v_tStockUOxFlux.reserve(this->NTimeStep); + // + for (unsigned short i =0; i < this->NTimeStep; i++){ + double cumin = (this->v_StockUOx_CIN[i])->GetThisComposition(tIV).GetTotalMass(); + double cumout = (this->v_StockUOx_COU[i])->GetThisComposition(tIV).GetTotalMass(); + v_tStockUOxFlux.push_back(cumin - cumout); + } + return v_tStockUOxFlux; +} + +vector<double> Scenar_t::GetStockUOxInvOfIsotope(unsigned short Z) +{ + vector<double> v_tStockInv; + v_tStockInv.reserve(this->NTimeStep); + // + for (unsigned short i =0; i < this->NTimeStep; i++){// + v_tStockInv.push_back((this->v_StockUOx)[i]->GetSpeciesComposition(Z).GetTotalMass()); + } + return v_tStockInv; +} + + +vector<double> Scenar_t::GetStockUOxInvOfIsotope(unsigned short Z, unsigned short A, unsigned short I) +{ + IsotopicVector tIV; tIV.Add(Z,A,I,1); + vector<double> v_tStockInv; + v_tStockInv.reserve(this->NTimeStep); + // + for (unsigned short i =0; i < this->NTimeStep; i++){// + + v_tStockInv.push_back((this->v_StockUOx)[i]->GetThisComposition(tIV).GetTotalMass()); + + } + return v_tStockInv; +} + + +double Scenar_t::GetStockUOxInvOfIsotopeAtTime(unsigned short Z, double time) +{ + if (this->v_StockUOx.empty()) return 0.0; + else{ + unsigned short i = (unsigned short) (time/(this->TimeStep)); + // + return (this->v_StockUOx)[i]->GetSpeciesComposition(Z).GetTotalMass(); + } +} + + +double Scenar_t::GetStockUOxInvOfIsotopeAtTime(unsigned short Z, unsigned short A, unsigned short I, double time) +{ + if (this->v_StockUOx.empty()) return 0.0; + else{ + unsigned short i = (unsigned short) (time/(this->TimeStep)); + IsotopicVector tIV; tIV.Add(Z,A,I,1); + // + return (this->v_StockUOx)[i]->GetThisComposition(tIV).GetTotalMass(); + } +} + +//// + +vector<double> Scenar_t::GetStockMOxFlux(unsigned short Z) +{ + vector<double> v_tStockMOxFlux; + v_tStockMOxFlux.reserve(this->NTimeStep); + // + for (unsigned short i =0; i < this->NTimeStep; i++){ + double cumin = (this->v_StockMOx_CIN[i])->GetSpeciesComposition(Z).GetTotalMass(); + double cumout = (this->v_StockMOx_COU[i])->GetSpeciesComposition(Z).GetTotalMass(); + + v_tStockMOxFlux.push_back(cumin - cumout); + } + return v_tStockMOxFlux; +} + +vector<double> Scenar_t::GetStockMOxFlux(unsigned short Z, unsigned short A, unsigned short I) +{ + IsotopicVector tIV; tIV.Add(Z,A,I,1); + vector<double> v_tStockMOxFlux; + v_tStockMOxFlux.reserve(this->NTimeStep); + // + for (unsigned short i =0; i < this->NTimeStep; i++){ + double cumin = (this->v_StockMOx_CIN[i])->GetThisComposition(tIV).GetTotalMass(); + double cumout = (this->v_StockMOx_COU[i])->GetThisComposition(tIV).GetTotalMass(); + v_tStockMOxFlux.push_back(cumin - cumout); + } + return v_tStockMOxFlux; +} + +vector<double> Scenar_t::GetStockMOxInvOfIsotope(unsigned short Z) +{ + vector<double> v_tStockInv; + v_tStockInv.reserve(this->NTimeStep); + // + for (unsigned short i =0; i < this->NTimeStep; i++){// + v_tStockInv.push_back((this->v_StockMOx)[i]->GetSpeciesComposition(Z).GetTotalMass()); + } + return v_tStockInv; +} + + +vector<double> Scenar_t::GetStockMOxInvOfIsotope(unsigned short Z, unsigned short A, unsigned short I) +{ + IsotopicVector tIV; tIV.Add(Z,A,I,1); + vector<double> v_tStockInv; + v_tStockInv.reserve(this->NTimeStep); + // + for (unsigned short i =0; i < this->NTimeStep; i++){// + + v_tStockInv.push_back((this->v_StockMOx)[i]->GetThisComposition(tIV).GetTotalMass()); + + } + return v_tStockInv; +} + + +double Scenar_t::GetStockMOxInvOfIsotopeAtTime(unsigned short Z, double time) +{ + if (this->v_StockMOx.empty()) return 0.0; + else{ + unsigned short i = (unsigned short) (time/(this->TimeStep)); + // + return (this->v_StockMOx)[i]->GetSpeciesComposition(Z).GetTotalMass(); + } +} + + +double Scenar_t::GetStockMOxInvOfIsotopeAtTime(unsigned short Z, unsigned short A, unsigned short I, double time) +{ + if (this->v_StockMOx.empty()) return 0.0; + else{ + unsigned short i = (unsigned short) (time/(this->TimeStep)); + IsotopicVector tIV; tIV.Add(Z,A,I,1); + // + return (this->v_StockMOx)[i]->GetThisComposition(tIV).GetTotalMass(); + } +} + diff --git a/Utils/ROOT2ROOT/GENERIC/SSENAR/Scenar_t.hxx b/Utils/ROOT2ROOT/GENERIC/SSENAR/Scenar_t.hxx new file mode 100644 index 000000000..010eca4b6 --- /dev/null +++ b/Utils/ROOT2ROOT/GENERIC/SSENAR/Scenar_t.hxx @@ -0,0 +1,127 @@ +#ifndef SCENAR_T_HXX +#define SCENAR_T_HXX + +#include <vector> + +#include "ZAI.hxx" +#include "IsotopicVector.hxx" +#include "TObject.h" + +// +// data structure +// + +// Length of the vector for the constuctor +// this is actually equivalent to NTimeStep + +class ZAI; +class Isotopicvector; +// +class Scenar_t: public TObject +{ + public: + // + Scenar_t(); + Scenar_t(const unsigned short vlength); + virtual ~Scenar_t(); + // + double BU_UOx; // [GWd/t] + double TC_UOx; + double BU_MOx; + double TC_MOx; + double TF_MOx; + bool IsMOxAm; // + double Fr_MOx; // [%] + double NB_Fuel; + double Ks_Fuel; + double LF_Fuel; + unsigned short StMMOx; + unsigned short NTimeStep; + double TimeStep; + vector<double> Time; // [Year] + vector<double> Power; + // + bool S_IsOk; + // + double UOx_CT; + double MOx_CT; + // Cumulated Number of Load + unsigned short UOx_NLOAD; + unsigned short MOx_NLOAD; + unsigned short MOx_MLOAD; + + // Facilities IV of Interest + // + vector<IsotopicVector*> v_PWR_UOx_BOC; + vector<IsotopicVector*> v_PWR_MOx_BOC; + + vector<IsotopicVector*> v_PWR_UOx_EOC; + vector<IsotopicVector*> v_PWR_MOx_EOC; + + vector<IsotopicVector*> v_StockUOx_CIN; + vector<IsotopicVector*> v_StockUOx_COU; + vector<IsotopicVector*> v_StockUOx; + + vector<IsotopicVector*> v_StockMOx_CIN; + vector<IsotopicVector*> v_StockMOx_COU; + vector<IsotopicVector*> v_StockMOx; + + // Isotopic vectors, CYCLE + // + vector<IsotopicVector*> v_IVTOTAL; + vector<IsotopicVector*> v_IVINCYCLE; + vector<IsotopicVector*> v_IVWASTE; + + // + // Methods + // + + void Print(); + void Clear(); + // + vector<double> GetPWR_UOx_BOC(unsigned short Z); // + vector<double> GetPWR_UOx_BOC(unsigned short Z, unsigned short A, unsigned short I=0); + vector<double> GetPWR_MOx_BOC(unsigned short Z); + vector<double> GetPWR_MOx_BOC(unsigned short Z, unsigned short A, unsigned short I=0); + + vector<double> GetPWR_UOx_EOC(unsigned short Z); // + vector<double> GetPWR_UOx_EOC(unsigned short Z, unsigned short A, unsigned short I=0); // + vector<double> GetPWR_MOx_EOC(unsigned short Z); + vector<double> GetPWR_MOx_EOC(unsigned short Z, unsigned short A, unsigned short I=0); + + // + vector<double> GetTotalInvOfIsotope(unsigned short Z); + vector<double> GetTotalInvOfIsotope(unsigned short Z, unsigned short A, unsigned short I=0); + double GetTotalInvOfIsotopeAtTime(unsigned short Z, double Time); + double GetTotalInvOfIsotopeAtTime(unsigned short Z, unsigned short A, unsigned short I, double time); + + vector<double> GetIncycleInvOfIsotope(unsigned short Z); + vector<double> GetIncycleInvOfIsotope(unsigned short Z, unsigned short A, unsigned short I=0); + double GetIncycleInvOfIsotopeAtTime(unsigned short Z, double Time); + double GetIncycleInvOfIsotopeAtTime(unsigned short Z, unsigned short A, unsigned short I, double time); + + vector<double> GetWasteInvOfIsotope(unsigned short Z); + vector<double> GetWasteInvOfIsotope(unsigned short Z, unsigned short A, unsigned short I=0); + double GetWasteInvOfIsotopeAtTime(unsigned short Z, double Time); + double GetWasteInvOfIsotopeAtTime(unsigned short Z, unsigned short A, unsigned short I, double time); + + vector<double> GetStockUOxFlux(unsigned short Z); + vector<double> GetStockUOxFlux(unsigned short Z, unsigned short A, unsigned short I=0); + vector<double> GetStockUOxInvOfIsotope(unsigned short Z); + vector<double> GetStockUOxInvOfIsotope(unsigned short Z, unsigned short A, unsigned short I=0); + double GetStockUOxInvOfIsotopeAtTime(unsigned short Z, double time); + double GetStockUOxInvOfIsotopeAtTime(unsigned short Z, unsigned short A, unsigned short I, double time); + + vector<double> GetStockMOxFlux(unsigned short Z); + vector<double> GetStockMOxFlux(unsigned short Z, unsigned short A, unsigned short I=0); + vector<double> GetStockMOxInvOfIsotope(unsigned short Z); + vector<double> GetStockMOxInvOfIsotope(unsigned short Z, unsigned short A, unsigned short I=0); + double GetStockMOxInvOfIsotopeAtTime(unsigned short Z, double time); + double GetStockMOxInvOfIsotopeAtTime(unsigned short Z, unsigned short A, unsigned short I, double time); + + // + ClassDef (Scenar_t,1) + +}; + +#endif diff --git a/Utils/ROOT2ROOT/GENERIC/SSENAR/ZAI.cxx b/Utils/ROOT2ROOT/GENERIC/SSENAR/ZAI.cxx new file mode 100644 index 000000000..945d5e355 --- /dev/null +++ b/Utils/ROOT2ROOT/GENERIC/SSENAR/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/Utils/ROOT2ROOT/GENERIC/SSENAR/ZAI.hxx b/Utils/ROOT2ROOT/GENERIC/SSENAR/ZAI.hxx new file mode 100644 index 000000000..1ddb39d43 --- /dev/null +++ b/Utils/ROOT2ROOT/GENERIC/SSENAR/ZAI.hxx @@ -0,0 +1,96 @@ +#ifndef _ZAI_ +#define _ZAI_ + +/*! + \file + \brief Header file for ZAI classes. + */ + +#include <string> +#include "TObject.h" +#include <iostream> + +using namespace std; + +//-----------------------------------------------------------------------------// +//! Defines a nucleus + +/*! + Define a nuclei as 3 integer Z,A,I: + \li its charge number : Z + \li its mass number : A + \li its Isomeric state : I (0 : fundamental, 1 : first isomeric state ,...) + + The aim of this class is to discribe each ZAI. + + @author BaM + @version 2.0 + */ +//________________________________________________________________________ + + + +class ZAI : public TObject +{ +public: + + +//********* Constructor/Destructor Method *********// + + /*! + \name Constructor/Desctructor + */ + //@{ + + ZAI(); ///< Default constructor + + //{ + ///< Normal Constructor. + /*! + Default: No parent + \param Z : number of protons + \param A : number of nucleons (A = 0 means natural isotopes) + */ + ZAI(int Z, int A, int I = 0); + //} + + + ~ZAI(); ///< Normal Destructor. + + +//********* ZAI main attributes Method *********// + + /*! + \name ZAI main attributes + */ + //@{ + int Z() const { return fZ; } //!< returns the number of protons + int A() const { return fA; } //!< returns the number of nucleons + int I() const { return fI; } //!< returns the Isomeric State + int N() const { return fA-fZ; } //!< returns the number of neutrons + + //@} + + + + ZAI operator = (ZAI IVa); //!< + bool operator <(const ZAI& zai) const { return (fZ != zai.Z())? + (fZ < zai.Z()) : ( (fA != zai.A())? + (fA < zai.A()) : (fI < zai.I()) ); }//!< Compartor operator + + bool operator != (const ZAI& zai) const { return ( fZ != zai.Z() ) || ( fA != zai.A() ) || ( fI != zai.I() ); }//!< Compartor operator + bool operator == (const ZAI& zai) const { return ( fZ == zai.Z() && fA == zai.A() && fI == zai.I()); }//!< Compartor operator + void Print() const { cout << fZ << " " << fA << " " << fI << endl;}//!< Print in standard output : Z A I + + +protected : + + short fZ; ///< number of protons + short fA; ///< number of nucleons (A = 0 means natural isotopes) + short fI; ///< Isomeric state + + ClassDef(ZAI,1); +}; + + +#endif -- GitLab