+ Add method in NPReaction to calculate total cross section

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;
     }
   }
   

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;

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;
 

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;
 }