From b6e9c724cd4fd33fb48d39f17c156a9f2bb0e27d Mon Sep 17 00:00:00 2001
From: Baptiste LENIAU <baptiste.leniau@subatech.in2p3.fr>
Date: Fri, 22 Apr 2016 14:49:26 +0000
Subject: [PATCH] copy Utils from 4.1
git-svn-id: svn+ssh://svn.in2p3.fr/class@856 0e7d625b-0364-4367-a6be-d5be4a48d228
---
.../EQM/FBR_MLP_Keff/GenerateRootFile.cxx | 347 ++++++++
.../BaM/Utils/EQM/FBR_MLP_Keff/Train_MLP.cxx | 177 ++++
.../FBR_MLP_Keff/include/GenerateRootFile.hxx | 51 ++
.../EQM/FBR_MLP_Keff/include/StringLine.hxx | 275 ++++++
.../Utils/EQM/FBR_MLP_Keff/include/ZAI.hxx | 42 +
.../Utils/EQM/MLP_Kinf/GenerateRootFile.cxx | 343 ++++++++
branches/BaM/Utils/EQM/MLP_Kinf/Train_MLP.cxx | 178 ++++
.../EQM/MLP_Kinf/include/GenerateRootFile.hxx | 51 ++
.../Utils/EQM/MLP_Kinf/include/StringLine.hxx | 275 ++++++
.../BaM/Utils/EQM/MLP_Kinf/include/ZAI.hxx | 42 +
.../BaM/Utils/EQM/PWR_MOX_MLP/Train_MLP.cxx | 163 ++++
.../BaM/Utils/MURE2CLASS/BinaryFormat2.hxx | 141 ++++
branches/BaM/Utils/MURE2CLASS/MURE2CLASS.cxx | 758 +++++++++++++++++
branches/BaM/Utils/MURE2CLASS/StringLine.hxx | 275 ++++++
branches/BaM/Utils/MURE2CLASS/ZAI.hxx | 50 ++
.../BaM/Utils/XSM/MLP/BuildInput/Gene.cxx | 497 +++++++++++
.../BaM/Utils/XSM/MLP/BuildInput/Gene.hxx | 53 ++
.../Utils/XSM/MLP/BuildInput/StringLine.hxx | 275 ++++++
branches/BaM/Utils/XSM/MLP/BuildInput/ZAI.hxx | 42 +
.../MLP/Train/EvaluateTrainingCommands.dat | 5 +
.../BaM/Utils/XSM/MLP/Train/LaunchTraining.sh | 21 +
branches/BaM/Utils/XSM/MLP/Train/Train_XS.cxx | 185 +++++
branches/BaM/Utils/XSM/MLP/Train/deviations.C | 209 +++++
branches/BaM/Utils/XSM/MLP/Train/tmvaglob.C | 785 ++++++++++++++++++
branches/BaM/Utils/rootlogon.C | 5 +
25 files changed, 5245 insertions(+)
create mode 100755 branches/BaM/Utils/EQM/FBR_MLP_Keff/GenerateRootFile.cxx
create mode 100755 branches/BaM/Utils/EQM/FBR_MLP_Keff/Train_MLP.cxx
create mode 100755 branches/BaM/Utils/EQM/FBR_MLP_Keff/include/GenerateRootFile.hxx
create mode 100755 branches/BaM/Utils/EQM/FBR_MLP_Keff/include/StringLine.hxx
create mode 100755 branches/BaM/Utils/EQM/FBR_MLP_Keff/include/ZAI.hxx
create mode 100755 branches/BaM/Utils/EQM/MLP_Kinf/GenerateRootFile.cxx
create mode 100755 branches/BaM/Utils/EQM/MLP_Kinf/Train_MLP.cxx
create mode 100755 branches/BaM/Utils/EQM/MLP_Kinf/include/GenerateRootFile.hxx
create mode 100755 branches/BaM/Utils/EQM/MLP_Kinf/include/StringLine.hxx
create mode 100755 branches/BaM/Utils/EQM/MLP_Kinf/include/ZAI.hxx
create mode 100755 branches/BaM/Utils/EQM/PWR_MOX_MLP/Train_MLP.cxx
create mode 100644 branches/BaM/Utils/MURE2CLASS/BinaryFormat2.hxx
create mode 100755 branches/BaM/Utils/MURE2CLASS/MURE2CLASS.cxx
create mode 100644 branches/BaM/Utils/MURE2CLASS/StringLine.hxx
create mode 100644 branches/BaM/Utils/MURE2CLASS/ZAI.hxx
create mode 100755 branches/BaM/Utils/XSM/MLP/BuildInput/Gene.cxx
create mode 100755 branches/BaM/Utils/XSM/MLP/BuildInput/Gene.hxx
create mode 100755 branches/BaM/Utils/XSM/MLP/BuildInput/StringLine.hxx
create mode 100755 branches/BaM/Utils/XSM/MLP/BuildInput/ZAI.hxx
create mode 100644 branches/BaM/Utils/XSM/MLP/Train/EvaluateTrainingCommands.dat
create mode 100755 branches/BaM/Utils/XSM/MLP/Train/LaunchTraining.sh
create mode 100755 branches/BaM/Utils/XSM/MLP/Train/Train_XS.cxx
create mode 100755 branches/BaM/Utils/XSM/MLP/Train/deviations.C
create mode 100755 branches/BaM/Utils/XSM/MLP/Train/tmvaglob.C
create mode 100755 branches/BaM/Utils/rootlogon.C
diff --git a/branches/BaM/Utils/EQM/FBR_MLP_Keff/GenerateRootFile.cxx b/branches/BaM/Utils/EQM/FBR_MLP_Keff/GenerateRootFile.cxx
new file mode 100755
index 000000000..373bd7037
--- /dev/null
+++ b/branches/BaM/Utils/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/branches/BaM/Utils/EQM/FBR_MLP_Keff/Train_MLP.cxx b/branches/BaM/Utils/EQM/FBR_MLP_Keff/Train_MLP.cxx
new file mode 100755
index 000000000..7ca138d39
--- /dev/null
+++ b/branches/BaM/Utils/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/branches/BaM/Utils/EQM/FBR_MLP_Keff/include/GenerateRootFile.hxx b/branches/BaM/Utils/EQM/FBR_MLP_Keff/include/GenerateRootFile.hxx
new file mode 100755
index 000000000..ef7b9ff21
--- /dev/null
+++ b/branches/BaM/Utils/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/branches/BaM/Utils/EQM/FBR_MLP_Keff/include/StringLine.hxx b/branches/BaM/Utils/EQM/FBR_MLP_Keff/include/StringLine.hxx
new file mode 100755
index 000000000..da8500370
--- /dev/null
+++ b/branches/BaM/Utils/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/branches/BaM/Utils/EQM/FBR_MLP_Keff/include/ZAI.hxx b/branches/BaM/Utils/EQM/FBR_MLP_Keff/include/ZAI.hxx
new file mode 100755
index 000000000..1ee9136f9
--- /dev/null
+++ b/branches/BaM/Utils/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/branches/BaM/Utils/EQM/MLP_Kinf/GenerateRootFile.cxx b/branches/BaM/Utils/EQM/MLP_Kinf/GenerateRootFile.cxx
new file mode 100755
index 000000000..8bd2d961c
--- /dev/null
+++ b/branches/BaM/Utils/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/branches/BaM/Utils/EQM/MLP_Kinf/Train_MLP.cxx b/branches/BaM/Utils/EQM/MLP_Kinf/Train_MLP.cxx
new file mode 100755
index 000000000..8f8aa47e8
--- /dev/null
+++ b/branches/BaM/Utils/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/branches/BaM/Utils/EQM/MLP_Kinf/include/GenerateRootFile.hxx b/branches/BaM/Utils/EQM/MLP_Kinf/include/GenerateRootFile.hxx
new file mode 100755
index 000000000..ef7b9ff21
--- /dev/null
+++ b/branches/BaM/Utils/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/branches/BaM/Utils/EQM/MLP_Kinf/include/StringLine.hxx b/branches/BaM/Utils/EQM/MLP_Kinf/include/StringLine.hxx
new file mode 100755
index 000000000..da8500370
--- /dev/null
+++ b/branches/BaM/Utils/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/branches/BaM/Utils/EQM/MLP_Kinf/include/ZAI.hxx b/branches/BaM/Utils/EQM/MLP_Kinf/include/ZAI.hxx
new file mode 100755
index 000000000..1ee9136f9
--- /dev/null
+++ b/branches/BaM/Utils/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_
diff --git a/branches/BaM/Utils/EQM/PWR_MOX_MLP/Train_MLP.cxx b/branches/BaM/Utils/EQM/PWR_MOX_MLP/Train_MLP.cxx
new file mode 100755
index 000000000..e701f0dc0
--- /dev/null
+++ b/branches/BaM/Utils/EQM/PWR_MOX_MLP/Train_MLP.cxx
@@ -0,0 +1,163 @@
+//
+// This program train and test a MLP using TMVA from a training data in a form of a TTRee
+//
+//@author Root_tmva_team modified by BaL
+#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!=1)
+ {
+ std::cout << "TrainMLP Usage :" << std::endl;
+ std::cout << "\tTrainMLP YourTrainingData.root" << std::endl;
+ }
+ //---------------------------------------------------------------
+ // This loads the library
+ TMVA::Tools::Instance();
+ // ---------------------------------------------------------------
+
+ std::cout << std::endl;
+ std::cout << "==> Start TMVARegression" << std::endl;
+
+
+ // Create a new root output file
+ TString outfileName( "TMVA_MOX_Equivalence.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" )]
+ factory->AddVariable( "BU" , "Burn Up" , "GWj_tMLi" , 'F' );
+ factory->AddVariable( "U5_enrichment","U5 enrichment", "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("teneur");
+
+ // 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_Equivalence"<<0;
+ // Neural network (MLP)
+ factory->BookMethod( TMVA::Types::kMLP, Name.str().c_str(), "!H:!V:VarTransform=Norm:NeuronType=tanh:NCycles=20000:HiddenLayers=N+10: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/branches/BaM/Utils/MURE2CLASS/BinaryFormat2.hxx b/branches/BaM/Utils/MURE2CLASS/BinaryFormat2.hxx
new file mode 100644
index 000000000..570a0010d
--- /dev/null
+++ b/branches/BaM/Utils/MURE2CLASS/BinaryFormat2.hxx
@@ -0,0 +1,141 @@
+#ifndef _BINARYFORMAT2_
+#define _BINARYFORMAT2_
+
+#include <fstream>
+using namespace std;
+
+/*!
+ \file
+ \brief Header file for binary output files management.
+
+ This file is to be included in any units which manipulate (read/void write) MURE
+ binary output files. These structures allow easy way to handle
+ binary output of MURE evolution data (BDATA_* files).
+
+ @author JW
+ @author PTO
+ @version 1.0
+*/
+
+//! Header of MURE output binary file
+struct FileHeader
+{
+ FileHeader() {Version=1;} //!< Constructor
+ short Version; //!< Binary gile version
+ float Time; //!< Time of MCNP step at which this file is printed
+ float K; //!< keff
+ float Kerr; //!< keff error
+ int NCells; //!< Number of evolving cells
+ void write(ofstream &out); //!< Write the header into a stream
+ void read(ifstream &in); //!< Read the header from a stream
+};
+
+//! Header of an evolving cell in a binary file
+struct CellHeader
+{
+ int CellNumber; //!< MCNP number of this cell
+ float Volume; //!< Volume of the cell
+ float Flux; //!< Flux in this cell
+ float FluxErr; //!< Flux error
+ int NNucleusRecords; //!< Number of nuclei records in this cell
+ void write(ofstream &out); //!< Write the header into a stream
+ void read(ifstream &in); //!< Read the header from a stream
+};
+
+//! Record of a nucleus in a binary file
+struct NucleusRecord
+{
+ short Z; //!< Proton number of the nucleus
+ short A; //!< Nucleon number of the nucleus
+ short I; //!< Isomeric state of the nucleus
+ float Mass; //!< Atomic mass
+ float Proportion; //!< Number of this nuclei in cell
+ short NReactionRecords; //!< Number of reaction records of this nucleus
+ void write(ofstream &out); //!< Write the record into a stream
+ void read(ifstream &in); //!< Read the record from a stream
+};
+
+//! Record of a reaction in a binary file
+struct ReactionRecord
+{
+ short Code; //!< Reaction code
+ float Sigma; //!< Reaction cross-section
+ float SigmaErr; //!< Cross-section error
+ void write(ofstream &out); //!< Write the record into a stream
+ void read(ifstream &in); //!< Read the record from a stream
+};
+//--------------------------------------------------------
+// Implementation of methods
+//-------------------------------------------------------
+void FileHeader::write(ofstream &out)
+{
+ out.write((char*)&Version,sizeof(Version));
+ out.write((char*)&Time,sizeof(Time));
+ out.write((char*)&K,sizeof(K));
+ out.write((char*)&Kerr,sizeof(Kerr));
+ out.write((char*)&NCells,sizeof(NCells));
+}
+
+void FileHeader::read(ifstream &in)
+{
+ in.read((char*)&Version,sizeof(Version));
+ in.read((char*)&Time,sizeof(Time));
+ in.read((char*)&K,sizeof(K));
+ in.read((char*)&Kerr,sizeof(Kerr));
+ in.read((char*)&NCells,sizeof(NCells));
+}
+
+void CellHeader::write(ofstream &out)
+{
+ out.write((char*)&CellNumber,sizeof(CellNumber));
+ out.write((char*)&Volume,sizeof(Volume));
+ out.write((char*)&Flux,sizeof(Flux));
+ out.write((char*)&FluxErr,sizeof(FluxErr));
+ out.write((char*)&NNucleusRecords,sizeof(NNucleusRecords));
+}
+
+void CellHeader::read(ifstream &in)
+{
+ in.read((char*)&CellNumber,sizeof(CellNumber));
+ in.read((char*)&Volume,sizeof(Volume));
+ in.read((char*)&Flux,sizeof(Flux));
+ in.read((char*)&FluxErr,sizeof(FluxErr));
+ in.read((char*)&NNucleusRecords,sizeof(NNucleusRecords));
+}
+
+void NucleusRecord::write(ofstream &out)
+{
+ out.write((char*)&Z,sizeof(Z));
+ out.write((char*)&A,sizeof(A));
+ out.write((char*)&I,sizeof(I));
+ out.write((char*)&Mass,sizeof(Mass));
+ out.write((char*)&Proportion,sizeof(Proportion));
+ out.write((char*)&NReactionRecords,sizeof(NReactionRecords));
+}
+
+void NucleusRecord::read(ifstream &in)
+{
+ in.read((char*)&Z,sizeof(Z));
+ in.read((char*)&A,sizeof(A));
+ in.read((char*)&I,sizeof(I));
+ in.read((char*)&Mass,sizeof(Mass));
+ in.read((char*)&Proportion,sizeof(Proportion));
+ in.read((char*)&NReactionRecords,sizeof(NReactionRecords));
+}
+
+void ReactionRecord::write(ofstream &out)
+{
+ out.write((char*)&Code,sizeof(Code));
+ out.write((char*)&Sigma,sizeof(Sigma));
+ out.write((char*)&SigmaErr,sizeof(SigmaErr));
+}
+
+void ReactionRecord::read(ifstream &in)
+{
+ in.read((char*)&Code,sizeof(Code));
+ in.read((char*)&Sigma,sizeof(Sigma));
+ in.read((char*)&SigmaErr,sizeof(SigmaErr));
+}
+
+
+#endif
diff --git a/branches/BaM/Utils/MURE2CLASS/MURE2CLASS.cxx b/branches/BaM/Utils/MURE2CLASS/MURE2CLASS.cxx
new file mode 100755
index 000000000..ff93a7615
--- /dev/null
+++ b/branches/BaM/Utils/MURE2CLASS/MURE2CLASS.cxx
@@ -0,0 +1,758 @@
+// DESCRIPTION
+// This program convert MURE output tp EvolutionData forma
+//
+//@author BaM
+#include "BinaryFormat2.hxx"
+#include "StringLine.hxx"
+#include "ZAI.hxx"
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <iostream>
+#include <fstream>
+#include <sstream>
+#include <string>
+#include <vector>
+#include <algorithm>
+
+using std::string;
+using namespace std;
+
+//----------------------------------------------------------------------------------------------------------------
+//----------------------------------------------------------------------------------------------------------------
+//------------------------------------- METHODS ------------------------------------------------------------------
+//----------------------------------------------------------------------------------------------------------------
+//----------------------------------------------------------------------------------------------------------------
+
+string itoa(int num)
+{
+ ostringstream os(ostringstream::out);
+ os << "" << num;
+ return os.str();
+}
+
+int ReadCommentVersion(ifstream& in, string filename)
+{
+ if(!in.good())
+ {
+ cerr << "ERROR in MureRead:: Cannot find ASCII file " << filename << " to Read" << endl;
+ return -1;
+ }
+ //
+ //Go throught the comments at the beginning of the file and Read the version number of Data Format
+ //
+ string TheComment;
+ int MureDataVersion;
+ bool bComment=true;
+ while (bComment)
+ {
+ in>>TheComment;
+ if (TheComment=="%")
+ {
+ getline(in,TheComment); // the "in>>" command keep the cursor at the end of line
+ }
+ else if (TheComment=="V")
+ {
+ in>>MureDataVersion;
+ bComment=false;
+ }
+ else
+ {
+ cout << "No comments or format version at the beginning of data files" << endl;
+ bComment=false;
+ MureDataVersion=0;
+ }
+ }
+ return MureDataVersion;
+}
+
+//----------------------------------------------------------------------------------------------------------------
+//----------------------------------------------------------------------------------------------------------------
+//------------------------------------- MAIN ---------------------------------------------------------------------
+//----------------------------------------------------------------------------------------------------------------
+//----------------------------------------------------------------------------------------------------------------
+
+int main(int argc, char** argv)
+{
+
+ // =========================================================================================
+ // VARIABLES
+ // =========================================================================================
+ string sPWD = getenv("PWD");
+
+ string s_tmp = ""; int i_tmp = 0; double tmp = 0;
+ string DataASCIIFile = "DATA_"; int NDataASCIIFile = 0;
+ string DataBinFile = "BDATA_"; int NDataBinFile = 0;
+ string DataFile = ""; int NDataFile = 0;
+ string LastDataFile = "";
+
+ string OutDataFile = "";
+ string OutDataFileInfo = "";
+ string OutLOG = "";
+
+ bool IsThereASCCIFile = false;
+ bool IsThereBinFile = false;
+
+ vector <int> vCellNumber;
+ vector <string> vCellComment;
+
+ int NumberOfCells = 1;
+
+ vector < double > vTime;
+ vector < double > vKeff;
+ vector < double > vFlux;
+ vector < vector < double > > vFlux1;
+
+ vector <ZAI > zai1;
+ vector <vector <ZAI > > zai2;
+ vector <vector <vector<ZAI > > > zai3;
+
+ vector <ZAI > zai_FS_1;
+ vector <vector <ZAI > > zai_FS_2;
+ vector <vector <vector<ZAI > > > zai_FS_3;
+
+
+ vector <ZAI > zai_SC_1;
+ vector <vector <ZAI > > zai_SC_2;
+ vector <vector <vector<ZAI > > > zai_SC_3;
+
+
+ vector <ZAI > zai_SN_1;
+ vector <vector <ZAI > > zai_SN_2;
+ vector <vector <vector<ZAI > > > zai_SN_3;
+
+
+ vector <bool > HasToBePrint;
+
+ // =========================================================================================
+ // WELCOME !
+ // =========================================================================================
+
+ // MURE output path
+ if(argc != 9)
+ {
+ cout << "Argument problem for MURE2CLASS... EXIT!" << endl;
+ cout << "Arg should be :" << endl << "\t1 Path," << endl << "\t2 OutName," << endl << "\t3 ReactorType," << endl << "\t4 FuelType," << endl;
+ cout << "\t5 Power," << endl << endl << "\t6 NormalizationFactor(Before Normalization)," << endl << "\t7 WantedCell (by default should be 0)" << endl << "\t8 Step to Skip," << endl;
+ exit(1);}
+ string DBPath = argv[1];
+ string OutName = argv[2];
+ string ReactorType = argv[3];
+ string FuelType = argv[4];
+
+ string Power = argv[5];
+ double NormalizationFactor = atof(argv[6]);
+ int WantedCell = atoi(argv[7]);
+ int StepToSkip = atoi(argv[8]);
+ // Name of the output file
+ OutDataFile = "./" + OutName + ".dat";
+ OutDataFileInfo = "./" + OutName + ".info";
+ OutLOG = "./" + OutName + ".log";
+
+ ofstream OutputLog(OutLOG.c_str());
+
+ // Check how many DATAxxx and BDATAxxx files there are
+ string Command = "find " + DBPath + " -name \"" + DataASCIIFile + "*\" > ASCII.tmp";
+ system(Command.c_str());
+ Command = "find " + DBPath + " -name \"" + DataBinFile + "*\" > BIN.tmp";
+ system(Command.c_str());
+
+ ifstream ASCIITMP("ASCII.tmp");
+ if (ASCIITMP.is_open())
+ {
+ while (!ASCIITMP.eof())
+ {
+ getline(ASCIITMP,s_tmp);
+ i_tmp++;
+ }
+ } ASCIITMP.close(); i_tmp--; NDataASCIIFile = i_tmp;
+
+ i_tmp=0;
+ ifstream BINTMP("BIN.tmp");
+ if (BINTMP.is_open())
+ {
+ while (!BINTMP.eof())
+ {
+ getline(BINTMP,s_tmp);
+ i_tmp++;
+ }
+ } BINTMP.close(); i_tmp--; NDataBinFile = i_tmp;
+
+ // Remove temporary files...
+ Command = "\\rm -f ASCII.tmp BIN.tmp";
+ system(Command.c_str());
+
+ // Say if there is Binary or ASCII file
+ if (NDataASCIIFile>0) {IsThereASCCIFile=true; DataFile = DataASCIIFile; NDataFile = NDataASCIIFile;}
+ else if (NDataBinFile>0) {IsThereBinFile =true; DataFile = DataBinFile; NDataFile = NDataBinFile;}
+ else {cout << endl << "There is no MURE DATA file in the given path... EXIT!" << endl << endl; exit(1);}
+
+ // =========================================================================================
+ // READ (B)DATA files...
+ // =========================================================================================
+
+ OutputLog << endl;
+ OutputLog << "===================================================" << endl;
+ OutputLog << "---------------------------------------------------" << endl;
+ if (IsThereBinFile) OutputLog << endl << "Binary file detected... " << endl << endl;
+ else cout << endl << "ASCII file detected... " << endl << endl;
+ // /*sleep(1)*/;
+
+ for (int t=0; t<NDataFile; t++)
+ {
+ string SFX = "";
+ if (t<=9) SFX = "00" + itoa(t);
+ else if (t>=10 && t<=99) SFX = "0" + itoa(t);
+ else if (t>=100 && t<=999) SFX = "" + itoa(t);
+ else {cout << endl << "there is more than 1000 DATA files, not yet implemented... EXIT!" << endl << endl; exit(1);}
+ string s_File = DBPath + "/" + DataFile + SFX;
+ if (t%10==0) {string ff = DataFile + SFX; cout << "Reading file " << ff << endl; /*sleep(1)*/;}
+ LastDataFile = DataFile + SFX;
+
+ // READ BINARY FILE
+ if (IsThereBinFile)
+ {
+ ifstream in(s_File.c_str(),ios::binary);
+ if(!in.good()){cout << "Cannot find Binary file " << s_File << endl; exit(0);}
+ //
+ //Read The version number of writing Mure evolving Data Format
+ //
+ char TextV;
+ int MureDataVersion=0;
+ in.read((char*)&TextV, sizeof(char));
+
+ if(TextV=='V') in.read((char*)&MureDataVersion, sizeof(int));
+ else {in.close(); in.clear(); in.open(s_File.c_str(),ios::binary);}
+
+ string TheComment;
+
+ // Read the File Header
+ FileHeader FH;
+ FH.read(in);
+ NumberOfCells = FH.NCells;
+
+ // --------------------------- TIME ----------------------
+ vTime.push_back(FH.Time);
+ vKeff.push_back(FH.K);
+
+ CellHeader CH;
+ NucleusRecord NR;
+ bool STOP2READ=false;
+
+ for(int c=0 ; c<NumberOfCells; c++)
+ {
+ if(STOP2READ)break;
+ CH.read(in);
+ vFlux.push_back(CH.Flux);
+
+ // read the Cell Comment
+ int StringSize; in.read((char*)&StringSize, sizeof(int)); TheComment.resize(StringSize);
+ char tmp[StringSize+1]; in.read(tmp, (StringSize+1)*sizeof(char)); tmp[StringSize]='\0';
+ TheComment=tmp;
+
+ bool SkipCell=false;
+ if(CH.CellNumber<0 && TheComment!="WasteStorage") SkipCell=true;
+ if(TheComment=="WasteStorage") STOP2READ=true;
+
+ if (t==0)
+ {
+ vCellNumber.push_back(CH.CellNumber);
+ vCellComment.push_back(TheComment);
+ }
+
+ //read the spatial variables
+ int NSpatialVariables; in.read((char*)&NSpatialVariables, sizeof(int));
+
+ vector<double> SpatialVariables(NSpatialVariables);
+ vector<string> SpatialVariableNames;
+
+ for(int j=0; j<NSpatialVariables; j++) in.read((char*)&SpatialVariables[j],sizeof(double));
+
+ //read the spatial variable names
+ for(int j=0; j<NSpatialVariables; j++)
+ {
+ int StringSize; in.read((char*)&StringSize,sizeof(int));
+ s_tmp=""; s_tmp.resize(StringSize); char tmp[StringSize+1];
+ in.read(tmp, (StringSize+1)*sizeof(char)); tmp[StringSize]='\0';
+ s_tmp=tmp; SpatialVariableNames.push_back(s_tmp);
+ }
+
+ for(int i=0; i<CH.NNucleusRecords; i++)
+ {
+ // Read the Nucleus Record
+ NR.read(in);
+
+ if(!SkipCell)
+ {
+ ZAI zai(NR.Z, NR.A, NR.I, NR.Proportion);
+ zai1.push_back(zai);
+ //cout<<zai.Z()<<" "<<zai.A()<<" "<<zai.I()<<" "<<zai.Prop()<<endl;
+ //if (t==0 && NR.Proportion>=2e+25) cout << NR.Z << " " << NR.A << " " << NR.I << " " << NR.Proportion << endl;
+ }
+
+ ReactionRecord RR;
+ for (int k=0; k<NR.NReactionRecords; k++)
+ {
+ RR.read(in);
+ if(RR.Sigma==0.0) RR.Sigma=1e-10;
+
+ if(!SkipCell)
+ {
+ ZAI zai(NR.Z, NR.A, NR.I, RR.Sigma);
+ if (RR.Code==18) zai_FS_1.push_back(zai);
+ if (RR.Code==102) zai_SC_1.push_back(zai);
+ if (RR.Code==16) zai_SN_1.push_back(zai);
+
+ }
+
+
+
+
+
+ }
+ }
+ if(!SkipCell)
+ {
+ zai2.push_back(zai1);
+ zai1.clear();
+ zai_FS_2.push_back(zai_FS_1);
+ zai_FS_1.clear();
+ zai_SC_2.push_back(zai_SC_1);
+ zai_SC_1.clear();
+ zai_SN_2.push_back(zai_SN_1);
+ zai_SN_1.clear();
+ }
+ }
+ in.close();
+ zai3.push_back(zai2); zai2.clear();
+ zai_FS_3.push_back(zai_FS_2); zai_FS_2.clear();
+ zai_SC_3.push_back(zai_SC_2); zai_SC_2.clear();
+ zai_SN_3.push_back(zai_SN_2); zai_SN_2.clear();
+
+ vFlux1.push_back(vFlux); vFlux.clear();
+
+ }
+ // READ ASCII FILE
+ else if (IsThereASCCIFile)
+ {
+ ifstream in(s_File.c_str());
+ int TranspCode=1;
+ if(!in.good()) {if(!in.good()){OutputLog << "Cannot find ASCII file " << s_File << endl; exit(0);}}
+ //
+ //Read The version number of writing Mure evolving Data Format
+ //
+ int MureDataVersion=ReadCommentVersion(in, s_File);
+ if(MureDataVersion==0) {in.close(); in.clear(); in.open(s_File.c_str());}
+ else if(MureDataVersion<0) {OutputLog << "Old version of MURE... Check!!! exit(1)"; exit(0);}
+
+ // Read out the TIME the KEFF and the KEFF ERROR
+ double Time,Keff,Keff_Err;
+ in>>Time>>Keff>>Keff_Err;
+ in>>NumberOfCells;
+
+ // --------------------------- TIME ----------------------
+ vTime.push_back(Time);
+ vKeff.push_back(Keff);
+
+ bool STOP2READ=false;
+
+ for(int c=0 ; c<NumberOfCells; c++)
+ {
+ if(STOP2READ)break;
+ int CellNumber;
+ double Volume,Flux,FluxErr;
+ in>>CellNumber>>Volume>>Flux>>FluxErr;
+
+ vFlux.push_back(Flux);
+
+ int NNucleusRecords;
+ in>>NNucleusRecords;
+
+ // read the Cell Comment and Spatial Variables
+ string TheComment;
+ getline(in,TheComment); // the "in>>" command keep the cursor at the end of line
+ getline(in,TheComment); // so two "getline" commands are requested
+
+ bool SkipCell=false;
+ if(CellNumber<0 && TheComment!="WasteStorage") SkipCell=true;
+ if(TheComment=="WasteStorage") STOP2READ=true;
+
+ if (t==0)
+ {
+ vCellNumber.push_back(CellNumber);
+ vCellComment.push_back(TheComment);
+ }
+
+ int NSpatialVariables;
+ in>>NSpatialVariables;
+
+ vector<double> SpatialVariables;
+ vector<string> SpatialVariableNames;
+
+ for(int j=0; j<NSpatialVariables; j++)
+ {
+ string TheString;
+ double val;
+ in>>TheString>>val;
+ SpatialVariableNames.push_back(TheString);
+ SpatialVariables.push_back(val);
+ }
+
+ for(int i=0; i<NNucleusRecords; i++)
+ {
+ // Read the Nucleus Record
+ int Z,A,I;
+ double Mass,Proportion;
+ in>>Z>>A>>I>>Mass>>Proportion;
+ cout << Z << endl;
+ if(!SkipCell)
+ {
+ ZAI zai(Z, A, I, Proportion);
+ zai1.push_back(zai);
+ }
+
+ int NReactionRecords;
+ in>>NReactionRecords;
+ for (int k=0; k<NReactionRecords; k++)
+ {
+ int Code;
+ double Sigma,SigmaErr;
+ in>>Code>>Sigma>>SigmaErr;
+ if(!SkipCell)
+ {
+ ZAI zai(Z, A, I, Proportion);
+ if (Code==18) zai_FS_1.push_back(zai);
+ if (Code==102) zai_SC_1.push_back(zai);
+ if (Code==16) zai_SN_1.push_back(zai);
+
+ }
+
+ if(Sigma==0.0) Sigma=1e-10;
+ }
+
+
+
+
+ }
+ if(!SkipCell)
+ {
+ zai2.push_back(zai1);
+ zai1.clear();
+ zai_FS_2.push_back(zai_FS_1);
+ zai_FS_1.clear();
+ zai_SC_2.push_back(zai_SC_1);
+ zai_SC_1.clear();
+ zai_SN_2.push_back(zai_SN_1);
+ zai_SN_1.clear();
+ }
+ }
+ in.close();
+ zai3.push_back(zai2); zai2.clear();
+ zai_FS_3.push_back(zai_FS_2); zai_FS_2.clear();
+ zai_SC_3.push_back(zai_SC_2); zai_SC_2.clear();
+ zai_SN_3.push_back(zai_SN_2); zai_SN_2.clear();
+ vFlux1.push_back(vFlux); vFlux.clear();
+
+
+ }
+ // else ERROR
+ else {OutputLog << endl << "DATA files are neither binary nor ASCII?? ... EXIT!" << endl; exit(1);}
+ }
+ OutputLog << endl << "Last file read : " << LastDataFile << endl << endl;
+ OutputLog << "All (B)DATA files have been read..." << endl << endl;
+ OutputLog << "---------------------------------------------------" << endl;
+ OutputLog << "===================================================" << endl;
+
+ // =========================================================================================
+ // Manage several Cells
+ // =========================================================================================
+
+ // Time bins
+ int Nt = zai3.size();
+ // Number of cells
+ int Nc = zai3[0].size();
+
+ // Number of nuclides
+ int Ni = zai3[0][0].size();
+
+ double CycleTime = (vTime[vTime.size()-1] - vTime[0]) / 3600. / 24 / 365.4;
+
+ // Manage the cells...
+ bool SumOfCell = false;
+
+ if (vCellNumber.size()>=2)
+ {
+ OutputLog << endl << "===================================================" << endl;
+ OutputLog << "-------------- WARNING ----------------------------" << endl;
+ OutputLog << "===================================================" << endl << endl;
+ OutputLog << "THERE IS MORE THAN ONE CELL... Cells are : " << endl << endl; /*sleep(1)*/;
+ for(int i=0; i<vCellNumber.size(); i++)
+ {
+ OutputLog << "index : " << i << " - Cell number : " << vCellNumber[i] << endl;
+ OutputLog << "Cell comment : " << vCellComment[i] << endl << endl;
+ }
+ if(WantedCell==0) SumOfCell = true;
+
+ if (!SumOfCell && WantedCell>Nc) {OutputLog << "The index you choose is not defined... EXIT!" << endl; exit(1);}
+
+ zai2.clear();
+ zai_FS_2.clear();
+ zai_SC_2.clear();
+ zai_SN_2.clear();
+ vFlux.clear();
+ for (int t=0; t<Nt; t++)
+ {
+ for (int i=0; i<Ni; i++)
+ {
+ double SUM = 0;
+ for (int c=0; c<Nc; c++) SUM += zai3[t][c][i].Prop();
+ ZAI zai(zai3[t][0][i].Z(),zai3[t][0][i].A(),zai3[t][0][i].I(),SUM);
+ zai1.push_back(zai);
+ }
+ for (int i=0; i<(int)zai_FS_3[0][0].size(); i++)
+ {
+ double SUM = 0;
+ for (int c=0; c<Nc; c++) SUM += zai_FS_3[t][c][i].Prop()/Nc;
+ ZAI zai(zai_FS_3[t][0][i].Z(),zai_FS_3[t][0][i].A(),zai_FS_3[t][0][i].I(),SUM);
+ zai_FS_1.push_back(zai);
+ }
+ for (int i=0; i<(int)zai_SC_3[0][0].size(); i++)
+ {
+ double SUM = 0;
+ for (int c=0; c<Nc; c++) SUM += zai_SC_3[t][c][i].Prop()/Nc;
+ ZAI zai(zai_SC_3[t][0][i].Z(),zai_SC_3[t][0][i].A(),zai_SC_3[t][0][i].I(),SUM);
+ zai_SC_1.push_back(zai);
+ }
+ for (int i=0; i<(int)zai_SN_3[0][0].size(); i++)
+ {
+ double SUM = 0;
+ for (int c=0; c<Nc; c++) SUM += zai_SN_3[t][c][i].Prop()/Nc;
+ ZAI zai(zai_SN_3[t][0][i].Z(),zai_SN_3[t][0][i].A(),zai_SN_3[t][0][i].I(),SUM);
+ zai_SN_1.push_back(zai);
+ }
+
+ double SUM = 0;
+ for (int j=0; j<(int)vFlux1[t].size(); j++) SUM += vFlux1[t][j];
+
+ vFlux.push_back(SUM);
+
+ zai2.push_back(zai1);
+ zai1.clear();
+ zai_FS_2.push_back(zai_FS_1);
+ zai_FS_1.clear();
+ zai_SC_2.push_back(zai_SC_1);
+ zai_SC_1.clear();
+ zai_SN_2.push_back(zai_SN_1);
+ zai_SN_1.clear();
+
+ }
+ OutputLog << endl << "---------------------------------------------------" << endl;
+ OutputLog << "-------------- Sum of cells done ------------------" << endl;
+ OutputLog << "---------------------------------------------------" << endl << endl;
+ }
+ else WantedCell=0;
+
+
+ // =========================================================================================
+ // Manage the cutoff and calculate total N and M at t=0
+ // =========================================================================================
+
+ double NTotal = 0;
+ double MTotalFissile = 0;
+ if (SumOfCell)
+ {
+ for (int i=0; i<Ni; i++)
+ {
+ double Ni = zai2[0][i].Prop();
+ NTotal += Ni;
+ if (zai2[0][i].Z()>=90) MTotalFissile += Ni*zai2[0][i].A()/6.023e23;
+ }
+ }
+ else
+ {
+ for (int i=0; i<Ni; i++)
+ {
+ double Ni = zai3[0][WantedCell][i].Prop();
+ NTotal += Ni;
+ if (zai3[0][WantedCell][i].Z()>=90) MTotalFissile += Ni*zai3[0][WantedCell][i].A()/6.023e23;
+ }
+ }
+
+ // =========================================================================================
+ // Write the DATABASE and convert Number to Mass - Manage Normalization Factor
+ // =========================================================================================
+
+ ofstream Output(OutDataFile.c_str());
+ Output.precision(16);
+ Output << "time";
+ for(int t=StepToSkip; t<vTime.size(); t++) Output << " " << vTime.at(t)-vTime.at(StepToSkip);
+ Output << endl;
+
+ Output << "keff";
+ for(int t=StepToSkip; t<vKeff.size(); t++) Output << " " << vKeff.at(t);
+ Output << endl;
+
+ Output << "flux";
+ if (SumOfCell)
+ {
+ for(int t=StepToSkip; t<vTime.size(); t++) Output << " " << vFlux.at(t);
+ Output << endl;
+ }
+ else
+ {
+ for(int t=StepToSkip; t<vTime.size(); t++) Output << " " << vFlux1[t][WantedCell];
+ Output << endl;
+ }
+ int NPrinted=0;
+ {
+ if (SumOfCell)
+ {
+ for(int i=0; i < zai2[0].size(); i++)
+ {
+ Output << "Inv " << zai2[0][i].Z() << " " << zai2[0][i].A() << " " << zai2[0][i].I() << " ";
+ for (int t=StepToSkip; t<vTime.size(); t++)
+ {
+ double Val = zai2[t][i].Prop() * NormalizationFactor;
+ Output << Val << " ";
+ if (t==StepToSkip) NPrinted++;
+
+ }
+ Output << endl;
+ }
+ }
+ else
+ {
+ for(int i=0; i < zai3[0][0].size(); i++)
+ {
+ Output << "Inv " << zai3[0][0][i].Z() << " " << zai3[0][0][i].A() << " " << zai3[0][0][i].I() << " ";
+ for (int t=StepToSkip; t<vTime.size(); t++)
+ {
+ double Val = zai3[t][WantedCell][i].Prop() * NormalizationFactor;
+ Output << Val << " ";
+ if (t==StepToSkip) NPrinted++;
+ }
+ Output << endl;
+ }
+ }
+ }
+
+ if (SumOfCell)
+ {
+ for(int i=0; i< (int)zai_FS_2[0].size(); i++)
+ {
+
+ Output << "XSFis " << zai_FS_2[0][i].Z() << " " << zai_FS_2[0][i].A() << " " << zai_FS_2[0][i].I() << " ";
+ for (int t=StepToSkip; t<vTime.size(); t++)
+ {
+ double Val = zai_FS_2[t][i].Prop();
+ Output << Val << " ";
+
+ }
+ Output << endl;
+ }
+ }
+ else
+ {
+ for(int i=0; i< (int)zai_FS_3[0][0].size(); i++)
+ {
+
+ Output << "XSFis " << zai_FS_3[0][0][i].Z() << " " << zai_FS_3[0][0][i].A() << " " << zai_FS_3[0][0][i].I() << " ";
+ for (int t=StepToSkip; t<vTime.size(); t++)
+ {
+ double Val = zai_FS_3[t][WantedCell][i].Prop();
+ Output << Val << " ";
+ }
+ Output << endl;
+ }
+ }
+
+ if (SumOfCell)
+ {
+ for(int i=0; i< (int)zai_SC_2[0].size(); i++)
+ {
+ Output << "XSCap " << zai_SC_2[0][i].Z() << " " << zai_SC_2[0][i].A() << " " << zai_SC_2[0][i].I() << " ";
+ for (int t=StepToSkip; t<vTime.size(); t++)
+ {
+ double Val = zai_SC_2[t][i].Prop();
+ Output << Val << " ";
+
+ }
+ Output << endl;
+ }
+ }
+ else
+ {
+ for(int i=0; i< (int)zai_SC_3[0][0].size(); i++)
+ {
+ Output << "XSCap " << zai_SC_3[0][0][i].Z() << " " << zai_SC_3[0][0][i].A() << " " << zai_SC_3[0][0][i].I() << " ";
+ for (int t=StepToSkip; t<vTime.size(); t++)
+ {
+ double Val = zai_SC_3[t][WantedCell][i].Prop();
+ Output << Val << " ";
+ }
+ Output << endl;
+ }
+ }
+
+ if (SumOfCell)
+ {
+ for(int i=0; i< (int)zai_SN_2[0].size(); i++)
+ {
+ Output << "XSn2n " << zai_SN_2[0][i].Z() << " " << zai_SN_2[0][i].A() << " " << zai_SN_2[0][i].I() << " ";
+ for (int t=StepToSkip; t<vTime.size(); t++)
+ {
+ double Val = zai_SN_2[t][i].Prop();
+ Output << Val << " ";
+
+ }
+ Output << endl;
+ }
+ }
+ else
+ {
+ for(int i=0; i< (int)zai_SN_3[0][0].size(); i++)
+ {
+ Output << "XSn2n " << zai_SN_3[0][0][i].Z() << " " << zai_SN_3[0][0][i].A() << " " << zai_SN_3[0][0][i].I() << " ";
+ for (int t=StepToSkip; t<vTime.size(); t++)
+ {
+ double Val = zai_SN_3[t][WantedCell][i].Prop();
+ Output << Val << " ";
+ }
+ Output << endl;
+ }
+ }
+ Output.close();
+
+ ofstream OutputInfo(OutDataFileInfo.c_str());
+ OutputInfo << "Reactor " << ReactorType << endl;
+ OutputInfo << "Fueltype " << FuelType << endl;
+ OutputInfo << "CycleTime " << CycleTime << endl;
+ OutputInfo << "AssemblyHeavyMetalMass " << MTotalFissile << " g" << endl;
+ OutputInfo << "ConstantPower " << Power << " W" << endl;
+
+ OutputInfo << "Nnuclei " << NPrinted << endl;
+ OutputInfo << "NormalizationFactor " << NormalizationFactor << endl;
+ OutputInfo << "FinalHeavyMetalMass " << MTotalFissile*NormalizationFactor << " g" << endl;
+
+ OutputInfo.close();
+
+ // =========================================================================================
+ // BYE!
+ // =========================================================================================
+
+ OutputLog << "===================================================" << endl;
+ OutputLog << "---------------------------------------------------" << endl;
+
+ OutputLog << endl << "The database " << OutDataFile << " has been generated..." << endl;
+ OutputLog << "The database information " << OutDataFileInfo << " has been generated..." << endl << endl;
+ OutputLog << NPrinted << " nuclides have been written" << endl << endl;
+
+ OutputLog << "---------------------------------------------------" << endl;
+ OutputLog << "===================================================" << endl;
+
+}
+
+
+/*
+ g++ -o MURE2CLASS MURE2CLASS.cxx
+ */
diff --git a/branches/BaM/Utils/MURE2CLASS/StringLine.hxx b/branches/BaM/Utils/MURE2CLASS/StringLine.hxx
new file mode 100644
index 000000000..17cbf14eb
--- /dev/null
+++ b/branches/BaM/Utils/MURE2CLASS/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 2.01
+*/
+
+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
+ */
+ static 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/branches/BaM/Utils/MURE2CLASS/ZAI.hxx b/branches/BaM/Utils/MURE2CLASS/ZAI.hxx
new file mode 100644
index 000000000..ad18e9652
--- /dev/null
+++ b/branches/BaM/Utils/MURE2CLASS/ZAI.hxx
@@ -0,0 +1,50 @@
+#ifndef _ZAI_
+#define _ZAI_
+
+#include <math.h>
+#include <vector>
+#include <iostream>
+#include <string>
+using namespace std;
+
+class ZAI
+{
+ public:
+
+ ZAI() {fZ=0;fA=0;fI=0;fProp=0;fHalfLifeTime=0;} //!< Default constructor
+ ZAI(int Z, int A, int I) {fZ=Z;fA=A;fI=I;fProp=0.0;fHalfLifeTime=0;}
+ ZAI(int Z, int A, int I, double Prop) {fZ=Z;fA=A;fI=I;fProp=Prop;fHalfLifeTime=0;}
+// ~ZAI();
+
+
+ int Z(){return fZ;} //!< returns the number of protons
+ int A(){return fA;} //!< returns the number of nucleons
+ int I(){return fI;} //!< returns the Isomeric state (Ground State, ith excited)
+ int N(){return fA-fZ;} //!< returns the number of neutrons
+ double Prop(){return fProp;}//!< get the proportion of a ZAI
+
+ void Set(int Z, int A, int I, double p) {fZ=Z;fA=A;fI=I;fProp=p;}
+ void SetZ(int Z){fZ=Z;} //!< Set the proton number
+ void SetA(int A){fA=A;} //!< Set the nucleon number
+ void SetI(int I){fI=I;} //!< Set the isomeric state (0=gs,...)
+ void SetProp(double p){fProp=p;}//!< set the proportion of a ZAI
+
+ void SetHalfLifeTime(double T_12){fHalfLifeTime=T_12;} //!< Set the Decay Half life [s]
+ double GetHalfLifeTime(){return fHalfLifeTime;} //!< returns the Decay Half life [s]
+ void SetDecayConstant(double lambda){fHalfLifeTime=log(2.)/lambda;} //!< Set the Decay Constant
+ double GetDecayConstant(){if (fHalfLifeTime)return log(2.)/fHalfLifeTime;return 0.;} //!< returns the Decay Half life
+
+ void Stable(){fStable=true;} //!< Set that a ZAI is stable to decay
+ bool IsStable(){return fStable;}//!< Return whether a ZAI is stable
+
+ protected:
+
+ int fZ; //!< number of protons
+ int fA; //!< number of nucleons (A=0 means natural isotopes)
+ int fI; //!< Isomeric state (Ground State, ith excited)
+ double fProp; //!< Mass of a ZAI (from the BaseSummary.dat file
+ double fHalfLifeTime; //!< Decay Half life constant
+ bool fStable; //!< whether a ZAI is stable
+
+};
+#endif
diff --git a/branches/BaM/Utils/XSM/MLP/BuildInput/Gene.cxx b/branches/BaM/Utils/XSM/MLP/BuildInput/Gene.cxx
new file mode 100755
index 000000000..1bf019439
--- /dev/null
+++ b/branches/BaM/Utils/XSM/MLP/BuildInput/Gene.cxx
@@ -0,0 +1,497 @@
+/**********************************************************/
+// 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 "Gene.hxx"
+#include <TH1F.h>
+#include <TH2D.h>
+#include <TFile.h>
+#include <TTree.h>
+#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();
+}
+//--------------------------------------------------------------------------------------------------
+//Give a name to a cross section
+//!!!!!!!!!!!!!!!!!!!!!!!!!!//
+// YOU MUST KEEP THE FORMAT :
+// XS_Z_A_I_fis (fission cross section)
+// XS_Z_A_I_cap (n,gamma cross section)
+// XS_Z_A_I_n2n (n,2n cross section)
+//!!!!!!!!!!!!!!!!!!!!!!!!!!//
+string NameXS(ZAI act,string xs)
+{
+ stringstream Name;
+ Name<<"XS_"<<act.Z<<"_"<<act.A<<"_"<<act.I<<"_"<<xs;
+
+ return Name.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;
+/**********************init map********************/
+ map < ZAI,vector<double> > mAllXS;
+ map < ZAI, vector<double> > mAllInventories;
+ for(int act=0;act<int(fAllNuclei.size());act++ )
+ {
+ vector<double> InitVect;
+ for(int Tstep=0 ;Tstep<fNOfTimeStep;Tstep++)
+ {
+ InitVect.push_back(0);
+ }
+ mAllInventories.insert( pair<ZAI,vector<double> >(fAllNuclei[act],InitVect) );
+ }
+
+ for(int act=0;act<int(fAllNuclei.size());act++ )
+ {
+ vector< double> InitVect;
+ for(int xs=0;xs<3;xs++)
+ {
+
+ InitVect.push_back(0);
+
+ }
+
+ mAllXS.insert(pair<ZAI,vector< double> >(fAllNuclei[act], InitVect) );
+ }
+
+/**********************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" );
+
+
+/**********************cross section**************************************************/
+
+ string XSType[3]={"fis","cap","n2n"}; //!!!!!!DO NOT TOUCH THIS
+
+ for(int act=0;act<int(fAllNuclei.size());act++ )
+ {
+ for(int xs=0;xs<3;xs++)
+ {
+ string NamedXS= NameXS(fAllNuclei[act],XSType[xs]);
+ string NameXSBis=NamedXS+"/D";
+
+ if( xs==0 && fXSFis[0][fAllNuclei[act]].size()!=0) //Check if the reaction we want to branch exists
+ if(fXSFis[0][fAllNuclei[act]][0]!=0)
+ fOutT->Branch( NamedXS.c_str() , &mAllXS[fAllNuclei[act]][xs], NameXSBis.c_str() );
+
+ if(xs==1 && fXSCap[0][fAllNuclei[act]].size()!=0)
+ if(fXSCap[0][fAllNuclei[act]][0]!=0)
+ fOutT->Branch( NamedXS.c_str() , &mAllXS[fAllNuclei[act]][xs], NameXSBis.c_str() );
+
+ if(xs==2 && fXSN2N[0][fAllNuclei[act]].size()!=0)
+ if(fXSN2N[0][fAllNuclei[act]][0]!=0)
+ fOutT->Branch( NamedXS.c_str() , &mAllXS[fAllNuclei[act]][xs], NameXSBis.c_str() );
+
+
+ }
+ }
+/**********************FILLING THE TTREE**************************************************/
+
+ //Fill containing all the output of the networks to train
+ ofstream InputNetwork("TrainingInput.cxx");
+
+ 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];
+ for(int act=0;act<int(fAllNuclei.size());act++)
+ {
+
+ if(fXSFis[b][fAllNuclei[act]].size()!=0) // reaction may not be present
+ {
+ if(fXSFis[b][fAllNuclei[act]][Tstep] !=0 )
+ { mAllXS[fAllNuclei[act]][0] = fXSFis[b][fAllNuclei[act]][Tstep];
+ if(b==0 && Tstep==0)
+ InputNetwork<<"OUTPUT.push_back(\""<<NameXS(fAllNuclei[act],XSType[0])<<"\");"<<endl;
+ }
+ }
+ if(fXSCap[b][fAllNuclei[act]].size()!=0)
+ {
+ if(fXSCap[b][fAllNuclei[act]][Tstep] !=0)
+ { mAllXS[fAllNuclei[act]][1] = fXSCap[b][fAllNuclei[act]][Tstep];
+ if(b==0 && Tstep==0)
+ InputNetwork<<"OUTPUT.push_back(\""<<NameXS(fAllNuclei[act],XSType[1])<<"\");"<<endl;
+ }
+ }
+ if( fXSN2N[b][fAllNuclei[act]].size()!=0)
+ {
+ if(fXSN2N[b][fAllNuclei[act]][Tstep] !=0)
+ { mAllXS[fAllNuclei[act]][2] = fXSN2N[b][fAllNuclei[act]][Tstep];
+ if(b==0 && Tstep==0)
+ InputNetwork<<"OUTPUT.push_back(\""<<NameXS(fAllNuclei[act],XSType[2])<<"\");"<<endl;
+ }
+ }
+
+ }
+ 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
+
+ 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());
+
+ /****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;
+
+ //XS
+ map<ZAI, vector <double> > mapFistmp;
+ map<ZAI, vector <double> > mapCaptmp;
+ map<ZAI, vector <double> > mapN2Ntmp;
+ do
+ {
+
+ start = 0;
+ int z;
+ string tmp2 = StringLine::NextWord(line, start, ' ');
+ if (tmp2 == "XSFis")
+ {
+ z = atoi(StringLine::NextWord(line, start, ' ').c_str());
+ int a = atoi(StringLine::NextWord(line, start, ' ').c_str());
+ int i = atoi(StringLine::NextWord(line, start, ' ').c_str());
+
+ if(a!=0 && z!=0)
+ {
+
+ ZAI zaitmp(z, a, i);
+
+ vector<double> XSTime;
+ for(int i = 0; i < (int)vT.size(); i++)
+ { long double q = atof(StringLine::NextWord(line, start, ' ').c_str());
+ XSTime.push_back(q);
+ }
+ mapFistmp.insert( pair<ZAI,vector<double> >(zaitmp,XSTime));
+ }
+ }
+ else if (tmp2 == "XSCap")
+ {
+ z = atoi(StringLine::NextWord(line, start, ' ').c_str());
+ int a = atoi(StringLine::NextWord(line, start, ' ').c_str());
+ int i = atoi(StringLine::NextWord(line, start, ' ').c_str());
+
+ if(a!=0 && z!=0)
+ {
+
+ ZAI zaitmp(z, a, i);
+ vector<double> XSTime;
+ for(int i = 0; i < (int)vT.size(); i++)
+ { long double q = atof(StringLine::NextWord(line, start, ' ').c_str());
+ XSTime.push_back(q);
+ }
+ mapCaptmp.insert( pair<ZAI,vector<double> >(zaitmp,XSTime));
+
+ }
+ }
+ else if (tmp2 == "XSn2n")
+ {
+ z = atoi(StringLine::NextWord(line, start, ' ').c_str());
+ int a = atoi(StringLine::NextWord(line, start, ' ').c_str());
+ int i = atoi(StringLine::NextWord(line, start, ' ').c_str());
+
+ if(a!=0 && z!=0)
+ {
+ ZAI zaitmp(z, a, i);
+ vector<double> XSTime;
+ for(int i = 0; i < (int)vT.size(); i++)
+ { long double q = atof(StringLine::NextWord(line, start, ' ').c_str());
+ XSTime.push_back(q);
+ }
+ mapN2Ntmp.insert( pair<ZAI,vector<double> >(zaitmp,XSTime));
+ }
+ }
+ getline(DecayDB, line);
+ start = 0;
+ tmp2 = StringLine::NextWord(line, start, ' ');
+
+ if(line == "") break;
+ }while (!DecayDB.eof() );
+
+ fXSFis.push_back(mapFistmp);
+ fXSCap.push_back(mapCaptmp);
+ fXSN2N.push_back(mapN2Ntmp);
+
+ 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 Gene Gene.cxx `root-config --cflags` `root-config --libs`
+
+*/
\ No newline at end of file
diff --git a/branches/BaM/Utils/XSM/MLP/BuildInput/Gene.hxx b/branches/BaM/Utils/XSM/MLP/BuildInput/Gene.hxx
new file mode 100755
index 000000000..11d75520f
--- /dev/null
+++ b/branches/BaM/Utils/XSM/MLP/BuildInput/Gene.hxx
@@ -0,0 +1,53 @@
+#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<ZAI> fAllNuclei; //All the nuclei present in the fuel
+
+vector< map < ZAI, vector<double> > > fXSFis; // map of fission cross section fXSFis[NumberOfTheEvolution][ZAI][TimeStep]
+vector< map < ZAI, vector<double> > > fXSCap;
+vector< map < ZAI, vector<double> > > fXSN2N;
+
+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 DumpForTestingNeuron(string filename);
+void DumpInputNeuron(string filename);
\ No newline at end of file
diff --git a/branches/BaM/Utils/XSM/MLP/BuildInput/StringLine.hxx b/branches/BaM/Utils/XSM/MLP/BuildInput/StringLine.hxx
new file mode 100755
index 000000000..da8500370
--- /dev/null
+++ b/branches/BaM/Utils/XSM/MLP/BuildInput/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/branches/BaM/Utils/XSM/MLP/BuildInput/ZAI.hxx b/branches/BaM/Utils/XSM/MLP/BuildInput/ZAI.hxx
new file mode 100755
index 000000000..1ee9136f9
--- /dev/null
+++ b/branches/BaM/Utils/XSM/MLP/BuildInput/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/branches/BaM/Utils/XSM/MLP/Train/EvaluateTrainingCommands.dat b/branches/BaM/Utils/XSM/MLP/Train/EvaluateTrainingCommands.dat
new file mode 100644
index 000000000..7f5793c64
--- /dev/null
+++ b/branches/BaM/Utils/XSM/MLP/Train/EvaluateTrainingCommands.dat
@@ -0,0 +1,5 @@
+in terminal type (with N the number of cross section trained:
+
+root
+.L deviations.C
+for(int i=0;i<N;i++) {stringstream ss;ss<<"Training_output_"<<i<<".root";deviations(ss.str().c_str(),0,kTRUE,kFALSE,kFALSE); }
\ No newline at end of file
diff --git a/branches/BaM/Utils/XSM/MLP/Train/LaunchTraining.sh b/branches/BaM/Utils/XSM/MLP/Train/LaunchTraining.sh
new file mode 100755
index 000000000..b41b38288
--- /dev/null
+++ b/branches/BaM/Utils/XSM/MLP/Train/LaunchTraining.sh
@@ -0,0 +1,21 @@
+#!/bin/bash
+
+# Script to train MLPs
+#@author BaL
+#
+
+echo "--------------------------"
+echo "--- Run Training from MLP $1 to $2 ---"
+echo "--------------------------"
+
+#LigneDeDepart=$((0))
+#NbreSimu=$((700))
+LigneDeDepart=$1
+LigneFinal=$2
+
+echo LigneDeDepart $LigneDeDepart LigneFinal $LigneFinal exclu
+
+for ((ligne=$LigneDeDepart; ligne<$LigneFinal; ligne=ligne+1 )) ;
+do
+ Train_XS $ligne
+done
diff --git a/branches/BaM/Utils/XSM/MLP/Train/Train_XS.cxx b/branches/BaM/Utils/XSM/MLP/Train/Train_XS.cxx
new file mode 100755
index 000000000..2d79051d7
--- /dev/null
+++ b/branches/BaM/Utils/XSM/MLP/Train/Train_XS.cxx
@@ -0,0 +1,185 @@
+/***********************************/
+// Train one MLP "INDICE" from the
+//file ../BuildInput/TrainingInput.root
+//
+//@author Root_tmva_Team modified by BaL
+/***********************************/
+
+#include <cstdlib>
+#include <vector>
+#include <string>
+#include <sstream>
+#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"
+
+#if not defined(__CINT__) || defined(__MAKECINT__)
+#include "TMVA/Tools.h"
+#include "TMVA/Factory.h"
+#endif
+
+using namespace TMVA;
+std::vector<std::string>OUTPUT;
+void LOAD_OUTPUT()
+{
+
+ #include "../BuildInput/TrainingInput.cxx"
+
+}
+void Train_XS_Time(int INDICE)
+{
+
+ //---------------------------------------------------------------
+ // This loads the library
+ TMVA::Tools::Instance();
+
+ // ---------------------------------------------------------------
+
+ std::cout << std::endl;
+ std::cout << "==> Start TMVARegression" << std::endl;
+
+// --------------------------------------------------------------------------------------------------
+// --- Here the preparation phase begins
+
+ // Create a new root OUTPUT file
+ std::stringstream OutTMVA;
+ OutTMVA<<"Training_output_"<<INDICE<<".root";
+ TString outfileName( OutTMVA.str().c_str() );
+ TFile* OUTPUTFile = TFile::Open( outfileName, "UPDATE" );//RECREATE
+
+ // Create the factory object. Later you can choose the methods
+ // whose performance you'd like to investigate. The factory will
+ // then run the performance analysis for you.
+ //
+ // 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" );
+
+//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
+//@@Change : change the each first arguments with the name used in your training sample
+ // factory->AddVariable( "TheNameOfYourNucleusInTheTrainingFile" , "whatever" , "whatever", 'F' );
+ factory->AddVariable( "U5" , "U 235" , "FractionIsotopic", 'F' );
+ factory->AddVariable( "U8" , "U 238" , "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" , "seconds" , 'F' );
+
+ // Add the variable carrying the regression target
+ factory->AddTarget( OUTPUT[INDICE].c_str() ); //The name of the MLP output
+
+ // 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 = "../BuildInput/TrainingInput.root";
+ 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 );
+
+ // 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:
+ // factory->PrepareTrainingAndTestTree( mycut, "SplitMode=random:!V" );
+
+ // ---- Book MVA methods
+ //
+ // please lookup the various method configuration options in the corresponding cxx files, eg:
+ // src/MethoCuts.cxx, etc, or here: http://tmva.sourceforge.net/optionRef.html
+ // it is possible to preset ranges in the option string in which the cut optimisation should be done:
+ // "...:CutRangeMin[2]=-1:CutRangeMax[2]=1"...", where [2] is the third input variable
+
+ std::stringstream Name;
+ Name<< OUTPUT[INDICE];
+ // Neural network (MLP)
+ factory->BookMethod( TMVA::Types::kMLP, Name.str().c_str(), "!H:!V:VarTransform=Norm:NeuronType=tanh:NCycles=20000:HiddenLayers=N,N: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;
+
+}
+int main(int argc, char const *argv[])
+{
+ LOAD_OUTPUT();
+ if(argc != 2 || argv[1] == "-h")
+ {
+ std::cout << "Usage : TrainXS i"<< std::endl;
+ std::cout << "With i the cross section index " << std::endl;
+ std::cout << "File ../BuildInput/TrainingInput.cxx indicates\n indice ranging from 0 to "<< OUTPUT.size()-1 << std::endl;
+ exit(0);
+ }
+
+ Train_XS_Time(atoi(argv[1])) ;
+
+ return 0;
+}
+ /*
+ COMPILATION :
+
+ g++ -o Train_XS `root-config --cflags` Train_XS.cxx `root-config --glibs` -lTMVA
+
+ */
diff --git a/branches/BaM/Utils/XSM/MLP/Train/deviations.C b/branches/BaM/Utils/XSM/MLP/Train/deviations.C
new file mode 100755
index 000000000..db3ab5d1d
--- /dev/null
+++ b/branches/BaM/Utils/XSM/MLP/Train/deviations.C
@@ -0,0 +1,209 @@
+/***********************/
+// Root macro
+//Plot & record distribution of mean (relative or absolute)
+// difference and the standard deviation of this difference
+// between the MLP exectution and a set of known events
+//Save the mean and the standard deviation in file :
+//XS_accuracy.dat
+// This file is a slighly modified version of the file
+// $ROOTSYS/tmva/test/deviations.C
+/***********************/
+#include "TLegend.h"
+#include "TText.h"
+#include "TH2.h"
+#include "tmvaglob.C"
+
+enum HistType { MVAType = 0, ProbaType = 1, RarityType = 2, CompareType = 3 };
+
+// input: - Input file (result from TMVA)
+// - use of TMVA plotting TStyle
+void deviations( TString fin = "TMVAReg.root",
+ HistType htype = MVAType, Bool_t showTarget, Bool_t useTMVAStyle = kTRUE ,Bool_t Absolute = kTRUE )
+{
+ // set style and remove existing canvas'
+ TMVAGlob::Initialize( useTMVAStyle );
+ gStyle->SetNumberContours(999);
+ ofstream Variance("XS_accuracy.dat",ios::app);
+
+ // switches
+ const Bool_t Save_Images = kTRUE;
+
+ // checks if file with name "fin" is already open, and if not opens one
+ TFile* file = TMVAGlob::OpenFile( fin );
+
+ // define Canvas layout here!
+ Int_t xPad = 1; // no of plots in x
+ Int_t yPad = 1; // no of plots in y
+ Int_t noPad = xPad * yPad ;
+ const Int_t width = 650; // size of canvas
+
+ // this defines how many canvases we need
+ TCanvas* c[100];
+
+ // counter variables
+ Int_t countCanvas = 0;
+
+ // search for the right histograms in full list of keys
+ // TList* methods = new TMap();
+
+ TIter next(file->GetListOfKeys());
+ TKey *key(0);
+ while ((key = (TKey*)next())) {
+
+ if (!TString(key->GetName()).BeginsWith("Method_")) continue;
+ if (!gROOT->GetClass(key->GetClassName())->InheritsFrom("TDirectory")) continue;
+
+ TString methodName;
+ TMVAGlob::GetMethodName(methodName,key);
+ cout << "--- Plotting deviation for method: " << methodName << endl;
+
+ TObjString *mN = new TObjString( methodName );
+ TDirectory* mDir = (TDirectory*)key->ReadObj();
+
+ TList* jobNames = new TList();
+
+ TIter keyIt(mDir->GetListOfKeys());
+ TKey *titkey;
+ while ((titkey = (TKey*)keyIt())) {
+
+ if (!gROOT->GetClass(titkey->GetClassName())->InheritsFrom("TDirectory")) continue;
+
+ TDirectory *titDir = (TDirectory *)titkey->ReadObj();
+
+ TObjString *jN = new TObjString( titDir->GetName() );
+ if (!jobNames->Contains( jN )) jobNames->Add( jN );
+ else delete jN;
+
+ TString methodTitle;
+ TMVAGlob::GetMethodTitle(methodTitle,titDir);
+
+ TString hname = "MVA_" + methodTitle;
+ TIter dirKeyIt( titDir->GetListOfKeys() );
+ TKey* dirKey;
+
+ Int_t countPlots = 0;
+ while ((dirKey = (TKey*)dirKeyIt())){
+ if (dirKey->ReadObj()->InheritsFrom("TH2F")) {
+ TString s(dirKey->ReadObj()->GetName());
+ if (s.Contains("_reg_") &&
+ ( (showTarget && s.Contains("_tgt")) || (!showTarget && !s.Contains("_tgt")) ) &&
+ s.Contains( (htype == CompareType ? "train" : "test" )))
+ {
+ c[countCanvas] = new TCanvas( Form("canvas%d", countCanvas+1),
+ Form( "Regression output deviation versus %s for method: %s",
+ (showTarget ? "target" : "input variables"), methodName.Data() ),
+ countCanvas*50+100, (countCanvas+1)*20, width, (Int_t)width*0.72 );
+ c[countCanvas]->SetRightMargin(0.10); // leave space for border
+
+ TH1* htmp = (TH1*)dirKey->ReadObj();
+
+ int NbinX = htmp->GetNbinsX();
+ int NbinY = htmp->GetNbinsY()*2;
+ double XMin = htmp->GetXaxis()->GetBinLowEdge(0);
+ double XMax = htmp->GetXaxis()->GetBinUpEdge(NbinX-1);
+
+ double YMax =-1e10;
+ double YMin =1e10;
+
+ for(int binx=0;binx < htmp->GetNbinsX() ;binx++)
+ {
+ double x = htmp->GetXaxis()->GetBinLowEdge(binx);
+
+ for(int biny=0;biny < htmp->GetNbinsY() ;biny++)
+ {
+ int Bin=htmp->GetBin(binx,biny);
+ double y = htmp->GetYaxis()->GetBinLowEdge(biny);
+
+ if( (y/x > YMax) && (htmp->GetBinContent(Bin) !=0) )
+ { YMax = y/x; }
+
+ if( (y/x < YMin) && (htmp->GetBinContent(Bin) !=0) )
+ { YMin = y/x ; }
+ }
+ }
+
+ double PerOneStep=0.004;
+ NbinY = int( (YMax-YMin)/PerOneStep );
+
+ // TH2D* h=new TH2D("toto", "titi",NbinX, XMin , XMax,NbinY, YMin, YMax );
+
+ TH2D* h=new TH2D("toto", "titi",50, XMin , XMax,20, YMin, YMax );
+
+ for(int binx=0;binx < htmp->GetNbinsX() ;binx++)
+ {
+ double x = htmp->GetXaxis()->GetBinLowEdge(binx);
+
+ for(int biny=0;biny < htmp->GetNbinsY() ;biny++)
+ {
+ int Bin=htmp->GetBin(binx,biny);
+
+ double y = htmp->GetYaxis()->GetBinLowEdge(biny);
+
+
+ h->Fill(x,y/x,htmp->GetBinContent(Bin));
+
+ }
+
+ }
+
+
+
+ h->SetTitle( Form("Relative Output deviation for method: %s (%s sample)",
+ hname.Data(), (htype == CompareType ? "training" : "test" )) );
+
+
+ htmp->SetTitle( Form("Output deviation for method: %s (%s sample)",
+ hname.Data(), (htype == CompareType ? "training" : "test" )) );
+
+
+ if(!Absolute)
+ {
+ h->GetXaxis()->SetTitle(htmp->GetXaxis()->GetTitle());
+ h->GetYaxis()->SetTitle(htmp->GetYaxis()->GetTitle());
+
+ TH1D* h_py = h->ProjectionY();
+ h_py->GetYaxis()->SetTitle(htmp->GetYaxis()->GetTitle());
+
+ h_py->Fit("gaus","M","",YMin,YMax);
+ TF1 *gfit= (TF1 *)h_py->GetFunction("gaus");
+ Variance<<htmp->GetXaxis()->GetTitle()<<" "<< gfit->GetParameter(1) <<" "<< gfit->GetParameter(2)<<endl;
+ // h_py->DrawCopy();
+
+ /* h->Draw("colz");
+ TLine* l = new TLine( h->GetXaxis()->GetXmin(), 0, h->GetXaxis()->GetXmax(), 0 );
+ l->SetLineStyle(2);
+ l->Draw();
+ */
+ }
+ else
+ {
+ htmp->Draw("colz");
+ TLine* l = new TLine( htmp->GetXaxis()->GetXmin(), 0, htmp->GetXaxis()->GetXmax(), 0 );
+ l->SetLineStyle(2);
+ l->Draw();
+ }
+
+ // update and print
+
+ cout << "plotting logo" << endl;
+ TMVAGlob::plot_logo(1.058);
+ c[countCanvas]->Update();
+
+ TString fname = Form( "plots/deviation_%s_%s_%s_c%i",
+ methodName.Data(),
+ (showTarget ? "target" : "vars"),
+ (htype == CompareType ? "training" : "test" ), countPlots );
+ TMVAGlob::imgconv( c[countCanvas], fname );
+
+ //countPlots++;
+ // countCanvas++;
+ }
+ }
+ }
+ }
+ }
+}
+
+
+
+
diff --git a/branches/BaM/Utils/XSM/MLP/Train/tmvaglob.C b/branches/BaM/Utils/XSM/MLP/Train/tmvaglob.C
new file mode 100755
index 000000000..a796dff84
--- /dev/null
+++ b/branches/BaM/Utils/XSM/MLP/Train/tmvaglob.C
@@ -0,0 +1,785 @@
+// global TMVA style settings
+#ifndef TMVA_TMVAGLOB
+#define TMVA_TMVAGLOB
+
+#include <iostream>
+#include <vector>
+
+#include "TPad.h"
+#include "TCanvas.h"
+#include "TColor.h"
+#include "TSystem.h"
+#include "TImage.h"
+#include "TKey.h"
+#include "TH1.h"
+#include "TROOT.h"
+#include "TStyle.h"
+#include "TFile.h"
+#include "TDirectory.h"
+#include "TObjArray.h"
+
+
+#include "RVersion.h"
+
+using std::cout;
+using std::endl;
+
+namespace TMVAGlob {
+
+ // --------- S t y l e ---------------------------
+ static Bool_t UsePaperStyle = 0;
+ // -----------------------------------------------
+
+ enum TypeOfPlot { kId = 0,
+ kNorm,
+ kDecorrelated,
+ kPCA,
+ kGaussDecorr,
+ kNumOfMethods };
+
+ static Int_t c_Canvas = TColor::GetColor( "#f0f0f0" );
+ static Int_t c_FrameFill = TColor::GetColor( "#fffffd" );
+ static Int_t c_TitleBox = TColor::GetColor( "#5D6B7D" );
+ static Int_t c_TitleBorder = TColor::GetColor( "#7D8B9D" );
+ static Int_t c_TitleText = TColor::GetColor( "#FFFFFF" );
+ static Int_t c_SignalLine = TColor::GetColor( "#0000ee" );
+ static Int_t c_SignalFill = TColor::GetColor( "#7d99d1" );
+ static Int_t c_BackgroundLine = TColor::GetColor( "#ff0000" );
+ static Int_t c_BackgroundFill = TColor::GetColor( "#ff0000" );
+ static Int_t c_NovelBlue = TColor::GetColor( "#2244a5" );
+
+ // set the style
+ void SetSignalAndBackgroundStyle( TH1* sig, TH1* bkg, TH1* all = 0 )
+ {
+ //signal
+ // const Int_t FillColor__S = 38 + 150; // change of Color Scheme in ROOT-5.16.
+ // convince yourself with gROOT->GetListOfColors()->Print()
+ Int_t FillColor__S = c_SignalFill;
+ Int_t FillStyle__S = 1001;
+ Int_t LineColor__S = c_SignalLine;
+ Int_t LineWidth__S = 2;
+
+ // background
+ //Int_t icolor = UsePaperStyle ? 2 + 100 : 2;
+ Int_t FillColor__B = c_BackgroundFill;
+ Int_t FillStyle__B = 3554;
+ Int_t LineColor__B = c_BackgroundLine;
+ Int_t LineWidth__B = 2;
+
+ if (sig != NULL) {
+ sig->SetLineColor( LineColor__S );
+ sig->SetLineWidth( LineWidth__S );
+ sig->SetFillStyle( FillStyle__S );
+ sig->SetFillColor( FillColor__S );
+ }
+
+ if (bkg != NULL) {
+ bkg->SetLineColor( LineColor__B );
+ bkg->SetLineWidth( LineWidth__B );
+ bkg->SetFillStyle( FillStyle__B );
+ bkg->SetFillColor( FillColor__B );
+ }
+
+ if (all != NULL) {
+ all->SetLineColor( LineColor__S );
+ all->SetLineWidth( LineWidth__S );
+ all->SetFillStyle( FillStyle__S );
+ all->SetFillColor( FillColor__S );
+ }
+ }
+
+ void SetMultiClassStyle( TObjArray* hists )
+ {
+ //signal
+ // const Int_t FillColor__S = 38 + 150; // change of Color Scheme in ROOT-5.16.
+ // convince yourself with gROOT->GetListOfColors()->Print()
+ //Int_t FillColor__S = c_SignalFill;
+ //Int_t FillStyle__S = 1001;
+ //Int_t LineColor__S = c_SignalLine;
+ //Int_t LineWidth__S = 2;
+
+ // background
+ //Int_t icolor = UsePaperStyle ? 2 + 100 : 2;
+ //Int_t FillColor__B = c_BackgroundFill;
+ //Int_t FillStyle__B = 3554;
+ //Int_t LineColor__B = c_BackgroundLine;
+ //Int_t LineWidth__B = 2;
+
+ Int_t FillColors[10] = {38,2,3,6,7,8,9,11};
+ Int_t LineColors[10] = {4,2,3,6,7,8,9,11};
+ Int_t FillStyles[5] = {1001,3554,3003,3545,0};
+
+ for(Int_t i=0; i<hists->GetEntriesFast(); ++i){
+ ((TH1*)(*hists)[i])->SetFillColor(FillColors[i%10]);
+ ((TH1*)(*hists)[i])->SetFillStyle(FillStyles[i%5]);
+ ((TH1*)(*hists)[i])->SetLineColor(LineColors[i%10]);
+ ((TH1*)(*hists)[i])->SetLineWidth(2);
+ }
+ }
+
+ // set frame styles
+ void SetFrameStyle( TH1* frame, Float_t scale = 1.0 )
+ {
+ frame->SetLabelOffset( 0.012, "X" );// label offset on x axis
+ frame->SetLabelOffset( 0.012, "Y" );// label offset on x axis
+ frame->GetXaxis()->SetTitleOffset( 1.25 );
+ frame->GetYaxis()->SetTitleOffset( 1.22 );
+ frame->GetXaxis()->SetTitleSize( 0.045*scale );
+ frame->GetYaxis()->SetTitleSize( 0.045*scale );
+ Float_t labelSize = 0.04*scale;
+ frame->GetXaxis()->SetLabelSize( labelSize );
+ frame->GetYaxis()->SetLabelSize( labelSize );
+
+ // global style settings
+ gPad->SetTicks();
+ gPad->SetLeftMargin ( 0.108*scale );
+ gPad->SetRightMargin ( 0.050*scale );
+ gPad->SetBottomMargin( 0.120*scale );
+ }
+
+ void SetTMVAStyle() {
+
+ TStyle *TMVAStyle = gROOT->GetStyle("TMVA");
+ if(TMVAStyle!=0) {
+ gROOT->SetStyle("TMVA");
+ return;
+ }
+
+ TMVAStyle = new TStyle(*gROOT->GetStyle("Plain")); // our style is based on Plain
+ TMVAStyle->SetName("TMVA");
+ TMVAStyle->SetTitle("TMVA style based on \"Plain\" with modifications defined in tmvaglob.C");
+ gROOT->GetListOfStyles()->Add(TMVAStyle);
+ gROOT->SetStyle("TMVA");
+
+ TMVAStyle->SetLineStyleString( 5, "[52 12]" );
+ TMVAStyle->SetLineStyleString( 6, "[22 12]" );
+ TMVAStyle->SetLineStyleString( 7, "[22 10 7 10]" );
+
+ // the pretty color palette of old
+ TMVAStyle->SetPalette((UsePaperStyle ? 18 : 1),0);
+
+ // use plain black on white colors
+ TMVAStyle->SetFrameBorderMode(0);
+ TMVAStyle->SetCanvasBorderMode(0);
+ TMVAStyle->SetPadBorderMode(0);
+ TMVAStyle->SetPadColor(0);
+ TMVAStyle->SetFillStyle(0);
+
+ TMVAStyle->SetLegendBorderSize(0);
+
+ // title properties
+ // TMVAStyle->SetTitleW(.4);
+ // TMVAStyle->SetTitleH(.10);
+ // MVAStyle->SetTitleX(.5);
+ // TMVAStyle->SetTitleY(.9);
+ TMVAStyle->SetTitleFillColor( c_TitleBox );
+ TMVAStyle->SetTitleTextColor( c_TitleText );
+ TMVAStyle->SetTitleBorderSize( 1 );
+ TMVAStyle->SetLineColor( c_TitleBorder );
+ if (!UsePaperStyle) {
+ TMVAStyle->SetFrameFillColor( c_FrameFill );
+ TMVAStyle->SetCanvasColor( c_Canvas );
+ }
+
+ // set the paper & margin sizes
+ TMVAStyle->SetPaperSize(20,26);
+ TMVAStyle->SetPadTopMargin(0.10);
+ TMVAStyle->SetPadRightMargin(0.05);
+ TMVAStyle->SetPadBottomMargin(0.11);
+ TMVAStyle->SetPadLeftMargin(0.12);
+
+ // use bold lines and markers
+ TMVAStyle->SetMarkerStyle(21);
+ TMVAStyle->SetMarkerSize(0.3);
+ TMVAStyle->SetHistLineWidth(2);
+ TMVAStyle->SetLineStyleString(2,"[12 12]"); // postscript dashes
+
+ // do not display any of the standard histogram decorations
+ TMVAStyle->SetOptTitle(1);
+ TMVAStyle->SetTitleH(0.052);
+
+ TMVAStyle->SetOptStat(0);
+ TMVAStyle->SetOptFit(0);
+
+ // put tick marks on top and RHS of plots
+ TMVAStyle->SetPadTickX(1);
+ TMVAStyle->SetPadTickY(1);
+
+ }
+
+ void DestroyCanvases()
+ {
+
+ TList* loc = (TList*)gROOT->GetListOfCanvases();
+ TListIter itc(loc);
+ TObject *o(0);
+ while ((o = itc())) delete o;
+ }
+
+ // set style and remove existing canvas'
+ void Initialize( Bool_t useTMVAStyle = kTRUE )
+ {
+ // destroy canvas'
+ DestroyCanvases();
+
+ // set style
+ if (!useTMVAStyle) {
+ gROOT->SetStyle("Plain");
+ gStyle->SetOptStat(0);
+ return;
+ }
+
+ SetTMVAStyle();
+ }
+
+ // checks if file with name "fin" is already open, and if not opens one
+ TFile* OpenFile( const TString& fin )
+ {
+ TFile* file = gDirectory->GetFile();
+ if (file==0 || fin != file->GetName()) {
+ if (file != 0) {
+ gROOT->cd();
+ file->Close();
+ }
+ cout << "--- Opening root file " << fin << " in read mode" << endl;
+ file = TFile::Open( fin, "READ" );
+ }
+ else {
+ file = gDirectory->GetFile();
+ }
+
+ file->cd();
+ return file;
+ }
+
+ // used to create output file for canvas
+ void imgconv( TCanvas* c, const TString & fname )
+ {
+ // return;
+ if (NULL == c) {
+ cout << "*** Error in TMVAGlob::imgconv: canvas is NULL" << endl;
+ }
+ else {
+ // create directory if not existing
+ TString f = fname;
+ TString dir = f.Remove( f.Last( '/' ), f.Length() - f.Last( '/' ) );
+ gSystem->mkdir( dir );
+
+ TString pngName = fname + ".png";
+ TString gifName = fname + ".gif";
+ TString epsName = fname + ".eps";
+ c->cd();
+
+ // create eps (other option: c->Print( epsName ))
+ if (UsePaperStyle) {
+ c->Print(epsName);
+ }
+ else {
+ cout << "--- --------------------------------------------------------------------" << endl;
+ cout << "--- If you want to save the image as eps, gif or png, please comment out " << endl;
+ cout << "--- the corresponding lines (line no. 239-241) in tmvaglob.C" << endl;
+ cout << "--- --------------------------------------------------------------------" << endl;
+ c->Print(epsName);
+ c->Print(pngName);
+ // c->Print(gifName);
+ }
+ }
+ }
+
+ TImage * findImage(const char * imageName)
+ {
+ // looks for the image in macropath
+ TString macroPath(gROOT->GetMacroPath()); // look for the image in here
+ Ssiz_t curIndex(0);
+ TImage *img(0);
+ while(1) {
+ Ssiz_t pathStart = curIndex;
+ curIndex = macroPath.Index(":",curIndex);
+ Ssiz_t pathEnd = (curIndex==-1)?macroPath.Length():curIndex;
+ TString path(macroPath(pathStart,pathEnd-pathStart));
+
+ gSystem->ExpandPathName(path);
+ const char* fullName = Form("%s/%s", path.Data(), imageName);
+
+ Bool_t fileFound = ! gSystem->AccessPathName(fullName);
+
+ if(fileFound) {
+ img = TImage::Open(fullName);
+ break;
+ }
+ if(curIndex==-1) break;
+ curIndex++;
+ }
+ return img;
+ }
+
+ void plot_logo( Float_t v_scale = 1.0, Float_t skew = 1.0 )
+ {
+
+ TImage *img = findImage("tmva_logo.gif");
+ if (!img) {
+ cout << "+++ Could not open image tmva_logo.gif" << endl;
+ return;
+ }
+
+ img->SetConstRatio(kFALSE);
+ UInt_t h_ = img->GetHeight();
+ UInt_t w_ = img->GetWidth();
+
+ Float_t r = w_/h_;
+ gPad->Update();
+ Float_t rpad = Double_t(gPad->VtoAbsPixel(0) - gPad->VtoAbsPixel(1))/(gPad->UtoAbsPixel(1) - gPad->UtoAbsPixel(0));
+ r *= rpad;
+
+ Float_t d = 0.055;
+ // absolute coordinates
+ Float_t x1R = 1 - gStyle->GetPadRightMargin();
+ Float_t y1B = 1 - gStyle->GetPadTopMargin()+.01; // we like the logo to sit a bit above the histo
+
+ Float_t x1L = x1R - d*r/skew;
+ Float_t y1T = y1B + d*v_scale*skew;
+ if (y1T>0.99) y1T = 0.99;
+
+ TPad *p1 = new TPad("imgpad", "imgpad", x1L, y1B, x1R, y1T );
+ p1->SetRightMargin(0);
+ p1->SetBottomMargin(0);
+ p1->SetLeftMargin(0);
+ p1->SetTopMargin(0);
+ p1->Draw();
+
+ Int_t xSizeInPixel = p1->UtoAbsPixel(1) - p1->UtoAbsPixel(0);
+ Int_t ySizeInPixel = p1->VtoAbsPixel(0) - p1->VtoAbsPixel(1);
+ if (xSizeInPixel<=25 || ySizeInPixel<=25) {
+ delete p1;
+ return; // ROOT doesn't draw smaller than this
+ }
+
+ p1->cd();
+ img->Draw();
+ }
+
+ void NormalizeHist( TH1* h )
+ {
+ if (h==0) return;
+ if (h->GetSumw2N() == 0) h->Sumw2();
+ if(h->GetSumOfWeights()!=0) {
+ Float_t dx = (h->GetXaxis()->GetXmax() - h->GetXaxis()->GetXmin())/h->GetNbinsX();
+ h->Scale( 1.0/h->GetSumOfWeights()/dx );
+ }
+ }
+ void NormalizeHists( TH1* sig, TH1* bkg = 0 )
+ {
+ if (sig->GetSumw2N() == 0) sig->Sumw2();
+ if (bkg && bkg->GetSumw2N() == 0) bkg->Sumw2();
+
+ if(sig->GetSumOfWeights()!=0) {
+ Float_t dx = (sig->GetXaxis()->GetXmax() - sig->GetXaxis()->GetXmin())/sig->GetNbinsX();
+ sig->Scale( 1.0/sig->GetSumOfWeights()/dx );
+ }
+ if (bkg != 0 && bkg->GetSumOfWeights()!=0) {
+ Float_t dx = (bkg->GetXaxis()->GetXmax() - bkg->GetXaxis()->GetXmin())/bkg->GetNbinsX();
+ bkg->Scale( 1.0/bkg->GetSumOfWeights()/dx );
+ }
+ }
+
+ // the following are tools to help handling different methods and titles
+
+
+ void GetMethodName( TString & name, TKey * mkey ) {
+ if (mkey==0) return;
+ name = mkey->GetName();
+ name.ReplaceAll("Method_","");
+ }
+
+ void GetMethodTitle( TString & name, TKey * ikey ) {
+ if (ikey==0) return;
+ name = ikey->GetName();
+ }
+
+ void GetMethodName( TString & name, TDirectory * mdir ) {
+ if (mdir==0) return;
+ name = mdir->GetName();
+ name.ReplaceAll("Method_","");
+ }
+
+ void GetMethodTitle( TString & name, TDirectory * idir ) {
+ if (idir==0) return;
+ name = idir->GetName();
+ }
+
+ TKey *NextKey( TIter & keyIter, TString className) {
+ TKey *key=(TKey *)keyIter.Next();
+ TKey *rkey=0;
+ Bool_t loop=(key!=0);
+ //
+ while (loop) {
+ TClass *cl = gROOT->GetClass(key->GetClassName());
+ if (cl->InheritsFrom(className.Data())) {
+ loop = kFALSE;
+ rkey = key;
+ } else {
+ key = (TKey *)keyIter.Next();
+ if (key==0) loop = kFALSE;
+ }
+ }
+ return rkey;
+ }
+
+ UInt_t GetListOfKeys( TList& keys, TString inherits, TDirectory *dir=0 )
+ {
+ // get a list of keys with a given inheritance
+ // the list contains TKey objects
+ if (dir==0) dir = gDirectory;
+ TIter mnext(dir->GetListOfKeys());
+ TKey *mkey;
+ keys.Clear();
+ keys.SetOwner(kFALSE);
+ UInt_t ni=0;
+ while ((mkey = (TKey*)mnext())) {
+ // make sure, that we only look at TDirectory with name Method_<xxx>
+ TClass *cl = gROOT->GetClass(mkey->GetClassName());
+ if (cl->InheritsFrom(inherits)) {
+ keys.Add(mkey);
+ ni++;
+ }
+ }
+ return ni;
+ }
+
+ Int_t GetNumberOfTargets( TDirectory *dir )
+ {
+ if (!dir) {
+ cout << "tmvaglob::GetNumberOfTargets is called with *dir==NULL :( " << endl;
+ return 0;
+ }
+ TIter next(dir->GetListOfKeys());
+ TKey* key = 0;
+ Int_t noTrgts = 0;
+
+ while ((key = (TKey*)next())) {
+ if (key->GetCycle() != 1) continue;
+ if (TString(key->GetName()).Contains("__Regression_target")) noTrgts++;
+ }
+ return noTrgts;
+ }
+
+ Int_t GetNumberOfInputVariables( TDirectory *dir )
+ {
+ TIter next(dir->GetListOfKeys());
+ TKey* key = 0;
+ Int_t noVars = 0;
+
+ while ((key = (TKey*)next())) {
+ if (key->GetCycle() != 1) continue;
+
+ // count number of variables (signal is sufficient), exclude target(s)
+ if (TString(key->GetName()).Contains("__Signal") || (TString(key->GetName()).Contains("__Regression") && !(TString(key->GetName()).Contains("__Regression_target")))) noVars++;
+ }
+
+ return noVars;
+ }
+
+ std::vector<TString> GetInputVariableNames(TDirectory *dir )
+ {
+ TIter next(dir->GetListOfKeys());
+ TKey* key = 0;
+ //set<std::string> varnames;
+ std::vector<TString> names;
+
+ while ((key = (TKey*)next())) {
+ if (key->GetCycle() != 1) continue;
+ TClass *cl = gROOT->GetClass(key->GetClassName());
+ if (!cl->InheritsFrom("TH1")) continue;
+ TString name(key->GetName());
+ Int_t pos = name.First("__");
+ name.Remove(pos);
+ Bool_t hasname = false;
+ std::vector<TString>::const_iterator iter = names.begin();
+ while(iter != names.end()){
+ if(name.CompareTo(*iter)==0)
+ hasname=true;
+ iter++;
+ }
+ if(!hasname)
+ names.push_back(name);
+ }
+ return names;
+ }
+
+ Int_t GetNumberOfInputVariablesMultiClass( TDirectory *dir ){
+ std::vector<TString> names(GetInputVariableNames(dir));
+ return names.end() - names.begin();
+ }
+
+ std::vector<TString> GetClassNames(TDirectory *dir )
+ {
+
+ TIter next(dir->GetListOfKeys());
+ TKey* key = 0;
+ //set<std::string> varnames;
+ std::vector<TString> names;
+
+ while ((key = (TKey*)next())) {
+ if (key->GetCycle() != 1) continue;
+ TClass *cl = gROOT->GetClass(key->GetClassName());
+ if (!cl->InheritsFrom("TH1")) continue;
+ TString name(key->GetName());
+ name.ReplaceAll("_Deco","");
+ name.ReplaceAll("_Gauss","");
+ name.ReplaceAll("_PCA","");
+ name.ReplaceAll("_Id","");
+ name.ReplaceAll("_vs_","");
+ char c = '_';
+ Int_t pos = name.Last(c);
+ name.Remove(0,pos+1);
+
+ /*Int_t pos = name.First("__");
+ name.Remove(0,pos+2);
+ char c = '_';
+ pos = name.Last(c);
+ name.Remove(pos);
+ if(name.Contains("Gauss")){
+ pos = name.Last(c);
+ name.Remove(pos);
+ }
+ pos = name.Last(c);
+ if(pos!=-1)
+ name.Remove(0,pos+1);
+ */
+ Bool_t hasname = false;
+ std::vector<TString>::const_iterator iter = names.begin();
+ while(iter != names.end()){
+ if(name.CompareTo(*iter)==0)
+ hasname=true;
+ iter++;
+ }
+ if(!hasname)
+ names.push_back(name);
+ }
+ return names;
+ }
+
+
+ TKey* FindMethod( TString name, TDirectory *dir=0 )
+ {
+ // find the key for a method
+ if (dir==0) dir = gDirectory;
+ TIter mnext(dir->GetListOfKeys());
+ TKey *mkey;
+ TKey *retkey=0;
+ Bool_t loop=kTRUE;
+ while (loop) {
+ mkey = (TKey*)mnext();
+ if (mkey==0) {
+ loop = kFALSE;
+ }
+ else {
+ TString clname = mkey->GetClassName();
+ TClass *cl = gROOT->GetClass(clname);
+ if (cl->InheritsFrom("TDirectory")) {
+ TString mname = mkey->GetName(); // method name
+ TString tname = "Method_"+name; // target name
+ if (mname==tname) { // target found!
+ loop = kFALSE;
+ retkey = mkey;
+ }
+ }
+ }
+ }
+ return retkey;
+ }
+
+ Bool_t ExistMethodName( TString name, TDirectory *dir=0 )
+ {
+ // find the key for a method
+ if (dir==0) dir = gDirectory;
+ TIter mnext(dir->GetListOfKeys());
+ TKey *mkey;
+ Bool_t loop=kTRUE;
+ while (loop) {
+ mkey = (TKey*)mnext();
+ if (mkey==0) {
+ loop = kFALSE;
+ }
+ else {
+ TString clname = mkey->GetClassName();
+ TString keyname = mkey->GetName();
+ TClass *cl = gROOT->GetClass(clname);
+ if (keyname.Contains("Method") && cl->InheritsFrom("TDirectory")) {
+
+ TDirectory* d_ = (TDirectory*)dir->Get( keyname );
+ if (!d_) {
+ cout << "HUUUGE TROUBLES IN TMVAGlob::ExistMethodName() --> contact authors" << endl;
+ return kFALSE;
+ }
+
+ TIter mnext_(d_->GetListOfKeys());
+ TKey *mkey_;
+ while ((mkey_ = (TKey*)mnext_())) {
+ TString clname_ = mkey_->GetClassName();
+ TClass *cl_ = gROOT->GetClass(clname_);
+ if (cl_->InheritsFrom("TDirectory")) {
+ TString mname = mkey_->GetName(); // method name
+ if (mname==name) { // target found!
+ return kTRUE;
+ }
+ }
+ }
+ }
+ }
+ }
+ return kFALSE;
+ }
+
+ UInt_t GetListOfMethods( TList & methods, TDirectory *dir=0 )
+ {
+ // get a list of methods
+ // the list contains TKey objects
+ if (dir==0) dir = gDirectory;
+ TIter mnext(dir->GetListOfKeys());
+ TKey *mkey;
+ methods.Clear();
+ methods.SetOwner(kFALSE);
+ UInt_t ni=0;
+ while ((mkey = (TKey*)mnext())) {
+ // make sure, that we only look at TDirectory with name Method_<xxx>
+ TString name = mkey->GetClassName();
+ TClass *cl = gROOT->GetClass(name);
+ if (cl->InheritsFrom("TDirectory")) {
+ if (TString(mkey->GetName()).BeginsWith("Method_")) {
+ methods.Add(mkey);
+ ni++;
+ }
+ }
+ }
+ cout << "--- Found " << ni << " classifier types" << endl;
+ return ni;
+ }
+
+ UInt_t GetListOfJobs( TFile* file, TList& jobdirs)
+ {
+ // get a list of all jobs in all method directories
+ // based on ideas by Peter and Joerg found in macro deviations.C
+ TIter next(file->GetListOfKeys());
+ TKey *key(0);
+ while ((key = (TKey*)next())) {
+
+ if (TString(key->GetName()).BeginsWith("Method_")) {
+ if (gROOT->GetClass(key->GetClassName())->InheritsFrom("TDirectory")) {
+
+ TDirectory* mDir = (TDirectory*)key->ReadObj();
+
+ TIter keyIt(mDir->GetListOfKeys());
+ TKey *jobkey;
+ while ((jobkey = (TKey*)keyIt())) {
+ if (!gROOT->GetClass(jobkey->GetClassName())->InheritsFrom("TDirectory")) continue;
+
+ TDirectory *jobDir = (TDirectory *)jobkey->ReadObj();
+ cout << "jobdir name " << jobDir->GetName() << endl;
+ jobdirs.Add(jobDir);
+ }
+ }
+ }
+ }
+ return jobdirs.GetSize();
+ }
+
+ UInt_t GetListOfTitles( TDirectory *rfdir, TList & titles )
+ {
+ // get a list of titles (i.e TDirectory) given a method dir
+ UInt_t ni=0;
+ if (rfdir==0) return 0;
+ TList *keys = rfdir->GetListOfKeys();
+ if (keys==0) {
+ cout << "+++ Directory '" << rfdir->GetName() << "' contains no keys" << endl;
+ return 0;
+ }
+ //
+ TIter rfnext(rfdir->GetListOfKeys());
+ TKey *rfkey;
+ titles.Clear();
+ titles.SetOwner(kFALSE);
+ while ((rfkey = (TKey*)rfnext())) {
+ // make sure, that we only look at histograms
+ TClass *cl = gROOT->GetClass(rfkey->GetClassName());
+ if (cl->InheritsFrom("TDirectory")) {
+ titles.Add(rfkey);
+ ni++;
+ }
+ }
+ cout << "--- Found " << ni << " instance(s) of the method " << rfdir->GetName() << endl;
+ return ni;
+ }
+
+ UInt_t GetListOfTitles( TString & methodName, TList & titles, TDirectory *dir=0 )
+ {
+ // get the list of all titles for a given method
+ // if the input dir is 0, gDirectory is used
+ // returns a list of keys
+ UInt_t ni=0;
+ if (dir==0) dir = gDirectory;
+ TDirectory* rfdir = (TDirectory*)dir->Get( methodName );
+ if (rfdir==0) {
+ cout << "+++ Could not locate directory '" << methodName << endl;
+ return 0;
+ }
+
+ return GetListOfTitles( rfdir, titles );
+
+ TList *keys = rfdir->GetListOfKeys();
+ if (keys==0) {
+ cout << "+++ Directory '" << methodName << "' contains no keys" << endl;
+ return 0;
+ }
+ //
+ TIter rfnext(rfdir->GetListOfKeys());
+ TKey *rfkey;
+ titles.Clear();
+ titles.SetOwner(kFALSE);
+ while ((rfkey = (TKey*)rfnext())) {
+ // make sure, that we only look at histograms
+ TClass *cl = gROOT->GetClass(rfkey->GetClassName());
+ if (cl->InheritsFrom("TDirectory")) {
+ titles.Add(rfkey);
+ ni++;
+ }
+ }
+ cout << "--- Found " << ni << " instance(s) of the method " << methodName << endl;
+ return ni;
+ }
+
+ TDirectory *GetInputVariablesDir( TMVAGlob::TypeOfPlot type, TDirectory *dir=0 )
+ {
+ // get the InputVariables directory
+ const TString directories[TMVAGlob::kNumOfMethods] = { "InputVariables_Id",
+ "InputVariables_Deco",
+ "InputVariables_PCA",
+ "InputVariables_Gauss_Deco" };
+ if (dir==0) dir = gDirectory;
+
+ // get top dir containing all hists of the variables
+ dir = (TDirectory*)gDirectory->Get( directories[type] );
+ if (dir==0) {
+ cout << "+++ Could not locate input variable directory '" << directories[type] << endl;
+ return 0;
+ }
+ return dir;
+ }
+
+ TDirectory *GetCorrelationPlotsDir( TMVAGlob::TypeOfPlot type, TDirectory *dir=0 )
+ {
+ // get the CorrelationPlots directory
+ if (dir==0) dir = GetInputVariablesDir( type, 0 );
+ if (dir==0) return 0;
+ //
+ TDirectory* corrdir = (TDirectory*)dir->Get( "CorrelationPlots" );
+ if (corrdir==0) {
+ cout << "+++ Could not find CorrelationPlots directory 'CorrelationPlots'" << endl;
+ return 0;
+ }
+ return corrdir;
+ }
+
+}
+
+#endif
diff --git a/branches/BaM/Utils/rootlogon.C b/branches/BaM/Utils/rootlogon.C
new file mode 100755
index 000000000..6bc2ed9c5
--- /dev/null
+++ b/branches/BaM/Utils/rootlogon.C
@@ -0,0 +1,5 @@
+{
+ #include "TMatrix"
+ gSystem.Load("$CLASS_lib/libCLASSpkg_root.so");
+
+}
--
GitLab