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

Merged muonspin/musrfit into master

Browse Source
This commit is contained in:
Zaher Salman 2016-05-03 12:39:29 +02:00
commit c9ddde8ea1
4 changed files with 37 additions and 49 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

30
src/classes/PRunAsymmetryRRF.cpp
View File

@ -187,15 +187,7 @@ Double_t PRunAsymmetryRRF::CalcChiSquare(const std::vector<Double_t>& par)
// calculate chi square // calculate chi square
Double_t time(1.0); Double_t time(1.0);
Int_t i, N(static_cast<Int_t>(fData.GetValue()->size())); Int_t i;
// 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()));
if (startTimeBin < 0)
startTimeBin = 0;
Int_t endTimeBin = static_cast<Int_t>(floor((fFitEndTime - fData.GetDataTimeStart())/fData.GetDataTimeStep())) + 1;
if (endTimeBin > N)
endTimeBin = N;
// Calculate the theory function once to ensure one function evaluation for the current set of parameters. // Calculate the theory function once to ensure one function evaluation for the current set of parameters.
// This is needed for the LF and user functions where some non-thread-save calculations only need to be calculated once // This is needed for the LF and user functions where some non-thread-save calculations only need to be calculated once
@ -204,12 +196,12 @@ Double_t PRunAsymmetryRRF::CalcChiSquare(const std::vector<Double_t>& par)
asymFcnValue = fTheory->Func(time, par, fFuncValues); asymFcnValue = fTheory->Func(time, par, fFuncValues);
#ifdef HAVE_GOMP #ifdef HAVE_GOMP
Int_t chunk = (endTimeBin - startTimeBin)/omp_get_num_procs(); Int_t chunk = (fEndTimeBin - fStartTimeBin)/omp_get_num_procs();
if (chunk < 10) if (chunk < 10)
chunk = 10; chunk = 10;
#pragma omp parallel for default(shared) private(i,time,diff,asymFcnValue,a,b,f) schedule(dynamic,chunk) reduction(+:chisq) #pragma omp parallel for default(shared) private(i,time,diff,asymFcnValue,a,b,f) schedule(dynamic,chunk) reduction(+:chisq)
#endif #endif
for (i=startTimeBin; i<endTimeBin; ++i) { for (i=fStartTimeBin; i<fEndTimeBin; ++i) {
time = fData.GetDataTimeStart() + (Double_t)i*fData.GetDataTimeStep(); time = fData.GetDataTimeStart() + (Double_t)i*fData.GetDataTimeStep();
switch (fAlphaBetaTag) { switch (fAlphaBetaTag) {
case 1: // alpha == 1, beta == 1 case 1: // alpha == 1, beta == 1
@ -387,15 +379,15 @@ void PRunAsymmetryRRF::SetFitRangeBin(const TString fitRange)
void PRunAsymmetryRRF::CalcNoOfFitBins() void PRunAsymmetryRRF::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 // 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<Int_t>(ceil((fFitStartTime - fData.GetDataTimeStart())/fData.GetDataTimeStep())); fStartTimeBin = static_cast<Int_t>(ceil((fFitStartTime - fData.GetDataTimeStart())/fData.GetDataTimeStep()));
if (startTimeBin < 0) if (fStartTimeBin < 0)
startTimeBin = 0; fStartTimeBin = 0;
Int_t endTimeBin = static_cast<Int_t>(floor((fFitEndTime - fData.GetDataTimeStart())/fData.GetDataTimeStep())) + 1; fEndTimeBin = static_cast<Int_t>(floor((fFitEndTime - fData.GetDataTimeStart())/fData.GetDataTimeStep())) + 1;
if (endTimeBin > static_cast<Int_t>(fData.GetValue()->size())) if (fEndTimeBin > static_cast<Int_t>(fData.GetValue()->size()))
endTimeBin = fData.GetValue()->size(); fEndTimeBin = fData.GetValue()->size();
if (endTimeBin > startTimeBin) if (fEndTimeBin > fStartTimeBin)
fNoOfFitBins = endTimeBin - startTimeBin; fNoOfFitBins = fEndTimeBin - fStartTimeBin;
else else
fNoOfFitBins = 0; fNoOfFitBins = 0;
} }

44
src/classes/PRunSingleHistoRRF.cpp
View File

@ -159,15 +159,7 @@ Double_t PRunSingleHistoRRF::CalcChiSquare(const std::vector<Double_t>& par)
// calculate chi square // calculate chi square
Double_t time(1.0); Double_t time(1.0);
Int_t i, N(static_cast<Int_t>(fData.GetValue()->size())); Int_t i;
// 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()));
if (startTimeBin < 0)
startTimeBin = 0;
Int_t endTimeBin = static_cast<Int_t>(floor((fFitEndTime - fData.GetDataTimeStart())/fData.GetDataTimeStep())) + 1;
if (endTimeBin > N)
endTimeBin = N;
// Calculate the theory function once to ensure one function evaluation for the current set of parameters. // Calculate the theory function once to ensure one function evaluation for the current set of parameters.
// This is needed for the LF and user functions where some non-thread-save calculations only need to be calculated once // This is needed for the LF and user functions where some non-thread-save calculations only need to be calculated once
@ -176,12 +168,12 @@ Double_t PRunSingleHistoRRF::CalcChiSquare(const std::vector<Double_t>& par)
time = fTheory->Func(time, par, fFuncValues); time = fTheory->Func(time, par, fFuncValues);
#ifdef HAVE_GOMP #ifdef HAVE_GOMP
Int_t chunk = (endTimeBin - startTimeBin)/omp_get_num_procs(); Int_t chunk = (fEndTimeBin - fStartTimeBin)/omp_get_num_procs();
if (chunk < 10) if (chunk < 10)
chunk = 10; chunk = 10;
#pragma omp parallel for default(shared) private(i,time,diff) schedule(dynamic,chunk) reduction(+:chisq) #pragma omp parallel for default(shared) private(i,time,diff) schedule(dynamic,chunk) reduction(+:chisq)
#endif #endif
for (i=startTimeBin; i<endTimeBin; ++i) { for (i=fStartTimeBin; i<fEndTimeBin; ++i) {
time = fData.GetDataTimeStart() + (Double_t)i*fData.GetDataTimeStep(); time = fData.GetDataTimeStart() + (Double_t)i*fData.GetDataTimeStep();
diff = fData.GetValue()->at(i) - fTheory->Func(time, par, fFuncValues); diff = fData.GetValue()->at(i) - fTheory->Func(time, par, fFuncValues);
chisq += diff*diff / (fData.GetError()->at(i)*fData.GetError()->at(i)); chisq += diff*diff / (fData.GetError()->at(i)*fData.GetError()->at(i));
@ -215,15 +207,7 @@ Double_t PRunSingleHistoRRF::CalcChiSquareExpected(const std::vector<Double_t>&
// calculate chi square // calculate chi square
Double_t time(1.0); Double_t time(1.0);
Int_t i, N(static_cast<Int_t>(fData.GetValue()->size())); Int_t i;
// 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()));
if (startTimeBin < 0)
startTimeBin = 0;
Int_t endTimeBin = static_cast<Int_t>(floor((fFitEndTime - fData.GetDataTimeStart())/fData.GetDataTimeStep())) + 1;
if (endTimeBin > N)
endTimeBin = N;
// Calculate the theory function once to ensure one function evaluation for the current set of parameters. // Calculate the theory function once to ensure one function evaluation for the current set of parameters.
// This is needed for the LF and user functions where some non-thread-save calculations only need to be calculated once // This is needed for the LF and user functions where some non-thread-save calculations only need to be calculated once
@ -232,12 +216,12 @@ Double_t PRunSingleHistoRRF::CalcChiSquareExpected(const std::vector<Double_t>&
time = fTheory->Func(time, par, fFuncValues); time = fTheory->Func(time, par, fFuncValues);
#ifdef HAVE_GOMP #ifdef HAVE_GOMP
Int_t chunk = (endTimeBin - startTimeBin)/omp_get_num_procs(); Int_t chunk = (fEndTimeBin - fStartTimeBin)/omp_get_num_procs();
if (chunk < 10) if (chunk < 10)
chunk = 10; chunk = 10;
#pragma omp parallel for default(shared) private(i,time,diff) schedule(dynamic,chunk) reduction(+:chisq) #pragma omp parallel for default(shared) private(i,time,diff) schedule(dynamic,chunk) reduction(+:chisq)
#endif #endif
for (i=startTimeBin; i < endTimeBin; ++i) { for (i=fStartTimeBin; i < fEndTimeBin; ++i) {
time = fData.GetDataTimeStart() + (Double_t)i*fData.GetDataTimeStep(); time = fData.GetDataTimeStart() + (Double_t)i*fData.GetDataTimeStep();
theo = fTheory->Func(time, par, fFuncValues); theo = fTheory->Func(time, par, fFuncValues);
diff = fData.GetValue()->at(i) - theo; diff = fData.GetValue()->at(i) - theo;
@ -412,15 +396,15 @@ void PRunSingleHistoRRF::SetFitRangeBin(const TString fitRange)
void PRunSingleHistoRRF::CalcNoOfFitBins() void PRunSingleHistoRRF::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 // 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<Int_t>(ceil((fFitStartTime - fData.GetDataTimeStart())/fData.GetDataTimeStep())); fStartTimeBin = static_cast<Int_t>(ceil((fFitStartTime - fData.GetDataTimeStart())/fData.GetDataTimeStep()));
if (startTimeBin < 0) if (fStartTimeBin < 0)
startTimeBin = 0; fStartTimeBin = 0;
Int_t endTimeBin = static_cast<Int_t>(floor((fFitEndTime - fData.GetDataTimeStart())/fData.GetDataTimeStep())) + 1; fEndTimeBin = static_cast<Int_t>(floor((fFitEndTime - fData.GetDataTimeStart())/fData.GetDataTimeStep())) + 1;
if (endTimeBin > static_cast<Int_t>(fData.GetValue()->size())) if (fEndTimeBin > static_cast<Int_t>(fData.GetValue()->size()))
endTimeBin = fData.GetValue()->size(); fEndTimeBin = fData.GetValue()->size();
if (endTimeBin > startTimeBin) if (fEndTimeBin > fStartTimeBin)
fNoOfFitBins = endTimeBin - startTimeBin; fNoOfFitBins = fEndTimeBin - fStartTimeBin;
else else
fNoOfFitBins = 0; fNoOfFitBins = 0;
} }

6
src/include/PRunAsymmetryRRF.h
View File

@ -52,6 +52,9 @@ class PRunAsymmetryRRF : public PRunBase
virtual void SetFitRangeBin(const TString fitRange); virtual void SetFitRangeBin(const TString fitRange);
virtual Int_t GetStartTimeBin() { return fStartTimeBin; }
virtual Int_t GetEndTimeBin() { return fEndTimeBin; }
protected: protected:
virtual void CalcNoOfFitBins(); virtual void CalcNoOfFitBins();
virtual Bool_t PrepareData(); virtual Bool_t PrepareData();
@ -70,6 +73,9 @@ class PRunAsymmetryRRF : public PRunBase
Int_t fGoodBins[4]; ///< keep first/last good bins. 0=fgb, 1=lgb (forward); 2=fgb, 3=lgb (backward) Int_t fGoodBins[4]; ///< keep first/last good bins. 0=fgb, 1=lgb (forward); 2=fgb, 3=lgb (backward)
Int_t fStartTimeBin; ///< bin at which the fit starts
Int_t fEndTimeBin; ///< bin at which the fit ends
Bool_t SubtractFixBkg(); Bool_t SubtractFixBkg();
Bool_t SubtractEstimatedBkg(); Bool_t SubtractEstimatedBkg();

6
src/include/PRunSingleHistoRRF.h
View File

@ -51,6 +51,9 @@ class PRunSingleHistoRRF : public PRunBase
virtual void SetFitRangeBin(const TString fitRange); virtual void SetFitRangeBin(const TString fitRange);
virtual Int_t GetStartTimeBin() { return fStartTimeBin; }
virtual Int_t GetEndTimeBin() { return fEndTimeBin; }
protected: protected:
virtual void CalcNoOfFitBins(); virtual void CalcNoOfFitBins();
virtual Bool_t PrepareData(); virtual Bool_t PrepareData();
@ -66,6 +69,9 @@ class PRunSingleHistoRRF : public PRunBase
Int_t fGoodBins[2]; ///< keep first/last good bins. 0=fgb, 1=lgb Int_t fGoodBins[2]; ///< keep first/last good bins. 0=fgb, 1=lgb
Int_t fStartTimeBin; ///< bin at which the fit starts
Int_t fEndTimeBin; ///< bin at which the fit ends
PDoubleVector fForward; ///< forward histo data PDoubleVector fForward; ///< forward histo data
PDoubleVector fM; ///< vector holding M(t) = [N(t)-N_bkg] exp(+t/tau). Needed to estimate N0. PDoubleVector fM; ///< vector holding M(t) = [N(t)-N_bkg] exp(+t/tau). Needed to estimate N0.
PDoubleVector fMerr; ///< vector holding the error of M(t): M_err = exp(+t/tau) sqrt(N(t)). PDoubleVector fMerr; ///< vector holding the error of M(t): M_err = exp(+t/tau) sqrt(N(t)).