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

some changes related to the bin data plotting

Browse Source
This commit is contained in:
nemu
2008-04-05 19:14:00 +00:00
parent b47e5a4672
commit c463cb2530
8 changed files with 111 additions and 72 deletions
Download Patch File Download Diff File Expand all files Collapse all files

46
src/classes/PRunAsymmetry.cpp
View File

@ -112,7 +112,7 @@ PRunAsymmetry::PRunAsymmetry(PMsrHandler *msrInfo, PRunDataHandler *rawData, uns
//cout << endl << ">> PRunAsymmetry::PRunAsymmetry(): fAlphaBetaTag = " << fAlphaBetaTag;
// calculate fData
// calculate fFitData
if (!PrepareData())
fValid = false;
}
@ -153,34 +153,34 @@ double PRunAsymmetry::CalcChiSquare(const std::vector<double>& par)
}
// calculate chisq
for (unsigned int i=0; i<fData.fValue.size(); i++) {
if ((fData.fTime[i]>=fFitStartTime) && (fData.fTime[i]<=fFitStopTime)) {
for (unsigned int i=0; i<fFitData.fValue.size(); i++) {
if ((fFitData.fTime[i]>=fFitStartTime) && (fFitData.fTime[i]<=fFitStopTime)) {
switch (fAlphaBetaTag) {
case 1: // alpha == 1, beta == 1
asymFcnValue = fTheory->Func(fData.fTime[i], par, fFuncValues);
asymFcnValue = fTheory->Func(fFitData.fTime[i], par, fFuncValues);
break;
case 2: // alpha != 1, beta == 1
a = par[fRunInfo->fAlphaParamNo-1];
f = fTheory->Func(fData.fTime[i], par, fFuncValues);
f = fTheory->Func(fFitData.fTime[i], par, fFuncValues);
asymFcnValue = (f*(a+1.0)-(a-1.0))/((a+1.0)-f*(a-1.0));
break;
case 3: // alpha == 1, beta != 1
b = par[fRunInfo->fBetaParamNo-1];
f = fTheory->Func(fData.fTime[i], par, fFuncValues);
f = fTheory->Func(fFitData.fTime[i], par, fFuncValues);
asymFcnValue = f*(b+1.0)/(2.0-f*(b-1.0));
break;
case 4: // alpha != 1, beta != 1
a = par[fRunInfo->fAlphaParamNo-1];
b = par[fRunInfo->fBetaParamNo-1];
f = fTheory->Func(fData.fTime[i], par, fFuncValues);
f = fTheory->Func(fFitData.fTime[i], par, fFuncValues);
asymFcnValue = (f*(a*b+1.0)-(a-1.0))/((a+1.0)-f*(a*b-1.0));
break;
default:
break;
}
//if (i==0) cout << endl << "A(0) = " << asymFcnValue;
diff = fData.fValue[i] - asymFcnValue;
chisq += diff*diff / (fData.fError[i]*fData.fError[i]);
diff = fFitData.fValue[i] - asymFcnValue;
chisq += diff*diff / (fFitData.fError[i]*fFitData.fError[i]);
}
}
@ -227,32 +227,32 @@ void PRunAsymmetry::CalcTheory()
// calculate asymmetry
double asymFcnValue = 0.0;
double a, b, f;
for (unsigned int i=0; i<fData.fTime.size(); i++) {
for (unsigned int i=0; i<fFitData.fTime.size(); i++) {
switch (fAlphaBetaTag) {
case 1: // alpha == 1, beta == 1
asymFcnValue = fTheory->Func(fData.fTime[i], par, fFuncValues);
asymFcnValue = fTheory->Func(fFitData.fTime[i], par, fFuncValues);
break;
case 2: // alpha != 1, beta == 1
a = par[fRunInfo->fAlphaParamNo-1];
f = fTheory->Func(fData.fTime[i], par, fFuncValues);
f = fTheory->Func(fFitData.fTime[i], par, fFuncValues);
asymFcnValue = (f*(a+1.0)-(a-1.0))/((a+1.0)-f*(a-1.0));
break;
case 3: // alpha == 1, beta != 1
b = par[fRunInfo->fBetaParamNo-1];
f = fTheory->Func(fData.fTime[i], par, fFuncValues);
f = fTheory->Func(fFitData.fTime[i], par, fFuncValues);
asymFcnValue = f*(b+1.0)/(2.0-f*(b-1.0));
break;
case 4: // alpha != 1, beta != 1
a = par[fRunInfo->fAlphaParamNo-1];
b = par[fRunInfo->fBetaParamNo-1];
f = fTheory->Func(fData.fTime[i], par, fFuncValues);
f = fTheory->Func(fFitData.fTime[i], par, fFuncValues);
asymFcnValue = (f*(a*b+1.0)-(a-1.0))/((a+1.0)-f*(a*b-1.0));
break;
default:
asymFcnValue = 0.0;
break;
}
fData.fTheory.push_back(asymFcnValue);
fFitData.fTheory.push_back(asymFcnValue);
}
// clean up
@ -272,7 +272,7 @@ void PRunAsymmetry::CalcTheory()
*/
bool PRunAsymmetry::PrepareData()
{
//cout << endl << "in PRunAsymmetry::PrepareData(): will feed fData";
//cout << endl << "in PRunAsymmetry::PrepareData(): will feed fFitData";
// get forward/backward histo from PRunDataHandler object ------------------------
// get the correct run
@ -459,7 +459,7 @@ bool PRunAsymmetry::PrepareData()
cout << endl << " forward/backward time are not equal! This cannot be handled";
return false;
}
fData.fTime.push_back(forwardPacked.fTime[i]);
fFitData.fTime.push_back(forwardPacked.fTime[i]);
// to make the formulae more readable
f = forwardPacked.fValue[i];
b = backwardPacked.fValue[i];
@ -470,19 +470,19 @@ bool PRunAsymmetry::PrepareData()
asym = (f-b) / (f+b);
else
asym = 0.0;
fData.fValue.push_back(asym);
fFitData.fValue.push_back(asym);
// calculate the error
if (f+b != 0.0)
error = 2.0/((f+b)*(f+b))*TMath::Sqrt(b*b*ef*ef+eb*eb*f*f);
else
error = 1.0;
fData.fError.push_back(error);
fFitData.fError.push_back(error);
}
/*
FILE *fp = fopen("asym.dat", "w");
for (unsigned int i=0; i<fData.fTime.size(); i++) {
fprintf(fp, "%lf, %lf, %lf\n", fData.fTime[i], fData.fValue[i], fData.fError[i]);
for (unsigned int i=0; i<fFitData.fTime.size(); i++) {
fprintf(fp, "%lf, %lf, %lf\n", fFitData.fTime[i], fFitData.fValue[i], fFitData.fError[i]);
}
fclose(fp);
return false;
@ -490,8 +490,8 @@ return false;
// count the number of bins to be fitted
fNoOfFitBins=0;
for (unsigned int i=0; i<fData.fValue.size(); i++) {
if ((fData.fTime[i] >= fFitStartTime) && (fData.fTime[i] <= fFitStopTime))
for (unsigned int i=0; i<fFitData.fValue.size(); i++) {
if ((fFitData.fTime[i] >= fFitStartTime) && (fFitData.fTime[i] <= fFitStopTime))
fNoOfFitBins++;
}

13
src/classes/PRunBase.cpp
View File

@ -122,10 +122,15 @@ PRunBase::~PRunBase()
{
fParamNo.clear();
fData.fTime.clear();
fData.fValue.clear();
fData.fError.clear();
fData.fTheory.clear();
fFitData.fTime.clear();
fFitData.fValue.clear();
fFitData.fError.clear();
fFitData.fTheory.clear();
fBinData.fTime.clear();
fBinData.fValue.clear();
fBinData.fError.clear();
fBinData.fTheory.clear();
fT0s.clear();

8
src/classes/PRunNonMusr.cpp
View File

@ -60,7 +60,7 @@ PRunNonMusr::PRunNonMusr(PMsrHandler *msrInfo, PRunDataHandler *rawData, unsigne
{
bool success;
// calculate fData
// calculate fFitData
if (success) {
success = PrepareData();
}
@ -130,12 +130,12 @@ bool PRunNonMusr::PrepareData()
{
bool success = true;
cout << endl << "in PRunNonMusr::PrepareData(): will feed fData";
cout << endl << "in PRunNonMusr::PrepareData(): will feed fFitData";
// count the number of bins to be fitted
fNoOfFitBins=0;
for (unsigned int i=0; i<fData.fValue.size(); i++) {
if ((fData.fTime[i] >= fFitStartTime) && (fData.fTime[i] <= fFitStopTime))
for (unsigned int i=0; i<fFitData.fValue.size(); i++) {
if ((fFitData.fTime[i] >= fFitStartTime) && (fFitData.fTime[i] <= fFitStopTime))
fNoOfFitBins++;
}

8
src/classes/PRunRRF.cpp
View File

@ -60,7 +60,7 @@ PRunRRF::PRunRRF(PMsrHandler *msrInfo, PRunDataHandler *rawData, unsigned int ru
{
bool success;
// calculate fData
// calculate fFitData
if (success) {
success = PrepareData();
}
@ -130,12 +130,12 @@ bool PRunRRF::PrepareData()
{
bool success = true;
cout << endl << "in PRunRRF::PrepareData(): will feed fData";
cout << endl << "in PRunRRF::PrepareData(): will feed fFitData";
// count the number of bins to be fitted
fNoOfFitBins=0;
for (unsigned int i=0; i<fData.fValue.size(); i++) {
if ((fData.fTime[i] >= fFitStartTime) && (fData.fTime[i] <= fFitStopTime))
for (unsigned int i=0; i<fFitData.fValue.size(); i++) {
if ((fFitData.fTime[i] >= fFitStartTime) && (fFitData.fTime[i] <= fFitStopTime))
fNoOfFitBins++;
}

74
src/classes/PRunSingleHisto.cpp
View File

@ -127,11 +127,11 @@ double PRunSingleHisto::CalcChiSquare(const std::vector<double>& par)
}
// calculate chi square
for (unsigned int i=0; i<fData.fValue.size(); i++) {
if ((fData.fTime[i]>=fFitStartTime) && (fData.fTime[i]<=fFitStopTime)) {
diff = fData.fValue[i] -
(N0*TMath::Exp(-fData.fTime[i]/tau)*(1+fTheory->Func(fData.fTime[i], par, fFuncValues))+bkg);
chisq += diff*diff / (fData.fError[i]*fData.fError[i]);
for (unsigned int i=0; i<fFitData.fValue.size(); i++) {
if ((fFitData.fTime[i]>=fFitStartTime) && (fFitData.fTime[i]<=fFitStopTime)) {
diff = fFitData.fValue[i] -
(N0*TMath::Exp(-fFitData.fTime[i]/tau)*(1+fTheory->Func(fFitData.fTime[i], par, fFuncValues))+bkg);
chisq += diff*diff / (fFitData.fError[i]*fFitData.fError[i]);
}
}
@ -139,15 +139,15 @@ double PRunSingleHisto::CalcChiSquare(const std::vector<double>& par)
// static int counter = 0;
// TString fln=fRunInfo->fRunName+"_"+(Long_t)fRunInfo->fForwardHistoNo+"_data.dat";
// ofstream f(fln.Data(),ios_base::out);
// for (unsigned int i=0; i<fData.fValue.size(); i++) {
// f << endl << fData.fTime[i] << " " << fData.fValue[i] << " " << fData.fError[i];
// for (unsigned int i=0; i<fFitData.fValue.size(); i++) {
// f << endl << fFitData.fTime[i] << " " << fFitData.fValue[i] << " " << fFitData.fError[i];
// }
// f.close();
//
// fln=fRunInfo->fRunName+"_"+(Long_t)fRunInfo->fForwardHistoNo+"_theo.dat";
// ofstream ft(fln.Data(),ios_base::out);
// for (unsigned int i=0; i<fData.fValue.size(); i++) {
// ft << endl << fData.fTime[i] << " " << N0*TMath::Exp(-fData.fTime[i]/tau)*(1+fTheory->Func(fData.fTime[i], par))+bkg;
// for (unsigned int i=0; i<fFitData.fValue.size(); i++) {
// ft << endl << fFitData.fTime[i] << " " << N0*TMath::Exp(-fFitData.fTime[i]/tau)*(1+fTheory->Func(fFitData.fTime[i], par))+bkg;
// }
// ft.close();
// counter++;
@ -200,13 +200,13 @@ double PRunSingleHisto::CalcMaxLikelihood(const std::vector<double>& par)
// calculate maximum log likelihood
double theo;
double data;
for (unsigned int i=0; i<fData.fValue.size(); i++) {
if ((fData.fTime[i]>=fFitStartTime) && (fData.fTime[i]<=fFitStopTime)) {
for (unsigned int i=0; i<fFitData.fValue.size(); i++) {
if ((fFitData.fTime[i]>=fFitStartTime) && (fFitData.fTime[i]<=fFitStopTime)) {
// calculate theory for the given parameter set
theo = N0*TMath::Exp(-fData.fTime[i]/tau)*(1+fTheory->Func(fData.fTime[i], par, fFuncValues))+bkg;
theo = N0*TMath::Exp(-fFitData.fTime[i]/tau)*(1+fTheory->Func(fFitData.fTime[i], par, fFuncValues))+bkg;
// check if data value is not too small
if (fData.fValue[i] > 1.0e-9)
data = fData.fValue[i];
if (fFitData.fValue[i] > 1.0e-9)
data = fFitData.fValue[i];
else
data = 1.0e-9;
// add maximum log likelihood contribution of bin i
@ -251,8 +251,8 @@ void PRunSingleHisto::CalcTheory()
}
// calculate theory
for (unsigned int i=0; i<fData.fTime.size(); i++) {
fData.fTheory.push_back(N0*TMath::Exp(-fData.fTime[i]/tau)*(1+fTheory->Func(fData.fTime[i], par, fFuncValues))+bkg);
for (unsigned int i=0; i<fFitData.fTime.size(); i++) {
fFitData.fTheory.push_back(N0*TMath::Exp(-fFitData.fTime[i]/tau)*(1+fTheory->Func(fFitData.fTime[i], par, fFuncValues))+bkg);
}
// clean up
@ -268,7 +268,7 @@ void PRunSingleHisto::CalcTheory()
*/
bool PRunSingleHisto::PrepareData()
{
// cout << endl << "in PRunSingleHisto::PrepareData(): will feed fData";
// cout << endl << "in PRunSingleHisto::PrepareData(): will feed fFitData";
// get the proper run
PRawRunData* runData = fRawData->GetRunData(fRunInfo->fRunName);
@ -333,7 +333,7 @@ bool PRunSingleHisto::PrepareData()
// first get start data, end data, and t0
unsigned int start = fRunInfo->fDataRange[0];
unsigned int end = fRunInfo->fDataRange[1];
double t0 = fT0s[0];
double t0 = fT0s[0];
// check if start, end, and t0 make any sense
// 1st check if start and end are in proper order
if (end < start) { // need to swap them
@ -356,7 +356,8 @@ bool PRunSingleHisto::PrepareData()
cout << endl << "PRunSingleHisto::PrepareData(): t0 data bin doesn't make any sense!";
return false;
}
// everything looks fine, hence fill data set
// everything looks fine, hence fill fit data set
double value = 0.0;
for (unsigned i=start; i<end; i++) {
if (((i-start) % fRunInfo->fPacking == 0) && (i != start)) { // fill data
@ -364,12 +365,12 @@ bool PRunSingleHisto::PrepareData()
// the value is normalize to per bin
value /= fRunInfo->fPacking;
// time shifted so that packing is included correctly, i.e. t0 == t0 after packing
fData.fTime.push_back(fTimeResolution*((double)i-(double)t0-(double)fRunInfo->fPacking));
fData.fValue.push_back(value);
fFitData.fTime.push_back(fTimeResolution*((double)i-(double)t0-(double)fRunInfo->fPacking));
fFitData.fValue.push_back(value);
if (value == 0.0)
fData.fError.push_back(1.0);
fFitData.fError.push_back(1.0);
else
fData.fError.push_back(TMath::Sqrt(value));
fFitData.fError.push_back(TMath::Sqrt(value));
value = 0.0;
}
value += runData->fDataBin[histoNo][i];
@ -377,11 +378,34 @@ bool PRunSingleHisto::PrepareData()
// count the number of bins to be fitted
fNoOfFitBins=0;
for (unsigned int i=0; i<fData.fValue.size(); i++) {
if ((fData.fTime[i] >= fFitStartTime) && (fData.fTime[i] <= fFitStopTime))
for (unsigned int i=0; i<fFitData.fValue.size(); i++) {
if ((fFitData.fTime[i] >= fFitStartTime) && (fFitData.fTime[i] <= fFitStopTime))
fNoOfFitBins++;
}
// fill the bin data set (used for plots and t0's search
start = (unsigned int)t0-fRunInfo->fPacking*((unsigned int)t0/fRunInfo->fPacking);
end = (unsigned int)t0+fRunInfo->fPacking*((runData->fDataBin[histoNo].size()-1-(unsigned int)t0)/fRunInfo->fPacking);
value = 0.0;
//cout << endl << ">> start = " << start << ", end = " << end << ", " << runData->fDataBin[histoNo].size();
for (unsigned i=start; i<end; i++) {
if (((i-start) % fRunInfo->fPacking == 0) && (i != start)) { // 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;
// time shifted so that packing is included correctly, i.e. t0 == t0 after packing
fBinData.fTime.push_back(fTimeResolution*((double)i-(double)t0-(double)fRunInfo->fPacking));
fBinData.fValue.push_back(value);
if (value == 0.0)
fBinData.fError.push_back(1.0);
else
fBinData.fError.push_back(TMath::Sqrt(value));
value = 0.0;
}
value += runData->fDataBin[histoNo][i];
}
//cout << endl << "fBinData.fValue.size() = " << fBinData.fValue.size();
return true;
}