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

(i) added mu minus. (ii) single histo is now fed with a proper error vector which covers all allowed background options.

Browse Source
This commit is contained in:
suter_a 2016-04-25 12:20:24 +02:00
parent 319fb45f97
commit da2f0a1049
2 changed files with 97 additions and 5 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

93
src/classes/PFitterFcnDKS.cpp
View File

@ -103,7 +103,7 @@ Double_t PFitterFcnDKS::operator()(const std::vector<Double_t>& par) const
// calc chisq/log-mlh
chisq = 0.0;
ierr += fDKS.callLaunchChiSquare(FITTYPE_SINGLE_HISTO, fMemDataSingleHisto[i], fMemDataSingleHisto[i], dksp.fNoOfFitBins,
ierr += fDKS.callLaunchChiSquare(FITTYPE_SINGLE_HISTO, fMemDataSingleHisto[i], fMemDataSingleHistoErr[i], dksp.fNoOfFitBins,
par.size(), dksp.fFun.size(), dksp.fMap.size(),
dksp.fStartTime , dksp.fPackedTimeResolution, chisq);
value += chisq;
@ -141,6 +141,30 @@ Double_t PFitterFcnDKS::operator()(const std::vector<Double_t>& par) const
}
}
// mu mius
for (UInt_t i=0; i<fRunListCollection->GetNoOfMuMinus(); i++) {
// get current values of N0, Nbkg, tau, the functions, the maps, the time resolution, the fit start time, etc.
ierr = fRunListCollection->GetMuMinusParams(i, par, dksp);
// set fun values
ierr += fDKS.writeFunctions(&dksp.fFun[0], dksp.fFun.size());
// set map values
ierr += fDKS.writeMaps(&dksp.fMap[0], dksp.fMap.size());
// calc chisq/log-mlh
chisq = 0.0;
ierr += fDKS.callLaunchChiSquare(FITTYPE_MU_MINUS, fMemDataMuMinus[i], fMemDataMuMinusErr[i], dksp.fNoOfFitBins,
par.size(), dksp.fFun.size(), dksp.fMap.size(),
dksp.fStartTime , dksp.fPackedTimeResolution, chisq);
value += chisq;
if (ierr != 0) {
cerr << "PFitterFcnDKS::operator(): **ERROR** Kernel launch for mu minus failed!" << endl;
exit (EXIT_FAILURE);
}
}
if (dksp.fScaleN0AndBkg)
norm = dksp.fPackedTimeResolution*1.0e3;
@ -169,7 +193,14 @@ void PFitterFcnDKS::CalcExpectedChiSquare(const std::vector<Double_t> &par, Doub
// single histo
for (UInt_t i=0; i<fRunListCollection->GetNoOfSingleHisto(); i++) {
value = fRunListCollection->GetSingleHistoChisqExpected(par, i); // calculate the expected chisq for single histo run block 'i'
value = fRunListCollection->GetSingleRunChisqExpected(par, i); // calculate the expected chisq for single histo run block 'i'
expectedChisqPerRun.push_back(value);
totalExpectedChisq += value;
}
// mu minus
for (UInt_t i=0; i<fRunListCollection->GetNoOfMuMinus(); i++) {
value = fRunListCollection->GetSingleRunChisqExpected(par, i); // calculate the expected chisq for mu minus run block 'i'
expectedChisqPerRun.push_back(value);
totalExpectedChisq += value;
}
@ -221,6 +252,11 @@ void PFitterFcnDKS::InitDKS(const UInt_t dksTag)
if (minSize > size)
minSize = size;
}
for (UInt_t i=0; i<fRunListCollection->GetNoOfMuMinus(); i++) {
size = fRunListCollection->GetNoOfBinsFitted(i);
if (minSize > size)
minSize = size;
}
if (minSize == 1e9) {
cerr << ">> PFitterFcnDKS::InitDKS: **ERROR** failed to get data size to be fitted." << endl;
fValid = false;
@ -261,6 +297,7 @@ void PFitterFcnDKS::InitDKS(const UInt_t dksTag)
// single histos
fMemDataSingleHisto.resize(fRunListCollection->GetNoOfSingleHisto());
fMemDataSingleHistoErr.resize(fRunListCollection->GetNoOfSingleHisto());
for (UInt_t i=0; i<fRunListCollection->GetNoOfSingleHisto(); i++) {
runData = fRunListCollection->GetSingleHisto(i);
if (runData == 0) {
@ -269,6 +306,7 @@ void PFitterFcnDKS::InitDKS(const UInt_t dksTag)
return;
}
fMemDataSingleHisto[i] = fDKS.allocateMemory<Double_t>(minSize, ierr);
fMemDataSingleHistoErr[i] = fDKS.allocateMemory<Double_t>(minSize, ierr);
if (ierr != 0) {
cerr << ">> PFitterFcnDKS::InitDKS: **ERROR** failed to allocated single histo data set memory (i=" << i << ") on the GPU" << endl;
fValid = false;
@ -281,6 +319,7 @@ void PFitterFcnDKS::InitDKS(const UInt_t dksTag)
return;
}
fDKS.writeData<double>(fMemDataSingleHisto[i], &runData->GetValue()->at(startTimeBin), minSize);
fDKS.writeData<double>(fMemDataSingleHistoErr[i], &runData->GetError()->at(startTimeBin), minSize);
}
if (fRunListCollection->GetNoOfSingleHisto() > 0)
fitType = FITTYPE_SINGLE_HISTO;
@ -322,6 +361,43 @@ void PFitterFcnDKS::InitDKS(const UInt_t dksTag)
}
}
// mu minus
fMemDataMuMinus.resize(fRunListCollection->GetNoOfMuMinus());
fMemDataMuMinusErr.resize(fRunListCollection->GetNoOfMuMinus());
for (UInt_t i=0; i<fRunListCollection->GetNoOfMuMinus(); i++) {
runData = fRunListCollection->GetMuMinus(i);
if (runData == 0) {
cerr << ">> PFitterFcnDKS::InitDKS: **ERROR** failed to get mu minus data set (i=" << i << ") from fRunListCollection." << endl;
fValid = false;
return;
}
fMemDataMuMinus[i] = fDKS.allocateMemory<Double_t>(minSize, ierr);
fMemDataMuMinusErr[i] = fDKS.allocateMemory<Double_t>(minSize, ierr);
if (ierr != 0) {
cerr << ">> PFitterFcnDKS::InitDKS: **ERROR** failed to allocated mu minus data set memory (i=" << i << ") on the GPU" << endl;
fValid = false;
return;
}
Int_t startTimeBin = fRunListCollection->GetStartTimeBin(MSR_FITTYPE_MU_MINUS, i);
if (startTimeBin < 0) {
cerr << ">> PFitterFcnDKS::InitDKS: **ERROR** startTimeBin undefind (mu minus fit)." << endl;
fValid = false;
return;
}
fDKS.writeData<double>(fMemDataMuMinus[i], &runData->GetValue()->at(startTimeBin), minSize);
fDKS.writeData<double>(fMemDataMuMinusErr[i], &runData->GetError()->at(startTimeBin), minSize);
}
if (fRunListCollection->GetNoOfMuMinus() > 0) {
if (fitType == FITTYPE_UNDEFINED) {
fitType = FITTYPE_MU_MINUS;
} else {
cerr << ">>PFitterFcnDKS::InitDKS: **ERROR** mixed fit types found. This is currently not supported!" << endl;
fValid = false;
return;
}
}
// set the function string and compile the program
ierr = fDKS.callCompileProgram(fTheoStr, !fUseChi2);
if (ierr != 0) {
@ -355,8 +431,10 @@ void PFitterFcnDKS::FreeDKS()
for (UInt_t i=0; i<fMemDataSingleHisto.size(); i++) {
runData = fRunListCollection->GetSingleHisto(i);
fDKS.freeMemory<Double_t>(fMemDataSingleHisto[i], runData->GetValue()->size());
fDKS.freeMemory<Double_t>(fMemDataSingleHistoErr[i], runData->GetValue()->size());
}
fMemDataSingleHisto.clear();
fMemDataSingleHistoErr.clear();
// asymmetry
for (UInt_t i=0; i<fMemDataAsymmetry.size(); i++) {
@ -364,4 +442,15 @@ void PFitterFcnDKS::FreeDKS()
fDKS.freeMemory<Double_t>(fMemDataAsymmetry[i], runData->GetValue()->size());
fDKS.freeMemory<Double_t>(fMemDataAsymmetryErr[i], runData->GetValue()->size());
}
fMemDataAsymmetry.clear();
fMemDataAsymmetryErr.clear();
// mu minus
for (UInt_t i=0; i<fMemDataMuMinus.size(); i++) {
runData = fRunListCollection->GetMuMinus(i);
fDKS.freeMemory<Double_t>(fMemDataMuMinus[i], runData->GetValue()->size());
fDKS.freeMemory<Double_t>(fMemDataMuMinusErr[i], runData->GetValue()->size());
}
fMemDataMuMinus.clear();
fMemDataMuMinusErr.clear();
}

9
src/include/PFitterFcnDKS.h
View File

@ -68,9 +68,12 @@ class PFitterFcnDKS : public ROOT::Minuit2::FCNBase
mutable DKSBaseMuSR fDKS;
vector<void *> fMemDataSingleHisto; ///< vector holding the initial addresses of the single histo data sets on the GPU
vector<void *> fMemDataAsymmetry; ///< vector holding the initial addresses of the asymmetry data sets on the GPU
vector<void *> fMemDataAsymmetryErr; ///< vector holding the initial addresses of the asymmetry error sets on the GPU
vector<void *> fMemDataSingleHisto; ///< vector holding the initial addresses of the single histo data sets on the GPU
vector<void *> fMemDataSingleHistoErr; ///< vector holding the initial addresses of the single histo error sets on the GPU
vector<void *> fMemDataAsymmetry; ///< vector holding the initial addresses of the asymmetry data sets on the GPU
vector<void *> fMemDataAsymmetryErr; ///< vector holding the initial addresses of the asymmetry error sets on the GPU
vector<void *> fMemDataMuMinus; ///< vector holding the initial addresses of the mu minus data sets on the GPU
vector<void *> fMemDataMuMinusErr; ///< vector holding the initial addresses of the mu minus error sets on the GPU
vector<PNidx> fNidx; ///< N0 / Nbkg parameter index vector
virtual void InitDKS(const UInt_t dksTag);