diff --git a/VeikonKone.dox b/VeikonKone.dox
index e4948149e075f646cb35e9cc759b6a5a1f0490c2..57425b784af33783518d4ad90fefe835e48d15fd 100644
--- a/VeikonKone.dox
+++ b/VeikonKone.dox
@@ -753,7 +753,7 @@ WARN_LOGFILE =
# spaces.
# Note: If this tag is empty the current directory is searched.
-INPUT = /home/oliskir/f20/VeikonKone/include/
+INPUT = /home/oliskir/Desktop/f20/VeikonKone/include/
# This tag can be used to specify the character encoding of the source files
# that doxygen parses. Internally doxygen uses the UTF-8 encoding. Doxygen uses
@@ -2347,4 +2347,4 @@ GENERATE_LEGEND = YES
# The default value is: YES.
# This tag requires that the tag HAVE_DOT is set to YES.
-DOT_CLEANUP = YES
+DOT_CLEANUP = YES
\ No newline at end of file
diff --git a/allowed.mac b/allowed.mac
index d6f64e9d6a5856a56f5a844392654be7fc05cfd9..501288e86316594cbbc8e8d3b0c9b4a3db8a153b 100644
--- a/allowed.mac
+++ b/allowed.mac
@@ -1,5 +1,4 @@
-####/control/execute vis.mac
/control/execute init.mac
# magnetic field strength
@@ -25,8 +24,23 @@
# output
/VK/output/zeroSuppression true
-/VK/output/openFile output/allowed_blinkers_thick.root
-/VK/storeTrajectoryData true
+/VK/storeTrajectoryData false
+
+# 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
-# run simulation
-/run/beamOn 1000000
diff --git a/analysis/VKAdd b/analysis/VKAdd
index 935ff94cf130e8a5efac9ea0a4cb62c4f9d5baca..d66d47a8cb2877bef4af59c662b4663a8945693e 100755
Binary files a/analysis/VKAdd and b/analysis/VKAdd differ
diff --git a/analysis/VKAnalysis.dox b/analysis/VKAnalysis.dox
index 13ece24dbfe71b74171fa65ef1ec5612404fdc3b..69bdf4683c4c6d5302eb5d3e1f53b51bffa49863 100644
--- a/analysis/VKAnalysis.dox
+++ b/analysis/VKAnalysis.dox
@@ -753,7 +753,7 @@ WARN_LOGFILE =
# spaces.
# Note: If this tag is empty the current directory is searched.
-INPUT = /home/oliskir/Desktop/F20-IGISOL/sim/VeikonKone/analysis/include/
+INPUT = /home/oliskir/Desktop/f20/VeikonKone/analysis/include/
# This tag can be used to specify the character encoding of the source files
# that doxygen parses. Internally doxygen uses the UTF-8 encoding. Doxygen uses
diff --git a/analysis/plot/CreateHisto.C b/analysis/plot/CreateHisto.C
new file mode 100644
index 0000000000000000000000000000000000000000..48d9df590ca39b137f16a0ce0d9d7e47458df973
--- /dev/null
+++ b/analysis/plot/CreateHisto.C
@@ -0,0 +1,114 @@
+
+/*
+ * 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
+{
+ // set some parameters
+ Double_t backgrHours = 189.4; // runs 192-199
+ Double_t seconds = hours*3600; // beam time
+ Double_t rate = 30e3; // Hz
+ Double_t vetoThres = 100; // keV
+ Double_t branching = 1e-6;
+ Double_t dtPileup = 1000; // ns
+
+ TString simdir = "../output/220517/";
+ 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> hname = {"backgr", "allowed", "forbidden", "pileup", "summing"};
+ vector<double> sf;
+ vector<TH1F*> hist, histVeto;
+
+ // canvas
+ TCanvas * c1 = new TCanvas("c1");
+
+ int i = 0;
+ for (auto& fn : fname) {
+
+ // output histogram names
+ TString hn = hname[i];
+ TString hnV = hn + "Veto";
+
+ if (fn == "backgr") {
+ // add files
+ TChain * chain = new TChain("data");
+ chain->Add(datadir+"*0192*");
+ chain->Add(datadir+"*0193*");
+ chain->Add(datadir+"*0194*");
+ chain->Add(datadir+"*0195*");
+ chain->Add(datadir+"*0196*");
+ chain->Add(datadir+"*0197*");
+ chain->Add(datadir+"*0198*");
+ chain->Add(datadir+"*0199*");
+ // histograms
+ hist.push_back(new TH1F(hn, hn, bins, emin, emax));
+ histVeto.push_back(new TH1F(hnV, hnV, bins, emin, emax));
+ // veto cut
+ TCut veto("M1 <= 1 && E1 < 140 + 0.21*E0");
+ // draw
+ chain->Draw("E0>>"+hn, "Z0 != 10");
+ chain->Draw("E0>>"+hnV, "Z0 != 10" && veto);
+ // scaling factor
+ sf.push_back(hours/backgrHours);
+ }
+ else {
+ TFile *f = TFile::Open(simdir+Form("%s.root", fn.c_str()));
+ TTree *t = (TTree*)f -> Get("Detector");
+ // histograms
+ hist.push_back(new TH1F(hn, hn, bins, emin, emax));
+ histVeto.push_back(new TH1F(hnV, hnV, bins, emin, emax));
+ // veto cut
+ TCut veto = Form("EDepVeto < %.1f", vetoThres);
+ // draw
+ t->Draw("EDepSignal>>"+hn);
+ 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;
+ sf.push_back(x);
+ }
+
+ // scale
+ hist[i]->Scale(sf[i]);
+ histVeto[i]->Scale(sf[i]);
+
+ i++;
+ }
+
+
+ // print info
+
+ cout << endl;
+ 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);
+ double y = hist[i]->Integral(binx1, binx2);
+ double yV = histVeto[i]->Integral(binx1, binx2);
+
+ cout << " " << fname[i] << " " << sf[i] << " " << Form("%.1f / %.1f",y,yV) << endl;
+ }
+ cout << "--------------------" << endl;
+
+
+ // "detach" the histograms from file
+ for (auto& h : hist) h->SetDirectory(0);
+ for (auto& h : histVeto) h->SetDirectory(0);
+
+
+ // write histograms to new file
+ TString no = "his.root";
+ TFile *fo = TFile::Open(no, "recreate");
+ for (auto& h : hist) h->Write();
+ for (auto& h : histVeto) h->Write();
+ fo -> Close();
+
+ cout << " Histograms written to: " << no << endl << endl;
+}
diff --git a/analysis/plot/MagnetScan_sim.C b/analysis/plot/MagnetScan_sim.C
new file mode 100644
index 0000000000000000000000000000000000000000..f12c20c9d67f25fb785b19f8ce32e1526123cb1c
--- /dev/null
+++ b/analysis/plot/MagnetScan_sim.C
@@ -0,0 +1,126 @@
+
+/*
+ * Magnet scan
+ */
+
+void MagnetScan_sim(Int_t bins = 100, Double_t emin = 500, Double_t emax = 2500, Bool_t vetoOn = false) {
+
+ TString dir = "../../output/scan/";
+ TString ss = "_c";
+
+ vector<TString> file;
+ vector<double> current;
+ double ymax = 0;
+
+ double calib = 0.264;
+
+ file.push_back(dir + "i26" + "_b" + ".root"); current.push_back(26.0);
+ file.push_back(dir + "i24" + "_b" + ".root"); current.push_back(24.0);
+ file.push_back(dir + "i22" + ss + ".root"); current.push_back(22.0);
+ file.push_back(dir + "i20" + ss + ".root"); current.push_back(20.0);
+ file.push_back(dir + "i18" + ss + ".root"); current.push_back(18.0);
+ file.push_back(dir + "i16" + ss + ".root"); current.push_back(16.0);
+
+ // source activity
+ double A0 = 398e3; // Start activity in Bq
+ double years = 8.0; // years since source was made
+ double t12 = 31.6; // half-life
+ 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;
+
+ // duration of measurement in minutes
+ double duration = 1;
+
+ // get histograms
+ vector<TH1F*> hist;
+ Int_t i=0;
+ for (auto& fname : file) {
+ TFile * f = new TFile(fname, "READ");
+ TString hname = Form("I%.0f", current[i]);
+ hist.push_back(new TH1F(hname, hname, bins, emin/calib, emax/calib));
+ TTree * t = (TTree*)f->Get("Detector");
+
+ // 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;
+
+ // veto cut
+ TCut cut("");
+ if (vetoOn) cut = TCut("PhotonsVeto == 0");
+
+ // draw
+ t->Draw("PhotonsSignal>>"+hname, cut);
+
+ // scale
+ hist[i]->Scale(sf);
+
+ // find max
+ double y = hist[i]->GetMaximum();
+ if (y > ymax) ymax = y;
+
+ i++;
+ }
+
+
+ // canvas
+ TCanvas * c1 = new TCanvas("c1");
+
+ // Legend
+ Double_t xl1=.7, yl1=0.6, xl2=xl1+0.15, yl2=yl1+0.25;
+ TLegend *leg = new TLegend(xl1, yl1, xl2, yl2);
+
+ // draw histograms
+ i=0;
+ for (auto& h : hist) {
+ cout << Form("I/I_max=%.0f %%", current[i])
+ << ", Peak integral=" << h->Integral(h->FindBin(1400/calib), h->FindBin(1900/calib)) << endl;
+
+ h->GetXaxis()->SetLimits(emin, emax);
+
+ h->SetMaximum(1.2*ymax);
+
+ int col = (i < 4)? i+1 : i+2;
+ h->SetLineColor(col);
+
+ int sty = i % 2 + 1;
+ h->SetLineStyle(sty);
+
+ int wid = i % 2 + 1;
+ h->SetLineWidth(wid);
+
+ if (i == 0) h->Draw();
+ else h->Draw("same");
+
+ leg -> AddEntry(h, Form("I/I_{max}=%.0f %%", current[i]), "l");
+
+ i++;
+ }
+
+ // make pretty
+ gPad->SetLogy(0);
+ gPad->SetGrid(0,1);
+ gPad->SetTicks();
+
+ gStyle->SetOptStat("");
+
+ hist[0]->GetXaxis()->SetTitle("Energy (keV)");
+ hist[0]->GetYaxis()->SetTitle(Form("Counts / %.0f keV / minute", (emax-emin)/bins));
+
+ TString title = "^{207}Bi magnet scan (Simulation)";
+ hist[0]->SetTitle(title);
+
+ // legend
+ leg -> SetTextSize(15);
+ leg -> SetTextFont(133);
+ leg -> SetFillColor(0);
+ leg -> Draw("same");
+
+ TString fname = "/home/oliskir/Dropbox/Work_ln-s/F20-IGISOL/analysis/figures/";
+ fname += "magnet_scan_sim" + ss + ".png";
+ c1->Print(fname);
+}
diff --git a/analysis/plot/MomentumWindow.C b/analysis/plot/MomentumWindow.C
new file mode 100644
index 0000000000000000000000000000000000000000..5cad25141d2ab0d369cd07e8c5323099df9ea5ba
--- /dev/null
+++ b/analysis/plot/MomentumWindow.C
@@ -0,0 +1,230 @@
+
+/*
+ * Momentum window. May 2017
+ * Combines simulated and experimental data.
+ */
+
+void MomentumWindow(Double_t xscale = 0.97, 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";
+
+ double calib = 147.2; // Bmax (T) => I/Imax (%)
+ calib *= xscale;
+
+ 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};
+
+ vector<TString> file;
+ for (auto& b : bmax) {
+ TString fname = dir + Form("bmax%.0f.root", b*1E4);
+ file.push_back(fname);
+ }
+
+ for (auto&b : bmax) b *= calib;
+
+ // source activity
+ double A0 = 398e3; // Start activity in Bq
+ double years = 8.0; // years since source was made
+ double t12 = 31.6; // half-life
+ 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;
+
+ // duration of measurement in minutes
+ double duration = 1;
+
+ // get histograms
+ vector<TH1F*> hist;
+ Int_t i=0;
+ for (auto& fname : file) {
+ TFile * f = new TFile(fname, "READ");
+ TString hname = Form("sim%i", i+1);
+ hist.push_back(new TH1F(hname, hname, bins, emin, emax));
+ TTree * t = (TTree*)f->Get("Detector");
+
+ // 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;
+
+ // veto cut
+ TCut cut("");
+ if (vetoOn) cut = TCut("EDepVeto == 0");
+
+ // draw
+ t->Draw("EDepSignal>>"+hname, cut);
+
+ // scale
+ hist[i]->Scale(sf);
+
+ i++;
+
+ cout << fname << " " << hname << endl;
+ }
+
+ // do integrals
+ vector<double> integrals[ncuts];
+ cout << " I/I_max (%) Counts/min Peak frac. (%)" << endl;
+ i=0;
+ for (auto& h : hist) {
+
+ vector<double> counts;
+ for (int j = 0; j < ncuts; j++) {
+ counts.push_back(h->Integral(h->FindBin(cuts[j]), h->FindBin(emax)));
+ integrals[j].push_back(counts[j]);
+ }
+
+ double frac = counts[2] / counts[0];
+
+ cout << Form(" %.1f", bmax[i])
+ << Form(" %.0f", counts[0])
+ << Form(" %.0f%%", frac*100.) << endl;
+
+
+ i++;
+ }
+
+
+ //----------------------------------------------------------------------//
+ 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);
+
+ double expduration = 5;
+
+ // get histograms
+ vector<TH1F*> exphist;
+ i=0;
+ for (auto& fname : expfile) {
+ TFile * f = new TFile(fname, "READ");
+ TString hname = Form("exp%i", i+1);
+ 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");
+ t->Draw("E0>>"+hname, cut);
+ exphist[i]->Scale(1./expduration);
+ i++;
+ cout << fname << " " << hname << endl;
+ }
+
+ // do exp integrals
+ vector<double> expintegrals[expncuts];
+ 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)));
+ expintegrals[j].push_back(counts[j]);
+ }
+
+ cout << Form(" %.0f%%", current[i])
+ << Form(" %.0f", counts[0]) << endl;
+
+
+ i++;
+ }
+ //----------------------------------------------------------------------//
+
+ // canvas
+ TCanvas * c1 = new TCanvas("c1");
+
+ // Legend
+ Double_t xl1=.63, yl1=0.5, xl2=xl1+0.22, yl2=yl1+0.35;
+ TLegend *leg = new TLegend(xl1, yl1, xl2, yl2);
+
+ // TGraphs
+ const int M = ncuts;
+ TGraph* gr[M];
+ for (int i = 0; i < M; i++) {
+ gr[i] = new TGraph(integrals[i].size(), bmax.data(), integrals[i].data());
+ }
+
+ // Multigraph
+ TMultiGraph *mg = new TMultiGraph();
+ for (int j = 0; j < M; j++) {
+ 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");
+ }
+
+ mg->SetTitle("Momentum window (1.69 MeV)");
+ mg->Draw("AC");
+
+ // experimental data
+ const int N = expncuts;
+ TGraphErrors * expgr[N];
+ for (int i = 0; i < N; i++) {
+ vector<double> ex, ey;
+ for (auto& y : expintegrals[i]) {
+ ex.push_back(0);
+ ey.push_back(sqrt(y));
+ }
+ expgr[i] = new TGraphErrors(expintegrals[i].size(), current.data(), expintegrals[i].data(), ex.data(), ey.data());
+
+ int col = 4;
+ if (i == 1) col = 6;
+ expgr[i]->SetLineColor(col);
+ expgr[i]->SetLineStyle(1+i);
+ expgr[i]->SetMarkerColor(col);
+ expgr[i]->SetMarkerStyle(20+i);
+ expgr[i]->Draw("samePL");
+
+ leg->AddEntry(expgr[i], Form("EXP [E > %.0f keV]", expcuts[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
+ mg->GetHistogram()->SetMinimum(0); // Y
+
+ leg->Draw("Same");
+
+ // text
+ TText *t = new TText(14.5, 1300, Form("x-scaling=%.3f", scale));
+ t->SetTextAlign(11);
+ t->SetTextColor(1);
+ t->SetTextFont(43);
+ t->SetTextSize(20);
+ t->Draw("same");
+
+ // make pretty
+ gPad->SetLogy(0);
+ gPad->SetGrid(0,1);
+ gPad->SetTicks();
+
+ gStyle->SetOptStat("");
+
+ TString fname = "/home/oliskir/Dropbox/Work_ln-s/F20-IGISOL/analysis/figures/";
+ fname += "momentum_window_sim" + ss + ".png";
+ c1->Print(fname);
+}
diff --git a/analysis/plot/SuperimpMay2017.C b/analysis/plot/SuperimpMay2017.C
new file mode 100644
index 0000000000000000000000000000000000000000..b8b8e9a07fa0db0e5f829c881b0e8555b42c99e5
--- /dev/null
+++ b/analysis/plot/SuperimpMay2017.C
@@ -0,0 +1,240 @@
+
+/*
+ * Superimpose spectra created with CreateHisto.C
+ * May 2017.
+ */
+
+
+void SuperimpMay2017(double xmin = 5000, double xmax = 7500, bool veto = true)
+{
+ // output
+ TString ofname = "vk.eps";
+
+ // input
+ TString ifname = "his.root";
+
+ // files
+ vector<TString> hname;
+ hname.push_back("forbidden");
+ hname.push_back("allowed");
+ hname.push_back("summing");
+ hname.push_back("pileup");
+ hname.push_back("backgr");
+ if (veto) {
+ for (auto& hn : hname) hn += "Veto";
+ }
+
+ // line width
+ vector<Int_t> lineWidth;
+ lineWidth.push_back(2);
+ lineWidth.push_back(1);
+ lineWidth.push_back(1);
+ lineWidth.push_back(1);
+ lineWidth.push_back(1);
+
+ // line colors
+ vector<Int_t> lineColor;
+ lineColor.push_back(1);
+ lineColor.push_back(1);
+ lineColor.push_back(1);
+ lineColor.push_back(1);
+ lineColor.push_back(1);
+
+ // line styles
+ vector<Int_t> lineStyle;
+ lineStyle.push_back(1);
+ lineStyle.push_back(1);
+ lineStyle.push_back(1);
+ lineStyle.push_back(1);
+ lineStyle.push_back(2);
+
+ // fill colors
+ vector<Int_t> fillColor;
+ fillColor.push_back(lineColor[0]);
+ fillColor.push_back(lineColor[1]);
+ fillColor.push_back(lineColor[2]);
+ fillColor.push_back(lineColor[3]);
+ fillColor.push_back(lineColor[4]);
+
+ // fill styles
+ vector<Int_t> fillStyle;
+ fillStyle.push_back(0);
+ fillStyle.push_back(3003);
+ fillStyle.push_back(3004);
+ fillStyle.push_back(3001);
+ fillStyle.push_back(0);
+
+ // legends
+ vector<TString> legend;
+ legend.push_back("2^{+}#rightarrow 0^{+} (br = 1#times10^{-6})");
+ legend.push_back("2^{+}#rightarrow 2^{+}");
+ legend.push_back("#beta#gamma summing");
+ legend.push_back("#beta#beta pile-up (#tau = 1 #mus)");
+ legend.push_back("CR background");
+
+ // y axis
+ Double_t ymin = 0;
+ Double_t ymax = 20;
+ Bool_t logy = true;
+ if (logy) {
+ ymin = 0.01;
+ ymax = 1e3;
+ }
+
+ // Prepare canvas
+ TCanvas *canvas = new TCanvas("c1", "canvas", 900, 700);
+ canvas -> Clear();
+ gROOT->SetStyle("Plain");
+ gPad->UseCurrentStyle();
+ gStyle->SetOptStat("");
+ gStyle->SetFillColor(-1);
+ gStyle->SetLegendBorderSize(0);
+
+ Int_t rf = 133; // times new roman, precision=3
+ Int_t itf = 13; // italic times new roman, precision=3
+ Float_t tsize = 30; // axes title size
+ Float_t lsize = 25; // legend and label size
+ Int_t lw = 1; // default line width
+
+ // Legend
+ Double_t xl1=.5, yl1=0.75, xl2=xl1+0.35, yl2=yl1+0.25;
+ TLegend *leg = new TLegend(xl1, yl1, xl2, yl2);
+
+ canvas->SetFrameBorderMode(0);
+ canvas->SetBorderMode(0);
+ canvas->SetFillColor(kWhite);
+ canvas->SetLeftMargin(0.1);
+
+ Float_t xmarg = 0.005; // percent of canvas size
+ Float_t ymarg = 0.005; // percent of canvas size
+ canvas -> Divide(1, 1, xmarg, ymarg);
+
+ canvas->cd(1);
+ canvas->cd(1)->SetBottomMargin(0.11);
+ canvas->cd(1)->SetTopMargin(0.28);
+ canvas->cd(1)->SetLeftMargin(0.15);
+ canvas->cd(1)->SetTicks(1,1);
+ canvas->cd(1)->SetLogy(1);
+
+
+ TFile *f = TFile::Open(ifname);
+
+ // loop over histograms
+ for (Int_t i=0; i<hname.size(); i++) {
+
+ // get histogram
+ TString hn = hname[i];
+ TH1F *h = (TH1F*)gDirectory->Get(hn);
+
+ // set line and fill properties
+ h -> SetLineWidth(lineWidth[i]);
+ h -> SetLineColor(lineColor[i]);
+ h -> SetLineStyle(lineStyle[i]);
+ h -> SetFillColorAlpha(fillColor[i], 0.60); // 60% transparency
+ h -> SetFillStyle(fillStyle[i]);
+
+ // draw
+ TString opt = "hist";
+ if (i>0) opt = "samehist";
+ h -> Draw(opt);
+
+ // update canvas
+ canvas -> Update();
+
+ // "detach" the histograms from file and close file
+ h -> SetDirectory(0);
+
+ // set log-y axis
+ gPad -> SetLogy(logy);
+
+ // legend
+ TString label = "l";
+ if (fillStyle[i] != 0) label = "f";
+ leg -> AddEntry(h, legend[i], label);
+ leg -> SetTextSize(lsize);
+ leg -> SetTextFont(rf);
+ leg -> SetFillColor(0);
+ leg -> Draw("same");
+
+ // title
+ h -> SetTitle("^{20}F #beta spectrum (GEANT4)");
+
+ // get rid of stats
+ h -> SetStats(kFALSE);
+
+ // axes
+
+ // X axis
+ TAxis * xax = h -> GetXaxis();
+ xax->SetTitleColor(kBlack);
+ xax->SetTitleFont(rf);
+ xax->SetTitleSize(tsize);
+ xax->SetLabelFont(rf);
+ xax->SetLabelSize(lsize);
+ TString s = Form("E (MeV)");
+ xax->SetTitle(s);
+ xax->SetTitleOffset(1);
+ xax->SetLimits(xax->GetXmin()/1000, xax->GetXmax()/1000); // keV -> MeV
+ xax->SetRangeUser(xmin/1000, xmax/1000);
+
+ // Y axis
+ TAxis * yax = h -> GetYaxis();
+ yax->SetTitleFont(rf);
+ yax->SetLabelFont(rf);
+ yax->SetTitleSize(tsize);
+ yax->SetTitleColor(kBlack);
+ yax->SetLabelSize(lsize);
+ s = Form("Counts / 50 keV");
+ yax->SetTitle(s);
+ yax->SetTitleOffset(1.0);
+ h->SetMinimum(ymin);
+ h->SetMaximum(ymax);
+ }
+
+ 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");
+ t1->SetTextSize(tsize);
+ t1->SetTextAngle(0);
+ t1->SetTextColor(1);
+ t1->SetTextFont(rf);
+ t1->Draw("same");
+
+/*
+ // TEXT
+ Double_t xt = xmin;
+ Double_t yt = 10*log10(ymax)/6.3 * ymax;
+ Double_t incr = 2.5*log10(ymax)/6.3;
+
+ TLatex *t1 = new TLatex(xt, yt, "B_{max} = 0.48 T");
+ t1->SetTextSize(tsize);
+ t1->SetTextAngle(0);
+ t1->SetTextColor(1);
+ t1->SetTextFont(rf);
+ t1->Draw("same");
+
+ yt *= incr;
+ t1 = new TLatex(xt, yt, "Duration = 1 week");
+ t1->SetTextSize(tsize);
+ t1->SetTextAngle(0);
+ t1->SetTextColor(1);
+ t1->SetTextFont(rf);
+ t1->Draw("same");
+
+ yt *= incr;
+ t1 = new TLatex(xt, yt, "^{20}F yield = 40 kHz");
+ t1->SetTextSize(tsize);
+ t1->SetTextAngle(0);
+ t1->SetTextColor(1);
+ t1->SetTextFont(rf);
+ t1->Draw("same");
+
+*/
+ f -> Close();
+
+ canvas -> Print(ofname);
+
+ TString fPng = "/home/oliskir/Desktop/F20Spectrum.png";
+ canvas->Print(fPng);
+}
diff --git a/analysis/plot/histo.C b/analysis/plot/histo.C
index 17c283f1564cc462880d12638efecdbe37641846..49198485bc3ccd1be533f68ad86550157fb96acc 100644
--- a/analysis/plot/histo.C
+++ b/analysis/plot/histo.C
@@ -1,5 +1,8 @@
/*
+ * This macro was used to generate histograms
+ * appearing in 2nd 20F proposal for IGISOL (submitted in Sep 2016).
+ *
* source options:
* - muon
* - forbidden
diff --git a/compile.sh b/compile.sh
index f8a0ae00ecf8241be0ff5fc04a8c35d9b497a8f0..0226cae77c0fdbf30306298548a8418b99f39a13 100755
--- a/compile.sh
+++ b/compile.sh
@@ -2,13 +2,12 @@ rm -f CMakeCache.txt
rm -rf CMakeFiles
#@carbonoli
-#cmake -DCMAKE_BUILD_TYPE=DEBUG -DGeant4_DIR=/home/oliskir/src/geant-4.10/geant4.10.3-install/lib/Geant4-10.3.0 $pwd
+cmake -DCMAKE_BUILD_TYPE=DEBUG -DGeant4_DIR=/home/oliskir/src/geant-4.10/geant4.10.3-install/lib/Geant4-10.3.0 $pwd
#@squamish
#cmake -DCMAKE_BUILD_TYPE=DEBUG -DGeant4_DIR=/home/oliskir/src/geant4/geant4.10.2-install/lib/Geant4-10.2.0 $pwd -DCMAKE_BUILD_TYPE=DEBUG
#@stkernfys
-cmake -DCMAKE_BUILD_TYPE=DEBUG -DGeant4_DIR=/home/kernadmin/geant4/geant4.10.2-install/lib/Geant4-10.2.2/ $pwd -DCMAKE_BUILD_TYPE=DEBUG
+#cmake -DCMAKE_BUILD_TYPE=DEBUG -DGeant4_DIR=/home/kernadmin/geant4/geant4.10.2-install/lib/Geant4-10.2.2/ $pwd -DCMAKE_BUILD_TYPE=DEBUG
make -j4
-
diff --git a/electron.mac b/electron.mac
index ab02b4805569fe067d2a1571a0108dd72af04cf1..333fc753ad1d91da573c42686e632f1c41ad3fec 100644
--- a/electron.mac
+++ b/electron.mac
@@ -20,7 +20,66 @@
/gps/energy 1682.232 keV
# magnetic field
-/VK/field/maxbfield 0.145 tesla
-
-/VK/output/openFile output/electron.root
+/VK/field/maxbfield 0.1250 tesla
+/VK/output/openFile output/test.root
/run/beamOn 10000
+
+
+
+/VK/eventVerbose 0
+
+/VK/field/maxbfield 0.1100 tesla
+##/VK/output/openFile output/scan/bmax1100.root
+##/run/beamOn 10000
+
+/VK/field/maxbfield 0.1150 tesla
+##/VK/output/openFile output/scan/bmax1150.root
+##/run/beamOn 10000
+
+/VK/field/maxbfield 0.1200 tesla
+##/VK/output/openFile output/scan/bmax1200.root
+##/run/beamOn 10000
+
+/VK/field/maxbfield 0.1250 tesla
+##/VK/output/openFile output/scan/bmax1250.root
+##/run/beamOn 10000
+
+/VK/field/maxbfield 0.1300 tesla
+##/VK/output/openFile output/scan/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/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 bc8c01abf7ef4871220736cc4559e4d372223ede..cc4d217401a44e47e8e41a606d9353812ce0c2a4 100644
--- a/forbidden.mac
+++ b/forbidden.mac
@@ -1,10 +1,8 @@
-###/control/execute vis.mac
/control/execute init.mac
# magnetic field strength
/VK/field/maxbfield 0.48 tesla
-/VK/detector/volumes/blinkers false
# source dimensions
/VK/source/diskShapedSource true
@@ -21,13 +19,13 @@
# energy spectrum
/VK/source/betaSpectrum true
/VK/source/betaEndPoint 7.00 MeV
-/VK/source/betaMinEnergy 5.00 MeV
+/VK/source/betaMinEnergy 0.00 MeV
/VK/source/betaMaxEnergy 7.00 MeV
# output
/VK/output/zeroSuppression true
-/VK/output/openFile output/forbidden.root
-/VK/storeTrajectoryData true
+/VK/storeTrajectoryData false
+/VK/output/openFile output/220517/forbidden.root
# run simulation
-/run/beamOn 100000
+/run/beamOn 1000000
diff --git a/gamma.mac b/gamma.mac
new file mode 100644
index 0000000000000000000000000000000000000000..c187a8d1b4f3f2ab0bb2952bcf68a3efd135dd2c
--- /dev/null
+++ b/gamma.mac
@@ -0,0 +1,30 @@
+
+/control/execute init.mac
+
+# magnetic field strength
+/VK/field/maxbfield 0.48 tesla
+
+# source dimensions
+/VK/source/diskShapedSource true
+/VK/source/innerDiameter 0 cm
+/VK/source/outerDiameter 4 mm
+
+# radiation
+/gps/particle gamma
+/gps/position 0 0 -30 cm
+/gps/ang/type iso
+/gps/ang/mintheta 90 deg
+/gps/ang/maxtheta 180 deg
+
+# energy spectrum
+/gps/energy 1.60 MeV
+
+# output
+/VK/output/zeroSuppression true
+/VK/storeTrajectoryData false
+/VK/output/openFile output/220517/gamma.root
+
+# run simulation
+/run/beamOn 5000000
+
+
diff --git a/include/VKEventAction.hh b/include/VKEventAction.hh
index ead84d7e30484d0852c589fe5c52e3fb3b3b406a..ab58288e26c7d5a12281a6c8c5de6ca2813ff56d 100644
--- a/include/VKEventAction.hh
+++ b/include/VKEventAction.hh
@@ -63,13 +63,14 @@ class VKEventAction : public G4UserEventAction
G4int fSaveThreshold;
G4int fScintCollID;
+ G4int fVetoCollID;
G4int fSiPMCollID;
G4int fVerbose;
G4int fSiPMThreshold;
G4bool fForcedrawphotons;
G4bool fForcenophotons;
- G4int fPrintModulo;
+ G4int fEvtPct;
};
diff --git a/include/VKMagnet.hh b/include/VKMagnet.hh
index 3dace04f323d29b3072bf76b9ef090a0b0c4cb09..abdc286de2720640f29e08910053167612d8f39f 100644
--- a/include/VKMagnet.hh
+++ b/include/VKMagnet.hh
@@ -100,6 +100,9 @@ class VKMagnet : public G4PVPlacement
G4Tubs* fBrassWalls_sol;
G4LogicalVolume* fBrassWalls_log;
+
+//*** G4Tubs* extra_sol;
+//*** G4LogicalVolume* extra_log;
};
#endif
diff --git a/include/VKRecorderRoot.hh b/include/VKRecorderRoot.hh
index eb8570699c66d270f1605067e81209ea28684dd3..25b3a5295c4bf870fc54c8adb6ca35dd8db9a2aa 100644
--- a/include/VKRecorderRoot.hh
+++ b/include/VKRecorderRoot.hh
@@ -103,7 +103,8 @@ class VKRecorderRoot : public VKRecorderBase {
void CreateDetBranches(); // detection branches
void CreateTraBranches(); // trajectory branches
- void WriteEventToTree();
+ void CreateScintBranches(); // scintillator branches
+ void FillTree();
void ClearRootBranches();
void CloseRootFile();
void WriteMetaDataToRootFile();
@@ -130,6 +131,8 @@ class VKRecorderRoot : public VKRecorderBase {
G4double fEDepSignal;
G4double fEDepVeto;
+ G4double fPhotonsSignal;
+ G4double fPhotonsVeto;
// ROOT output
TFile *fTFile;
diff --git a/include/VKScintillatorSensitiveDetector.hh b/include/VKScintillatorSensitiveDetector.hh
index 0748457fb219aa4c25b93ee618eff35647c3ab0a..877bc6b67a92fb498d05ca46bba0fc5cd88b576e 100644
--- a/include/VKScintillatorSensitiveDetector.hh
+++ b/include/VKScintillatorSensitiveDetector.hh
@@ -49,9 +49,8 @@ class VKScintillatorSensitiveDetector : public G4VSensitiveDetector
virtual void PrintAll();
private:
-
+ G4int fHitsCID;
VKScintillatorHitsCollection* fScintCollection; // <-- typedef in VKScintillatorHit
-
};
#endif
diff --git a/include/VKSiPMScintInfo.hh b/include/VKSiPMScintInfo.hh
new file mode 100644
index 0000000000000000000000000000000000000000..735935ad2b08abb74b82bf10e19c27c5c46be042
--- /dev/null
+++ b/include/VKSiPMScintInfo.hh
@@ -0,0 +1,101 @@
+//
+// ********************************************************************
+// * License and Disclaimer *
+// * *
+// * The Geant4 software is copyright of the Copyright Holders of *
+// * the Geant4 Collaboration. It is provided under the terms and *
+// * conditions of the Geant4 Software License, included in the file *
+// * LICENSE and available at http://cern.ch/geant4/license . These *
+// * include a list of copyright holders. *
+// * *
+// * Neither the authors of this software system, nor their employing *
+// * institutes,nor the agencies providing financial support for this *
+// * work make any representation or warranty, express or implied, *
+// * regarding this software system or assume any liability for its *
+// * use. Please see the license in the file LICENSE and URL above *
+// * for the full disclaimer and the limitation of liability. *
+// * *
+// * This code implementation is the result of the scientific and *
+// * technical work of the GEANT4 collaboration. *
+// * By using, copying, modifying or distributing the software (or *
+// * any work based on the software) you agree to acknowledge its *
+// * use in resulting scientific publications, and indicate your *
+// * acceptance of all terms of the Geant4 Software license. *
+// ********************************************************************
+//
+
+#ifndef VKSiPMScintInfo_h
+#define VKSiPMScintInfo_h 1
+
+#include "VKSiPMScintInfo.hh"
+
+#include "G4VUserEventInformation.hh"
+#include "G4ThreeVector.hh"
+#include "globals.hh"
+
+#include <vector>
+
+class VKSiPMScintInfo : public G4VUserEventInformation
+{
+ public:
+
+ VKSiPMScintInfo();
+ virtual ~VKSiPMScintInfo();
+
+ inline virtual void Print() const {};
+
+ void IncPhotonCount_Scint() {fPhotonCount_Scint++;}
+ void IncPhotonCount_Ceren() {fPhotonCount_Ceren++;}
+ void IncEDep(G4double dep) {fTotE+=dep;}
+ void IncAbsorption() {fAbsorptionCount++;}
+ void IncBoundaryAbsorption() {fBoundaryAbsorptionCount++;}
+ void IncHitCount(G4int i=1) {fHitCount+=i;}
+
+ void SetEWeightPos(const G4ThreeVector& p) {fEWeightPos=p;}
+ void SetReconPos(const G4ThreeVector& p) {fReconPos=p;}
+ void SetConvPos(const G4ThreeVector& p) {fConvPos=p; fConvPosSet=true;}
+ void SetPosMax(const G4ThreeVector& p, G4double edep) {fPosMax=p; fEdepMax=edep;}
+
+ G4int GetPhotonCount_Scint() const {return fPhotonCount_Scint;}
+ G4int GetPhotonCount_Ceren() const {return fPhotonCount_Ceren;}
+ G4int GetHitCount() const {return fHitCount;}
+ G4double GetEDep() const {return fTotE;}
+ G4int GetAbsorptionCount() const {return fAbsorptionCount;}
+ G4int GetBoundaryAbsorptionCount() const {return fBoundaryAbsorptionCount;}
+
+ G4ThreeVector GetEWeightPos() {return fEWeightPos;}
+ G4ThreeVector GetReconPos() {return fReconPos;}
+ G4ThreeVector GetConvPos() {return fConvPos;}
+ G4ThreeVector GetPosMax() {return fPosMax;}
+ G4double GetEDepMax() {return fEdepMax;}
+ G4double IsConvPosSet() {return fConvPosSet;}
+
+ // Gets the total photon count produced
+ G4int GetPhotonCount() {return fPhotonCount_Scint + fPhotonCount_Ceren;}
+
+ void IncSiPMsAboveThreshold() {fSiPMsAboveThreshold++;}
+ G4int GetSiPMsAboveThreshold() {return fSiPMsAboveThreshold;}
+
+ private:
+
+ G4int fHitCount;
+ G4int fPhotonCount_Scint;
+ G4int fPhotonCount_Ceren;
+ G4int fAbsorptionCount;
+ G4int fBoundaryAbsorptionCount;
+
+ G4double fTotE;
+
+ // These only have meaning if totE > 0
+ // If totE = 0 then these wont be set by EndOfEventAction
+ G4ThreeVector fEWeightPos;
+ G4ThreeVector fReconPos; // Also relies on hitCount>0
+ G4ThreeVector fConvPos; // true (initial) conversion position
+ G4bool fConvPosSet;
+ G4ThreeVector fPosMax;
+ G4double fEdepMax;
+
+ G4int fSiPMsAboveThreshold;
+};
+
+#endif
diff --git a/include/VKSiPMSensitiveDetector.hh b/include/VKSiPMSensitiveDetector.hh
index b7f552eda0d7208f8e7840f05889618eeca5a394..3e716818a4587d92194d68159babb617b0202d79 100644
--- a/include/VKSiPMSensitiveDetector.hh
+++ b/include/VKSiPMSensitiveDetector.hh
@@ -64,7 +64,7 @@ class VKSiPMSensitiveDetector : public G4VSensitiveDetector
void SetSiPMPositions(const std::vector<G4ThreeVector>& positions);
private:
-
+ G4int fHitsCID;
VKSiPMHitsCollection* fSiPMHitCollection; // <-- typedef in VKSiPMHit.hh
G4DataVector* fSiPMPositionsX;
diff --git a/include/VKUserEventInformation.hh b/include/VKUserEventInformation.hh
index a0f1840d6cb59f480f234ba8c88b32e052687fa4..fe049185e4b65b12c1eb451a0bb269cf9c1570a3 100644
--- a/include/VKUserEventInformation.hh
+++ b/include/VKUserEventInformation.hh
@@ -27,10 +27,14 @@
#ifndef VKUserEventInformation_h
#define VKUserEventInformation_h 1
+#include "VKSiPMScintInfo.hh"
+
#include "G4VUserEventInformation.hh"
#include "G4ThreeVector.hh"
#include "globals.hh"
+#include <vector>
+
class VKUserEventInformation : public G4VUserEventInformation
{
public:
@@ -38,60 +42,51 @@ class VKUserEventInformation : public G4VUserEventInformation
VKUserEventInformation();
virtual ~VKUserEventInformation();
- inline virtual void Print()const{};
-
- void IncPhotonCount_Scint() {fPhotonCount_Scint++;}
- void IncPhotonCount_Ceren() {fPhotonCount_Ceren++;}
- void IncEDep(G4double dep) {fTotE+=dep;}
- void IncAbsorption() {fAbsorptionCount++;}
- void IncBoundaryAbsorption() {fBoundaryAbsorptionCount++;}
- void IncHitCount(G4int i=1) {fHitCount+=i;}
-
- void SetEWeightPos(const G4ThreeVector& p) {fEWeightPos=p;}
- void SetReconPos(const G4ThreeVector& p) {fReconPos=p;}
- void SetConvPos(const G4ThreeVector& p) {fConvPos=p; fConvPosSet=true;}
- void SetPosMax(const G4ThreeVector& p, G4double edep) {fPosMax=p; fEdepMax=edep;}
-
- G4int GetPhotonCount_Scint() const {return fPhotonCount_Scint;}
- G4int GetPhotonCount_Ceren() const {return fPhotonCount_Ceren;}
- G4int GetHitCount() const {return fHitCount;}
- G4double GetEDep() const {return fTotE;}
- G4int GetAbsorptionCount() const {return fAbsorptionCount;}
- G4int GetBoundaryAbsorptionCount() const {return fBoundaryAbsorptionCount;}
-
- G4ThreeVector GetEWeightPos() {return fEWeightPos;}
- G4ThreeVector GetReconPos() {return fReconPos;}
- G4ThreeVector GetConvPos() {return fConvPos;}
- G4ThreeVector GetPosMax() {return fPosMax;}
- G4double GetEDepMax() {return fEdepMax;}
- G4double IsConvPosSet() {return fConvPosSet;}
+ inline virtual void Print() const{};
+
+ void IncPhotonCount_Scint(size_t id);
+ void IncPhotonCount_Ceren(size_t id);
+ void IncEDep(size_t id, G4double dep);
+ void IncAbsorption(size_t id);
+ void IncBoundaryAbsorption(size_t id);
+ void IncHitCount(size_t id, G4int i=1);
+
+ void SetEWeightPos(size_t id, const G4ThreeVector& p);
+ void SetReconPos(size_t id, const G4ThreeVector& p);
+ void SetPosMax(size_t id, const G4ThreeVector& p, G4double edep);
+
+ G4int GetPhotonCount_Scint(size_t id);
+ G4int GetPhotonCount_Ceren(size_t id);
+ G4int GetHitCount(size_t id);
+ G4double GetEDep(size_t id);
+ G4int GetAbsorptionCount(size_t id);
+ G4int GetBoundaryAbsorptionCount(size_t id);
+
+ G4ThreeVector GetEWeightPos(size_t id);
+ G4ThreeVector GetReconPos(size_t id);
+ G4ThreeVector GetPosMax(size_t id);
+ G4double GetEDepMax(size_t id);
// Gets the total photon count produced
- G4int GetPhotonCount() {return fPhotonCount_Scint + fPhotonCount_Ceren;}
+ G4int GetPhotonCount(size_t id);
+ G4int GetPhotonCount();
- void IncSiPMsAboveThreshold() {fSiPMsAboveThreshold++;}
- G4int GetSiPMsAboveThreshold() {return fSiPMsAboveThreshold;}
+ void IncSiPMsAboveThreshold(size_t id);
+ G4int GetSiPMsAboveThreshold(size_t id);
- private:
- G4int fHitCount;
- G4int fPhotonCount_Scint;
- G4int fPhotonCount_Ceren;
- G4int fAbsorptionCount;
- G4int fBoundaryAbsorptionCount;
+ void SetConvPos(const G4ThreeVector& p) {fConvPos=p; fConvPosSet=true;}
+ G4ThreeVector GetConvPos() {return fConvPos;}
+ G4double IsConvPosSet() {return fConvPosSet;}
- G4double fTotE;
+ private:
+ void CheckCapacity(size_t id) const;
+ void EnsureCapacity(size_t id);
- // These only have meaning if totE > 0
- // If totE = 0 then these wont be set by EndOfEventAction
- G4ThreeVector fEWeightPos;
- G4ThreeVector fReconPos; // Also relies on hitCount>0
G4ThreeVector fConvPos; // true (initial) conversion position
G4bool fConvPosSet;
- G4ThreeVector fPosMax;
- G4double fEdepMax;
-
- G4int fSiPMsAboveThreshold;
+
+ std::vector<VKSiPMScintInfo*> info;
};
diff --git a/init.mac b/init.mac
index be4e23e30b7d1dbc0fd15d8112fac7a4819350c9..8b1326855ee5d85a1fd39eab89bea558e9178d7b 100644
--- a/init.mac
+++ b/init.mac
@@ -4,6 +4,7 @@
/event/verbose 0
/tracking/verbose 0
#/process/em/verbose 0 # Command available for G4EmLivermorePhysics, but not for StandardEM
+/VK/eventVerbose 0
# physics
/VK/physics/opticalPhysics false
diff --git a/scinttest.mac b/scinttest.mac
new file mode 100644
index 0000000000000000000000000000000000000000..b17a398fbd1a7ca78ab6f97e117153d3482c1e62
--- /dev/null
+++ b/scinttest.mac
@@ -0,0 +1,38 @@
+
+# verbosity
+/run/verbose 0
+/event/verbose 0
+/tracking/verbose 0
+#/process/em/verbose 0 # Command available for G4EmLivermorePhysics, but not for StandardEM
+/VK/eventVerbose 0
+
+# physics
+/VK/physics/opticalPhysics true
+
+# setup
+/VK/detector/volumes/betaDetector true
+/VK/detector/volumes/magnet false
+/VK/detector/volumes/catcherFoil false
+/VK/detector/volumes/shield false
+/VK/detector/volumes/blinkers false
+
+# field
+/VK/field/enable false
+
+# initialize
+/run/initialize
+
+# output
+/VK/output/zeroSuppression true
+/VK/storeTrajectoryData true
+
+# radiation generator
+###/gps/position 0 5 33 cm
+###/gps/direction 0 -1 0
+/gps/position 0 0 0 cm
+/gps/direction 0 0 1
+
+/gps/energy 500 keV
+
+/VK/output/openFile output/scinttest.root
+/run/beamOn 100
diff --git a/src/VKBetaDetector.cpp b/src/VKBetaDetector.cpp
index 7744a5e3715a40c395dcde693edf6a4fa94c00df..bd7881e770fb0db546bd000687e85f3e374f41e5 100644
--- a/src/VKBetaDetector.cpp
+++ b/src/VKBetaDetector.cpp
@@ -119,7 +119,7 @@ VKBetaDetector::VKBetaDetector(G4RotationMatrix *pRot,
// ----------------------------------------------- //
- // SIGNAL detector //
+ // SIGNAL scintillator //
// ----------------------------------------------- //
// solid
@@ -136,7 +136,7 @@ VKBetaDetector::VKBetaDetector(G4RotationMatrix *pRot,
// ----------------------------------------------- //
- // Active VETO shield //
+ // VETO scintillator //
// ----------------------------------------------- //
innerRadius = fSignalDiameter/2 + fReflectorThickness;
@@ -198,39 +198,45 @@ VKBetaDetector::VKBetaDetector(G4RotationMatrix *pRot,
// ----------------------------------------------- //
- // SiPMs (modelled as PMTs for now!) //
+ // SiPMs (modelled as PMTs) //
// ----------------------------------------------- //
// SiPM solid and logical volumes
fSiPM_sol = new G4Tubs("SiPM_sol", 0., fSiPMDiameter/2, fSiPMLength/2, 0.*deg, 360.*deg);
fSiPM_log = new G4LogicalVolume(fSiPM_sol, fGlass, "SiPM_log");
- // Place one SiPM on the back side of the SIGNAL detector
+ // Photocathode volumes
+ //
+ // The "photocathode" is a metal slab at the back of the glass that
+ // is only a very rough approximation of the real thing since it only
+ // absorbs or detects the photons based on the efficiency set below.
+ //
+ fPhotocathode_sol = new G4Tubs("photocathode_sol", 0., fSiPMDiameter/2, fSiPMLength/4, 0.*deg, 360.*deg);
+ fPhotocathode_log = new G4LogicalVolume(fPhotocathode_sol, fAluminium, "photocathode_log");
+ pos = G4ThreeVector(0, 0, fSiPMLength/4);
+ new G4PVPlacement(0, pos, fPhotocathode_log, "photocathode", fSiPM_log, false, fCheckOverlaps);
+
+ // Place one SiPM on the back side of the SIGNAL scintillator
displ = -(fContainerLength - fSiPMLength)/2 + fReflectorThickness + fSignalLength;
pos = G4ThreeVector(0, 0, displ);
copyNo = 0;
new G4PVPlacement(0, pos, fSiPM_log, "SiPM", fContainer_log, false, copyNo, fCheckOverlaps);
fSiPMPositions.push_back(pos);
- // Place another SiPM on the back side of the VETO shield
+ // Place another SiPM on the back side of the VETO scintillator
displ = -(fContainerLength - fSiPMLength)/2 + fReflectorThickness + fVetoLength;
pos = G4ThreeVector(0, 0, displ);
copyNo = 1;
new G4PVPlacement(0, pos, fSiPM_log, "SiPM", fContainer_log, false, copyNo, fCheckOverlaps);
fSiPMPositions.push_back(pos);
- // Photocathode solid and logical volumes
- //
- // *** The "photocathode" is a metal slab at the back of the glass that
- // *** is only a very rough approximation of the real thing since it only
- // *** absorbs or detects the photons based on the efficiency set below.
- // *** Note that we do *not* create a physical volume for the photocathode
- //
- fPhotocathode_sol = new G4Tubs("photocathode_sol", 0., fSiPMDiameter/2,
- fSiPMLength/4, 0.*deg, 360.*deg);
- fPhotocathode_log = new G4LogicalVolume(fPhotocathode_sol, fAluminium, "photocathode_log");
-
-
+ // ***IMPORTANT***
+ // VKEventAction::EndOfEventAction() assumes that there are only two SiPMs, with copyNo=0
+ // attached to SIGNAL and copyNo=1 attached to VETO. If you wish to add more SiPMs, you
+ // must make the necessary changes in VKEventAction::EndOfEventAction().
+ // Also check if changes need to be made to VKStackingAction::ClassifyNewTrack() and in
+ // VKRecorderRoot.cpp.
+
// Visualization attributes
VisAttributes();
@@ -332,13 +338,27 @@ void VKBetaDetector::SurfaceProperties()
new G4OpticalSurface("photocathOptSurf", glisur, polished, dielectric_metal);
photocathOptSurf -> SetMaterialPropertiesTable(photocath_mpt);
+ //
+ // 6um Mylar window (same properties as Reflective coating)
+ //
+ G4double mylarRefl[] = {fRefl, fRefl}; //{0.3, 0.3};
+ G4double mylarEff[] = {0.0, 0.0};
+ G4MaterialPropertiesTable* mylarWindow_mpt = new G4MaterialPropertiesTable();
+ mylarWindow_mpt -> AddProperty("REFLECTIVITY", ephoton, mylarRefl, num);
+ mylarWindow_mpt -> AddProperty("EFFICIENCY", ephoton, mylarEff, num);
+ // Create optical surface
+ G4OpticalSurface* mylarWindowOptSurf =
+ new G4OpticalSurface("mylarWindowOptSurf", unified, polished, dielectric_metal);
+ mylarWindowOptSurf -> SetMaterialPropertiesTable(mylarWindow_mpt);
+
+
//
// Create logical skin surfaces
//
- // *** NB: Mylar window gets same properties has reflective coating
+ // *** NB: Mylar window gets same properties as reflective coating
//
new G4LogicalSkinSurface("coating_surf", fCoating_log, reflCoatOptSurf);
- new G4LogicalSkinSurface("window_surf", fWindow_log, reflCoatOptSurf);
+ new G4LogicalSkinSurface("window_surf", fWindow_log, mylarWindowOptSurf);
new G4LogicalSkinSurface("photocathode_surf", fPhotocathode_log, photocathOptSurf);
}
diff --git a/src/VKBlinkers.cpp b/src/VKBlinkers.cpp
index b851833e51463384d3e75d13be302a476a214485..20d97d4b37899e830c261e77e7328526d8616252 100644
--- a/src/VKBlinkers.cpp
+++ b/src/VKBlinkers.cpp
@@ -63,8 +63,10 @@ VKBlinkers::VKBlinkers(G4RotationMatrix *pRot,
// Get materials
G4Material* pmma = fMaterials -> GetMaterial("G4_PLEXIGLASS");
G4Material* labVacuum = fMaterials -> GetMaterial("LABVACUUM");
+ G4Material* alu = fMaterials -> GetMaterial("ALUMINIUM");
assert(pmma != nullptr); // check that material exists
assert(labVacuum != nullptr);
+ assert(alu != nullptr);
// rotation matrices, translation and transformation
G4RotationMatrix doNotRotate;
@@ -125,6 +127,29 @@ VKBlinkers::VKBlinkers(G4RotationMatrix *pRot,
new G4PVPlacement(0, G4ThreeVector(0, 0, -fConeLength/2), fBlinkers_log,
"blinkers", fContainer_log, false, 0, fCheckOverlaps);
+
+/*
+ //-------------------------------------------------------------//
+ // Simple Al blinkers used for tests in April 2017.
+ // - fCutOffAngle = 60*deg;
+ // - fThickness = 3*mm;
+ //
+ // Tube
+ G4Tubs* tube_sol = new G4Tubs("blinkersTube_sol",
+ fInnerDiameter/2,
+ fOuterDiameter/2,
+ fTotalLength/2,
+ startAngle=0.*deg,
+ spanningAngle=360.*deg);
+ // logical volume
+ fBlinkers_log = new G4LogicalVolume(tube_sol, alu, "blinkers_log");
+ // placement
+ new G4PVPlacement(0, G4ThreeVector(0, 0, fTotalLength/2), fBlinkers_log,
+ "blinkers", fContainer_log, false, 0, fCheckOverlaps);
+ //-------------------------------------------------------------//
+*/
+
+
// visualization attributes
VisAttributes();
diff --git a/src/VKDetectorConstruction.cpp b/src/VKDetectorConstruction.cpp
index d59d9439aac7803b50f915d1764766ee5df0c922..54ee73add0fd4b91ef868ecaedc04b07dfcc6833 100644
--- a/src/VKDetectorConstruction.cpp
+++ b/src/VKDetectorConstruction.cpp
@@ -273,17 +273,24 @@ void VKDetectorConstruction::ConstructSDandField()
}
if (fSetSDs) {
+ // set verbosity
+ G4SDManager::GetSDMpointer()->SetVerboseLevel(0);
+
// register sensitive detectors with SD manager
G4SDManager* SDman = G4SDManager::GetSDMpointer();
SDman -> AddNewDetector(fSiPM_SD.Get());
SDman -> AddNewDetector(fBetaDetSignal_SD.Get());
SDman -> AddNewDetector(fBetaDetVeto_SD.Get());
- // register sensitive detectors to logical volumes
- SetSensitiveDetector(fBetaDetector->GetPhotocathodeLogVol(), fSiPM_SD.Get());
- SetSensitiveDetector(fBetaDetector->GetSignalLogVol(), fBetaDetSignal_SD.Get());
- SetSensitiveDetector(fBetaDetector->GetVetoLogVol(), fBetaDetVeto_SD.Get());
+
+ // Only add detectors once
fSetSDs = false;
}
+
+ // register sensitive detectors to logical volumes
+ // (repeat every time geometry is reinitialized)
+ SetSensitiveDetector(fBetaDetector->GetPhotocathodeLogVol(), fSiPM_SD.Get());
+ SetSensitiveDetector(fBetaDetector->GetSignalLogVol(), fBetaDetSignal_SD.Get());
+ SetSensitiveDetector(fBetaDetector->GetVetoLogVol(), fBetaDetVeto_SD.Get());
}
}
diff --git a/src/VKEventAction.cpp b/src/VKEventAction.cpp
index 59f5cb377d6b086cbdf293391a4c066b317f52f3..96bc38babaf6b96e575d7c587ffc296896f914f0 100644
--- a/src/VKEventAction.cpp
+++ b/src/VKEventAction.cpp
@@ -31,6 +31,8 @@
#include "VKTrajectory.hh"
#include "VKRecorderBase.hh"
+#include "VKScintillatorSensitiveDetector.hh"
+
#include "G4EventManager.hh"
#include "G4SDManager.hh"
#include "G4RunManager.hh"
@@ -43,17 +45,22 @@
#include "G4UImanager.hh"
#include "G4SystemOfUnits.hh"
+#include <vector>
+
+using namespace std;
+
VKEventAction::VKEventAction(VKRecorderBase* r)
: fRecorder(r),
fSaveThreshold(0),
fScintCollID(-1),
+ fVetoCollID(-1),
fSiPMCollID(-1),
fVerbose(0),
fSiPMThreshold(1),
fForcedrawphotons(false),
fForcenophotons(false),
- fPrintModulo(100)
+ fEvtPct(0)
{
fEventMessenger = new VKEventMessenger(this);
}
@@ -72,8 +79,9 @@ void VKEventAction::BeginOfEventAction(const G4Event* anEvent)
G4SDManager* SDman = G4SDManager::GetSDMpointer();
if (fScintCollID < 0) fScintCollID = SDman -> GetCollectionID("Signal_SD_HitCollection");
+ if (fVetoCollID < 0) fVetoCollID = SDman -> GetCollectionID("Veto_SD_HitCollection");
if (fSiPMCollID < 0) fSiPMCollID = SDman -> GetCollectionID("SiPM_SD_HitCollection");
-
+
if (fRecorder) fRecorder -> RecordBeginOfEvent(anEvent);
}
@@ -81,7 +89,7 @@ void VKEventAction::BeginOfEventAction(const G4Event* anEvent)
void VKEventAction::EndOfEventAction(const G4Event* anEvent)
{
// event info
- VKUserEventInformation* eventInformation
+ VKUserEventInformation* eventInfo
= (VKUserEventInformation*)anEvent -> GetUserInformation();
// trajectories
@@ -103,54 +111,63 @@ void VKEventAction::EndOfEventAction(const G4Event* anEvent)
}
}
- VKScintillatorHitsCollection* scintHC = 0;
+ VKScintillatorHitsCollection * scintHC = 0;
+ VKScintillatorHitsCollection * vetoHC = 0;
VKSiPMHitsCollection* sipmHC = 0;
G4HCofThisEvent* hitsCE = anEvent -> GetHCofThisEvent();
// Get hit collections
if (hitsCE) {
if (fScintCollID >= 0) scintHC = (VKScintillatorHitsCollection*)(hitsCE->GetHC(fScintCollID));
+ if (fVetoCollID >= 0) vetoHC = (VKScintillatorHitsCollection*)(hitsCE->GetHC(fVetoCollID));
if (fSiPMCollID >= 0) sipmHC = (VKSiPMHitsCollection*)(hitsCE->GetHC(fSiPMCollID));
}
- // Hits in scintillator
- if (scintHC) {
- int n_hit = scintHC -> entries(); // number of hits
- G4ThreeVector eWeightPos(0.); // energy-weighted position of hit
- G4double edep;
- G4double edepMax = 0;
-
- //gather info on hits in scintillator
- for (int i=0; i<n_hit; i++) {
- edep = (*scintHC)[i] -> GetEdep();
- eventInformation -> IncEDep(edep); // add up the energy deposition
- eWeightPos += (*scintHC)[i]->GetPos() * edep; // calculate energy-weighted position
- if (edep > edepMax) {
- edepMax = edep; // store max energy deposition
- G4ThreeVector posMax = (*scintHC)[i]->GetPos();
- eventInformation -> SetPosMax(posMax, edep);
+ // Hits in scintillators (signal and veto)
+ vector<VKScintillatorHitsCollection*> scintHCs = {scintHC, vetoHC};
+ G4int id = 0;
+ for (auto& hc : scintHCs) {
+ if (hc) {
+ int n_hit = hc -> entries(); // number of hits
+ G4ThreeVector eWeightPos(0.); // energy-weighted position of hit
+ G4double edep;
+ G4double edepMax = 0;
+
+ //gather info on hits in scintillator
+ for (int i=0; i<n_hit; i++) {
+ edep = (*hc)[i] -> GetEdep();
+ eventInfo -> IncEDep(id, edep); // add up the energy deposition
+ eWeightPos += (*hc)[i]->GetPos() * edep; // calculate energy-weighted position
+ if (edep > edepMax) {
+ edepMax = edep; // store max energy deposition
+ G4ThreeVector posMax = (*hc)[i]->GetPos();
+ eventInfo -> SetPosMax(id, posMax, edep);
+ }
}
- }
- if (eventInformation->GetEDep() == 0.) {
- if (fVerbose>0) G4cout << "No hits in the scintillator this event." << G4endl;
- }
- else{
- // Finish calculation of energy-weighted position
- eWeightPos /= eventInformation->GetEDep();
- eventInformation -> SetEWeightPos(eWeightPos);
+ if (eventInfo->GetEDep(id) == 0.) {
+ if (fVerbose>0) G4cout << "No hits in " << hc->GetName() << " this event." << G4endl;
+ }
+ else{
+ // Finish calculation of energy-weighted position
+ eWeightPos /= eventInfo->GetEDep(id);
+ eventInfo -> SetEWeightPos(id, eWeightPos);
+ if (fVerbose > 0) {
+ G4cout << "\tEnergy weighted position of hits in " << hc->GetName() << " : "
+ << eWeightPos/mm << G4endl;
+ }
+ }
+
+ // print total energy deposition
if (fVerbose > 0) {
- G4cout << "\tEnergy weighted position of hits in scintillator : "
- << eWeightPos/mm << G4endl;
+ G4cout << "\tTotal energy deposition in " << hc->GetName() << " : "
+ << eventInfo->GetEDep(id) / keV << " (keV)" << G4endl;
}
}
-
- // print total energy deposition
- if (fVerbose > 0) {
- G4cout << "\tTotal energy deposition in scintillator : "
- << eventInformation->GetEDep() / keV << " (keV)" << G4endl;
- }
+ id++;
}
+ G4String scintName[2] = {"Signal", "Veto"};
+
// Hits in SiPM
if (sipmHC) {
G4ThreeVector reconPos(0.,0.,0.);
@@ -158,55 +175,62 @@ void VKEventAction::EndOfEventAction(const G4Event* anEvent)
//Gather info from all SiPMs
for (G4int i=0; i<sipms; i++) {
+ G4int id = (*sipmHC)[i]->GetSiPMNumber();
+ G4String name = (id < 2)? scintName[id] : "Unknown";
+ if(fVerbose > 0) {
+ G4cout << "\t======= SiPM data for: " << name << " =======" << G4endl;
+ }
+
G4int photons = (*sipmHC)[i]->GetPhotonCount();
- eventInformation -> IncHitCount(photons);
+ eventInfo -> IncHitCount(id, photons);
reconPos += (*sipmHC)[i]->GetSiPMPos() * photons;
if (photons >= fSiPMThreshold) {
- eventInformation -> IncSiPMsAboveThreshold();
+ eventInfo -> IncSiPMsAboveThreshold(id);
}
else { // if it wasnt above the threshold, turn it back off
(*sipmHC)[i]->SetDrawit(false);
}
- }
- // if there were any hits, determine position
- if (eventInformation->GetHitCount() > 0) {
- reconPos /= eventInformation->GetHitCount();
- if(fVerbose > 0) {
- G4cout << "\tReconstructed position of hits in SiPM : "
- << reconPos/mm << G4endl;
+ // if there were any hits, determine position
+ if (eventInfo->GetHitCount(id) > 0) {
+ reconPos /= eventInfo->GetHitCount(id);
+ if(fVerbose > 0) {
+ G4cout << "\tReconstructed position of hits in SiPMs: "
+ << reconPos/mm << G4endl;
+ }
+ eventInfo -> SetReconPos(id, reconPos);
+ }
+
+ //End of event output. later to be controlled by a verbose level
+ if (fVerbose > 0) {
+ G4cout << "\tNumber of photons that hit SiPMs in this event : "
+ << eventInfo->GetHitCount(id) << G4endl;
+ G4cout << "\tNumber of SiPMs above threshold ("<<fSiPMThreshold<<") : "
+ << eventInfo->GetSiPMsAboveThreshold(id) << G4endl;
+ G4cout << "\tNumber of photons produced by scintillation in this event : "
+ << eventInfo->GetPhotonCount_Scint(id) << G4endl;
+ G4cout << "\tNumber of photons produced by cerenkov in this event : "
+ << eventInfo->GetPhotonCount_Ceren(id) << G4endl;
+ G4cout << "\tNumber of photons absorbed (OpAbsorption) in this event : "
+ << eventInfo->GetAbsorptionCount(id) << G4endl;
+ G4cout << "\tNumber of photons absorbed at boundaries (OpBoundary) in "
+ << "this event : " << eventInfo->GetBoundaryAbsorptionCount(id)
+ << G4endl;
+ G4cout << "Unacounted for photons in this event : "
+ << (eventInfo->GetPhotonCount_Scint(id) +
+ eventInfo->GetPhotonCount_Ceren(id) -
+ eventInfo->GetAbsorptionCount(id) -
+ eventInfo->GetHitCount(id) -
+ eventInfo->GetBoundaryAbsorptionCount(id))
+ << G4endl;
}
- eventInformation -> SetReconPos(reconPos);
}
- sipmHC -> DrawAllHits();
- }
- //End of event output. later to be controlled by a verbose level
- if (fVerbose > 0) {
- G4cout << "\tNumber of photons that hit SiPMs in this event : "
- << eventInformation->GetHitCount() << G4endl;
- G4cout << "\tNumber of SiPMs above threshold ("<<fSiPMThreshold<<") : "
- << eventInformation->GetSiPMsAboveThreshold() << G4endl;
- G4cout << "\tNumber of photons produced by scintillation in this event : "
- << eventInformation->GetPhotonCount_Scint() << G4endl;
- G4cout << "\tNumber of photons produced by cerenkov in this event : "
- << eventInformation->GetPhotonCount_Ceren() << G4endl;
- G4cout << "\tNumber of photons absorbed (OpAbsorption) in this event : "
- << eventInformation->GetAbsorptionCount() << G4endl;
- G4cout << "\tNumber of photons absorbed at boundaries (OpBoundary) in "
- << "this event : " << eventInformation->GetBoundaryAbsorptionCount()
- << G4endl;
- G4cout << "Unacounted for photons in this event : "
- << (eventInformation->GetPhotonCount_Scint() +
- eventInformation->GetPhotonCount_Ceren() -
- eventInformation->GetAbsorptionCount() -
- eventInformation->GetHitCount() -
- eventInformation->GetBoundaryAbsorptionCount())
- << G4endl;
+ sipmHC -> DrawAllHits();
}
// If we have set the flag to save 'special' events, save here
- if (fSaveThreshold&&eventInformation->GetPhotonCount() <= fSaveThreshold)
+ if (fSaveThreshold && eventInfo->GetPhotonCount() <= fSaveThreshold)
G4RunManager::GetRunManager()->rndmSaveThisEvent();
// Print progress
@@ -233,16 +257,19 @@ void VKEventAction::SetSaveThreshold(G4int save)
void VKEventAction::PrintProgress()
{
// number of events to be processed
- G4int evtNoMax = G4RunManager::GetRunManager()
+ double evtMax = G4RunManager::GetRunManager()
->GetCurrentRun()->GetNumberOfEventToBeProcessed();
// current event no.
- G4int evtNo = G4RunManager::GetRunManager()
+ double evt = G4RunManager::GetRunManager()
->GetCurrentEvent()->GetEventID();
- if (evtNo % fPrintModulo == 0 && evtNo > 0)
+ if (evt == 0) fEvtPct = 0;
+
+ G4int evtPct = (evt + 0.5) / evtMax * 100;
+ if (evtPct > fEvtPct && evt > 0)
{
- float f = (float)(evtNo)/(float)evtNoMax *100.;
- printf(" %i%% of %i events simulated ...", (int)f, evtNoMax);
+ fEvtPct = evtPct;
+ printf(" %i%% of %i events simulated ...", evtPct, (int)evtMax);
fflush(stdout);
printf(" \r");
}
diff --git a/src/VKGeneralPhysics.cpp b/src/VKGeneralPhysics.cpp
index 21f44cfa580941ce75cd08fc726a45ab0ffd3c06..5b4e4a951f6821cae5f65364543bce7db920d020 100644
--- a/src/VKGeneralPhysics.cpp
+++ b/src/VKGeneralPhysics.cpp
@@ -29,7 +29,7 @@
#include "G4ios.hh"
#include <iomanip>
#include "G4Decay.hh"
-
+#include <G4Proton.hh>
VKGeneralPhysics::VKGeneralPhysics(const G4String& name)
: G4VPhysicsConstructor(name)
@@ -61,6 +61,9 @@ void VKGeneralPhysics::ConstructParticle()
G4Geantino::GeantinoDefinition();
G4ChargedGeantino::ChargedGeantinoDefinition();
*/
+ // proton
+ G4Proton::ProtonDefinition();
+
// generic ion
G4GenericIon::GenericIonDefinition();
}
diff --git a/src/VKMagnet.cpp b/src/VKMagnet.cpp
index 8ae5aa0aa8b4abdcaf583296f06adfbb80f34e7f..6025a16d0b1a50579b8c68466a32a460477b564e 100644
--- a/src/VKMagnet.cpp
+++ b/src/VKMagnet.cpp
@@ -64,10 +64,12 @@ VKMagnet::VKMagnet(G4RotationMatrix *pRot,
G4Material* copper = fMaterials -> GetMaterial("COPPER");
G4Material* iron = fMaterials -> GetMaterial("IRON");
G4Material* brass = fMaterials -> GetMaterial("BRASS");
+ G4Material* pmma = fMaterials -> GetMaterial("G4_PLEXIGLASS");
assert(copper != nullptr); // check that material exists
assert(iron != nullptr);
assert(brass != nullptr);
assert(labVacuum != nullptr);
+ assert(pmma != nullptr);
// rotation matrices, translation and transformation
G4RotationMatrix doNotRotate, rotate180DegAroundY;
@@ -195,6 +197,19 @@ VKMagnet::VKMagnet(G4RotationMatrix *pRot,
"brassWalls", fContainer_log, false, 0, fCheckOverlaps);
+/*-------------------------------------------------------------
+ extra_sol = new G4Tubs("extra_sol",
+ fBrassInnerDiameter/2 - 1.*cm,
+ fBrassInnerDiameter/2,
+ 3./2.*cm,
+ startAngle=0*deg,
+ spanningAngle=360*deg);
+ extra_log = new G4LogicalVolume(extra_sol, pmma, "extra_log");
+ new G4PVPlacement(0, G4ThreeVector(0,0,0), extra_log,
+ "extra", fContainer_log, false, 0, fCheckOverlaps);
+--------------------------------------------------------------*/
+
+
// visualization attributes
VisAttributes();
diff --git a/src/VKPhysicsList.cpp b/src/VKPhysicsList.cpp
index 2c54d6c7f8657aad873ea56cf8e48199f01a2e84..2e904f75989e0b59ee226ba15dd326031f9ebcab 100644
--- a/src/VKPhysicsList.cpp
+++ b/src/VKPhysicsList.cpp
@@ -29,6 +29,7 @@
#include "VKEMPhysics.hh"
#include "VKMuonPhysics.hh"
#include "VKPhysicsListMessenger.hh"
+#include "G4EmStandardPhysics_option4.hh"
#include "G4EmLowEPPhysics.hh"
#include "G4EmLivermorePhysics.hh"
@@ -50,13 +51,14 @@ VKPhysicsList::VKPhysicsList()
// default cut value (1.0mm)
defaultCutValue = 1.0*mm;
+ SetDefaultCutValue(defaultCutValue);
// General Physics
RegisterPhysics(new VKGeneralPhysics("general"));
// EM Physics
- RegisterPhysics(new VKEMPhysics("standard EM"));
- // RegisterPhysics(new G4EmStandardPhysics_option4()); // not yet tested
+ // RegisterPhysics(new VKEMPhysics("standard EM"));
+ RegisterPhysics(new G4EmStandardPhysics_option4());
// RegisterPhysics(new G4EmLivermorePhysics()); // overpredicts muon energy loss by a factor of 2!
// RegisterPhysics(new G4EmLowEPPhysics());
// RegisterPhysics(new G4EmPenelopePhysics()); // not yet tested
diff --git a/src/VKRecorderRoot.cpp b/src/VKRecorderRoot.cpp
index 72035b429681fa82d7db8537fb9537f3c522eb5f..8cb1e113f4eda56580736f92e3c326dcd3fec819 100644
--- a/src/VKRecorderRoot.cpp
+++ b/src/VKRecorderRoot.cpp
@@ -29,9 +29,13 @@
#include "VKBetaDetector.hh"
#include "VKRecorderRootMessenger.hh"
#include "VKPrimaryGeneratorAction.hh"
+#include "VKUserEventInformation.hh"
+#include "VKScintillatorHit.hh"
+#include "VKSiPMHit.hh"
#include "G4VPhysicalVolume.hh"
#include "G4RunManager.hh"
+#include "G4EventManager.hh"
#include <TNamed.h>
#include <TParameter.h>
@@ -146,9 +150,22 @@ void VKRecorderRoot::RecordBeginOfEvent(const G4Event* anEvent)
}
-void VKRecorderRoot::RecordEndOfEvent(const G4Event*)
+void VKRecorderRoot::RecordEndOfEvent(const G4Event* anEvent)
{
- WriteEventToTree();
+ // event info
+ VKUserEventInformation* eventInfo
+ = (VKUserEventInformation*)anEvent -> GetUserInformation();
+
+ // deposited energy
+ fEDepSignal = eventInfo->GetEDep(0);
+ fEDepVeto = eventInfo->GetEDep(1);
+
+ // detected optical photons
+ fPhotonsSignal = eventInfo->GetHitCount(0);
+ fPhotonsVeto = eventInfo->GetHitCount(1);
+
+ // fill ROOT tree
+ FillTree();
}
@@ -206,17 +223,6 @@ void VKRecorderRoot::RecordStep(const G4Step* aStep)
}
}
}
-
-
- // ------------------------------------------ //
- // ENERGY DEPOSITION IN SIGNAL AND VETO //
- // ------------------------------------------ //
-
- if (fBetaOn) {
- G4double EdepStep = aStep -> GetTotalEnergyDeposit();
- if (prePV == fSignalPV) fEDepSignal += EdepStep;
- if (prePV == fVetoPV) fEDepVeto += EdepStep;
- }
}
@@ -261,21 +267,24 @@ void VKRecorderRoot::CreateTraBranches()
}
-void VKRecorderRoot::WriteEventToTree()
+void VKRecorderRoot::FillTree()
{
if (fTFile == NULL) return;
// zero suppression
G4bool skip = fZeroSuppression && fEDepSignal==0 && fEDepVeto==0;
+
// beta detector
- _eDepSignal = fEDepSignal /keV;
- _eDepVeto = fEDepVeto /keV;
- _photonsSignal = 0;
- _photonsVeto = 0;
+ _eDepSignal = fEDepSignal /keV;
+ _eDepVeto = fEDepVeto /keV;
+ _photonsSignal = fPhotonsSignal;
+ _photonsVeto = fPhotonsVeto;
+
if (fBetaOn && !skip) fDetTree -> Fill();
+
// virtual detectors
// focal plane
@@ -285,6 +294,7 @@ void VKRecorderRoot::WriteEventToTree()
_Efp = fEkin_fp /keV;
_THfp = fTheta_fp /deg;
_PHIfp = fPhi_fp /deg;
+
// midway plane
_Xmp = fX_mp /mm;
_Ymp = fY_mp /mm;
@@ -292,6 +302,7 @@ void VKRecorderRoot::WriteEventToTree()
_Emp = fEkin_mp /keV;
_THmp = fTheta_mp /deg;
_PHImp = fPhi_mp /deg;
+
// source
_Xso = fX_so /mm;
_Yso = fY_so /mm;
@@ -299,8 +310,8 @@ void VKRecorderRoot::WriteEventToTree()
_Eso = fEkin_so /keV;
_THso = fTheta_so /deg;
_PHIso = fPhi_so /deg;
- // write to tree
- if (fVirtualOn && !skip) fTraTree -> Fill();
+
+ if (fVirtualOn) fTraTree -> Fill();
}
@@ -308,6 +319,7 @@ void VKRecorderRoot::ClearRootBranches()
{
_eDepSignal = _eDepVeto = 0;
_photonsSignal = _photonsVeto = 0;
+
_Xfp = _Yfp = _Zfp = _Efp = _THfp = _PHIfp = 0;
_Xmp = _Ymp = _Zmp = _Emp = _THmp = _PHImp = 0;
_Xso = _Yso = _Zso = _Eso = _THso = _PHIso = 0;
diff --git a/src/VKRunAction.cpp b/src/VKRunAction.cpp
index e3658a0483003c6918205942a7a0a866b366bcad..e2ee23233d15775a31e65ec04afe2f4a3f61029f 100644
--- a/src/VKRunAction.cpp
+++ b/src/VKRunAction.cpp
@@ -65,7 +65,8 @@ void VKRunAction::EndOfRunAction(const G4Run* aRun)
{
// run duration
fStopwatch -> Stop();
- G4cout << " ---> End run " << G4endl;
+ G4cout << " ---> End run" << G4endl;
+ G4cout << " ---> " << G4RunManager::GetRunManager()->GetCurrentRun()->GetNumberOfEventToBeProcessed() << " events simulated" << G4endl;
G4cout << " ---> Run duration: " << Form("%.2f", fStopwatch->RealTime()) << " seconds";
double rate = aRun->GetNumberOfEventToBeProcessed()/(fStopwatch->RealTime()/3600);
G4cout << " (" << Form("%.2e", rate) << " events/hour)" << G4endl;
diff --git a/src/VKScintillatorSensitiveDetector.cpp b/src/VKScintillatorSensitiveDetector.cpp
index 44e166209c0abd351c9d3d61df686f931c27c457..142a564a880729f76fa30fabda0ad2d74d9efcad 100644
--- a/src/VKScintillatorSensitiveDetector.cpp
+++ b/src/VKScintillatorSensitiveDetector.cpp
@@ -36,10 +36,12 @@
#include "G4TouchableHistory.hh"
#include "G4ios.hh"
#include "G4VProcess.hh"
+#include "G4SDManager.hh"
VKScintillatorSensitiveDetector::VKScintillatorSensitiveDetector(G4String name)
- : G4VSensitiveDetector(name)
+ : G4VSensitiveDetector(name),
+ fHitsCID(-1)
{
fScintCollection = NULL;
collectionName.insert(name+"_HitCollection");
@@ -54,20 +56,23 @@ VKScintillatorSensitiveDetector::~VKScintillatorSensitiveDetector()
void VKScintillatorSensitiveDetector::Initialize(G4HCofThisEvent* hitsCE)
{
+ // create hits collection
fScintCollection = new VKScintillatorHitsCollection(SensitiveDetectorName, collectionName[0]);
- // A way to keep all the hits of this event in one place if needed
- static G4int hitsCID = -1;
- if (hitsCID<0) {
- hitsCID = GetCollectionID(0); // "CID" = Collection ID
+ // get collection ID
+ if (fHitsCID < 0) {
+ fHitsCID = G4SDManager::GetSDMpointer()->GetCollectionID(collectionName[0]);
}
- hitsCE -> AddHitsCollection(hitsCID, fScintCollection);
+
+ // register collection
+ hitsCE -> AddHitsCollection(fHitsCID, fScintCollection);
}
-G4bool VKScintillatorSensitiveDetector::ProcessHits(G4Step* aStep,G4TouchableHistory* )
+G4bool VKScintillatorSensitiveDetector::ProcessHits(G4Step* aStep, G4TouchableHistory* )
{
G4double edep = aStep -> GetTotalEnergyDeposit();
+
if (edep==0.) return false; // No energy deposited so do not count as hit
// pre-step point and volume
diff --git a/src/VKSiPMScintInfo.cpp b/src/VKSiPMScintInfo.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..080069ee9eb87a551b75abb827774d6f2b0683ce
--- /dev/null
+++ b/src/VKSiPMScintInfo.cpp
@@ -0,0 +1,52 @@
+//
+// ********************************************************************
+// * License and Disclaimer *
+// * *
+// * The Geant4 software is copyright of the Copyright Holders of *
+// * the Geant4 Collaboration. It is provided under the terms and *
+// * conditions of the Geant4 Software License, included in the file *
+// * LICENSE and available at http://cern.ch/geant4/license . These *
+// * include a list of copyright holders. *
+// * *
+// * Neither the authors of this software system, nor their employing *
+// * institutes,nor the agencies providing financial support for this *
+// * work make any representation or warranty, express or implied, *
+// * regarding this software system or assume any liability for its *
+// * use. Please see the license in the file LICENSE and URL above *
+// * for the full disclaimer and the limitation of liability. *
+// * *
+// * This code implementation is the result of the scientific and *
+// * technical work of the GEANT4 collaboration. *
+// * By using, copying, modifying or distributing the software (or *
+// * any work based on the software) you agree to acknowledge its *
+// * use in resulting scientific publications, and indicate your *
+// * acceptance of all terms of the Geant4 Software license. *
+// ********************************************************************
+//
+
+#include "VKSiPMScintInfo.hh"
+
+
+VKSiPMScintInfo::VKSiPMScintInfo()
+ : fHitCount(0),
+ fPhotonCount_Scint(0),
+ fPhotonCount_Ceren(0),
+ fAbsorptionCount(0),
+ fBoundaryAbsorptionCount(0),
+ fTotE(0.),
+ fEWeightPos(0.),
+ fReconPos(0.),
+ fConvPos(0.),
+ fConvPosSet(false),
+ fPosMax(0.),
+ fEdepMax(0.),
+ fSiPMsAboveThreshold(0)
+{
+ // do nothing
+}
+
+
+VKSiPMScintInfo::~VKSiPMScintInfo()
+{
+ // do nothing
+}
diff --git a/src/VKSiPMSensitiveDetector.cpp b/src/VKSiPMSensitiveDetector.cpp
index 603fdba9e551e8f1fedcd4588101367227e8da9d..a0aa24007cf185c9fad840571efbfd29e296cb3b 100644
--- a/src/VKSiPMSensitiveDetector.cpp
+++ b/src/VKSiPMSensitiveDetector.cpp
@@ -39,9 +39,12 @@
#include "G4ParticleTypes.hh"
#include "G4ParticleDefinition.hh"
+#include <iostream> // std::cin, std::cout
+
VKSiPMSensitiveDetector::VKSiPMSensitiveDetector(G4String name)
: G4VSensitiveDetector(name),
+ fHitsCID(-1),
fSiPMHitCollection(0),
fSiPMPositionsX(0),
fSiPMPositionsY(0),
@@ -70,12 +73,14 @@ void VKSiPMSensitiveDetector::SetSiPMPositions(const std::vector<G4ThreeVector>&
void VKSiPMSensitiveDetector::Initialize(G4HCofThisEvent* hitsCE)
{
fSiPMHitCollection = new VKSiPMHitsCollection(SensitiveDetectorName, collectionName[0]);
- // Store collection with event and keep ID
- static G4int hitCID = -1;
- if (hitCID<0) {
- hitCID = GetCollectionID(0);
+
+ // get collection ID
+ if (fHitsCID < 0) {
+ fHitsCID = G4SDManager::GetSDMpointer()->GetCollectionID(collectionName[0]);
}
- hitsCE -> AddHitsCollection(hitCID, fSiPMHitCollection );
+
+ // register collection
+ hitsCE -> AddHitsCollection(fHitsCID, fSiPMHitCollection);
}
@@ -91,12 +96,11 @@ G4bool VKSiPMSensitiveDetector::ProcessHits(G4Step* ,G4TouchableHistory* )
G4bool VKSiPMSensitiveDetector::ProcessHits_constStep(const G4Step* aStep, G4TouchableHistory*)
{
-
// Only do something if this is an optical photon
G4ParticleDefinition * partDef = aStep->GetTrack()->GetDefinition();
if (partDef != G4OpticalPhoton::OpticalPhotonDefinition()) return false;
- // User replica number 1 since photocathode is a daughter volume <--- HMMM... IS THIS STILL THE CASE ??!!
+ // User replica number 1 since photocathode is a daughter volume
// to the SiPM which was replicated
G4int SiPMNumber =
aStep->GetPostStepPoint()->GetTouchable()->GetReplicaNumber(1);
diff --git a/src/VKStackingAction.cpp b/src/VKStackingAction.cpp
index 3dc02b5cfc9441c6afb28401638702946b3533a5..61a7f8af2fd8817414090edd1923506278a3a843 100644
--- a/src/VKStackingAction.cpp
+++ b/src/VKStackingAction.cpp
@@ -52,7 +52,7 @@ VKStackingAction::~VKStackingAction()
G4ClassificationOfNewTrack VKStackingAction::ClassifyNewTrack(const G4Track * aTrack)
{
- VKUserEventInformation* eventInformation =
+ VKUserEventInformation* eventInfo =
(VKUserEventInformation*)G4EventManager::GetEventManager()
-> GetConstCurrentEvent() -> GetUserInformation();
@@ -61,10 +61,15 @@ G4ClassificationOfNewTrack VKStackingAction::ClassifyNewTrack(const G4Track * aT
{
if (aTrack->GetParentID() > 0) // particle is secondary
{
+ G4int scintId = 0;
+ auto name = aTrack->GetVolume()->GetName();
+ if (name == "signal") scintId = 0;
+ else if (name == "veto") scintId = 1;
+
if (aTrack->GetCreatorProcess()->GetProcessName() == "Scintillation")
- eventInformation->IncPhotonCount_Scint();
+ eventInfo->IncPhotonCount_Scint(scintId);
else if (aTrack->GetCreatorProcess()->GetProcessName() == "Cerenkov")
- eventInformation->IncPhotonCount_Ceren();
+ eventInfo->IncPhotonCount_Ceren(scintId);
}
}
diff --git a/src/VKSteppingAction.cpp b/src/VKSteppingAction.cpp
index 7bec1e35f614f072cc7a38607232faebfade2050..44ef59487d01aa20353ff476fd7ed9eafc859f76 100644
--- a/src/VKSteppingAction.cpp
+++ b/src/VKSteppingAction.cpp
@@ -34,6 +34,8 @@
#include "VKSteppingMessenger.hh"
#include "VKRecorderBase.hh"
+#include "VKScintillatorSensitiveDetector.hh"
+
#include "G4SteppingManager.hh"
#include "G4SDManager.hh"
#include "G4EventManager.hh"
@@ -70,18 +72,18 @@ void VKSteppingAction::UserSteppingAction(const G4Step * theStep)
if (theTrack->GetCurrentStepNumber() == 1) fExpectedNextStatus = Undefined;
// track info
- VKUserTrackInformation* trackInformation
+ VKUserTrackInformation* trackInfo
= (VKUserTrackInformation*)theTrack -> GetUserInformation();
// event info
- VKUserEventInformation* eventInformation
+ VKUserEventInformation* eventInfo
= (VKUserEventInformation*)G4EventManager::GetEventManager()
-> GetConstCurrentEvent() -> GetUserInformation();
// pre-step info
G4StepPoint* thePrePoint = theStep -> GetPreStepPoint();
G4VPhysicalVolume* thePrePV = thePrePoint -> GetPhysicalVolume();
-
+
// post-step info
G4StepPoint* thePostPoint = theStep -> GetPostStepPoint();
G4VPhysicalVolume* thePostPV = thePostPoint -> GetPhysicalVolume();
@@ -120,7 +122,7 @@ void VKSteppingAction::UserSteppingAction(const G4Step * theStep)
// If we havent already found the conversion position and there were
// secondaries generated, then search for it
- if (!eventInformation->IsConvPosSet() && tN2ndariesTot>0 ) {
+ if (!eventInfo->IsConvPosSet() && tN2ndariesTot>0 ) {
size_t lp1min = (*fSecondary).size()-tN2ndariesTot;
size_t lp1max = (*fSecondary).size();
// loop over secondaries
@@ -133,13 +135,14 @@ void VKSteppingAction::UserSteppingAction(const G4Step * theStep)
if (creatorName=="phot"||creatorName=="compt"||creatorName=="conv") {
// since this is happening before the secondary is being tracked
// the Vertex position has not been set yet(set in initial step)
- eventInformation -> SetConvPos((*fSecondary)[lp1]->GetPosition());
+ eventInfo -> SetConvPos((*fSecondary)[lp1]->GetPosition());
}
}
}
}
- // *** What does this do !??
+ // If user has selected to only allow primaries to go one step
+ // before being killed, kill track now
if (fOneStepPrimaries && (thePrePV->GetName()=="signal" || thePrePV->GetName()=="veto")) {
theTrack -> SetTrackStatus(fStopAndKill);
}
@@ -158,14 +161,21 @@ void VKSteppingAction::UserSteppingAction(const G4Step * theStep)
// kill optical photons entering `world`
if (thePostPV->GetName() == "world") theTrack -> SetTrackStatus(fStopAndKill);
+ // Check which SiPM/Scintillator the optical photon is in
+ auto preName = thePrePV->GetName();
+ auto preCopyNo = thePrePV->GetCopyNo();
+ G4int SiPMScintId = 0;
+ if (preName == "signal" || (preName == "SiPM" && preCopyNo == 0)) SiPMScintId = 0;
+ else if (preName == "veto" || (preName == "SiPM" && preCopyNo == 1)) SiPMScintId = 1;
+
// If photon was absorbed by the absorption process
G4String procName = thePostPoint->GetProcessDefinedStep()->GetProcessName();
if (procName == "OpAbsorption") {
- eventInformation -> IncAbsorption();
- trackInformation -> AddTrackStatusFlag(absorbed);
+ eventInfo -> IncAbsorption(SiPMScintId);
+ trackInfo -> AddTrackStatusFlag(absorbed);
}
- boundaryStatus = boundary -> GetStatus();
+ boundaryStatus = boundary -> GetStatus();
// Check to see if the partcile was actually at a boundary
// Otherwise the boundary status may not be valid
@@ -185,8 +195,8 @@ void VKSteppingAction::UserSteppingAction(const G4Step * theStep)
// Depending on boundary status, take different course of action
switch (boundaryStatus) {
case Absorption:
- trackInformation->AddTrackStatusFlag(boundaryAbsorbed);
- eventInformation->IncBoundaryAbsorption();
+ trackInfo->AddTrackStatusFlag(boundaryAbsorbed);
+ eventInfo->IncBoundaryAbsorption(SiPMScintId);
break;
case Detection: {
//
@@ -197,10 +207,10 @@ void VKSteppingAction::UserSteppingAction(const G4Step * theStep)
// absorbed and status is Detection
//
G4SDManager* SDman = G4SDManager::GetSDMpointer();
- G4String sdName = "/VKDet/SiPMSD";
+ G4String sdName = "SiPM_SD";
VKSiPMSensitiveDetector* sipmSD = (VKSiPMSensitiveDetector*)SDman -> FindSensitiveDetector(sdName);
if (sipmSD) sipmSD -> ProcessHits_constStep(theStep, NULL);
- trackInformation -> AddTrackStatusFlag(hitSiPM);
+ trackInfo -> AddTrackStatusFlag(hitSiPM);
break;
}
case FresnelReflection:
@@ -209,7 +219,7 @@ void VKSteppingAction::UserSteppingAction(const G4Step * theStep)
case LobeReflection:
case SpikeReflection:
case BackScattering:
- trackInformation -> IncReflections();
+ trackInfo -> IncReflections(); // this gets executed for all types of reflections
fExpectedNextStatus = StepTooSmall;
break;
default:
@@ -223,3 +233,20 @@ void VKSteppingAction::UserSteppingAction(const G4Step * theStep)
// kill track if step number exceeds 10,000
if (theTrack->GetCurrentStepNumber() > 1E4) theTrack->SetTrackStatus(fStopAndKill);
}
+
+/*
+ * 1: Undefined
+ * 2: FresnelRefraction
+ * 3: FresnelReflection
+ * 4: TotalInternalReflection
+ * 5: LambertianReflection
+ * 6: LobeReflection
+ * 7: SpikeReflection
+ * 8: BackScattering
+ * 9: Absorption
+ * 10: Detection
+ * 11: NotAtBoundary
+ * 12: SameMaterial
+ * 13: StepTooSmall
+ * 14: NoRINDEX
+*/
diff --git a/src/VKUserEventInformation.cpp b/src/VKUserEventInformation.cpp
index 657289b0fd89bfd88587800a66e0248617f0cf55..d322e7c50e4c055731904674032f33e0178b8fd9 100644
--- a/src/VKUserEventInformation.cpp
+++ b/src/VKUserEventInformation.cpp
@@ -28,19 +28,8 @@
VKUserEventInformation::VKUserEventInformation()
- : fHitCount(0),
- fPhotonCount_Scint(0),
- fPhotonCount_Ceren(0),
- fAbsorptionCount(0),
- fBoundaryAbsorptionCount(0),
- fTotE(0.),
- fEWeightPos(0.),
- fReconPos(0.),
- fConvPos(0.),
- fConvPosSet(false),
- fPosMax(0.),
- fEdepMax(0.),
- fSiPMsAboveThreshold(0)
+ : fConvPos(0.),
+ fConvPosSet(false)
{
// do nothing
}
@@ -48,5 +37,163 @@ VKUserEventInformation::VKUserEventInformation()
VKUserEventInformation::~VKUserEventInformation()
{
- // do nothing
+ for (auto& i : info) delete i;
+}
+
+
+void VKUserEventInformation::IncPhotonCount_Scint(size_t id)
+{
+ EnsureCapacity(id);
+ info[id]->IncPhotonCount_Scint();
+}
+
+void VKUserEventInformation::IncPhotonCount_Ceren(size_t id)
+{
+ EnsureCapacity(id);
+ info[id]->IncPhotonCount_Ceren();
+}
+
+void VKUserEventInformation::IncEDep(size_t id, G4double dep)
+{
+ EnsureCapacity(id);
+ info[id]->IncEDep(dep);
+}
+
+void VKUserEventInformation::IncAbsorption(size_t id)
+{
+ EnsureCapacity(id);
+ info[id]->IncAbsorption();
+}
+
+void VKUserEventInformation::IncBoundaryAbsorption(size_t id)
+{
+ EnsureCapacity(id);
+ info[id]->IncBoundaryAbsorption();
+}
+
+void VKUserEventInformation::IncHitCount(size_t id, G4int i)
+{
+ EnsureCapacity(id);
+ info[id]->IncHitCount(i);
+}
+
+void VKUserEventInformation::SetEWeightPos(size_t id, const G4ThreeVector& p)
+{
+ EnsureCapacity(id);
+ info[id]->SetEWeightPos(p);
+}
+
+void VKUserEventInformation::SetReconPos(size_t id, const G4ThreeVector& p)
+{
+ EnsureCapacity(id);
+ info[id]->SetReconPos(p);
+}
+
+void VKUserEventInformation::SetPosMax(size_t id, const G4ThreeVector& p, G4double edep)
+{
+ EnsureCapacity(id);
+ info[id]->SetPosMax(p, edep);
+}
+
+void VKUserEventInformation::IncSiPMsAboveThreshold(size_t id)
+{
+ EnsureCapacity(id);
+ info[id]->IncSiPMsAboveThreshold();
+}
+
+G4int VKUserEventInformation::GetPhotonCount_Scint(size_t id)
+{
+ EnsureCapacity(id);
+ return info[id]->GetPhotonCount_Scint();
+}
+
+G4int VKUserEventInformation::GetPhotonCount_Ceren(size_t id)
+{
+ EnsureCapacity(id);
+ return info[id]->GetPhotonCount_Ceren();
+}
+
+G4int VKUserEventInformation::GetHitCount(size_t id)
+{
+ EnsureCapacity(id);
+ return info[id]->GetHitCount();
+}
+
+G4double VKUserEventInformation::GetEDep(size_t id)
+{
+ EnsureCapacity(id);
+ return info[id]->GetEDep();
+}
+
+G4int VKUserEventInformation::GetAbsorptionCount(size_t id)
+{
+ EnsureCapacity(id);
+ return info[id]->GetAbsorptionCount();
+}
+
+G4int VKUserEventInformation::GetBoundaryAbsorptionCount(size_t id)
+{
+ EnsureCapacity(id);
+ return info[id]->GetBoundaryAbsorptionCount();
+}
+
+G4ThreeVector VKUserEventInformation::GetEWeightPos(size_t id)
+{
+ EnsureCapacity(id);
+ return info[id]->GetEWeightPos();
+}
+
+G4ThreeVector VKUserEventInformation::GetReconPos(size_t id)
+{
+ EnsureCapacity(id);
+ return info[id]->GetReconPos();
+}
+
+G4ThreeVector VKUserEventInformation::GetPosMax(size_t id)
+{
+ EnsureCapacity(id);
+ return info[id]->GetPosMax();
+}
+
+G4double VKUserEventInformation::GetEDepMax(size_t id)
+{
+ EnsureCapacity(id);
+ return info[id]->GetEDepMax();
+}
+
+G4int VKUserEventInformation::GetPhotonCount(size_t id)
+{
+ EnsureCapacity(id);
+ return info[id]->GetPhotonCount();
+}
+
+G4int VKUserEventInformation::GetPhotonCount()
+{
+ G4int photons = 0;
+ for (auto& i : info) photons += i->GetPhotonCount();
+ return photons;
+}
+
+G4int VKUserEventInformation::GetSiPMsAboveThreshold(size_t id)
+{
+ EnsureCapacity(id);
+ return info[id]->GetSiPMsAboveThreshold();
+}
+
+void VKUserEventInformation::EnsureCapacity(size_t id)
+{
+ if (id >= info.size()) info.emplace_back(new VKSiPMScintInfo());
+}
+
+void VKUserEventInformation::CheckCapacity(size_t id) const
+{
+ if (id >= info.size())
+ {
+ G4ExceptionDescription ed;
+ ed << "VKUserEventInformation::CheckCapacity(size_t): "
+ << "Attempt to access information for non-existing SiPM-Scintillator unit (id: " << id << ")"
+ << G4endl;
+ G4Exception("VKUserEventInformation::CheckCapacity(size_t):", "VK",
+ FatalException, ed, "Something is wrong with number of SiPM-Scintillator unit");
+ }
}