LMU/musrfit
LMU/musrfit
5
0
Fork 0
Code Issues 6 Pull Requests Actions Packages Projects Releases 3 Activity

corrected data handling (see MUSR-58)

Browse Source
This commit is contained in:
nemu 2009-05-19 14:47:24 +00:00
parent 2a12a17881
commit 0b2e989141
8 changed files with 117 additions and 66 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
src/classes/PMusrCanvas.cpp
View File

@ -1382,7 +1382,7 @@ void PMusrCanvas::HandleDataSet(unsigned int plotNo, unsigned int runNo, PRunDat
end = data->fDataTimeStart + data->fValue.size()*data->fDataTimeStep + data->fDataTimeStep/2.0;
// invoke histo
dataHisto = new TH1F(name, name, data->fValue.size()+2, start, end);
dataHisto = new TH1F(name, name, data->fValue.size()+1, start, end);
// fill histogram
for (unsigned int i=0; i<data->fValue.size(); i++) {
@ -1422,7 +1422,7 @@ void PMusrCanvas::HandleDataSet(unsigned int plotNo, unsigned int runNo, PRunDat
//cout << endl << ">> start = " << start << ", end = " << end << endl;
// invoke histo
theoHisto = new TH1F(name, name, data->fTheory.size()+2, start, end);
theoHisto = new TH1F(name, name, data->fTheory.size()+1, start, end);
// fill histogram
for (unsigned int i=0; i<data->fTheory.size(); i++) {

93
src/classes/PRunAsymmetry.cpp
View File

@ -699,7 +699,7 @@ bool PRunAsymmetry::PrepareFitData(PRawRunData* runData, unsigned int histoNo[2]
double f, b, ef, eb;
// fill data time start, and step
// data start at data_start-t0
fData.fDataTimeStart = fTimeResolution*(((double)start[0]-t0[0])+(double)fRunInfo->fPacking/2.0);
fData.fDataTimeStart = fTimeResolution*((double)start[0]-t0[0]);
fData.fDataTimeStep = fTimeResolution*(double)fRunInfo->fPacking;
for (unsigned int i=0; i<noOfBins; i++) {
// to make the formulae more readable
@ -761,11 +761,22 @@ bool PRunAsymmetry::PrepareViewData(PRawRunData* runData, unsigned int histoNo[2
// transform raw histo data. This is done the following way (for details see the manual):
// first rebin the data, than calculate the asymmetry
// first get start data, end data, and t0
int start[2] = {fRunInfo->fDataRange[0]-fRunInfo->fPacking*(fRunInfo->fDataRange[0]/fRunInfo->fPacking),
fRunInfo->fDataRange[2]-fRunInfo->fPacking*(fRunInfo->fDataRange[2]/fRunInfo->fPacking)};
int val = fRunInfo->fDataRange[0]-fRunInfo->fPacking*(fRunInfo->fDataRange[0]/fRunInfo->fPacking);
do {
if (fRunInfo->fDataRange[2] - fRunInfo->fDataRange[0] < 0)
val += fRunInfo->fPacking;
} while (val + fRunInfo->fDataRange[2] - fRunInfo->fDataRange[0] < 0);
int start[2] = {val, val + fRunInfo->fDataRange[2] - fRunInfo->fDataRange[0]};
int end[2];
double t0[2] = {fT0s[0], fT0s[1]};
/*
cout << endl << ">> start[0]=" << start[0] << ", end[0]=" << end[0];
cout << endl << ">> start[1]=" << start[1] << ", end[1]=" << end[1];
cout << endl;
*/
// make sure that there are equal number of rebinned bins in forward and backward
unsigned int noOfBins0 = (runData->fDataBin[histoNo[0]].size()-start[0])/fRunInfo->fPacking;
unsigned int noOfBins1 = (runData->fDataBin[histoNo[1]].size()-start[1])/fRunInfo->fPacking;
@ -806,34 +817,44 @@ bool PRunAsymmetry::PrepareViewData(PRawRunData* runData, unsigned int histoNo[2
double error = 0.0;
// forward
for (int i=start[0]; i<end[0]; i++) {
if (((i-start[0]) % fRunInfo->fPacking == 0) && (i != start[0])) { // fill data
// in order that after rebinning the fit does not need to be redone (important for plots)
// the value is normalize to per bin
value /= fRunInfo->fPacking;
forwardPacked.fValue.push_back(value);
if (value == 0.0)
forwardPacked.fError.push_back(1.0);
else
forwardPacked.fError.push_back(TMath::Sqrt(error)/fRunInfo->fPacking);
value = 0.0;
error = 0.0;
if (fRunInfo->fPacking == 1) {
forwardPacked.fValue.push_back(fForward[i]);
forwardPacked.fError.push_back(fForwardErr[i]);
} else { // packed data, i.e. fRunInfo->fPacking > 1
if (((i-start[0]) % fRunInfo->fPacking == 0) && (i != start[0])) { // fill data
// in order that after rebinning the fit does not need to be redone (important for plots)
// the value is normalize to per bin
value /= fRunInfo->fPacking;
forwardPacked.fValue.push_back(value);
if (value == 0.0)
forwardPacked.fError.push_back(1.0);
else
forwardPacked.fError.push_back(TMath::Sqrt(error)/fRunInfo->fPacking);
value = 0.0;
error = 0.0;
}
}
value += fForward[i];
error += fForwardErr[i]*fForwardErr[i];
}
// backward
for (int i=start[1]; i<end[1]; i++) {
if (((i-start[1]) % fRunInfo->fPacking == 0) && (i != start[1])) { // fill data
// in order that after rebinning the fit does not need to be redone (important for plots)
// the value is normalize to per bin
value /= fRunInfo->fPacking;
backwardPacked.fValue.push_back(value);
if (value == 0.0)
backwardPacked.fError.push_back(1.0);
else
backwardPacked.fError.push_back(TMath::Sqrt(error)/fRunInfo->fPacking);
value = 0.0;
error = 0.0;
if (fRunInfo->fPacking == 1) {
backwardPacked.fValue.push_back(fBackward[i]);
backwardPacked.fError.push_back(fBackwardErr[i]);
} else { // packed data, i.e. fRunInfo->fPacking > 1
if (((i-start[1]) % fRunInfo->fPacking == 0) && (i != start[1])) { // fill data
// in order that after rebinning the fit does not need to be redone (important for plots)
// the value is normalize to per bin
value /= fRunInfo->fPacking;
backwardPacked.fValue.push_back(value);
if (value == 0.0)
backwardPacked.fError.push_back(1.0);
else
backwardPacked.fError.push_back(TMath::Sqrt(error)/fRunInfo->fPacking);
value = 0.0;
error = 0.0;
}
}
value += fBackward[i];
error += fBackwardErr[i]*fBackwardErr[i];
@ -851,9 +872,14 @@ bool PRunAsymmetry::PrepareViewData(PRawRunData* runData, unsigned int histoNo[2
double f, b, ef, eb, alpha = 1.0, beta = 1.0;
// fill data time start, and step
// data start at data_start-t0
fData.fDataTimeStart = fTimeResolution*(((double)start[0]-t0[0])+(double)fRunInfo->fPacking/2.0);
fData.fDataTimeStart = fTimeResolution*((double)start[0]-t0[0]);
fData.fDataTimeStep = fTimeResolution*(double)fRunInfo->fPacking;
/*
cout << endl << ">> start time = " << fData.fDataTimeStart << ", step = " << fData.fDataTimeStep;
cout << endl << "--------------------------------" << endl;
*/
// get the proper alpha and beta
switch (fAlphaBetaTag) {
case 1: // alpha == 1, beta == 1
@ -877,7 +903,7 @@ bool PRunAsymmetry::PrepareViewData(PRawRunData* runData, unsigned int histoNo[2
}
//cout << endl << ">> alpha = " << alpha << ", beta = " << beta;
for (unsigned int i=0; i<noOfBins; i++) {
for (unsigned int i=0; i<forwardPacked.fValue.size(); i++) {
// to make the formulae more readable
f = forwardPacked.fValue[i];
b = backwardPacked.fValue[i];
@ -924,11 +950,16 @@ bool PRunAsymmetry::PrepareViewData(PRawRunData* runData, unsigned int histoNo[2
// calculate theory
unsigned int size = runData->fDataBin[histoNo[0]].size();
double startTime = -fT0s[0]*fTimeResolution;
fData.fTheoryTimeStart = startTime;
fData.fTheoryTimeStep = fTimeResolution;
double factor = 1.0;
if (fData.fValue.size() * 10 > runData->fDataBin[histoNo[0]].size()) {
size = fData.fValue.size() * 10;
factor = (double)runData->fDataBin[histoNo[0]].size() / (double)size;
}
//cout << endl << ">> runData->fDataBin[histoNo[0]].size() = " << runData->fDataBin[histoNo[0]].size() << ", fData.fValue.size() * 10 = " << fData.fValue.size() * 10 << ", size = " << size << ", factor = " << factor << endl;
fData.fTheoryTimeStart = fData.fDataTimeStart;
fData.fTheoryTimeStep = fTimeResolution*factor;
for (unsigned int i=0; i<size; i++) {
time = startTime + (double)i*fTimeResolution;
time = fData.fTheoryTimeStart + (double)i*fTimeResolution*factor;
value = fTheory->Func(time, par, fFuncValues);
if (fabs(value) > 10.0) { // dirty hack needs to be fixed!!
value = 0.0;

1
src/classes/PRunDataHandler.cpp
View File

@ -865,6 +865,7 @@ cout << endl;
for (int j=0; j<psiBin.get_histoLength_bin(); j++) {
histoData.push_back(histo[j]);
}
delete histo;
runData.fDataBin.push_back(histoData);
histoData.clear();
}

36
src/classes/PRunSingleHisto.cpp
View File

@ -501,11 +501,13 @@ bool PRunSingleHisto::PrepareFitData(PRawRunData* runData, const unsigned int hi
double normalizer = 1.0;
// data start at data_start-t0
// time shifted so that packing is included correctly, i.e. t0 == t0 after packing
fData.fDataTimeStart = fTimeResolution*(((double)start-(double)t0)+(double)(fRunInfo->fPacking-1)/2.0);
fData.fDataTimeStart = fTimeResolution*((double)start-(double)t0);
fData.fDataTimeStep = fTimeResolution*fRunInfo->fPacking;
for (int i=start; i<end; i++) {
if (fRunInfo->fPacking == 1) {
value = runData->fDataBin[histoNo][i];
normalizer = fRunInfo->fPacking * (fTimeResolution * 1e3); // fTimeResolution us->ns
value /= normalizer;
fData.fValue.push_back(value);
if (value == 0.0)
fData.fError.push_back(1.0);
@ -579,8 +581,9 @@ bool PRunSingleHisto::PrepareRawViewData(PRawRunData* runData, const unsigned in
double value = 0.0;
// data start at data_start-t0
// time shifted so that packing is included correctly, i.e. t0 == t0 after packing
fData.fDataTimeStart = fTimeResolution*((double)start-(double)t0+(double)(fRunInfo->fPacking-1)/2.0);
fData.fDataTimeStart = fTimeResolution*((double)start-(double)t0);
fData.fDataTimeStep = fTimeResolution*fRunInfo->fPacking;
/*
cout << endl << ">> time resolution = " << fTimeResolution;
cout << endl << ">> start = " << start << ", t0 = " << t0 << ", packing = " << fRunInfo->fPacking;
@ -662,11 +665,17 @@ cout << endl << ">> data start time = " << fData.fDataTimeStart;
// calculate theory
unsigned int size = runData->fDataBin[histoNo].size();
double factor = 1.0;
if (fData.fValue.size() * 10 > runData->fDataBin[histoNo].size()) {
size = fData.fValue.size() * 10;
factor = (double)runData->fDataBin[histoNo].size() / (double)size;
}
//cout << endl << ">> runData->fDataBin[histoNo].size() = " << runData->fDataBin[histoNo].size() << ", fData.fValue.size() * 10 = " << fData.fValue.size() * 10 << ", size = " << size << ", factor = " << factor << endl;
double theoryValue;
fData.fTheoryTimeStart = fData.fDataTimeStart;
fData.fTheoryTimeStep = fTimeResolution;
fData.fTheoryTimeStep = fTimeResolution*factor;
for (unsigned int i=0; i<size; i++) {
time = fData.fTheoryTimeStart + i*fTimeResolution;
time = fData.fTheoryTimeStart + i*fTimeResolution*factor;
theoryValue = fTheory->Func(time, par, fFuncValues);
if (fabs(theoryValue) > 10.0) { // dirty hack needs to be fixed!!
theoryValue = 0.0;
@ -773,14 +782,15 @@ bool PRunSingleHisto::PrepareViewData(PRawRunData* runData, const unsigned int h
double time;
// data start at data_start-t0
// time shifted so that packing is included correctly, i.e. t0 == t0 after packing
fData.fDataTimeStart = fTimeResolution*((double)start-(double)t0+(double)(fRunInfo->fPacking-1)/2.0);
fData.fDataTimeStart = fTimeResolution*((double)start-(double)t0);
fData.fDataTimeStep = fTimeResolution*fRunInfo->fPacking;
/*
cout << endl << ">> start time = " << fData.fDataTimeStart << ", step = " << fData.fDataTimeStep;
cout << endl << ">> start = " << start << ", end = " << end;
cout << endl << "--------------------------------";
cout << endl << "--------------------------------" << endl;
*/
double normalizer = 1.0;
for (int i=start; i<end; i++) {
if (((i-start) % fRunInfo->fPacking == 0) && (i != start)) { // fill data
@ -815,13 +825,19 @@ cout << endl << "--------------------------------";
}
// calculate theory
unsigned int size = runData->fDataBin[histoNo].size();
double theoryValue;
unsigned int size = runData->fDataBin[histoNo].size();
double factor = 1.0;
if (fData.fValue.size() * 10 > runData->fDataBin[histoNo].size()) {
size = fData.fValue.size() * 10;
factor = (double)runData->fDataBin[histoNo].size() / (double)size;
}
//cout << endl << ">> runData->fDataBin[histoNo].size() = " << runData->fDataBin[histoNo].size() << ", fData.fValue.size() * 10 = " << fData.fValue.size() * 10 << ", size = " << size << ", factor = " << factor << endl;
fData.fTheoryTimeStart = fData.fDataTimeStart;
fData.fTheoryTimeStep = fTimeResolution;
fData.fTheoryTimeStep = fTimeResolution*factor;
//cout << endl << ">> size=" << size << ", startTime=" << startTime << ", fTimeResolution=" << fTimeResolution;
for (unsigned int i=0; i<size; i++) {
time = fData.fTheoryTimeStart + (double)i*fTimeResolution;
time = fData.fTheoryTimeStart + (double)i*fTimeResolution*factor;
theoryValue = fTheory->Func(time, par, fFuncValues);
if (fabs(theoryValue) > 10.0) { // dirty hack needs to be fixed!!
theoryValue = 0.0;

30
src/tests/skewedGaussianTest/fakeData.cpp
View File

@ -83,7 +83,8 @@ void fakeDataSyntax()
int main(int argc, char *argv[])
{
const Double_t gamma_mu = TMath::TwoPi() * 0.1355; // in (rad/ns/T)
const Double_t tau_mu = 2197.147; // muon life time in ns
// muon life time in (us), see PRL99, 032001 (2007)
const Double_t tau_mu = 2197.019; // muon life time in ns
if ((argc != 4) && (argc !=6)) {
fakeDataSyntax();
@ -285,10 +286,11 @@ int main(int argc, char *argv[])
for (UInt_t i=0; i<asym.size(); i++) { // loop over all histos
for (Int_t j=0; j<noOfChannels; j++) { // loop over time
for (UInt_t k=0; k<pB.size(); k++) { // loop over p(B)
if (j < t0[i])
if (j < t0[i]) {
dval += 0.0;
else
} else {
dval += pB[k]*TMath::Cos(gamma_mu*B[k]*timeResolution*(j-t0[i])+phase[i]);
}
}
ddata.push_back(asym[i]*dval);
dval = 0.0;
@ -353,10 +355,10 @@ int main(int argc, char *argv[])
name += i;
title = "asym";
title += i;
hh = new TH1F(name.Data(), title.Data(), noOfChannels,
-timeResolution/2.0, (noOfChannels+0.5)*timeResolution);
hh = new TH1F(name.Data(), title.Data(), noOfChannels+1,
-timeResolution*0.5, (noOfChannels+0.5)*timeResolution);
for (Int_t j=0; j<noOfChannels; j++) {
hh->SetBinContent(j, asymmetry[i][j]);
hh->SetBinContent(j+1, asymmetry[i][j]);
}
hh->Draw("*H HIST");
@ -380,10 +382,10 @@ int main(int argc, char *argv[])
name += i;
title = "histo";
title += i;
hh = new TH1F(name.Data(), title.Data(), noOfChannels,
-timeResolution/2.0, (noOfChannels+0.5)*timeResolution);
hh = new TH1F(name.Data(), title.Data(), noOfChannels+1,
-timeResolution*0.5, (noOfChannels+0.5)*timeResolution);
for (Int_t j=0; j<noOfChannels; j++) {
hh->SetBinContent(j, histo[i][j]);
hh->SetBinContent(j+1, histo[i][j]);
}
hh->Draw("*H HIST");
@ -411,18 +413,18 @@ int main(int argc, char *argv[])
// create histos
name = "theoHisto";
name += i;
theoHisto = new TH1F(name.Data(), name.Data(), noOfChannels,
-timeResolution/2.0, (noOfChannels+0.5)*timeResolution);
theoHisto = new TH1F(name.Data(), name.Data(), noOfChannels+1,
-timeResolution*0.5, (noOfChannels+0.5)*timeResolution);
if (i < 10)
name = "hDecay0";
else
name = "hDecay";
name += i;
fakeHisto = new TH1F(name.Data(), name.Data(), noOfChannels,
-timeResolution/2.0, (noOfChannels+0.5)*timeResolution);
fakeHisto = new TH1F(name.Data(), name.Data(), noOfChannels+1,
-timeResolution*0.5, (noOfChannels+0.5)*timeResolution);
// fill theoHisto
for (Int_t j=0; j<noOfChannels; j++)
theoHisto->SetBinContent(j, histo[i][j]);
theoHisto->SetBinContent(j+1, histo[i][j]);
// fill fakeHisto
fakeHisto->FillRandom(theoHisto, (Int_t)theoHisto->Integral());

3
src/tests/skewedGaussianTest/fakeDataNemu.cpp
View File

@ -82,7 +82,8 @@ void fakeDataNemuSyntax()
int main(int argc, char *argv[])
{
const Double_t gamma_mu = TMath::TwoPi() * 0.1355; // in (rad/ns/T)
const Double_t tau_mu = 2197.147; // muon life time in ns
// muon life time in (us), see PRL99, 032001 (2007)
const Double_t tau_mu = 2197.019; // muon life time in ns
if ((argc != 4) && (argc !=6)) {
fakeDataNemuSyntax();

2
src/tests/skewedGaussianTest/paramInput.dat
View File

@ -13,7 +13,7 @@ t0, 3300, 3300, 3300, 3300
asym, 0.24, 0.24, 0.24, 0.24
#------------------------------------------------------
# phases in (°)
phase, 10.0, 90.0, 170.0, 270.0
phase, 0.0, 90.0, 180.0, 270.0
#------------------------------------------------------
# N0's
N0, 1100, 1000, 1002, 990

14
src/tests/skewedGaussianTest/skewedGaussian.C
View File

@ -41,18 +41,18 @@ void skewedGaussian()
FILE *fp;
char fln[256];
const Double_t w = 0.8; // weight of the skewed Gaussian
const Double_t B0 = 142.0; // skewed Gaussian B0 (G)
const Double_t sm = 2.5; // skewed Gaussian sigma- (G)
const Double_t sp = 4.5; // skewed Gaussian sigma+ (G)
const Double_t w = 1.0; // weight of the skewed Gaussian
const Double_t B0 = 30.0; // skewed Gaussian B0 (G)
const Double_t sm = 2.0; // skewed Gaussian sigma- (G)
const Double_t sp = 2.0; // skewed Gaussian sigma+ (G)
const Double_t B0ext = 150.0; // external field Gaussian B0 (G)
const Double_t sext = 1.2; // external field Gaussian sigma (G)
const Double_t B0ext = 30.0; // external field Gaussian B0 (G)
const Double_t sext = 10; // external field Gaussian sigma (G)
sprintf(fln, "skewedGauss-B%0.2lf-sm%0.2lf-sp%0.2lf-w%0.1lf-Bext%0.2lf-sext%0.2lf.dat",
B0, sm, sp, w, B0ext, sext);
const Int_t noOfPoints = 2000;
const Int_t noOfPoints = 8000;
const Double_t res = 0.1;