diff --git a/src/ToDo.txt b/src/ToDo.txt index ca239a29..83bec8e7 100644 --- a/src/ToDo.txt +++ b/src/ToDo.txt @@ -41,11 +41,20 @@ short term: * implement table based theory functions (LF stuff) static GKT LF **DONE** 08-03-12 +* setup startup handler with infos like: data path, symbol list, color list, ... + **DONE** 08-04-04 + +* at the moment that startup handler is looked for in the directory where musrfit/musrview + is executed, this is stupid. Define a place where to look for it. + +* check concept to get bin data for musrview on single histo (fBinData in PRunBase). + If it is ok, implement it for the others. + * do I need to cleanup AddText() objects from TPaveText etc myself? **CHECK** * something is strange with the coordinate system in TPaveText! **CHECK** -* setup startup handler with infos like: data path, symbol list, color list, ... +* check problem with the mathmore implementation of 1F1(m;n;z) --------------------- intermediate term: diff --git a/src/classes/PRunAsymmetry.cpp b/src/classes/PRunAsymmetry.cpp index 21ba1a04..391da7fe 100644 --- a/src/classes/PRunAsymmetry.cpp +++ b/src/classes/PRunAsymmetry.cpp @@ -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& par) } // calculate chisq - for (unsigned int i=0; i=fFitStartTime) && (fData.fTime[i]<=fFitStopTime)) { + for (unsigned int i=0; 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; iFunc(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= fFitStartTime) && (fData.fTime[i] <= fFitStopTime)) + for (unsigned int i=0; i= fFitStartTime) && (fFitData.fTime[i] <= fFitStopTime)) fNoOfFitBins++; } diff --git a/src/classes/PRunBase.cpp b/src/classes/PRunBase.cpp index dc57afab..b52e6662 100644 --- a/src/classes/PRunBase.cpp +++ b/src/classes/PRunBase.cpp @@ -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(); diff --git a/src/classes/PRunNonMusr.cpp b/src/classes/PRunNonMusr.cpp index b360d2b0..f1c81ea8 100644 --- a/src/classes/PRunNonMusr.cpp +++ b/src/classes/PRunNonMusr.cpp @@ -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= fFitStartTime) && (fData.fTime[i] <= fFitStopTime)) + for (unsigned int i=0; i= fFitStartTime) && (fFitData.fTime[i] <= fFitStopTime)) fNoOfFitBins++; } diff --git a/src/classes/PRunRRF.cpp b/src/classes/PRunRRF.cpp index 6b79bc94..db4a6c77 100644 --- a/src/classes/PRunRRF.cpp +++ b/src/classes/PRunRRF.cpp @@ -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= fFitStartTime) && (fData.fTime[i] <= fFitStopTime)) + for (unsigned int i=0; i= fFitStartTime) && (fFitData.fTime[i] <= fFitStopTime)) fNoOfFitBins++; } diff --git a/src/classes/PRunSingleHisto.cpp b/src/classes/PRunSingleHisto.cpp index 1a1ea820..99f90778 100644 --- a/src/classes/PRunSingleHisto.cpp +++ b/src/classes/PRunSingleHisto.cpp @@ -127,11 +127,11 @@ double PRunSingleHisto::CalcChiSquare(const std::vector& par) } // calculate chi square - for (unsigned int i=0; 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=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& 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; ifRunName+"_"+(Long_t)fRunInfo->fForwardHistoNo+"_theo.dat"; // ofstream ft(fln.Data(),ios_base::out); -// for (unsigned int i=0; iFunc(fData.fTime[i], par))+bkg; +// for (unsigned int i=0; iFunc(fFitData.fTime[i], par))+bkg; // } // ft.close(); // counter++; @@ -200,13 +200,13 @@ double PRunSingleHisto::CalcMaxLikelihood(const std::vector& par) // calculate maximum log likelihood double theo; double data; - for (unsigned int i=0; i=fFitStartTime) && (fData.fTime[i]<=fFitStopTime)) { + for (unsigned int i=0; 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; iFunc(fData.fTime[i], par, fFuncValues))+bkg); + for (unsigned int i=0; iFunc(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; ifPacking == 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= fFitStartTime) && (fData.fTime[i] <= fFitStopTime)) + for (unsigned int i=0; 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; ifPacking == 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; } diff --git a/src/include/PMusr.h b/src/include/PMusr.h index 50e9eed6..7151925f 100644 --- a/src/include/PMusr.h +++ b/src/include/PMusr.h @@ -98,7 +98,7 @@ using namespace std; */ typedef vector PIntVector; -//------------------------------------------------------------------------ +//------------------------------------------------------------- /** *

*/ @@ -116,7 +116,7 @@ typedef vector PComplexVector; */ typedef vector PStringVector; -//------------------------------------------------------------------------------------------ +//------------------------------------------------------------- /** *

Predominantly used in PRunBase. */ @@ -127,7 +127,7 @@ typedef struct { vector fTheory; } PRunData; -//------------------------------------------------------------------------ +//------------------------------------------------------------- /** *

*/ @@ -142,7 +142,7 @@ typedef struct { vector fDataBin; ///< vector of all histos of a run } PRawRunData; -//------------------------------------------------------------------------ +//------------------------------------------------------------- /** *

*/ diff --git a/src/include/PRunBase.h b/src/include/PRunBase.h index d192e854..e5b55d6a 100644 --- a/src/include/PRunBase.h +++ b/src/include/PRunBase.h @@ -70,7 +70,7 @@ 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 &fData; } + virtual PRunData* GetData() { return &fFitData; } virtual void CleanUp(); virtual bool IsValid() { return fValid; } @@ -82,16 +82,17 @@ class PRunBase PMsrRunStructure *fRunInfo; ///< run info used to filter out needed infos for the run PRunDataHandler *fRawData; ///< holds the raw run data - vector fParamNo; ///< vector of parameter numbers for the specifc run + PIntVector fParamNo; ///< vector of parameter numbers for the specifc run - PRunData fData; ///< data to be fitted - double fTimeResolution; ///< time resolution - vector fT0s; ///< all t0's of a run! The derived classes will handle it + 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 + double fTimeResolution; ///< time resolution + PDoubleVector fT0s; ///< all t0's of a run! The derived classes will handle it virtual bool PrepareData() = 0; // pure virtual, i.e. needs to be implemented by the deriving class!! - vector fFuncValues; ///< is keeping the values of the functions from the FUNCTIONS block - PTheory *fTheory; ///< theory needed to calculate chi-square + PDoubleVector fFuncValues; ///< is keeping the values of the functions from the FUNCTIONS block + PTheory *fTheory; ///< theory needed to calculate chi-square }; #endif // _PRUNBASE_H_