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

more work on the global section

Browse Source
This commit is contained in:
suter_a 2014-12-09 15:47:13 +01:00
parent 3f021c112c
commit 918f1eb110
6 changed files with 188 additions and 37 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

66
src/classes/PMsrHandler.cpp
View File

@ -604,7 +604,37 @@ Int_t PMsrHandler::WriteMsrLogFile(const Bool_t messages)
}
break;
case MSR_TAG_GLOBAL:
// STILL QUITE A FEW ITEMS MISSING !!!!!!
sstr = str;
if (sstr.BeginsWith("fittype")) {
fout.width(16);
switch (fGlobal.GetFitType()) {
case MSR_FITTYPE_SINGLE_HISTO:
fout << left << "fittype" << MSR_FITTYPE_SINGLE_HISTO << " (single histogram fit)" << endl;
break;
case MSR_FITTYPE_ASYM:
fout << left << "fittype" << MSR_FITTYPE_ASYM << " (asymmetry fit)" << endl ;
break;
case MSR_FITTYPE_MU_MINUS:
fout << left << "fittype" << MSR_FITTYPE_MU_MINUS << " (mu minus fit)" << endl ;
break;
case MSR_FITTYPE_NON_MUSR:
fout << left << "fittype" << MSR_FITTYPE_NON_MUSR << " (non muSR fit)" << endl ;
break;
default:
break;
}
} else if (sstr.BeginsWith("data")) {
// STILL MISSING
} else if (sstr.BeginsWith("t0")) {
// STILL MISSING
} else if (sstr.BeginsWith("fit")) {
// STILL MISSING
} else if (sstr.BeginsWith("packing")) {
// STILL MISSING
} else {
fout << str.Data() << endl;
}
break;
case MSR_TAG_RUN:
sstr = str;
@ -2645,6 +2675,41 @@ Bool_t PMsrHandler::HandleGlobalEntry(PMsrLines &lines)
} else if (iter->fLine.BeginsWith("fit", TString::kIgnoreCase)) { // fit range
if (tokens->GetEntries() < 3) {
error = true;
} else { // fit given in time, i.e. fit <start> <end>, where <start>, <end> are given as doubles
if (iter->fLine.Contains("fgb", TString::kIgnoreCase)) { // fit given in bins, i.e. fit fgb+n0 lgb-n1
// check 1st entry, i.e. fgb[+n0]
ostr = dynamic_cast<TObjString*>(tokens->At(1));
str = ostr->GetString();
Ssiz_t idx = str.First("+");
TString numStr = str;
if (idx > -1) { // '+' present hence extract n0
numStr.Remove(0,idx+1);
if (numStr.IsFloat()) {
global.SetFitRangeOffset(numStr.Atoi(), 0);
} else {
error = true;
}
} else { // n0 == 0
global.SetFitRangeOffset(0, 0);
}
// check 2nd entry, i.e. lgb[-n1]
ostr = dynamic_cast<TObjString*>(tokens->At(2));
str = ostr->GetString();
idx = str.First("-");
numStr = str;
if (idx > -1) { // '-' present hence extract n1
numStr.Remove(0,idx+1);
if (numStr.IsFloat()) {
global.SetFitRangeOffset(numStr.Atoi(), 1);
} else {
error = true;
}
} else { // n0 == 0
global.SetFitRangeOffset(0, 0);
}
if (!error)
global.SetFitRangeInBins(true);
} else { // fit given in time, i.e. fit <start> <end>, where <start>, <end> are given as doubles
for (Int_t i=1; i<3; i++) {
ostr = dynamic_cast<TObjString*>(tokens->At(i));
@ -2655,6 +2720,7 @@ Bool_t PMsrHandler::HandleGlobalEntry(PMsrLines &lines)
error = true;
}
}
}
} else if (iter->fLine.BeginsWith("packing", TString::kIgnoreCase)) { // packing
if (tokens->GetEntries() < 2) {
error = true;

48
src/classes/PMusr.cpp
View File

@ -709,8 +709,11 @@ PMsrGlobalBlock::PMsrGlobalBlock()
for (UInt_t i=0; i<4; i++) {
fDataRange[i] = -1; // undefined data bin range
}
fFitRangeInBins = false; // default is that fit range is given in time NOT bins
fFitRange[0] = PMUSR_UNDEFINED; // undefined start fit range
fFitRange[1] = PMUSR_UNDEFINED; // undefined end fit range
fFitRangeOffset[0] = -1; // undefined start fit range offset
fFitRangeOffset[1] = -1; // undefined end fit range offset
fPacking = -1; // undefined packing/rebinning
}
@ -833,6 +836,43 @@ void PMsrGlobalBlock::SetFitRange(Double_t dval, UInt_t idx)
fFitRange[idx] = dval;
}
//--------------------------------------------------------------------------
// GetFitRangeOffset (public)
//--------------------------------------------------------------------------
/**
* <p> get fit range offset value at position idx. idx: 0=fit range offset start, 1=fit range offset end.
*
* <b>return:</b>
* - fit range offset value, if idx is within proper boundaries
* - -1, otherwise
*
* \param idx index of the fit range value to be returned
*/
Int_t PMsrGlobalBlock::GetFitRangeOffset(UInt_t idx)
{
if (idx >= 2)
return -1;
return fFitRangeOffset[idx];
}
//--------------------------------------------------------------------------
// SetFitRangeOffset (public)
//--------------------------------------------------------------------------
/**
* <p> set fit range offset value at position idx. Illegale values will be ignored.
*
* \param ival value to be set
* \param idx index of the fit range value to be set
*/
void PMsrGlobalBlock::SetFitRangeOffset(Int_t ival, UInt_t idx)
{
if (idx >= 2)
return;
fFitRangeOffset[idx] = ival;
}
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// implementation PMsrRunBlock
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
@ -863,8 +903,8 @@ PMsrRunBlock::PMsrRunBlock()
fFitRangeInBins = false; // default is that fit range is given in time NOT bins
fFitRange[0] = PMUSR_UNDEFINED; // undefined start fit range
fFitRange[1] = PMUSR_UNDEFINED; // undefined end fit range
fFitRangeOffset[0] = 0; // undefined start fit range offset
fFitRangeOffset[1] = 0; // undefined end fit range offset
fFitRangeOffset[0] = -1; // undefined start fit range offset
fFitRangeOffset[1] = -1; // undefined end fit range offset
fPacking = -1; // undefined packing
fXYDataIndex[0] = -1; // undefined x data index (NonMusr)
fXYDataIndex[1] = -1; // undefined y data index (NonMusr)
@ -916,8 +956,8 @@ void PMsrRunBlock::CleanUp()
fFitRangeInBins = false; // default is that fit range is given in time NOT bins
fFitRange[0] = PMUSR_UNDEFINED; // undefined start fit range
fFitRange[1] = PMUSR_UNDEFINED; // undefined end fit range
fFitRangeOffset[0] = 0; // undefined start fit range offset
fFitRangeOffset[1] = 0; // undefined end fit range offset
fFitRangeOffset[0] = -1; // undefined start fit range offset
fFitRangeOffset[1] = -1; // undefined end fit range offset
fPacking = -1; // undefined packing
fXYDataIndex[0] = -1; // undefined x data index (NonMusr)
fXYDataIndex[1] = -1; // undefined y data index (NonMusr)

4
src/classes/PRunBase.cpp
View File

@ -54,8 +54,8 @@ PRunBase::PRunBase()
fRawData = 0;
fTimeResolution = -1.0;
fFitStartTime = 0.0;
fFitEndTime = 0.0;
fFitStartTime = PMUSR_UNDEFINED;
fFitEndTime = PMUSR_UNDEFINED;
fValid = true;
fHandleTag = kEmpty;

76
src/classes/PRunSingleHisto.cpp
View File

@ -58,6 +58,7 @@ PRunSingleHisto::PRunSingleHisto() : PRunBase()
fScaleN0AndBkg = true;
fNoOfFitBins = 0;
fBackground = 0;
fPacking = -1;
// the 2 following variables are need in case fit range is given in bins, and since
// the fit range can be changed in the command block, these variables need to be accessible
@ -81,6 +82,18 @@ PRunSingleHisto::PRunSingleHisto(PMsrHandler *msrInfo, PRunDataHandler *rawData,
fScaleN0AndBkg = IsScaleN0AndBkg();
fNoOfFitBins = 0;
fPacking = fRunInfo->GetPacking();
if (fPacking == -1) { // i.e. packing is NOT given in the RUN-block, it must be given in the GLOBAL-block
fPacking = fMsrInfo->GetMsrGlobal()->GetPacking();
}
if (fPacking == -1) { // this should NOT happen, somethin is severely wrong
cerr << endl << ">> PRunSingleHisto::PRunSingleHisto: **SEVERE ERROR**: Couldn't find any packing information!";
cerr << endl << ">> This is very bad :-(, will quit ...";
cerr << endl;
fValid = false;
return;
}
// the 2 following variables are need in case fit range is given in bins, and since
// the fit range can be changed in the command block, these variables need to be accessible
fGoodBins[0] = -1;
@ -192,7 +205,7 @@ Double_t PRunSingleHisto::CalcChiSquare(const std::vector<Double_t>& par)
// the correction factor is need since the data scales like pack*t_res,
// whereas the error scales like sqrt(pack*t_res)
if (fScaleN0AndBkg)
chisq *= fRunInfo->GetPacking() * (fTimeResolution * 1.0e3);
chisq *= fPacking * (fTimeResolution * 1.0e3);
return chisq;
}
@ -285,7 +298,7 @@ Double_t PRunSingleHisto::CalcChiSquareExpected(const std::vector<Double_t>& par
// the correction factor is need since the data scales like pack*t_res,
// whereas the error scales like sqrt(pack*t_res)
if (fScaleN0AndBkg)
chisq *= fRunInfo->GetPacking() * (fTimeResolution * 1.0e3);
chisq *= fPacking * (fTimeResolution * 1.0e3);
return chisq;
}
@ -352,7 +365,7 @@ Double_t PRunSingleHisto::CalcMaxLikelihood(const std::vector<Double_t>& par)
Double_t normalizer = 1.0;
if (fScaleN0AndBkg)
normalizer = fRunInfo->GetPacking() * (fTimeResolution * 1.0e3);
normalizer = fPacking * (fTimeResolution * 1.0e3);
// In order not to have an IF in the next loop, determine the start and end bins for the fit range now
Int_t startTimeBin = static_cast<Int_t>(ceil((fFitStartTime - fData.GetDataTimeStart())/fData.GetDataTimeStep()));
@ -608,6 +621,9 @@ Bool_t PRunSingleHisto::PrepareData()
{
Bool_t success = true;
// keep the Global block info
PMsrGlobalBlock *globalBlock = fMsrInfo->GetMsrGlobal();
// get the proper run
PRawRunData* runData = fRawData->GetRunData(*fRunInfo->GetRunName());
if (!runData) { // couldn't get run
@ -644,14 +660,22 @@ Bool_t PRunSingleHisto::PrepareData()
fT0s[i] = fRunInfo->GetT0Bin(i);
}
// fill in the T0's from the GLOBAL block section (if present)
for (UInt_t i=0; i<globalBlock->GetT0BinSize(); i++) {
if (fT0s[i] == -1) { // i.e. not given in the RUN block section
fT0s[i] = globalBlock->GetT0Bin(i);
}
}
// fill in the T0's from the data file, if not already present in the msr-file
for (UInt_t i=0; i<histoNo.size(); i++) {
if (fT0s[i] == -1.0) // i.e. not present in the msr-file, try the data file
if (fT0s[i] == -1.0) { // i.e. not present in the msr-file, try the data file
if (runData->GetT0Bin(histoNo[i]) > 0.0) {
fT0s[i] = runData->GetT0Bin(histoNo[i]);
fRunInfo->SetT0Bin(fT0s[i], i); // keep value for the msr-file
}
}
}
// fill in the T0's gaps, i.e. in case the T0's are NOT in the msr-file and NOT in the data file
for (UInt_t i=0; i<histoNo.size(); i++) {
@ -827,6 +851,15 @@ Bool_t PRunSingleHisto::PrepareFitData(PRawRunData* runData, const UInt_t histoN
Int_t end;
start = fRunInfo->GetDataRange(0);
end = fRunInfo->GetDataRange(1);
// check if data range has been provided, and if not try to get it from the GLOBAL block section
if (start < 0) {
start = fMsrInfo->GetMsrGlobal()->GetDataRange(0);
}
if (end < 0) {
end = fMsrInfo->GetMsrGlobal()->GetDataRange(1);
}
// check if data range has been provided, and if not try to estimate them
if (start < 0) {
Int_t offset = (Int_t)(10.0e-3/fTimeResolution);
@ -864,7 +897,7 @@ Bool_t PRunSingleHisto::PrepareFitData(PRawRunData* runData, const UInt_t histoN
return false;
}
// keep good bins for potential latter use
// keep good bins for potential later use
fGoodBins[0] = start;
fGoodBins[1] = end;
@ -877,6 +910,11 @@ Bool_t PRunSingleHisto::PrepareFitData(PRawRunData* runData, const UInt_t histoN
fRunInfo->SetFitRange(fFitEndTime, 1);
}
// check if fit range is given in the run block. If not, i.e. it has to be found in the global section!
if (fRunInfo->GetFitRange(0) == PMUSR_UNDEFINED) {
// ANYTHING NEEDED AT THIS POINT??? as35
}
// check how the background shall be handled
if (fRunInfo->GetBkgFitParamNo() == -1) { // bkg shall **NOT** be fitted
// subtract background from histogramms ------------------------------------------
@ -908,13 +946,13 @@ Bool_t PRunSingleHisto::PrepareFitData(PRawRunData* runData, const UInt_t histoN
// in order that after rebinning the fit does not need to be redone (important for plots)
// the value is normalize to per 1 nsec if scaling is whished
if (fScaleN0AndBkg)
normalizer = fRunInfo->GetPacking() * (fTimeResolution * 1.0e3); // fTimeResolution us->ns
normalizer = fPacking * (fTimeResolution * 1.0e3); // fTimeResolution us->ns
// data start at data_start-t0
// time shifted so that packing is included correctly, i.e. t0 == t0 after packing
fData.SetDataTimeStart(fTimeResolution*((Double_t)start-(Double_t)t0+(Double_t)(fRunInfo->GetPacking()-1)/2.0));
fData.SetDataTimeStep(fTimeResolution*fRunInfo->GetPacking());
fData.SetDataTimeStart(fTimeResolution*((Double_t)start-(Double_t)t0+(Double_t)(fPacking-1)/2.0));
fData.SetDataTimeStep(fTimeResolution*fPacking);
for (Int_t i=start; i<end; i++) {
if (fRunInfo->GetPacking() == 1) {
if (fPacking == 1) {
value = fForward[i];
value /= normalizer;
fData.AppendValue(value);
@ -922,8 +960,8 @@ Bool_t PRunSingleHisto::PrepareFitData(PRawRunData* runData, const UInt_t histoN
fData.AppendErrorValue(1.0/normalizer);
else
fData.AppendErrorValue(TMath::Sqrt(value));
} else { // packed data, i.e. fRunInfo->GetPacking() > 1
if (((i-start) % fRunInfo->GetPacking() == 0) && (i != start)) { // fill data
} else { // packed data, i.e. fPacking > 1
if (((i-start) % fPacking == 0) && (i != start)) { // fill data
value /= normalizer;
fData.AppendValue(value);
if (value == 0.0)
@ -964,7 +1002,7 @@ Bool_t PRunSingleHisto::PrepareFitData(PRawRunData* runData, const UInt_t histoN
Bool_t PRunSingleHisto::PrepareRawViewData(PRawRunData* runData, const UInt_t histoNo)
{
// check if view_packing is wished
Int_t packing = fRunInfo->GetPacking();
Int_t packing = fPacking;
if (fMsrInfo->GetMsrPlotList()->at(0).fViewPacking > 0) {
packing = fMsrInfo->GetMsrPlotList()->at(0).fViewPacking;
}
@ -974,7 +1012,7 @@ Bool_t PRunSingleHisto::PrepareRawViewData(PRawRunData* runData, const UInt_t hi
if (fScaleN0AndBkg) {
dataNorm = 1.0/ (packing * (fTimeResolution * 1.0e3)); // fTimeResolution us->ns
} else if (!fScaleN0AndBkg && (fMsrInfo->GetMsrPlotList()->at(0).fViewPacking > 0)) {
theoryNorm = (Double_t)fMsrInfo->GetMsrPlotList()->at(0).fViewPacking/(Double_t)fRunInfo->GetPacking();
theoryNorm = (Double_t)fMsrInfo->GetMsrPlotList()->at(0).fViewPacking/(Double_t)packing;
}
// raw data, since PMusrCanvas is doing ranging etc.
@ -1157,7 +1195,7 @@ Bool_t PRunSingleHisto::PrepareRawViewData(PRawRunData* runData, const UInt_t hi
Bool_t PRunSingleHisto::PrepareViewData(PRawRunData* runData, const UInt_t histoNo)
{
// check if view_packing is wished. This is a global option for all PLOT blocks!
Int_t packing = fRunInfo->GetPacking();
Int_t packing = fPacking;
if (fMsrInfo->GetMsrPlotList()->at(0).fViewPacking > 0) {
packing = fMsrInfo->GetMsrPlotList()->at(0).fViewPacking;
}
@ -1171,7 +1209,7 @@ Bool_t PRunSingleHisto::PrepareViewData(PRawRunData* runData, const UInt_t histo
if (fScaleN0AndBkg) {
dataNorm = 1.0/ (packing * (fTimeResolution * 1.0e3)); // fTimeResolution us->ns
} else if (!fScaleN0AndBkg && (fMsrInfo->GetMsrPlotList()->at(0).fViewPacking > 0)) {
theoryNorm = (Double_t)fMsrInfo->GetMsrPlotList()->at(0).fViewPacking/(Double_t)fRunInfo->GetPacking();
theoryNorm = (Double_t)fMsrInfo->GetMsrPlotList()->at(0).fViewPacking/(Double_t)packing;
}
// transform raw histo data. This is done the following way (for details see the manual):
@ -1461,7 +1499,7 @@ void PRunSingleHisto::EstimateN0()
if (fScaleN0AndBkg) {
N0 /= fTimeResolution*1.0e3;
} else {
N0 *= fRunInfo->GetPacking();
N0 *= fPacking;
}
cout << ">> PRunSingleHisto::EstimateN0: found N0=" << param->at(paramNo-1).fValue << ", will set it to N0=" << N0 << endl;
@ -1507,10 +1545,10 @@ Bool_t PRunSingleHisto::EstimateBkg(UInt_t histoNo)
// calculate proper background range
if (beamPeriod != 0.0) {
Double_t timeBkg = (Double_t)(end-start)*(fTimeResolution*fRunInfo->GetPacking()); // length of the background intervall in time
Double_t timeBkg = (Double_t)(end-start)*(fTimeResolution*fPacking); // length of the background intervall in time
UInt_t fullCycles = (UInt_t)(timeBkg/beamPeriod); // how many proton beam cylces can be placed within the proposed background intervall
// correct the end of the background intervall such that the background is as close as possible to a multiple of the proton cylce
end = start + (UInt_t) ((fullCycles*beamPeriod)/(fTimeResolution*fRunInfo->GetPacking()));
end = start + (UInt_t) ((fullCycles*beamPeriod)/(fTimeResolution*fPacking));
cout << endl << "PRunSingleHisto::EstimatBkg(): Background " << start << ", " << end;
if (end == start)
end = fRunInfo->GetBkgRange(1);
@ -1545,7 +1583,7 @@ Bool_t PRunSingleHisto::EstimateBkg(UInt_t histoNo)
if (fScaleN0AndBkg)
fBackground = bkg / (fTimeResolution * 1e3); // keep background (per 1 nsec) for chisq, max.log.likelihood, fTimeResolution us->ns
else
fBackground = bkg * fRunInfo->GetPacking(); // keep background (per bin)
fBackground = bkg * fPacking; // keep background (per bin)
fRunInfo->SetBkgEstimated(fBackground, 0);

6
src/include/PMusr.h
View File

@ -537,20 +537,26 @@ class PMsrGlobalBlock {
virtual Int_t GetDataRange(UInt_t idx);
virtual UInt_t GetT0BinSize() { return fT0.size(); }
virtual Double_t GetT0Bin(UInt_t idx=0);
virtual Bool_t IsFitRangeInBin() { return fFitRangeInBins; }
virtual Double_t GetFitRange(UInt_t idx);
virtual Int_t GetFitRangeOffset(UInt_t idx);
virtual Int_t GetPacking() { return fPacking; }
virtual void SetFitType(Int_t ival) { fFitType = ival; }
virtual void SetDataRange(Int_t ival, Int_t idx);
virtual void SetT0Bin(Double_t dval, Int_t idx=-1);
virtual void SetFitRangeInBins(Bool_t bval) { fFitRangeInBins = bval; }
virtual void SetFitRange(Double_t dval, UInt_t idx);
virtual void SetFitRangeOffset(Int_t ival, UInt_t idx);
virtual void SetPacking(Int_t ival) { fPacking = ival; }
private:
Int_t fFitType; ///< fit type: 0=single histo fit, 2=asymmetry fit, 4=mu^- single histo fit, 8=non muSR fit
Int_t fDataRange[4]; ///< data bin range (fit type 0, 2, 4)
PDoubleVector fT0; ///< t0 bins (fit type 0, 2, 4). if fit type 0 -> f0, f1, f2, ...; if fit type
Bool_t fFitRangeInBins; ///< flag telling if fit range is given in time or in bins
Double_t fFitRange[2]; ///< fit range in (us)
Int_t fFitRangeOffset[2]; ///< if fit range is given in bins it can have the form fit fgb+n0 lgb-n1. This variable holds the n0 and n1.
Int_t fPacking; ///< packing/rebinning
};

1
src/include/PRunSingleHisto.h
View File

@ -62,6 +62,7 @@ class PRunSingleHisto : public PRunBase
Bool_t fScaleN0AndBkg; ///< true=scale N0 and background to 1/ns, otherwise 1/bin
UInt_t fNoOfFitBins; ///< number of bins to be fitted
Double_t fBackground; ///< needed if background range is given (units: 1/bin)
Int_t fPacking; ///< packing for this particular run. Either given in the RUN- or GLOBAL-block.
Int_t fGoodBins[2]; ///< keep first/last good bins. 0=fgb, 1=lgb