Skip to content
Snippets Groups Projects
Commit 9454c78b authored by Oliver Kirsebom's avatar Oliver Kirsebom
Browse files

bla

parent 514cb33e
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
CMakeLists.txt
+ 11
− 11
View file @ 9454c78b
......@@ -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})
......
VeikonKone.dox
+ 1
− 1
View file @ 9454c78b
......@@ -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
......
compile.sh
+ 2
− 2
View file @ 9454c78b
......@@ -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
......
LXe.cpp → mac/LXe.cpp
+ 0
− 0
View file @ 9454c78b
File moved
LXe.in → mac/LXe.in
+ 0
− 0
View file @ 9454c78b
File moved
LXe.out → mac/LXe.out
+ 0
− 0
View file @ 9454c78b
File moved
run.png → mac/run.png
+ 0
− 0
View file @ 9454c78b
File moved
output/f20scan/TransmissionExp.C
+ 105
− 55
View file @ 9454c78b
......@@ -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;
}
}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment