diff --git a/CMakeLists.txt b/CMakeLists.txt
index 7de020ac0443f61604ea0b6e6cc934299d48ec4f..f6562e78e2a562a1f337df439265ea6c22a90978 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -45,6 +45,18 @@ if (ROOT_FOUND)
list(APPEND LIBRARIES ${ROOT_LIBRARIES} Spectrum MathMore)
endif(ROOT_FOUND)
+# Include GSL
+find_package(GSL)
+if (GSL_FOUND)
+ message(STATUS "Found GSL")
+ include_directories(${GSL_INCLUDES})
+ add_definitions("-DHAVE_GSL")
+ if (GSL_ODE_V2)
+ add_definitions("-DGSL_ODE_V2")
+ endif()
+ list(APPEND LIBRARIES ${GSL_LIBRARIES} gslcblas)
+endif (GSL_FOUND)
+
# Include Boost
find_package(Boost REQUIRED)
if (Boost_FOUND)
diff --git a/analysis/VKAdd b/analysis/VKAdd
index 366f7e96aef89b01e88ebc30e27da3eb4d208783..95c13f6c34272a49ed789a761e9c0b386a807ac5 100755
Binary files a/analysis/VKAdd and b/analysis/VKAdd differ
diff --git a/analysis/include/VKBaseData.hh b/analysis/include/VKBaseData.hh
index 5cfc95dd535162587d638661553548fa10839e60..29b09752bc0a971555c1e5e994be86a5be864231 100644
--- a/analysis/include/VKBaseData.hh
+++ b/analysis/include/VKBaseData.hh
@@ -5,6 +5,7 @@
#include <TFile.h>
#include <TTree.h>
#include <TRandom3.h>
+#include <TF1.h>
#include <string>
@@ -34,6 +35,7 @@ class VKBaseData {
virtual double GetNormalization() const = 0;
virtual double GetArea() const = 0;
virtual const VKEvent& GetEvent(size_t id) = 0;
+ virtual bool HasV2Branches() const = 0;
inline void AddEnergyResolution(bool enable) {fAddRes=enable;}
@@ -43,10 +45,11 @@ class VKBaseData {
void SampleResolution(Float_t& e);
void CreateOutputBranches(TTree * t);
+ TF1 * fRes;
TRandom3 * fRNG;
bool fAddRes;
- double fP0, fP1;
+ double fP0, fP1, fTauRel;
TFile *fOutputFile;
TTree *fOutputTree;
diff --git a/analysis/include/VKPileupData.hh b/analysis/include/VKPileupData.hh
index f584a95f9eca4dd4b167fee46e3ed0ec654ad011..9cc9e23d3c0de67c6a89e3683495c8a523c4761e 100644
--- a/analysis/include/VKPileupData.hh
+++ b/analysis/include/VKPileupData.hh
@@ -15,6 +15,7 @@ class VKPileupData : public VKBaseData {
inline size_t GetEntries() const override {return fEntries;}
inline double GetNormalization() const override {return fNormalization;}
inline double GetArea() const override {return fSd->GetArea();}
+ inline bool HasV2Branches() const override {return fHasV2Branches;}
const VKEvent& GetEvent(size_t id) override;
private:
@@ -22,6 +23,7 @@ class VKPileupData : public VKBaseData {
size_t fEntries;
double fNormalization;
+ bool fHasV2Branches;
VKEvent fEvent;
};
diff --git a/analysis/include/VKSimData.hh b/analysis/include/VKSimData.hh
index fabe0d617d92d935f8e473dbe4b716db774d8c38..97c416a4ffcdf72731ea125e9d1dffe65e3d6615 100644
--- a/analysis/include/VKSimData.hh
+++ b/analysis/include/VKSimData.hh
@@ -14,6 +14,7 @@ class VKSimData : public VKBaseData {
inline size_t GetEntries() const override {return fEntries;}
inline double GetNormalization() const override {return fNormalization;}
inline double GetArea() const override {return fArea;}
+ inline bool HasV2Branches() const override {return fHasV2Branches;}
const VKEvent& GetEvent(size_t id) override;
private:
@@ -22,6 +23,7 @@ class VKSimData : public VKBaseData {
size_t fEntries;
double fNormalization;
double fArea;
+ bool fHasV2Branches;
VKEvent fEvent;
TFile *fInputFile;
diff --git a/analysis/include/VKSummingData.hh b/analysis/include/VKSummingData.hh
index 95f535c04fe40620a02299ee0b9fae066310efcd..673f954cd9bc0f1a7d9ee5b186f06d742f5cc926 100644
--- a/analysis/include/VKSummingData.hh
+++ b/analysis/include/VKSummingData.hh
@@ -16,6 +16,7 @@ class VKSummingData : public VKBaseData {
inline size_t GetEntries() const override {return fEntries;}
inline double GetNormalization() const override {return fNormalization;}
inline double GetArea() const override {return fArea;}
+ inline bool HasV2Branches() const override {return fHasV2Branches;}
const VKEvent& GetEvent(size_t id) override;
private:
@@ -25,6 +26,7 @@ class VKSummingData : public VKBaseData {
size_t fEntries;
double fNormalization;
double fArea;
+ bool fHasV2Branches;
VKEvent fEvent;
};
diff --git a/analysis/loop.sh b/analysis/loop.sh
index 3fa9c01944d88caa6a8979740e4b1eafee22e449..d1ac3efa588ce8719b2c7dcc80bba0bfb23ce5f7 100755
--- a/analysis/loop.sh
+++ b/analysis/loop.sh
@@ -1,6 +1,6 @@
-for b in 10 15 20 25 30 35 40 44 46 48 50
+for b in 1500 2000 2500 3000 3500 4000
do
- ./VKAdd --r -o ../output/f20scan/200617/res/b$b.root ../output/f20scan/200617/b$b.root
+ ./VKAdd -r ../output/f20scan/15mm_y5mm_v1/b$b.root
done
diff --git a/analysis/main.cpp b/analysis/main.cpp
index 71034387b29438d77330becf5b8112ada8bdf1f6..ceab07b01ba4423a1097e294a5496cb5f672d36d 100644
--- a/analysis/main.cpp
+++ b/analysis/main.cpp
@@ -33,7 +33,7 @@ int main( int argc, char* argv[] ){
// get data
VKBaseData * sd;
- if (cml.input.size()==1) {
+ if (cml.input.size() == 1) {
if (!cml.pileup) sd = new VKSimData(cml.input[0]);
else sd = new VKPileupData(cml.input[0]);
}
@@ -41,8 +41,28 @@ int main( int argc, char* argv[] ){
sd = new VKSummingData(cml.input[0], cml.input[1]);
}
+ // auto-generate output name
+ if (cml.output.empty()) {
+ cml.output = "output";
+ for (int j = 0; j < cml.input.size(); j++) {
+ auto& i = cml.input[j];
+ size_t p1 = i.rfind("/");
+ size_t p2 = i.rfind(".root");
+ string fname = i.substr(p1, p2-p1);
+ cml.output += fname;
+ if (j < cml.input.size()-1) cml.output += "_";
+ }
+ if (cml.pileup) cml.output += "_pileup";
+ else if (cml.input.size() == 2) cml.output += "_summing";
+ if (cml.resolution) cml.output += "_res";
+ cml.output += ".root";
+ }
+
// set experimental resolution
sd -> AddEnergyResolution(cml.resolution);
+ if (cml.resolution) {
+ cout << " Folding with experimental resolution" << endl;
+ }
// create output file
sd -> CreateOutput(cml.output);
@@ -82,7 +102,7 @@ CommandlineOptions parseoptions(char *const *argv, int argc) {
/* getopt_long stores the option index here. */
int option_index = 0;
- c = getopt_long (argc, argv, "h:o:",
+ c = getopt_long (argc, argv, "hrpo:",
long_options, &option_index);
/* Detect the end of the options. */
@@ -142,9 +162,6 @@ CommandlineOptions parseoptions(char *const *argv, int argc) {
// The input file should come last with no switch
for (int i = optind; i < argc; i++) cml.input.push_back(argv[i]);
- // Default output file
- if (cml.output.empty()) cml.output = "output/o.root";
-
return cml;
}
diff --git a/analysis/src/VKBaseData.cpp b/analysis/src/VKBaseData.cpp
index b17b26f9822c1bdec2d99898d2d7e87d98073285..f7e1a81b46717842c8a7bd0bb4bc78f7524901b4 100644
--- a/analysis/src/VKBaseData.cpp
+++ b/analysis/src/VKBaseData.cpp
@@ -11,11 +11,25 @@ VKBaseData::VKBaseData()
{
// coefficients of the energy-resolution function
- //fP0 = 30.44;
- //fP1 = 17.30;
-
- fP0 = 0.0;
- fP1 = 15; //25.2;
+ //Old parameters:
+ //fP0 = 30.44;
+ //fP1 = 17.30;
+ //fTauRel = 0;
+
+ fP0 = 0.0;
+ fP1 = 5.5;
+ fTauRel = 1.93;
+
+ double emin = -3000;
+ double emax = 100;
+
+ // experimental resolution function
+ if (fTauRel == 0) { // Gaussian
+ fRes = new TF1("fRes", "TMath::Gaus(x, [0], [1], true)", emin, emax);
+ }
+ else { // Gaussian with low-energy tail
+ fRes = new TF1("fRes", "1./[2] * TMath::Exp((x-[0])/[2] + 0.5 * TMath::Power([1]/[2], 2)) * TMath::Erfc(1./sqrt(2) * ((x-[0])/[1] + [1]/[2]))", emin, emax);
+ }
}
@@ -87,24 +101,40 @@ void VKBaseData::CreateOutput(string filename)
void VKBaseData::CreateOutputBranches(TTree * t)
{
- t->Branch("EDepSignal", &_ESout, "_ESout/F");
- t->Branch("EDepVeto" , &_EVout, "_EVout/F");
t->Branch("PhotonsSignal", &_PSout, "_PSout/i");
t->Branch("PhotonsVeto" , &_PVout , "_PVout/i");
- t->Branch("EDepV2Front", &_EF2out, "_EF2out/F");
- t->Branch("EDepV2Signal", &_ES2out, "_ES2out/F");
- t->Branch("EDepV2Veto" , &_EV2out, "_EV2out/F");
+
+ if (this->HasV2Branches()) {
+ t->Branch("EDepV2Front", &_EF2out, "_EF2out/F");
+ t->Branch("EDepV2Signal", &_ES2out, "_ES2out/F");
+ t->Branch("EDepV2Veto" , &_EV2out, "_EV2out/F");
+ }
+ else {
+ t->Branch("EDepSignal", &_ESout, "_ESout/F");
+ t->Branch("EDepVeto" , &_EVout, "_EVout/F");
+ }
}
void VKBaseData::SampleResolution(Float_t& e)
{
if (!fAddRes) return;
-
+
double fwhm = sqrt(pow(fP0,2) + fP1 * e);
double sigma = fwhm / 2.35;
+ double tau = sigma * fTauRel;
+
+ // set parameter values
+ fRes->SetParameter(0, 0);
+ fRes->SetParameter(1, sigma);
+ if (fTauRel > 0) {
+ fRes->SetParameter(0, tau); // shift function so mean is at zero
+ fRes->SetParameter(2, tau);
+ }
+
+ if (fwhm > 50 && e > 50) {
+ e += fRes->GetRandom();
+ }
- if (e > 0) e += fRNG -> Gaus(0, sigma);
-
if (e<=0) e = 0;
}
diff --git a/analysis/src/VKPileupData.cpp b/analysis/src/VKPileupData.cpp
index 6dab791ce01349f4c3593560a94af147074fdb34..62ee33f6949e27c0b7d68ebe78256d78c0060147 100644
--- a/analysis/src/VKPileupData.cpp
+++ b/analysis/src/VKPileupData.cpp
@@ -12,6 +12,9 @@ VKPileupData::VKPileupData(string fname)
if (n % 2 != 0) n -= 1;
fEntries = n/2;
+ // check if tree has V1 or V2 branches
+ fHasV2Branches = fSd->HasV2Branches();
+
// calculate normalization
double norm = fSd -> GetNormalization();
double frac = (double)fEntries / (double)fSd->GetEntries();
diff --git a/analysis/src/VKSimData.cpp b/analysis/src/VKSimData.cpp
index 6e1e3d26a14cfe1ded6e369d3b917e609a4ceee1..8541c896622830319319c9820b1708990c0cbeed 100644
--- a/analysis/src/VKSimData.cpp
+++ b/analysis/src/VKSimData.cpp
@@ -19,6 +19,9 @@ VKSimData::VKSimData(string filename)
fInputTree = (TTree*)fInputFile -> Get("Detector");
SetInputBranches(fInputTree);
+ // check if tree has V1 or V2 branches
+ fHasV2Branches = fInputTree->GetBranch("EDepV2Signal");
+
// number of events going to output file
fEntries = fInputTree -> GetEntries();
diff --git a/analysis/src/VKSummingData.cpp b/analysis/src/VKSummingData.cpp
index a41dd04d62f01ee87f05c555a20daedfd2d26f9b..41770ff1cc5db963fdfc0edac6e56e44f2ebe03f 100644
--- a/analysis/src/VKSummingData.cpp
+++ b/analysis/src/VKSummingData.cpp
@@ -15,6 +15,9 @@ VKSummingData::VKSummingData(string fname1, string fname2)
size_t n2 = fSd2 -> GetEntries();
fEntries = TMath::Min(n1, n2);
+ // check if tree has V1 or V2 branches
+ fHasV2Branches = fSd1->HasV2Branches();
+
// calculate normalization
double nmax = (double)TMath::Max(n1, n2);
double norm1 = fSd1 -> GetNormalization();
diff --git a/cmake/FindGSL.cmake b/cmake/FindGSL.cmake
new file mode 100644
index 0000000000000000000000000000000000000000..64a578fa6904af76e113c7b5d2b1e2ae74e630d7
--- /dev/null
+++ b/cmake/FindGSL.cmake
@@ -0,0 +1,32 @@
+# - Find GSL
+# Find the native GSL includes and library
+#
+# GSL_INCLUDES - where to find gsl/gsl_*.h, etc.
+# GSL_LIBRARIES - List of libraries when using GSL.
+# GSL_FOUND - True if GSL found.
+
+
+if (GSL_INCLUDES)
+ # Already in cache, be silent
+ set (GSL_FIND_QUIETLY TRUE)
+endif (GSL_INCLUDES)
+
+find_path (GSL_INCLUDES gsl/gsl_math.h)
+
+find_library (GSL_LIB NAMES gsl)
+
+set (GSL_CBLAS_LIB "" CACHE FILEPATH "If your program fails to link
+(usually because GSL is not automatically linking a CBLAS and no other
+component of your project provides a CBLAS) then you may need to point
+this variable to a valid CBLAS. Usually GSL is distributed with
+libgslcblas.{a,so} (next to GSL_LIB) which you may use if an optimized
+CBLAS is unavailable.")
+
+set (GSL_LIBRARIES "${GSL_LIB}" "${GSL_CBLAS_LIB}")
+
+# handle the QUIETLY and REQUIRED arguments and set GSL_FOUND to TRUE if
+# all listed variables are TRUE
+include (FindPackageHandleStandardArgs)
+find_package_handle_standard_args (GSL DEFAULT_MSG GSL_LIBRARIES GSL_INCLUDES)
+
+mark_as_advanced (GSL_LIB GSL_CBLAS_LIB GSL_INCLUDES)
diff --git a/include/VKBetaSpectrum.hh b/include/VKBetaSpectrum.hh
new file mode 100644
index 0000000000000000000000000000000000000000..0b109445f656b42b28b74c1fa2764f3620f8cc21
--- /dev/null
+++ b/include/VKBetaSpectrum.hh
@@ -0,0 +1,43 @@
+
+#ifndef VKBetaSpectrum_H
+#define VKBetaSpectrum_H 1
+
+#include "TF1.h"
+
+
+class VKBetaSpectrum
+{
+ public:
+
+ /*
+ * Allowed beta spectrum
+ * @param Z_daughter Charge of daughter nucleus (use negative value for positron emission)
+ * @param A Nucleon numbers
+ * @param endPointKinEnergy End-point kinetic energy
+ * @param shapeFunc Shape function (x = total beta energy in MeV)
+ */
+ VKBetaSpectrum(int Z_daughter, int A, G4double endPointKinEnergy, TF1 * shapeFunc = nullptr);
+ virtual ~VKBetaSpectrum();
+
+ /*
+ * Compute spectrum
+ * @param x Kinetic energy in MeV
+ * @param p Parameters (not used)
+ */
+ double Evaluate(double *x, double *par);
+
+ private:
+ double bremsstrahlungFactor(double W, double W0);
+
+ TF1 * shapeFunc;
+ double Z, A;
+ double W0;
+ double me;
+ double alpha;
+ double gamma;
+ double R;
+ double d1;
+ double M;
+};
+
+#endif
diff --git a/include/VKCatcherFoil.hh b/include/VKCatcherFoil.hh
index 30d16dab15f9973fff34108de685ed8cfa9819ce..873e32759a8485bc5e35c9df5fe824990c3fa93e 100644
--- a/include/VKCatcherFoil.hh
+++ b/include/VKCatcherFoil.hh
@@ -70,6 +70,7 @@ class VKCatcherFoil : public G4PVPlacement
G4double fFrameThickness;
G4double fContainerDiameter;
G4double fContainerLength;
+ G4String fFoilMaterial;
// Geometry
bool fCheckOverlaps;
diff --git a/include/VKDetectorConstruction.hh b/include/VKDetectorConstruction.hh
index f21ccc6088047fa686c0acc85420e32dc831e177..350152766ba5c06686e3820dfceac28d0de94045 100644
--- a/include/VKDetectorConstruction.hh
+++ b/include/VKDetectorConstruction.hh
@@ -150,6 +150,30 @@ class VKDetectorConstruction : public G4VUserDetectorConstruction
*/
inline void SetCheckOverlaps(G4bool b) { fCheckOverlaps=b; }
+ /*
+ * Set position of catcher foil
+ */
+ void SetCatcherFoilPosition(const G4ThreeVector& pos);
+
+ /*
+ * Set thickness of catcher foil
+ */
+ void SetCatcherFoilThickness(G4double thickness);
+
+ /*
+ * Set material of catcher foil
+ */
+ void SetCatcherFoilMaterial(G4String name);
+
+ /*
+ * Set position of beta detector
+ */
+ void SetBetaDetectorPosition(const G4ThreeVector& pos);
+
+ /*
+ * Set position of beta detector V2.0
+ */
+ void SetBetaDetectorV2Position(const G4ThreeVector& pos);
//
// Methods to obtain information about the setup:
@@ -205,6 +229,16 @@ class VKDetectorConstruction : public G4VUserDetectorConstruction
*/
inline const G4ThreeVector&