478 lines
16 KiB
C++
478 lines
16 KiB
C++
/***************************************************************************
|
|
|
|
PRunNonMusr.cpp
|
|
|
|
Author: Andreas Suter
|
|
e-mail: andreas.suter@psi.ch
|
|
|
|
$Id$
|
|
|
|
***************************************************************************/
|
|
|
|
/***************************************************************************
|
|
* Copyright (C) 2007 by Andreas Suter *
|
|
* andreas.suter@psi.ch *
|
|
* *
|
|
* This program is free software; you can redistribute it and/or modify *
|
|
* it under the terms of the GNU General Public License as published by *
|
|
* the Free Software Foundation; either version 2 of the License, or *
|
|
* (at your option) any later version. *
|
|
* *
|
|
* This program is distributed in the hope that it will be useful, *
|
|
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
|
|
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
|
|
* GNU General Public License for more details. *
|
|
* *
|
|
* You should have received a copy of the GNU General Public License *
|
|
* along with this program; if not, write to the *
|
|
* Free Software Foundation, Inc., *
|
|
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
|
|
***************************************************************************/
|
|
|
|
#include <iostream>
|
|
|
|
#include "PRunNonMusr.h"
|
|
|
|
//--------------------------------------------------------------------------
|
|
// Constructor
|
|
//--------------------------------------------------------------------------
|
|
/**
|
|
* <p>Constructor.
|
|
*/
|
|
PRunNonMusr::PRunNonMusr() : PRunBase()
|
|
{
|
|
fFitStartTime = 0.0;
|
|
fFitStopTime = 0.0;
|
|
fNoOfFitBins = 0;
|
|
|
|
fHandleTag = kEmpty;
|
|
|
|
fRawRunData = 0;
|
|
}
|
|
|
|
//--------------------------------------------------------------------------
|
|
// Constructor
|
|
//--------------------------------------------------------------------------
|
|
/**
|
|
* <p>Constructor
|
|
*
|
|
* \param msrInfo pointer to the msr-file handler
|
|
* \param rawData raw run data
|
|
* \param runNo number of the run within the msr-file
|
|
* \param tag tag showing what shall be done: kFit == fitting, kView == viewing
|
|
*/
|
|
PRunNonMusr::PRunNonMusr(PMsrHandler *msrInfo, PRunDataHandler *rawData, UInt_t runNo, EPMusrHandleTag tag) : PRunBase(msrInfo, rawData, runNo, tag)
|
|
{
|
|
// get the proper run
|
|
fRawRunData = fRawData->GetRunData(*(fRunInfo->GetRunName()));
|
|
if (!fRawRunData) { // couldn't get run
|
|
cerr << endl << "PRunNonMusr::PRunNonMusr(): **ERROR** Couldn't get raw run data!";
|
|
cerr << endl;
|
|
fValid = false;
|
|
}
|
|
|
|
// calculate fData
|
|
if (!PrepareData())
|
|
fValid = false;
|
|
}
|
|
|
|
//--------------------------------------------------------------------------
|
|
// Destructor
|
|
//--------------------------------------------------------------------------
|
|
/**
|
|
* <p>Destructor
|
|
*/
|
|
PRunNonMusr::~PRunNonMusr()
|
|
{
|
|
}
|
|
|
|
//--------------------------------------------------------------------------
|
|
// CalcChiSquare
|
|
//--------------------------------------------------------------------------
|
|
/**
|
|
* <p>Calculate chi-square.
|
|
*
|
|
* <b>return:</b>
|
|
* - chisq value
|
|
*
|
|
* \param par parameter vector iterated by minuit2
|
|
*/
|
|
Double_t PRunNonMusr::CalcChiSquare(const std::vector<Double_t>& par)
|
|
{
|
|
Double_t chisq = 0.0;
|
|
Double_t diff = 0.0;
|
|
|
|
// calculate functions
|
|
for (Int_t i=0; i<fMsrInfo->GetNoOfFuncs(); i++) {
|
|
fFuncValues[i] = fMsrInfo->EvalFunc(fMsrInfo->GetFuncNo(i), *fRunInfo->GetMap(), par);
|
|
}
|
|
|
|
// calculate chi square
|
|
Double_t x;
|
|
for (UInt_t i=0; i<fData.GetValue()->size(); i++) {
|
|
x = fData.GetX()->at(i);
|
|
if ((x>=fFitStartTime) && (x<=fFitStopTime)) {
|
|
diff = fData.GetValue()->at(i) - fTheory->Func(x, par, fFuncValues);
|
|
chisq += diff*diff / (fData.GetError()->at(i)*fData.GetError()->at(i));
|
|
}
|
|
}
|
|
|
|
return chisq;
|
|
}
|
|
|
|
//--------------------------------------------------------------------------
|
|
// CalcMaxLikelihood
|
|
//--------------------------------------------------------------------------
|
|
/**
|
|
* <p>NOT IMPLEMENTED!!
|
|
*
|
|
* \param par parameter vector iterated by minuit2
|
|
*/
|
|
Double_t PRunNonMusr::CalcMaxLikelihood(const std::vector<Double_t>& par)
|
|
{
|
|
cout << endl << "PRunSingleHisto::CalcMaxLikelihood(): not implemented yet ..." << endl;
|
|
|
|
return 1.0;
|
|
}
|
|
|
|
//--------------------------------------------------------------------------
|
|
// CalcTheory
|
|
//--------------------------------------------------------------------------
|
|
/**
|
|
* <p>Calculate theory for a given set of fit-parameters.
|
|
*/
|
|
void PRunNonMusr::CalcTheory()
|
|
{
|
|
}
|
|
|
|
//--------------------------------------------------------------------------
|
|
// PrepareData
|
|
//--------------------------------------------------------------------------
|
|
/**
|
|
* <p>Prepare data for fitting or viewing.
|
|
*
|
|
* <b>return:</b>
|
|
* - true if everthing went smooth
|
|
* - false, otherwise.
|
|
*/
|
|
Bool_t PRunNonMusr::PrepareData()
|
|
{
|
|
Bool_t success = true;
|
|
|
|
if (fRunInfo->GetRunNameSize() > 1) { // ADDRUN present which is not supported for NonMusr
|
|
cerr << endl << ">> PRunNonMusr::PrepareData(): **WARNING** ADDRUN NOT SUPPORTED FOR THIS FIT TYPE, WILL IGNORE IT." << endl;
|
|
}
|
|
|
|
if (fHandleTag == kFit)
|
|
success = PrepareFitData();
|
|
else if (fHandleTag == kView)
|
|
success = PrepareViewData();
|
|
else
|
|
success = false;
|
|
|
|
return success;
|
|
}
|
|
|
|
//--------------------------------------------------------------------------
|
|
// PrepareFitData
|
|
//--------------------------------------------------------------------------
|
|
/**
|
|
* <p>Prepare data for fitting.
|
|
*
|
|
* <b>return:</b>
|
|
* - true if everthing went smooth
|
|
* - false, otherwise.
|
|
*/
|
|
Bool_t PRunNonMusr::PrepareFitData()
|
|
{
|
|
Bool_t success = true;
|
|
|
|
// keep start/stop time for fit: here the meaning is of course start x, stop x
|
|
fFitStartTime = fRunInfo->GetFitRange(0);
|
|
fFitStopTime = fRunInfo->GetFitRange(1);
|
|
|
|
// get x-, y-index
|
|
UInt_t xIndex = GetXIndex();
|
|
UInt_t yIndex = GetYIndex();
|
|
|
|
// pack the raw data
|
|
Double_t value = 0.0;
|
|
Double_t err = 0.0;
|
|
for (UInt_t i=0; i<fRawRunData->fDataNonMusr.GetData()->at(xIndex).size(); i++) {
|
|
if (fRunInfo->GetPacking() == 1) {
|
|
fData.AppendXValue(fRawRunData->fDataNonMusr.GetData()->at(xIndex).at(i));
|
|
fData.AppendValue(fRawRunData->fDataNonMusr.GetData()->at(yIndex).at(i));
|
|
fData.AppendErrorValue(fRawRunData->fDataNonMusr.GetErrData()->at(yIndex).at(i));
|
|
} else { // packed data, i.e. fRunInfo->GetPacking() > 1
|
|
if ((i % fRunInfo->GetPacking() == 0) && (i != 0)) { // fill data
|
|
fData.AppendXValue(fRawRunData->fDataNonMusr.GetData()->at(xIndex).at(i)-(fRawRunData->fDataNonMusr.GetData()->at(xIndex).at(i)-fRawRunData->fDataNonMusr.GetData()->at(xIndex).at(i-fRunInfo->GetPacking()))/2.0);
|
|
fData.AppendValue(value);
|
|
fData.AppendErrorValue(TMath::Sqrt(err));
|
|
value = 0.0;
|
|
err = 0.0;
|
|
}
|
|
// sum raw data values
|
|
value += fRawRunData->fDataNonMusr.GetData()->at(yIndex).at(i);
|
|
err += fRawRunData->fDataNonMusr.GetErrData()->at(yIndex).at(i)*fRawRunData->fDataNonMusr.GetErrData()->at(yIndex).at(i);
|
|
}
|
|
}
|
|
|
|
// count the number of bins to be fitted
|
|
fNoOfFitBins=0;
|
|
Double_t x;
|
|
for (UInt_t i=0; i<fData.GetValue()->size(); i++) {
|
|
x = fData.GetX()->at(i);
|
|
if ((x >= fFitStartTime) && (x <= fFitStopTime))
|
|
fNoOfFitBins++;
|
|
}
|
|
|
|
return success;
|
|
}
|
|
|
|
//--------------------------------------------------------------------------
|
|
// PrepareViewData
|
|
//--------------------------------------------------------------------------
|
|
/**
|
|
* <p>Prepare data for viewing.
|
|
*
|
|
* <b>return:</b>
|
|
* - true if everthing went smooth
|
|
* - false, otherwise.
|
|
*/
|
|
Bool_t PRunNonMusr::PrepareViewData()
|
|
{
|
|
Bool_t success = true;
|
|
|
|
// get x-, y-index
|
|
UInt_t xIndex = GetXIndex();
|
|
UInt_t yIndex = GetYIndex();
|
|
|
|
// fill data histo
|
|
// pack the raw data
|
|
Double_t value = 0.0;
|
|
Double_t err = 0.0;
|
|
for (UInt_t i=0; i<fRawRunData->fDataNonMusr.GetData()->at(xIndex).size(); i++) {
|
|
if (fRunInfo->GetPacking() == 1) {
|
|
fData.AppendXValue(fRawRunData->fDataNonMusr.GetData()->at(xIndex).at(i));
|
|
fData.AppendValue(fRawRunData->fDataNonMusr.GetData()->at(yIndex).at(i));
|
|
fData.AppendErrorValue(fRawRunData->fDataNonMusr.GetErrData()->at(yIndex).at(i));
|
|
} else { // packed data, i.e. fRunInfo->GetPacking() > 1
|
|
if ((i % fRunInfo->GetPacking() == 0) && (i != 0)) { // fill data
|
|
fData.AppendXValue(fRawRunData->fDataNonMusr.GetData()->at(xIndex).at(i)-(fRawRunData->fDataNonMusr.GetData()->at(xIndex).at(i)-fRawRunData->fDataNonMusr.GetData()->at(xIndex).at(i-fRunInfo->GetPacking()))/2.0);
|
|
fData.AppendValue(value);
|
|
fData.AppendErrorValue(TMath::Sqrt(err));
|
|
value = 0.0;
|
|
err = 0.0;
|
|
}
|
|
// sum raw data values
|
|
value += fRawRunData->fDataNonMusr.GetData()->at(yIndex).at(i);
|
|
err += fRawRunData->fDataNonMusr.GetErrData()->at(yIndex).at(i)*fRawRunData->fDataNonMusr.GetErrData()->at(yIndex).at(i);
|
|
}
|
|
}
|
|
|
|
// count the number of bins to be fitted
|
|
fNoOfFitBins = fData.GetValue()->size();
|
|
|
|
// fill theory histo
|
|
// feed the parameter vector
|
|
std::vector<Double_t> par;
|
|
PMsrParamList *paramList = fMsrInfo->GetMsrParamList();
|
|
for (UInt_t i=0; i<paramList->size(); i++)
|
|
par.push_back((*paramList)[i].fValue);
|
|
// calculate functions
|
|
for (Int_t i=0; i<fMsrInfo->GetNoOfFuncs(); i++) {
|
|
fFuncValues[i] = fMsrInfo->EvalFunc(fMsrInfo->GetFuncNo(i), *fRunInfo->GetMap(), par);
|
|
}
|
|
|
|
// get plot range
|
|
PMsrPlotList *plotList;
|
|
PMsrPlotStructure plotBlock;
|
|
plotList = fMsrInfo->GetMsrPlotList();
|
|
// find the proper plot block
|
|
// Here a small complication to be handled: there are potentially multiple
|
|
// run blocks and the run might be present in various of these run blocks. In
|
|
// order to get a nice resolution on the theory the following procedure will be
|
|
// followed: the smallest x-interval found will be used to for the fXTheory resolution
|
|
// which is 1000 function points. The function will be calculated from the smallest
|
|
// xmin found up to the largest xmax found.
|
|
Double_t xMin = 0.0, xMax = 0.0;
|
|
|
|
// init xMin/xMax, xAbsMin/xAbsMax
|
|
plotBlock = plotList->at(0);
|
|
if (plotBlock.fTmin.size() == 0) { // check if no range information is present
|
|
PMsrRunList *runList = fMsrInfo->GetMsrRunList();
|
|
xMin = runList->at(0).GetFitRange(0);
|
|
xMax = runList->at(0).GetFitRange(1);
|
|
for (UInt_t i=1; i<runList->size(); i++) {
|
|
if (runList->at(i).GetFitRange(0) < xMin)
|
|
xMin = runList->at(i).GetFitRange(0);
|
|
if (runList->at(i).GetFitRange(1) > xMax)
|
|
xMax = runList->at(i).GetFitRange(1);
|
|
}
|
|
} else if (plotBlock.fTmin.size() == 1) { // check if 'range' information is present
|
|
xMin = plotBlock.fTmin[0];
|
|
xMax = plotBlock.fTmax[0];
|
|
} else if (plotBlock.fTmin.size() > 1) { // check if 'sub_ranges' information is present
|
|
xMin = plotBlock.fTmin[0];
|
|
xMax = plotBlock.fTmax[0];
|
|
for (UInt_t i=1; i<plotBlock.fTmin.size(); i++) {
|
|
if (plotBlock.fTmin[i] < xMin)
|
|
xMin = plotBlock.fTmin[i];
|
|
if (plotBlock.fTmax[i] > xMax)
|
|
xMax = plotBlock.fTmax[i];
|
|
}
|
|
}
|
|
|
|
if (plotBlock.fUseFitRanges) { // check if 'use_fit_ranges' information is present
|
|
PMsrRunList *runList = fMsrInfo->GetMsrRunList();
|
|
xMin = runList->at(0).GetFitRange(0);
|
|
xMax = runList->at(0).GetFitRange(1);
|
|
for (UInt_t i=1; i<runList->size(); i++) {
|
|
if (runList->at(i).GetFitRange(0) < xMin)
|
|
xMin = runList->at(i).GetFitRange(0);
|
|
if (runList->at(i).GetFitRange(1) > xMax)
|
|
xMax = runList->at(i).GetFitRange(1);
|
|
}
|
|
}
|
|
|
|
for (UInt_t i=1; i<plotList->size(); i++) { // go through all the plot blocks
|
|
plotBlock = plotList->at(i);
|
|
|
|
if (plotBlock.fTmin.size() == 0) { // check if no range information is present
|
|
PMsrRunList *runList = fMsrInfo->GetMsrRunList();
|
|
for (UInt_t i=0; i<runList->size(); i++) {
|
|
if (runList->at(i).GetFitRange(0) < xMin)
|
|
xMin = runList->at(i).GetFitRange(0);
|
|
if (runList->at(i).GetFitRange(1) > xMax)
|
|
xMax = runList->at(i).GetFitRange(1);
|
|
}
|
|
} else if (plotBlock.fTmin.size() == 1) { // check if 'range' information is present
|
|
if (plotBlock.fTmin[0] < xMin)
|
|
xMin = plotBlock.fTmin[0];
|
|
if (plotBlock.fTmax[0] > xMax)
|
|
xMax = plotBlock.fTmax[0];
|
|
} else if (plotBlock.fTmin.size() > 1) { // check if 'sub_ranges' information is present
|
|
for (UInt_t i=0; i<plotBlock.fTmin.size(); i++) {
|
|
if (plotBlock.fTmin[i] < xMin)
|
|
xMin = plotBlock.fTmin[i];
|
|
if (plotBlock.fTmax[i] > xMax)
|
|
xMax = plotBlock.fTmax[i];
|
|
}
|
|
}
|
|
|
|
if (plotBlock.fUseFitRanges) { // check if 'use_fit_ranges' information is present
|
|
PMsrRunList *runList = fMsrInfo->GetMsrRunList();
|
|
for (UInt_t i=0; i<runList->size(); i++) {
|
|
if (runList->at(i).GetFitRange(0) < xMin)
|
|
xMin = runList->at(i).GetFitRange(0);
|
|
if (runList->at(i).GetFitRange(1) > xMax)
|
|
xMax = runList->at(i).GetFitRange(1);
|
|
}
|
|
}
|
|
}
|
|
|
|
// typically take 1000 points to calculate the theory, except if there are more data points, than take that number
|
|
Double_t xStep;
|
|
if (fData.GetX()->size() > 1000.0)
|
|
xStep = (xMax-xMin)/fData.GetX()->size();
|
|
else
|
|
xStep = (xMax-xMin)/1000.0;
|
|
|
|
Double_t xx = xMin;
|
|
do {
|
|
// fill x-vector
|
|
fData.AppendXTheoryValue(xx);
|
|
// fill y-vector
|
|
fData.AppendTheoryValue(fTheory->Func(xx, par, fFuncValues));
|
|
// calculate next xx
|
|
xx += xStep;
|
|
} while (xx < xMax);
|
|
|
|
// clean up
|
|
par.clear();
|
|
|
|
return success;
|
|
}
|
|
|
|
//--------------------------------------------------------------------------
|
|
// GetXIndex
|
|
//--------------------------------------------------------------------------
|
|
/**
|
|
* <p>Returns the x-axis data index.
|
|
*
|
|
* <b>return:</b>
|
|
* - x-index
|
|
*/
|
|
UInt_t PRunNonMusr::GetXIndex()
|
|
{
|
|
UInt_t index = 0;
|
|
Bool_t found = false;
|
|
|
|
if (fRawRunData->fDataNonMusr.FromAscii()) { // ascii-file format
|
|
index = 0;
|
|
found = true;
|
|
} else { // db-file format
|
|
if (fRunInfo->GetXDataIndex() > 0) { // xy-data already indices
|
|
index = fRunInfo->GetXDataIndex()-1; // since xy-data start with 1 ...
|
|
found = true;
|
|
} else { // xy-data data tags which needs to be converted to an index
|
|
for (UInt_t i=0; i<fRawRunData->fDataNonMusr.GetDataTags()->size(); i++) {
|
|
if (fRawRunData->fDataNonMusr.GetDataTags()->at(i).CompareTo(*fRunInfo->GetXDataLabel()) == 0) {
|
|
index = i;
|
|
found = true;
|
|
break;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
if (!found) {
|
|
cerr << endl << "PRunNonMusr::GetXIndex(): **ERROR** Couldn't obtain x-data index!";
|
|
cerr << endl;
|
|
assert(0);
|
|
}
|
|
|
|
return index;
|
|
}
|
|
|
|
//--------------------------------------------------------------------------
|
|
// GetYIndex
|
|
//--------------------------------------------------------------------------
|
|
/**
|
|
* <p>Returns the y-axis data index.
|
|
*
|
|
* <b>return:</b>
|
|
* - y-index
|
|
*/
|
|
UInt_t PRunNonMusr::GetYIndex()
|
|
{
|
|
UInt_t index = 0;
|
|
Bool_t found = false;
|
|
|
|
if (fRawRunData->fDataNonMusr.FromAscii()) { // ascii-file format
|
|
index = 1;
|
|
found = true;
|
|
} else { // db-file format
|
|
if (fRunInfo->GetYDataIndex() > 0) { // xy-data already indices
|
|
index = fRunInfo->GetYDataIndex()-1; // since xy-data start with 1 ...
|
|
found = true;
|
|
} else { // xy-data data tags which needs to be converted to an index
|
|
for (UInt_t i=0; i<fRawRunData->fDataNonMusr.GetDataTags()->size(); i++) {
|
|
if (fRawRunData->fDataNonMusr.GetDataTags()->at(i).CompareTo(*fRunInfo->GetYDataLabel()) == 0) {
|
|
index = i;
|
|
found = true;
|
|
break;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
if (!found) {
|
|
cerr << endl << "PRunNonMusr::GetYIndex(): **ERROR** Couldn't obtain y-data index!";
|
|
cerr << endl;
|
|
assert(0);
|
|
}
|
|
|
|
return index;
|
|
}
|