diff --git a/NPLib/Physics/NPFunction.cxx b/NPLib/Physics/NPFunction.cxx
index 22680f365fb2f49a3451974277eb87c5971464f1..e7a03c66b2a6bed16c3f58736c0058d177934708 100644
--- a/NPLib/Physics/NPFunction.cxx
+++ b/NPLib/Physics/NPFunction.cxx
@@ -70,21 +70,26 @@ TH1F* Read1DProfile(string filename,string HistName)
     while (getline(ASCII, LineBuffer)) {
        stringstream iss(LineBuffer);
        if (!(iss >> xb >> wb)) {continue;}   // skip comment lines 
-       cout << xb << "\t" << wb << endl;
        // fill vectors
        x.push_back(xb);
        w.push_back(wb);
+//       cout << xb << "\t" << wb << endl;
        // compute xmin / xmax / size of x array
        if (xb > xmax) xmax = xb;
        if (xb < xmin) xmin = xb;
        size++;
     }
+    Double_t dx = (xmax - xmin) / (size - 1);
+//    cout << xmin << "\t" << xmax << "\t" << size << "\t" << dx << endl;
 
     // fill histo
-    h = new TH1F(HistName.c_str(), HistName.c_str(), size, xmin, xmax);
+//    h = new TH1F(HistName.c_str(), HistName.c_str(), size, xmin-dx/2, xmax+dx/2);
+    h = new TH1F(HistName.c_str(), HistName.c_str(), size, xmin, xmax+dx);
+//    h = new TH1F(HistName.c_str(), HistName.c_str(), size, xmin, xmax);
     for (unsigned int i = 0; i < size; 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;
     }
   }
   
diff --git a/NPLib/Physics/NPReaction.cxx b/NPLib/Physics/NPReaction.cxx
index e48b763b5d411efbcd105d23bf11e7014cb43b39..710692eba9fccf644f5fd2320385de3fc838a5b5 100644
--- a/NPLib/Physics/NPReaction.cxx
+++ b/NPLib/Physics/NPReaction.cxx
@@ -697,6 +697,17 @@ TGraph* Reaction::GetELabVersusThetaCM(double AngleStep_CM){
 }
 
 
+////////////////////////////////////////////////////////////////////////////////////////////
+Double_t Reaction::GetTotalCrossSection() const {
+   Double_t stot = fCrossSectionHist->Integral("width"); // take bin width into account (in deg!)
+   stot *= M_PI/180; // correct so that bin width is in rad
+   stot *= 2*M_PI;   // integration over phi
+
+   return stot;
+}
+
+
+
 ////////////////////////////////////////////////////////////////////////////////////////////
 void Reaction::PrintKinematic(){
 	int size = 360;
diff --git a/NPLib/Physics/NPReaction.h b/NPLib/Physics/NPReaction.h
index ac713aaf3879d0094ee7fb2d93ad089855d77521..7f93db29684a2e902f9dbb8e12a36d1729320684 100644
--- a/NPLib/Physics/NPReaction.h
+++ b/NPLib/Physics/NPReaction.h
@@ -220,6 +220,9 @@ namespace NPL{
       double	GetE_CM_3()	{return ECM_3;}
       double	GetE_CM_4()	{return ECM_4;}
 
+      // calculate total cross section
+      Double_t GetTotalCrossSection() const;
+
       // Print private paremeter
       void Print() const;
 
diff --git a/Projects/MUGAST/ShowResults.C b/Projects/MUGAST/ShowResults.C
index c39df127273eff5bd3c63ec0cba353dfa7e9dd5d..e81c790ea5f7a0d8591321310fc472dafcb10af4 100644
--- a/Projects/MUGAST/ShowResults.C
+++ b/Projects/MUGAST/ShowResults.C
@@ -78,13 +78,14 @@ void ShowResults()
 }
 
 
-
 void CountingRates(Double_t ibeam, Double_t ubt)
 {
    // load event generator file
    NPL::Reaction *reaction = new NPL::Reaction();
-   myReaction->ReadConfigurationFile("30Pdp.reaction");
+   reaction->ReadConfigurationFile("30Pdp.reaction");
+//   reaction->ReadConfigurationFile("11Be_d3He.reaction");
    // get angular distribution
-   TH1F *dsig = myReaction->GetCrossSectionHist();
+   TH1F *dsig = reaction->GetCrossSectionHist();
    dsig->Draw();
+   cout << "total cross section = " << reaction->GetTotalCrossSection() << " mb" << endl;
 }