diff --git a/ChangeLog b/ChangeLog index bde8abee..efb4cc68 100644 --- a/ChangeLog +++ b/ChangeLog @@ -6,6 +6,8 @@ changes since 0.8.0 =================================== +NEW added the command SCALE_N0_BKG TRUE | FALSE to the command-block. This can be used to force a single histogram fit + to use either 1/ns scaling for N0 and background or 1/bins one. NEW any2many: some more work, including the PSI-BIN write routines which are officially not released yet. NEW extended global mode in msr2data NEW any2many: an attempt to write the universial musr-data-file converter. Just started, needs still some work. diff --git a/src/classes/PFitter.cpp b/src/classes/PFitter.cpp index a6909a43..44d8a150 100644 --- a/src/classes/PFitter.cpp +++ b/src/classes/PFitter.cpp @@ -304,6 +304,8 @@ Bool_t PFitter::CheckCommands() continue; } else if (it->fLine.Contains("MAX_LIKELIHOOD", TString::kIgnoreCase)) { continue; + } else if (it->fLine.Contains("SCALE_N0_BKG", TString::kIgnoreCase)) { + continue; } else if (it->fLine.Contains("INTERACTIVE", TString::kIgnoreCase)) { cmd.first = PMN_INTERACTIVE; cmd.second = cmdLineNo; diff --git a/src/classes/PMusrCanvas.cpp b/src/classes/PMusrCanvas.cpp index 8a1c7a2b..4b5915bd 100644 --- a/src/classes/PMusrCanvas.cpp +++ b/src/classes/PMusrCanvas.cpp @@ -115,6 +115,7 @@ ClassImpQ(PMusrCanvas) */ PMusrCanvas::PMusrCanvas() { + fScaleN0AndBkg = true; fValid = false; fDifferenceView = false; fCurrentPlotView = PV_DATA; @@ -344,6 +345,8 @@ void PMusrCanvas::SetMsrHandler(PMsrHandler *msrHandler) { fMsrHandler = msrHandler; + fScaleN0AndBkg = IsScaleN0AndBkg(); + // check if a fourier block is present in the msr-file, and if yes extract the given values if (fMsrHandler->GetMsrFourierList()->fFourierBlockPresent) { fFourier.fFourierBlockPresent = true; @@ -1298,6 +1301,7 @@ void PMusrCanvas::InitFourier() */ void PMusrCanvas::InitMusrCanvas(const Char_t* title, Int_t wtopx, Int_t wtopy, Int_t ww, Int_t wh) { + fScaleN0AndBkg = true; fValid = false; fDifferenceView = false; fCurrentPlotView = PV_DATA; @@ -3027,7 +3031,10 @@ void PMusrCanvas::PlotData(Bool_t unzoom) if (runList->at(0).IsLifetimeCorrected()) { // lifetime correction yAxisTitle = "asymmetry"; } else { // no liftime correction - yAxisTitle = "N(t) per nsec"; + if (fScaleN0AndBkg) + yAxisTitle = "N(t) per nsec"; + else + yAxisTitle = "N(t) per bin"; } break; case MSR_PLOT_ASYM: @@ -5165,3 +5172,46 @@ void PMusrCanvas::SaveDataAscii() cout << endl << ">> Data windows saved in ascii format ..." << endl; } + +//-------------------------------------------------------------------------- +// IsScaleN0AndBkg (private) +//-------------------------------------------------------------------------- +/** + *

Checks if N0/Bkg normalization to 1/ns is whished. The default is yes, since most of the users want to have it that way. + * To overwrite this, one should add the line 'SCALE_N0_BKG FALSE' to the command block of the msr-file. + * + * return: + * - true, if scaling of N0 and Bkg to 1/ns is whished + * - false, otherwise + * + * \param histoNo forward histogram number of the run + */ +Bool_t PMusrCanvas::IsScaleN0AndBkg() +{ + Bool_t willScale = true; + + PMsrLines *cmd = fMsrHandler->GetMsrCommands(); + for (UInt_t i=0; isize(); i++) { + if (cmd->at(i).fLine.Contains("SCALE_N0_BKG", TString::kIgnoreCase)) { + TObjArray *tokens = 0; + TObjString *ostr = 0; + TString str; + tokens = cmd->at(i).fLine.Tokenize(" \t"); + if (tokens->GetEntries() != 2) { + cerr << endl << ">> PRunSingleHisto::IsScaleN0AndBkg(): **WARNING** Found uncorrect 'SCALE_N0_BKG' command, will ignore it."; + cerr << endl << ">> Allowed commands: SCALE_N0_BKG TRUE | FALSE" << endl; + return willScale; + } + ostr = dynamic_cast(tokens->At(1)); + str = ostr->GetString(); + if (!str.CompareTo("FALSE", TString::kIgnoreCase)) { + willScale = false; + } + // clean up + if (tokens) + delete tokens; + } + } + + return willScale; +} diff --git a/src/classes/PRunSingleHisto.cpp b/src/classes/PRunSingleHisto.cpp index 9cd995ee..241ff0a7 100644 --- a/src/classes/PRunSingleHisto.cpp +++ b/src/classes/PRunSingleHisto.cpp @@ -33,6 +33,10 @@ #include using namespace std; +#include +#include +#include + #include "PMusr.h" #include "PRunSingleHisto.h" @@ -44,6 +48,7 @@ using namespace std; */ PRunSingleHisto::PRunSingleHisto() : PRunBase() { + fScaleN0AndBkg = true; fNoOfFitBins = 0; } @@ -60,6 +65,8 @@ PRunSingleHisto::PRunSingleHisto() : PRunBase() */ PRunSingleHisto::PRunSingleHisto(PMsrHandler *msrInfo, PRunDataHandler *rawData, UInt_t runNo, EPMusrHandleTag tag) : PRunBase(msrInfo, rawData, runNo, tag) { + fScaleN0AndBkg = IsScaleN0AndBkg(); + if (!PrepareData()) { cerr << endl << ">> PRunSingleHisto::PRunSingleHisto: **SEVERE ERROR**: Couldn't prepare data for fitting!"; cerr << endl << ">> This is very bad :-(, will quit ..."; @@ -79,7 +86,7 @@ PRunSingleHisto::~PRunSingleHisto() } //-------------------------------------------------------------------------- -// CalcChiSquare +// CalcChiSquare (public) //-------------------------------------------------------------------------- /** *

Calculate chi-square. @@ -94,7 +101,7 @@ Double_t PRunSingleHisto::CalcChiSquare(const std::vector& par) Double_t chisq = 0.0; Double_t diff = 0.0; - Double_t N0; + Double_t N0 = 0.0; // check if norm is a parameter or a function if (fRunInfo->GetNormParamNo() < MSR_PARAM_FUN_OFFSET) { // norm is a parameter @@ -144,11 +151,14 @@ Double_t PRunSingleHisto::CalcChiSquare(const std::vector& par) // the correction factor is need since the data scales like pack*t_res, // whereas the error scales like sqrt(pack*t_res) - return chisq * fRunInfo->GetPacking() * (fTimeResolution * 1.0e3); + if (fScaleN0AndBkg) + chisq *= fRunInfo->GetPacking() * (fTimeResolution * 1.0e3); + + return chisq; } //-------------------------------------------------------------------------- -// CalcMaxLikelihood +// CalcMaxLikelihood (public) //-------------------------------------------------------------------------- /** *

Calculate log maximum-likelihood. @@ -204,7 +214,9 @@ Double_t PRunSingleHisto::CalcMaxLikelihood(const std::vector& par) Double_t data; Double_t time; // norm is needed since there is no simple scaling like in chisq case to get the correct Max.Log.Likelihood value when normlizing N(t) to 1/ns - Double_t normalizer = fRunInfo->GetPacking() * (fTimeResolution * 1.0e3); + Double_t normalizer = 1.0; + if (fScaleN0AndBkg) + normalizer = fRunInfo->GetPacking() * (fTimeResolution * 1.0e3); for (UInt_t i=0; isize(); i++) { time = fData.GetDataTimeStart() + (Double_t)i*fData.GetDataTimeStep(); if ((time>=fFitStartTime) && (time<=fFitEndTime)) { @@ -225,7 +237,7 @@ Double_t PRunSingleHisto::CalcMaxLikelihood(const std::vector& par) } //-------------------------------------------------------------------------- -// CalcTheory +// CalcTheory (public) //-------------------------------------------------------------------------- /** *

Calculate theory for a given set of fit-parameters. @@ -275,9 +287,9 @@ void PRunSingleHisto::CalcTheory() } // calculate theory - UInt_t size = fData.GetValue()->size(); - Double_t start = fData.GetDataTimeStart(); - Double_t resolution = fData.GetDataTimeStep(); + UInt_t size = fData.GetValue()->size(); + Double_t start = fData.GetDataTimeStart(); + Double_t resolution = fData.GetDataTimeStep(); Double_t time; for (UInt_t i=0; iPrepare data for fitting or viewing. What is already processed at this stage: @@ -501,7 +513,7 @@ Bool_t PRunSingleHisto::PrepareData() // keep the time resolution in (us) fTimeResolution = runData->GetTimeResolution()/1.0e3; cout.precision(8); - cout << endl << ">> PRunSingleHisto::PrepareData(): time resolution=" << fixed << runData->GetTimeResolution() << "(ns)" << endl; + cout << endl << ">> PRunSingleHisto::PrepareData(): time resolution=" << fixed << runData->GetTimeResolution() << "(ns)"; if (fHandleTag == kFit) success = PrepareFitData(runData, histoNo[0]); @@ -519,7 +531,7 @@ Bool_t PRunSingleHisto::PrepareData() } //-------------------------------------------------------------------------- -// PrepareFitData +// PrepareFitData (private) //-------------------------------------------------------------------------- /** *

Take the pre-processed data (i.e. grouping and addrun are preformed) and form the histogram for fitting. @@ -607,9 +619,11 @@ Bool_t PRunSingleHisto::PrepareFitData(PRawRunData* runData, const UInt_t histoN // everything looks fine, hence fill data set Int_t t0 = fT0s[0]; Double_t value = 0.0; + Double_t normalizer = 1.0; // in order that after rebinning the fit does not need to be redone (important for plots) - // the value is normalize to per 1 nsec - Double_t normalizer = fRunInfo->GetPacking() * (fTimeResolution * 1.0e3); // fTimeResolution us->ns + // the value is normalize to per 1 nsec if scaling is whished + if (fScaleN0AndBkg) + normalizer = fRunInfo->GetPacking() * (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)); @@ -651,7 +665,7 @@ Bool_t PRunSingleHisto::PrepareFitData(PRawRunData* runData, const UInt_t histoN } //-------------------------------------------------------------------------- -// PrepareRawViewData +// PrepareRawViewData (private) //-------------------------------------------------------------------------- /** *

Take the pre-processed data (i.e. grouping and addrun are preformed) and form the histogram for viewing @@ -723,7 +737,8 @@ Bool_t PRunSingleHisto::PrepareRawViewData(PRawRunData* runData, const UInt_t hi if (((i-start) % packing == 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 1 nsec - normalizer = packing * (fTimeResolution * 1e3); // fTimeResolution us->ns + if (fScaleN0AndBkg) + normalizer = packing * (fTimeResolution * 1e3); // fTimeResolution us->ns value /= normalizer; fData.AppendValue(value); if (value == 0.0) @@ -828,7 +843,7 @@ Bool_t PRunSingleHisto::PrepareRawViewData(PRawRunData* runData, const UInt_t hi } //-------------------------------------------------------------------------- -// PrepareViewData +// PrepareViewData (private) //-------------------------------------------------------------------------- /** *

Take the pre-processed data (i.e. grouping and addrun are preformed) and form the histogram for viewing @@ -975,7 +990,8 @@ Bool_t PRunSingleHisto::PrepareViewData(PRawRunData* runData, const UInt_t histo if (((i-start) % packing == 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 1 nsec - normalizer = packing * (fTimeResolution * 1.0e3); // fTimeResolution us->ns + if (fScaleN0AndBkg) + normalizer = packing * (fTimeResolution * 1.0e3); // fTimeResolution us->ns value /= normalizer; time = (((Double_t)i-(Double_t)(packing-1)/2.0)-t0)*fTimeResolution; expval = TMath::Exp(+time/tau)/N0; @@ -1106,7 +1122,7 @@ Bool_t PRunSingleHisto::PrepareViewData(PRawRunData* runData, const UInt_t histo } //-------------------------------------------------------------------------- -// EstimatBkg +// EstimatBkg (private) //-------------------------------------------------------------------------- /** *

Estimate the background for a given interval. @@ -1179,11 +1195,58 @@ Bool_t PRunSingleHisto::EstimateBkg(UInt_t histoNo) // forward for (UInt_t i=start; iGetDataBin(histoNo)->at(i); +cout << endl << "debug> bkg=" << bkg << ", end=" << end << ", start=" << start; bkg /= static_cast(end - start + 1); - fBackground = bkg / (fTimeResolution * 1e3); // keep background (per 1 nsec) for chisq, max.log.likelihood, fTimeResolution us->ns + 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) cout << endl << ">> fRunInfo->fRunName=" << fRunInfo->GetRunName()->Data() << ", histNo=" << histoNo << ", fBackground=" << fBackground; return true; } + +//-------------------------------------------------------------------------- +// IsScaleN0AndBkg (private) +//-------------------------------------------------------------------------- +/** + *

Checks if N0/Bkg normalization to 1/ns is whished. The default is yes, since most of the users want to have it that way. + * To overwrite this, one should add the line 'SCALE_N0_BKG FALSE' to the command block of the msr-file. + * + * return: + * - true, if scaling of N0 and Bkg to 1/ns is whished + * - false, otherwise + * + * \param histoNo forward histogram number of the run + */ +Bool_t PRunSingleHisto::IsScaleN0AndBkg() +{ + Bool_t willScale = true; + + PMsrLines *cmd = fMsrInfo->GetMsrCommands(); + for (UInt_t i=0; isize(); i++) { + if (cmd->at(i).fLine.Contains("SCALE_N0_BKG", TString::kIgnoreCase)) { + TObjArray *tokens = 0; + TObjString *ostr = 0; + TString str; + tokens = cmd->at(i).fLine.Tokenize(" \t"); + if (tokens->GetEntries() != 2) { + cerr << endl << ">> PRunSingleHisto::IsScaleN0AndBkg(): **WARNING** Found uncorrect 'SCALE_N0_BKG' command, will ignore it."; + cerr << endl << ">> Allowed commands: SCALE_N0_BKG TRUE | FALSE" << endl; + return willScale; + } + ostr = dynamic_cast(tokens->At(1)); + str = ostr->GetString(); + if (!str.CompareTo("FALSE", TString::kIgnoreCase)) { + willScale = false; + } + // clean up + if (tokens) + delete tokens; + } + } + + return willScale; +} diff --git a/src/include/PMusrCanvas.h b/src/include/PMusrCanvas.h index e32f05ca..114364cd 100644 --- a/src/include/PMusrCanvas.h +++ b/src/include/PMusrCanvas.h @@ -231,6 +231,7 @@ class PMusrCanvas : public TObject, public TQObject virtual void SaveGraphicsAndQuit(Char_t *fileName, Char_t *graphicsFormat); private: + Bool_t fScaleN0AndBkg; ///< true=N0 and background is scaled to (1/ns), otherwise (1/bin) for the single histogram case Bool_t fBatchMode; ///< musrview in ROOT batch mode Bool_t fValid; ///< if true, everything looks OK Bool_t fDifferenceView; ///< tag showing that the shown data, fourier, are the difference between data and theory @@ -318,6 +319,8 @@ class PMusrCanvas : public TObject, public TQObject virtual void SaveDataAscii(); + virtual Bool_t IsScaleN0AndBkg(); + ClassDef(PMusrCanvas, 1) }; diff --git a/src/include/PRunSingleHisto.h b/src/include/PRunSingleHisto.h index 38b01f6b..196ddb62 100644 --- a/src/include/PRunSingleHisto.h +++ b/src/include/PRunSingleHisto.h @@ -57,10 +57,12 @@ class PRunSingleHisto : public PRunBase virtual Bool_t PrepareViewData(PRawRunData* runData, const UInt_t histoNo); private: + 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) - Bool_t EstimateBkg(UInt_t histoNo); + virtual Bool_t EstimateBkg(UInt_t histoNo); + virtual Bool_t IsScaleN0AndBkg(); }; #endif // _PRUNSINGLEHISTO_H_