From aff483b7d50ff511b93e53b550d50213f775a5e2 Mon Sep 17 00:00:00 2001
From: Oliver Kirsebom <oliver.kirsebom@gmail.com>
Date: Wed, 31 May 2017 21:38:18 +0300
Subject: [PATCH] macros
---
allowed.mac | 19 +---
analysis/plot/CreateHisto.C | 16 ++--
analysis/plot/MomentumWindow.C | 148 +++++++++++++++++++-----------
analysis/plot/SuperimpMay2017.C | 30 ++++--
anygamma.mac | 36 ++++++++
betaplus.mac | 30 ++++++
electron.mac | 91 +++++++++----------
forbidden.mac | 8 +-
gamma.mac | 4 +-
output/f20scan/ExpScanF20.C | 156 ++++++++++++++++++++++++++++++++
output/f20scan/ScanF20.C | 110 ++++++++++++++++++++++
src/VKCatcherFoil.cpp | 15 ++-
12 files changed, 523 insertions(+), 140 deletions(-)
create mode 100644 anygamma.mac
create mode 100644 betaplus.mac
create mode 100644 output/f20scan/ExpScanF20.C
create mode 100644 output/f20scan/ScanF20.C
diff --git a/allowed.mac b/allowed.mac
index 501288e..0e86e18 100644
--- a/allowed.mac
+++ b/allowed.mac
@@ -19,7 +19,7 @@
# energy spectrum
/VK/source/betaSpectrum true
/VK/source/betaEndPoint 5.40 MeV
-/VK/source/betaMinEnergy 0.00 MeV
+/VK/source/betaMinEnergy 5.00 MeV
/VK/source/betaMaxEnergy 5.40 MeV
# output
@@ -28,19 +28,6 @@
# run simulations
-/VK/field/maxbfield 0.44 tesla
-/VK/output/openFile output/230517/allowed_bmax4400.root
-/run/beamOn 50000000
-
-/VK/field/maxbfield 0.46 tesla
-/VK/output/openFile output/230517/allowed_bmax4600.root
-/run/beamOn 50000000
-
-/VK/field/maxbfield 0.48 tesla
-/VK/output/openFile output/230517/allowed_bmax4800.root
-/run/beamOn 50000000
-
/VK/field/maxbfield 0.50 tesla
-/VK/output/openFile output/230517/allowed_bmax5000.root
-/run/beamOn 50000000
-
+/VK/output/openFile output/f20scan/b50/allowed_cut.root
+/run/beamOn 20000000
diff --git a/analysis/plot/CreateHisto.C b/analysis/plot/CreateHisto.C
index 48d9df5..a2baf73 100644
--- a/analysis/plot/CreateHisto.C
+++ b/analysis/plot/CreateHisto.C
@@ -3,7 +3,7 @@
* This macro was used to generate histograms for test run May 2017.
*/
-void CreateHisto(double hours = 168, int binWidth = 100, double ROImin = 5600, double ROImax = 7000) // 168 hours = 7 days
+void CreateHisto(double hours = 24, int binWidth = 100, double ROImin = 5505, double ROImax = 6695) // 168 hours = 7 days
{
// set some parameters
Double_t backgrHours = 189.4; // runs 192-199
@@ -13,13 +13,13 @@ void CreateHisto(double hours = 168, int binWidth = 100, double ROImin = 5600, d
Double_t branching = 1e-6;
Double_t dtPileup = 1000; // ns
- TString simdir = "../output/220517/";
+ TString simdir = "../../output/f20scan/b50/";
TString datadir = "/home/oliskir/Desktop/f20/data/sorted/";
const double emin = 0;
const double emax = 1E4;
const int bins = (emax - emin) / binWidth;
- vector<string> fname = {"backgr", "allowed", "forbidden", "pileup", "summing"};
+ vector<string> fname = {"backgr", "allowed_cut_res", "forbidden_res", "pileup_res", "summing_res"};
vector<string> hname = {"backgr", "allowed", "forbidden", "pileup", "summing"};
vector<double> sf;
vector<TH1F*> hist, histVeto;
@@ -34,7 +34,7 @@ void CreateHisto(double hours = 168, int binWidth = 100, double ROImin = 5600, d
TString hn = hname[i];
TString hnV = hn + "Veto";
- if (fn == "backgr") {
+ if (hn == "backgr") {
// add files
TChain * chain = new TChain("data");
chain->Add(datadir+"*0192*");
@@ -69,8 +69,8 @@ void CreateHisto(double hours = 168, int binWidth = 100, double ROImin = 5600, d
t->Draw("EDepSignal>>"+hnV, veto);
// scaling factor
double norm = ((TParameter<double>*)(t->GetUserInfo()->FindObject("NORMALIZATION")))->GetVal(); double x = rate * seconds / norm;
- if (fn == "forbidden") x *= branching;
- else if (fn == "pileup") x *= dtPileup * 1E-9 * rate;
+ if (hn == "forbidden") x *= branching;
+ else if (hn == "pileup") x *= dtPileup * 1E-9 * rate;
sf.push_back(x);
}
@@ -88,8 +88,8 @@ void CreateHisto(double hours = 168, int binWidth = 100, double ROImin = 5600, d
cout << " Name Scaling Counts_ROI" << endl;
cout << "---------------------------------" << endl;
for (int i = 0; i < hist.size(); i++) {
- binx1 = hist[i]->GetXaxis()->FindBin(ROImin);
- binx2 = hist[i]->GetXaxis()->FindBin(ROImax);
+ int binx1 = hist[i]->GetXaxis()->FindBin(ROImin);
+ int binx2 = hist[i]->GetXaxis()->FindBin(ROImax);
double y = hist[i]->Integral(binx1, binx2);
double yV = histVeto[i]->Integral(binx1, binx2);
diff --git a/analysis/plot/MomentumWindow.C b/analysis/plot/MomentumWindow.C
index 5cad251..ef84f18 100644
--- a/analysis/plot/MomentumWindow.C
+++ b/analysis/plot/MomentumWindow.C
@@ -2,26 +2,68 @@
/*
* Momentum window. May 2017
* Combines simulated and experimental data.
+ *
+ * peaks: 1) 0.482 MeV, 2) 0.976 MeV, 3) 1.682 MeV
*/
-void MomentumWindow(Double_t xscale = 0.97, Int_t bins = 2000, Double_t emin = 0, Double_t emax = 2000, Bool_t vetoOn = false) {
+void MomentumWindow(int peak = 3, Double_t xscale = 0.925, Int_t bins = 2000, Double_t emin = 0, Double_t emax = 2000, Bool_t vetoOn = false) {
- TString dir = "../../output/scan/";
- TString expdir = "/home/oliskir/Desktop/F20-IGISOL/data/sorted/";
- TString ss = "_c";
+ TString dir = Form("../../output/scan/e%i/", peak);
+ TString expdir = "/home/oliskir/Desktop/f20/data/sorted/new/";
+
+ double xmin;
+ double xmax;
+ if (peak == 1) {
+ xmin = 3;
+ xmax = 14;
+ }
+ else if (peak == 2) {
+ xmin = 9;
+ xmax = 18;
+ }
+ else if (peak == 3) {
+ xmin = 14;
+ xmax = 30;
+ }
- double calib = 147.2; // Bmax (T) => I/Imax (%)
- calib *= xscale;
+ double ampsPerTesla = 1067;
+ double maxAmps = 700; // Danfysik supply
+ double calib = ampsPerTesla / maxAmps * 100; // Bmax (T) => I/Imax (%)
+ calib *= xscale; // small correction
- const int ncuts = 3;
- double cuts[] = {10, 1200, 1670}; // keV
-
- // simulated data
- vector<double> bmax = {0.110, 0.115, 0.120, 0.125, 0.130, 0.135, 0.140, 0.145, 0.155, 0.160, 0.165, 0.170, 0.175};
+ cout << "Calibration coeffiecient: " << calib << " %/T" << endl;
+
+ const int ncuts = 2;
+ double cuts_min[ncuts], cuts_max[ncuts];
+ if (peak == 1) {
+ cuts_min[0] = 200; cuts_max[0] = 550;
+ cuts_min[1] = 350; cuts_max[1] = 550;
+ }
+ else if (peak == 2) {
+ cuts_min[0] = 600; cuts_max[0] = 1150;
+ cuts_min[1] = 850; cuts_max[1] = 1150;
+ }
+ else if (peak == 3) {
+ cuts_min[0] = 1200; cuts_max[0] = 1950;
+ cuts_min[1] = 1550; cuts_max[1] = 1950;
+ }
+
+ // simulated data
+ vector<double> bmax;
+ if (peak == 1) {
+ bmax = vector<double>{0.030, 0.035, 0.040, 0.045, 0.050, 0.055, 0.060, 0.065, 0.070, 0.075, 0.080, 0.085, 0.090, 0.095};
+ }
+ else if (peak == 2) {
+ bmax = vector<double>{0.070, 0.075, 0.080, 0.085, 0.090, 0.095, 0.100, 0.105, 0.110, 0.115, 0.120, 0.125, 0.130, 0.135};
+ }
+ else if (peak == 3) {
+ bmax = vector<double>{0.110, 0.115, 0.120, 0.125, 0.130, 0.135, 0.140, 0.145, 0.150, 0.155, 0.160, 0.165, 0.170, 0.175};
+ }
+
vector<TString> file;
for (auto& b : bmax) {
- TString fname = dir + Form("bmax%.0f.root", b*1E4);
+ TString fname = dir + Form("bmax%04d.root", (int)(b*1E4));
file.push_back(fname);
}
@@ -34,12 +76,12 @@ void MomentumWindow(Double_t xscale = 0.97, Int_t bins = 2000, Double_t emin = 0
double A = A0 * TMath::Exp(- years * TMath::Log(2) / t12);
cout << "Activity: " << A << " Bq" << endl;
- // intensity of 1.6-MeV line
- double intensity = (0.2522 + 0.0403)/1E3;
+ // intensities
+ double intensity;
+ if (peak == 1) intensity = (15.62 + 4.38)/1E3;
+ else if (peak == 2) intensity = (72.43 + 18.38)/1E3;
+ else if (peak == 3) intensity = (0.2522 + 0.0403)/1E3;
- // duration of measurement in minutes
- double duration = 1;
-
// get histograms
vector<TH1F*> hist;
Int_t i=0;
@@ -51,9 +93,9 @@ void MomentumWindow(Double_t xscale = 0.97, Int_t bins = 2000, Double_t emin = 0
// Normalization
Double_t norm = ((TParameter<double>*)(t->GetUserInfo()->FindObject("NORMALIZATION")))->GetVal();
- // Calculate scaling factor
- Double_t sf = intensity * 60 * duration * A / norm;
- if (i==0) cout << "Scaling factor: " << sf << endl;
+ // Calculate scaling factor (for 1 minute measurement)
+ Double_t sf = intensity * 60 * A / norm;
+ if (i==0) cout << "Simulation normalization factor: " << sf << endl;
// veto cut
TCut cut("");
@@ -71,6 +113,7 @@ void MomentumWindow(Double_t xscale = 0.97, Int_t bins = 2000, Double_t emin = 0
}
// do integrals
+ double ymax = 0;
vector<double> integrals[ncuts];
cout << " I/I_max (%) Counts/min Peak frac. (%)" << endl;
i=0;
@@ -78,8 +121,10 @@ void MomentumWindow(Double_t xscale = 0.97, Int_t bins = 2000, Double_t emin = 0
vector<double> counts;
for (int j = 0; j < ncuts; j++) {
- counts.push_back(h->Integral(h->FindBin(cuts[j]), h->FindBin(emax)));
+ counts.push_back(h->Integral(h->FindBin(cuts_min[j]), h->FindBin(cuts_max[j])));
integrals[j].push_back(counts[j]);
+
+ if (counts[j] > ymax) ymax = counts[j];
}
double frac = counts[2] / counts[0];
@@ -94,24 +139,18 @@ void MomentumWindow(Double_t xscale = 0.97, Int_t bins = 2000, Double_t emin = 0
//----------------------------------------------------------------------//
- const int expncuts = 2;
- double expcuts[] = {1200, 1550}; // keV
-
// experimental data
vector<TString> expfile;
vector<double> current;
- expfile.push_back(expdir+"s0169_00.root"); current.push_back(28.0);
- expfile.push_back(expdir+"s0170_00.root"); current.push_back(26.0);
- expfile.push_back(expdir+"s0171_00.root"); current.push_back(24.0);
- expfile.push_back(expdir+"s0181_00.root"); current.push_back(23.0);
- expfile.push_back(expdir+"s0172_00.root"); current.push_back(22.0);
- expfile.push_back(expdir+"s0176_00.root"); current.push_back(21.0);
- expfile.push_back(expdir+"s0173_00.root"); current.push_back(20.0);
- expfile.push_back(expdir+"s0182_00.root"); current.push_back(19.0);
- expfile.push_back(expdir+"s0174_00.root"); current.push_back(18.0);
- expfile.push_back(expdir+"s0175_00.root"); current.push_back(16.0);
+ vector<double> duration;
- double expduration = 5;
+ for (int i = 0; i <= 28; i++) {
+ int run = 272 + i;
+ expfile.push_back(expdir + Form("s0%i_00.root", run));
+ current.push_back(i);
+ if (i <= 15) duration.push_back(1);
+ else duration.push_back(5);
+ }
// get histograms
vector<TH1F*> exphist;
@@ -122,29 +161,29 @@ void MomentumWindow(Double_t xscale = 0.97, Int_t bins = 2000, Double_t emin = 0
exphist.push_back(new TH1F(hname, hname, bins, emin, emax));
TTree * t = (TTree*)f->Get("data");
TCut cut("");
- if (vetoOn) cut = TCut("M1<=1 && E1<140+0.21*E0");
+ if (vetoOn) cut = TCut("M0 == 1 && M1 <= 1 && E1 < 140 + 0.21*E0");
t->Draw("E0>>"+hname, cut);
- exphist[i]->Scale(1./expduration);
+ exphist[i]->Scale(1./duration[i]);
i++;
cout << fname << " " << hname << endl;
}
// do exp integrals
- vector<double> expintegrals[expncuts];
+ vector<double> expintegrals[ncuts];
cout << " I/I_max (%) Counts/min" << endl;
i=0;
for (auto& h : exphist) {
vector<double> counts;
- for (int j = 0; j < expncuts; j++) {
- counts.push_back(h->Integral(h->FindBin(expcuts[j]), h->FindBin(emax)));
+ for (int j = 0; j < ncuts; j++) {
+ counts.push_back(h->Integral(h->FindBin(cuts_min[j]), h->FindBin(cuts_max[j])));
expintegrals[j].push_back(counts[j]);
+
}
cout << Form(" %.0f%%", current[i])
<< Form(" %.0f", counts[0]) << endl;
-
i++;
}
//----------------------------------------------------------------------//
@@ -153,7 +192,7 @@ void MomentumWindow(Double_t xscale = 0.97, Int_t bins = 2000, Double_t emin = 0
TCanvas * c1 = new TCanvas("c1");
// Legend
- Double_t xl1=.63, yl1=0.5, xl2=xl1+0.22, yl2=yl1+0.35;
+ Double_t xl1=.63, yl1=0.5, xl2=xl1+0.22, yl2=yl1+0.25;
TLegend *leg = new TLegend(xl1, yl1, xl2, yl2);
// TGraphs
@@ -169,14 +208,19 @@ void MomentumWindow(Double_t xscale = 0.97, Int_t bins = 2000, Double_t emin = 0
gr[j]->SetLineColor(j+1);
gr[j]->SetLineWidth(2);
mg->Add(gr[j]);
- leg->AddEntry(gr[j], Form("SIM [E > %.0f keV}", cuts[j]), "l");
+ leg->AddEntry(gr[j], Form("SIM [%.0f<E<%.0f]", cuts_min[j], cuts_max[j]), "l");
}
- mg->SetTitle("Momentum window (1.69 MeV)");
+ string sEner;
+ if (peak == 1) sEner = "0.482";
+ else if (peak == 2) sEner = "0.976";
+ else if (peak == 3) sEner = "1.682";
+
+ mg->SetTitle(Form("Momentum window (%s)", sEner.c_str()));
mg->Draw("AC");
// experimental data
- const int N = expncuts;
+ const int N = ncuts;
TGraphErrors * expgr[N];
for (int i = 0; i < N; i++) {
vector<double> ex, ey;
@@ -194,27 +238,25 @@ void MomentumWindow(Double_t xscale = 0.97, Int_t bins = 2000, Double_t emin = 0
expgr[i]->SetMarkerStyle(20+i);
expgr[i]->Draw("samePL");
- leg->AddEntry(expgr[i], Form("EXP [E > %.0f keV]", expcuts[i]), "p");
+ leg->AddEntry(expgr[i], Form("EXP [%.0f<E<%.0f]", cuts_min[i], cuts_max[i]), "p");
}
-
-
auto yaxis = mg->GetYaxis();
yaxis->SetTitle("Counts/minute");
auto xaxis = mg->GetXaxis();
xaxis->SetTitle("I/I_{max} (%)");
- xaxis->SetLimits(14, 30); // along X
- mg->GetHistogram()->SetMaximum(1500); // along
+ xaxis->SetLimits(xmin, xmax); // along X
+ mg->GetHistogram()->SetMaximum(1.3*ymax); // along
mg->GetHistogram()->SetMinimum(0); // Y
leg->Draw("Same");
// text
- TText *t = new TText(14.5, 1300, Form("x-scaling=%.3f", scale));
+ TText *t = new TText(0.8*xmax, 1.15*ymax, Form("multiplier: %.3f", xscale));
t->SetTextAlign(11);
t->SetTextColor(1);
t->SetTextFont(43);
- t->SetTextSize(20);
+ t->SetTextSize(16);
t->Draw("same");
// make pretty
@@ -225,6 +267,6 @@ void MomentumWindow(Double_t xscale = 0.97, Int_t bins = 2000, Double_t emin = 0
gStyle->SetOptStat("");
TString fname = "/home/oliskir/Dropbox/Work_ln-s/F20-IGISOL/analysis/figures/";
- fname += "momentum_window_sim" + ss + ".png";
+ fname += Form("momentum_window_e%i.png", peak);
c1->Print(fname);
}
diff --git a/analysis/plot/SuperimpMay2017.C b/analysis/plot/SuperimpMay2017.C
index b8b8e9a..ab463e6 100644
--- a/analysis/plot/SuperimpMay2017.C
+++ b/analysis/plot/SuperimpMay2017.C
@@ -5,7 +5,7 @@
*/
-void SuperimpMay2017(double xmin = 5000, double xmax = 7500, bool veto = true)
+void SuperimpMay2017(double ROImin = 5505, double ROImax = 6695, double xmin = 5000, double xmax = 7500, bool veto = true)
{
// output
TString ofname = "vk.eps";
@@ -118,7 +118,7 @@ void SuperimpMay2017(double xmin = 5000, double xmax = 7500, bool veto = true)
TFile *f = TFile::Open(ifname);
-
+
// loop over histograms
for (Int_t i=0; i<hname.size(); i++) {
@@ -146,23 +146,26 @@ void SuperimpMay2017(double xmin = 5000, double xmax = 7500, bool veto = true)
// set log-y axis
gPad -> SetLogy(logy);
+
+ // integrals
+ int binx1 = h->GetXaxis()->FindBin(ROImin);
+ int binx2 = h->GetXaxis()->FindBin(ROImax);
+ double counts = h->Integral(binx1, binx2);
// legend
TString label = "l";
if (fillStyle[i] != 0) label = "f";
- leg -> AddEntry(h, legend[i], label);
+ leg -> AddEntry(h, legend[i] + Form(" [%.1f]", counts), label);
leg -> SetTextSize(lsize);
leg -> SetTextFont(rf);
leg -> SetFillColor(0);
leg -> Draw("same");
// title
- h -> SetTitle("^{20}F #beta spectrum (GEANT4)");
+ h -> SetTitle("VETO spectrum");
// get rid of stats
h -> SetStats(kFALSE);
-
- // axes
// X axis
TAxis * xax = h -> GetXaxis();
@@ -194,13 +197,24 @@ void SuperimpMay2017(double xmin = 5000, double xmax = 7500, bool veto = true)
Double_t xt = xmin + 0.0*(xmax-xmin);
xt /= 1000;
Double_t yt = 1.2*ymax;
- TLatex *t1 = new TLatex(xt, yt, "1 week @ 30 kHz");
+ TLatex *t1 = new TLatex(xt, yt, "24 hours @ 30 kHz");
t1->SetTextSize(tsize);
t1->SetTextAngle(0);
t1->SetTextColor(1);
t1->SetTextFont(rf);
t1->Draw("same");
+ xt = xmin + 0.0*(xmax-xmin);
+ xt /= 1000;
+ yt = 3.0*ymax;
+ TLatex *t2 = new TLatex(xt, yt, Form("ROI: %.1f-%.1f MeV", ROImin/1000, ROImax/1000));
+ t2->SetTextSize(tsize);
+ t2->SetTextAngle(0);
+ t2->SetTextColor(1);
+ t2->SetTextFont(rf);
+ t2->Draw("same");
+
+
/*
// TEXT
Double_t xt = xmin;
@@ -235,6 +249,6 @@ void SuperimpMay2017(double xmin = 5000, double xmax = 7500, bool veto = true)
canvas -> Print(ofname);
- TString fPng = "/home/oliskir/Desktop/F20Spectrum.png";
+ TString fPng = "/home/oliskir/Desktop/specs/F20Spectrum.png";
canvas->Print(fPng);
}
diff --git a/anygamma.mac b/anygamma.mac
new file mode 100644
index 0000000..6b87865
--- /dev/null
+++ b/anygamma.mac
@@ -0,0 +1,36 @@
+
+/control/execute init.mac
+
+# source dimensions
+/VK/source/diskShapedSource true
+/VK/source/innerDiameter 0 cm
+/VK/source/outerDiameter 4 mm
+
+# radiation
+/gps/particle gamma
+/gps/position 0.1 0 -29.9 cm
+/gps/ang/type iso
+/gps/ang/mintheta 90 deg
+/gps/ang/maxtheta 180 deg
+
+# energy spectrum
+/gps/energy 1.770 MeV
+
+# output
+/VK/output/zeroSuppression true
+/VK/storeTrajectoryData false
+
+
+# run simulation
+/VK/field/maxbfield 0.125 tesla
+/VK/output/openFile output/g1.root
+/run/beamOn 280000
+
+/VK/field/maxbfield 0.135 tesla
+/VK/output/openFile output/g2.root
+/run/beamOn 280000
+
+/VK/field/maxbfield 0.145 tesla
+/VK/output/openFile output/g3.root
+/run/beamOn 280000
+
diff --git a/betaplus.mac b/betaplus.mac
new file mode 100644
index 0000000..60eb87f
--- /dev/null
+++ b/betaplus.mac
@@ -0,0 +1,30 @@
+
+/control/execute init.mac
+
+# output
+/VK/output/zeroSuppression true
+/VK/storeTrajectoryData true
+
+# source dimensions
+/VK/source/diskShapedSource true
+/VK/source/innerDiameter 0 cm
+/VK/source/outerDiameter 3 mm
+
+# radiation
+/gps/particle e+
+/gps/position 0 0 -30 cm
+/gps/ang/type iso
+/gps/ang/mintheta 90 deg
+/gps/ang/maxtheta 180 deg
+
+# energy spectrum
+/VK/source/betaSpectrum true
+/VK/source/betaEndPoint 0.805 MeV
+/VK/source/betaMinEnergy 0.00 MeV
+/VK/source/betaMaxEnergy 0.805 MeV
+
+# magnetic field
+/VK/field/maxbfield 0.05 tesla
+/VK/output/openFile output/test.root
+/run/beamOn 10000
+
diff --git a/electron.mac b/electron.mac
index 333fc75..0081187 100644
--- a/electron.mac
+++ b/electron.mac
@@ -1,5 +1,6 @@
/control/execute init.mac
+/VK/eventVerbose 0
# output
/VK/output/zeroSuppression true
@@ -16,70 +17,64 @@
/gps/ang/mintheta 90 deg
/gps/ang/maxtheta 180 deg
-#/gps/energy 975.657 keV
-/gps/energy 1682.232 keV
+#/gps/energy 481.697 keV
+/gps/energy 975.657 keV
+#/gps/energy 1682.232 keV
-# magnetic field
-/VK/field/maxbfield 0.1250 tesla
-/VK/output/openFile output/test.root
+/VK/field/maxbfield 0.0700 tesla
+/VK/output/openFile output/scan/e2/bmax0700.root
/run/beamOn 10000
+/VK/field/maxbfield 0.0750 tesla
+/VK/output/openFile output/scan/e2/bmax0750.root
+/run/beamOn 10000
+/VK/field/maxbfield 0.0800 tesla
+/VK/output/openFile output/scan/e2/bmax0800.root
+/run/beamOn 10000
-/VK/eventVerbose 0
+/VK/field/maxbfield 0.0850 tesla
+/VK/output/openFile output/scan/e2/bmax0850.root
+/run/beamOn 10000
+
+/VK/field/maxbfield 0.0900 tesla
+/VK/output/openFile output/scan/e2/bmax0900.root
+/run/beamOn 10000
+
+/VK/field/maxbfield 0.0950 tesla
+/VK/output/openFile output/scan/e2/bmax0950.root
+/run/beamOn 10000
+
+/VK/field/maxbfield 0.1000 tesla
+/VK/output/openFile output/scan/e2/bmax1000.root
+/run/beamOn 10000
+
+/VK/field/maxbfield 0.1050 tesla
+/VK/output/openFile output/scan/e2/bmax1050.root
+/run/beamOn 10000
/VK/field/maxbfield 0.1100 tesla
-##/VK/output/openFile output/scan/bmax1100.root
-##/run/beamOn 10000
+/VK/output/openFile output/scan/e2/bmax1100.root
+/run/beamOn 10000
/VK/field/maxbfield 0.1150 tesla
-##/VK/output/openFile output/scan/bmax1150.root
-##/run/beamOn 10000
+/VK/output/openFile output/scan/e2/bmax1150.root
+/run/beamOn 10000
/VK/field/maxbfield 0.1200 tesla
-##/VK/output/openFile output/scan/bmax1200.root
-##/run/beamOn 10000
+/VK/output/openFile output/scan/e2/bmax1200.root
+/run/beamOn 10000
/VK/field/maxbfield 0.1250 tesla
-##/VK/output/openFile output/scan/bmax1250.root
-##/run/beamOn 10000
+/VK/output/openFile output/scan/e2/bmax1250.root
+/run/beamOn 10000
/VK/field/maxbfield 0.1300 tesla
-##/VK/output/openFile output/scan/bmax1300.root
-##/run/beamOn 10000
+/VK/output/openFile output/scan/e2/bmax1300.root
+/run/beamOn 10000
/VK/field/maxbfield 0.1350 tesla
-##/VK/output/openFile output/scan/bmax1350.root
-##/run/beamOn 10000
-
-/VK/field/maxbfield 0.1400 tesla
-##/VK/output/openFile output/scan/bmax1400.root
-##/run/beamOn 10000
-
-/VK/field/maxbfield 0.1450 tesla
-##/VK/output/openFile output/scan/bmax1450.root
-##/run/beamOn 10000
-
-/VK/field/maxbfield 0.1500 tesla
-##/VK/output/openFile output/scan/bmax1500.root
-##/run/beamOn 10000
-
-/VK/field/maxbfield 0.1550 tesla
-##/VK/output/openFile output/scan/bmax1550.root
-##/run/beamOn 10000
-
-/VK/field/maxbfield 0.1600 tesla
-##/VK/output/openFile output/scan/bmax1600.root
-##/run/beamOn 10000
-
-/VK/field/maxbfield 0.1650 tesla
-##/VK/output/openFile output/scan/bmax1650.root
-##/run/beamOn 10000
+/VK/output/openFile output/scan/e2/bmax1350.root
+/run/beamOn 10000
-/VK/field/maxbfield 0.1700 tesla
-##/VK/output/openFile output/scan/bmax1700.root
-##/run/beamOn 10000
-/VK/field/maxbfield 0.1750 tesla
-##/VK/output/openFile output/scan/bmax1750.root
-##/run/beamOn 10000
diff --git a/forbidden.mac b/forbidden.mac
index cc4d217..64f7eab 100644
--- a/forbidden.mac
+++ b/forbidden.mac
@@ -2,7 +2,7 @@
/control/execute init.mac
# magnetic field strength
-/VK/field/maxbfield 0.48 tesla
+/VK/field/maxbfield 0.49 tesla
# source dimensions
/VK/source/diskShapedSource true
@@ -19,13 +19,13 @@
# energy spectrum
/VK/source/betaSpectrum true
/VK/source/betaEndPoint 7.00 MeV
-/VK/source/betaMinEnergy 0.00 MeV
+/VK/source/betaMinEnergy 5.00 MeV
/VK/source/betaMaxEnergy 7.00 MeV
# output
/VK/output/zeroSuppression true
/VK/storeTrajectoryData false
-/VK/output/openFile output/220517/forbidden.root
+/VK/output/openFile output/f20scan/b49/forbidden.root
# run simulation
-/run/beamOn 1000000
+/run/beamOn 100000
diff --git a/gamma.mac b/gamma.mac
index c187a8d..0aec575 100644
--- a/gamma.mac
+++ b/gamma.mac
@@ -2,7 +2,7 @@
/control/execute init.mac
# magnetic field strength
-/VK/field/maxbfield 0.48 tesla
+/VK/field/maxbfield 0.46 tesla
# source dimensions
/VK/source/diskShapedSource true
@@ -22,7 +22,7 @@
# output
/VK/output/zeroSuppression true
/VK/storeTrajectoryData false
-/VK/output/openFile output/220517/gamma.root
+/VK/output/openFile output/f20scan/b46/gamma.root
# run simulation
/run/beamOn 5000000
diff --git a/output/f20scan/ExpScanF20.C b/output/f20scan/ExpScanF20.C
new file mode 100644
index 0000000..12a5386
--- /dev/null
+++ b/output/f20scan/ExpScanF20.C
@@ -0,0 +1,156 @@
+
+/*
+ * Experiment beta spectrum for various B-field settings
+ * with simulated beta spectrum + measured background superimposed.
+ * 1.6-MeV gamma-ray line is used for normalisation.
+ */
+
+void ExpScanF20(int binWidth = 50, double ROImin = 1500, double ROImax = 7000, double gEff = 0.00296E-2)
+{
+ const double emin = 0;
+ const double emax = 8E3;
+ const int bins = (emax - emin) / binWidth;
+
+ vector<int> bmax = {15, 20, 25, 30, 35, 40, 44, 46, 48, 50};
+ vector<int> runs = {366, 332, 367, 333, 368, 334, 369, 360, 361, 358};
+ vector<double> minutes = {15, 30, 15, 30, 15, 30, 30, 60, 60, 60};
+ vector<double> gCutLow = {1480, 1436, 1480, 1436, 1480, 1436, 1480, 1500, 1500, 1500};
+ vector<double> gCutHigh = {1580, 1534, 1580, 1534, 1580, 1534, 1580, 1600, 1600, 1600};
+/*
+ vector<int> bmax = {20, 30, 40, 46, 48, 50};
+ vector<int> runs = {332, 333, 334, 360, 361, 358};
+ vector<double> minutes = {30, 30, 30, 60, 60, 60};
+ vector<double> gCutLow = {1436, 1436, 1436, 1500, 1500, 1500};
+ vector<double> gCutHigh = {1534, 1534, 1534, 1600, 1600, 1600};
+*/
+ TString datadir = "/home/oliskir/Desktop/f20/data/sorted/new/";
+
+ // canvas
+ TCanvas * c1 = new TCanvas("c1");
+
+
+ int i = 0;
+ for (auto& b : bmax) {
+
+ Double_t hours = minutes[i]/60;
+ Double_t seconds = minutes[i]*60;
+
+ // ----------- background measurement --------------- //
+ double backgrHours = 11.5;
+ // add files
+ TChain * chain = new TChain("data");
+ chain->Add(datadir+"*0311*");
+ // histogram
+ TH1F * hbackgr = new TH1F("backgr", "backgr", bins, emin, emax);
+ // draw
+ chain->Draw("E0>>backgr", "Z0 != 10");
+ // scaling factor
+ hbackgr->Scale(hours/backgrHours);
+ // -------------------------------------------------- //
+
+ // --------- beta-spectrum measurement ------- //
+ TChain * cexp = new TChain("data");
+ cexp->Add(datadir+Form("*0%i*", runs[i]));
+ // histogram
+ TH1F * hbeta = new TH1F("beta", "beta", bins, emin, emax);
+ TH1F * hgamma = new TH1F("gamma", "gamma", 1000, 0, 3000);
+ // draw
+ cexp->Draw("E0>>beta", "Z0 != 10");
+ cexp->Draw("E2>>gamma", "Z2 != 10");
+ // integrate 1.6-MeV peak
+ int binx1 = hgamma->GetXaxis()->FindBin(gCutLow[i]);
+ int binx2 = hgamma->GetXaxis()->FindBin(gCutHigh[i]);
+ double gN = hgamma->Integral(binx1, binx2);
+ // -------------------------------------------------- //
+
+ // -------------- beta-spectrum simulation -------------- //
+ // histogram
+ TH1F * h = new TH1F("dummy", "dummy", bins, emin, emax);
+ double decays = 1;
+ if (i != 0 && i != 2 && i != 4 && i != 6) { // **** TEMPORARY FIX ****
+ TString fname = Form("b%i/allowed_bmax%i_res.root", b, 100*b);
+ TString hname = Form("h%i", b);
+ TFile * f = new TFile(fname, "READ");
+ TTree * t = (TTree*)f->Get("Detector");
+ // histogram
+ h = new TH1F(hname, hname, bins, emin, emax);
+ // draw
+ t->Draw("EDepSignal>>"+hname, "EDepSignal > 100");
+ // scaling factor
+ double norm = ((TParameter<double>*)(t->GetUserInfo()->FindObject("NORMALIZATION")))->GetVal();
+ decays = gN / gEff;
+ double x = decays / norm;
+ // scale
+ h->Scale(x);
+ // add background
+ h->Add(hbackgr);
+ }
+ // ---------------------------------------------------- //
+
+ // draw
+ hbeta->SetMaximum(5*hbeta->GetMaximum());
+ hbeta->SetLineColor(1);
+ hbeta->SetLineWidth(1);
+ hbeta->SetLineStyle(1);
+ hbeta->Draw("hist");
+
+ h->SetLineColor(2);
+ h->Draw("histsame");
+
+ hbackgr->SetLineColor(4);
+ hbackgr->SetLineStyle(2);
+ hbackgr->Draw("histsame");
+
+
+ // integrate
+ binx1 = h->GetXaxis()->FindBin(ROImin);
+ binx2 = h->GetXaxis()->FindBin(ROImax);
+ double simCounts = h->Integral(binx1, binx2);
+ double expCounts = hbeta->Integral(binx1, binx2);
+
+
+ // make pretty
+ hbeta->SetTitle(Form("%.0f min at B_{max} = %.2f T", minutes[i], 0.01*(double)b));
+ hbeta->SetMinimum(0.5);
+ hbeta->GetXaxis()->SetTitle("Energy (keV)");
+ hbeta->GetYaxis()->SetTitle(Form("Counts / %i keV", binWidth));
+
+ double yt = 0.7*h->GetMaximum();
+ TLatex *tx = new TLatex(5500, yt, Form("#splitline{Counts %.1f-%.1f MeV:}{#splitline{SIM: %.1E counts*}{EXP: %.1E counts}}", ROImin/1000, ROImax/1000, simCounts, expCounts));
+ tx->SetTextAlign(11);
+ tx->SetTextColor(1);
+ tx->SetTextFont(43);
+ tx->SetTextSize(16);
+ tx->Draw("same");
+
+ tx = new TLatex(5500, 0.1*yt, "* Normalized to 1.6-MeV #gamma-ray");
+ tx->SetTextAlign(11);
+ tx->SetTextColor(1);
+ tx->SetTextFont(43);
+ tx->SetTextSize(12);
+ tx->Draw("same");
+
+ gPad->SetLogy();
+ gPad->SetGrid();
+ gPad->Update();
+ gStyle->SetOptStat("");
+
+ Double_t xl1=.67, yl1=0.5, xl2=xl1+0.22, yl2=yl1+0.11;
+ TLegend * leg = new TLegend(xl1, yl1, xl2, yl2);
+ leg -> AddEntry(h, "Sim #beta + Backgr", "l");
+ leg -> AddEntry(hbackgr, "Backgr", "l");
+ leg -> AddEntry(hbeta, "Exp #beta", "l");
+ leg -> SetTextSize(15);
+ leg -> SetTextFont(133);
+ leg -> Draw("same");
+
+ cout << " B_max: " << 0.01*b << " T" << endl;
+ cout << "Rate (based on gammas): " << Form("%.1f", decays / seconds / 1E3) << " kHz" << endl;
+ cout << "Rate (based on betas): " << Form("%.1f", expCounts / simCounts * decays / seconds / 1E3) << " kHz" << endl;
+
+ c1->Print(Form("/home/oliskir/Desktop/specs/expf20spec_bmax%i.png", 100*b));
+
+ cin.get();
+ i++;
+ }
+}
diff --git a/output/f20scan/ScanF20.C b/output/f20scan/ScanF20.C
new file mode 100644
index 0000000..40cd60d
--- /dev/null
+++ b/output/f20scan/ScanF20.C
@@ -0,0 +1,110 @@
+
+/*
+ * Simulated beta spectrum for various B-field settings with
+ * experimental background spectrum superimposed.
+ */
+
+void ScanF20(double hours = 1, double rate = 30E3, int binWidth = 25, double ROImin = 1500, double ROImax = 7000)
+{
+ Double_t seconds = hours*3600;
+
+ const double emin = 0;
+ const double emax = 8E3;
+ const int bins = (emax - emin) / binWidth;
+
+ vector<int> bmax = {10, 20, 30, 40, 44, 46, 48, 50};
+
+ // canvas
+ TCanvas * c1 = new TCanvas("c1");
+
+
+ // ----------- background measurement --------------- //
+ TString datadir = "/home/oliskir/Desktop/f20/data/sorted/new/";
+ double backgrHours = 11.5;
+ // add files
+ TChain * chain = new TChain("data");
+ chain->Add(datadir+"*0311*");
+ // histogram
+ TH1F * hbackgr = new TH1F("backgr", "backgr", bins, emin, emax);
+ // draw
+ chain->Draw("E0>>backgr", "Z0 != 10");
+ // scaling factor
+ hbackgr->Scale(hours/backgrHours);
+ // ----------- background measurement --------------- //
+
+
+ int i = 0;
+ for (auto& b : bmax) {
+
+ TString fname = Form("b%i/allowed_bmax%i_res.root", b, 100*b);
+ TString hname = Form("h%i", b);
+
+ cout << " Input: " << fname << endl;
+
+ TFile * f = new TFile(fname, "READ");
+ TTree * t = (TTree*)f->Get("Detector");
+
+ // histograms
+ TH1F * h = new TH1F(hname, hname, bins, emin, emax);
+
+ // draw
+ t->Draw("EDepSignal>>"+hname, "EDepSignal > 100");
+
+ // scaling factor
+ double norm = ((TParameter<double>*)(t->GetUserInfo()->FindObject("NORMALIZATION")))->GetVal();
+ double x = rate * seconds / norm;
+
+ // scale
+ h->Scale(x);
+
+ // draw
+ h->SetMaximum(5*h->GetMaximum());
+ h->SetLineColor(4);
+ h->Draw("hist");
+
+ hbackgr->SetLineColor(2);
+ hbackgr->SetLineStyle(2);
+ hbackgr->Draw("histsame");
+
+ // make pretty
+ h->SetTitle(Form("1 hour at 30 kHz rate and B_{max} = %.2f T", 0.01*(double)b));
+ h->SetMinimum(0.5);
+ int binx1 = h->GetXaxis()->FindBin(ROImin);
+ int binx2 = h->GetXaxis()->FindBin(ROImax);
+ double y = h->Integral(binx1, binx2);
+ h->GetXaxis()->SetTitle("Energy (keV)");
+ h->GetYaxis()->SetTitle(Form("Counts / %i keV", binWidth));
+
+ double ybackgr = hbackgr->Integral(binx1, binx2);
+
+ double yt = 0.3*h->GetMaximum();
+ TText *tx = new TText(4000, yt, Form("Counts/min %.1f-%.1f MeV: %.1E", ROImin/1000, ROImax/1000, (y+ybackgr)/60.));
+ tx->SetTextAlign(11);
+ tx->SetTextColor(1);
+ tx->SetTextFont(43);
+ tx->SetTextSize(16);
+ tx->Draw("same");
+
+ cout << " B_max: " << 0.01*b << " T" << endl;
+ cout << Form(" Counts %.0f-%.0f keV: %.0f", ROImin, ROImax, y) << endl;
+ cout << " Scaling factor: " << x << endl;
+
+ gPad->SetLogy();
+ gPad->SetGrid();
+ gPad->Update();
+ gStyle->SetOptStat("");
+
+ Double_t xl1=.7, yl1=0.5, xl2=xl1+0.18, yl2=yl1+0.11;
+ TLegend * leg = new TLegend(xl1, yl1, xl2, yl2);
+ leg -> AddEntry(h, "^{20}F #beta spectrum", "l");
+ leg -> AddEntry(hbackgr, "Background", "l");
+ leg -> SetTextSize(15);
+ leg -> SetTextFont(133);
+ leg -> Draw("same");
+
+ c1->Print(Form("/home/oliskir/Desktop/specs/f20spec_bmax%i.png", 100*b));
+
+ cin.get();
+ i++;
+ }
+}
diff --git a/src/VKCatcherFoil.cpp b/src/VKCatcherFoil.cpp
index 3f82496..0209399 100644
--- a/src/VKCatcherFoil.cpp
+++ b/src/VKCatcherFoil.cpp
@@ -64,9 +64,11 @@ VKCatcherFoil::VKCatcherFoil(G4RotationMatrix *pRot,
G4Material* labVacuum = fMaterials -> GetMaterial("LABVACUUM");
G4Material* carbon = fMaterials -> GetMaterial("CARBON");
G4Material* alu = fMaterials -> GetMaterial("ALUMINIUM");
+ G4Material* copper = fMaterials -> GetMaterial("COPPER"); // only used for testing
assert(carbon != nullptr); // check that material exists
assert(alu != nullptr);
assert(labVacuum != nullptr);
+ assert(copper != nullptr);
// Initialize variables
G4double startAngle, spanningAngle;
@@ -89,7 +91,7 @@ VKCatcherFoil::VKCatcherFoil(G4RotationMatrix *pRot,
// ------------------------ //
// Foil //
// ------------------------ //
-
+
// solid
fFoil_sol = new G4Tubs("foil_sol", 0., fFoilDiameter/2, fFoilThickness/2,
startAngle=0.*deg, spanningAngle=360.*deg);
@@ -100,6 +102,17 @@ VKCatcherFoil::VKCatcherFoil(G4RotationMatrix *pRot,
new G4PVPlacement(0, G4ThreeVector(), fFoil_log, "foil", fContainer_log, false,
0, fCheckOverlaps);
+/*-----------------just for test purposes ----------
+ // solid
+ fFoil_sol = new G4Tubs("foil_sol", 0., fFoilDiameter/2, 2*mm/2,
+ startAngle=0.*deg, spanningAngle=360.*deg);
+ // logical
+ fFoil_log = new G4LogicalVolume(fFoil_sol, copper, "foil_log");
+
+ // physical
+ new G4PVPlacement(0, G4ThreeVector(), fFoil_log, "foil", fContainer_log, false,
+ 0, fCheckOverlaps);
+*/--------------------------------------------------
// ------------------------ //
// Frame //
--
GitLab