From b5955839701a892098247bfdac52c0d901c4ffc5 Mon Sep 17 00:00:00 2001
From: Nicolas de Sereville <deserevi@ipno.in2p3.fr>
Date: Wed, 27 Jan 2016 17:41:33 +0100
Subject: [PATCH] + Add method in NPReaction to calculate total cross section
---
NPLib/Physics/NPFunction.cxx | 9 +++++++--
NPLib/Physics/NPReaction.cxx | 11 +++++++++++
NPLib/Physics/NPReaction.h | 3 +++
Projects/MUGAST/ShowResults.C | 7 ++++---
4 files changed, 25 insertions(+), 5 deletions(-)
diff --git a/NPLib/Physics/NPFunction.cxx b/NPLib/Physics/NPFunction.cxx
index 22680f365..e7a03c66b 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 e48b763b5..710692eba 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 ac713aaf3..7f93db296 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 c39df1272..e81c790ea 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;
}
--
GitLab