diff --git a/src/classes/PFitterFcnDKS.cpp b/src/classes/PFitterFcnDKS.cpp index 305be307..e68cf0f6 100644 --- a/src/classes/PFitterFcnDKS.cpp +++ b/src/classes/PFitterFcnDKS.cpp @@ -81,6 +81,17 @@ Double_t PFitterFcnDKS::operator()(const std::vector& par) const { Double_t value = 0.0; + // write parameter to GPU + Int_t ierr = fDKS.writeParams(&par[0], par.size()); + + // loop over all data sets + for (UInt_t i=0; iGetNoOfSingleHisto(); i++) { + // set N0, Nbkg + // set fun values + // set map values + // calc chisq/log-mlh + } + return value; } @@ -122,6 +133,8 @@ void PFitterFcnDKS::CalcExpectedChiSquare(const std::vector &par, Doub */ void PFitterFcnDKS::InitDKS(const UInt_t dksTag) { + Int_t ierr = 0; + // if any device was allocated before, free the device resources FreeDKS(); @@ -142,16 +155,89 @@ void PFitterFcnDKS::InitDKS(const UInt_t dksTag) // init chisq buffer on the GPU + // 1) calculated the maximum size for the data needed. + Int_t maxSize = -1, size = -1; + Int_t parSize = -1, mapSize = -1, funSize = -1; + for (UInt_t i=0; iGetNoOfSingleHisto(); i++) { + size = fRunListCollection->GetSingleHisto(i)->GetValue()->size(); + if (maxSize < size) + maxSize = size; + } + for (UInt_t i=0; iGetNoOfAsymmetry(); i++) { + size = fRunListCollection->GetAsymmetry(i)->GetValue()->size(); + if (maxSize < size) + maxSize = size; + } + if (maxSize == -1) { + cerr << ">> PFitterFcnDKS::InitDKS: **ERROR** failed to get data size to be fitted." << endl; + fValid = false; + return; + } + + // 2) get number of parameters / functions / maps + parSize = fRunListCollection->GetNoOfParameters(); + if (parSize == -1) { + cerr << ">> PFitterFcnDKS::InitDKS: **ERROR** failed to get number of fit parameters." << endl; + fValid = false; + return; + } + funSize = fRunListCollection->GetNoOfFunctions(); + if (funSize == -1) { + cerr << ">> PFitterFcnDKS::InitDKS: **ERROR** failed to get number of functions." << endl; + fValid = false; + return; + } + mapSize = fRunListCollection->GetNoOfMaps(); + if (mapSize == -1) { + cerr << ">> PFitterFcnDKS::InitDKS: **ERROR** failed to get number of maps." << endl; + fValid = false; + return; + } + cout << "debug> parSize=" << parSize << ", funSize=" << funSize << ", mapSize=" << mapSize << endl; + + // now ready to init the chisq buffer on the GPU + ierr = fDKS.initChiSquare(maxSize, parSize, funSize, mapSize); + if (ierr != 0) { + cerr << ">> PFitterFcnDKS::InitDKS: **ERROR** failed to allocate the necessary chisq buffer on the GPU." << endl; + fValid = false; + return; + } + // allocate memory for the data on the GPU/CPU and transfer the data sets + fMemData.resize(fRunListCollection->GetNoOfSingleHisto()); + PRunData *runData=0; + for (UInt_t i=0; iGetNoOfSingleHisto(); i++) { + runData = fRunListCollection->GetSingleHisto(i); + if (runData == 0) { + cerr << ">> PFitterFcnDKS::InitDKS: **ERROR** failed to get data set (i=" << i << ") from fRunListCollection." << endl; + fValid = false; + return; + } + fMemData[i] = fDKS.allocateMemory(runData->GetValue()->size(), ierr); + if (ierr != 0) { + cerr << ">> PFitterFcnDKS::InitDKS: **ERROR** failed to allocated data set memory (i=" << i << ") on the GPU" << endl; + fValid = false; + return; + } + fDKS.writeData(fMemData[i], runData->GetValue(), runData->GetValue()->size()); + } // set the function string and compile the program - Int_t ierr = fDKS.callCompileProgram(fTheoStr, !fUseChi2); + ierr = fDKS.callCompileProgram(fTheoStr, !fUseChi2); if (ierr != 0) { cerr << ">> PFitterFcnDKS::InitDKS: **ERROR** failed to compile theory!" << endl; fValid = false; 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) { @@ -172,5 +258,10 @@ void PFitterFcnDKS::InitDKS(const UInt_t dksTag) */ void PFitterFcnDKS::FreeDKS() { - + PRunData *runData=0; + for (UInt_t i=0; iGetSingleHisto(i); + fDKS.freeMemory(fMemData[i], runData->GetValue()->size()); + } + fMemData.clear(); } diff --git a/src/classes/PMsrHandler.cpp b/src/classes/PMsrHandler.cpp index c33a3295..9793b0a0 100644 --- a/src/classes/PMsrHandler.cpp +++ b/src/classes/PMsrHandler.cpp @@ -69,6 +69,8 @@ PMsrHandler::PMsrHandler(const Char_t *fileName, PStartupOptions *startupOptions fFuncHandler = 0; + fNoOfMaps = -1; + // check if the file name given is a path-file-name, and if yes, split it into path and file name. if (fFileName.Contains("/")) { Int_t idx = -1; @@ -5664,6 +5666,17 @@ Bool_t PMsrHandler::CheckMaps() } } + if (result == true) { + PIntVector *map = fRuns[0].GetMap(); + fNoOfMaps = 0; + for (UInt_t i=0; isize(); i++) + if (map->at(i) != 0) + fNoOfMaps++; + if (fNoOfMaps == 0) + fNoOfMaps = -1; + } + cout << "debug> fNoOfMaps=" << fNoOfMaps << endl; + // clean up mapVec.clear(); mapBlock.clear(); diff --git a/src/classes/PRunListCollection.cpp b/src/classes/PRunListCollection.cpp index c9c5e56f..1a902721 100644 --- a/src/classes/PRunListCollection.cpp +++ b/src/classes/PRunListCollection.cpp @@ -1062,3 +1062,44 @@ const Char_t* PRunListCollection::GetYAxisTitle(const TString &runName, const UI return result; } +//-------------------------------------------------------------------------- +// GetN0Idx (public) +//-------------------------------------------------------------------------- +/** + * \brief PRunListCollection::GetN0Idx + * \param idx + * \return + */ +Int_t PRunListCollection::GetN0Idx(UInt_t idx) +{ + Int_t N0idx = -1; + + // make sure idx is within proper bounds + if (idx >= fRunSingleHistoList.size()) + return N0idx; + + N0idx = fMsrInfo->GetMsrRunList()->at(idx).GetNormParamNo(); + + return N0idx; +} + +//-------------------------------------------------------------------------- +// GetNbkgIdx (public) +//-------------------------------------------------------------------------- +/** + * \brief PRunListCollection::GetNbkgIdx + * \param idx + * \return + */ +Int_t PRunListCollection::GetNbkgIdx(UInt_t idx) +{ + 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; +} diff --git a/src/include/PFitterFcnDKS.h b/src/include/PFitterFcnDKS.h index 8cf2e67f..b6392b29 100644 --- a/src/include/PFitterFcnDKS.h +++ b/src/include/PFitterFcnDKS.h @@ -36,6 +36,11 @@ #include "DKSBaseMuSR.h" #include "PRunListCollection.h" +typedef struct { + UInt_t fN0; ///< N0 parameter index + UInt_t fNbkg; ///< Nbkg parameter index +} PNidx; + /** *

This is the minuit2 interface function class providing the function to be optimized (chisq or log max-likelihood). */ @@ -63,6 +68,9 @@ class PFitterFcnDKS : public ROOT::Minuit2::FCNBase mutable DKSBaseMuSR fDKS; + vector fMemData; ///< vector holding the initial addresses of the data sets on the GPU + vector fNidx; ///< N0 / Nbkg parameter index vector + virtual void InitDKS(const UInt_t dksTag); virtual void FreeDKS(); }; diff --git a/src/include/PMsrHandler.h b/src/include/PMsrHandler.h index 95e3bf76..b0ecfc97 100644 --- a/src/include/PMsrHandler.h +++ b/src/include/PMsrHandler.h @@ -113,6 +113,7 @@ class PMsrHandler virtual std::string GetDKSTheoryString(); virtual UInt_t GetDKSTag(); + virtual Int_t GetNoOfMaps() { return fNoOfMaps; } private: Bool_t fFourierOnly; ///< flag indicating if Fourier transform only is wished. If yes, some part of the msr-file blocks are not needed. @@ -139,6 +140,8 @@ class PMsrHandler Bool_t fCopyStatisticsBlock; ///< flag, if true: just copy to old statistics block (musrt0), otherwise write a new one (musrfit) + Int_t fNoOfMaps; + virtual Bool_t HandleFitParameterEntry(PMsrLines &line); virtual Bool_t HandleTheoryEntry(PMsrLines &line); virtual Bool_t HandleFunctionsEntry(PMsrLines &line); diff --git a/src/include/PRunListCollection.h b/src/include/PRunListCollection.h index 00b1112e..649a745f 100644 --- a/src/include/PRunListCollection.h +++ b/src/include/PRunListCollection.h @@ -100,6 +100,12 @@ class PRunListCollection virtual const Char_t* GetXAxisTitle(const TString &runName, const UInt_t idx) const; virtual const Char_t* GetYAxisTitle(const TString &runName, const UInt_t idx) const; + 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); + private: PMsrHandler *fMsrInfo; ///< pointer to the msr-file handler PRunDataHandler *fData; ///< pointer to the run-data handler