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

made startTimeBin and endTimeBin class variables.

Browse Source
This commit is contained in:
suter_a
2016-04-25 12:18:27 +02:00
parent 71fbce110b
commit 319fb45f97
2 changed files with 30 additions and 41 deletions
Download Patch File Download Diff File Expand all files Collapse all files

64
src/classes/PRunMuMinus.cpp
View File

@ -52,6 +52,9 @@ PRunMuMinus::PRunMuMinus() : PRunBase()
fGoodBins[0] = -1; fGoodBins[0] = -1;
fGoodBins[1] = -1; fGoodBins[1] = -1;
fStartTimeBin = -1;
fEndTimeBin = -1;
fHandleTag = kEmpty; fHandleTag = kEmpty;
} }
@ -87,6 +90,9 @@ PRunMuMinus::PRunMuMinus(PMsrHandler *msrInfo, PRunDataHandler *rawData, UInt_t
fGoodBins[0] = -1; fGoodBins[0] = -1;
fGoodBins[1] = -1; fGoodBins[1] = -1;
fStartTimeBin = -1;
fEndTimeBin = -1;
if (!PrepareData()) { if (!PrepareData()) {
cerr << endl << ">> PRunMuMinus::PRunMuMinus: **SEVERE ERROR**: Couldn't prepare data for fitting!"; cerr << endl << ">> PRunMuMinus::PRunMuMinus: **SEVERE ERROR**: Couldn't prepare data for fitting!";
cerr << endl << ">> This is very bad :-(, will quit ..."; cerr << endl << ">> This is very bad :-(, will quit ...";
@ -130,15 +136,7 @@ Double_t PRunMuMinus::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
@ -147,12 +145,12 @@ Double_t PRunMuMinus::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));
@ -186,15 +184,7 @@ Double_t PRunMuMinus::CalcChiSquareExpected(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
@ -203,12 +193,12 @@ Double_t PRunMuMinus::CalcChiSquareExpected(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();
theo = fTheory->Func(time, par, fFuncValues); theo = fTheory->Func(time, par, fFuncValues);
diff = fData.GetValue()->at(i) - theo; diff = fData.GetValue()->at(i) - theo;
@ -243,15 +233,7 @@ Double_t PRunMuMinus::CalcMaxLikelihood(const std::vector<Double_t>& par)
Double_t theo; Double_t theo;
Double_t data; Double_t data;
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
@ -260,12 +242,12 @@ Double_t PRunMuMinus::CalcMaxLikelihood(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,theo,data) schedule(dynamic,chunk) reduction(-:mllh) #pragma omp parallel for default(shared) private(i,time,theo,data) schedule(dynamic,chunk) reduction(-:mllh)
#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();
// calculate theory for the given parameter set // calculate theory for the given parameter set
theo = fTheory->Func(time, par, fFuncValues); theo = fTheory->Func(time, par, fFuncValues);
@ -401,15 +383,15 @@ void PRunMuMinus::SetFitRangeBin(const TString fitRange)
void PRunMuMinus::CalcNoOfFitBins() void PRunMuMinus::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;
} }

7
src/include/PRunMuMinus.h
View File

@ -51,6 +51,10 @@ class PRunMuMinus : 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; }
virtual Int_t GetPacking() { return fPacking; }
protected: protected:
virtual void CalcNoOfFitBins(); virtual void CalcNoOfFitBins();
virtual Bool_t PrepareData(); virtual Bool_t PrepareData();
@ -65,6 +69,9 @@ class PRunMuMinus : public PRunBase
PDoubleVector fForward; ///< forward histo data PDoubleVector fForward; ///< forward histo data
Int_t fStartTimeBin; ///< bin at which the fit starts
Int_t fEndTimeBin; ///< bin at which the fit ends
virtual Bool_t GetProperT0(PRawRunData* runData, PMsrGlobalBlock *globalBlock, PUIntVector &histoNo); virtual Bool_t GetProperT0(PRawRunData* runData, PMsrGlobalBlock *globalBlock, PUIntVector &histoNo);
virtual Bool_t GetProperDataRange(); virtual Bool_t GetProperDataRange();
virtual void GetProperFitRange(PMsrGlobalBlock *globalBlock); virtual void GetProperFitRange(PMsrGlobalBlock *globalBlock);