diff --git a/ChangeLog b/ChangeLog index b82d1c67..bd77ee06 100644 --- a/ChangeLog +++ b/ChangeLog @@ -8,6 +8,7 @@ changes since 0.9.0 =================================== NEW the chi^2 calculation in single-histogram and asymmetry fits is parallelized if musrfit is built using a compiler supporting OpenMP (e.g. GCC >= 4.2) + Using --disable-omp this feature can be disabled on the configure level. NEW any2many: force the user to define the exact NeXus ouput format (HDF4,HDF5,XML) FIXED the evaluation of the LF depolarization functions (before numerous unnecessary calculations were performed) FIXED casting problem between uint32 and time_t for some compilers (MUSR-185) diff --git a/configure.ac b/configure.ac index a658e2b9..a37f17fd 100644 --- a/configure.ac +++ b/configure.ac @@ -691,13 +691,20 @@ fi AC_SUBST(FFTW3_LIBS) AC_SUBST(FFTW3_CFLAGS) -# Check for gomp library. If available musrfit is also linked against it (used for parallel chisq calculation and in libFitPofB). -SAVED_CXXFLAGS="$CXXFLAGS" -CXXFLAGS="$CXXFLAGS -fopenmp" -SAVED_LIBSS="$LIBS" -LIBS="$LIBS -fopenmp -lgomp" -AC_SEARCH_LIBS([omp_get_num_procs], [gomp], [AC_DEFINE([HAVE_GOMP], [1], [Define to 1 if gomp is available])], - [CXXFLAGS="$SAVED_CXXFLAGS" LIBS="$SAVED_LIBS"], []) +dnl ----------------------------------------------- +dnl Ask user if OpenMP support should be disabled (used for parallel chisq calculation and in libFitPofB) +dnl ----------------------------------------------- + +AC_ARG_ENABLE([omp], [AC_HELP_STRING([--enable-omp],[build musrfit with OpenMP support [default=yes]])]) + +if test x"$enable_omp" != xno; then + SAVED_CXXFLAGS="$CXXFLAGS" + CXXFLAGS="$CXXFLAGS -fopenmp" + SAVED_LIBSS="$LIBS" + LIBS="$LIBS -fopenmp -lgomp" + AC_SEARCH_LIBS([omp_get_num_procs], [gomp], [AC_DEFINE([HAVE_GOMP], [1], [Define to 1 if gomp is available])], + [CXXFLAGS="$SAVED_CXXFLAGS" LIBS="$SAVED_LIBS"], []) +fi dnl ----------------------------------------------- dnl Set host specific compiler and linker flags diff --git a/src/classes/PRunAsymmetry.cpp b/src/classes/PRunAsymmetry.cpp index dd2986b2..f3ccfc74 100644 --- a/src/classes/PRunAsymmetry.cpp +++ b/src/classes/PRunAsymmetry.cpp @@ -209,6 +209,7 @@ Double_t PRunAsymmetry::CalcChiSquare(const std::vector& par) asymFcnValue = (f*(a*b+1.0)-(a-1.0))/((a+1.0)-f*(a*b-1.0)); break; default: + asymFcnValue = 0.0; break; } diff = fData.GetValue()->at(i) - asymFcnValue; @@ -243,17 +244,40 @@ Double_t PRunAsymmetry::CalcMaxLikelihood(const std::vector& par) */ UInt_t PRunAsymmetry::GetNoOfFitBins() { - Double_t time; - fNoOfFitBins=0; - for (UInt_t i=0; isize(); i++) { - time = fData.GetDataTimeStart() + (Double_t)i * fData.GetDataTimeStep(); - if ((time >= fFitStartTime) && (time <= fFitEndTime)) - fNoOfFitBins++; - } +// Double_t time; +// fNoOfFitBins=0; +// for (UInt_t i=0; isize(); i++) { +// time = fData.GetDataTimeStart() + (Double_t)i * fData.GetDataTimeStep(); +// if ((time >= fFitStartTime) && (time <= fFitEndTime)) +// fNoOfFitBins++; +// } + CalcNoOfFitBins(); return fNoOfFitBins; } +//-------------------------------------------------------------------------- +// CalcNoOfFitBins (private) +//-------------------------------------------------------------------------- +/** + *

Calculate the number of fitted bins for the current fit range. + */ +void PRunAsymmetry::CalcNoOfFitBins() +{ + // In order not having to loop over all bins and to stay consistent with the chisq method, calculate the start and end bins explicitly + Int_t startTimeBin = static_cast(ceil((fFitStartTime - fData.GetDataTimeStart())/fData.GetDataTimeStep())); + if (startTimeBin < 0) + startTimeBin = 0; + Int_t endTimeBin = static_cast(floor((fFitEndTime - fData.GetDataTimeStart())/fData.GetDataTimeStep())) + 1; + if (endTimeBin > static_cast(fData.GetValue()->size())) + endTimeBin = fData.GetValue()->size(); + + if (endTimeBin > startTimeBin) + fNoOfFitBins = endTimeBin - startTimeBin; + else + fNoOfFitBins = 0; +} + //-------------------------------------------------------------------------- // CalcTheory //-------------------------------------------------------------------------- @@ -662,7 +686,7 @@ Bool_t PRunAsymmetry::PrepareData() * */ Bool_t PRunAsymmetry::SubtractFixBkg() -{ +{ Double_t dval; for (UInt_t i=0; iGetBkgFix(1) * fTimeResolution * 1.0e3; // bkg per ns -> bkg per bin; 1.0e3: us -> ns } - + return true; } @@ -970,14 +994,15 @@ Bool_t PRunAsymmetry::PrepareFitData(PRawRunData* runData, UInt_t histoNo[2]) fData.AppendErrorValue(error); } - // count the number of bins to be fitted - Double_t time; - fNoOfFitBins=0; - for (UInt_t i=0; isize(); i++) { - time = fData.GetDataTimeStart() + (Double_t)i * fData.GetDataTimeStep(); - if ((time >= fFitStartTime) && (time <= fFitEndTime)) - fNoOfFitBins++; - } +// // count the number of bins to be fitted +// Double_t time; +// fNoOfFitBins=0; +// for (UInt_t i=0; isize(); i++) { +// time = fData.GetDataTimeStart() + (Double_t)i * fData.GetDataTimeStep(); +// if ((time >= fFitStartTime) && (time <= fFitEndTime)) +// fNoOfFitBins++; +// } + CalcNoOfFitBins(); // clean up fForward.clear(); @@ -1199,14 +1224,15 @@ Bool_t PRunAsymmetry::PrepareViewData(PRawRunData* runData, UInt_t histoNo[2]) fData.AppendErrorValue(error); } - // count the number of bins to be fitted - Double_t time; - fNoOfFitBins=0; - for (UInt_t i=0; isize(); i++) { - time = fData.GetDataTimeStart() + (Double_t)i * fData.GetDataTimeStep(); - if ((time >= fFitStartTime) && (time <= fFitEndTime)) - fNoOfFitBins++; - } +// // count the number of bins to be fitted +// Double_t time; +// fNoOfFitBins=0; +// for (UInt_t i=0; isize(); i++) { +// time = fData.GetDataTimeStart() + (Double_t)i * fData.GetDataTimeStep(); +// if ((time >= fFitStartTime) && (time <= fFitEndTime)) +// fNoOfFitBins++; +// } + CalcNoOfFitBins(); // clean up fForward.clear(); @@ -1221,6 +1247,7 @@ Bool_t PRunAsymmetry::PrepareViewData(PRawRunData* runData, UInt_t histoNo[2]) } // calculate theory + Double_t time; UInt_t size = runData->GetDataBin(histoNo[0])->size(); Double_t factor = 1.0; if (fData.GetValue()->size() * 10 > runData->GetDataBin(histoNo[0])->size()) { diff --git a/src/classes/PRunSingleHisto.cpp b/src/classes/PRunSingleHisto.cpp index f1454f3e..aa269fee 100644 --- a/src/classes/PRunSingleHisto.cpp +++ b/src/classes/PRunSingleHisto.cpp @@ -360,17 +360,41 @@ void PRunSingleHisto::CalcTheory() */ UInt_t PRunSingleHisto::GetNoOfFitBins() { - fNoOfFitBins=0; - Double_t time; - for (UInt_t i=0; isize(); i++) { - time = fData.GetDataTimeStart() + (Double_t)i*fData.GetDataTimeStep(); - if ((time >= fFitStartTime) && (time <= fFitEndTime)) - fNoOfFitBins++; - } +// fNoOfFitBins=0; +// +// Double_t time; +// for (UInt_t i=0; isize(); i++) { +// time = fData.GetDataTimeStart() + (Double_t)i*fData.GetDataTimeStep(); +// if ((time >= fFitStartTime) && (time <= fFitEndTime)) +// fNoOfFitBins++; +// } + CalcNoOfFitBins(); return fNoOfFitBins; } +//-------------------------------------------------------------------------- +// CalcNoOfFitBins (private) +//-------------------------------------------------------------------------- +/** + *

Calculate the number of fitted bins for the current fit range. + */ +void PRunSingleHisto::CalcNoOfFitBins() +{ + // In order not having to loop over all bins and to stay consistent with the chisq method, calculate the start and end bins explicitly + Int_t startTimeBin = static_cast(ceil((fFitStartTime - fData.GetDataTimeStart())/fData.GetDataTimeStep())); + if (startTimeBin < 0) + startTimeBin = 0; + Int_t endTimeBin = static_cast(floor((fFitEndTime - fData.GetDataTimeStart())/fData.GetDataTimeStep())) + 1; + if (endTimeBin > static_cast(fData.GetValue()->size())) + endTimeBin = fData.GetValue()->size(); + + if (endTimeBin > startTimeBin) + fNoOfFitBins = endTimeBin - startTimeBin; + else + fNoOfFitBins = 0; +} + //-------------------------------------------------------------------------- // PrepareData (private) //-------------------------------------------------------------------------- @@ -702,14 +726,15 @@ Bool_t PRunSingleHisto::PrepareFitData(PRawRunData* runData, const UInt_t histoN } } - // count the number of bins to be fitted - fNoOfFitBins=0; - Double_t time; - for (UInt_t i=0; isize(); i++) { - time = fData.GetDataTimeStart() + (Double_t)i*fData.GetDataTimeStep(); - if ((time >= fFitStartTime) && (time <= fFitEndTime)) - fNoOfFitBins++; - } +// // count the number of bins to be fitted +// fNoOfFitBins=0; +// Double_t time; +// for (UInt_t i=0; isize(); i++) { +// time = fData.GetDataTimeStart() + (Double_t)i*fData.GetDataTimeStep(); +// if ((time >= fFitStartTime) && (time <= fFitEndTime)) +// fNoOfFitBins++; +// } + CalcNoOfFitBins(); return true; } @@ -801,15 +826,16 @@ Bool_t PRunSingleHisto::PrepareRawViewData(PRawRunData* runData, const UInt_t hi value += runData->GetDataBin(histoNo)->at(i); } - // count the number of bins - fNoOfFitBins=0; - - Double_t time; - for (UInt_t i=0; isize(); i++) { - time = fData.GetDataTimeStart() + (Double_t)i*fData.GetDataTimeStep(); - if ((time >= fFitStartTime) && (time <= fFitEndTime)) - fNoOfFitBins++; - } +// // count the number of bins +// fNoOfFitBins=0; +// +// Double_t time; +// for (UInt_t i=0; isize(); i++) { +// time = fData.GetDataTimeStart() + (Double_t)i*fData.GetDataTimeStep(); +// if ((time >= fFitStartTime) && (time <= fFitEndTime)) +// fNoOfFitBins++; +// } + CalcNoOfFitBins(); // fill theory vector for kView // feed the parameter vector @@ -874,6 +900,7 @@ Bool_t PRunSingleHisto::PrepareRawViewData(PRawRunData* runData, const UInt_t hi size = fData.GetValue()->size() * 10; factor = (Double_t)runData->GetDataBin(histoNo)->size() / (Double_t)size; } + Double_t time; Double_t theoryValue; fData.SetTheoryTimeStart(fData.GetDataTimeStart()); fData.SetTheoryTimeStep(fTimeResolution*factor); @@ -1096,13 +1123,14 @@ Bool_t PRunSingleHisto::PrepareViewData(PRawRunData* runData, const UInt_t histo } } - // count the number of bins to be fitted - fNoOfFitBins=0; - for (UInt_t i=0; isize(); i++) { - time = fData.GetDataTimeStart() + (Double_t)i*fData.GetDataTimeStep(); - if ((time >= fFitStartTime) && (time <= fFitEndTime)) - fNoOfFitBins++; - } +// // count the number of bins to be fitted +// fNoOfFitBins=0; +// for (UInt_t i=0; isize(); i++) { +// time = fData.GetDataTimeStart() + (Double_t)i*fData.GetDataTimeStep(); +// if ((time >= fFitStartTime) && (time <= fFitEndTime)) +// fNoOfFitBins++; +// } + CalcNoOfFitBins(); // calculate functions for (Int_t i=0; iGetNoOfFuncs(); i++) { diff --git a/src/include/PRunAsymmetry.h b/src/include/PRunAsymmetry.h index 9f77cf69..bb3ce942 100644 --- a/src/include/PRunAsymmetry.h +++ b/src/include/PRunAsymmetry.h @@ -52,6 +52,7 @@ class PRunAsymmetry : public PRunBase virtual UInt_t GetNoOfFitBins(); protected: + virtual void CalcNoOfFitBins(); virtual Bool_t PrepareData(); virtual Bool_t PrepareFitData(PRawRunData* runData, UInt_t histoNo[2]); virtual Bool_t PrepareViewData(PRawRunData* runData, UInt_t histoNo[2]); diff --git a/src/include/PRunSingleHisto.h b/src/include/PRunSingleHisto.h index 196ddb62..db2c0f96 100644 --- a/src/include/PRunSingleHisto.h +++ b/src/include/PRunSingleHisto.h @@ -51,6 +51,7 @@ class PRunSingleHisto : public PRunBase virtual UInt_t GetNoOfFitBins(); protected: + virtual void CalcNoOfFitBins(); virtual Bool_t PrepareData(); virtual Bool_t PrepareFitData(PRawRunData* runData, const UInt_t histoNo); virtual Bool_t PrepareRawViewData(PRawRunData* runData, const UInt_t histoNo);