From 126e3f62ea040d63cb9265e42bd4576c7a718578 Mon Sep 17 00:00:00 2001
From: Baptiste LENIAU <baptiste.leniau@subatech.in2p3.fr>
Date: Wed, 16 Sep 2015 17:00:26 +0000
Subject: [PATCH] adding Utils to build .xml files for MLP prediction
git-svn-id: svn+ssh://svn.in2p3.fr/class@766 0e7d625b-0364-4367-a6be-d5be4a48d228
---
.../EQM/FBR_MLP_Keff/GenerateRootFile.cxx | 347 ++++++++++++++++++
Utils/trunk/EQM/FBR_MLP_Keff/Train_MLP.cxx | 177 +++++++++
.../FBR_MLP_Keff/include/GenerateRootFile.hxx | 51 +++
.../EQM/FBR_MLP_Keff/include/StringLine.hxx | 275 ++++++++++++++
Utils/trunk/EQM/FBR_MLP_Keff/include/ZAI.hxx | 42 +++
Utils/trunk/EQM/MLP_Kinf/GenerateRootFile.cxx | 343 +++++++++++++++++
Utils/trunk/EQM/MLP_Kinf/Train_MLP.cxx | 178 +++++++++
.../EQM/MLP_Kinf/include/GenerateRootFile.hxx | 51 +++
.../trunk/EQM/MLP_Kinf/include/StringLine.hxx | 275 ++++++++++++++
Utils/trunk/EQM/MLP_Kinf/include/ZAI.hxx | 42 +++
10 files changed, 1781 insertions(+)
create mode 100755 Utils/trunk/EQM/FBR_MLP_Keff/GenerateRootFile.cxx
create mode 100755 Utils/trunk/EQM/FBR_MLP_Keff/Train_MLP.cxx
create mode 100755 Utils/trunk/EQM/FBR_MLP_Keff/include/GenerateRootFile.hxx
create mode 100755 Utils/trunk/EQM/FBR_MLP_Keff/include/StringLine.hxx
create mode 100755 Utils/trunk/EQM/FBR_MLP_Keff/include/ZAI.hxx
create mode 100755 Utils/trunk/EQM/MLP_Kinf/GenerateRootFile.cxx
create mode 100755 Utils/trunk/EQM/MLP_Kinf/Train_MLP.cxx
create mode 100755 Utils/trunk/EQM/MLP_Kinf/include/GenerateRootFile.hxx
create mode 100755 Utils/trunk/EQM/MLP_Kinf/include/StringLine.hxx
create mode 100755 Utils/trunk/EQM/MLP_Kinf/include/ZAI.hxx
diff --git a/Utils/trunk/EQM/FBR_MLP_Keff/GenerateRootFile.cxx b/Utils/trunk/EQM/FBR_MLP_Keff/GenerateRootFile.cxx
new file mode 100755
index 000000000..373bd7037
--- /dev/null
+++ b/Utils/trunk/EQM/FBR_MLP_Keff/GenerateRootFile.cxx
@@ -0,0 +1,347 @@
+/**********************************************************/
+// Make the input file for the MLPs training
+//
+// This programs reads a set of .dat files which are the
+// results of a depletion calculation (see manual and
+// looks for XS_CLOSEST). From the reading it fills a
+// TTree (Data) and write it in a file named
+// TrainingInput.root .
+// The file TrainingInput.cxx is the list of MLP outputs
+// (cross sections)
+//
+//@author BaM, BaL
+/**********************************************************/
+#include "include/GenerateRootFile.hxx"
+#include <TH1F.h>
+#include <TH2D.h>
+#include <TFile.h>
+#include <TTree.h>
+#include "include/StringLine.hxx"
+#include <TString.h>
+#include <string>
+#include <cmath>
+#include <iostream>
+#include <fstream>
+#include <algorithm>
+#include <map>
+#include <sstream>
+
+using namespace std;
+//--------------------------------------------------------------------------------------------------
+/*************************
+ MAIN
+*************************/
+int main(int argc, char ** argv){
+
+ if(argc!=2)
+ {
+ cout << "Usage : BuildInputTree Path" << endl;
+ cout << " Where Path is the path to the folder containing the Evolution Datas" << endl;
+ cout << " i.e the (.dat) files" << endl;
+ exit(0);
+ }
+
+ fEvolutionDataFolder = argv[1];
+
+ InitMass(); //Load nuclei masses
+ CheckJob(); // looks fot the .dat files in the fEvolutionDataFolder
+
+ FilePath = "DB_TMP/";
+ DataPath = FilePath + "Data/";
+ string Command = "mkdir -p " + DataPath;
+ system(Command.c_str());
+
+ cout << "Reading .dat files ..." << endl;
+ for(int i = 0; i < (int)JobName.size(); i++)
+ {
+ ReadAndFill(JobName[i]);
+ if (i%100 == 0)
+ cout << "\r" << i << " .dat files read" <<flush;
+ }
+ cout << "Filling the TTree ..." << endl;
+ DumpInputNeuron("TrainingInput.root");
+ cout << "Training input generated in file : TrainingInput.root " << endl;
+ cout << "Names of MLP outputs in file : TrainingInput.cxx " << endl;
+ system("rm -r DB_TMP");
+
+}
+//--------------------------------------------------------------------------------------------------
+// Function definitions
+//--------------------------------------------------------------------------------------------------
+//Convert int to string
+string itoa(int num)
+{
+ ostringstream os(ostringstream::out);
+ os << num;
+ return os.str();
+}
+//--------------------------------------------------------------------------------------------------
+void DumpInputNeuron(string filename)
+{
+ TFile* fOutFile = new TFile(filename.c_str(),"RECREATE");
+ TTree* fOutT = new TTree("Data", "Data");
+
+/**********************INITIALISATIONNN********************/
+ //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@//
+ //@@@Change the input value according to your fresh fuel compo
+ //-> here MOX FUEL
+ //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@//
+ double U5 = 0;
+ double U8 = 0;
+ double Pu8 = 0;
+ double Pu9 = 0;
+ double Pu10 = 0;
+ double Pu11 = 0;
+ double Pu12 = 0;
+ double Am1 = 0;
+
+ //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@//
+ double k_eff = 0;
+
+
+/**********************BRANCHING**************************************************/
+/**********************Fresh fuel**************************************************/
+ //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@//
+ //@@@Change the input value according to your fresh fuel compo
+ //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@//
+ fOutT->Branch( "U5" ,&U5 ,"U5/D" );
+ fOutT->Branch( "U8" ,&U8 ,"U8/D" );
+ fOutT->Branch( "Pu8" ,&Pu8 ,"Pu8/D" );
+ fOutT->Branch( "Pu9" ,&Pu9 ,"Pu9/D" );
+ fOutT->Branch( "Pu10" ,&Pu10 ,"Pu10/D" );
+ fOutT->Branch( "Pu11" ,&Pu11 ,"Pu11/D" );
+ fOutT->Branch( "Pu12" ,&Pu12 ,"Pu12/D" );
+ fOutT->Branch( "Am1" ,&Am1 ,"Am1/D" );
+ //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@//
+ fOutT->Branch( "k_eff" ,&k_eff ,"k_eff/D" );
+
+/**********************FILLING THE TTREE**************************************************/
+
+ //Fill containing all the output of the networks to train
+
+ int NumOfBase=fActinideCompoInit.size();
+ for(int b=0;b<NumOfBase;b++)
+ {
+ //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@//
+ //@@@Change the input value according to your fresh fuel compo
+ //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@//
+ // (Z , A ,I)
+ U5 = fActinideCompoInit[b].GetZAIIsotopicQuantity(ZAI(92,235,0));
+ U8 = fActinideCompoInit[b].GetZAIIsotopicQuantity(ZAI(92,238,0));
+ Pu8 = fActinideCompoInit[b].GetZAIIsotopicQuantity(ZAI(94,238,0));
+ Pu9 = fActinideCompoInit[b].GetZAIIsotopicQuantity(ZAI(94,239,0));
+ Pu10 = fActinideCompoInit[b].GetZAIIsotopicQuantity(ZAI(94,240,0));
+ Pu11 = fActinideCompoInit[b].GetZAIIsotopicQuantity(ZAI(94,241,0));
+ Pu12 = fActinideCompoInit[b].GetZAIIsotopicQuantity(ZAI(94,242,0));
+ Am1 = fActinideCompoInit[b].GetZAIIsotopicQuantity(ZAI(95,241,0));
+ //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@//
+ //@@@Change the TStep according the step you want
+ //@@@ I have assumed all your depletion calculations have
+ //@@@ the same time binning
+ //@@@ hoping is what you have done ....
+ //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@//
+ int Tstep = 0; //at begining of cycle if 0
+
+ //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@//
+
+ k_eff=fkeff[b][Tstep];
+
+ fOutT->Fill();
+ }
+
+ fOutFile->Write();
+ delete fOutT;
+ fOutFile-> Close();
+ delete fOutFile;
+
+}
+//--------------------------------------------------------------------------------------------------
+void InitMass()
+{
+//Set tghe mass of the nuceli
+
+//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@//
+//@@@Change : ADD THE MASS OF THE NUCLEI PRESENT IN YOUR FRESH FUEL
+//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@//
+
+ ZAI U238 = ZAI(92,238,0);
+ ZAImass.insert(pair<ZAI,double>(U238, 238050788.247e-6));
+
+ ZAI U234 = ZAI(92,234,0);
+ ZAImass.insert(pair<ZAI,double>(U234, 234041000.000e-6));
+
+ ZAI U235 = ZAI(92,235,0);
+ ZAImass.insert(pair<ZAI,double>(U235, 235043929.918e-6));
+
+ ZAI Pu238 = ZAI(94,238,0);
+ ZAImass.insert(pair<ZAI,double>(Pu238,238049559.894e-6));
+
+ ZAI Pu239 = ZAI(94,239,0);
+ ZAImass.insert(pair<ZAI,double>(Pu239,239052163.381e-6));
+
+ ZAI Pu240 = ZAI(94,240,0);
+ ZAImass.insert(pair<ZAI,double>(Pu240,240053813.545e-6));
+
+ ZAI Pu241 = ZAI(94,241,0);
+ ZAImass.insert(pair<ZAI,double>(Pu241,241056851.456e-6));
+
+ ZAI Pu242 = ZAI(94,242,0);
+ ZAImass.insert(pair<ZAI,double>(Pu242,242058742.611e-6));
+
+ ZAI Am241 = ZAI(95,241,0);
+ ZAImass.insert(pair<ZAI,double>(Am241,241056829.144e-6));
+
+ ZAI O16 = ZAI(8,16,0);
+ ZAImass.insert(pair<ZAI,double>(O16, 0.));
+ //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@//
+
+}
+//--------------------------------------------------------------------------------------------------
+void CheckJob()
+{ //LOAD THE LIST OF EvolutionData
+ cout << "Scanning " << fEvolutionDataFolder << " for .dat files ..." << endl;
+ cout << "Please wait ..."<< endl;
+
+ string Command = "find "+ fEvolutionDataFolder + " -name \"*.dat\" > JOB.tmp";
+ system(Command.c_str());
+
+ ifstream JOB("JOB.tmp");
+ if (JOB.is_open())
+ {
+ while (!JOB.eof())
+ {
+ string tmp;
+ getline(JOB,tmp);
+ JobName.push_back(tmp);
+ }
+ } JOB.close(); JobName.pop_back();
+
+ // Remove temporary files...
+ Command = "\\rm -f JOB.tmp";
+ system(Command.c_str());
+ random_shuffle(JobName.begin(), JobName.end());
+ cout << "Scan complete" <<endl;
+}
+
+//--------------------------------------------------------------------------------------------------
+void ReadAndFill(string jobname)
+{ //Read a .dat file and fill XS maps and the fuel initial composition
+//Read a .dat file and fill XS maps and the fuel initial composition
+
+ vector<double> vT;
+ vector<int> Z;
+ vector<int> A;
+ vector<int> I;
+ vector<double> Q;
+
+ ifstream DecayDB(jobname.c_str()); // Open the File
+ if(!DecayDB)
+ {
+ cout << "!!Warning!! !!!EvolutiveProduct!!! \n Can't open \"" << jobname << "\"\n" << endl;
+ }
+
+ string line;
+ int start = 0;
+
+ getline(DecayDB, line);
+
+ /******Getting Time vecotr ....******/
+ if( StringLine::NextWord(line, start, ' ') != "time")
+ {
+ cout << "!!Bad Trouble!! !!!EvolutiveProduct!!! Bad Database file : " << jobname << endl;
+ exit (1);
+ }
+
+ while(start < (int)line.size())
+ vT.push_back(atof(StringLine::NextWord(line, start, ' ').c_str()));
+
+ fNOfTimeStep=int(vT.size());
+
+
+ start = 0;
+ getline(DecayDB, line);
+ string tmp = StringLine::NextWord(line, start, ' ');
+ vector<double> vKeff;
+ if ( tmp == "keff" || tmp == "Keff" )
+ {
+ while(start < (int)line.size())
+ vKeff.push_back(atof(StringLine::NextWord(line, start, ' ').c_str()));
+ }
+ fkeff.push_back(vKeff);
+
+
+ /****Getting Inventories***/
+ getline(DecayDB, line);
+ do
+ {
+ start = 0;
+ int z;
+ string tmp2 = StringLine::NextWord(line, start, ' ');
+ if (tmp2 == "Inv") {
+ z = atoi(StringLine::NextWord(line, start, ' ').c_str());
+ }
+ else z = atoi(tmp2.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);
+ Z.push_back(z);
+ A.push_back(a);
+ I.push_back(i);
+ if(!fIsAllNucleiAlreadyFill)
+ {
+ fAllNuclei.push_back(zaitmp);
+ fTime=vT;
+ }
+
+ long double q = atof(StringLine::NextWord(line, start, ' ').c_str());
+ Q.push_back(q);
+
+ }
+
+ getline(DecayDB, line);
+ start = 0;
+ tmp2 = StringLine::NextWord(line, start, ' ');
+
+ if(line == "" || line == "CrossSection" || tmp2 == "XSFis" || tmp2 == "XSCap" || tmp2 == "XSn2n") break;
+ }while (!DecayDB.eof() );
+
+ if(fAllNuclei.size()!=0)
+ fIsAllNucleiAlreadyFill=true;
+
+ DecayDB.close();
+
+ double N = 0;
+ for(int i=0; i < (int)Z.size()-2; i++)
+ {
+ if( Z[i]>89 )
+ N += Q[i];
+ }
+
+
+ IsotopicVector CompoBasei;
+
+ for(int i=0; i < (int)Z.size()-2; i++)
+ {
+ if(Z[i]>89)
+ {
+ ZAI zai = ZAI(Z[i], A[i], I[i]);
+ CompoBasei.IVquantity.insert(pair<ZAI,double>(zai,Q[i]/N));
+ }
+ }
+
+ fActinideCompoInit.push_back(CompoBasei);
+
+GoodJobName.push_back(jobname);
+
+
+}
+/*-------------------------------------------------------------------------------------------------
+COMPILATION :
+
+g++ -o GenerateRootFile GenerateRootFile.cxx `root-config --cflags` `root-config --libs`
+
+*/
\ No newline at end of file
diff --git a/Utils/trunk/EQM/FBR_MLP_Keff/Train_MLP.cxx b/Utils/trunk/EQM/FBR_MLP_Keff/Train_MLP.cxx
new file mode 100755
index 000000000..7ca138d39
--- /dev/null
+++ b/Utils/trunk/EQM/FBR_MLP_Keff/Train_MLP.cxx
@@ -0,0 +1,177 @@
+//
+// This program train and test a MLP using TMVA from a training data in a form of a TTRee
+//
+#include <cstdlib>
+#include <iostream>
+#include <map>
+#include <string>
+#include <sstream>
+#include "TChain.h"
+#include "TFile.h"
+#include "TTree.h"
+#include "TString.h"
+#include "TObjString.h"
+#include "TSystem.h"
+#include "TROOT.h"
+
+#include "TMVARegGui.C"
+
+#if not defined(__CINT__) || defined(__MAKECINT__)
+#include "TMVA/Tools.h"
+#include "TMVA/Factory.h"
+#endif
+
+using namespace TMVA;
+
+int main(int argc, char const *argv[])
+{
+ if(argc!=2)
+ {
+ std::cout << "TrainMLP Usage :" << std::endl;
+ std::cout << "\tTrainMLP YourTrainingData.root" << std::endl;
+ }
+ //---------------------------------------------------------------
+ // This loads the library
+ TMVA::Tools::Instance();
+ // ---------------------------------------------------------------
+ // if(argc !=3)
+ // {
+ // cout<<"error : "<<endl;
+ // cout<<"Usage : Train_MLP rootfile.root targetname : "<<endl;
+ // cout<<"\t targetname : keff or k_0 or k_1 ..."<<endl;
+ // exit(1);
+ // }
+ std::string TargetName = "k_eff";//argv[2];
+
+
+ std::cout << std::endl;
+ std::cout << "==> Start TMVARegression" << std::endl;
+
+
+ // Create a new root output file
+ TString outfileName( "TMVA_"+TargetName+".root" );
+ TFile* outputFile = TFile::Open( outfileName, "RECREATE" );//UPDATE
+
+
+ // The first argument is the base of the name of all the
+ // weightfiles in the directory weight/
+ //
+ // The second argument is the output file for the training results
+ // All TMVA output can be suppressed by removing the "!" (not) in
+ // front of the "Silent" argument in the option string
+
+ TMVA::Factory *factory = new TMVA::Factory( "TMVARegression", outputFile,
+ "!V:!Silent:Color:DrawProgressBar" );
+
+ // If you wish to modify default settings
+ // (please check "src/Config.h" to see all available global options)
+ // (TMVA::gConfig().GetVariablePlotting()).fTimesRMS = 8.0;
+ // (TMVA::gConfig().GetIONames()).fWeightFileDir = "myWeightDirectory";
+
+ // Define the input variables that shall be used for the TMVA training
+ // note that you may also use variable expressions, such as: "3*var1/var2*abs(var3)"
+ // [all types of expressions that can also be parsed by TTree::Draw( "expression" )]
+
+ //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@//
+ //@@@Change/add variables if needed
+ //@@@ the variable name are the one you use in GenerateRootFile.cxx
+ //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@//
+ factory->AddVariable( "U8" , "U 238 " , "FractionIsotopic", 'F' );
+ factory->AddVariable( "U5" , "U 235 " , "FractionIsotopic", 'F' );
+ factory->AddVariable( "Pu8" , "Pu 238 " , "FractionIsotopic", 'F' );
+ factory->AddVariable( "Pu9" , "Pu 239 " , "FractionIsotopic", 'F' );
+ factory->AddVariable( "Pu10" , "Pu 240 " , "FractionIsotopic", 'F' );
+ factory->AddVariable( "Pu11" , "Pu 241 " , "FractionIsotopic", 'F' );
+ factory->AddVariable( "Pu12" , "Pu 242 " , "FractionIsotopic", 'F' );
+ factory->AddVariable( "Am1" , "Am 241 " , "FractionIsotopic", 'F' );
+
+
+ // You can add so-called "Spectator variables", which are not used in the MVA training,
+ // but will appear in the final "TestTree" produced by TMVA. This TestTree will contain the
+ // input variables, the response values of all trained MVAs, and the spectator variables
+ //factory->AddSpectator( "spec1:=var1*2", "Spectator 1", "units", 'F' );
+ //factory->AddSpectator( "spec2:=var1*3", "Spectator 2", "units", 'F' );
+
+ // Add the variable carrying the regression target
+
+ factory->AddTarget(TargetName);
+
+ // It is also possible to declare additional targets for multi-dimensional regression, ie:
+ // -- factory->AddTarget( "fvalue2" );
+ // BUT: this is currently ONLY implemented for MLP
+
+ // Read training and test data (see TMVAClassification for reading ASCII files)
+ // load the signal and background event samples from ROOT trees
+ TFile *input(0);
+ TString fname = argv[1]; //Training Data input file
+ if (!gSystem->AccessPathName( fname ))
+ input = TFile::Open( fname ); // check if file in local directory exists
+
+ if (!input) {
+ std::cout << "ERROR: could not open data file" << std::endl;
+ exit(1);
+ }
+ std::cout << "--- TMVARegression : Using input file: " << input->GetName() << std::endl;
+
+ // --- Register the regression tree
+
+ TTree *regTree = (TTree*)input->Get("Data");
+
+ // global event weights per tree (see below for setting event-wise weights)
+ Double_t regWeight = 1.0;
+
+ // You can add an arbitrary number of regression trees
+ factory->AddRegressionTree( regTree, regWeight );
+
+ // This would set individual event weights (the variables defined in the
+ // expression need to exist in the original TTree)
+ // factory->SetWeightExpression( "var1", "Regression" );
+
+ // Apply additional cuts on the signal and background samples (can be different)
+ TCut mycut = ""; // for example: TCut mycut = "abs(var1)<0.5 && abs(var2-0.5)<1";
+
+ // tell the factory to use all remaining events in the trees after training for testing:
+ factory->PrepareTrainingAndTestTree( mycut,
+ "nTrain_Regression=0:nTest_Regression=0:SplitMode=Random:NormMode=NumEvents:!V" );
+
+ // If no numbers of events are given, half of the events in the tree are used
+ // for training, and the other half for testing:
+
+ // ---- Book MVA methods
+ //
+ std::stringstream Name;
+ Name<<"MLP_"<<TargetName;
+ // Neural network (MLP)
+ factory->BookMethod( TMVA::Types::kMLP, Name.str().c_str(), "!H:V:VarTransform=Norm:NeuronType=tanh:NCycles=16000:HiddenLayers=8,:TestRate=6:TrainingMethod=BFGS:Sampling=0.3:SamplingEpoch=0.8:ConvergenceImprove=1e-6:ConvergenceTests=15:!UseRegulator" );
+
+
+ // --------------------------------------------------------------------------------------------------
+
+ // ---- Now you can tell the factory to train, test, and evaluate the MVAs
+
+ // Train MVAs using the set of training events
+ factory->TrainAllMethods();
+
+ // ---- Evaluate all MVAs using the set of test events
+ factory->TestAllMethods();
+
+ // ----- Evaluate and compare performance of all configured MVAs
+ factory->EvaluateAllMethods();
+
+ // --------------------------------------------------------------
+
+ // Save the output
+ outputFile->Close();
+
+ std::cout << "==> Wrote root file: " << outputFile->GetName() << std::endl;
+ std::cout << "==> TMVARegression is done!" << std::endl;
+
+ delete factory;
+
+ // Launch the GUI for the root macros
+if (!gROOT->IsBatch()) TMVARegGui( outfileName );
+}
+/*
+
+g++ -o Train_MLP `root-config --cflags` Train_MLP.cxx `root-config --glibs` -lTMVA -I$ROOTSYS/tmva/test/
+*/
\ No newline at end of file
diff --git a/Utils/trunk/EQM/FBR_MLP_Keff/include/GenerateRootFile.hxx b/Utils/trunk/EQM/FBR_MLP_Keff/include/GenerateRootFile.hxx
new file mode 100755
index 000000000..ef7b9ff21
--- /dev/null
+++ b/Utils/trunk/EQM/FBR_MLP_Keff/include/GenerateRootFile.hxx
@@ -0,0 +1,51 @@
+#include "ZAI.hxx"
+#include <TGraph.h>
+#include <TGraph.h>
+
+#include <vector>
+#include <stdio.h>
+#include <stdlib.h>
+#include <iostream>
+#include <fstream>
+#include <sstream>
+#include <string>
+#include <vector>
+#include <sstream>
+#include <cmath>
+#include <map>
+#include <iostream>
+#include <iomanip>
+using namespace std;
+
+string dtoa(double num)
+{
+ ostringstream os(ostringstream::out);
+ os<<setprecision(3)<<num;
+ return os.str();
+}
+
+string FilePath;
+string DataPath;
+
+map<ZAI, double> ZAImass;
+vector<string> JobName;
+vector<string> GoodJobName;
+
+vector<double> fTime; //Time vector of the depletion calculation (second)
+vector< vector<double> > fkeff;
+
+vector<ZAI> fAllNuclei; //All the nuclei present in the fuel
+
+
+vector<IsotopicVector> fActinideCompoInit; //Fresh fuel composition
+
+int fNOfTimeStep=0; //number of time step in the Evolution
+
+string fEvolutionDataFolder = "";
+
+bool fIsAllNucleiAlreadyFill=false;
+
+void InitMass();
+void CheckJob();
+void ReadAndFill(string jobname);
+void DumpInputNeuron(string filename);
\ No newline at end of file
diff --git a/Utils/trunk/EQM/FBR_MLP_Keff/include/StringLine.hxx b/Utils/trunk/EQM/FBR_MLP_Keff/include/StringLine.hxx
new file mode 100755
index 000000000..da8500370
--- /dev/null
+++ b/Utils/trunk/EQM/FBR_MLP_Keff/include/StringLine.hxx
@@ -0,0 +1,275 @@
+#ifndef _STRINGLINE_
+#define _STRINGLINE_
+
+#include <string>
+#include <sstream>
+#include <iostream>
+#include <algorithm>
+#include <cctype>
+using namespace std;
+/*!
+ \file
+ \brief Header file for StingLine class.
+*/
+
+// Class extracting fields from a string / line.
+/*!
+ The aim of this class is to provide tools to extract fields ("word") from
+ a string and convert a string in Upper/Lower case.
+ All methods are static so that it is not necessary to create object to use them
+
+ example:
+ \code
+ string line="The temperature is : 300.6 K";
+ int start;
+
+ 1st method: creation of StringLine
+
+ start=0;
+ StringLine SL;
+ string the=SL.NextWord(line,start);
+ string temperature_is=SL.NextWord(line,start,':');
+ string colon=SL.NextWord(line,start);
+ double T=atof(SL.NextWord(line,start).c_str());
+ cout<<the<<endl<<temperature_is<<endl<<T<<endl;
+
+ 2nd method: "using" the static methods
+
+ start=0;
+ the=StringLine::NextWord(line,start);
+ temperature_is=StringLine::NextWord(line,start,':');
+ colon=StringLine::NextWord(line,start);
+ T=atof(StringLine::NextWord(line,start).c_str());
+ cout<<the<<endl<<temperature_is<<endl<<T<<endl;
+ \endcode
+ @author PTO
+ @version 0.1
+*/
+
+class StringLine
+{
+ public:
+ // Find the next word in a line.
+ /*!
+ Find Next word in a line starting from position "start" in the line. If an alternative
+ separator is given, the word length is defined by the first position of sep or alt_sep found.
+ The first value of start is in general 0 (i.e. the beginning of the Line)
+ \param Line : a line containing words
+ \param start : from where to start to find the begining of a word
+ \param sep : the separator between 2 words (default=space)
+ \param alt_sep : the alternative separator between 2 words (default='')
+ */
+ static string NextWord(string Line,int &start,char sep=' ', char alt_sep='\0');
+ // Find the previous word in a line.
+ /*!
+ Find Previous word in a line starting from position "start" in the line. If an alternative
+ separator is given, the word length is defined by the first position of sep or alt_sep found.
+ The first value of start is in general the end of the Line.
+ \param Line : a line containing words
+ \param start : from where to start to find the begining of a word
+ \param sep : the separator between 2 words (default=space)
+ \param alt_sep : the alternative separator between 2 words (default='')
+ */
+ static string PreviousWord(string Line,int &start,char sep=' ', char alt_sep='\0');
+ static void ToLower(string &Line); // convert a string to Lower case
+ static void ToUpper(string &Line); // convert a string to Upper case
+
+ // Find \p search in \p Line from the begining.
+ /*!
+ returns the position, starting from the begenning of the first occurence
+ of \p search in \p Line if it is found, else returns -1
+ \param search : a string to find
+ \param Line : where to search
+ */
+ static int Find(const char *search,string Line);
+ // Find \p search in \p Line from the end.
+ /*!
+ returns the position, starting from the end of the first occurence
+ of \p search in \p Line if it is found, else returns -1
+ \param search : a string to find
+ \param Line : where to search
+ */
+ static int rFind(const char *search,string Line);
+ // convert a input type (\p in_T) to another (\p out_T).
+ /*!
+ Example:
+ \code
+ string s="32.12";
+ double t=StringLine::convert<double>(s);
+ string temperature=StringLine::convert<string>(300.);
+ \endcode
+ \param t : the input value
+ */
+ template <class out_T, class in_T> static out_T convert(const in_T & t);
+ // Find the start of a word in a line.
+ /*!
+ \param Line : a line containing words
+ \param CurrentPosition : from where to start to find the begining of a word
+ \param sep : the separator between 2 words (default=space)
+ \param alt_sep : the alternative separator between 2 words (default='')
+ */
+ static int GetStartWord(string Line,int CurrentPosition,char sep=' ', char alt_sep='\0');
+ // Find the end of a word in a line.
+ /*!
+ \param Line : a line containing words
+ \param CurrentPosition : from where to start to find the end of a word
+ \param sep : the separator between 2 words (default=space)
+ \param alt_sep : the alternative separator between 2 words (default='')
+ */
+ static int GetEndWord(string Line,int CurrentPosition,char sep=' ', char alt_sep='\0');
+ // Replace a sub-string by an other in a string.
+ /*!
+ \param InLine : the string which contains the sub-string to replace
+ \param ToReplace : the sub-string to replace
+ \param By : the sub-string ToReplace is replaced by the sub-string By in Inline
+ */
+ string ReplaceAll(string InLine, string ToReplace, string By);
+};
+
+
+//_________________________________________________________________________________
+inline string StringLine::NextWord(string Line,int &start,char sep, char alt_sep)
+{
+ string Word="";
+ if(start>=int(Line.size()))
+ {
+ return Word;
+ }
+ start=GetStartWord(Line,start,sep,alt_sep);
+ int wordlength=GetEndWord(Line,start,sep,alt_sep)-start;
+
+ Word=Line.substr(start,wordlength);
+
+ start+=wordlength;
+ return Word;
+}
+//_________________________________________________________________________________
+inline string StringLine::PreviousWord(string Line,int &start,char sep, char alt_sep)
+{
+ string Word="";
+ if(start<=0)
+ {
+ return Word;
+ }
+ int pos=Line.rfind(sep,start);
+ int alt_pos=-1;
+ int real_pos=pos;
+ char real_sep=sep;
+ if(alt_sep!='\0')
+ {
+ alt_pos=Line.rfind(alt_sep,start);
+ real_pos=max(pos,alt_pos);
+ if(real_pos!=pos)
+ real_sep=alt_sep;
+ }
+ int wordlength=start-Line.rfind(real_sep,real_pos);
+ if(real_pos<=0)
+ {
+ Word=Line.substr(0,start+1);
+ start=0;
+ return Word;
+ }
+ Word=Line.substr(real_pos+1,wordlength);
+
+ start-=wordlength+1;
+ return Word;
+}
+
+//_________________________________________________________________________________
+inline void StringLine::ToLower(string &Line)
+{
+ transform (Line.begin(), Line.end(), // source
+ Line.begin(), // destination
+ (int(*)(int))tolower); // operation
+}
+
+//_________________________________________________________________________________
+inline void StringLine::ToUpper(string &Line)
+{
+ transform (Line.begin(), Line.end(), // source
+ Line.begin(), // destination
+ (int(*)(int))toupper); // operation
+}
+
+//_________________________________________________________________________________
+inline int StringLine::GetStartWord(string Line,int CurrentPosition,char sep, char alt_sep)
+{
+ int pos=Line.find(sep,CurrentPosition);
+ int alt_pos=-1;
+ if(alt_sep!='\0')
+ alt_pos=Line.find(alt_sep,CurrentPosition);
+ int real_pos=pos;
+ char real_sep=sep;
+ if(alt_pos>=0)
+ {
+ real_pos=min(pos,alt_pos);
+ if(pos==int(string::npos))real_pos=alt_pos;
+ if(real_pos!=pos)
+ real_sep=alt_sep;
+ }
+ if(real_pos==int(string::npos)) return CurrentPosition;
+ while(CurrentPosition<int(Line.size()) && Line[CurrentPosition]==real_sep)
+ CurrentPosition++;
+ return CurrentPosition;
+}
+
+//_________________________________________________________________________________
+inline int StringLine::GetEndWord(string Line,int CurrentPosition,char sep, char alt_sep)
+{
+ int pos=Line.find(sep,CurrentPosition);
+ int alt_pos=-1;
+ if(alt_sep!='\0')
+ alt_pos=Line.find(alt_sep,CurrentPosition);
+ int real_pos=pos;
+ if(alt_pos>=0)
+ {
+ real_pos=min(pos,alt_pos);
+ if(pos==int(string::npos))real_pos=alt_pos;
+ }
+ if(real_pos==int(string::npos))
+ return Line.size();
+ return real_pos;
+}
+
+//_________________________________________________________________________________
+inline int StringLine::Find(const char *search,string Line)
+{
+ size_t Pos=Line.find(search);
+ if(Pos != string::npos ) return Pos;
+ return -1;
+}
+
+//_________________________________________________________________________________
+inline int StringLine::rFind(const char *search,string Line)
+{
+ size_t Pos=Line.rfind(search);
+ if(Pos != string::npos) return Pos;
+ return -1;
+}
+
+//_________________________________________________________________________________
+template <class out_T, class in_T>
+inline out_T StringLine::convert(const in_T & t)
+{
+ stringstream stream;
+ stream << t; // insert value to stream
+ out_T result; // store conversion's result here
+ stream >> result; // write value to result
+ return result;
+}
+
+//_________________________________________________________________________________
+inline string StringLine::ReplaceAll(string InLine, string ToReplace, string By)
+{
+ int start=0;
+ int pos=InLine.find(ToReplace,start);
+ while(pos!=int(string::npos))
+ {
+ InLine.replace(pos,ToReplace.size(),By);
+ start=0;
+ pos=InLine.find(ToReplace,start);
+ }
+ return InLine;
+
+}
+#endif
diff --git a/Utils/trunk/EQM/FBR_MLP_Keff/include/ZAI.hxx b/Utils/trunk/EQM/FBR_MLP_Keff/include/ZAI.hxx
new file mode 100755
index 000000000..1ee9136f9
--- /dev/null
+++ b/Utils/trunk/EQM/FBR_MLP_Keff/include/ZAI.hxx
@@ -0,0 +1,42 @@
+#ifndef MN_ZAI_H_
+#define MN_ZAI_H_
+
+#include <map>
+using namespace std;
+
+class ZAI {
+ public :
+ int Z;
+ int A;
+ int I;
+
+
+ ZAI(int z, int a, int i=0) { Z=z;A=a; I=i;}
+
+ bool operator <(const ZAI& zai) const {return (Z != zai.Z)?
+ (Z < zai.Z) : ( (A != zai.A)?
+ (A < zai.A) : (I < zai.I) );}
+};
+
+class IsotopicVector {
+ public :
+
+ IsotopicVector() {} ///< Normal Constructor.
+ ~IsotopicVector() {}; ///< Normal Destructor.
+ double GetZAIIsotopicQuantity(const ZAI& zai) const
+ {
+ map<ZAI ,double> IsotopicQuantity = IVquantity;
+
+ map<ZAI ,double>::iterator it;
+ it = IsotopicQuantity.find(zai);
+
+ if ( it != IsotopicQuantity.end() )
+ return it->second;
+ else
+ return 0;
+ }
+ map<ZAI,double> IVquantity;
+ double MASS;
+};
+
+#endif //MN_ZAI_H_
diff --git a/Utils/trunk/EQM/MLP_Kinf/GenerateRootFile.cxx b/Utils/trunk/EQM/MLP_Kinf/GenerateRootFile.cxx
new file mode 100755
index 000000000..8bd2d961c
--- /dev/null
+++ b/Utils/trunk/EQM/MLP_Kinf/GenerateRootFile.cxx
@@ -0,0 +1,343 @@
+/**********************************************************/
+// Make the input file for the MLPs training
+//
+// This programs reads a set of .dat files which are the
+// results of a depletion calculation (see manual and
+// looks for XS_CLOSEST). From the reading it fills a
+// TTree (Data) and write it in a file named
+// TrainingInput.root .
+// The file TrainingInput.cxx is the list of MLP outputs
+// (cross sections)
+//
+//@author BaM, BaL
+/**********************************************************/
+#include "include/GenerateRootFile.hxx"
+#include <TH1F.h>
+#include <TH2D.h>
+#include <TFile.h>
+#include <TTree.h>
+#include "include/StringLine.hxx"
+#include <TString.h>
+#include <string>
+#include <cmath>
+#include <iostream>
+#include <fstream>
+#include <algorithm>
+#include <map>
+#include <sstream>
+
+using namespace std;
+//--------------------------------------------------------------------------------------------------
+/*************************
+ MAIN
+*************************/
+int main(int argc, char ** argv){
+
+ if(argc!=2)
+ {
+ cout << "Usage : BuildInputTree Path" << endl;
+ cout << " Where Path is the path to the folder containing the Evolution Datas" << endl;
+ cout << " i.e the (.dat) files" << endl;
+ exit(0);
+ }
+
+ fEvolutionDataFolder = argv[1];
+
+ InitMass(); //Load nuclei masses
+ CheckJob(); // looks fot the .dat files in the fEvolutionDataFolder
+
+ FilePath = "DB_TMP/";
+ DataPath = FilePath + "Data/";
+ string Command = "mkdir -p " + DataPath;
+ system(Command.c_str());
+
+ cout << "Reading .dat files ..." << endl;
+ for(int i = 0; i < (int)JobName.size(); i++)
+ {
+ ReadAndFill(JobName[i]);
+ if (i%100 == 0)
+ cout << "\r" << i << " .dat files read" <<flush;
+ }
+ cout << "Filling the TTree ..." << endl;
+ DumpInputNeuron("TrainingInput.root");
+ cout << "Training input generated in file : TrainingInput.root " << endl;
+ cout << "Names of MLP outputs in file : TrainingInput.cxx " << endl;
+ system("rm -r DB_TMP");
+
+}
+//--------------------------------------------------------------------------------------------------
+// Function definitions
+//--------------------------------------------------------------------------------------------------
+//Convert int to string
+string itoa(int num)
+{
+ ostringstream os(ostringstream::out);
+ os << num;
+ return os.str();
+}
+//--------------------------------------------------------------------------------------------------
+void DumpInputNeuron(string filename)
+{
+ TFile* fOutFile = new TFile(filename.c_str(),"RECREATE");
+ TTree* fOutT = new TTree("Data", "Data");
+
+/**********************INITIALISATIONNN********************/
+ //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@//
+ //@@@Change the input value according to your fresh fuel compo
+ //-> here MOX FUEL
+ //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@//
+ double U5 = 0;
+ double U8 = 0;
+ double Pu8 = 0;
+ double Pu9 = 0;
+ double Pu10 = 0;
+ double Pu11 = 0;
+ double Pu12 = 0;
+ double Am1 = 0;
+
+ //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@//
+ double Time = 0;
+ double k_inf = 0;
+
+
+/**********************BRANCHING**************************************************/
+/**********************Fresh fuel**************************************************/
+ //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@//
+ //@@@Change the input value according to your fresh fuel compo
+ //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@//
+ fOutT->Branch( "U5" ,&U5 ,"U5/D" );
+ fOutT->Branch( "U8" ,&U8 ,"U8/D" );
+ fOutT->Branch( "Pu8" ,&Pu8 ,"Pu8/D" );
+ fOutT->Branch( "Pu9" ,&Pu9 ,"Pu9/D" );
+ fOutT->Branch( "Pu10" ,&Pu10 ,"Pu10/D" );
+ fOutT->Branch( "Pu11" ,&Pu11 ,"Pu11/D" );
+ fOutT->Branch( "Pu12" ,&Pu12 ,"Pu12/D" );
+ fOutT->Branch( "Am1" ,&Am1 ,"Am1/D" );
+ //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@//
+ fOutT->Branch( "Time" ,&Time ,"Time/D" );
+ fOutT->Branch( "k_inf" ,&k_inf ,"k_inf/D" );
+
+/**********************FILLING THE TTREE**************************************************/
+
+ //Fill containing all the output of the networks to train
+
+ int NumOfBase=fActinideCompoInit.size();
+ for(int b=0;b<NumOfBase;b++)
+ {
+ //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@//
+ //@@@Change the input value according to your fresh fuel compo
+ //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@//
+ // (Z , A ,I)
+ U5 = fActinideCompoInit[b].GetZAIIsotopicQuantity(ZAI(92,235,0));
+ U8 = fActinideCompoInit[b].GetZAIIsotopicQuantity(ZAI(92,238,0));
+ Pu8 = fActinideCompoInit[b].GetZAIIsotopicQuantity(ZAI(94,238,0));
+ Pu9 = fActinideCompoInit[b].GetZAIIsotopicQuantity(ZAI(94,239,0));
+ Pu10 = fActinideCompoInit[b].GetZAIIsotopicQuantity(ZAI(94,240,0));
+ Pu11 = fActinideCompoInit[b].GetZAIIsotopicQuantity(ZAI(94,241,0));
+ Pu12 = fActinideCompoInit[b].GetZAIIsotopicQuantity(ZAI(94,242,0));
+ Am1 = fActinideCompoInit[b].GetZAIIsotopicQuantity(ZAI(95,241,0));
+ //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@//
+ for(int Tstep=0 ;Tstep<fNOfTimeStep;Tstep++ )
+ {
+ Time=fTime[Tstep];
+ k_inf=fkeff[b][Tstep];
+ fOutT->Fill();
+ }
+ }
+
+ fOutFile->Write();
+ delete fOutT;
+ fOutFile-> Close();
+ delete fOutFile;
+
+}
+//--------------------------------------------------------------------------------------------------
+void InitMass()
+{
+//Set tghe mass of the nuceli
+
+//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@//
+//@@@Change : ADD THE MASS OF THE NUCLEI PRESENT IN YOUR FRESH FUEL
+//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@//
+
+ ZAI U238 = ZAI(92,238,0);
+ ZAImass.insert(pair<ZAI,double>(U238, 238050788.247e-6));
+
+ ZAI U234 = ZAI(92,234,0);
+ ZAImass.insert(pair<ZAI,double>(U234, 234041000.000e-6));
+
+ ZAI U235 = ZAI(92,235,0);
+ ZAImass.insert(pair<ZAI,double>(U235, 235043929.918e-6));
+
+ ZAI Pu238 = ZAI(94,238,0);
+ ZAImass.insert(pair<ZAI,double>(Pu238,238049559.894e-6));
+
+ ZAI Pu239 = ZAI(94,239,0);
+ ZAImass.insert(pair<ZAI,double>(Pu239,239052163.381e-6));
+
+ ZAI Pu240 = ZAI(94,240,0);
+ ZAImass.insert(pair<ZAI,double>(Pu240,240053813.545e-6));
+
+ ZAI Pu241 = ZAI(94,241,0);
+ ZAImass.insert(pair<ZAI,double>(Pu241,241056851.456e-6));
+
+ ZAI Pu242 = ZAI(94,242,0);
+ ZAImass.insert(pair<ZAI,double>(Pu242,242058742.611e-6));
+
+ ZAI Am241 = ZAI(95,241,0);
+ ZAImass.insert(pair<ZAI,double>(Am241,241056829.144e-6));
+
+ ZAI O16 = ZAI(8,16,0);
+ ZAImass.insert(pair<ZAI,double>(O16, 0.));
+ //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@//
+
+}
+//--------------------------------------------------------------------------------------------------
+void CheckJob()
+{ //LOAD THE LIST OF EvolutionData
+ cout << "Scanning " << fEvolutionDataFolder << " for .dat files ..." << endl;
+ cout << "Please wait ..."<< endl;
+
+ string Command = "find "+ fEvolutionDataFolder + " -name \"*.dat\" > JOB.tmp";
+ system(Command.c_str());
+
+ ifstream JOB("JOB.tmp");
+ if (JOB.is_open())
+ {
+ while (!JOB.eof())
+ {
+ string tmp;
+ getline(JOB,tmp);
+ JobName.push_back(tmp);
+ }
+ } JOB.close(); JobName.pop_back();
+
+ // Remove temporary files...
+ Command = "\\rm -f JOB.tmp";
+ system(Command.c_str());
+ random_shuffle(JobName.begin(), JobName.end());
+ cout << "Scan complete" <<endl;
+}
+
+//--------------------------------------------------------------------------------------------------
+void ReadAndFill(string jobname)
+{ //Read a .dat file and fill XS maps and the fuel initial composition
+//Read a .dat file and fill XS maps and the fuel initial composition
+
+ vector<double> vT;
+ vector<int> Z;
+ vector<int> A;
+ vector<int> I;
+ vector<double> Q;
+
+ ifstream DecayDB(jobname.c_str()); // Open the File
+ if(!DecayDB)
+ {
+ cout << "!!Warning!! !!!EvolutiveProduct!!! \n Can't open \"" << jobname << "\"\n" << endl;
+ }
+
+ string line;
+ int start = 0;
+
+ getline(DecayDB, line);
+
+ /******Getting Time vecotr ....******/
+ if( StringLine::NextWord(line, start, ' ') != "time")
+ {
+ cout << "!!Bad Trouble!! !!!EvolutiveProduct!!! Bad Database file : " << jobname << endl;
+ exit (1);
+ }
+
+ while(start < (int)line.size())
+ vT.push_back(atof(StringLine::NextWord(line, start, ' ').c_str()));
+
+ fNOfTimeStep=int(vT.size());
+
+
+ start = 0;
+ getline(DecayDB, line);
+ string tmp = StringLine::NextWord(line, start, ' ');
+ vector<double> vKeff;
+ if ( tmp == "keff" || tmp == "Keff" )
+ {
+ while(start < (int)line.size())
+ vKeff.push_back(atof(StringLine::NextWord(line, start, ' ').c_str()));
+ }
+ fkeff.push_back(vKeff);
+
+
+ /****Getting Inventories***/
+ getline(DecayDB, line);
+ do
+ {
+ start = 0;
+ int z;
+ string tmp2 = StringLine::NextWord(line, start, ' ');
+ if (tmp2 == "Inv") {
+ z = atoi(StringLine::NextWord(line, start, ' ').c_str());
+ }
+ else z = atoi(tmp2.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);
+ Z.push_back(z);
+ A.push_back(a);
+ I.push_back(i);
+ if(!fIsAllNucleiAlreadyFill)
+ {
+ fAllNuclei.push_back(zaitmp);
+ fTime=vT;
+ }
+
+ long double q = atof(StringLine::NextWord(line, start, ' ').c_str());
+ Q.push_back(q);
+
+ }
+
+ getline(DecayDB, line);
+ start = 0;
+ tmp2 = StringLine::NextWord(line, start, ' ');
+
+ if(line == "" || line == "CrossSection" || tmp2 == "XSFis" || tmp2 == "XSCap" || tmp2 == "XSn2n") break;
+ }while (!DecayDB.eof() );
+
+ if(fAllNuclei.size()!=0)
+ fIsAllNucleiAlreadyFill=true;
+
+ DecayDB.close();
+
+ double N = 0;
+ for(int i=0; i < (int)Z.size()-2; i++)
+ {
+ if( Z[i]>89 )
+ N += Q[i];
+ }
+
+
+ IsotopicVector CompoBasei;
+
+ for(int i=0; i < (int)Z.size()-2; i++)
+ {
+ if(Z[i]>89)
+ {
+ ZAI zai = ZAI(Z[i], A[i], I[i]);
+ CompoBasei.IVquantity.insert(pair<ZAI,double>(zai,Q[i]/N));
+ }
+ }
+
+ fActinideCompoInit.push_back(CompoBasei);
+
+GoodJobName.push_back(jobname);
+
+
+}
+/*-------------------------------------------------------------------------------------------------
+COMPILATION :
+
+g++ -o GenerateRootFile GenerateRootFile.cxx `root-config --cflags` `root-config --libs`
+
+*/
\ No newline at end of file
diff --git a/Utils/trunk/EQM/MLP_Kinf/Train_MLP.cxx b/Utils/trunk/EQM/MLP_Kinf/Train_MLP.cxx
new file mode 100755
index 000000000..8f8aa47e8
--- /dev/null
+++ b/Utils/trunk/EQM/MLP_Kinf/Train_MLP.cxx
@@ -0,0 +1,178 @@
+//
+// This program train and test a MLP using TMVA from a training data in a form of a TTRee
+//
+#include <cstdlib>
+#include <iostream>
+#include <map>
+#include <string>
+#include <sstream>
+#include "TChain.h"
+#include "TFile.h"
+#include "TTree.h"
+#include "TString.h"
+#include "TObjString.h"
+#include "TSystem.h"
+#include "TROOT.h"
+
+#include "TMVARegGui.C"
+
+#if not defined(__CINT__) || defined(__MAKECINT__)
+#include "TMVA/Tools.h"
+#include "TMVA/Factory.h"
+#endif
+
+using namespace TMVA;
+
+int main(int argc, char const *argv[])
+{
+ if(argc!=2)
+ {
+ std::cout << "TrainMLP Usage :" << std::endl;
+ std::cout << "\tTrainMLP YourTrainingData.root" << std::endl;
+ }
+ //---------------------------------------------------------------
+ // This loads the library
+ TMVA::Tools::Instance();
+ // ---------------------------------------------------------------
+ // if(argc !=3)
+ // {
+ // cout<<"error : "<<endl;
+ // cout<<"Usage : Train_MLP rootfile.root targetname : "<<endl;
+ // cout<<"\t targetname : keff or k_0 or k_1 ..."<<endl;
+ // exit(1);
+ // }
+ std::string TargetName = "k_inf";//argv[2];
+
+
+ std::cout << std::endl;
+ std::cout << "==> Start TMVARegression" << std::endl;
+
+
+ // Create a new root output file
+ TString outfileName( "TMVA_"+TargetName+".root" );
+ TFile* outputFile = TFile::Open( outfileName, "RECREATE" );//UPDATE
+
+
+ // The first argument is the base of the name of all the
+ // weightfiles in the directory weight/
+ //
+ // The second argument is the output file for the training results
+ // All TMVA output can be suppressed by removing the "!" (not) in
+ // front of the "Silent" argument in the option string
+
+ TMVA::Factory *factory = new TMVA::Factory( "TMVARegression", outputFile,
+ "!V:!Silent:Color:DrawProgressBar" );
+
+ // If you wish to modify default settings
+ // (please check "src/Config.h" to see all available global options)
+ // (TMVA::gConfig().GetVariablePlotting()).fTimesRMS = 8.0;
+ // (TMVA::gConfig().GetIONames()).fWeightFileDir = "myWeightDirectory";
+
+ // Define the input variables that shall be used for the TMVA training
+ // note that you may also use variable expressions, such as: "3*var1/var2*abs(var3)"
+ // [all types of expressions that can also be parsed by TTree::Draw( "expression" )]
+
+ //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@//
+ //@@@Change/add variables if needed
+ //@@@ the variable name are the one you use in GenerateRootFile.cxx
+ //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@//
+ factory->AddVariable( "U8" , "U 238 " , "FractionIsotopic", 'F' );
+ factory->AddVariable( "U5" , "U 235 " , "FractionIsotopic", 'F' );
+ factory->AddVariable( "Pu8" , "Pu 238 " , "FractionIsotopic", 'F' );
+ factory->AddVariable( "Pu9" , "Pu 239 " , "FractionIsotopic", 'F' );
+ factory->AddVariable( "Pu10" , "Pu 240 " , "FractionIsotopic", 'F' );
+ factory->AddVariable( "Pu11" , "Pu 241 " , "FractionIsotopic", 'F' );
+ factory->AddVariable( "Pu12" , "Pu 242 " , "FractionIsotopic", 'F' );
+ factory->AddVariable( "Am1" , "Am 241 " , "FractionIsotopic", 'F' );
+ factory->AddVariable( "Time" , "Time" , "temps(s)", 'F' );
+
+
+ // You can add so-called "Spectator variables", which are not used in the MVA training,
+ // but will appear in the final "TestTree" produced by TMVA. This TestTree will contain the
+ // input variables, the response values of all trained MVAs, and the spectator variables
+ //factory->AddSpectator( "spec1:=var1*2", "Spectator 1", "units", 'F' );
+ //factory->AddSpectator( "spec2:=var1*3", "Spectator 2", "units", 'F' );
+
+ // Add the variable carrying the regression target
+
+ factory->AddTarget(TargetName);
+
+ // It is also possible to declare additional targets for multi-dimensional regression, ie:
+ // -- factory->AddTarget( "fvalue2" );
+ // BUT: this is currently ONLY implemented for MLP
+
+ // Read training and test data (see TMVAClassification for reading ASCII files)
+ // load the signal and background event samples from ROOT trees
+ TFile *input(0);
+ TString fname = argv[1]; //Training Data input file
+ if (!gSystem->AccessPathName( fname ))
+ input = TFile::Open( fname ); // check if file in local directory exists
+
+ if (!input) {
+ std::cout << "ERROR: could not open data file" << std::endl;
+ exit(1);
+ }
+ std::cout << "--- TMVARegression : Using input file: " << input->GetName() << std::endl;
+
+ // --- Register the regression tree
+
+ TTree *regTree = (TTree*)input->Get("Data");
+
+ // global event weights per tree (see below for setting event-wise weights)
+ Double_t regWeight = 1.0;
+
+ // You can add an arbitrary number of regression trees
+ factory->AddRegressionTree( regTree, regWeight );
+
+ // This would set individual event weights (the variables defined in the
+ // expression need to exist in the original TTree)
+ // factory->SetWeightExpression( "var1", "Regression" );
+
+ // Apply additional cuts on the signal and background samples (can be different)
+ TCut mycut = ""; // for example: TCut mycut = "abs(var1)<0.5 && abs(var2-0.5)<1";
+
+ // tell the factory to use all remaining events in the trees after training for testing:
+ factory->PrepareTrainingAndTestTree( mycut,
+ "nTrain_Regression=0:nTest_Regression=0:SplitMode=Random:NormMode=NumEvents:!V" );
+
+ // If no numbers of events are given, half of the events in the tree are used
+ // for training, and the other half for testing:
+
+ // ---- Book MVA methods
+ //
+ std::stringstream Name;
+ Name<<"MLP_"<<TargetName;
+ // Neural network (MLP)
+ factory->BookMethod( TMVA::Types::kMLP, Name.str().c_str(), "!H:V:VarTransform=Norm:NeuronType=tanh:NCycles=16000:HiddenLayers=8,:TestRate=6:TrainingMethod=BFGS:Sampling=0.3:SamplingEpoch=0.8:ConvergenceImprove=1e-6:ConvergenceTests=15:!UseRegulator" );
+
+
+ // --------------------------------------------------------------------------------------------------
+
+ // ---- Now you can tell the factory to train, test, and evaluate the MVAs
+
+ // Train MVAs using the set of training events
+ factory->TrainAllMethods();
+
+ // ---- Evaluate all MVAs using the set of test events
+ factory->TestAllMethods();
+
+ // ----- Evaluate and compare performance of all configured MVAs
+ factory->EvaluateAllMethods();
+
+ // --------------------------------------------------------------
+
+ // Save the output
+ outputFile->Close();
+
+ std::cout << "==> Wrote root file: " << outputFile->GetName() << std::endl;
+ std::cout << "==> TMVARegression is done!" << std::endl;
+
+ delete factory;
+
+ // Launch the GUI for the root macros
+if (!gROOT->IsBatch()) TMVARegGui( outfileName );
+}
+/*
+
+g++ -o Train_MLP `root-config --cflags` Train_MLP.cxx `root-config --glibs` -lTMVA -I$ROOTSYS/tmva/test/
+*/
\ No newline at end of file
diff --git a/Utils/trunk/EQM/MLP_Kinf/include/GenerateRootFile.hxx b/Utils/trunk/EQM/MLP_Kinf/include/GenerateRootFile.hxx
new file mode 100755
index 000000000..ef7b9ff21
--- /dev/null
+++ b/Utils/trunk/EQM/MLP_Kinf/include/GenerateRootFile.hxx
@@ -0,0 +1,51 @@
+#include "ZAI.hxx"
+#include <TGraph.h>
+#include <TGraph.h>
+
+#include <vector>
+#include <stdio.h>
+#include <stdlib.h>
+#include <iostream>
+#include <fstream>
+#include <sstream>
+#include <string>
+#include <vector>
+#include <sstream>
+#include <cmath>
+#include <map>
+#include <iostream>
+#include <iomanip>
+using namespace std;
+
+string dtoa(double num)
+{
+ ostringstream os(ostringstream::out);
+ os<<setprecision(3)<<num;
+ return os.str();
+}
+
+string FilePath;
+string DataPath;
+
+map<ZAI, double> ZAImass;
+vector<string> JobName;
+vector<string> GoodJobName;
+
+vector<double> fTime; //Time vector of the depletion calculation (second)
+vector< vector<double> > fkeff;
+
+vector<ZAI> fAllNuclei; //All the nuclei present in the fuel
+
+
+vector<IsotopicVector> fActinideCompoInit; //Fresh fuel composition
+
+int fNOfTimeStep=0; //number of time step in the Evolution
+
+string fEvolutionDataFolder = "";
+
+bool fIsAllNucleiAlreadyFill=false;
+
+void InitMass();
+void CheckJob();
+void ReadAndFill(string jobname);
+void DumpInputNeuron(string filename);
\ No newline at end of file
diff --git a/Utils/trunk/EQM/MLP_Kinf/include/StringLine.hxx b/Utils/trunk/EQM/MLP_Kinf/include/StringLine.hxx
new file mode 100755
index 000000000..da8500370
--- /dev/null
+++ b/Utils/trunk/EQM/MLP_Kinf/include/StringLine.hxx
@@ -0,0 +1,275 @@
+#ifndef _STRINGLINE_
+#define _STRINGLINE_
+
+#include <string>
+#include <sstream>
+#include <iostream>
+#include <algorithm>
+#include <cctype>
+using namespace std;
+/*!
+ \file
+ \brief Header file for StingLine class.
+*/
+
+// Class extracting fields from a string / line.
+/*!
+ The aim of this class is to provide tools to extract fields ("word") from
+ a string and convert a string in Upper/Lower case.
+ All methods are static so that it is not necessary to create object to use them
+
+ example:
+ \code
+ string line="The temperature is : 300.6 K";
+ int start;
+
+ 1st method: creation of StringLine
+
+ start=0;
+ StringLine SL;
+ string the=SL.NextWord(line,start);
+ string temperature_is=SL.NextWord(line,start,':');
+ string colon=SL.NextWord(line,start);
+ double T=atof(SL.NextWord(line,start).c_str());
+ cout<<the<<endl<<temperature_is<<endl<<T<<endl;
+
+ 2nd method: "using" the static methods
+
+ start=0;
+ the=StringLine::NextWord(line,start);
+ temperature_is=StringLine::NextWord(line,start,':');
+ colon=StringLine::NextWord(line,start);
+ T=atof(StringLine::NextWord(line,start).c_str());
+ cout<<the<<endl<<temperature_is<<endl<<T<<endl;
+ \endcode
+ @author PTO
+ @version 0.1
+*/
+
+class StringLine
+{
+ public:
+ // Find the next word in a line.
+ /*!
+ Find Next word in a line starting from position "start" in the line. If an alternative
+ separator is given, the word length is defined by the first position of sep or alt_sep found.
+ The first value of start is in general 0 (i.e. the beginning of the Line)
+ \param Line : a line containing words
+ \param start : from where to start to find the begining of a word
+ \param sep : the separator between 2 words (default=space)
+ \param alt_sep : the alternative separator between 2 words (default='')
+ */
+ static string NextWord(string Line,int &start,char sep=' ', char alt_sep='\0');
+ // Find the previous word in a line.
+ /*!
+ Find Previous word in a line starting from position "start" in the line. If an alternative
+ separator is given, the word length is defined by the first position of sep or alt_sep found.
+ The first value of start is in general the end of the Line.
+ \param Line : a line containing words
+ \param start : from where to start to find the begining of a word
+ \param sep : the separator between 2 words (default=space)
+ \param alt_sep : the alternative separator between 2 words (default='')
+ */
+ static string PreviousWord(string Line,int &start,char sep=' ', char alt_sep='\0');
+ static void ToLower(string &Line); // convert a string to Lower case
+ static void ToUpper(string &Line); // convert a string to Upper case
+
+ // Find \p search in \p Line from the begining.
+ /*!
+ returns the position, starting from the begenning of the first occurence
+ of \p search in \p Line if it is found, else returns -1
+ \param search : a string to find
+ \param Line : where to search
+ */
+ static int Find(const char *search,string Line);
+ // Find \p search in \p Line from the end.
+ /*!
+ returns the position, starting from the end of the first occurence
+ of \p search in \p Line if it is found, else returns -1
+ \param search : a string to find
+ \param Line : where to search
+ */
+ static int rFind(const char *search,string Line);
+ // convert a input type (\p in_T) to another (\p out_T).
+ /*!
+ Example:
+ \code
+ string s="32.12";
+ double t=StringLine::convert<double>(s);
+ string temperature=StringLine::convert<string>(300.);
+ \endcode
+ \param t : the input value
+ */
+ template <class out_T, class in_T> static out_T convert(const in_T & t);
+ // Find the start of a word in a line.
+ /*!
+ \param Line : a line containing words
+ \param CurrentPosition : from where to start to find the begining of a word
+ \param sep : the separator between 2 words (default=space)
+ \param alt_sep : the alternative separator between 2 words (default='')
+ */
+ static int GetStartWord(string Line,int CurrentPosition,char sep=' ', char alt_sep='\0');
+ // Find the end of a word in a line.
+ /*!
+ \param Line : a line containing words
+ \param CurrentPosition : from where to start to find the end of a word
+ \param sep : the separator between 2 words (default=space)
+ \param alt_sep : the alternative separator between 2 words (default='')
+ */
+ static int GetEndWord(string Line,int CurrentPosition,char sep=' ', char alt_sep='\0');
+ // Replace a sub-string by an other in a string.
+ /*!
+ \param InLine : the string which contains the sub-string to replace
+ \param ToReplace : the sub-string to replace
+ \param By : the sub-string ToReplace is replaced by the sub-string By in Inline
+ */
+ string ReplaceAll(string InLine, string ToReplace, string By);
+};
+
+
+//_________________________________________________________________________________
+inline string StringLine::NextWord(string Line,int &start,char sep, char alt_sep)
+{
+ string Word="";
+ if(start>=int(Line.size()))
+ {
+ return Word;
+ }
+ start=GetStartWord(Line,start,sep,alt_sep);
+ int wordlength=GetEndWord(Line,start,sep,alt_sep)-start;
+
+ Word=Line.substr(start,wordlength);
+
+ start+=wordlength;
+ return Word;
+}
+//_________________________________________________________________________________
+inline string StringLine::PreviousWord(string Line,int &start,char sep, char alt_sep)
+{
+ string Word="";
+ if(start<=0)
+ {
+ return Word;
+ }
+ int pos=Line.rfind(sep,start);
+ int alt_pos=-1;
+ int real_pos=pos;
+ char real_sep=sep;
+ if(alt_sep!='\0')
+ {
+ alt_pos=Line.rfind(alt_sep,start);
+ real_pos=max(pos,alt_pos);
+ if(real_pos!=pos)
+ real_sep=alt_sep;
+ }
+ int wordlength=start-Line.rfind(real_sep,real_pos);
+ if(real_pos<=0)
+ {
+ Word=Line.substr(0,start+1);
+ start=0;
+ return Word;
+ }
+ Word=Line.substr(real_pos+1,wordlength);
+
+ start-=wordlength+1;
+ return Word;
+}
+
+//_________________________________________________________________________________
+inline void StringLine::ToLower(string &Line)
+{
+ transform (Line.begin(), Line.end(), // source
+ Line.begin(), // destination
+ (int(*)(int))tolower); // operation
+}
+
+//_________________________________________________________________________________
+inline void StringLine::ToUpper(string &Line)
+{
+ transform (Line.begin(), Line.end(), // source
+ Line.begin(), // destination
+ (int(*)(int))toupper); // operation
+}
+
+//_________________________________________________________________________________
+inline int StringLine::GetStartWord(string Line,int CurrentPosition,char sep, char alt_sep)
+{
+ int pos=Line.find(sep,CurrentPosition);
+ int alt_pos=-1;
+ if(alt_sep!='\0')
+ alt_pos=Line.find(alt_sep,CurrentPosition);
+ int real_pos=pos;
+ char real_sep=sep;
+ if(alt_pos>=0)
+ {
+ real_pos=min(pos,alt_pos);
+ if(pos==int(string::npos))real_pos=alt_pos;
+ if(real_pos!=pos)
+ real_sep=alt_sep;
+ }
+ if(real_pos==int(string::npos)) return CurrentPosition;
+ while(CurrentPosition<int(Line.size()) && Line[CurrentPosition]==real_sep)
+ CurrentPosition++;
+ return CurrentPosition;
+}
+
+//_________________________________________________________________________________
+inline int StringLine::GetEndWord(string Line,int CurrentPosition,char sep, char alt_sep)
+{
+ int pos=Line.find(sep,CurrentPosition);
+ int alt_pos=-1;
+ if(alt_sep!='\0')
+ alt_pos=Line.find(alt_sep,CurrentPosition);
+ int real_pos=pos;
+ if(alt_pos>=0)
+ {
+ real_pos=min(pos,alt_pos);
+ if(pos==int(string::npos))real_pos=alt_pos;
+ }
+ if(real_pos==int(string::npos))
+ return Line.size();
+ return real_pos;
+}
+
+//_________________________________________________________________________________
+inline int StringLine::Find(const char *search,string Line)
+{
+ size_t Pos=Line.find(search);
+ if(Pos != string::npos ) return Pos;
+ return -1;
+}
+
+//_________________________________________________________________________________
+inline int StringLine::rFind(const char *search,string Line)
+{
+ size_t Pos=Line.rfind(search);
+ if(Pos != string::npos) return Pos;
+ return -1;
+}
+
+//_________________________________________________________________________________
+template <class out_T, class in_T>
+inline out_T StringLine::convert(const in_T & t)
+{
+ stringstream stream;
+ stream << t; // insert value to stream
+ out_T result; // store conversion's result here
+ stream >> result; // write value to result
+ return result;
+}
+
+//_________________________________________________________________________________
+inline string StringLine::ReplaceAll(string InLine, string ToReplace, string By)
+{
+ int start=0;
+ int pos=InLine.find(ToReplace,start);
+ while(pos!=int(string::npos))
+ {
+ InLine.replace(pos,ToReplace.size(),By);
+ start=0;
+ pos=InLine.find(ToReplace,start);
+ }
+ return InLine;
+
+}
+#endif
diff --git a/Utils/trunk/EQM/MLP_Kinf/include/ZAI.hxx b/Utils/trunk/EQM/MLP_Kinf/include/ZAI.hxx
new file mode 100755
index 000000000..1ee9136f9
--- /dev/null
+++ b/Utils/trunk/EQM/MLP_Kinf/include/ZAI.hxx
@@ -0,0 +1,42 @@
+#ifndef MN_ZAI_H_
+#define MN_ZAI_H_
+
+#include <map>
+using namespace std;
+
+class ZAI {
+ public :
+ int Z;
+ int A;
+ int I;
+
+
+ ZAI(int z, int a, int i=0) { Z=z;A=a; I=i;}
+
+ bool operator <(const ZAI& zai) const {return (Z != zai.Z)?
+ (Z < zai.Z) : ( (A != zai.A)?
+ (A < zai.A) : (I < zai.I) );}
+};
+
+class IsotopicVector {
+ public :
+
+ IsotopicVector() {} ///< Normal Constructor.
+ ~IsotopicVector() {}; ///< Normal Destructor.
+ double GetZAIIsotopicQuantity(const ZAI& zai) const
+ {
+ map<ZAI ,double> IsotopicQuantity = IVquantity;
+
+ map<ZAI ,double>::iterator it;
+ it = IsotopicQuantity.find(zai);
+
+ if ( it != IsotopicQuantity.end() )
+ return it->second;
+ else
+ return 0;
+ }
+ map<ZAI,double> IVquantity;
+ double MASS;
+};
+
+#endif //MN_ZAI_H_
--
GitLab