some more work towards DKS/GPU integration

src/classes/PFitterFcnDKS.cpp

@@ -79,17 +79,38 @@ PFitterFcnDKS::~PFitterFcnDKS()
  */
 Double_t PFitterFcnDKS::operator()(const std::vector<Double_t>& 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; i<fRunListCollection->GetNoOfSingleHisto(); 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; i<fRunListCollection->GetNoOfSingleHisto(); 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) {

src/classes/PRunListCollection.cpp

@@ -1063,43 +1063,96 @@ const Char_t* PRunListCollection::GetYAxisTitle(const TString &runName, const UI
 }
 
 //--------------------------------------------------------------------------
-// GetN0Idx (public)
+//  GetSingleHistoParams (public)
 //--------------------------------------------------------------------------
-/**
+Int_t PRunListCollection::GetSingleHistoParams(UInt_t idx, const std::vector<Double_t>& par, PSingleHistoParams &shp)
- * \brief PRunListCollection::GetN0Idx
- * \param idx
- * \return
- */
-Int_t PRunListCollection::GetN0Idx(UInt_t idx)
 {
-  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; i<fMsrInfo->GetNoOfFuncs(); i++) {
+    funcNo = fMsrInfo->GetFuncNo(i);
+    shp.fFun.push_back(fMsrInfo->EvalFunc(funcNo, *runInfo.GetMap(), par));
+  }
+  for (UInt_t i=0; i<shp.fFun.size(); i++)
+    cout << "debug> 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<shp.fMap.size(); i++)
+    cout << "debug> map" << i+1 << " = " << shp.fMap[i] << endl;
+
+  return ierr;
 }
 
 //--------------------------------------------------------------------------
-// GetNbkgIdx (public)
+//  InitSingleHistoParams (private)
 //--------------------------------------------------------------------------
 /**
- * \brief PRunListCollection::GetNbkgIdx
+ * \brief PRunListCollection::InitSingleHistoParams
- * \param idx
+ * \param param
- * \return
  */
-Int_t PRunListCollection::GetNbkgIdx(UInt_t idx)
+void PRunListCollection::InitSingleHistoParams(PSingleHistoParams &param)
 {
-  Int_t NbkgIdx = -1;
+  param.fN0   = -1.0;
-
+  param.fNbkg = -1.0;
-  // make sure idx is within proper bounds
+  param.fTau  = -1.0;
-  if (idx >= fRunSingleHistoList.size())
+  param.fFun.clear();
-    return NbkgIdx;
+  param.fMap.clear();
-
+}
-  NbkgIdx = fMsrInfo->GetMsrRunList()->at(idx).GetBkgFitParamNo();
+
-
+//--------------------------------------------------------------------------
-  return NbkgIdx;
+//  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();
 }

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

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;
+
 /**
  * <p>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 GetSingleHistoParams(UInt_t idx, const std::vector<Double_t>& par, PSingleHistoParams &shp);
-    virtual Int_t GetNbkgIdx(UInt_t idx);
 
   private:
     PMsrHandler *fMsrInfo;  ///< pointer to the msr-file handler
@@ -116,6 +123,9 @@ class PRunListCollection
     vector<PRunAsymmetryRRF*>   fRunAsymmetryRRFList;   ///< stores all processed asymmetry RRF data
     vector<PRunMuMinus*>        fRunMuMinusList;        ///< stores all processed mu-minus data
     vector<PRunNonMusr*>        fRunNonMusrList;        ///< stores all processed non-muSR data
+
+    virtual void InitSingleHistoParams(PSingleHistoParams &param);
+    virtual Double_t GetBackground(Int_t idx);
 };
 
 #endif // _PRUNLISTCOLLECTION_H_

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();