From a29b790e042c98461cb65b457cfdd726372233f1 Mon Sep 17 00:00:00 2001 From: Andreas Suter Date: Fri, 1 Apr 2016 17:44:49 +0200 Subject: [PATCH] some more work towards DKS/GPU integration --- src/classes/PFitterFcnDKS.cpp | 31 ++++++--- src/classes/PRunListCollection.cpp | 105 ++++++++++++++++++++++------- src/classes/PRunSingleHisto.cpp | 1 + src/include/PRunListCollection.h | 14 +++- src/include/PRunSingleHisto.h | 2 + 5 files changed, 116 insertions(+), 37 deletions(-) diff --git a/src/classes/PFitterFcnDKS.cpp b/src/classes/PFitterFcnDKS.cpp index e68cf0f6..546a9f09 100644 --- a/src/classes/PFitterFcnDKS.cpp +++ b/src/classes/PFitterFcnDKS.cpp @@ -79,17 +79,38 @@ PFitterFcnDKS::~PFitterFcnDKS() */ Double_t PFitterFcnDKS::operator()(const std::vector& par) const { - Double_t value = 0.0; + Double_t value = 0.0, chisq = 0.0; // write parameter to GPU Int_t ierr = fDKS.writeParams(&par[0], par.size()); // loop over all data sets + PSingleHistoParams shp; for (UInt_t i=0; iGetNoOfSingleHisto(); i++) { + // get current values of N0, Nbkg, tau, the functions, and the maps + ierr = fRunListCollection->GetSingleHistoParams(i, par, shp); // set N0, Nbkg + ierr += fDKS.callSetConsts(shp.fN0, shp.fTau, shp.fNbkg); + // set fun values + ierr += fDKS.writeFunctions(&shp.fFun[0], shp.fFun.size()); + // set map values + ierr += fDKS.writeMaps(&shp.fMap[0], shp.fMap.size()); + // calc chisq/log-mlh +/* + chisq = 0.0; + ierr += fDKS.callLauchChiSquare(fMemData[i], fMemData[i], fData.size(), + par.size(), shp.fFun.size(), shp.fMap.size(), + dt0 , dt, 0, chisq); + value += chisq; +*/ + if (ierr != 0) { + cerr << "PFitterFcnDKS::operator(): **ERROR** Kernel launch failed!" << endl; + exit (EXIT_FAILURE); + } + } return value; @@ -230,14 +251,6 @@ void PFitterFcnDKS::InitDKS(const UInt_t dksTag) return; } - // keep the N0 and Nbkg for each data set - fNidx.resize(fRunListCollection->GetNoOfSingleHisto()); - for (UInt_t i=0; iGetNoOfSingleHisto(); i++) { - fNidx[i].fN0 = fRunListCollection->GetN0Idx(i); - fNidx[i].fNbkg = fRunListCollection->GetNbkgIdx(i); - cout << "debug> runNo=" << i << ", N0 idx=" << fNidx[i].fN0 << ", Nbkg idx=" << fNidx[i].fNbkg << endl; - } - // checks device properties if openCL ierr = fDKS.checkMuSRKernels(); if (ierr != 0) { diff --git a/src/classes/PRunListCollection.cpp b/src/classes/PRunListCollection.cpp index 1a902721..4ac1cbdf 100644 --- a/src/classes/PRunListCollection.cpp +++ b/src/classes/PRunListCollection.cpp @@ -1063,43 +1063,96 @@ const Char_t* PRunListCollection::GetYAxisTitle(const TString &runName, const UI } //-------------------------------------------------------------------------- -// GetN0Idx (public) +// GetSingleHistoParams (public) //-------------------------------------------------------------------------- -/** - * \brief PRunListCollection::GetN0Idx - * \param idx - * \return - */ -Int_t PRunListCollection::GetN0Idx(UInt_t idx) +Int_t PRunListCollection::GetSingleHistoParams(UInt_t idx, const std::vector& par, PSingleHistoParams &shp) { - Int_t N0idx = -1; - + Int_t ierr = 0; // make sure idx is within proper bounds if (idx >= fRunSingleHistoList.size()) - return N0idx; + return 1; - N0idx = fMsrInfo->GetMsrRunList()->at(idx).GetNormParamNo(); + // init param + InitSingleHistoParams(shp); - return N0idx; + // check if norm is a parameter or a function + PMsrRunBlock runInfo = fMsrInfo->GetMsrRunList()->at(idx); + if (runInfo.GetNormParamNo() < MSR_PARAM_FUN_OFFSET) { // norm is a parameter + shp.fN0 = par[runInfo.GetNormParamNo()-1]; + } else { // norm is a function + // get function number + UInt_t funNo = runInfo.GetNormParamNo()-MSR_PARAM_FUN_OFFSET; + // evaluate function + shp.fN0 = fMsrInfo->EvalFunc(funNo, *runInfo.GetMap(), par); + } + cout << "debug> shp.fN0 = " << shp.fN0 << endl; + + // get tau + if (runInfo.GetLifetimeParamNo() != -1) + shp.fTau = par[runInfo.GetLifetimeParamNo()-1]; + else + shp.fTau = PMUON_LIFETIME; + cout << "debug> shp.fTau = " << shp.fTau << endl; + + // get background + if (runInfo.GetBkgFitParamNo() == -1) { // bkg not fitted + if (runInfo.GetBkgFix(0) == PMUSR_UNDEFINED) { // no fixed background given (background interval) + shp.fNbkg = GetBackground(idx); + } else { // fixed bkg given + shp.fNbkg = runInfo.GetBkgFix(0); + } + } else { // bkg fitted + shp.fNbkg = par[runInfo.GetBkgFitParamNo()-1]; + } + cout << "debug> shp.fNbkg = " << shp.fNbkg << endl; + + // calculate functions + Int_t funcNo = 0; + for (Int_t i=0; iGetNoOfFuncs(); i++) { + funcNo = fMsrInfo->GetFuncNo(i); + shp.fFun.push_back(fMsrInfo->EvalFunc(funcNo, *runInfo.GetMap(), par)); + } + for (UInt_t i=0; i fun" << i+1 << " = " << shp.fFun[i] << endl; + + // get map vector + shp.fMap = *runInfo.GetMap(); + shp.fMap.erase(shp.fMap.begin()+GetNoOfMaps(), shp.fMap.end()); + for (UInt_t i=0; i map" << i+1 << " = " << shp.fMap[i] << endl; + + return ierr; } //-------------------------------------------------------------------------- -// GetNbkgIdx (public) +// InitSingleHistoParams (private) //-------------------------------------------------------------------------- /** - * \brief PRunListCollection::GetNbkgIdx - * \param idx - * \return + * \brief PRunListCollection::InitSingleHistoParams + * \param param */ -Int_t PRunListCollection::GetNbkgIdx(UInt_t idx) +void PRunListCollection::InitSingleHistoParams(PSingleHistoParams ¶m) { - Int_t NbkgIdx = -1; - - // make sure idx is within proper bounds - if (idx >= fRunSingleHistoList.size()) - return NbkgIdx; - - NbkgIdx = fMsrInfo->GetMsrRunList()->at(idx).GetBkgFitParamNo(); - - return NbkgIdx; + param.fN0 = -1.0; + param.fNbkg = -1.0; + param.fTau = -1.0; + param.fFun.clear(); + param.fMap.clear(); +} + +//-------------------------------------------------------------------------- +// GetBackground (private) +//-------------------------------------------------------------------------- +/** + * @brief PRunListCollection::GetBackground + * @param idx + * @return + */ +Double_t PRunListCollection::GetBackground(Int_t idx) +{ + // make sure idx is within proper bounds + if (idx >= (Int_t)fRunSingleHistoList.size()) + return 0.0; + + return fRunSingleHistoList[idx]->GetBackground(); } diff --git a/src/classes/PRunSingleHisto.cpp b/src/classes/PRunSingleHisto.cpp index 2881a32d..28be62f2 100644 --- a/src/classes/PRunSingleHisto.cpp +++ b/src/classes/PRunSingleHisto.cpp @@ -81,6 +81,7 @@ PRunSingleHisto::PRunSingleHisto(PMsrHandler *msrInfo, PRunDataHandler *rawData, { fScaleN0AndBkg = IsScaleN0AndBkg(); fNoOfFitBins = 0; + fBackground = 0; fPacking = fRunInfo->GetPacking(); if (fPacking == -1) { // i.e. packing is NOT given in the RUN-block, it must be given in the GLOBAL-block diff --git a/src/include/PRunListCollection.h b/src/include/PRunListCollection.h index 649a745f..7f1d878f 100644 --- a/src/include/PRunListCollection.h +++ b/src/include/PRunListCollection.h @@ -43,6 +43,14 @@ using namespace std; #include "PRunMuMinus.h" #include "PRunNonMusr.h" +typedef struct { + Double_t fN0; + Double_t fNbkg; + Double_t fTau; + PDoubleVector fFun; + PIntVector fMap; +} PSingleHistoParams; + /** *

Handler class handling all processed data of an msr-file. All calls of minuit2 are going through this class. */ @@ -103,8 +111,7 @@ class PRunListCollection virtual Int_t GetNoOfParameters() { return fMsrInfo->GetNoOfParams(); } virtual Int_t GetNoOfFunctions() { return fMsrInfo->GetNoOfFuncs(); } virtual Int_t GetNoOfMaps() { return fMsrInfo->GetNoOfMaps(); } - virtual Int_t GetN0Idx(UInt_t idx); - virtual Int_t GetNbkgIdx(UInt_t idx); + virtual Int_t GetSingleHistoParams(UInt_t idx, const std::vector& par, PSingleHistoParams &shp); private: PMsrHandler *fMsrInfo; ///< pointer to the msr-file handler @@ -116,6 +123,9 @@ class PRunListCollection vector fRunAsymmetryRRFList; ///< stores all processed asymmetry RRF data vector fRunMuMinusList; ///< stores all processed mu-minus data vector fRunNonMusrList; ///< stores all processed non-muSR data + + virtual void InitSingleHistoParams(PSingleHistoParams ¶m); + virtual Double_t GetBackground(Int_t idx); }; #endif // _PRUNLISTCOLLECTION_H_ diff --git a/src/include/PRunSingleHisto.h b/src/include/PRunSingleHisto.h index 3980cff1..2e736e82 100644 --- a/src/include/PRunSingleHisto.h +++ b/src/include/PRunSingleHisto.h @@ -51,6 +51,8 @@ class PRunSingleHisto : public PRunBase virtual void SetFitRangeBin(const TString fitRange); + virtual Double_t GetBackground() { return fBackground; } + protected: virtual void CalcNoOfFitBins(); virtual Bool_t PrepareData();