From 20665da9ebfb52588edd94718aa4573ad101d9f5 Mon Sep 17 00:00:00 2001 From: nemu Date: Sun, 6 Feb 2011 14:30:07 +0000 Subject: [PATCH] fix of fix concerning cyclotron background correction for asymmetry fits (MUSR-175). For single histogram fits, the chisq/maxLH is now estimated correctly (see ChangeLog for more details) --- ChangeLog | 3 +++ src/classes/PRunAsymmetry.cpp | 2 +- src/classes/PRunSingleHisto.cpp | 23 +++++++++++++---------- 3 files changed, 17 insertions(+), 11 deletions(-) diff --git a/ChangeLog b/ChangeLog index d74fd7ae..40ef4078 100644 --- a/ChangeLog +++ b/ChangeLog @@ -11,6 +11,9 @@ NEW any2many: an attempt to write the universial musr-data-file converter. Just NEW musrt0: added the possibility to show the t0 saved in the data file 's'. Furthermore added the option --getT0FromPromptPeak, -g with : will, in non-interactive mode estimate the t0's from the prompt peak and write it into the msr-file (MUSR-133). +FIXED for single histogram fits, the chisq given was wrong (not the fit-result though). I missed some needed scaling + when normalizing to 1/ns rather than bins. The same is true for log max. likelihood. I am rather depressed that + nobody so far found this but only R. Scheuermann pointed out there might be a problem. FIXED background range correction, such that it is a multiple of the proton cycle time FIXED linking of BMWlibs on Cygwin FIXED various bugs in msr2data diff --git a/src/classes/PRunAsymmetry.cpp b/src/classes/PRunAsymmetry.cpp index cefdc364..7815bc98 100644 --- a/src/classes/PRunAsymmetry.cpp +++ b/src/classes/PRunAsymmetry.cpp @@ -693,7 +693,7 @@ Bool_t PRunAsymmetry::SubtractEstimatedBkg() // calculate proper background range for (UInt_t i=0; i<2; i++) { 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[i]-start[i])*(fTimeResolution*fRunInfo->GetPacking()); // 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[i] = start[i] + (UInt_t) ((fullCycles*beamPeriod)/(fTimeResolution*fRunInfo->GetPacking())); diff --git a/src/classes/PRunSingleHisto.cpp b/src/classes/PRunSingleHisto.cpp index a0748e86..9cd995ee 100644 --- a/src/classes/PRunSingleHisto.cpp +++ b/src/classes/PRunSingleHisto.cpp @@ -142,7 +142,9 @@ Double_t PRunSingleHisto::CalcChiSquare(const std::vector& par) } } - return chisq; + // 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); } //-------------------------------------------------------------------------- @@ -201,14 +203,17 @@ Double_t PRunSingleHisto::CalcMaxLikelihood(const std::vector& par) Double_t theo; 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); for (UInt_t i=0; isize(); i++) { time = fData.GetDataTimeStart() + (Double_t)i*fData.GetDataTimeStep(); if ((time>=fFitStartTime) && (time<=fFitEndTime)) { // calculate theory for the given parameter set theo = N0*TMath::Exp(-time/tau)*(1+fTheory->Func(time, par, fFuncValues))+bkg; + theo *= normalizer; // check if data value is not too small if (fData.GetValue()->at(i) > 1.0e-9) - data = fData.GetValue()->at(i); + data = normalizer*fData.GetValue()->at(i); else data = 1.0e-9; // add maximum log likelihood contribution of bin i @@ -371,7 +376,7 @@ Bool_t PRunSingleHisto::PrepareData() if (fRunInfo->GetT0Size() <= i) { // t0 for i not present in the msr-file, i.e. #t0's != #forward histos if (static_cast(runData->GetT0Size()) > fRunInfo->GetForwardHistoNo(i)-1) { // t0 for i present in the data file fT0s.push_back(runData->GetT0(fRunInfo->GetForwardHistoNo(i)-1)); - } else { // t0 is neither in the run data nor in the msr-file -> will try estimated ones! + } else { // t0 is neither in the run data nor in the msr-file -> will try estimated one! fT0s.push_back(runData->GetT0Estimated(fRunInfo->GetForwardHistoNo(i)-1)); cerr << endl << ">> PRunSingleHisto::PrepareData(): **WARNING** NO t0's found, neither in the run data nor in the msr-file!"; cerr << endl << ">> run: " << fRunInfo->GetRunName()->Data(); @@ -555,7 +560,7 @@ Bool_t PRunSingleHisto::PrepareFitData(PRawRunData* runData, const UInt_t histoN cerr << endl; } - // check if start, end, and t0 make any sense + // check if start and end make any sense // 1st check if start and end are in proper order if (end < start) { // need to swap them Int_t keep = end; @@ -602,7 +607,9 @@ 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 // 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)); @@ -610,7 +617,6 @@ Bool_t PRunSingleHisto::PrepareFitData(PRawRunData* runData, const UInt_t histoN for (Int_t i=start; iGetPacking() == 1) { value = runData->GetDataBin(histoNo)->at(i); - normalizer = fRunInfo->GetPacking() * (fTimeResolution * 1e3); // fTimeResolution us->ns value /= normalizer; fData.AppendValue(value); if (value == 0.0) @@ -619,15 +625,12 @@ Bool_t PRunSingleHisto::PrepareFitData(PRawRunData* runData, const UInt_t histoN fData.AppendErrorValue(TMath::Sqrt(value)); } else { // packed data, i.e. fRunInfo->GetPacking() > 1 if (((i-start) % fRunInfo->GetPacking() == 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 = fRunInfo->GetPacking() * (fTimeResolution * 1e3); // fTimeResolution us->ns value /= normalizer; fData.AppendValue(value); if (value == 0.0) fData.AppendErrorValue(1.0); else - fData.AppendErrorValue(TMath::Sqrt(value/normalizer)); + fData.AppendErrorValue(TMath::Sqrt(value)); // reset values value = 0.0; }