diff --git a/VeikonKone.dox b/VeikonKone.dox
index 9feef7ef8c5006c9581a813f4808a78d988ac8e5..00957094131304ba0bd8087b1eb1e3ce3d4579a9 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/Desktop/f20/VeikonKone/include/
+INPUT = /home/oliskir/Desktop/f20/sim/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
diff --git a/compile.sh b/compile.sh
index 0226cae77c0fdbf30306298548a8418b99f39a13..a161863d50ab5a09fe0f210b2e565c6c2de53059 100755
--- a/compile.sh
+++ b/compile.sh
@@ -2,10 +2,10 @@ 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
+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
diff --git a/mac/f20scan.mac b/mac/f20scan.mac
index ab1783d8bb0fc06e5efb27eb27a899bce0246f2d..37ddcdc0dfdb5e268b5432d55e2ac5d8cd63a9fb 100644
--- a/mac/f20scan.mac
+++ b/mac/f20scan.mac
@@ -28,46 +28,46 @@
# run simulations
-/VK/field/maxbfield 0.10 tesla
-/VK/output/openFile output/f20scan/070617/b10.root
-/run/beamOn 10000
+###/VK/field/maxbfield 0.10 tesla
+###/VK/output/openFile output/f20scan/080617/b10.root
+###/run/beamOn 100000
-/VK/field/maxbfield 0.15 tesla
-/VK/output/openFile output/f20scan/070617/b15.root
-/run/beamOn 10000
+###/VK/field/maxbfield 0.15 tesla
+###/VK/output/openFile output/f20scan/080617/b15.root
+###/run/beamOn 100000
-/VK/field/maxbfield 0.20 tesla
-/VK/output/openFile output/f20scan/070617/b20.root
-/run/beamOn 10000
+###/VK/field/maxbfield 0.20 tesla
+###/VK/output/openFile output/f20scan/080617/b20.root
+###/run/beamOn 100000
-/VK/field/maxbfield 0.25 tesla
-/VK/output/openFile output/f20scan/070617/b25.root
-/run/beamOn 10000
+###/VK/field/maxbfield 0.25 tesla
+###/VK/output/openFile output/f20scan/080617/b25.root
+###/run/beamOn 100000
-/VK/field/maxbfield 0.30 tesla
-/VK/output/openFile output/f20scan/070617/b30.root
-/run/beamOn 10000
+###/VK/field/maxbfield 0.30 tesla
+###/VK/output/openFile output/f20scan/080617/b30.root
+###/run/beamOn 100000
-/VK/field/maxbfield 0.35 tesla
-/VK/output/openFile output/f20scan/070617/b35.root
-/run/beamOn 10000
+###/VK/field/maxbfield 0.35 tesla
+###/VK/output/openFile output/f20scan/080617/b35.root
+###/run/beamOn 100000
-/VK/field/maxbfield 0.40 tesla
-/VK/output/openFile output/f20scan/070617/b40.root
-/run/beamOn 10000
+###/VK/field/maxbfield 0.40 tesla
+###/VK/output/openFile output/f20scan/080617/b40.root
+###/run/beamOn 100000
-/VK/field/maxbfield 0.44 tesla
-/VK/output/openFile output/f20scan/070617/b44.root
-/run/beamOn 10000
+###/VK/field/maxbfield 0.44 tesla
+###/VK/output/openFile output/f20scan/080617/b44.root
+###/run/beamOn 100000
/VK/field/maxbfield 0.46 tesla
-/VK/output/openFile output/f20scan/070617/b46.root
-/run/beamOn 10000
+/VK/output/openFile output/f20scan/080617/b46.root
+/run/beamOn 1000000
/VK/field/maxbfield 0.48 tesla
-/VK/output/openFile output/f20scan/070617/b48.root
-/run/beamOn 10000
+/VK/output/openFile output/f20scan/080617/b48.root
+/run/beamOn 1000000
/VK/field/maxbfield 0.50 tesla
-/VK/output/openFile output/f20scan/070617/b50.root
-/run/beamOn 10000
+/VK/output/openFile output/f20scan/080617/b50.root
+/run/beamOn 1000000
diff --git a/output/f20scan/TransmissionExp.C b/output/f20scan/TransmissionExp.C
index 6ec42d11e45e98fe20fee528f05f6b2f512b316c..b5b66b1afc9af3a38fb87aa5e7722b48b61ca51f 100644
--- a/output/f20scan/TransmissionExp.C
+++ b/output/f20scan/TransmissionExp.C
@@ -1,4 +1,114 @@
+
+
+/*
+ * 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.
*/
@@ -41,10 +151,6 @@ double GetDuration(TChain * chain)
void FitGammaPeak(TH1D * h0, double& counts, double& e0, double& sigma, double& dcounts, double& de0, double& dsigma)
{
- // TSpectrum
- int maxPeaks = 1;
- TSpectrum * tspec = new TSpectrum(maxPeaks);
-
// make a copy
TH1D * h = (TH1D*)h0->Clone("hcopy");
h->GetXaxis()->SetRangeUser(1350, 1700);
@@ -52,38 +158,31 @@ void FitGammaPeak(TH1D * h0, double& counts, double& e0, double& sigma, double&
// get bin width
double binWidth = h->GetBinWidth(1);
- // estimate background
- TH1D * hbackgr = (TH1D*)tspec->Background(h);
-
// search for 1.6-MeV peak
- tspec->Search(h);
- double * xPos = tspec->GetPositionX();
- double * yPos = tspec->GetPositionX();
-
- // subtract background
- TH1D * hsub = (TH1D*)h->Clone("hsub");
- TH1D * hinv = (TH1D*)hbackgr->Clone("hinv");
- hinv->Scale(-1);
-//*** hsub->Add(hinv);
+ 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[0] - 50;
- double fitmax = xPos[0] + 50;
+ double fitmin = xPos - 75;
+ double fitmax = xPos + 75;
// start values
- fitf->SetParameter(0, yPos[0]);
- fitf->SetParameter(1, xPos[0]);
+ fitf->SetParameter(0, yPos);
+ fitf->SetParameter(1, xPos);
fitf->SetParameter(2, 20);
fitf->SetParameter(3, 0);
// fit
- hsub->Fit("fitf", "EBQ", "", fitmin, fitmax);
+ h->Fit("fitf", "EBQ", "", fitmin, fitmax);
fitmin = fitf->GetParameter(1) - 3*fitf->GetParameter(2);
fitmax = fitf->GetParameter(1) + 3*fitf->GetParameter(2);
- hsub->Fit("fitf", "ES+BQ", "", fitmin, fitmax);
+ h->Fit("fitf", "ES+BQ", "", fitmin, fitmax);
// get fitted values
counts = fitf->GetParameter(0) / binWidth;
@@ -95,26 +194,129 @@ void FitGammaPeak(TH1D * h0, double& counts, double& e0, double& sigma, double&
}
+/*
+ * 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, double emin = 100, double emax = 5400)
+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 = "070617/";
+ TString simdir = "080617/";
// gamma detection efficiency
- const double gammaEff = 0.54E-4;
- const double dgammaEff = 0.10E-4;
+ 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};//{15, 20, 25, 30, 35, 40, 44, 46, 48, 50};
+ vector<int> bmax = {15, 20, 25, 30, 35, 40, 44, 46, 48, 50};
// runs at each setting
vector<int> runs[N];
@@ -134,7 +336,7 @@ void TransmissionExp(bool veto = false, double emin = 100, double emax = 5400)
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");
+ cBackgr->Draw("E0 >> hBackgr0", "Z0 != 10" && vetoCut);
double hoursBackgr = GetDuration(cBackgr);
cout << endl;
@@ -146,7 +348,18 @@ void TransmissionExp(bool veto = false, double emin = 100, double emax = 5400)
// 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)
@@ -163,50 +376,63 @@ void TransmissionExp(bool veto = false, double emin = 100, double emax = 5400)
// integrate background spectrum
TH1D * hBackgr = (TH1D*)hBackgr0->Clone("hBackgr");
- int binx1 = hBackgr->GetXaxis()->FindBin(emin);
- int binx2 = hBackgr->GetXaxis()->FindBin(emax);
- double backgr = hBackgr->Integral(binx1, binx2);
- double dbackgr = sqrt(backgr);
+ n = hBackgr->Integral(binx1, binx2);
+ DWU backgr(n, sqrt(n));
// scale
- backgr *= hours/hoursBackgr;
- dbackgr *= hours/hoursBackgr;
+ 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, 1000, 0, 3000);
+ TH1D * hg = new TH1D(ng, ng, 600, 0, 3000);
// draw
- cExp->Draw("E0 >> " + nb, "M0==1 && Z0!=10");
- cExp->Draw("E2 >> " + ng, "M2==1 && Z2!=10");
+ cExp->Draw("E0 >> " + nb, "M0==1 && Z0!=10" && vetoCut);
+ cExp->Draw("E2 >> " + ng, "M2==1 && Z2!=10" && vetoCut);
// integrate beta spectrum
- double betas = hb->Integral(binx1, binx2);
- double dbetas = sqrt(betas);
+ n = hb->Integral(binx1, binx2);
+ DWU betas(n, sqrt(n));
// subtract background
betas -= backgr;
- dbetas = sqrt(pow(dbetas, 2) + pow(dbackgr, 2));
// integrate 1.6-MeV peak in singles spectrum
- double gammas, dgammas, centroid, dcentroid, sigma, dsigma;
- FitGammaPeak(hg, gammas, centroid, sigma, dgammas, dcentroid, dsigma);
+ 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, 1000, 0, 3000);
+ TH1F * hbg = new TH1F(nbg, nbg, bins, emin, emax);
// draw
- double egmin = centroid - 3*sigma;
- double egmax = centroid + 3*sigma;
- cExp->Draw("E0 >> " + nbg, Form("M0==1 && M2==1 && Z0 != 10 && Z2 != 10 && E2 >= %f && E2 <= %f", egmin, egmax));
+ cExp->Draw("E0 >> " + nbg, Form("M0==1 && M2==1 && Z0 != 10 && Z2 != 10 && E2 >= %f && E2 <= %f", egmin, egmax) && vetoCut);
// integrate coincidence beta spectrum
- double coinc = hbg->Integral(binx1, binx2);
- double dcoinc = sqrt(coinc);
+ n = hbg->Integral(binx1, binx2);
+ DWU coinc(n, sqrt(n));
// simulated beta spectrum
TString fname = simdir + Form("b%i.root", b);
@@ -217,30 +443,168 @@ void TransmissionExp(bool veto = false, double emin = 100, double emax = 5400)
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, "");
+ tSim->Draw("EDepSignal >> " + nbSim, cutSim);
// integrate
- double betasSim = hbSim->Integral(binx1, binx2);
- double dbetasSim = sqrt(betasSim);
+ n = hbSim->Integral(binx1, binx2);
+ DWU betasSim(n, sqrt(n));
- // scaling factor
- double decays = gammas / gammaEff;
- double norm = ((TParameter<double>*)(tSim->GetUserInfo()->FindObject("NORMALIZATION")))->GetVal();
- double x = decays / norm;
+ // 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);
- double relUncert = sqrt(pow(dbetasSim/betasSim, 2) + pow(dgammas/gammas, 2) + pow(dgammaEff/gammaEff, 2));
- betasSim *= x;
- dbetasSim = betasSim * relUncert;
+ 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, dbackgr/backgr*100.) << endl;
- cout << Form(" Betas: %.2E (%.1f%%)", betas, dbetas/betas*100.) << endl;
- cout << Form(" Gammas: %.0f (%.1f%%)", gammas, dgammas/gammas*100.) << endl;
- cout << Form(" Coinc: %.0f (%.1f%%)", coinc, dcoinc/coinc*100.) << endl;
- cout << Form(" Sim betas: %.2E (%.1f%%)", betasSim, dbetasSim/betasSim*100.) << 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);
}