diff --git a/analysis/VKAdd b/analysis/VKAdd
index d66d47a8cb2877bef4af59c662b4663a8945693e..260f6faa9c359584481ef64f9f129254a80a57e4 100755
Binary files a/analysis/VKAdd and b/analysis/VKAdd differ
diff --git a/analysis/VKAnalysis.dox b/analysis/VKAnalysis.dox
index 69bdf4683c4c6d5302eb5d3e1f53b51bffa49863..276ac49501b062d0a98692434cca60d9cebf3bae 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/VeikonKone/analysis/include/
+INPUT = /home/oliskir/Desktop/f20/sim/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/loop.sh b/analysis/loop.sh
index 78c1085d556fa1e9995bc4f3c3004353a64ee451..3fa9c01944d88caa6a8979740e4b1eafee22e449 100755
--- a/analysis/loop.sh
+++ b/analysis/loop.sh
@@ -1,6 +1,6 @@
-./VKAdd -o output/bi207_482_K.root --r ../output/bi207_482_K.root
-./VKAdd -o output/bi207_482_L.root --r ../output/bi207_482_L.root
-./VKAdd -o output/bi207_976_K.root --r ../output/bi207_976_K.root
-./VKAdd -o output/bi207_976_L.root --r ../output/bi207_976_L.root
-./VKAdd -o output/bi207_1682_K.root --r ../output/bi207_1682_K.root
-./VKAdd -o output/bi207_1682_L.root --r ../output/bi207_1682_L.root
+
+for b in 10 15 20 25 30 35 40 44 46 48 50
+do
+ ./VKAdd --r -o ../output/f20scan/200617/res/b$b.root ../output/f20scan/200617/b$b.root
+done
+
diff --git a/output/f20scan/BetaSpec.C b/output/f20scan/BetaSpec.C
deleted file mode 100644
index 1f3e3125080467fcd14602d6e879a1fbda566e8a..0000000000000000000000000000000000000000
--- a/output/f20scan/BetaSpec.C
+++ /dev/null
@@ -1,110 +0,0 @@
-
-/*
- * Simulated beta spectrum for various B-field settings with
- * experimental background spectrum superimposed.
- */
-
-void BetaSpec(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/output/f20scan/BetaSpecExp.C b/output/f20scan/BetaSpecExp.C
deleted file mode 100644
index 8def1b5fa663ca2190a32f3b65192d5e20072447..0000000000000000000000000000000000000000
--- a/output/f20scan/BetaSpecExp.C
+++ /dev/null
@@ -1,156 +0,0 @@
-
-/*
- * 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 BetaSpecExp(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/TransmissionExp.C b/output/f20scan/TransmissionExp.C
deleted file mode 100644
index b5b66b1afc9af3a38fb87aa5e7722b48b61ca51f..0000000000000000000000000000000000000000
--- a/output/f20scan/TransmissionExp.C
+++ /dev/null
@@ -1,610 +0,0 @@
-
-
-
-/*
- * Class that holds a double with uncertainty.
- */
-class DWU {
- public:
- DWU(double v = 0, double u = 0) : val(v), unc(u) {}
- virtual ~DWU(){}
-
- double val;
- double unc;
-
- // Assignment
-
- // =
- DWU& operator= (const DWU& d) {
- if (this != &d) { // check for self-assignment
- val = d.val;
- unc = d.unc;
- }
- return *this;
- }
-
- // Compound assignments
-
- // +=
- DWU& operator+=(const DWU& d) {
- double sum = val + d.val;
- double dsum = sqrt(pow(unc,2)+pow(d.unc,2));
- val = sum;
- unc = dsum;
- return *this;
- }
-
- // -=
- DWU& operator-=(const DWU& d) {
- double diff = val - d.val;
- double ddiff = sqrt(pow(unc,2)+pow(d.unc,2));
- val = diff;
- unc = ddiff;
- return *this;
- }
-
- // *=
- DWU& operator*=(const DWU& d) {
- double prod = val * d.val;
- double dprod = sqrt(pow(d.val*unc,2)+pow(val*d.unc,2));
- val = prod;
- unc = dprod;
- return *this;
- }
-
- // /=
- DWU& operator/=(const DWU& d) {
- double ratio = 0, dratio = 0;
- if (d.val != 0) {
- ratio = val / d.val;
- dratio = sqrt(pow(unc/d.val,2)+pow(val*d.unc/pow(d.val,2),2));
- }
- val = ratio;
- unc = dratio;
- return *this;
- }
-
- // Binary Arithmetic Operators
-
- // +
- DWU operator+ (const DWU& d) {
- DWU result = *this;
- result += d;
- return result;
- }
-
- // -
- DWU operator- (const DWU& d) {
- DWU result = *this;
- result -= d;
- return result;
- }
-
- // *
- DWU operator* (const DWU& d) {
- DWU result = *this;
- result *= d;
- return result;
- }
-
- // /
- DWU operator/ (const DWU& d) {
- DWU result = *this;
- result /= d;
- return result;
- }
-
- // scale
- void scale(const double& x) {
- val *= x;
- unc *= TMath::Abs(x);
- }
-
- // shift
- void shift(const double& x) {
- val += x;
- }
-
- private:
-};
-
-
-/*
- * Get duration in hours.
- */
-
-double GetDuration(TChain * chain)
-{
- // branches
- ULong64_t tms;
- chain->SetBranchAddress("tms", &tms);
-
- // find tmin and tmax
- int shifts = 0;
- double shift = 0;
- double Tp = 0, TFirst = 0, TLast = 0;
- for (int i=0; i<chain->GetEntries(); i++)
- {
- chain->GetEntry(i);
- double T = tms;
- if (Tp > T) {
- shift += Tp;
- shifts++;
- }
- if (i == 0) TFirst = T + shift;
- if (i == chain->GetEntries() - 1) TLast = T + shift;
- Tp = T;
- }
-
- double seconds = (TLast - TFirst) * 1E-3;
- double minutes = seconds / 60.;
- double hours = minutes / 60.;
-
- chain->ResetBranchAddresses();
-
- return hours;
-}
-
-/*
- * Locate and fit 1.6-MeV gamma peak
- */
-
-void FitGammaPeak(TH1D * h0, double& counts, double& e0, double& sigma, double& dcounts, double& de0, double& dsigma)
-{
- // make a copy
- TH1D * h = (TH1D*)h0->Clone("hcopy");
- h->GetXaxis()->SetRangeUser(1350, 1700);
-
- // get bin width
- double binWidth = h->GetBinWidth(1);
-
- // search for 1.6-MeV peak
- int b = h->GetMaximumBin();
- double xPos = h->GetBinCenter(b);
- double yPos = h->GetBinContent(b);
-
- // fit function
- TF1 * fitf = new TF1("fitf", "[0]*TMath::Gaus(x,[1],[2],true)+[3]", 0, 1E4);
-
- cout << " {Gamma peak found at: " << xPos << " keV}" << endl;
-
- // fit range
- double fitmin = xPos - 75;
- double fitmax = xPos + 75;
-
- // start values
- fitf->SetParameter(0, yPos);
- fitf->SetParameter(1, xPos);
- fitf->SetParameter(2, 20);
- fitf->SetParameter(3, 0);
-
- // fit
- h->Fit("fitf", "EBQ", "", fitmin, fitmax);
- fitmin = fitf->GetParameter(1) - 3*fitf->GetParameter(2);
- fitmax = fitf->GetParameter(1) + 3*fitf->GetParameter(2);
- h->Fit("fitf", "ES+BQ", "", fitmin, fitmax);
-
- // get fitted values
- counts = fitf->GetParameter(0) / binWidth;
- dcounts = fitf->GetParError(0) / binWidth;
- e0 = fitf->GetParameter(1);
- de0 = fitf->GetParError(1);
- sigma = fitf->GetParameter(2);
- dsigma = fitf->GetParError(2);
-}
-
-
-/*
- * Divide data into short runs to reduce gain drift.
- * For each run locate and fit 1.6-MeV gamma peak.
- * Finally, add up the counts from each.
- */
-
-void FitSlidingGammaPeak(int id, double minutes, TChain * chain, bool veto, double& counts, double& dcounts, double& sigma, double& x0min, double& x0max)
-{
- const double twindow = minutes * 60 * 1E3;
- double tstart = 0;
-
- ULong64_t entries = chain->GetEntries();
-
- // branches
- ULong64_t tms;
- chain->SetBranchAddress("tms", &tms);
- UInt_t M0, M1, M2;
- vector<double> E0(10), E1(10), E2(10);
- vector<UInt_t> Z0(10), Z1(10), Z2(10);
- chain->SetBranchAddress("M0", &M0);
- chain->SetBranchAddress("M1", &M1);
- chain->SetBranchAddress("M2", &M2);
- chain->SetBranchAddress("E0", E0.data());
- chain->SetBranchAddress("E1", E1.data());
- chain->SetBranchAddress("E2", E2.data());
- chain->SetBranchAddress("Z0", Z0.data());
- chain->SetBranchAddress("Z1", Z1.data());
- chain->SetBranchAddress("Z2", Z2.data());
-
- DWU nsum(0,0);
- DWU sigsum(0,0);
- x0min = 3000;
- x0max = 0;
-
- // histogram
- TString hname = Form("hgamma_%i", id);
- TH1D * h = new TH1D(hname, hname, 600, 0, 3000);
-
- // loop over data
- ULong64_t tms_prev = 0;
- ULong64_t toffset = 0;
- int fits = 0;
- for (int j = 0; j < entries; j++)
- {
- chain->GetEntry(j);
-
- // clock gets reset at every new file
- tms += toffset;
- if (tms < tms_prev) {
- tms -= toffset;
- toffset = tms_prev;
- tms += toffset;
- }
-
- if (tms >= tstart + twindow || j == entries-1) {
- // perform fit
- double n, dn, x0, dx0, sig, dsig;
- FitGammaPeak(h, n, x0, sig, dn, dx0, dsig);
- fits++;
-
- nsum += DWU(n, dn);
- sigsum += DWU(sig, dsig);
-
- x0min = TMath::Min(x0, x0min);
- x0max = TMath::Max(x0, x0max);
-
- // start of next time window
- tstart += twindow;
-
- // clear histogram
- h->Reset();
- }
-
- // fill histogram
- if (M2==1 && Z2[0]!=10) {
- bool vetoCut = true;
- if (veto) {
- vetoCut = M1==0 || (M1==1 && M0==1 && E1[0] < 180 + 0.21*E0[0]);
- }
- if (vetoCut) h->Fill(E2[0]);
- }
-
- tms_prev = tms;
- }
-
- counts = nsum.val;
- dcounts = nsum.unc;
-
- sigsum.scale(1. / (double)fits);
- sigma = sigsum.val;
-
- chain->ResetBranchAddresses();
-}
-
-
-/*
- * Plot transmission as a function of magnetic-field setting.
- */
-
-void TransmissionExp(bool veto = false, bool logy = false, double emin = 100, double emax = 5500)
-{
- // veto cut
- TCut vetoCut("");
- if (veto) vetoCut = TCut("M1==0 || (M1==1 && M0==1 && Sum$(E1) < 180 + 0.21*Sum$(E0))");
-
- // directories
- TString expdir = "/home/oliskir/Desktop/f20/data/sorted/new/";
- TString simdir = "080617/";
-
- // gamma detection efficiency
- const DWU gammaEff0(0.52E-4, 0.05E-4);
- vector<DWU> gammaEff;
-
- // how wide should the gate on the 1.6-MeV peak be? (expressed in units of sigma)
- // the gate is used to find beta-gamma coincidences.
- const double stdDevs = 2.5;
-
- // beta spectrum binning
- const int bins = (emax - emin) / 10;
-
- // number of magnetic field settings
- const int N = 10;
- vector<int> bmax = {15, 20, 25, 30, 35, 40, 44, 46, 48, 50};
-
- // runs at each setting
- vector<int> runs[N];
- runs[0] = {366};
- runs[1] = {332};
- runs[2] = {367};
- runs[3] = {333};
- runs[4] = {368};
- runs[5] = {334};
- runs[6] = {369};
- runs[7] = {360};
- runs[8] = {361, 362, 363, 364, 365, 370, 371, 372, 373, 375, 376, 379, 380};
- runs[9] = {358};
-
- // background run (337)
- int rBackgr = 337;
- TChain * cBackgr = new TChain("data");
- cBackgr->Add(expdir+Form("*0%i*", rBackgr));
- TH1F * hBackgr0 = new TH1F("hBackgr0", "hBackgr0", bins, emin, emax);
- cBackgr->Draw("E0 >> hBackgr0", "Z0 != 10" && vetoCut);
- double hoursBackgr = GetDuration(cBackgr);
-
- cout << endl;
- cout << " Background run:" << endl;
- cout << Form(" * Run: %i", rBackgr) << endl;
- cout << Form(" * Entries: %lli", cBackgr->GetEntries()) << endl;
- cout << Form(" * Duration: %.1f hours", hoursBackgr) << endl;
- cout << endl;
-
- // canvas
- TCanvas * c1 = new TCanvas("c1");
-
- // transmission
- vector<double> gammaxvec, transxvec;
- vector<DWU> transmission, transmissionSim;
-
- // beta spectrum integration limits
- int binx1 = hBackgr0->GetXaxis()->FindBin(emin);
- int binx2 = hBackgr0->GetXaxis()->FindBin(emax);
-
- // variables to store number of counts
- double n, dn;
-
- // loop over magnetic field settings
- int i = 0;
- for (auto& b : bmax)
- {
- // add files
- TChain * cExp = new TChain("data");
- for (auto& r : runs[i]) {
- TString fname = expdir + Form("*0%i*", r);
- cExp->Add(fname);
- }
-
- // duration of experimental run
- double hours = GetDuration(cExp);
-
- // integrate background spectrum
- TH1D * hBackgr = (TH1D*)hBackgr0->Clone("hBackgr");
- n = hBackgr->Integral(binx1, binx2);
- DWU backgr(n, sqrt(n));
-
- // scale
- backgr.scale(hours/hoursBackgr);
- hBackgr->Scale(hours/hoursBackgr);
-
- // singles histograms
- TString nb = Form("hb%i", i);
- TH1D * hb = new TH1D(nb, nb, bins, emin, emax);
- TString ng = Form("hg%i", i);
- TH1D * hg = new TH1D(ng, ng, 600, 0, 3000);
-
- // draw
- cExp->Draw("E0 >> " + nb, "M0==1 && Z0!=10" && vetoCut);
- cExp->Draw("E2 >> " + ng, "M2==1 && Z2!=10" && vetoCut);
-
- // integrate beta spectrum
- n = hb->Integral(binx1, binx2);
- DWU betas(n, sqrt(n));
-
- // subtract background
- betas -= backgr;
-
- // integrate 1.6-MeV peak in singles spectrum
- double centroid, dcentroid, sigma, dsigma;
-//*** FitGammaPeak(hg, n, centroid, sigma, dn, dcentroid, dsigma);
- double centroidMin = centroid;
- double centroidMax = centroid;
- double split = (int)(hours/0.5) + 1;
- double timeWindow = 1.01 * hours / split * 60;
- FitSlidingGammaPeak(i, timeWindow, cExp, veto, n, dn, sigma, centroidMin, centroidMax);
- centroid = (centroidMin + centroidMax) / 2;
- DWU gammas(n, dn);
-
- // number of decays
- DWU decays = gammas;
- decays.scale(1./gammaEff0.val);
-
- // transmission
- transmission.push_back(betas / decays);
- transxvec.push_back(141.0*1E-2*b);
-
- // coincidence gamma gate
- double egmin = centroidMin - stdDevs*sigma;
- double egmax = centroidMax + stdDevs*sigma;
-
- // coincidence histogram
- TString nbg = Form("hbg%i", i);
- TH1F * hbg = new TH1F(nbg, nbg, bins, emin, emax);
-
- // draw
- cExp->Draw("E0 >> " + nbg, Form("M0==1 && M2==1 && Z0 != 10 && Z2 != 10 && E2 >= %f && E2 <= %f", egmin, egmax) && vetoCut);
-
- // integrate coincidence beta spectrum
- n = hbg->Integral(binx1, binx2);
- DWU coinc(n, sqrt(n));
-
- // simulated beta spectrum
- TString fname = simdir + Form("b%i.root", b);
- TFile * fSim = new TFile(fname, "READ");
- TTree * tSim = (TTree*)fSim->Get("Detector");
-
- // histogram
- TString nbSim = Form("hbSim%i", i);
- TH1D * hbSim = new TH1D(nbSim, nbSim, bins, emin, emax);
-
- // apply veto?
- TCut cutSim("");
- if (veto) cutSim = TCut("EDepVeto < 10");
-
- // draw
- tSim->Draw("EDepSignal >> " + nbSim, cutSim);
-
- // integrate
- n = hbSim->Integral(binx1, binx2);
- DWU betasSim(n, sqrt(n));
-
- // normalisation
- double norm = ((TParameter<double>*)(tSim->GetUserInfo()->FindObject("NORMALIZATION")))->GetVal();
-
- // transmission
- betasSim.scale(1./norm);
- transmissionSim.push_back(betasSim);
-
- // scale spectrum and integral
- double x = decays.val / norm;
- hbSim->Scale(x);
- betasSim *= gammas;
- betasSim /= gammaEff0;
-
- // print data
- cout << Form(" Bmax: %.2f T", 0.01*b) << endl;
- cout << Form(" Duration: %.2f hours", hours) << endl;
- cout << Form(" Backgr: %.2E (%.1f%%)", backgr.val, backgr.unc/backgr.val*100.) << endl;
- cout << Form(" Betas: %.2E (%.1f%%)", betas.val, betas.unc/betas.val*100.) << endl;
- cout << Form(" Gammas: %.0f (%.1f%%) [centr: %.1f keV, sigma: %.2f keV]", gammas.val, gammas.unc/gammas.val*100., centroid, sigma) << endl;
- cout << Form(" Coinc: %.0f (%.1f%%)", coinc.val, coinc.unc/coinc.val*100.) << endl;
- cout << Form(" Sim betas: %.2E (%.1f%%)", betasSim.val, betasSim.unc/betasSim.val*100.) << endl;
-
- // gamma detection efficiency
- if (coinc.val >= 3) {
- gammaxvec.push_back(141.0*1E-2*b);
- gammaEff.push_back(coinc / betas);
- }
-
- // increment counter
- i++;
- }
-
-
- // --------- make graph of gamma efficiency ----------
- const int M = gammaEff.size();
- double x[M], y[M], dy[M];
- for (int i = 0; i < M; i++) {
- x[i] = gammaxvec[i];
- y[i] = gammaEff[i].val;
- dy[i] = gammaEff[i].unc;
- }
- TGraphErrors * gr = new TGraphErrors(M, x, y, 0, dy);
-
- // fit const
- TF1 * fconst = new TF1("fconst", "[0]", 0, 100);
- TFitResultPtr fitRes = gr->Fit("fconst", "ESBQ", "", 0, 100);
-
- // make pretty plot
- gr->SetMarkerStyle(21);
- gr->SetMarkerColor(4);
- gr->Draw("AP");
- gPad->SetTicks();
- gr->SetTitle("#gamma-ray detection efficiency");
- gr->GetXaxis()->SetTitle("I/I_{max} (%)");
- gr->GetYaxis()->SetTitle("#varepsilon_{#gamma}");
-
- // +- 1 sigma lines
- double x1 = gr->GetXaxis()->GetXmin();
- double x2 = gr->GetXaxis()->GetXmax();
- double ylow = fconst->GetParameter(0) - fconst->GetParError(0);
- TLine * low = new TLine(x1, ylow, x2, ylow);
- double yhigh = fconst->GetParameter(0) + fconst->GetParError(0);
- TLine * high = new TLine(x1, yhigh, x2, yhigh);
- low->SetLineColor(2);
- high->SetLineColor(2);
- low->SetLineStyle(2);
- high->SetLineStyle(2);
- low->Draw("same");
- high->Draw("same");
-
- // text box with average value
- TLatex *tx = new TLatex(32, ylow+2*(yhigh-ylow), Form("#varepsilon_{#gamma} = %.2f(%i)#times 10^{-4}", fconst->GetParameter(0)*1E4, (int)(fconst->GetParError(0)*1E6)));
- tx->SetTextAlign(11);
- tx->SetTextColor(1);
- tx->SetTextFont(43);
- tx->SetTextSize(24);
- tx->Draw("same");
-
- c1->Print("/home/oliskir/Dropbox/Work_ln-s/F20-IGISOL/analysis/figures/gammaEff.png");
-
- // print result
- cout << endl;
- cout << Form(" Gamma efficiency: %.2E (%.1f%%)", fconst->GetParameter(0), fconst->GetParError(0) / fconst->GetParameter(0)*100.) << endl;
- cout << Form(" Chi^2/dof: %.2f", fitRes->Chi2() / fitRes->Ndf()) << endl;
-
-
- // --------- make graph of transmission ----------
- const int Mt = transmission.size();
- double xt[Mt], yt[Mt], dyt[Mt], ytSim[Mt], ytSimLow[Mt], ytSimHigh[Mt];
- for (int i = 0; i < Mt; i++) {
- xt[i] = transxvec[i];
- yt[i] = 100*transmission[i].val;
- dyt[i] = 100*transmission[i].unc;
- ytSim[i] = 100*transmissionSim[i].val;
- double d = gammaEff0.unc / gammaEff0.val;
- ytSimLow[i] = (1.-d) * ytSim[i];
- ytSimHigh[i] = (1.+d) * ytSim[i];
- }
- TGraphErrors * grExp = new TGraphErrors(Mt, xt, yt, 0, dyt);
- TGraph * grSim = new TGraph(Mt, xt, ytSim);
- TGraph * grSimLow = new TGraph(Mt, xt, ytSimLow);
- TGraph * grSimHigh = new TGraph(Mt, xt, ytSimHigh);
-
- if (logy) gPad->SetLogy();
-
- grExp->SetMarkerStyle(21);
- grExp->SetMarkerColor(4);
- grExp->SetTitle("Exp.");
- grExp->SetFillStyle(0);
-
- grSim->SetLineColor(2);
- grSim->SetTitle("Sim.");
- grSim->SetFillStyle(0);
-
- grSimLow->SetLineStyle(2);
- grSimLow->SetLineColor(2);
- grSimLow->SetTitle("Sim. -1#sigma");
- grSimLow->SetFillStyle(0);
-
- grSimHigh->SetLineStyle(2);
- grSimHigh->SetLineColor(2);
- grSimHigh->SetTitle("Sim. +1#sigma");
- grSimHigh->SetFillStyle(0);
-
- TMultiGraph *mg = new TMultiGraph();
- mg->Add(grExp, "p");
- mg->Add(grSim, "c");
- mg->Add(grSimLow, "c");
- mg->Add(grSimHigh, "c");
- mg->Draw("a");
-
- gPad->SetTicks();
-
- TString mgTitle = Form("#beta transmission (%.1f-%.1f MeV)", emin/1E3, emax/1E3);
- if (veto) mgTitle += ", VETO";
- mg->SetTitle(mgTitle);
- mg->GetXaxis()->SetTitle("I/I_{max} (%)");
- mg->GetYaxis()->SetTitle("#varepsilon_{#beta} (%)");
- mg->Draw("a");
-
- auto legend = new TLegend(0.71, 0.70, 0.88, 0.88);
- legend->AddEntry(grExp, "Exp.", "pe");
- legend->AddEntry(grSim, "Sim.", "l");
- legend->AddEntry(grSimLow, "#pm1#sigma", "l");
- legend->SetTextSize(0.04);
- legend->Draw("same");
-
- TString figname = "transm";
- if (veto) figname += "_veto";
- if (logy) figname += "_log";
- figname += ".png";
-
- c1->Print("/home/oliskir/Dropbox/Work_ln-s/F20-IGISOL/analysis/figures/" + figname);
-}