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

slightly more general handling of data. Needed to implement musrview, ...

Browse Source
This commit is contained in:
nemu 2008-04-07 10:57:42 +00:00
parent c463cb2530
commit a5d079d0a1
14 changed files with 330 additions and 259 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

346
src/classes/PRunAsymmetry.cpp
View File

@ -59,7 +59,7 @@ PRunAsymmetry::PRunAsymmetry() : PRunBase()
* \param msrInfo pointer to the msr info structure
* \param runNo number of the run of the msr-file
*/
PRunAsymmetry::PRunAsymmetry(PMsrHandler *msrInfo, PRunDataHandler *rawData, unsigned int runNo) : PRunBase(msrInfo, rawData, runNo)
PRunAsymmetry::PRunAsymmetry(PMsrHandler *msrInfo, PRunDataHandler *rawData, unsigned int runNo, EPMusrHandleTag tag) : PRunBase(msrInfo, rawData, runNo, tag)
{
// check if alpha and/or beta is fixed --------------------
@ -112,7 +112,7 @@ PRunAsymmetry::PRunAsymmetry(PMsrHandler *msrInfo, PRunDataHandler *rawData, uns
//cout << endl << ">> PRunAsymmetry::PRunAsymmetry(): fAlphaBetaTag = " << fAlphaBetaTag;
// calculate fFitData
// calculate fData
if (!PrepareData())
fValid = false;
}
@ -153,34 +153,36 @@ double PRunAsymmetry::CalcChiSquare(const std::vector<double>& par)
}
// calculate chisq
for (unsigned int i=0; i<fFitData.fValue.size(); i++) {
if ((fFitData.fTime[i]>=fFitStartTime) && (fFitData.fTime[i]<=fFitStopTime)) {
double time;
for (unsigned int i=0; i<fData.fValue.size(); i++) {
time = fData.fDataTimeStart + (double)i*fData.fDataTimeStep;
if ((time>=fFitStartTime) && (time<=fFitStopTime)) {
switch (fAlphaBetaTag) {
case 1: // alpha == 1, beta == 1
asymFcnValue = fTheory->Func(fFitData.fTime[i], par, fFuncValues);
asymFcnValue = fTheory->Func(time, par, fFuncValues);
break;
case 2: // alpha != 1, beta == 1
a = par[fRunInfo->fAlphaParamNo-1];
f = fTheory->Func(fFitData.fTime[i], par, fFuncValues);
f = fTheory->Func(time, 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(fFitData.fTime[i], par, fFuncValues);
f = fTheory->Func(time, 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(fFitData.fTime[i], par, fFuncValues);
f = fTheory->Func(time, 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 = fFitData.fValue[i] - asymFcnValue;
chisq += diff*diff / (fFitData.fError[i]*fFitData.fError[i]);
diff = fData.fValue[i] - asymFcnValue;
chisq += diff*diff / (fData.fError[i]*fData.fError[i]);
}
}
@ -227,32 +229,34 @@ void PRunAsymmetry::CalcTheory()
// calculate asymmetry
double asymFcnValue = 0.0;
double a, b, f;
for (unsigned int i=0; i<fFitData.fTime.size(); i++) {
double time;
for (unsigned int i=0; i<fData.fValue.size(); i++) {
time = fData.fDataTimeStart + (double)i*fData.fDataTimeStep;
switch (fAlphaBetaTag) {
case 1: // alpha == 1, beta == 1
asymFcnValue = fTheory->Func(fFitData.fTime[i], par, fFuncValues);
asymFcnValue = fTheory->Func(time, par, fFuncValues);
break;
case 2: // alpha != 1, beta == 1
a = par[fRunInfo->fAlphaParamNo-1];
f = fTheory->Func(fFitData.fTime[i], par, fFuncValues);
f = fTheory->Func(time, 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(fFitData.fTime[i], par, fFuncValues);
f = fTheory->Func(time, 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(fFitData.fTime[i], par, fFuncValues);
f = fTheory->Func(time, 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;
}
fFitData.fTheory.push_back(asymFcnValue);
fData.fTheory.push_back(asymFcnValue);
}
// clean up
@ -272,7 +276,7 @@ void PRunAsymmetry::CalcTheory()
*/
bool PRunAsymmetry::PrepareData()
{
//cout << endl << "in PRunAsymmetry::PrepareData(): will feed fFitData";
//cout << endl << "in PRunAsymmetry::PrepareData(): will feed fData";
// get forward/backward histo from PRunDataHandler object ------------------------
// get the correct run
@ -366,148 +370,21 @@ bool PRunAsymmetry::PrepareData()
return false;
}
// 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
unsigned int start[2] = {fRunInfo->fDataRange[0], fRunInfo->fDataRange[2]};
unsigned int end[2] = {fRunInfo->fDataRange[1], fRunInfo->fDataRange[3]};
double t0[2] = {fT0s[0], fT0s[1]};
// check if start, end, and t0 make any sense
// 1st check if start and end are in proper order
for (unsigned int i=0; i<2; i++) {
if (end[i] < start[i]) { // need to swap them
int keep = end[i];
end[i] = start[i];
start[i] = keep;
}
// 2nd check if start is within proper bounds
if ((start[i] < 0) || (start[i] > runData->fDataBin[histoNo[i]].size())) {
cout << endl << "PRunAsymmetry::PrepareData(): start data bin doesn't make any sense!";
return false;
}
// 3rd check if end is within proper bounds
if ((end[i] < 0) || (end[i] > runData->fDataBin[histoNo[i]].size())) {
cout << endl << "PRunAsymmetry::PrepareData(): end data bin doesn't make any sense!";
return false;
}
// 4th check if t0 is within proper bounds
if ((t0[i] < 0) || (t0[i] > runData->fDataBin[histoNo[i]].size())) {
cout << endl << "PRunAsymmetry::PrepareData(): t0 data bin doesn't make any sense!";
return false;
}
// everything looks fine, hence fill data set
bool status;
switch(fHandleTag) {
case kFit:
status = PrepareFitData(runData, histoNo);
break;
case kView:
status = PrepareViewData(runData, histoNo);
break;
default:
status = false;
break;
}
// everything looks fine, hence fill packed forward and backward histo
PRunData forwardPacked;
PRunData backwardPacked;
double value = 0.0;
double error = 0.0;
// forward
for (unsigned 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;
// time shifted so that packing is included correctly, i.e. t0 == t0 after packing
forwardPacked.fTime.push_back(fTimeResolution*((double)i-(double)t0[0]-(double)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 (unsigned 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;
// time shifted so that packing is included correctly, i.e. t0 == t0 after packing
backwardPacked.fTime.push_back(fTimeResolution*((double)i-(double)t0[1]-(double)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];
}
// check if packed forward and backward hist have the same size, otherwise something is wrong
if (forwardPacked.fTime.size() != backwardPacked.fTime.size()) {
cout << endl << "PRunAsymmetry::PrepareData(): **PANIC ERROR**:";
cout << endl << " packed forward and backward histo should have the same number of bins!";
cout << endl << " however found (f/b) : " << forwardPacked.fTime.size() << "/" << backwardPacked.fTime.size();
return false;
}
// form asymmetry including error propagation
double asym;
double f, b, ef, eb;
for (unsigned int i=0; i<forwardPacked.fTime.size(); i++) {
// check that forward time == backward time!!
if (forwardPacked.fTime[i] != backwardPacked.fTime[i]) {
cout << endl << "PRunAsymmetry::PrepareData(): **PANIC ERROR**:";
cout << endl << " forward/backward time are not equal! This cannot be handled";
return false;
}
fFitData.fTime.push_back(forwardPacked.fTime[i]);
// to make the formulae more readable
f = forwardPacked.fValue[i];
b = backwardPacked.fValue[i];
ef = forwardPacked.fError[i];
eb = backwardPacked.fError[i];
// check that there are indeed bins
if (f+b != 0.0)
asym = (f-b) / (f+b);
else
asym = 0.0;
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;
fFitData.fError.push_back(error);
}
/*
FILE *fp = fopen("asym.dat", "w");
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;
*/
// count the number of bins to be fitted
fNoOfFitBins=0;
for (unsigned int i=0; i<fFitData.fValue.size(); i++) {
if ((fFitData.fTime[i] >= fFitStartTime) && (fFitData.fTime[i] <= fFitStopTime))
fNoOfFitBins++;
}
// clean up
forwardPacked.fTime.clear();
forwardPacked.fValue.clear();
forwardPacked.fError.clear();
backwardPacked.fTime.clear();
backwardPacked.fValue.clear();
backwardPacked.fError.clear();
fForward.clear();
fForwardErr.clear();
fBackward.clear();
fBackwardErr.clear();
return true;
return status;
}
//--------------------------------------------------------------------------
@ -640,4 +517,159 @@ bool PRunAsymmetry::SubtractEstimatedBkg()
return true;
}
//--------------------------------------------------------------------------
// PrepareFitData
//--------------------------------------------------------------------------
/**
* <p>
*
*/
bool PRunAsymmetry::PrepareFitData(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
unsigned int start[2] = {fRunInfo->fDataRange[0], fRunInfo->fDataRange[2]};
unsigned int end[2] = {fRunInfo->fDataRange[1], fRunInfo->fDataRange[3]};
double t0[2] = {fT0s[0], fT0s[1]};
// check if start, end, and t0 make any sense
// 1st check if start and end are in proper order
for (unsigned int i=0; i<2; i++) {
if (end[i] < start[i]) { // need to swap them
int keep = end[i];
end[i] = start[i];
start[i] = keep;
}
// 2nd check if start is within proper bounds
if ((start[i] < 0) || (start[i] > runData->fDataBin[histoNo[i]].size())) {
cout << endl << "PRunAsymmetry::PrepareData(): start data bin doesn't make any sense!";
return false;
}
// 3rd check if end is within proper bounds
if ((end[i] < 0) || (end[i] > runData->fDataBin[histoNo[i]].size())) {
cout << endl << "PRunAsymmetry::PrepareData(): end data bin doesn't make any sense!";
return false;
}
// 4th check if t0 is within proper bounds
if ((t0[i] < 0) || (t0[i] > runData->fDataBin[histoNo[i]].size())) {
cout << endl << "PRunAsymmetry::PrepareData(): t0 data bin doesn't make any sense!";
return false;
}
}
// everything looks fine, hence fill packed forward and backward histo
PRunData forwardPacked;
PRunData backwardPacked;
double value = 0.0;
double error = 0.0;
// forward
for (unsigned 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;
}
value += fForward[i];
error += fForwardErr[i]*fForwardErr[i];
}
// backward
for (unsigned 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;
}
value += fBackward[i];
error += fBackwardErr[i]*fBackwardErr[i];
}
// check if packed forward and backward hist have the same size, otherwise something is wrong
if (forwardPacked.fValue.size() != backwardPacked.fValue.size()) {
cout << endl << "PRunAsymmetry::PrepareData(): **PANIC ERROR**:";
cout << endl << " packed forward and backward histo should have the same number of bins!";
cout << endl << " however found (f/b) : " << forwardPacked.fValue.size() << "/" << backwardPacked.fValue.size();
return false;
}
// form asymmetry including error propagation
double asym;
double f, b, ef, eb;
// fill data time start, and step
fData.fDataTimeStart = fTimeResolution*((double)fRunInfo->fPacking/2.0);
fData.fDataTimeStep = fTimeResolution*(double)fRunInfo->fPacking;
for (unsigned int i=0; i<forwardPacked.fValue.size(); i++) {
// to make the formulae more readable
f = forwardPacked.fValue[i];
b = backwardPacked.fValue[i];
ef = forwardPacked.fError[i];
eb = backwardPacked.fError[i];
// check that there are indeed bins
if (f+b != 0.0)
asym = (f-b) / (f+b);
else
asym = 0.0;
fData.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);
}
/*
FILE *fp = fopen("asym.dat", "w");
for (unsigned int i=0; i<fData.fValue.size(); i++) {
fprintf(fp, "%lf, %lf, %lf\n", fData.fDataTimeStart+(double)i*fData.fDataTimeStep, fData.fValue[i], fData.fError[i]);
}
fclose(fp);
return false;
*/
// count the number of bins to be fitted
double time;
fNoOfFitBins=0;
for (unsigned int i=0; i<fData.fValue.size(); i++) {
time = fData.fDataTimeStart + (double)i * fData.fDataTimeStep;
if ((time >= fFitStartTime) && (time <= fFitStopTime))
fNoOfFitBins++;
}
// clean up
forwardPacked.fValue.clear();
forwardPacked.fError.clear();
backwardPacked.fValue.clear();
backwardPacked.fError.clear();
fForward.clear();
fForwardErr.clear();
fBackward.clear();
fBackwardErr.clear();
return true;
}
//--------------------------------------------------------------------------
// PrepareViewData
//--------------------------------------------------------------------------
/**
* <p>
*
*/
bool PRunAsymmetry::PrepareViewData(PRawRunData* runData, unsigned int histoNo[2])
{
return true;
}

16
src/classes/PRunBase.cpp
View File

@ -58,6 +58,7 @@ PRunBase::PRunBase()
fTimeResolution = -1.0;
fValid = true;
fHandleTag = kEmpty;
}
//--------------------------------------------------------------------------
@ -69,9 +70,10 @@ PRunBase::PRunBase()
* \param msrInfo pointer to the msr info structure
* \param runNo number of the run of the msr-file
*/
PRunBase::PRunBase(PMsrHandler *msrInfo, PRunDataHandler *rawData, unsigned int runNo)
PRunBase::PRunBase(PMsrHandler *msrInfo, PRunDataHandler *rawData, unsigned int runNo, EPMusrHandleTag tag)
{
fValid = true;
fHandleTag = tag;
fRunNo = static_cast<int>(runNo);
if ((runNo < 0) || (runNo > msrInfo->GetMsrRunList()->size())) {
@ -122,15 +124,9 @@ PRunBase::~PRunBase()
{
fParamNo.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();
fData.fValue.clear();
fData.fError.clear();
fData.fTheory.clear();
fT0s.clear();

11
src/classes/PRunListCollection.cpp
View File

@ -86,8 +86,9 @@ PRunListCollection::~PRunListCollection()
* <p>
*
* \param runNo
* \param tag
*/
bool PRunListCollection::Add(int runNo)
bool PRunListCollection::Add(int runNo, EPMusrHandleTag tag)
{
bool success = true;
@ -101,22 +102,22 @@ bool PRunListCollection::Add(int runNo)
switch (fitType) {
case PRUN_SINGLE_HISTO:
fRunSingleHistoList.push_back(new PRunSingleHisto(fMsrInfo, fData, runNo));
fRunSingleHistoList.push_back(new PRunSingleHisto(fMsrInfo, fData, runNo, tag));
if (!fRunSingleHistoList[fRunSingleHistoList.size()-1]->IsValid())
success = false;
break;
case PRUN_ASYMMETRY:
fRunAsymmetryList.push_back(new PRunAsymmetry(fMsrInfo, fData, runNo));
fRunAsymmetryList.push_back(new PRunAsymmetry(fMsrInfo, fData, runNo, tag));
if (!fRunAsymmetryList[fRunAsymmetryList.size()-1]->IsValid())
success = false;
break;
case PRUN_RRF:
fRunRRFList.push_back(new PRunRRF(fMsrInfo, fData, runNo));
fRunRRFList.push_back(new PRunRRF(fMsrInfo, fData, runNo, tag));
if (!fRunRRFList[fRunRRFList.size()-1]->IsValid())
success = false;
break;
case PRUN_NON_MUSR:
fRunNonMusrList.push_back(new PRunNonMusr(fMsrInfo, fData, runNo));
fRunNonMusrList.push_back(new PRunNonMusr(fMsrInfo, fData, runNo, tag));
if (!fRunNonMusrList[fRunNonMusrList.size()-1]->IsValid())
success = false;
break;

11
src/classes/PRunNonMusr.cpp
View File

@ -45,6 +45,8 @@ PRunNonMusr::PRunNonMusr() : PRunBase()
fFitStartTime = 0.0;
fFitStopTime = 0.0;
fNoOfFitBins = 0;
fHandleTag = kEmpty;
}
//--------------------------------------------------------------------------
@ -56,7 +58,7 @@ PRunNonMusr::PRunNonMusr() : PRunBase()
* \param msrInfo pointer to the msr info structure
* \param runNo number of the run of the msr-file
*/
PRunNonMusr::PRunNonMusr(PMsrHandler *msrInfo, PRunDataHandler *rawData, unsigned int runNo) : PRunBase(msrInfo, rawData, runNo)
PRunNonMusr::PRunNonMusr(PMsrHandler *msrInfo, PRunDataHandler *rawData, unsigned int runNo, EPMusrHandleTag tag) : PRunBase(msrInfo, rawData, runNo, tag)
{
bool success;
@ -132,13 +134,6 @@ bool PRunNonMusr::PrepareData()
cout << endl << "in PRunNonMusr::PrepareData(): will feed fFitData";
// count the number of bins to be fitted
fNoOfFitBins=0;
for (unsigned int i=0; i<fFitData.fValue.size(); i++) {
if ((fFitData.fTime[i] >= fFitStartTime) && (fFitData.fTime[i] <= fFitStopTime))
fNoOfFitBins++;
}
return success;
}

13
src/classes/PRunRRF.cpp
View File

@ -45,6 +45,8 @@ PRunRRF::PRunRRF() : PRunBase()
fFitStartTime = 0.0;
fFitStopTime = 0.0;
fNoOfFitBins = 0;
fHandleTag = kEmpty;
}
//--------------------------------------------------------------------------
@ -56,7 +58,7 @@ PRunRRF::PRunRRF() : PRunBase()
* \param msrInfo pointer to the msr info structure
* \param runNo number of the run of the msr-file
*/
PRunRRF::PRunRRF(PMsrHandler *msrInfo, PRunDataHandler *rawData, unsigned int runNo) : PRunBase(msrInfo, rawData, runNo)
PRunRRF::PRunRRF(PMsrHandler *msrInfo, PRunDataHandler *rawData, unsigned int runNo, EPMusrHandleTag tag) : PRunBase(msrInfo, rawData, runNo, tag)
{
bool success;
@ -130,14 +132,7 @@ bool PRunRRF::PrepareData()
{
bool success = true;
cout << endl << "in PRunRRF::PrepareData(): will feed fFitData";
// count the number of bins to be fitted
fNoOfFitBins=0;
for (unsigned int i=0; i<fFitData.fValue.size(); i++) {
if ((fFitData.fTime[i] >= fFitStartTime) && (fFitData.fTime[i] <= fFitStopTime))
fNoOfFitBins++;
}
cout << endl << "in PRunRRF::PrepareData(): will feed fData";
return success;
}

134
src/classes/PRunSingleHisto.cpp
View File

@ -58,7 +58,7 @@ PRunSingleHisto::PRunSingleHisto() : PRunBase()
* \param msrInfo pointer to the msr info structure
* \param runNo number of the run of the msr-file
*/
PRunSingleHisto::PRunSingleHisto(PMsrHandler *msrInfo, PRunDataHandler *rawData, unsigned int runNo) : PRunBase(msrInfo, rawData, runNo)
PRunSingleHisto::PRunSingleHisto(PMsrHandler *msrInfo, PRunDataHandler *rawData, unsigned int runNo, EPMusrHandleTag tag) : PRunBase(msrInfo, rawData, runNo, tag)
{
if (!PrepareData()) {
cout << endl << "**SEVERE ERROR**: PRunSingleHisto::PRunSingleHisto: Couldn't prepare data for fitting!";
@ -127,11 +127,13 @@ double PRunSingleHisto::CalcChiSquare(const std::vector<double>& par)
}
// calculate chi square
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]);
double time;
for (unsigned int i=0; i<fData.fValue.size(); i++) {
time = fData.fDataTimeStart + (double)i*fData.fDataTimeStep;
if ((time>=fFitStartTime) && (time<=fFitStopTime)) {
diff = fData.fValue[i] -
(N0*TMath::Exp(-time/tau)*(1+fTheory->Func(time, par, fFuncValues))+bkg);
chisq += diff*diff / (fData.fError[i]*fData.fError[i]);
}
}
@ -139,15 +141,17 @@ 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<fFitData.fValue.size(); i++) {
// f << endl << fFitData.fTime[i] << " " << fFitData.fValue[i] << " " << fFitData.fError[i];
// for (unsigned int i=0; i<fData.fValue.size(); i++) {
// time = fData.fDataTimeStart + (double)i*fData.fDataTimeStep;
// f << endl << time << " " << fData.fValue[i] << " " << fData.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<fFitData.fValue.size(); i++) {
// ft << endl << fFitData.fTime[i] << " " << N0*TMath::Exp(-fFitData.fTime[i]/tau)*(1+fTheory->Func(fFitData.fTime[i], par))+bkg;
// for (unsigned int i=0; i<fData.fValue.size(); i++) {
// time = fData.fDataTimeStart + (double)i*fData.fDataTimeStep;
// ft << endl << time << " " << N0*TMath::Exp(-time/tau)*(1.0+fTheory->Func(time, par, fFuncValues))+bkg;
// }
// ft.close();
// counter++;
@ -200,13 +204,15 @@ double PRunSingleHisto::CalcMaxLikelihood(const std::vector<double>& par)
// calculate maximum log likelihood
double theo;
double data;
for (unsigned int i=0; i<fFitData.fValue.size(); i++) {
if ((fFitData.fTime[i]>=fFitStartTime) && (fFitData.fTime[i]<=fFitStopTime)) {
double time;
for (unsigned int i=0; i<fData.fValue.size(); i++) {
time = fData.fDataTimeStart + (double)i*fData.fDataTimeStep;
if ((time>=fFitStartTime) && (time<=fFitStopTime)) {
// calculate theory for the given parameter set
theo = N0*TMath::Exp(-fFitData.fTime[i]/tau)*(1+fTheory->Func(fFitData.fTime[i], par, fFuncValues))+bkg;
theo = N0*TMath::Exp(-time/tau)*(1+fTheory->Func(time, par, fFuncValues))+bkg;
// check if data value is not too small
if (fFitData.fValue[i] > 1.0e-9)
data = fFitData.fValue[i];
if (fData.fValue[i] > 1.0e-9)
data = fData.fValue[i];
else
data = 1.0e-9;
// add maximum log likelihood contribution of bin i
@ -251,8 +257,10 @@ void PRunSingleHisto::CalcTheory()
}
// calculate theory
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);
double time;
for (unsigned int i=0; i<fData.fValue.size(); i++) {
time = fData.fDataTimeStart + (double)i*fData.fDataTimeStep;
fData.fTheory.push_back(N0*TMath::Exp(-time/tau)*(1+fTheory->Func(time, par, fFuncValues))+bkg);
}
// clean up
@ -268,7 +276,7 @@ void PRunSingleHisto::CalcTheory()
*/
bool PRunSingleHisto::PrepareData()
{
// cout << endl << "in PRunSingleHisto::PrepareData(): will feed fFitData";
// cout << endl << "in PRunSingleHisto::PrepareData(): will feed fData";
// get the proper run
PRawRunData* runData = fRawData->GetRunData(fRunInfo->fRunName);
@ -357,20 +365,51 @@ bool PRunSingleHisto::PrepareData()
return false;
}
// everything looks fine, hence fill fit data set
// everything looks fine, hence fill data set
bool status;
switch(fHandleTag) {
case kFit:
status = PrepareFitData(start, end, t0, runData, histoNo);
break;
case kView:
status = PrepareViewData(start, end, t0, runData, histoNo);
break;
default:
status = false;
break;
}
return status;
}
//--------------------------------------------------------------------------
// PrepareFitData
//--------------------------------------------------------------------------
/**
* <p>
*
* \param start
* \param end
* \param runData
* \param histoNo
*/
bool PRunSingleHisto::PrepareFitData(unsigned int start, unsigned int end, double t0,
PRawRunData* runData, unsigned int histoNo)
{
double value = 0.0;
// time shifted so that packing is included correctly, i.e. t0 == t0 after packing
fData.fDataTimeStart = fTimeResolution*((double)fRunInfo->fPacking/2.0);
fData.fDataTimeStep = fTimeResolution*fRunInfo->fPacking;
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
fFitData.fTime.push_back(fTimeResolution*((double)i-(double)t0-(double)fRunInfo->fPacking));
fFitData.fValue.push_back(value);
fData.fValue.push_back(value);
if (value == 0.0)
fFitData.fError.push_back(1.0);
fData.fError.push_back(1.0);
else
fFitData.fError.push_back(TMath::Sqrt(value));
fData.fError.push_back(TMath::Sqrt(value));
value = 0.0;
}
value += runData->fDataBin[histoNo][i];
@ -378,35 +417,30 @@ bool PRunSingleHisto::PrepareData()
// count the number of bins to be fitted
fNoOfFitBins=0;
for (unsigned int i=0; i<fFitData.fValue.size(); i++) {
if ((fFitData.fTime[i] >= fFitStartTime) && (fFitData.fTime[i] <= fFitStopTime))
double time;
for (unsigned int i=0; i<fData.fValue.size(); i++) {
time = fData.fDataTimeStart + (double)i*fData.fDataTimeStep;
if ((time >= fFitStartTime) && (time <= 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;
}
//--------------------------------------------------------------------------
// PrepareViewData
//--------------------------------------------------------------------------
/**
* <p>
*
* \param start
* \param end
* \param runData
* \param histoNo
*/
bool PRunSingleHisto::PrepareViewData(unsigned int start, unsigned int end, double t0,
PRawRunData* runData, unsigned int histoNo)
{
// to be implemented yet
return true;
}

15
src/include/PMusr.h
View File

@ -118,12 +118,23 @@ typedef vector<TString> PStringVector;
//-------------------------------------------------------------
/**
* <p>Predominantly used in PRunBase.
* <p> data handling tag
*/
enum EPMusrHandleTag { kEmpty, kFit, kView };
//-------------------------------------------------------------
/**
* <p> Predominantly used in PRunBase.
*/
typedef struct {
vector<double> fTime;
// data related info
double fDataTimeStart;
double fDataTimeStep;
vector<double> fValue;
vector<double> fError;
// theory related info
double fTheoryTimeStart;
double fTheoryTimeStep;
vector<double> fTheory;
} PRunData;

4
src/include/PRunAsymmetry.h
View File

@ -38,7 +38,7 @@ class PRunAsymmetry : public PRunBase
{
public:
PRunAsymmetry();
PRunAsymmetry(PMsrHandler *msrInfo, PRunDataHandler *rawData, unsigned int runNo);
PRunAsymmetry(PMsrHandler *msrInfo, PRunDataHandler *rawData, unsigned int runNo, EPMusrHandleTag tag);
virtual ~PRunAsymmetry();
virtual double CalcChiSquare(const std::vector<double>& par);
@ -49,6 +49,8 @@ class PRunAsymmetry : public PRunBase
protected:
virtual bool PrepareData();
virtual bool PrepareFitData(PRawRunData* runData, unsigned int histoNo[2]);
virtual bool PrepareViewData(PRawRunData* runData, unsigned int histoNo[2]);
private:
unsigned int fAlphaBetaTag; ///< 1-> alpha = beta = 1; 2-> alpha != 1, beta = 1; 3-> alpha = 1, beta != 1; 4-> alpha != 1, beta != 1

9
src/include/PRunBase.h
View File

@ -61,7 +61,7 @@ class PRunBase
{
public:
PRunBase();
PRunBase(PMsrHandler *msrInfo, PRunDataHandler *rawData, unsigned int runNo);
PRunBase(PMsrHandler *msrInfo, PRunDataHandler *rawData, unsigned int runNo, EPMusrHandleTag tag);
virtual ~PRunBase();
virtual double CalcChiSquare(const vector<double>& par) = 0; // pure virtual, i.e. needs to be implemented by the deriving class!!
@ -70,13 +70,15 @@ class PRunBase
virtual void CalcTheory() = 0; // pure virtual, i.e. needs to be implemented by the deriving class!!
virtual unsigned int GetRunNo() { return fRunNo; }
virtual PRunData* GetData() { return &fFitData; }
virtual PRunData* GetData() { return &fData; }
virtual void CleanUp();
virtual bool IsValid() { return fValid; }
protected:
bool fValid;
EPMusrHandleTag fHandleTag; ///< tag telling whether this is used for fit, view, ...
unsigned int fRunNo; ///< number of the run within the msr file
PMsrHandler *fMsrInfo; ///< msr-file handler
PMsrRunStructure *fRunInfo; ///< run info used to filter out needed infos for the run
@ -84,8 +86,7 @@ class PRunBase
PIntVector fParamNo; ///< vector of parameter numbers for the specifc run
PRunData fFitData; ///< data to be fitted, i.e. binned data on the fit time range
PRunData fBinData; ///< binned data set, starting at raw data 0, i.e. at negative times used to plot and determine t0's
PRunData fData; ///< data to be fitted, viewed, i.e. binned data
double fTimeResolution; ///< time resolution
PDoubleVector fT0s; ///< all t0's of a run! The derived classes will handle it

2
src/include/PRunListCollection.h
View File

@ -51,7 +51,7 @@ class PRunListCollection
enum EDataSwitch { kIndex, kRunNo };
virtual bool Add(int runNo);
virtual bool Add(int runNo, EPMusrHandleTag tag);
virtual double GetSingleHistoChisq(const std::vector<double>& par);
virtual double GetAsymmetryChisq(const std::vector<double>& par);

2
src/include/PRunNonMusr.h
View File

@ -38,7 +38,7 @@ class PRunNonMusr : public PRunBase
{
public:
PRunNonMusr();
PRunNonMusr(PMsrHandler *msrInfo, PRunDataHandler *rawData, unsigned int runNo);
PRunNonMusr(PMsrHandler *msrInfo, PRunDataHandler *rawData, unsigned int runNo, EPMusrHandleTag tag);
virtual ~PRunNonMusr();
virtual double CalcChiSquare(const std::vector<double>& par);

2
src/include/PRunRRF.h
View File

@ -38,7 +38,7 @@ class PRunRRF : public PRunBase
{
public:
PRunRRF();
PRunRRF(PMsrHandler *msrInfo, PRunDataHandler *rawData, unsigned int runNo);
PRunRRF(PMsrHandler *msrInfo, PRunDataHandler *rawData, unsigned int runNo, EPMusrHandleTag tag);
virtual ~PRunRRF();
virtual double CalcChiSquare(const std::vector<double>& par);

4
src/include/PRunSingleHisto.h
View File

@ -38,7 +38,7 @@ class PRunSingleHisto : public PRunBase
{
public:
PRunSingleHisto();
PRunSingleHisto(PMsrHandler *msrInfo, PRunDataHandler *rawData, unsigned int runNo);
PRunSingleHisto(PMsrHandler *msrInfo, PRunDataHandler *rawData, unsigned int runNo, EPMusrHandleTag tag);
virtual ~PRunSingleHisto();
virtual double CalcChiSquare(const std::vector<double>& par);
@ -49,6 +49,8 @@ class PRunSingleHisto : public PRunBase
protected:
virtual bool PrepareData();
virtual bool PrepareFitData(unsigned int start, unsigned int end, double t0, PRawRunData* runData, unsigned int histoNo);
virtual bool PrepareViewData(unsigned int start, unsigned int end, double t0, PRawRunData* runData, unsigned int histoNo);
private:
double fFitStartTime;

20
src/musrfit.cpp
View File

@ -248,10 +248,12 @@ void musrfit_write_ascii(TString fln, PRunData *data, int runCounter)
}
// dump data
f << "% number of data values = " << data->fTime.size() << endl;
f << "% number of data values = " << data->fValue.size() << endl;
f << "% time (us), value, error, theory" << endl;
for (unsigned int i=0; i<data->fTime.size(); i++) {
f << data->fTime[i] << ", " << data->fValue[i] << ", " << data->fError[i] << ", " << data->fTheory[i] << endl;
double time;
for (unsigned int i=0; i<data->fValue.size(); i++) {
time = data->fDataTimeStart + (double)i*data->fDataTimeStep;
f << time << ", " << data->fValue[i] << ", " << data->fError[i] << ", " << data->fTheory[i] << endl;
}
// close file
@ -349,14 +351,14 @@ void musrfit_write_root(TFile &f, TString fln, PRunData *data, int runCounter)
TCanvas *c = new TCanvas(name, title.Data(), 10, 10, 800, 600);
// create histos
Double_t diff = data->fTime[1]-data->fTime[0];
Double_t diff = data->fDataTimeStep;
Double_t start = -diff/2.0;
Double_t end = data->fTime[data->fTime.size()-1]+diff/2.0;
TH1F *hdata = new TH1F("hdata", "run data", data->fTime.size(), start, end);
TH1F *htheo = new TH1F("htheo", "run theory", data->fTime.size(), start, end);
Double_t end = data->fDataTimeStep*(data->fValue.size()-1)+diff/2.0;
TH1F *hdata = new TH1F("hdata", "run data", data->fValue.size(), start, end);
TH1F *htheo = new TH1F("htheo", "run theory", data->fValue.size(), start, end);
// fill data
for (unsigned int i=0; i<data->fTime.size(); i++) {
for (unsigned int i=0; i<data->fValue.size(); i++) {
hdata->SetBinContent(i, data->fValue[i]);
hdata->SetBinError(i, data->fError[i]);
htheo->SetBinContent(i, data->fTheory[i]);
@ -555,7 +557,7 @@ int main(int argc, char *argv[])
// feed all the necessary histogramms for the fit
runListCollection = new PRunListCollection(msrHandler, dataHandler);
for (unsigned int i=0; i < msrHandler->GetMsrRunList()->size(); i++) {
success = runListCollection->Add(i);
success = runListCollection->Add(i, kFit);
if (!success) {
cout << endl << "Couldn't handle run no " << i << " ";
cout << (*msrHandler->GetMsrRunList())[i].fRunName.Data();