diff --git a/src/classes/PFitterFcnDKS.cpp b/src/classes/PFitterFcnDKS.cpp index 546a9f09..1a3d6f59 100644 --- a/src/classes/PFitterFcnDKS.cpp +++ b/src/classes/PFitterFcnDKS.cpp @@ -87,8 +87,9 @@ Double_t PFitterFcnDKS::operator()(const std::vector& par) const // 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 + // get current values of N0, Nbkg, tau, the functions, the maps, the time resolution, the fit start time, etc. ierr = fRunListCollection->GetSingleHistoParams(i, par, shp); + // set N0, Nbkg ierr += fDKS.callSetConsts(shp.fN0, shp.fTau, shp.fNbkg); @@ -99,13 +100,12 @@ Double_t PFitterFcnDKS::operator()(const std::vector& par) const 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); + ierr += fDKS.callLaunchChiSquare(fMemData[i], fMemData[i], shp.fNoOfFitBins, + par.size(), shp.fFun.size(), shp.fMap.size(), + shp.fStartTime , shp.fPackedTimeResolution, chisq); value += chisq; -*/ + if (ierr != 0) { cerr << "PFitterFcnDKS::operator(): **ERROR** Kernel launch failed!" << endl; exit (EXIT_FAILURE); @@ -113,7 +113,11 @@ Double_t PFitterFcnDKS::operator()(const std::vector& par) const } - return value; + Double_t norm = 1.0; + if (shp.fScaleN0AndBkg) + norm = shp.fPackedTimeResolution*1.0e3; + + return norm*value; } //-------------------------------------------------------------------------- @@ -176,24 +180,25 @@ 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; + // 1) calculated the minimum size for the data needed. + Int_t minSize = 1e9, 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; + size = fRunListCollection->GetNoOfBinsFitted(i); + if (minSize > size) + minSize = size; } for (UInt_t i=0; iGetNoOfAsymmetry(); i++) { - size = fRunListCollection->GetAsymmetry(i)->GetValue()->size(); - if (maxSize < size) - maxSize = size; + size = fRunListCollection->GetNoOfBinsFitted(i); + if (minSize > size) + minSize = size; } - if (maxSize == -1) { + if (minSize == 1e9) { cerr << ">> PFitterFcnDKS::InitDKS: **ERROR** failed to get data size to be fitted." << endl; fValid = false; return; } + cout << "debug> minSize = " << minSize << endl; // 2) get number of parameters / functions / maps parSize = fRunListCollection->GetNoOfParameters(); @@ -217,7 +222,7 @@ void PFitterFcnDKS::InitDKS(const UInt_t dksTag) 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); + ierr = fDKS.initChiSquare(minSize, parSize, funSize, mapSize); if (ierr != 0) { cerr << ">> PFitterFcnDKS::InitDKS: **ERROR** failed to allocate the necessary chisq buffer on the GPU." << endl; fValid = false; @@ -234,13 +239,21 @@ void PFitterFcnDKS::InitDKS(const UInt_t dksTag) fValid = false; return; } - fMemData[i] = fDKS.allocateMemory(runData->GetValue()->size(), ierr); + fMemData[i] = fDKS.allocateMemory(minSize, 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()); + Int_t startTimeBin = fRunListCollection->GetStartTimeBin(i); + if (startTimeBin < 0) { + cerr << ">> PFitterFcnDKS::InitDKS: **ERROR** startTimeBin undefind." << endl; + fValid = false; + return; + } + cout << "debug> InitDKS: startTimeBin = " << startTimeBin << endl; + cout << "debug> InitDKS: runData->GetValue()->size() - startTimeBin = " << runData->GetValue()->size() - startTimeBin << endl; + fDKS.writeData(fMemData[i], &runData->GetValue()->at(startTimeBin), minSize); } // set the function string and compile the program diff --git a/src/classes/PRunListCollection.cpp b/src/classes/PRunListCollection.cpp index 4ac1cbdf..c849bac9 100644 --- a/src/classes/PRunListCollection.cpp +++ b/src/classes/PRunListCollection.cpp @@ -605,7 +605,7 @@ UInt_t PRunListCollection::GetNoOfBinsFitted(const UInt_t idx) const } Int_t type = fMsrInfo->GetMsrRunList()->at(idx).GetFitType(); - if (type == -1) { // i.e. not forun in the RUN block, try the GLOBAL block + if (type == -1) { // i.e. not found in the RUN block, try the GLOBAL block type = fMsrInfo->GetMsrGlobal()->GetFitType(); } @@ -1062,9 +1062,50 @@ const Char_t* PRunListCollection::GetYAxisTitle(const TString &runName, const UI return result; } +//-------------------------------------------------------------------------- +// GetStartTimeBin (public) +//-------------------------------------------------------------------------- +/** + * @brief PRunListCollection::GetStartTimeBin + * @param idx + * @return + */ +Int_t PRunListCollection::GetStartTimeBin(UInt_t idx) +{ + // make sure idx is within proper bounds + if (idx >= fRunSingleHistoList.size()) + return -1; + + return fRunSingleHistoList[idx]->GetStartTimeBin(); +} + +//-------------------------------------------------------------------------- +// GetEndTimeBin (public) +//-------------------------------------------------------------------------- +/** + * @brief PRunListCollection::GetEndTimeBin + * @param idx + * @return + */ +Int_t PRunListCollection::GetEndTimeBin(UInt_t idx) +{ + // make sure idx is within proper bounds + if (idx >= fRunSingleHistoList.size()) + return -1; + + return fRunSingleHistoList[idx]->GetEndTimeBin(); +} + //-------------------------------------------------------------------------- // GetSingleHistoParams (public) //-------------------------------------------------------------------------- +/** + * @brief PRunListCollection::GetSingleHistoParams + * @param idx + * @param par + * @param shp + * @return + */ Int_t PRunListCollection::GetSingleHistoParams(UInt_t idx, const std::vector& par, PSingleHistoParams &shp) { Int_t ierr = 0; @@ -1075,6 +1116,9 @@ Int_t PRunListCollection::GetSingleHistoParams(UInt_t idx, const std::vectorGetScaleN0AndBkg(); + // 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 @@ -1085,14 +1129,14 @@ Int_t PRunListCollection::GetSingleHistoParams(UInt_t idx, const std::vectorEvalFunc(funNo, *runInfo.GetMap(), par); } - cout << "debug> shp.fN0 = " << shp.fN0 << endl; +// 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; +// cout << "debug> shp.fTau = " << shp.fTau << endl; // get background if (runInfo.GetBkgFitParamNo() == -1) { // bkg not fitted @@ -1104,7 +1148,24 @@ Int_t PRunListCollection::GetSingleHistoParams(UInt_t idx, const std::vector shp.fNbkg = " << shp.fNbkg << endl; +// cout << "debug> shp.fNbkg = " << shp.fNbkg << endl; + + // get packed time resolution + shp.fPackedTimeResolution = fRunSingleHistoList[idx]->GetData()->GetDataTimeStep(); + +// cout << "debug> fPackedTimeResolution = " << shp.fPackedTimeResolution*1.0e3 << " (ns)" << endl; + + // get start time + // fRunSingleHistoList[idx]->GetData()->GetDataTimeStart() : time of fgb, which is 0-bin of the fit-data-set + // fRunSingleHistoList[idx]->GetStartTimeBin() * shp.fPackedTimeResolution : time-offset from fgb-time to fit start time + shp.fStartTime = fRunSingleHistoList[idx]->GetData()->GetDataTimeStart() + fRunSingleHistoList[idx]->GetStartTimeBin() * shp.fPackedTimeResolution; + +// cout << "debug> startTime = " << shp.fStartTime << " (us)" << endl; + + // get number of bins fitted + shp.fNoOfFitBins = fRunSingleHistoList[idx]->GetNoOfFitBins(); + +// cout << "debug> startTimeBin = " << fRunSingleHistoList[idx]->GetStartTimeBin() << ", endTimeBin = " << fRunSingleHistoList[idx]->GetEndTimeBin() << ", noOfFittedBins = " << shp.fNoOfFitBins << endl; // calculate functions Int_t funcNo = 0; @@ -1112,14 +1173,21 @@ Int_t PRunListCollection::GetSingleHistoParams(UInt_t idx, const std::vectorGetFuncNo(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()); + // need to reduce map indexes by 1 since in C/C++ arrays start at 0 + for (UInt_t i=0; i map" << i+1 << " = " << shp.fMap[i] << endl; +*/ return ierr; } @@ -1133,9 +1201,13 @@ Int_t PRunListCollection::GetSingleHistoParams(UInt_t idx, const std::vector> PRunSingleHisto::PRunSingleHisto: **SEVERE ERROR**: Couldn't prepare data for fitting!"; cerr << endl << ">> This is very bad :-(, will quit ..."; @@ -174,15 +180,7 @@ Double_t PRunSingleHisto::CalcChiSquare(const std::vector& par) // calculate chi square Double_t time(1.0); - Int_t i, N(static_cast(fData.GetValue()->size())); - - // In order not to have an IF in the next loop, determine the start and end bins for the fit range now - Int_t startTimeBin = static_cast(ceil((fFitStartTime - fData.GetDataTimeStart())/fData.GetDataTimeStep())); - if (startTimeBin < 0) - startTimeBin = 0; - Int_t endTimeBin = static_cast(floor((fFitEndTime - fData.GetDataTimeStart())/fData.GetDataTimeStep())) + 1; - if (endTimeBin > N) - endTimeBin = N; + Int_t i; // Calculate the theory function once to ensure one function evaluation for the current set of parameters. // This is needed for the LF and user functions where some non-thread-save calculations only need to be calculated once @@ -191,12 +189,12 @@ Double_t PRunSingleHisto::CalcChiSquare(const std::vector& par) time = fTheory->Func(time, par, fFuncValues); #ifdef HAVE_GOMP - Int_t chunk = (endTimeBin - startTimeBin)/omp_get_num_procs(); + Int_t chunk = (fEndTimeBin - fStartTimeBin)/omp_get_num_procs(); if (chunk < 10) chunk = 10; #pragma omp parallel for default(shared) private(i,time,diff) schedule(dynamic,chunk) reduction(+:chisq) #endif - for (i=startTimeBin; i < endTimeBin; ++i) { + for (i=fStartTimeBin; iat(i) - (N0*TMath::Exp(-time/tau)*(1.0+fTheory->Func(time, par, fFuncValues))+bkg); @@ -267,15 +265,7 @@ Double_t PRunSingleHisto::CalcChiSquareExpected(const std::vector& par // calculate chi square Double_t time(1.0); - Int_t i, N(static_cast(fData.GetValue()->size())); - - // In order not to have an IF in the next loop, determine the start and end bins for the fit range now - Int_t startTimeBin = static_cast(ceil((fFitStartTime - fData.GetDataTimeStart())/fData.GetDataTimeStep())); - if (startTimeBin < 0) - startTimeBin = 0; - Int_t endTimeBin = static_cast(floor((fFitEndTime - fData.GetDataTimeStart())/fData.GetDataTimeStep())) + 1; - if (endTimeBin > N) - endTimeBin = N; + Int_t i; // Calculate the theory function once to ensure one function evaluation for the current set of parameters. // This is needed for the LF and user functions where some non-thread-save calculations only need to be calculated once @@ -284,12 +274,12 @@ Double_t PRunSingleHisto::CalcChiSquareExpected(const std::vector& par time = fTheory->Func(time, par, fFuncValues); #ifdef HAVE_GOMP - Int_t chunk = (endTimeBin - startTimeBin)/omp_get_num_procs(); + Int_t chunk = (fEndTimeBin - fStartTimeBin)/omp_get_num_procs(); if (chunk < 10) chunk = 10; #pragma omp parallel for default(shared) private(i,time,diff) schedule(dynamic,chunk) reduction(+:chisq) #endif - for (i=startTimeBin; i < endTimeBin; ++i) { + for (i=fStartTimeBin; i < fEndTimeBin; ++i) { time = fData.GetDataTimeStart() + (Double_t)i*fData.GetDataTimeStep(); theo = N0*TMath::Exp(-time/tau)*(1.0+fTheory->Func(time, par, fFuncValues))+bkg; diff = fData.GetValue()->at(i) - theo; @@ -360,7 +350,7 @@ Double_t PRunSingleHisto::CalcMaxLikelihood(const std::vector& par) Double_t theo; Double_t data; Double_t time(1.0); - Int_t i, N(static_cast(fData.GetValue()->size())); + Int_t i; // norm is needed since there is no simple scaling like in chisq case to get the correct Max.Log.Likelihood value when normlizing N(t) to 1/ns Double_t normalizer = 1.0; @@ -368,27 +358,28 @@ Double_t PRunSingleHisto::CalcMaxLikelihood(const std::vector& par) if (fScaleN0AndBkg) normalizer = fPacking * (fTimeResolution * 1.0e3); - // In order not to have an IF in the next loop, determine the start and end bins for the fit range now - Int_t startTimeBin = static_cast(ceil((fFitStartTime - fData.GetDataTimeStart())/fData.GetDataTimeStep())); - if (startTimeBin < 0) - startTimeBin = 0; - Int_t endTimeBin = static_cast(floor((fFitEndTime - fData.GetDataTimeStart())/fData.GetDataTimeStep())) + 1; - if (endTimeBin > N) - endTimeBin = N; - // Calculate the theory function once to ensure one function evaluation for the current set of parameters. // This is needed for the LF and user functions where some non-thread-save calculations only need to be calculated once // for a given set of parameters---which should be done outside of the parallelized loop. // For all other functions it means a tiny and acceptable overhead. time = fTheory->Func(time, par, fFuncValues); +/* +cout << "debug> normalizer=" << normalizer << endl; +for (i=fStartTimeBin; iFunc(time, par, fFuncValues))+bkg; + cout << "debug> " << i << ": time=" << time << ", data=" << fData.GetValue()->at(i) << ", theo=" << theo << endl; +} +*/ + #ifdef HAVE_GOMP - Int_t chunk = (endTimeBin - startTimeBin)/omp_get_num_procs(); + Int_t chunk = (fEndTimeBin - fStartTimeBin)/omp_get_num_procs(); if (chunk < 10) chunk = 10; #pragma omp parallel for default(shared) private(i,time,theo,data) schedule(dynamic,chunk) reduction(-:mllh) #endif - for (i=startTimeBin; iFunc(time, par, fFuncValues))+bkg; @@ -406,6 +397,13 @@ Double_t PRunSingleHisto::CalcMaxLikelihood(const std::vector& par) } else { mllh += (theo-data); } + + if (i<20) { + if (data > 1.0e-9) + cout << "CPU> " << i << ": time=" << time << ", data=" << data << ", theo=" << theo << ", mlh=" << 2.0*((theo-data) + data*log(data/theo)) << endl; + else + cout << "CPU> " << i << ": time=" << time << ", mlh=" << 2.0*(theo-data) << endl; + } } return 2.0*mllh; @@ -589,15 +587,15 @@ void PRunSingleHisto::SetFitRangeBin(const TString fitRange) void PRunSingleHisto::CalcNoOfFitBins() { // In order not having to loop over all bins and to stay consistent with the chisq method, calculate the start and end bins explicitly - Int_t startTimeBin = static_cast(ceil((fFitStartTime - fData.GetDataTimeStart())/fData.GetDataTimeStep())); - if (startTimeBin < 0) - startTimeBin = 0; - Int_t endTimeBin = static_cast(floor((fFitEndTime - fData.GetDataTimeStart())/fData.GetDataTimeStep())) + 1; - if (endTimeBin > static_cast(fData.GetValue()->size())) - endTimeBin = fData.GetValue()->size(); + fStartTimeBin = static_cast(ceil((fFitStartTime - fData.GetDataTimeStart())/fData.GetDataTimeStep())); + if (fStartTimeBin < 0) + fStartTimeBin = 0; + fEndTimeBin = static_cast(floor((fFitEndTime - fData.GetDataTimeStart())/fData.GetDataTimeStep())) + 1; + if (fEndTimeBin > static_cast(fData.GetValue()->size())) + fEndTimeBin = fData.GetValue()->size(); - if (endTimeBin > startTimeBin) - fNoOfFitBins = endTimeBin - startTimeBin; + if (fEndTimeBin > fStartTimeBin) + fNoOfFitBins = fEndTimeBin - fStartTimeBin; else fNoOfFitBins = 0; } @@ -1280,6 +1278,8 @@ Bool_t PRunSingleHisto::PrepareViewData(PRawRunData* runData, const UInt_t histo // clean up par.clear(); +//cout << "debug> fStartTimeBin = " << fStartTimeBin << ", fEndTimeBin = " << fEndTimeBin << endl; + return true; } diff --git a/src/include/PRunBase.h b/src/include/PRunBase.h index 37fad18b..6fb6f529 100644 --- a/src/include/PRunBase.h +++ b/src/include/PRunBase.h @@ -57,6 +57,7 @@ class PRunBase virtual void CalcTheory() = 0; ///< pure virtual, i.e. needs to be implemented by the deriving class!! + virtual Double_t GetTimeResolution() { return fTimeResolution; } ///< returns the native raw data time resolution virtual UInt_t GetRunNo() { return fRunNo; } ///< returns the number of runs of the msr-file virtual PRunData* GetData() { return &fData; } ///< returns the data to be fitted virtual void CleanUp(); diff --git a/src/include/PRunListCollection.h b/src/include/PRunListCollection.h index 7f1d878f..e7bed341 100644 --- a/src/include/PRunListCollection.h +++ b/src/include/PRunListCollection.h @@ -43,14 +43,23 @@ using namespace std; #include "PRunMuMinus.h" #include "PRunNonMusr.h" +//---------------------------------------------------------------------------------------------------------------- +/** + * + */ typedef struct { + Bool_t fScaleN0AndBkg; Double_t fN0; Double_t fNbkg; Double_t fTau; + Double_t fPackedTimeResolution; + Double_t fStartTime; + Int_t fNoOfFitBins; PDoubleVector fFun; PIntVector fMap; } PSingleHistoParams; +//---------------------------------------------------------------------------------------------------------------- /** *

Handler class handling all processed data of an msr-file. All calls of minuit2 are going through this class. */ @@ -111,6 +120,8 @@ 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 GetStartTimeBin(UInt_t idx); + virtual Int_t GetEndTimeBin(UInt_t idx); virtual Int_t GetSingleHistoParams(UInt_t idx, const std::vector& par, PSingleHistoParams &shp); private: diff --git a/src/include/PRunSingleHisto.h b/src/include/PRunSingleHisto.h index 2e736e82..15423636 100644 --- a/src/include/PRunSingleHisto.h +++ b/src/include/PRunSingleHisto.h @@ -53,6 +53,11 @@ class PRunSingleHisto : public PRunBase virtual Double_t GetBackground() { return fBackground; } + virtual Int_t GetStartTimeBin() { return fStartTimeBin; } + virtual Int_t GetEndTimeBin() { return fEndTimeBin; } + virtual Int_t GetPacking() { return fPacking; } + virtual Bool_t GetScaleN0AndBkg() { return fScaleN0AndBkg; } + protected: virtual void CalcNoOfFitBins(); virtual Bool_t PrepareData(); @@ -66,9 +71,12 @@ class PRunSingleHisto : public PRunBase Double_t fBackground; ///< needed if background range is given (units: 1/bin) Int_t fPacking; ///< packing for this particular run. Either given in the RUN- or GLOBAL-block. - Int_t fGoodBins[2]; ///< keep first/last good bins. 0=fgb, 1=lgb + Int_t fGoodBins[2]; ///< keep first/last good bins. 0=fgb, 1=lgb - PDoubleVector fForward; ///< forward histo data + PDoubleVector fForward; ///< forward histo data + + Int_t fStartTimeBin; ///< bin at which the fit starts + Int_t fEndTimeBin; ///< bin at which the fit ends virtual Bool_t GetProperT0(PRawRunData* runData, PMsrGlobalBlock *globalBlock, PUIntVector &histoNo); virtual Bool_t GetProperDataRange();