From da2f0a1049c01fe75514e327bcde39e490f00e1f Mon Sep 17 00:00:00 2001 From: Andreas Suter Date: Mon, 25 Apr 2016 12:20:24 +0200 Subject: [PATCH] (i) added mu minus. (ii) single histo is now fed with a proper error vector which covers all allowed background options. --- src/classes/PFitterFcnDKS.cpp | 93 ++++++++++++++++++++++++++++++++++- src/include/PFitterFcnDKS.h | 9 ++-- 2 files changed, 97 insertions(+), 5 deletions(-) diff --git a/src/classes/PFitterFcnDKS.cpp b/src/classes/PFitterFcnDKS.cpp index e028f340..d105d2c3 100644 --- a/src/classes/PFitterFcnDKS.cpp +++ b/src/classes/PFitterFcnDKS.cpp @@ -103,7 +103,7 @@ Double_t PFitterFcnDKS::operator()(const std::vector& 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& par) const } } + // mu mius + for (UInt_t i=0; iGetNoOfMuMinus(); 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 &par, Doub // single histo for (UInt_t i=0; iGetNoOfSingleHisto(); 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; iGetNoOfMuMinus(); 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; iGetNoOfMuMinus(); 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; iGetNoOfSingleHisto(); i++) { runData = fRunListCollection->GetSingleHisto(i); if (runData == 0) { @@ -269,6 +306,7 @@ void PFitterFcnDKS::InitDKS(const UInt_t dksTag) return; } fMemDataSingleHisto[i] = fDKS.allocateMemory(minSize, ierr); + fMemDataSingleHistoErr[i] = fDKS.allocateMemory(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(fMemDataSingleHisto[i], &runData->GetValue()->at(startTimeBin), minSize); + fDKS.writeData(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; iGetNoOfMuMinus(); 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(minSize, ierr); + fMemDataMuMinusErr[i] = fDKS.allocateMemory(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(fMemDataMuMinus[i], &runData->GetValue()->at(startTimeBin), minSize); + fDKS.writeData(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; iGetSingleHisto(i); fDKS.freeMemory(fMemDataSingleHisto[i], runData->GetValue()->size()); + fDKS.freeMemory(fMemDataSingleHistoErr[i], runData->GetValue()->size()); } fMemDataSingleHisto.clear(); + fMemDataSingleHistoErr.clear(); // asymmetry for (UInt_t i=0; i(fMemDataAsymmetry[i], runData->GetValue()->size()); fDKS.freeMemory(fMemDataAsymmetryErr[i], runData->GetValue()->size()); } + fMemDataAsymmetry.clear(); + fMemDataAsymmetryErr.clear(); + + // mu minus + for (UInt_t i=0; iGetMuMinus(i); + fDKS.freeMemory(fMemDataMuMinus[i], runData->GetValue()->size()); + fDKS.freeMemory(fMemDataMuMinusErr[i], runData->GetValue()->size()); + } + fMemDataMuMinus.clear(); + fMemDataMuMinusErr.clear(); } diff --git a/src/include/PFitterFcnDKS.h b/src/include/PFitterFcnDKS.h index e5508e1c..b2325b97 100644 --- a/src/include/PFitterFcnDKS.h +++ b/src/include/PFitterFcnDKS.h @@ -68,9 +68,12 @@ class PFitterFcnDKS : public ROOT::Minuit2::FCNBase mutable DKSBaseMuSR fDKS; - vector fMemDataSingleHisto; ///< vector holding the initial addresses of the single histo data sets on the GPU - vector fMemDataAsymmetry; ///< vector holding the initial addresses of the asymmetry data sets on the GPU - vector fMemDataAsymmetryErr; ///< vector holding the initial addresses of the asymmetry error sets on the GPU + vector fMemDataSingleHisto; ///< vector holding the initial addresses of the single histo data sets on the GPU + vector fMemDataSingleHistoErr; ///< vector holding the initial addresses of the single histo error sets on the GPU + vector fMemDataAsymmetry; ///< vector holding the initial addresses of the asymmetry data sets on the GPU + vector fMemDataAsymmetryErr; ///< vector holding the initial addresses of the asymmetry error sets on the GPU + vector fMemDataMuMinus; ///< vector holding the initial addresses of the mu minus data sets on the GPU + vector fMemDataMuMinusErr; ///< vector holding the initial addresses of the mu minus error sets on the GPU vector fNidx; ///< N0 / Nbkg parameter index vector virtual void InitDKS(const UInt_t dksTag);