diff --git a/CMakeLists.txt b/CMakeLists.txt
index 220c0494b518f349abba4556453aa5297510eac7..7de020ac0443f61604ea0b6e6cc934299d48ec4f 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -134,17 +134,17 @@ endif()
# relies on these scripts being in the current working directory.
#
set(VeikonKone_SCRIPTS
- LXe.out
- LXe.in
- defaults.mac
- cerenkov.mac
- wls.mac
- photon.mac
- reviewEvent.mac
- gui.mac
- icons.mac
- run.png
- vis.mac
+ mac/LXe.out
+ mac/LXe.in
+ mac/defaults.mac
+ mac/cerenkov.mac
+ mac/wls.mac
+ mac/photon.mac
+ mac/reviewEvent.mac
+ mac/gui.mac
+ mac/icons.mac
+ mac/run.png
+ mac/vis.mac
)
foreach(_script ${VeikonKone_SCRIPTS})
diff --git a/VeikonKone.dox b/VeikonKone.dox
index 00957094131304ba0bd8087b1eb1e3ce3d4579a9..9feef7ef8c5006c9581a813f4808a78d988ac8e5 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/sim/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
diff --git a/compile.sh b/compile.sh
index a161863d50ab5a09fe0f210b2e565c6c2de53059..0226cae77c0fdbf30306298548a8418b99f39a13 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/LXe.cpp b/mac/LXe.cpp
similarity index 100%
rename from LXe.cpp
rename to mac/LXe.cpp
diff --git a/LXe.in b/mac/LXe.in
similarity index 100%
rename from LXe.in
rename to mac/LXe.in
diff --git a/LXe.out b/mac/LXe.out
similarity index 100%
rename from LXe.out
rename to mac/LXe.out
diff --git a/run.png b/mac/run.png
similarity index 100%
rename from run.png
rename to mac/run.png
diff --git a/output/f20scan/TransmissionExp.C b/output/f20scan/TransmissionExp.C
index a785330833456abc37410af0ab81ce50d539e83a..6ec42d11e45e98fe20fee528f05f6b2f512b316c 100644
--- a/output/f20scan/TransmissionExp.C
+++ b/output/f20scan/TransmissionExp.C
@@ -6,14 +6,8 @@
double GetDuration(TChain * chain)
{
// branches
- UInt_t M0, M1, M2;
- chain->SetBranchAddress("M0", &M0);
- chain->SetBranchAddress("M1", &M1);
- chain->SetBranchAddress("M2", &M2);
- vector<ULong64_t> T0(10), T1(10), T2(10);
- chain->SetBranchAddress("T0", T0.data());
- chain->SetBranchAddress("T1", T1.data());
- chain->SetBranchAddress("T2", T2.data());
+ ULong64_t tms;
+ chain->SetBranchAddress("tms", &tms);
// find tmin and tmax
int shifts = 0;
@@ -22,14 +16,7 @@ double GetDuration(TChain * chain)
for (int i=0; i<chain->GetEntries(); i++)
{
chain->GetEntry(i);
-
- double T00 = 0, T10 = 0, T20 = 0;
- if (M0 > 0) T00 = T0[0];
- if (M1 > 0) T10 = T1[0];
- if (M2 > 0) T20 = T2[0];
-
- double T = TMath::Max(T00, TMath::Max(T10, T20));
-
+ double T = tms;
if (Tp > T) {
shift += Tp;
shifts++;
@@ -39,9 +26,11 @@ double GetDuration(TChain * chain)
Tp = T;
}
- double seconds = (TLast - TFirst) * 1E-8;
+ double seconds = (TLast - TFirst) * 1E-3;
double minutes = seconds / 60.;
double hours = minutes / 60.;
+
+ chain->ResetBranchAddresses();
return hours;
}
@@ -50,8 +39,59 @@ double GetDuration(TChain * chain)
* Locate and fit 1.6-MeV gamma peak
*/
-void FitGammaPeak(TH1D * h, double& counts, double& e0, double& sigma)
+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);
+
+ // 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);
+
+ // fit function
+ TF1 * fitf = new TF1("fitf", "[0]*TMath::Gaus(x,[1],[2],true)+[3]", 0, 1E4);
+
+ // fit range
+ double fitmin = xPos[0] - 50;
+ double fitmax = xPos[0] + 50;
+
+ // start values
+ fitf->SetParameter(0, yPos[0]);
+ fitf->SetParameter(1, xPos[0]);
+ fitf->SetParameter(2, 20);
+ fitf->SetParameter(3, 0);
+
+ // fit
+ hsub->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);
+
+ // 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);
}
@@ -59,23 +99,22 @@ void FitGammaPeak(TH1D * h, double& counts, double& e0, double& sigma)
* Plot transmission as a function of magnetic-field setting.
*/
-void TransmissionExp(bool veto = false)
+void TransmissionExp(bool veto = false, double emin = 100, double emax = 5400)
{
// directories
TString expdir = "/home/oliskir/Desktop/f20/data/sorted/new/";
TString simdir = "070617/";
// gamma detection efficiency
- const double gammaEff = 0.54E-6;
+ const double gammaEff = 0.54E-4;
+ const double dgammaEff = 0.10E-4;
- // beta spectrum integration limits
- const double emin = 100;
- const double emax = 5400;
- const int bins = 530;
+ // 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};
+ vector<int> bmax = {15};//{15, 20, 25, 30, 35, 40, 44, 46, 48, 50};
// runs at each setting
vector<int> runs[N];
@@ -87,7 +126,7 @@ void TransmissionExp(bool veto = false)
runs[5] = {334};
runs[6] = {369};
runs[7] = {360};
- runs[8] = {361, 362, 363, 364, 365, 370, 371, 372, 373, 374, 375, 376, 379, 380};
+ runs[8] = {361, 362, 363, 364, 365, 370, 371, 372, 373, 375, 376, 379, 380};
runs[9] = {358};
// background run (337)
@@ -112,8 +151,6 @@ void TransmissionExp(bool veto = false)
int i = 0;
for (auto& b : bmax)
{
- cout << " Bmax: " << 0.01*b << " T" << endl;
-
// add files
TChain * cExp = new TChain("data");
for (auto& r : runs[i]) {
@@ -123,17 +160,19 @@ void TransmissionExp(bool veto = false)
// duration of experimental run
double hours = GetDuration(cExp);
- cout << Form(" Duration: %.1f", hours) << endl;
- // scale background spectrum
+ // 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);
+
+ // scale
+ backgr *= hours/hoursBackgr;
+ dbackgr *= hours/hoursBackgr;
hBackgr->Scale(hours/hoursBackgr);
- // integrate background spectrum
- int binx1 = hBackgr->GetXaxis()->FindBin(emin);
- int binx2 = hBackgr->GetXaxis()->FindBin(emax);
- double backgr = hBackgr->Integral(binx1, binx2);
-
// singles histograms
TString nb = Form("hb%i", i);
TH1D * hb = new TH1D(nb, nb, bins, emin, emax);
@@ -141,16 +180,20 @@ void TransmissionExp(bool veto = false)
TH1D * hg = new TH1D(ng, ng, 1000, 0, 3000);
// draw
- cExp->Draw("E0 >> " + nb, "Z0 != 10");
- cExp->Draw("E2 >> " + ng, "Z2 != 10");
-
+ cExp->Draw("E0 >> " + nb, "M0==1 && Z0!=10");
+ cExp->Draw("E2 >> " + ng, "M2==1 && Z2!=10");
+
// integrate beta spectrum
- double betas = hb->Integral(binx1, binx2);
+ double betas = hb->Integral(binx1, binx2);
+ double dbetas = sqrt(betas);
+
+ // subtract background
betas -= backgr;
+ dbetas = sqrt(pow(dbetas, 2) + pow(dbackgr, 2));
// integrate 1.6-MeV peak in singles spectrum
- double gammas, centroid, sigma;
- FitGammaPeak(hg, gammas, centroid, sigma);
+ double gammas, dgammas, centroid, dcentroid, sigma, dsigma;
+ FitGammaPeak(hg, gammas, centroid, sigma, dgammas, dcentroid, dsigma);
// coincidence histogram
TString nbg = Form("hbg%i", i);
@@ -159,10 +202,11 @@ void TransmissionExp(bool veto = false)
// draw
double egmin = centroid - 3*sigma;
double egmax = centroid + 3*sigma;
- cExp->Draw("E0 >> " + nbg, Form("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));
// integrate coincidence beta spectrum
- double coincidences = hbg->Integral(binx1, binx2);
+ double coinc = hbg->Integral(binx1, binx2);
+ double dcoinc = sqrt(coinc);
// simulated beta spectrum
TString fname = simdir + Form("b%i.root", b);
@@ -176,21 +220,27 @@ void TransmissionExp(bool veto = false)
// draw
tSim->Draw("EDepSignal >> " + nbSim, "");
+ // integrate
+ double betasSim = hbSim->Integral(binx1, binx2);
+ double dbetasSim = sqrt(betasSim);
+
// scaling factor
double decays = gammas / gammaEff;
double norm = ((TParameter<double>*)(tSim->GetUserInfo()->FindObject("NORMALIZATION")))->GetVal();
- hbSim->Scale(decays/norm);
+ double x = decays / norm;
- // integrate
- double betasSim = hbSim->Integral(binx1, binx2);
-
- cout << Form(" Background: %.0f", backgr) << endl;
- cout << Form(" Betas: %.0f", betas) << endl;
- cout << Form(" Gammas: %.0f", gammas) << endl;
- cout << Form(" Coincidences: %.0f", coincidences) << endl;
- cout << Form(" Simulated betas: %.0f", betasSim) << endl;
-
- cin.get();
- i++;
+ // scale spectrum and integral
+ hbSim->Scale(x);
+ double relUncert = sqrt(pow(dbetasSim/betasSim, 2) + pow(dgammas/gammas, 2) + pow(dgammaEff/gammaEff, 2));
+ betasSim *= x;
+ dbetasSim = betasSim * relUncert;
+
+ 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;
}
}