diff --git a/NPLib/Physics/NPBeam.cxx b/NPLib/Physics/NPBeam.cxx
index 3d6fe698bd3e2c731a848d6697c4f9e6e2442ede..ff5f5d840cc2b1b1870fd63bf717d5ccf57ca3e2 100644
--- a/NPLib/Physics/NPBeam.cxx
+++ b/NPLib/Physics/NPBeam.cxx
@@ -77,7 +77,7 @@ Beam::Beam(){
   while(gDirectory->FindObjectAny(Form("EnergyHist_%i",offset))!=0)
     ++offset;
 
-  fEnergyHist  = new TH1F(Form("EnergyHist_%i",offset),"EnergyHist",1,0,1);
+  fEnergyHist  = new TH1D(Form("EnergyHist_%i",offset),"EnergyHist",1,0,1);
   fXThetaXHist = new TH2F(Form("XThetaXHis_%i",offset),"XThetaXHis",1,0,1,1,0,1);
   fYPhiYHist   = new TH2F(Form("YPhiYHist_%i",offset),"YPhiYHist",1,0,1,1,0,1);
 }
@@ -112,7 +112,7 @@ Beam::Beam(string name){
   while(gDirectory->FindObjectAny(Form("EnergyHist_%i",offset))!=0)
     ++offset;
 
-  fEnergyHist  = new TH1F(Form("EnergyHist_%i",offset),"EnergyHist",1,0,1);
+  fEnergyHist  = new TH1D(Form("EnergyHist_%i",offset),"EnergyHist",1,0,1);
   fXThetaXHist = new TH2F(Form("XThetaXHis_%i",offset),"XThetaXHis",1,0,1,1,0,1);
   fYPhiYHist   = new TH2F(Form("YPhiYHist_%i",offset),"YPhiYHist",1,0,1,1,0,1);
 }
diff --git a/NPLib/Physics/NPBeam.h b/NPLib/Physics/NPBeam.h
index d7b8fb5d8d9a4026662c8e6f0bfb0d2e364c3678..5d01cba8c20fa01a28b697f1599fa586d959837f 100644
--- a/NPLib/Physics/NPBeam.h
+++ b/NPLib/Physics/NPBeam.h
@@ -27,7 +27,7 @@
 #include <string>
 
 // ROOT header
-#include "TH1F.h"
+#include "TH1D.h"
 #include "TH2F.h"
 #include "TLorentzVector.h"
 #include "TRandom.h"
@@ -69,7 +69,7 @@ namespace NPL{
     double fZEmission;
     double fZProfile;
     // case of user given distribution
-    TH1F* fEnergyHist;
+    TH1D* fEnergyHist;
     TH2F* fXThetaXHist;
     TH2F* fYPhiYHist;
     
@@ -88,7 +88,7 @@ namespace NPL{
     void SetMeanPhiY    (double MeanPhiY)       {fMeanPhiY=MeanPhiY;}
     void SetSigmaThetaX (double SigmaThetaX)    {fSigmaThetaX=SigmaThetaX;}
     void SetSigmaPhiY   (double SigmaPhiY)      {fSigmaPhiY=SigmaPhiY;}
-    void SetEnergyHist  (TH1F*  EnergyHist)     {delete fEnergyHist; fEnergyHist   = EnergyHist;}
+    void SetEnergyHist  (TH1D*  EnergyHist)     {delete fEnergyHist; fEnergyHist   = EnergyHist;}
     void SetXThetaXHist (TH2F*  XThetaXHist)    {delete fXThetaXHist; fXThetaXHist = XThetaXHist;}
     void SetYPhiYHist   (TH2F*  YPhiYHist)      {delete fYPhiYHist; fYPhiYHist     = YPhiYHist;}
     void SetVerboseLevel(int verbose)           {fVerboseLevel = verbose;}
@@ -106,7 +106,7 @@ namespace NPL{
     double    GetMeanPhiY    () const {return fMeanPhiY;}
     double    GetSigmaThetaX () const {return fSigmaThetaX;}
     double    GetSigmaPhiY   () const {return fSigmaPhiY;}
-    TH1F*     GetEnergyHist  () const {return fEnergyHist;}
+    TH1D*     GetEnergyHist  () const {return fEnergyHist;}
     TH2F*     GetXThetaXHist () const {return fXThetaXHist;}
     TH2F*     GetYPhiYHist   () const {return fYPhiYHist;}
     int       GetVerboseLevel() const {return fVerboseLevel;}
diff --git a/NPLib/Physics/NPDecay.cxx b/NPLib/Physics/NPDecay.cxx
index eab03205eea8f25c48a14c7cf84669182959f2bf..bcd26e5433bb1caa251f7ddd685919fce8982b70 100644
--- a/NPLib/Physics/NPDecay.cxx
+++ b/NPLib/Physics/NPDecay.cxx
@@ -37,7 +37,7 @@
 //////////////////////////// Single Decay //////////////////////////////////////
 ////////////////////////////////////////////////////////////////////////////////
 
-NPL::SingleDecay::SingleDecay(std::string mother, double threshold, std::vector<std::string> daughter,std::vector<double> Ex, TH1F* CrossSection){
+NPL::SingleDecay::SingleDecay(std::string mother, double threshold, std::vector<std::string> daughter,std::vector<double> Ex, TH1D* CrossSection){
 
   m_MotherName = mother;
   m_DaughterName = daughter;
@@ -277,7 +277,7 @@ void NPL::Decay::ReadConfiguration(std::string MotherName, NPL::InputParser pars
       double LT = blocks[i]->GetDouble("LifeTime","ns");
       m_Shoot = blocks[i]->GetVectorInt("Shoot");
 
-      TH1F* h = NULL; 
+      TH1D* h = NULL; 
       std::vector<std::string> cs ;
 
 
diff --git a/NPLib/Physics/NPDecay.h b/NPLib/Physics/NPDecay.h
index 24f8646be1ad8bac91abea34f7a5d2a8ef12659c..b13a04f47a73ca8f8863a4bfcf785015515d606e 100644
--- a/NPLib/Physics/NPDecay.h
+++ b/NPLib/Physics/NPDecay.h
@@ -47,7 +47,7 @@ namespace NPL{
     public: 
         SingleDecay(){};
         SingleDecay(std::string mother,double threshold, std::vector<std::string> daughter, 
-        std::vector<double> Ex, TH1F* CrossSection);
+        std::vector<double> Ex, TH1D* CrossSection);
         ~SingleDecay(){};
     
     private:
@@ -60,7 +60,7 @@ namespace NPL{
         std::vector<double> m_DaughterMasses;
         std::vector<double> m_ExcitationEnergies;
         double m_Threshold;
-        TH1F* m_CrossSectionHist; 
+        TH1D* m_CrossSectionHist; 
         TGenPhaseSpace m_TGenPhaseSpace; 
  
    public:
diff --git a/NPLib/Physics/NPFunction.cxx b/NPLib/Physics/NPFunction.cxx
index 5d59b7e8e9983a5df5243a8b92f35ce040a69298..73e96d53cbcdcc7051cca3d8eb1987db4d473157 100644
--- a/NPLib/Physics/NPFunction.cxx
+++ b/NPLib/Physics/NPFunction.cxx
@@ -30,10 +30,10 @@ using namespace NPL;
 namespace NPL{
 
   // Check the type of Filename (root or ASCII) and extract build/extract a 1D histogramm
-  TH1F* Read1DProfile(string filename,string HistName) {  
+  TH1D* Read1DProfile(string filename,string HistName) {  
     ifstream ASCII;
     TFile ROOT;
-    TH1F* h;
+    TH1D* h;
 
     // test whether file format is ASCII or ROOT
     bool type = OpenASCIIorROOTFile(filename, ASCII , ROOT);
@@ -68,17 +68,18 @@ namespace NPL{
       //    cout << xmin << "\t" << xmax << "\t" << mysize << "\t" << dx << endl;
 
       // fill histo
-      h = new TH1F(HistName.c_str(), HistName.c_str(), mysize, xmin, xmax+dx);
+      h = new TH1D(HistName.c_str(), HistName.c_str(), mysize, xmin, xmax+dx);
       for (unsigned int i = 0; i < mysize; i++) {
         int bin = h->FindBin(x[i]);
         h->SetBinContent(bin,w[i]);
-        //      cout << i << "\t" << x[i] << "\t" << bin << "\t" << w[i] << "\t" << h->GetBinContent(bin) << endl;
+	//	cout << i << "\t" << x[i] << "\t" << bin << "\t" << w[i] << "\t" << h->GetBinContent(bin) << endl;
       }
     }
 
     // ROOT File case
     else{
-      h = (TH1F*) gDirectory->FindObjectAny(HistName.c_str());
+      h = (TH1D*) gDirectory->FindObjectAny(HistName.c_str());
+      cout << "ROOT file case" << endl;
       if(!h){
         cout << "Error: Histogram " << HistName << " not found in file " << filename << endl;
         exit(1);
diff --git a/NPLib/Physics/NPFunction.h b/NPLib/Physics/NPFunction.h
index 8150a644e5c321bb6c80aa2eee3bc893d993b851..2b1a2a30814f64c3f12458502aeb7594a6abf32e 100644
--- a/NPLib/Physics/NPFunction.h
+++ b/NPLib/Physics/NPFunction.h
@@ -45,7 +45,7 @@ namespace NPL{
 
 
   // Check the type of Filename (root or ASCII) and extract build/extract a 1D histogramm
-  TH1F* Read1DProfile(string filename,string HistName);
+  TH1D* Read1DProfile(string filename,string HistName);
   
   // Check the type of Filename (root or ASCII) and extract build/extract a 2D histogramm
   TH2F* Read2DProfile(string filename,string HistName);
diff --git a/NPLib/Physics/NPQFS.cxx b/NPLib/Physics/NPQFS.cxx
index 587919787ea6cbb1859a616b0989b37aabe9ad31..17710ab3d2cdf51d5d2e96e7c5ed66809324a08a 100644
--- a/NPLib/Physics/NPQFS.cxx
+++ b/NPLib/Physics/NPQFS.cxx
@@ -170,12 +170,12 @@ void QFS::ReadConfigurationFile(NPL::InputParser parser){
       }
       if(blocks[i]->HasToken("PerpMomentumPath")){
           vector<string> file_perp = blocks[i]->GetVectorString("PerpMomentumPath");
-          TH1F* Perptemp = Read1DProfile(file_perp[0], file_perp[1]);
+          TH1D* Perptemp = Read1DProfile(file_perp[0], file_perp[1]);
           SetPerpMomentumHist(Perptemp);
       }
       if(blocks[i]->HasToken("ParMomentumPath")){
           vector<string> file_par = blocks[i]->GetVectorString("ParMomentumPath");
-          TH1F* Partemp = Read1DProfile(file_par[0], file_par[1]);
+          TH1D* Partemp = Read1DProfile(file_par[0], file_par[1]);
           SetParMomentumHist(Partemp);
       }
       if(blocks[i]->HasToken("Deexcitation")){
diff --git a/NPLib/Physics/NPQFS.h b/NPLib/Physics/NPQFS.h
index c5954ef76ee1cd8eb3dfde21efaf66651ef3c4a1..5637d14714c4d1af87b228a9d25093b750d96ebe 100644
--- a/NPLib/Physics/NPQFS.h
+++ b/NPLib/Physics/NPQFS.h
@@ -53,7 +53,7 @@ using namespace NPL;
 #include "TVector3.h"
 #include "TGraph.h"
 #include "TCanvas.h"
-#include "TH1F.h"
+#include "TH1.h"
 
 using namespace std;
 
@@ -81,8 +81,8 @@ namespace NPL{
     double   fExcitationB;             // Excitation energy in MeV of beam-like ejectile
     double   fMomentumSigma;           // Width of the momentum distribution  (sigma)             
     TVector3 fInternalMomentum;        // Internal momentum of the removed cluster            
-    TH1F*    fPerpMomentumHist;        // Perpendicular momentum distribution in beam at rest frame
-    TH1F*    fParMomentumHist;         // Parallel momentum distribution in beam at rest frame
+    TH1D*    fPerpMomentumHist;        // Perpendicular momentum distribution in beam at rest frame
+    TH1D*    fParMomentumHist;         // Parallel momentum distribution in beam at rest frame
     double   fIsotropic;               
 
     bool fCondEcmPos;
@@ -173,9 +173,9 @@ namespace NPL{
     void SetInternalMomentum(const TVector3& mom) {fInternalMomentum = mom;}
     void SetMomentumSigma(const double& sigma) {fMomentumSigma = sigma;}
     void SetIsAllowed(const bool& val) {fIsAllowed = val;}
-    void SetPerpMomentumHist  (TH1F*  PerpMomentumHist)
+    void SetPerpMomentumHist  (TH1D*  PerpMomentumHist)
         {delete fPerpMomentumHist; fPerpMomentumHist   = PerpMomentumHist;}
-    void SetParMomentumHist  (TH1F*  ParMomentumHist)
+    void SetParMomentumHist  (TH1D*  ParMomentumHist)
         {delete fParMomentumHist; fParMomentumHist   = ParMomentumHist;}
 
     //GETTERS
diff --git a/NPLib/Physics/NPReaction.cxx b/NPLib/Physics/NPReaction.cxx
index e79397f08617ee70ee792f15b88e9bb090840e76..25bd9ea5a78a6e6e4f698d13edab10add892f526 100644
--- a/NPLib/Physics/NPReaction.cxx
+++ b/NPLib/Physics/NPReaction.cxx
@@ -157,7 +157,7 @@ Reaction::Reaction(string reaction) {
   while (gDirectory->FindObjectAny(Form("EnergyHist_%i", offset)) != 0)
     ++offset;
 
-  fCrossSectionHist = new TH1F(Form("EnergyHist_%i", offset), "Reaction_CS", 1, 0, 180);
+  fCrossSectionHist = new TH1D(Form("EnergyHist_%i", offset), "Reaction_CS", 1, 0, 180);
   fDoubleDifferentialCrossSectionHist = NULL;
 
   fshoot3 = true;
@@ -523,7 +523,7 @@ void Reaction::ReadConfigurationFile(NPL::InputParser parser) {
 
     if (blocks[i]->HasToken("CrossSectionPath")) {
       vector<string> file = blocks[i]->GetVectorString("CrossSectionPath");
-      TH1F* CStemp = Read1DProfile(file[0], file[1]);
+      TH1D* CStemp = Read1DProfile(file[0], file[1]);
 
       // multiply CStemp by sin(theta)
       TF1* fsin = new TF1("sin", Form("1/(sin(x*%f/180.))", M_PI), 0, 180);
@@ -536,7 +536,7 @@ void Reaction::ReadConfigurationFile(NPL::InputParser parser) {
       fLabCrossSection = true;
 
       vector<string> file = blocks[i]->GetVectorString("LabCrossSectionPath");
-      TH1F* CStemp = Read1DProfile(file[0], file[1]);
+      TH1D* CStemp = Read1DProfile(file[0], file[1]);
 
       // multiply CStemp by sin(theta)
       TF1* fsin = new TF1("sin", Form("1/(sin(x*%f/180.))", M_PI), 0, 180);
diff --git a/NPLib/Physics/NPReaction.h b/NPLib/Physics/NPReaction.h
index 71388278fdde2cbc1408432c84c3071436e447df..e7f3d8b583ef7781f81345ab39c4d1d9dfe5c2b5 100644
--- a/NPLib/Physics/NPReaction.h
+++ b/NPLib/Physics/NPReaction.h
@@ -40,7 +40,7 @@ using namespace NPL;
 // ROOT header
 #include "TCanvas.h"
 #include "TGraph.h"
-#include "TH1F.h"
+#include "TH1.h"
 #include "TLorentzRotation.h"
 #include "TLorentzVector.h"
 #include "TRandom.h"
@@ -95,9 +95,9 @@ namespace NPL {
     double fExcitation1;                       // Excitation energy in MeV of the beam, useful for isomers
     double fExcitation3;                       // Excitation energy in MeV
     double fExcitation4;                       // Excitation energy in MeV
-    TH1F* fCrossSectionHist;                   // Differential cross section in CM frame
+    TH1D* fCrossSectionHist;                   // Differential cross section in CM frame
     TH2F* fDoubleDifferentialCrossSectionHist; // Diff. CS CM frame vs Beam E
-    TH1F* fExcitationEnergyHist;               // Distribution of Excitation energy
+    TH1D* fExcitationEnergyHist;               // Distribution of Excitation energy
    public:
     // Getters and Setters
     void SetBeamEnergy(const double& eBeam) {
@@ -140,7 +140,7 @@ namespace NPL {
       initializePrecomputeVariable();
     }
     void SetVerboseLevel(const int& verbose) { fVerboseLevel = verbose; }
-    void SetCrossSectionHist(TH1F* CrossSectionHist) {
+    void SetCrossSectionHist(TH1D* CrossSectionHist) {
       delete fCrossSectionHist;
       fCrossSectionHist = CrossSectionHist;
     }
@@ -164,7 +164,7 @@ namespace NPL {
     Particle* GetNucleus3() { return GetParticle3(); }
     Particle* GetNucleus4() { return GetParticle4(); }
 
-    TH1F* GetCrossSectionHist() const { return fCrossSectionHist; }
+    TH1D* GetCrossSectionHist() const { return fCrossSectionHist; }
     int GetVerboseLevel() const { return fVerboseLevel; }
     bool GetShoot3() const { return fshoot3; }
     bool GetShoot4() const { return fshoot4; }