diff --git a/NPLib/CMakeLists.txt b/NPLib/CMakeLists.txt
index dae59391ea14715418ae5877ced3730b23deaf82..6f39ad17552faa35094855162c026ad87f51dbd0 100644
--- a/NPLib/CMakeLists.txt
+++ b/NPLib/CMakeLists.txt
@@ -11,11 +11,11 @@ include("${CMAKE_BINARY_DIR}/ressources/CMake/NPTool_CMake_Preamble.cmake")
 
 
 # Major change in the Core/Physics (affecting the detector/analysis/simulation)
-set(NPLIB_VERSION_MAJOR 2)
+set(NPLIB_VERSION_MAJOR 3)
 # Minor change in the Core/Physics (not affecting any other part)
-set(NPLIB_VERSION_MINOR 2)
+set(NPLIB_VERSION_MINOR 0)
 # Change any of the detector in NPA
-set(NPLIB_VERSION_DETA 45)
+set(NPLIB_VERSION_DETA 0)
 
 #activate Multithreading (on by default)
 if(NOT DEFINED NPMULTITHREADING)
diff --git a/NPLib/Core/NPDetectorManager.cxx b/NPLib/Core/NPDetectorManager.cxx
index 490684d236f9e0154f5470e7c5c3f78ae931845c..0b000bc0b6c24a3aa33b0be6b732b639c41145a1 100644
--- a/NPLib/Core/NPDetectorManager.cxx
+++ b/NPLib/Core/NPDetectorManager.cxx
@@ -88,6 +88,7 @@ NPL::DetectorManager::DetectorManager(){
   m_ShieldFrontRadius = 0 ; 
   m_ShieldBackRadius = 0 ;
   m_ShieldMaterial = "" ;
+  m_RootOutput=RootOutput::getInstance();
 }
 
 
@@ -266,7 +267,6 @@ void NPL::DetectorManager::BuildPhysicalEvent(){
   for (it =begin; it != end; ++it) {
     (it->second->*m_ClearEventPhysicsPtr)();
     (it->second->*m_BuildPhysicalPtr)();
-    if(m_FillSpectra){
       (it->second->*m_FillSpectra)();
       if(m_CheckSpectra)
         (it->second->*m_CheckSpectra)();
@@ -498,6 +498,12 @@ void NPL::DetectorManager::CheckSpectraServer(){
     std::cout <<"WARNING: requesting to check spectra server, which is not started" << std::endl; 
 
 }
+
+////////////////////////////////////////////////////////////////////////////////
+void  NPL::DetectorManager::FillOutputTree(){
+  m_RootOutput->Fill();
+}
+
 ////////////////////////////////////////////////////////////////////////////////
 bool NPL::DetectorManager::IsCryogenic() const {return  m_CryoTarget;};
 double NPL::DetectorManager::GetTargetThickness() const {return m_TargetThickness;};
diff --git a/NPLib/Core/NPDetectorManager.h b/NPLib/Core/NPDetectorManager.h
index 07e280bcc371132daff707f6cfedbc7dbcb63e97..f6f94d0c0f59fd47e66bddc300e5e7b5e4835766 100644
--- a/NPLib/Core/NPDetectorManager.h
+++ b/NPLib/Core/NPDetectorManager.h
@@ -23,6 +23,7 @@
 //   NPL
 #include "NPVDetector.h"
 #include "NPSpectraServer.h"
+#include "RootOutput.h"
 // ROOT
 #include "TH1.h"
 
@@ -67,6 +68,13 @@ namespace NPL{
     private :
       NPL::SpectraServer* m_SpectraServer;
 
+
+    private:// root outpu
+      RootOutput* m_RootOutput;
+
+    public:
+      void       FillOutputTree();
+
     private:   
       // The std::map containning all detectors
       // Using a Map one can access to any detector using its name
diff --git a/NPLib/Core/NPOptionManager.cxx b/NPLib/Core/NPOptionManager.cxx
index a38c8569ea4d021e0d6b2a4762ca956a6a86f817..450c4ddf7274b1bd5129a3c9720d45221b128d93 100644
--- a/NPLib/Core/NPOptionManager.cxx
+++ b/NPLib/Core/NPOptionManager.cxx
@@ -106,6 +106,8 @@ void NPOptionManager::ReadTheInputArgument(int argc, char** argv){
   fLastPhyFile = false;
   fLastResFile = false;
   fLastAnyFile = false;
+  fIsAnalysis  = false;
+  fIsSimulation= false;
   fVerboseLevel               = 1;
   fNumberOfEntryToAnalyse     = -1;
 	fFirstEntryToAnalyse        = 0;
@@ -114,7 +116,6 @@ void NPOptionManager::ReadTheInputArgument(int argc, char** argv){
   fDisableAllBranchOption = false;
   fInputPhysicalTreeOption = false;
   fGenerateHistoOption = false ;
-  fPROOFMode = false;
   fCircularTree = false;
   fOnline = false;
   fG4BatchMode = false;
@@ -158,6 +159,8 @@ void NPOptionManager::ReadTheInputArgument(int argc, char** argv){
     else if (argument == "-T" && argc >= i + 2)                   { std::string file = argv[++i] ; std::string tree = argv[++i]; CreateRunToTreatFile(file,tree);}
 
     else if (argument == "--cal" && argc >= i + 1)                fCalibrationFileName = argv[++i] ;
+    
+    else if (argument == "-S" && argc >= i + 1)                   fIsSplit=true; 
 
     else if (argument == "-C" && argc >= i + 1)                   fCalibrationFileName = argv[++i] ;
 
@@ -185,8 +188,6 @@ void NPOptionManager::ReadTheInputArgument(int argc, char** argv){
 
     else if (argument == "--generate-histo")                      fGenerateHistoOption = true ;
 
-    else if (argument == "--proof")                               fPROOFMode = true ;
-
     else if (argument == "-L")                                    fNumberOfEntryToAnalyse = atoi(argv[++i]) ;
 
     else if (argument == "--random-seed")                         fRandomSeed = atoi(argv[++i]) ;
@@ -444,12 +445,12 @@ void NPOptionManager::DisplayHelp(){
   std::cout << "\t--output -O <arg>\t\tSet arg as the Output File Name (output tree)" << std::endl ;
   std::cout << "\t--tree-name <arg>\t\tSet arg as the Output Tree Name " << std::endl ;
   std::cout << "\t--definition <definition> \tAdd <definition> to the list of definition" << std::endl  ;
+  std::cout << "\t-S \t\t\t\tOne tree output per detector" << std::endl  ;
   std::cout << "\t--verbose -V <arg>\t\tSet the verbose level, 0 for nothing, 1 for normal printout."<<std::endl;
 	std::cout  << "\t\t\t\t\tError and warning are not affected" << std::endl ;
   std::cout << std::endl << "NPAnalysis only:"<<std::endl;
   std::cout << "\t--run -R <arg>\t\t\tSet arg as the run to read file list" << std::endl  ;
   std::cout << "\t-T <name> <file>\t\tTree <name> from root file <file>" << std::endl  ;
-
   std::cout << "\t--cal -C <arg>\t\t\tSet arg as the calibration file list" << std::endl ;
   std::cout << "\t--disable-branch\t\tDisable of branch of Input tree except the one of the detector (faster)" << std::endl  ;
   std::cout << "\t--generate-histo -GH\t\tInstantiate the T*Spectra class of each detector" << std::endl ;
diff --git a/NPLib/Core/NPOptionManager.h b/NPLib/Core/NPOptionManager.h
index a793188b715dc4efbb37daae2c75411c2b092900..8268a4ba99d5c3bbbaccb0d2e6cd10386843ac48 100644
--- a/NPLib/Core/NPOptionManager.h
+++ b/NPLib/Core/NPOptionManager.h
@@ -103,9 +103,12 @@ class NPOptionManager{
       bool   GetGenerateHistoOption()      {return fGenerateHistoOption;}
       bool   GetCheckHistoOption()         {return fCheckHistoOption;}
       bool   GetOnline()                   {return fOnline;}
-      bool   GetPROOF()                    {return fPROOFMode;}
       bool   GetG4BatchMode()              {return fG4BatchMode;}
-      bool   GetCircularTree()                 {return fCircularTree;}
+      bool   GetCircularTree()             {return fCircularTree;}
+      bool   IsAnalysis()                  {return fIsAnalysis;};
+      bool   IsSimulation()                {return fIsSimulation;}
+      bool   IsSplit()                     {return fIsSplit;}
+
       int    GetVerboseLevel()             {return fVerboseLevel;}
       int    GetNumberOfEntryToAnalyse()   {return fNumberOfEntryToAnalyse;} 
       int    GetFirstEntryToAnalyse()      {return fFirstEntryToAnalyse;} 
@@ -121,7 +124,8 @@ class NPOptionManager{
       void SetDetectorFile(const std::string& name)  {fDetectorFileName = name;CheckDetectorConfiguration();}
       void SetRunToReadFile(const std::string& name) {fRunToReadFileName = name;}
       void SetVerboseLevel(int VerboseLevel)         {fVerboseLevel = VerboseLevel;}
- 
+      void SetIsAnalysis(bool val=true){fIsAnalysis=val;};
+      void SetIsSimulation(bool val=true){fIsSimulation=val;}
 
    public: // user definition
       bool HasDefinition(std::string def) {return(fDefinition.find(def)!=fDefinition.end());}
@@ -143,17 +147,19 @@ class NPOptionManager{
       std::string fCalibrationFileName;
       std::string fOutputFileName;
       std::string fOutputTreeName;
+      bool   fIsSplit; // One tree per detector
       bool   fDisableAllBranchOption;
       bool   fInputPhysicalTreeOption;
       bool   fGenerateHistoOption;
       bool   fCheckHistoOption;
       bool   fOnline; // true if spectra server is started
-      bool   fPROOFMode; // if true, the system run in a pROOF environment
       bool   fLastSimFile;
       bool   fLastPhyFile;
       bool   fLastResFile;
       bool   fLastAnyFile;
       bool   fCircularTree;
+      bool   fIsAnalysis;
+      bool   fIsSimulation;
       int    fVerboseLevel; // 0 for not talk, 1 for talking
       int    fNumberOfEntryToAnalyse; // use to limit the number of analysed in NPA
       int    fFirstEntryToAnalyse; // use to set the first event analysed in NPA (total: fFirstEntryToAnalyse -> fFirstEntryToAnalyse + fNumberOfEntryToAnalyse)
diff --git a/NPLib/Core/RootInput.cxx b/NPLib/Core/RootInput.cxx
index e14e60357e879a8ed0d0953e44e4f9b359695042..b92a4f76682c5d827a7de9080f7f45505ec230df 100644
--- a/NPLib/Core/RootInput.cxx
+++ b/NPLib/Core/RootInput.cxx
@@ -1,5 +1,5 @@
 /*****************************************************************************
- * Copyright (C) 2009-2016   this file is part of the NPTool Project         *
+ * Copyright (C) 2009-2021   this file is part of the NPTool Project         *
  *                                                                           *
  * For the licensing terms see $NPTOOL/Licence/NPTool_Licence                *
  * For the list of contributors see $NPTOOL/Licence/Contributors             *
@@ -9,7 +9,7 @@
  * Original Author: N. de Sereville  contact address: deserevi@ipno.in2p3.fr *
  *                                                                           *
  * Creation Date  : 21/07/09                                                 *
- * Last update    :                                                          *
+ * Last update    : 10/05/21 Add support for Split tree input (A. Matta)     *
  *---------------------------------------------------------------------------*
  * Decription: This class is a singleton class which deals with the ROOT     *
  *             input file and tree both for NPSimulation and NPAnalysis.     *
diff --git a/NPLib/Core/RootInput.h b/NPLib/Core/RootInput.h
index 51c4799566829323cc2d5ed4278d7cc1847f4d3d..ab64a9a73475effe0b848e075f050c715f61df4d 100644
--- a/NPLib/Core/RootInput.h
+++ b/NPLib/Core/RootInput.h
@@ -1,7 +1,7 @@
 #ifndef ROOTINPUT_HH
 #define ROOTINPUT_HH
 /*****************************************************************************
- * Copyright (C) 2009-2016    this file is part of the NPTool Project        *
+ * Copyright (C) 2009-2021   this file is part of the NPTool Project         *
  *                                                                           *
  * For the licensing terms see $NPTOOL/Licence/NPTool_Licence                *
  * For the list of contributors see $NPTOOL/Licence/Contributors             *
@@ -11,7 +11,7 @@
  * Original Author: N. de Sereville  contact address: deserevi@ipno.in2p3.fr *
  *                                                                           *
  * Creation Date  : 21/07/09                                                 *
- * Last update    :                                                          *
+ * Last update    : 10/05/21 Add support for Split tree input (A. Matta)     *
  *---------------------------------------------------------------------------*
  * Decription: This class is a singleton class which deals with the ROOT     *
  *             input file and tree both for NPSimulation and NPAnalysis.     *
diff --git a/NPLib/Core/RootOutput.cxx b/NPLib/Core/RootOutput.cxx
index 80400e84fe0546ab6966ad179f6d68c7761d6f83..d5fab575b67d7f910a5c0988eb78a68dac694efc 100644
--- a/NPLib/Core/RootOutput.cxx
+++ b/NPLib/Core/RootOutput.cxx
@@ -1,5 +1,5 @@
 /*****************************************************************************
- * Copyright (C) 2009-2016   this file is part of the NPTool Project         *
+ * Copyright (C) 2009-2021   this file is part of the NPTool Project         *
  *                                                                           *
  * For the licensing terms see $NPTOOL/Licence/NPTool_Licence                *
  * For the list of contributors see $NPTOOL/Licence/Contributors             *
@@ -9,13 +9,14 @@
  * Original Author: N. de Sereville  contact address: deserevi@ipno.in2p3.fr *
  *                                                                           *
  * Creation Date  : 21/07/09                                                 *
- * Last update    : 03/02/11                                                 *
+ * Last update    : 10/05/21                                                 *
  *---------------------------------------------------------------------------*
  * Decription: This class is a singleton class which deals with the ROOT     *
  *             output file and tree both for NPSimulation and NPAnalysis.    *
  *---------------------------------------------------------------------------*
  * Comment:                                                                  *
  *   + 03/02/11: Add support for TAsciiFile objects (N. de Sereville)        *
+ *   + 10/05/21: Add support for split tree output (A. Matta)                *
  *                                                                           *
  *                                                                           *
  *****************************************************************************/
@@ -31,10 +32,10 @@ using namespace std;
 
 RootOutput* RootOutput::instance = 0;
 ////////////////////////////////////////////////////////////////////////////////
-RootOutput* RootOutput::getInstance(std::string fileNameBase, std::string treeNameBase){
+RootOutput* RootOutput::getInstance(std::string fileNameBase, std::string treeNameBase,bool split){
   // A new instance of RootOutput is created if it does not exist:
   if (instance == 0) {
-    instance = new RootOutput(fileNameBase.c_str(), treeNameBase.c_str());
+    instance = new RootOutput(fileNameBase.c_str(), treeNameBase.c_str(),split);
   }
 
   // The instance of RootOutput is returned:
@@ -50,103 +51,160 @@ void RootOutput::Destroy(){
 }
 
 ////////////////////////////////////////////////////////////////////////////////
-RootOutput::RootOutput(std::string fileNameBase, std::string treeNameBase){
-  TDirectory* currentPath= gDirectory;
-
+RootOutput::RootOutput(std::string fileNameBase, std::string treeNameBase,bool split){
+
+  pSplit = (split || NPOptionManager::getInstance()->IsSplit());
+  cout << endl << "/////////// ROOT Output files ///////////" << endl;
+  cout << "Initializing ouput trees and files ";
+  if(pSplit)
+    cout << "in split mode (one tree per detector)" << endl;
+  else
+    cout << endl;
+  pTreeName=treeNameBase;
+  pCurrentDirectory= gDirectory;
   bool analysis=false;
   bool simulation=false;
-  if(fileNameBase.find("Analysis/")!=std::string::npos){
+
+  if(NPOptionManager::getInstance()->IsAnalysis()){
     analysis = true;
-    fileNameBase.erase(0,8);
   }
-  else if(fileNameBase.find("Simulation/")!=std::string::npos){
+  else if(NPOptionManager::getInstance()->IsAnalysis()){
     simulation= true;
-    fileNameBase.erase(0,10);
   }
 
-  // The ROOT file is created
-  if(!NPOptionManager::getInstance()->GetPROOF()){
-    std::string fileName;
-    if(analysis)
-      fileName = NPOptionManager::getInstance()->GetAnalysisOutputPath();
-    else if(simulation)
-      fileName = NPOptionManager::getInstance()->GetSimulationOutputPath();
-    else
-      fileName="./";
-
-    if (fileNameBase.find("root")!=std::string::npos) fileName += fileNameBase;
-    else fileName += fileNameBase + ".root";
-
-    pRootFile = new TFile(fileName.c_str(), "RECREATE");
-
-    if(treeNameBase=="SimulatedTree"){
-      string path = getenv("NPTOOL");
-        path+="/.last_sim_file";
-      ofstream last_sim_file(path.c_str());
-      last_sim_file << "TTreeName" << endl 
-       << "  " << treeNameBase <<endl
-       << "RootFileName" <<endl
-       << "  " << fileName<<endl;
-      last_sim_file.close();
-    }
-
-    else if(treeNameBase=="PhysicsTree"){
-      string path = getenv("NPTOOL");
-        path+="/.last_phy_file";
-      ofstream last_phy_file(path.c_str());
-      last_phy_file << "TTreeName" << endl 
-       << "  " << treeNameBase <<endl
-       << "RootFileName" <<endl
-       << "  " << fileName<<endl;
-      last_phy_file.close();
-    }
-
-
-    else if(treeNameBase=="ResultTree"){
-      string path = getenv("NPTOOL");
-      path+="/.last_res_file";
-      ofstream last_res_file(path.c_str());
-      last_res_file << "TTreeName" << endl 
-       << "  " << treeNameBase <<endl
-       << "RootFileName" <<endl
-       << "  " << fileName<<endl;
-      last_res_file.close();
-    }
-  
-    else{
-      string path = getenv("NPTOOL");
-        path+="/.last_any_file";
-      ofstream last_any_file(path.c_str());
-      last_any_file << "TTreeName" << endl 
-       << "  " << treeNameBase <<endl
-       << "RootFileName" <<endl
-       << "  " << fileName << endl;
-      last_any_file.close();
+  // Setup the base name
+  if(analysis)
+    pBaseName = NPOptionManager::getInstance()->GetAnalysisOutputPath();
+  else if(simulation)
+    pBaseName = NPOptionManager::getInstance()->GetSimulationOutputPath();
+  else
+    pBaseName="./";
+
+  pBaseName += "/"+fileNameBase;
+
+  if (fileNameBase.find("root")==std::string::npos) 
+    pBaseName += ".root";
+
+  if(pSplit){
+    // Create a folder for all the trees
+    string stripname = pBaseName;
+    stripname.erase(stripname.find(".root"),5);
+    string path = stripname.substr(0,stripname.rfind("/"));
+    string filebase = stripname.substr(stripname.rfind("/")+1);
+    string command = "mkdir -p "+stripname;
+    int res = system(command.c_str());
+    if(res!=0){
+      std::cout << "Error creating folder " << stripname << std::endl;
+      exit(1);
     }
+    // create the master file 
+    pMasterFile=stripname+"/"+filebase+".tree";
+    ofstream master(pMasterFile.c_str(),std::ofstream::trunc);
+    master.close();
+  }
+  else{
+    CreateTreeAndFile("global");
+  }
+  /////
+  // Create the last file 
+  if(treeNameBase=="SimulatedTree"){
+    string path = getenv("NPTOOL");
+    path+="/.last_sim_file";
+    ofstream last_sim_file(path.c_str());
+    last_sim_file << "Tree "<< pTreeName <<endl
+      << " " << pBaseName<<endl;
+    last_sim_file.close();
+  }
 
+  else if(treeNameBase=="PhysicsTree"){
+    string path = getenv("NPTOOL");
+    path+="/.last_phy_file";
+    ofstream last_phy_file(path.c_str());
+    last_phy_file << "Tree "<< pTreeName<<endl
+      << " " << pBaseName <<endl;
+    last_phy_file.close();
+  }
 
-
+  else if(treeNameBase=="ResultTree"){
+    string path = getenv("NPTOOL");
+    path+="/.last_res_file";
+    ofstream last_res_file(path.c_str());
+    last_res_file << "Tree " << pTreeName << endl 
+      << " " << pBaseName<<endl;
+    last_res_file.close();
   }
 
-  else{ // the file path must be the current directory
-    // Does not create the Output file at instantiation
-    pRootFile = 0 ;
+  else{
+    string path = getenv("NPTOOL");
+    path+="/.last_any_file";
+    ofstream last_any_file(path.c_str());
+    last_any_file << "Tree " << pTreeName <<endl
+      << " " << pBaseName<< endl;
+    last_any_file.close();
   }
 
-  pRootTree = new TTree(treeNameBase.c_str(), "Data created / analysed with the NPTool package");
-  pRootFile->SetCompressionLevel(1);
-  pRootList = new TList();
 
-  // Init TAsciiFile objects
   InitAsciiFiles();
-  gDirectory->cd(currentPath->GetPath()); 
-
-  if(NPOptionManager::getInstance()->GetCircularTree()){
-    cout << "Information: Output tree is set to circular mode" << endl;
-    pRootTree->SetCircular(1000); 
-  }
 }
 
+////////////////////////////////////////////////////////////////////////////////
+void RootOutput::CreateTreeAndFile(std::string name){
+  // Create the tree only if does not exist already
+  string file_name=pBaseName;
+
+  if(pRootFiles.find(name)==pRootFiles.end()){
+    if(pSplit){
+      string  strip= pBaseName.substr(pBaseName.rfind("/"));
+      strip = strip.substr(0,strip.rfind(".root"));
+      string  insertion= "_"+name;
+      file_name.insert(file_name.rfind(".root"),insertion);
+      file_name.insert(file_name.rfind("/"),strip);
+    }
+    cout << " - Creating output file " << file_name.c_str() << endl;
+    pRootFiles[name] = new TFile(file_name.c_str(), "RECREATE");
+    pRootTrees[name] = new TTree(pTreeName.c_str(), "Data created / analysed with the nptool package");
+    pRootFiles[name]->SetCompressionLevel(1);
+
+    // Init TAsciiFile objects
+    gDirectory->cd(pCurrentDirectory->GetPath()); 
+
+    if(NPOptionManager::getInstance()->GetCircularTree()){
+      cout << "Information: Output tree is set to circular mode" << endl;
+      pRootTrees[name]->SetCircular(1000); 
+    }
+
+    if(pSplit){// Add the tree to the .tree file
+      ofstream master(pMasterFile.c_str(),std::ofstream::app);
+      file_name = file_name.substr(file_name.rfind("/")+1);
+      master << file_name.c_str() << endl;
+      master.close();
+    }
+  } 
+  //
+  /*
+     for(auto it = m_DetectorMap.begin() ; it!=m_DetectorMap.end() ;++it){
+     string insertion = "_"+it->first;
+     master << filebase << insertion << ".root" << std::endl;
+     string filename=path+"/"+filebase+"/"+filebase+insertion+".root";
+     auto file = new TFile(filename.c_str(),"RECREATE");
+     string treename = "RawTree_"+it->first;
+     auto tree = new TTree("RawTree",treename.c_str());
+     m_TreeMap[it->first]=tree;
+     m_FileMap[it->first]=file;
+     tree->SetDirectory(file);
+     std::cout << "Splitting tree: " << filename << std::endl;
+     it->second->InitBranch(tree);
+     }
+     master.close();
+
+*/
+}
+////////////////////////////////////////////////////////////////////////////////
+void RootOutput::Fill(){
+ for(auto it = pRootTrees.begin();it!=pRootTrees.end();it++){
+  it->second->Fill();
+ } 
+}
 
 ////////////////////////////////////////////////////////////////////////////////
 void RootOutput::InitAsciiFiles(){
@@ -159,15 +217,15 @@ void RootOutput::InitAsciiFiles(){
   pEventGenerator = new TAsciiFile();
   pEventGenerator->SetNameTitle("EventGenerator", fileNameEG.c_str());
   pEventGenerator->Append(fileNameEG.c_str());
-  pEventGenerator->Write(0,TAsciiFile::kOverwrite);
-  
+  //pEventGenerator->Write(0,TAsciiFile::kOverwrite);
+
   // Detector configuration 
   // Get file name from NPOptionManager
   std::string fileNameDC = OptionManager->GetDetectorFile();
   pDetectorConfiguration = new TAsciiFile();
   pDetectorConfiguration->SetNameTitle("DetectorConfiguration", fileNameDC.c_str());
   pDetectorConfiguration->Append(fileNameDC.c_str());
-  pDetectorConfiguration->Write(0,TAsciiFile::kOverwrite);
+  //pDetectorConfiguration->Write(0,TAsciiFile::kOverwrite);
 
   // Run to treat file
   // Get file name from NPOptionManager
@@ -176,7 +234,7 @@ void RootOutput::InitAsciiFiles(){
     std::string fileNameRT = OptionManager->GetRunToReadFile();
     pRunToTreatFile->SetNameTitle("RunToTreat", fileNameRT.c_str());
     pRunToTreatFile->Append(fileNameRT.c_str());
-    pRunToTreatFile->Write(0,TAsciiFile::kOverwrite);
+    //pRunToTreatFile->Write(0,TAsciiFile::kOverwrite);
   }
 
   // Calibration files
@@ -184,61 +242,44 @@ void RootOutput::InitAsciiFiles(){
   if (!OptionManager->IsDefault("Calibration")) {
     std::string fileNameCal = OptionManager->GetCalibrationFile();
     pCalibrationFile->SetNameTitle("Calibration", fileNameCal.c_str());
-    pCalibrationFile->Write(0,TAsciiFile::kOverwrite);
+    //pCalibrationFile->Write(0,TAsciiFile::kOverwrite);
   }
 
   // Analysis configuration files
   pAnalysisConfigFile = new TAsciiFile();
   pAnalysisConfigFile->SetNameTitle("AnalysisConfig", "AnalysisConfig");
-  pAnalysisConfigFile->Write(0,TAsciiFile::kOverwrite);
+  //pAnalysisConfigFile->Write(0,TAsciiFile::kOverwrite);
 }
 
-
-
 ////////////////////////////////////////////////////////////////////////////////
 RootOutput::~RootOutput(){ 
   // The data is written to the file and the tree is closed:
-  if (pRootFile && !NPOptionManager::getInstance()->GetPROOF()) {
-    TDirectory* currentPath= gDirectory;
-    gDirectory->cd(pRootFile->GetPath());
-   
-    // write TAsciiFile if used
-    // EventGenerator
-    if (!pEventGenerator->IsEmpty()) pEventGenerator->Write(0,TAsciiFile::kOverwrite);
-    // DetectorConfiguration
-    if (!pDetectorConfiguration->IsEmpty()) pDetectorConfiguration->Write(0,TAsciiFile::kOverwrite);
-    // CalibrationFile
-    if (!pCalibrationFile->IsEmpty()) pCalibrationFile->Write(0,TAsciiFile::kOverwrite);
-    // RunToTreatFile
-    if (!pRunToTreatFile->IsEmpty()) pRunToTreatFile->Write(0,TAsciiFile::kOverwrite);
-    // Analysis ConfigFile
-    if (!pAnalysisConfigFile->IsEmpty()) pAnalysisConfigFile->Write(0,TAsciiFile::kOverwrite);
-   
-    cout << endl;
-    cout << endl << "Root Output summary" << endl;
-    cout << "  - Number of entries in the Tree: " << pRootTree->GetEntries() << endl;
-    cout << "  - Number of bites written to file: " << pRootTree->Write(0, TObject::kOverwrite) << endl;
-    pRootFile->Flush();
-    pRootFile->Purge(1);
-
-    gDirectory->cd(currentPath->GetPath());
-    pRootFile->Close();
-  }
-
-  else if (pRootFile && NPOptionManager::getInstance()->GetPROOF()){
-    if (!pEventGenerator->IsEmpty()) pEventGenerator->Write(0,TAsciiFile::kOverwrite);
-    // DetectorConfiguration
-    if (!pDetectorConfiguration->IsEmpty()) pDetectorConfiguration->Write(0,TAsciiFile::kOverwrite);
-    // CalibrationFile
-    if (!pCalibrationFile->IsEmpty()) pCalibrationFile->Write(0,TAsciiFile::kOverwrite);
-    // RunToTreatFile
-    if (!pRunToTreatFile->IsEmpty()) pRunToTreatFile->Write(0,TAsciiFile::kOverwrite);
-    // Analysis ConfigFile
-    if (!pAnalysisConfigFile->IsEmpty()) pAnalysisConfigFile->Write(0,TAsciiFile::kOverwrite);
-  }
-
-  else if(!pRootFile && NPOptionManager::getInstance()->GetPROOF()){
-
+  if (pRootFiles.size()>0) {
+    cout << endl << endl << "Root Output summary" << endl;
+    TDirectory* pCurrentDirectory= gDirectory;
+    for(auto it = pRootFiles.begin(); it!=pRootFiles.end();it++){
+      cout << " - " <<it->first << " tree and file " << endl;
+      gDirectory->cd(it->second->GetPath());
+      // write TAsciiFile if used
+      // EventGenerator
+      if (!pEventGenerator->IsEmpty()) pEventGenerator->Write(0,TAsciiFile::kOverwrite);
+      // DetectorConfiguration
+      if (!pDetectorConfiguration->IsEmpty()) pDetectorConfiguration->Write(0,TAsciiFile::kOverwrite);
+      // CalibrationFile
+      if (!pCalibrationFile->IsEmpty()) pCalibrationFile->Write(0,TAsciiFile::kOverwrite);
+      // RunToTreatFile
+      if (!pRunToTreatFile->IsEmpty()) pRunToTreatFile->Write(0,TAsciiFile::kOverwrite);
+      // Analysis ConfigFile
+      if (!pAnalysisConfigFile->IsEmpty()) pAnalysisConfigFile->Write(0,TAsciiFile::kOverwrite);
+
+      cout << "  -> Number of entries in the " << it->first << " Tree: " << pRootTrees[it->first]->GetEntries() << endl;
+      cout << "  -> Number of bites written to file: " << pRootTrees[it->first]->Write(0, TObject::kOverwrite) << endl;
+      it->second->Flush();
+      it->second->Purge(1);
+
+      gDirectory->cd(pCurrentDirectory->GetPath());
+      it->second->Close();
+    }
   }
 
   else {
@@ -247,20 +288,27 @@ RootOutput::~RootOutput(){
 }
 
 ////////////////////////////////////////////////////////////////////////////////
-TFile* RootOutput::InitFile(std::string fileNameBase){
-
-  if(NPOptionManager::getInstance()->GetPROOF()){
-    std::string GlobalPath = getenv("NPTOOL");
-    std::string fileName = GlobalPath + "/Outputs/Analysis/";
-    if (fileNameBase.find("root")!=std::string::npos) fileName += fileNameBase;
-    else fileName += fileNameBase + ".root";
-    pRootFile = new TFile(fileName.c_str(), "RECREATE");
-    pRootFile->Flush();
-    return pRootFile;
-  }
+TFile* RootOutput::GetFile(std::string name)  {
+  if(!pSplit)
+    name="global";
 
+  if(pRootFiles.find(name)!=pRootFiles.end())
+    return pRootFiles[name];
   else{
-    cout << "ERROR: Do not use RootOutput::InitFile without a proof environment (use --proof option to NPTool)" << endl ;
+    std::cout << "Error: Requested file for detector " << name << " does not exist" << std::endl;
     exit(1);
   }
+  return 0;
 }
+
+////////////////////////////////////////////////////////////////////////////////
+TTree* RootOutput::GetTree(std::string name) {
+  if(!pSplit)
+    name="global";
+
+  if(pRootTrees.find(name)==pRootTrees.end())
+    CreateTreeAndFile(name);
+  
+  return pRootTrees[name];
+}
+
diff --git a/NPLib/Core/RootOutput.h b/NPLib/Core/RootOutput.h
index 2ac1db6018403851367dd5f29259c46c369bbc9f..d2605aa0fbad743a6e7a3f398dd287de8f916a85 100644
--- a/NPLib/Core/RootOutput.h
+++ b/NPLib/Core/RootOutput.h
@@ -1,7 +1,7 @@
 #ifndef ROOTOUTPUT_HH
 #define ROOTOUTPUT_HH
 /*****************************************************************************
- * Copyright (C) 2009-2016    this file is part of the NPTool Project        *
+ * Copyright (C) 2009-2021   this file is part of the NPTool Project         *
  *                                                                           *
  * For the licensing terms see $NPTOOL/Licence/NPTool_Licence                *
  * For the list of contributors see $NPTOOL/Licence/Contributors             *
@@ -11,18 +11,20 @@
  * Original Author: N. de Sereville  contact address: deserevi@ipno.in2p3.fr *
  *                                                                           *
  * Creation Date  : 21/07/09                                                 *
- * Last update    : 03/02/11                                                 *
+ * Last update    : 10/05/21                                                 *
  *---------------------------------------------------------------------------*
  * Decription: This class is a singleton class which deals with the ROOT     *
  *             output file and tree both for NPSimulation and NPAnalysis.    *
  *---------------------------------------------------------------------------*
  * Comment:                                                                  *
  *   + 03/02/11: Add support for TAsciiFile objects (N. de Sereville)        *
+ *   + 10/05/21: Add support for split tree output (A. Matta)                *
  *                                                                           *
  *                                                                           *
  *****************************************************************************/
 
-
+// STL
+#include <map>
 
 // NPL headers
 #include "TAsciiFile.h"
@@ -41,7 +43,8 @@ public:
   // it does not yet exist:
   // (see the constructor for an explanation of the arguments)
   static RootOutput* getInstance(std::string fileNameBase = "Simulation",
-                                 std::string treeNameBase = "SimulatedTree");
+                                 std::string treeNameBase = "SimulatedTree",
+                                 bool split = false);
   
   // The analysis class instance can be deleted by calling the Destroy
   // method (NOTE: The class destructor is protected, and can thus not be
@@ -50,8 +53,8 @@ public:
   
 protected:
   // Constructor (protected)
-  RootOutput(std::string fileNameBase, std::string treeNameBase);
-  
+  RootOutput(std::string fileNameBase, std::string treeNameBase, bool split=false);
+ 
   // Destructor (protected)
   virtual ~RootOutput();
   
@@ -65,27 +68,32 @@ private:
   
 private:
   void InitAsciiFiles();
-  
+  void CreateTreeAndFile(std::string name);
 public:
-  TFile*      GetFile()                           {return pRootFile;}
-  TTree*      GetTree()                           {return pRootTree;}
-  TList*      GetList()                           {return pRootList;}
+  TFile*      GetFile(std::string name="global") ; 
+  TTree*      GetTree(std::string name="global") ;
+  TList*      GetList(std::string name="global") ;
   TAsciiFile* GetAsciiFileEventGenerator()        {return pEventGenerator;}
   TAsciiFile* GetAsciiFileDetectorConfiguration() {return pDetectorConfiguration;}
   TAsciiFile* GetAsciiFileCalibration()           {return pCalibrationFile;}
   TAsciiFile* GetAsciiFileRunToTreat()            {return pRunToTreatFile;}
   TAsciiFile* GetAsciiFileAnalysisConfig()        {return pAnalysisConfigFile;}
-  TFile*      InitFile(std::string fileNameBase); // use only for proof environment
+  void        Fill();
   
 private:
-  TFile      *pRootFile;
-  TTree      *pRootTree;
-  TList      *pRootList;
-  TAsciiFile *pEventGenerator;
-  TAsciiFile *pDetectorConfiguration;
-  TAsciiFile *pCalibrationFile;
-  TAsciiFile *pRunToTreatFile;
-  TAsciiFile *pAnalysisConfigFile;
+  std::string                   pBaseName;
+  std::string                   pTreeName;
+  std::string                   pMasterFile;
+  TDirectory*                   pCurrentDirectory;
+  bool                          pSplit;
+  std::map<std::string, TFile*> pRootFiles;
+  std::map<std::string, TTree*> pRootTrees;
+  std::map<std::string, TList*> pRootLists;
+  TAsciiFile* pEventGenerator;
+  TAsciiFile* pDetectorConfiguration;
+  TAsciiFile* pCalibrationFile;
+  TAsciiFile* pRunToTreatFile;
+  TAsciiFile* pAnalysisConfigFile;
   
 };
 
diff --git a/NPLib/Detectors/Minos/TMinosPhysics.cxx b/NPLib/Detectors/Minos/TMinosPhysics.cxx
index d23cedc5d7225810c8a514001da404948095153d..3fa3a22a15403b09f0c2192de9d41a075ae99c4a 100644
--- a/NPLib/Detectors/Minos/TMinosPhysics.cxx
+++ b/NPLib/Detectors/Minos/TMinosPhysics.cxx
@@ -306,7 +306,7 @@ void TMinosPhysics::InitializeRootInputPhysics() {
 
 ///////////////////////////////////////////////////////////////////////////
 void TMinosPhysics::InitializeRootOutput() {
-  TTree* outputTree = RootOutput::getInstance()->GetTree();
+  TTree* outputTree = RootOutput::getInstance()->GetTree("Minos");
   outputTree->Branch("Minos", "TMinosPhysics", &m_EventPhysics);
 }
 
diff --git a/NPLib/Detectors/Nebula/TNebulaPhysics.cxx b/NPLib/Detectors/Nebula/TNebulaPhysics.cxx
index 0035728f376a6bf045e7c837a2b0dc8789128dbc..7a74ce45addf621cee2e0907f0ca094fb5c1d14a 100644
--- a/NPLib/Detectors/Nebula/TNebulaPhysics.cxx
+++ b/NPLib/Detectors/Nebula/TNebulaPhysics.cxx
@@ -313,7 +313,7 @@ void TNebulaPhysics::InitializeRootInputPhysics() {
 
 ///////////////////////////////////////////////////////////////////////////
 void TNebulaPhysics::InitializeRootOutput() {
-  TTree* outputTree = RootOutput::getInstance()->GetTree();
+  TTree* outputTree = RootOutput::getInstance()->GetTree("Nebula");
   outputTree->Branch("Nebula", "TNebulaPhysics", &m_EventPhysics);
 }
 
diff --git a/NPLib/Detectors/Samurai/CMakeLists.txt b/NPLib/Detectors/Samurai/CMakeLists.txt
index 018365593859b6f7c9e87e6cd8cf76f8c381dca7..639fd771a6cacc88c3dfd0f1bcc7e981089be24f 100644
--- a/NPLib/Detectors/Samurai/CMakeLists.txt
+++ b/NPLib/Detectors/Samurai/CMakeLists.txt
@@ -2,7 +2,9 @@ add_custom_command(OUTPUT TSamuraiHodoscopeDataDict.cxx COMMAND ${CMAKE_BINARY_D
 
 add_custom_command(OUTPUT TSamuraiHodoscopePhysicsDict.cxx COMMAND ${CMAKE_BINARY_DIR}/scripts/build_dict.sh TSamuraiHodoscopePhysics.h TSamuraiHodoscopePhysicsDict.cxx TSamuraiHodoscopePhysics.rootmap libNPSamurai.dylib DEPENDS TSamuraiHodoscopePhysics.h)
 
+add_custom_command(OUTPUT TSamuraiBDCDataDict.cxx COMMAND ${CMAKE_BINARY_DIR}/scripts/build_dict.sh TSamuraiBDCData.h TSamuraiBDCDataDict.cxx TSamuraiBDCData.rootmap libNPSamurai.dylib DEPENDS TSamuraiBDCData.h)
 
+add_custom_command(OUTPUT TSamuraiBDCPhysicsDict.cxx COMMAND ${CMAKE_BINARY_DIR}/scripts/build_dict.sh TSamuraiBDCPhysics.h TSamuraiBDCPhysicsDict.cxx TSamuraiBDCPhysics.rootmap libNPSamurai.dylib DEPENDS TSamuraiBDCPhysics.h)
 
 add_custom_command(OUTPUT TSamuraiFDC2DataDict.cxx COMMAND ${CMAKE_BINARY_DIR}/scripts/build_dict.sh TSamuraiFDC2Data.h TSamuraiFDC2DataDict.cxx TSamuraiFDC2Data.rootmap libNPSamurai.dylib DEPENDS TSamuraiFDC2Data.h)
 
@@ -12,8 +14,8 @@ add_custom_command(OUTPUT TSamuraiFDC0DataDict.cxx COMMAND ${CMAKE_BINARY_DIR}/s
 
 add_custom_command(OUTPUT TSamuraiFDC0PhysicsDict.cxx COMMAND ${CMAKE_BINARY_DIR}/scripts/build_dict.sh TSamuraiFDC0Physics.h TSamuraiFDC0PhysicsDict.cxx TSamuraiFDC0Physics.rootmap libNPSamurai.dylib DEPENDS TSamuraiFDC0Physics.h)
 
-add_library(NPSamurai SHARED TSamuraiHodoscopeData.cxx TSamuraiFDC2Data.cxx TSamuraiFDC2DataDict.cxx TSamuraiHodoscopePhysics.cxx TSamuraiHodoscopePhysicsDict.cxx TSamuraiFDC2Physics.cxx TSamuraiHodoscopeDataDict.cxx TSamuraiFDC2PhysicsDict.cxx TSamuraiFDC0Data.cxx TSamuraiFDC0DataDict.cxx TSamuraiFDC0Physics.cxx TSamuraiFDC0PhysicsDict.cxx)
+add_library(NPSamurai SHARED TSamuraiHodoscopeData.cxx TSamuraiBDCData.cxx TSamuraiBDCPhysics.cxx TSamuraiBDCDataDict.cxx TSamuraiBDCPhysicsDict.cxx TSamuraiFDC2Data.cxx TSamuraiFDC2DataDict.cxx TSamuraiHodoscopePhysics.cxx TSamuraiHodoscopePhysicsDict.cxx TSamuraiFDC2Physics.cxx TSamuraiHodoscopeDataDict.cxx TSamuraiFDC2PhysicsDict.cxx TSamuraiFDC0Data.cxx TSamuraiFDC0DataDict.cxx TSamuraiFDC0Physics.cxx TSamuraiFDC0PhysicsDict.cxx)
 
 target_link_libraries(NPSamurai ${ROOT_LIBRARIES} NPCore NPTrackReconstruction) 
-install(FILES TSamuraiHodoscopeData.h TSamuraiHodoscopePhysics.h TSamuraiFDC2Data.h TSamuraiFDC2Physics.h  TSamuraiFDC0Data.h TSamuraiFDC0Physics.h SamuraiDCIndex.h DESTINATION ${CMAKE_INCLUDE_OUTPUT_DIRECTORY})
+install(FILES TSamuraiBDCData.h TSamuraiBDCPhysics.h TSamuraiHodoscopeData.h TSamuraiHodoscopePhysics.h TSamuraiFDC2Data.h TSamuraiFDC2Physics.h  TSamuraiFDC0Data.h TSamuraiFDC0Physics.h SamuraiDCIndex.h DESTINATION ${CMAKE_INCLUDE_OUTPUT_DIRECTORY})
 
diff --git a/NPLib/Detectors/Samurai/TSamuraiBDCData.cxx b/NPLib/Detectors/Samurai/TSamuraiBDCData.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..862fdde7635fbf08bdda51f0670f0ca0d8cb1ac4
--- /dev/null
+++ b/NPLib/Detectors/Samurai/TSamuraiBDCData.cxx
@@ -0,0 +1,59 @@
+#include "TSamuraiBDCData.h"
+#include <iostream>
+
+TSamuraiBDCData::TSamuraiBDCData(){};
+TSamuraiBDCData::~TSamuraiBDCData(){};
+
+////////////////////////////////////////////////////////////////////////////////
+void TSamuraiBDCData::SetData(const int& Det, const int& Layer, const int& Wire, const double& Time, const int& Edge){
+  fBDC_DetectorNbr.push_back(Det);
+  fBDC_LayerNbr.push_back(Layer); 
+  fBDC_WireNbr.push_back(Wire); 
+  fBDC_Time.push_back(Time); 
+  fBDC_Edge.push_back(Edge); 
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void TSamuraiBDCData::Clear(){
+  fBDC_DetectorNbr.clear();
+  fBDC_LayerNbr.clear();
+  fBDC_WireNbr.clear();
+  fBDC_Time.clear();
+  fBDC_Edge.clear();
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void TSamuraiBDCData::Print(){
+  using namespace std;
+
+  cout << " -- Event:" << endl;
+  cout << "   - Multiplicity: " << Mult() << endl;
+
+}
+////////////////////////////////////////////////////////////////////////////////
+unsigned int TSamuraiBDCData::MultLayer(unsigned int det , unsigned int layer, int edge){
+  unsigned int mult=0;
+  unsigned int size = fBDC_DetectorNbr.size();
+  for(unsigned int i = 0 ; i< size ; i++ ){
+    if(fBDC_DetectorNbr[i]==det)
+      if(fBDC_LayerNbr[i]==layer)
+        if(fBDC_Edge[i]==edge && edge!=-1) // edge type is specified (0 or 1)
+          mult++;
+        else if(edge==-1)// edge type is not specified
+          mult++;
+  }
+  return mult;
+
+}
+////////////////////////////////////////////////////////////////////////////////
+std::vector<int> TSamuraiBDCData::GetWire(unsigned int det , unsigned int layer){
+  std::vector<int> wires;
+  unsigned int size = fBDC_DetectorNbr.size();
+  for(unsigned int i = 0 ; i< size ; i++ ){
+    if(fBDC_DetectorNbr[i]==det)
+      if(fBDC_LayerNbr[i]==layer)
+        wires.push_back(fBDC_WireNbr[i]);
+  }
+  return wires;
+}
+ClassImp(TSamuraiBDCData); 
diff --git a/NPLib/Detectors/Samurai/TSamuraiBDCData.h b/NPLib/Detectors/Samurai/TSamuraiBDCData.h
new file mode 100644
index 0000000000000000000000000000000000000000..50e87d81d3b3fdaeea34c9263216b09a58577a61
--- /dev/null
+++ b/NPLib/Detectors/Samurai/TSamuraiBDCData.h
@@ -0,0 +1,38 @@
+#ifndef TSamuraiBDCData_H
+#define TSamuraiBDCData_H
+#include "TObject.h"
+#include <vector>
+
+class TSamuraiBDCData: public TObject{
+  public:
+    TSamuraiBDCData();
+    ~TSamuraiBDCData();
+
+  private:
+    std::vector<int> fBDC_DetectorNbr;
+    std::vector<int> fBDC_LayerNbr;
+    std::vector<int> fBDC_WireNbr;
+    std::vector<double> fBDC_Time;
+    std::vector<int> fBDC_Edge;
+  
+  public:
+    void Clear();
+    void Print();
+    void Clear(const Option_t*) {};
+    void Dump() const{};
+  
+  public:
+    void SetData(const int& Det, const int& Layer, const int& Wire, const double& Time, const int& Edge);
+    unsigned int Mult(){return fBDC_DetectorNbr.size();};
+    unsigned int MultLayer(unsigned int det , unsigned int layer, int edge=-1);
+    std::vector<int> GetWire(unsigned int det , unsigned int layer);
+    int const GetDetectorNbr(const unsigned int& i){return fBDC_DetectorNbr[i];};
+    int const GetLayerNbr(const unsigned int& i){return fBDC_LayerNbr[i];};
+    int const GetWireNbr(const unsigned int& i){return fBDC_WireNbr[i];};
+    double const GetTime(const unsigned int& i){return fBDC_Time[i];};
+    int const GetEdge(const unsigned int& i){return fBDC_Edge[i];};
+
+    ClassDef(TSamuraiBDCData,1); 
+};
+
+#endif
diff --git a/NPLib/Detectors/Samurai/TSamuraiBDCPhysics.cxx b/NPLib/Detectors/Samurai/TSamuraiBDCPhysics.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..6a73e365eda3e9b1b8f7ebf37d22882cdc87daba
--- /dev/null
+++ b/NPLib/Detectors/Samurai/TSamuraiBDCPhysics.cxx
@@ -0,0 +1,479 @@
+/*****************************************************************************
+ * Copyright (C) 2009-2020   this file is part of the NPTool Project         *
+ *                                                                           *
+ * For the licensing terms see $NPTOOL/Licence/NPTool_Licence                *
+ * For the list of contributors see $NPTOOL/Licence/Contributors             *
+ *****************************************************************************/
+
+/*****************************************************************************
+ * Original Author: Adrien MATTA  contact address: matta@lpccaen.in2p3.fr    *
+ *                                                                           *
+ * Creation Date  : May 2021                                                 *
+ * Last update    :                                                          *
+ *---------------------------------------------------------------------------*
+ * Decription:                                                               *
+ *  This class hold SamuraiBDC treated data                                  *
+ *                                                                           *
+ *---------------------------------------------------------------------------*
+ * Comment:                                                                  *
+ *                                                                           *
+ *****************************************************************************/
+#include "TSamuraiBDCPhysics.h"
+
+//   STL
+#include <sstream>
+#include <iostream>
+#include <cmath>
+#include <stdlib.h>
+#include <limits>
+using namespace std;
+
+//   NPL
+#include "RootInput.h"
+#include "RootOutput.h"
+#include "TAsciiFile.h"
+#include "NPOptionManager.h"
+#include "NPDetectorFactory.h"
+#include "NPSystemOfUnits.h"
+//   ROOT
+using namespace NPUNITS;
+///////////////////////////////////////////////////////////////////////////
+
+ClassImp(TSamuraiBDCPhysics)
+  ///////////////////////////////////////////////////////////////////////////
+  TSamuraiBDCPhysics::TSamuraiBDCPhysics(){
+    m_EventData         = new TSamuraiBDCData ;
+    m_PreTreatedData    = new TSamuraiBDCData ;
+    m_EventPhysics      = this ;
+    //m_Spectra           = NULL;
+    ToTThreshold_L = 0;
+    ToTThreshold_H = 180;
+    DriftLowThreshold=0.1 ;
+    DriftUpThreshold=2.4;
+    PowerThreshold=5;
+  }
+
+///////////////////////////////////////////////////////////////////////////
+void TSamuraiBDCPhysics::BuildSimplePhysicalEvent(){
+  BuildPhysicalEvent();
+}
+
+///////////////////////////////////////////////////////////////////////////
+void TSamuraiBDCPhysics::BuildPhysicalEvent(){
+  PreTreat();
+//  RemoveNoise();
+
+  // Map[plane angle, vector of spatial information]
+  static map<double, vector<double> > X ; 
+  static map<double, vector<double> > Z ; 
+  static map<double, vector<double> > R ; 
+  static int det,layer,wire;
+  X.clear();Z.clear();R.clear();
+  unsigned int size = Detector.size();
+  for(unsigned int i = 0 ; i < size ; i++){
+    if(DriftLength[i] > DriftLowThreshold && DriftLength[i] < DriftUpThreshold){
+      det = Detector[i];
+      layer = Layer[i];
+      wire = Wire[i]; 
+      SamuraiDCIndex idx(det,layer,wire);
+      X[Wire_Angle[idx]].push_back(Wire_X[idx]); 
+      Z[Wire_Angle[idx]].push_back(Wire_Z[idx]); 
+      R[Wire_Angle[idx]].push_back(DriftLength[i]); 
+    }
+  }
+  
+  // Reconstruct the vector for each of the plane of each of the detector
+  static double X0,X100,a,b; // store the BuildTrack2D results
+  static map<double, TVector3 > VX0 ;  
+  static map<double, TVector3 > VX100 ;  
+  static map<double, double > D ;// the minimum distance  
+  static unsigned int uid; uid=0;
+  VX0.clear();VX100.clear(),D.clear();
+
+  for(auto it = X.begin();it!=X.end();++it){
+#if __cplusplus > 199711L && NPMULTITHREADING 
+    m_reconstruction.AddPlan(uid++,X[it->first],Z[it->first],R[it->first]); 
+#else
+    D[it->first]=m_reconstruction.BuildTrack2D(X[it->first],Z[it->first],R[it->first],X0,X100,a,b); 
+#endif 
+    }
+#if __cplusplus > 199711L && NPMULTITHREADING 
+  // do all plan at once in parallele, return when all plan are done
+  m_reconstruction.BuildTrack2D();
+  uid=0;
+#endif
+  
+  for(auto it = X.begin();it!=X.end();++it){
+#if __cplusplus > 199711L && NPMULTITHREADING 
+  D[it->first]=m_reconstruction.GetResults(uid++,X0,X100,a,b); 
+#endif
+
+   // for Debug, write a file of 
+/*   { std::ofstream f("distance.txt", std::ios::app);
+   f<< D[it->first] << endl;
+   f.close();
+   }
+  */ 
+    // very large a means track perpendicular to the chamber, what happen when there is pile up
+    if(abs(a)>5000)
+      PileUp++;
+
+    Mult+=X[it->first].size();
+    // Position at z=0
+    TVector3 P(X0,0,0);
+    P.RotateZ(it->first);
+    VX0[it->first]=P;
+    // Position at z=100
+    TVector3 P100= TVector3(X100,0,0);
+    P100.RotateZ(it->first);
+    VX100[it->first]=P100;
+
+  }
+  
+  // Reconstruct the central position (z=0) for each detector
+  static vector<TVector3> C ;  
+  static vector<double  > W ; // weight based on D  
+  C.clear(),W.clear();
+  TVector3 P;
+  for(auto it1 = VX0.begin();it1!=VX0.end();++it1){
+    for(auto it2 = it1;it2!=VX0.end();++it2){
+      if(it1!=it2){// different plane
+        //cout << "BDC" << endl;
+        m_reconstruction.ResolvePlane(it1->second,it1->first,it2->second,it2->first,P);
+        // cout << "done " << D[it1->first] << " " << D[it2->first] << endl;
+        if(P.X()!=-10000 && D[it1->first]<PowerThreshold&& D[it2->first]<PowerThreshold){
+          C.push_back(P);
+          // Mean pos are weighted based on the the sum of distance from track
+          // to hit obtained during the minimisation
+          W.push_back(1./sqrt(D[it1->first]*D[it2->first]));
+          }
+      }
+    }
+  }
+  // Reconstruct the position at z=100 for each detector
+  static vector<TVector3> C100 ;  
+  C100.clear();
+  for(auto it1 = VX100.begin();it1!=VX100.end();++it1){
+    for(auto it2 = it1;it2!=VX100.end();++it2){
+      if(it1!=it2 ){// different plane, same detector
+        m_reconstruction.ResolvePlane(it1->second,it1->first,it2->second,it2->first,P);
+
+        if(P.X()!=-10000&& D[it1->first]<PowerThreshold && D[it2->first]<PowerThreshold)
+          C100.push_back(P);
+      }
+    }
+  }
+  
+  // Build the Reference position by averaging all possible pair 
+  size = C.size();
+  static double PosX100,PosY100,norm;
+  if(size){
+    PosX=0;
+    PosY=0;
+    PosX100=0;
+    PosY100=0;
+    norm=0;
+    for(unsigned int i = 0 ; i < size ; i++){
+      PosX+= C[i].X()*W[i]; 
+      PosY+= C[i].Y()*W[i]; 
+      PosX100+= C100[i].X()*W[i]; 
+      PosY100+= C100[i].Y()*W[i]; 
+      norm+=W[i];
+    } 
+    MultMean=size;
+    // Mean position at Z=0
+    PosX=PosX/norm; 
+    PosY=PosY/norm; 
+    // Mean position at Z=100
+    PosX100=PosX100/norm; 
+    PosY100=PosY100/norm; 
+
+    devX=0;
+    devY=0;
+    for(unsigned int i = 0 ; i < size ; i++){
+      devX+=W[i]*(C[i].X()-PosX)*(C[i].X()-PosX);
+      devY+=W[i]*(C[i].Y()-PosY)*(C[i].Y()-PosY);
+    }
+    devX=sqrt(devX/((size-1)*norm));
+    devY=sqrt(devY/((size-1)*norm));
+    // Compute ThetaX, angle between the Direction vector projection in XZ with
+    // the Z axis
+    //ThetaX=atan((PosX100-PosX)/100.);
+    ThetaX = (PosX100-PosX)/100.;
+    // Compute PhiY, angle between the Direction vector projection in YZ with
+    // the Z axis
+    //PhiY=atan((PosY100-PosY)/100.);
+    PhiY=(PosY100-PosY)/100.;
+    Dir=TVector3(PosX100-PosX,PosY100-PosY,100).Unit();
+  }
+
+  return;
+}
+///////////////////////////////////////////////////////////////////////////
+void TSamuraiBDCPhysics::PreTreat(){
+  ClearPreTreatedData();
+  static CalibrationManager* Cal = CalibrationManager::getInstance();
+  static string channel;
+  // one map per detector
+  map<unsigned int, vector<double> > X ; 
+  map<unsigned int, vector<double> > Z ; 
+  map<unsigned int, vector<double> > R ; 
+
+  unsigned int size = m_EventData->Mult();
+  for(unsigned int i = 0 ; i < size ; i++){
+    // EDGE=1 is the leading edge, IE, the real time.
+    // EDGE=0 is the trailing edge, so it helps build Tot
+    if(m_EventData->GetEdge(i)==1){
+      int det   = m_EventData->GetDetectorNbr(i); 
+      int layer = m_EventData->GetLayerNbr(i); 
+      int wire  = m_EventData->GetWireNbr(i); 
+      double time = m_EventData->GetTime(i);
+      double etime = 0;
+      // look for matching trailing edge   
+      for(unsigned int j = 0 ; j < size ; j++){
+        if(m_EventData->GetEdge(j)==0){
+          int edet   = m_EventData->GetDetectorNbr(j); 
+          int elayer = m_EventData->GetLayerNbr(j); 
+          int ewire  = m_EventData->GetWireNbr(j); 
+          // same wire
+          if(wire==ewire && layer==elayer && det==edet){
+            etime = m_EventData->GetTime(j); 
+          }    
+        }
+        if(etime && etime>time)
+          break;
+        else
+          etime=0;
+      }
+      // a valid wire must have an edge
+      if(etime && time && etime-time>ToTThreshold_L && etime-time<ToTThreshold_H){
+        Detector.push_back(det);
+        Layer.push_back(layer);       
+        Wire.push_back(wire);
+        Time.push_back(time);
+        ToT.push_back(etime-time);
+        channel="SamuraiBDC/L" + NPL::itoa(layer);
+        DriftLength.push_back(2.5-Cal->ApplySigmoid(channel,etime));
+      }
+    }
+
+  }
+  return;
+}
+///////////////////////////////////////////////////////////////////////////
+void TSamuraiBDCPhysics::RemoveNoise(){
+  // Remove the noise by looking if a matching wire exist in the adjacent layer
+  // this done by looking at the closest plane with the same orientation
+  unsigned int size = Detector.size(); 
+  for(unsigned int i = 0 ; i < size ; i++){
+    bool match=false;
+    int det = Detector[i];
+    int layer = Layer[i];
+    int wire = Wire[i];
+    // look for matching adjacent wire   
+
+    for(unsigned int j = 0 ; j < size ; j++){
+      int adet = Detector[j];
+      int alayer = Layer[j];
+      int awire = Wire[j];
+      bool blayer = false;
+      if(layer%2==0){
+        if(layer+1==alayer)
+          blayer=true;
+      }
+
+      else{
+        if(layer-1==alayer)
+          blayer=true;
+      }
+
+      if(det==adet && blayer && abs(wire-awire)<=1){
+        match=true;
+        break;
+      }
+    }
+
+    if(match)
+      Matched.push_back(true);
+    else
+      Matched.push_back(false);
+  }
+  return;
+}
+///////////////////////////////////////////////////////////////////////////
+void TSamuraiBDCPhysics::Clear(){
+  MultMean=0;
+  PileUp=0;
+  Mult=0;
+  PosX=PosY=-10000;
+  devX=devY=-10000;
+  DriftLength.clear();
+  Detector.clear();
+  Layer.clear();
+  Wire.clear();
+  Time.clear();
+  ToT.clear();
+  ParticleDirection.clear();
+  MiddlePosition.clear();
+  Matched.clear();
+}
+///////////////////////////////////////////////////////////////////////////
+
+////   Innherited from VDetector Class   ////
+
+///////////////////////////////////////////////////////////////////////////
+void TSamuraiBDCPhysics::ReadConfiguration(NPL::InputParser parser){
+  vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithToken("SAMURAIBDC");
+  if(NPOptionManager::getInstance()->GetVerboseLevel())
+    cout << "//// " << blocks.size() << " detector(s) found " << endl; 
+
+  vector<string> token= {"XML"};
+
+  for(unsigned int i = 0 ; i < blocks.size() ; i++){
+    cout << endl << "////  Samurai BDC (" << i+1 << ")" << endl;
+    string xmlpath = blocks[i]->GetString("XML");
+    NPL::XmlParser xml;
+    xml.LoadFile(xmlpath);
+    AddDC("SAMURAIBDC",xml);
+  }
+
+#if __cplusplus > 199711L && NPMULTITHREADING 
+  if(blocks.size()){ // if a detector is found, init the thread pool
+    // one thread for each plan X,Y = 2
+    // ! more than that this will not help !
+    m_reconstruction.SetNumberOfThread(2);
+    m_reconstruction.InitThreadPool();
+  }
+#endif 
+}
+
+///////////////////////////////////////////////////////////////////////////
+void TSamuraiBDCPhysics::AddDC(string name, NPL::XmlParser& xml){
+  std::vector<NPL::XML::block*> b = xml.GetAllBlocksWithName(name);  
+  // BDC case
+  if(name=="SAMURAIBDC"){
+    unsigned int det=0;
+    unsigned int sizeB = b.size();
+    for(unsigned int i = 0 ; i < sizeB ; i++){
+      unsigned int layer = b[i]->AsInt("layer"); 
+      unsigned int wire  = b[i]->AsInt("wireid"); 
+      double X = b[i]->AsDouble("wirepos");  
+      double Z = b[i]->AsDouble("wirez");  
+      string sDir = b[i]->AsString("anodedir");
+      double T=0;
+      if(sDir=="X")
+        T= 0*deg;
+      else if(sDir=="Y")
+        T= -90*deg;
+      else if(sDir=="U")
+        T=-30*deg;
+      else if(sDir=="V")
+        T=+30*deg;
+      else{
+        cout << "ERROR: Unknown layer orientation for Samurai BDC"<< endl;
+        exit(1);
+      }
+      SamuraiDCIndex idx(det,layer,wire);
+      Wire_X[idx]=X;
+      Wire_Z[idx]=Z;
+      Wire_Angle[idx]=T;
+    }
+  }
+}
+
+///////////////////////////////////////////////////////////////////////////
+void TSamuraiBDCPhysics::InitSpectra(){  
+  //m_Spectra = new TSamuraiBDCSpectra(m_NumberOfDetector);
+}
+
+///////////////////////////////////////////////////////////////////////////
+void TSamuraiBDCPhysics::FillSpectra(){  
+  //  m_Spectra -> FillRawSpectra(m_EventData);
+  //  m_Spectra -> FillPreTreatedSpectra(m_PreTreatedData);
+  //  m_Spectra -> FillPhysicsSpectra(m_EventPhysics);
+}
+///////////////////////////////////////////////////////////////////////////
+void TSamuraiBDCPhysics::CheckSpectra(){  
+  //  m_Spectra->CheckSpectra();  
+}
+///////////////////////////////////////////////////////////////////////////
+void TSamuraiBDCPhysics::ClearSpectra(){  
+  // To be done
+}
+///////////////////////////////////////////////////////////////////////////
+map< string , TH1*> TSamuraiBDCPhysics::GetSpectra() {
+  /*  if(m_Spectra)
+      return m_Spectra->GetMapHisto();
+      else{
+      map< string , TH1*> empty;
+      return empty;
+      }*/
+  map< string , TH1*> empty;
+  return empty;
+
+} 
+
+///////////////////////////////////////////////////////////////////////////
+void TSamuraiBDCPhysics::WriteSpectra(){
+  // m_Spectra->WriteSpectra();
+}
+///////////////////////////////////////////////////////////////////////////
+void TSamuraiBDCPhysics::AddParameterToCalibrationManager(){
+  CalibrationManager* Cal = CalibrationManager::getInstance();
+
+  // each layer
+  for( int l = 0 ; l < 8 ; ++l){
+    Cal->AddParameter("SamuraiBDC", "L"+ NPL::itoa(l),"BDC_L"+ NPL::itoa(l));
+  }
+
+}
+
+///////////////////////////////////////////////////////////////////////////
+void TSamuraiBDCPhysics::InitializeRootInputRaw(){
+  TChain* inputChain = RootInput::getInstance()->GetChain()   ;
+  inputChain->SetBranchStatus( "SamuraiBDC" , true );
+  // The following line is necessary only for system were the tree is splitted
+  // (older root version). The found argument silenced the Branches not found
+  // warning for non splitted tree.
+  if(inputChain->FindBranch("fDC_*"))
+    inputChain->SetBranchStatus( "fDC_*",true);
+  inputChain->SetBranchAddress( "SamuraiBDC" , &m_EventData );
+
+}
+
+///////////////////////////////////////////////////////////////////////////
+void TSamuraiBDCPhysics::InitializeRootInputPhysics(){
+}
+
+///////////////////////////////////////////////////////////////////////////
+void TSamuraiBDCPhysics::InitializeRootOutput(){
+  TTree* outputTree = RootOutput::getInstance()->GetTree("SamuraiBDC");
+  outputTree->Branch( "SamuraiBDC" , "TSamuraiBDCPhysics" , &m_EventPhysics );
+}
+
+////////////////////////////////////////////////////////////////////////////////
+//            Construct Method to be pass to the DetectorFactory              //
+////////////////////////////////////////////////////////////////////////////////
+NPL::VDetector* TSamuraiBDCPhysics::Construct(){
+  return (NPL::VDetector*) new TSamuraiBDCPhysics();
+}
+
+////////////////////////////////////////////////////////////////////////////////
+//            Registering the construct method to the factory                 //
+////////////////////////////////////////////////////////////////////////////////
+extern "C"{
+  class proxy_samuraiBDC{
+    public:
+      proxy_samuraiBDC(){
+        NPL::DetectorFactory::getInstance()->AddToken("SAMURAIBDC","Samurai");
+        NPL::DetectorFactory::getInstance()->AddToken("SAMURAIBDC1","Samurai");
+        NPL::DetectorFactory::getInstance()->AddToken("SAMURAIBDC2","Samurai");
+        NPL::DetectorFactory::getInstance()->AddDetector("SAMURAIBDC",TSamuraiBDCPhysics::Construct);
+        NPL::DetectorFactory::getInstance()->AddDetector("SAMURAIBDC1",TSamuraiBDCPhysics::Construct);
+        NPL::DetectorFactory::getInstance()->AddDetector("SAMURAIBDC2",TSamuraiBDCPhysics::Construct);
+      }
+  };
+
+  proxy_samuraiBDC p_samuraiBDC;
+}
+
diff --git a/NPLib/Detectors/Samurai/TSamuraiBDCPhysics.h b/NPLib/Detectors/Samurai/TSamuraiBDCPhysics.h
new file mode 100644
index 0000000000000000000000000000000000000000..b3ea3f7dfc4569d2580f4f565b46782007f4859f
--- /dev/null
+++ b/NPLib/Detectors/Samurai/TSamuraiBDCPhysics.h
@@ -0,0 +1,196 @@
+#ifndef TSamuraiBDCPhysics_H
+#define TSamuraiBDCPhysics_H
+/*****************************************************************************
+ * Copyright (C) 2009-2020    this file is part of the NPTool Project        *
+ *                                                                           *
+ * For the licensing terms see $NPTOOL/Licence/NPTool_Licence                *
+ * For the list of contributors see $NPTOOL/Licence/Contributors             *
+ *****************************************************************************/
+
+/*****************************************************************************
+ * Original Author: Adrien MATTA  contact address: matta@lpccaen.in2p3.fr    *
+ *                                                                           *
+ * Creation Date  : May 2021                                                 *
+ * Last update    :                                                          *
+ *---------------------------------------------------------------------------*
+ * Decription:                                                               *
+ *  This class hold SamuraiBDC treated data                                  *
+ *                                                                           *
+ *---------------------------------------------------------------------------*
+ * Comment:                                                                  *
+ *                                                                           *  
+ *                                                                           *
+ *                                                                           *
+ *****************************************************************************/
+// STL
+#include <vector>
+#include <map>
+#include <string>
+
+// NPL
+#include "TSamuraiBDCData.h"
+#include "SamuraiDCIndex.h"
+//#include "TSamuraiBDCSpectra.h"
+#include "NPCalibrationManager.h"
+#include "NPVDetector.h"
+#include "NPInputParser.h"
+#include "NPXmlParser.h"
+
+#if __cplusplus > 199711L && NPMULTITHREADING 
+#include "NPDCReconstructionMT.h"
+#else
+#include "NPDCReconstruction.h"
+#endif
+
+// ROOT 
+#include "TVector3.h" 
+
+// Forward declaration
+//class TSamuraiBDCSpectra;
+
+
+
+class TSamuraiBDCPhysics : public TObject, public NPL::VDetector{
+  public:
+    TSamuraiBDCPhysics();
+    ~TSamuraiBDCPhysics() {};
+
+  public: 
+    void Clear();   
+    void Clear(const Option_t*) {};
+
+  public:
+    //   Provide Physical Multiplicity
+    std::vector<double> DriftLength;
+    std::vector<int> Detector;
+    std::vector<int> Layer;
+    std::vector<int> Wire;
+    std::vector<double> Time;
+    std::vector<double> ToT;
+    std::vector<bool>   Matched;
+    // Computed variable
+    std::vector<TVector3> ParticleDirection;
+    std::vector<TVector3> MiddlePosition;
+
+    double PosX;
+    double PosY;
+    double ThetaX;
+    double PhiY;
+    double devX,devY;
+    TVector3 Dir;
+    int Mult;
+    int MultMean;
+    int PileUp;
+
+  public:
+    // Projected position at given Z plan
+    TVector3 ProjectedPosition(double Z);
+
+  private: // Charateristic of the DC 
+    void AddDC(std::string name, NPL::XmlParser&);//! take the XML file and fill in Wire_X and Layer_Angle
+    std::map<SamuraiDCIndex,double> Wire_X;//! X position of the wires
+    std::map<SamuraiDCIndex,double> Wire_Z;//! Z position of the wires
+    std::map<SamuraiDCIndex,double> Wire_Angle;//! Wire Angle (0 for X, 90 for Y, U and V are typically at +/-30)
+  
+  private: // Analysis
+    double ToTThreshold_H;//! a ToT Low threshold to remove noise
+    double ToTThreshold_L;//! a ToT High threshold to remove noise
+    // since the calibration is a sigmoid there quite a few event at the edge 
+    double DriftLowThreshold;//! Minimum Drift length to keep the hit 
+    double DriftUpThreshold;//! Maximum Drift length to keep the hit
+    double PowerThreshold;//! Maximum P2 minimisation value to keep the track   
+    void RemoveNoise();
+    // Construct the 2D track and ref position at Z=0 and Z=100 based on X,Z and Radius provided
+
+    // Object use to perform the DC reconstruction
+    #if __cplusplus > 199711L && NPMULTITHREADING 
+    NPL::DCReconstructionMT m_reconstruction;//!
+    #else
+    NPL::DCReconstruction m_reconstruction;//!
+    #endif
+
+  public: //   Innherited from VDetector Class
+
+    // Read stream at ConfigFile to pick-up parameters of detector (Position,...) using Token
+    void ReadConfiguration(NPL::InputParser) ;
+
+
+    // Add Parameter to the CalibrationManger
+    void AddParameterToCalibrationManager() ;      
+
+    // Activated associated Branches and link it to the private member DetectorData address
+    // In this method mother Branches (Detector) AND daughter leaf (fDetector_parameter) have to be activated
+    void InitializeRootInputRaw() ;
+
+    // Activated associated Branches and link it to the private member DetectorPhysics address
+    // In this method mother Branches (Detector) AND daughter leaf (parameter) have to be activated
+    void InitializeRootInputPhysics() ;
+
+    // Create associated branches and associated private member DetectorPhysics address
+    void InitializeRootOutput() ;
+
+    // This method is called at each event read from the Input Tree. Aime is to build treat Raw dat in order to extract physical parameter. 
+    void BuildPhysicalEvent() ;
+
+    // Same as above, but only the simplest event and/or simple method are used (low multiplicity, faster algorythm but less efficient ...).
+    // This method aimed to be used for analysis performed during experiment, when speed is requiered.
+    // NB: This method can eventually be the same as BuildPhysicalEvent.
+    void BuildSimplePhysicalEvent() ;
+
+    // Same as above but for online analysis
+    void BuildOnlinePhysicalEvent()  {BuildPhysicalEvent();};
+
+    // Those two method all to clear the Event Physics or Data
+    void ClearEventPhysics() {Clear();}      
+    void ClearEventData()    {m_EventData->Clear();}   
+
+    // Method related to the TSpectra classes, aimed at providing a framework for online applications
+    // Instantiate the Spectra class and the histogramm throught it
+    void InitSpectra();
+    // Fill the spectra hold by the spectra class
+    void FillSpectra();
+    // Used for Online mainly, perform check on the histo and for example change their color if issues are found
+    void CheckSpectra();
+    // Used for Online only, clear all the spectra hold by the Spectra class
+    void ClearSpectra();
+    // Write Spectra to file
+    void WriteSpectra();
+
+  public:      //   Specific to SamuraiBDC Array
+
+    //   Clear The PreTeated object
+    void ClearPreTreatedData()   {m_PreTreatedData->Clear();}
+
+    //   Remove bad channel, calibrate the data and apply threshold
+    void PreTreat();
+
+    // Retrieve raw and pre-treated data
+    TSamuraiBDCData* GetRawData()        const {return m_EventData;}
+    TSamuraiBDCData* GetPreTreatedData() const {return m_PreTreatedData;}
+  
+    double GetPosX(){return PosX;}
+    double GetPosY(){return PosY;}
+    double GetThetaX(){return ThetaX;}
+    double GetPhiY(){return PhiY;}
+    double GetDevX(){return devX;}
+    double GetDevY(){return devY;}
+    int GetPileUp(){return PileUp;}
+
+  private:   //   Root Input and Output tree classes
+    TSamuraiBDCData*         m_EventData;//!
+    TSamuraiBDCData*         m_PreTreatedData;//!
+    TSamuraiBDCPhysics*      m_EventPhysics;//!
+
+
+  private: // Spectra Class
+   // TSamuraiBDCSpectra* m_Spectra; // !
+
+  public: // Spectra Getter
+    std::map< std::string , TH1*> GetSpectra(); 
+
+  public: // Static constructor to be passed to the Detector Factory
+    static NPL::VDetector* Construct();
+    ClassDef(TSamuraiBDCPhysics,1)  // SamuraiBDCPhysics structure
+};
+
+#endif
diff --git a/NPLib/Detectors/Samurai/TSamuraiFDC0Physics.cxx b/NPLib/Detectors/Samurai/TSamuraiFDC0Physics.cxx
index 9c049f7d39684b509bf5f8b8ed57a94eacc07de8..d06456908363c118871c42eccab30ac456da80e9 100644
--- a/NPLib/Detectors/Samurai/TSamuraiFDC0Physics.cxx
+++ b/NPLib/Detectors/Samurai/TSamuraiFDC0Physics.cxx
@@ -447,7 +447,7 @@ void TSamuraiFDC0Physics::InitializeRootInputPhysics(){
 
 ///////////////////////////////////////////////////////////////////////////
 void TSamuraiFDC0Physics::InitializeRootOutput(){
-  TTree* outputTree = RootOutput::getInstance()->GetTree();
+  TTree* outputTree = RootOutput::getInstance()->GetTree("SamuraiFDC0");
   outputTree->Branch( "SamuraiFDC0" , "TSamuraiFDC0Physics" , &m_EventPhysics );
 }
 
diff --git a/NPLib/Detectors/Samurai/TSamuraiFDC2Physics.cxx b/NPLib/Detectors/Samurai/TSamuraiFDC2Physics.cxx
index cd72f64f4ebbf7f0f1c23e2841e931c2300f623d..a443bad1a8f5776e4bc894771be662dbafcc4bf0 100644
--- a/NPLib/Detectors/Samurai/TSamuraiFDC2Physics.cxx
+++ b/NPLib/Detectors/Samurai/TSamuraiFDC2Physics.cxx
@@ -471,7 +471,7 @@ void TSamuraiFDC2Physics::InitializeRootInputPhysics(){
 
 ///////////////////////////////////////////////////////////////////////////
 void TSamuraiFDC2Physics::InitializeRootOutput(){
-  TTree* outputTree = RootOutput::getInstance()->GetTree();
+  TTree* outputTree = RootOutput::getInstance()->GetTree("SamuraiFDC2");
   outputTree->Branch( "SamuraiFDC2" , "TSamuraiFDC2Physics" , &m_EventPhysics );
 }
 
diff --git a/NPLib/Detectors/Samurai/TSamuraiHodoscopePhysics.cxx b/NPLib/Detectors/Samurai/TSamuraiHodoscopePhysics.cxx
index 52de6e132b1ea805d7c5c2b878d409ff9139cff6..12741636a8431642dda9d1aea17a7f22ffa9456a 100644
--- a/NPLib/Detectors/Samurai/TSamuraiHodoscopePhysics.cxx
+++ b/NPLib/Detectors/Samurai/TSamuraiHodoscopePhysics.cxx
@@ -313,7 +313,7 @@ void TSamuraiHodoscopePhysics::InitializeRootInputPhysics() {
 
 ///////////////////////////////////////////////////////////////////////////
 void TSamuraiHodoscopePhysics::InitializeRootOutput() {
-  TTree* outputTree = RootOutput::getInstance()->GetTree();
+  TTree* outputTree = RootOutput::getInstance()->GetTree("SamuraiHodoscope");
   outputTree->Branch("SamuraiHodoscope", "TSamuraiHodoscopePhysics", &m_EventPhysics);
 }
 
diff --git a/NPLib/Utility/npanalysis.cxx b/NPLib/Utility/npanalysis.cxx
index 9ad671f30a3b2cc5009d85c6833332e9d0945aa2..c33c82ffe0f7e52454dfb56256146aa8aab2ea02 100644
--- a/NPLib/Utility/npanalysis.cxx
+++ b/NPLib/Utility/npanalysis.cxx
@@ -24,6 +24,7 @@ void ProgressDisplay(struct timeval& begin, struct timeval& end, unsigned long&
 int main(int argc , char** argv){
   // command line parsing
   NPOptionManager* myOptionManager = NPOptionManager::getInstance(argc,argv);
+  myOptionManager->SetIsAnalysis();
   std::string inputfilename = myOptionManager->GetRunToReadFile();
   // if input files are not given, use those from TAsciiFile
   if (myOptionManager->IsDefault("DetectorConfiguration")) {
@@ -68,8 +69,8 @@ int main(int argc , char** argv){
       OutputfileName="ResultTree";
   }
 
-  RootOutput::getInstance("Analysis/"+OutputfileName,TreeName);
-  TTree* tree= RootOutput::getInstance()->GetTree();
+  RootOutput::getInstance(OutputfileName,TreeName);
+  //TTree* tree= RootOutput::getInstance()->GetTree();
 
   // Instantiate the detector using a file
   NPL::DetectorManager* myDetector = new NPL::DetectorManager();
@@ -137,7 +138,7 @@ int main(int argc , char** argv){
         // Build the current event
         myDetector->BuildPhysicalEvent();
         // Fill the tree
-        tree->Fill();
+        myDetector->FillOutputTree();
 
         current_tree = Chain->GetTreeNumber()+1;
         //ProgressDisplay(begin,end,treated,inter,nentries,mean_rate,displayed,current_tree,total_tree);
@@ -182,7 +183,7 @@ int main(int argc , char** argv){
         // User Analysis
         UserAnalysis->TreatEvent();
         // Fill the tree      
-        tree->Fill();
+        myDetector->FillOutputTree();
       
         current_tree = Chain->GetTreeNumber()+1;
         //ProgressDisplay(begin,end,treated,inter,nentries,mean_rate,displayed,current_tree,total_tree);
@@ -216,7 +217,7 @@ int main(int argc , char** argv){
         // User Analysis
         UserAnalysis->TreatEvent();
         // Fill the tree      
-        tree->Fill();
+        myDetector->FillOutputTree();
 
         current_tree = Chain->GetTreeNumber()+1;
         //ProgressDisplay(begin,end,treated,inter,nentries,mean_rate,displayed,current_tree,total_tree);
diff --git a/NPSimulation/Simulation.cc b/NPSimulation/Simulation.cc
index ebf11bc194f820e31d21d1d568663580b413b1fd..76bb70d9d24ca41baae4e17fa6eb529270122601 100644
--- a/NPSimulation/Simulation.cc
+++ b/NPSimulation/Simulation.cc
@@ -39,6 +39,7 @@
 int main(int argc, char** argv){
     // Initialize NPOptionManager object
     NPOptionManager* OptionManager  = NPOptionManager::getInstance(argc, argv);
+    OptionManager->SetIsSimulation();
     if(OptionManager->GetVerboseLevel() > 0){
         string line;
         line.resize(80,'*');
@@ -76,7 +77,7 @@ int main(int argc, char** argv){
     ///////////////////////////////////////////////////////////////
     ///////////////// Initializing the Root Output ////////////////
     ///////////////////////////////////////////////////////////////
-    RootOutput::getInstance("Simulation/" + OptionManager->GetOutputFile());
+    RootOutput::getInstance(OptionManager->GetOutputFile());
     
     // Construct the default run manager
     G4RunManager* runManager = new G4RunManager;