LMU/musrfit
LMU/musrfit
5
0
Fork 0
Code Issues 6 Pull Requests Actions Packages Projects Releases 3 Activity

added DKS asymmetry fit support

Browse Source
This commit is contained in:
suter_a 2016-04-15 16:47:40 +02:00
parent 27eb664686
commit aedacbcd34
7 changed files with 249 additions and 82 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

112
src/classes/PFitterFcnDKS.cpp
View File

@ -85,37 +85,64 @@ Double_t PFitterFcnDKS::operator()(const std::vector<Double_t>& par) const
Int_t ierr = fDKS.writeParams(&par[0], par.size()); Int_t ierr = fDKS.writeParams(&par[0], par.size());
// loop over all data sets // loop over all data sets
PSingleHistoParams shp; PDKSParams dksp;
Double_t norm = 1.0;
// single histos
for (UInt_t i=0; i<fRunListCollection->GetNoOfSingleHisto(); i++) { for (UInt_t i=0; i<fRunListCollection->GetNoOfSingleHisto(); i++) {
// get current values of N0, Nbkg, tau, the functions, the maps, the time resolution, the fit start time, etc. // 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); ierr = fRunListCollection->GetSingleHistoParams(i, par, dksp);
// set N0, Nbkg // set N0, Nbkg
ierr += fDKS.callSetConsts(shp.fN0, shp.fTau, shp.fNbkg); ierr += fDKS.callSetConsts(dksp.fN0, dksp.fTau, dksp.fNbkg);
// set fun values // set fun values
ierr += fDKS.writeFunctions(&shp.fFun[0], shp.fFun.size()); ierr += fDKS.writeFunctions(&dksp.fFun[0], dksp.fFun.size());
// set map values // set map values
ierr += fDKS.writeMaps(&shp.fMap[0], shp.fMap.size()); ierr += fDKS.writeMaps(&dksp.fMap[0], dksp.fMap.size());
// calc chisq/log-mlh // calc chisq/log-mlh
chisq = 0.0; chisq = 0.0;
ierr += fDKS.callLaunchChiSquare(fMemData[i], fMemData[i], shp.fNoOfFitBins, ierr += fDKS.callLaunchChiSquare(FITTYPE_SINGLE_HISTO, fMemDataSingleHisto[i], fMemDataSingleHisto[i], dksp.fNoOfFitBins,
par.size(), shp.fFun.size(), shp.fMap.size(), par.size(), dksp.fFun.size(), dksp.fMap.size(),
shp.fStartTime , shp.fPackedTimeResolution, chisq); dksp.fStartTime , dksp.fPackedTimeResolution, chisq);
value += chisq; value += chisq;
if (ierr != 0) { if (ierr != 0) {
cerr << "PFitterFcnDKS::operator(): **ERROR** Kernel launch failed!" << endl; cerr << "PFitterFcnDKS::operator(): **ERROR** Kernel launch for single histo failed!" << endl;
exit (EXIT_FAILURE); exit (EXIT_FAILURE);
} }
} }
Double_t norm = 1.0; // asymmetries
if (shp.fScaleN0AndBkg) for (UInt_t i=0; i<fRunListCollection->GetNoOfAsymmetry(); i++) {
norm = shp.fPackedTimeResolution*1.0e3; // get current values of alpha and beta, the functions, the maps, the time resolution, the fit start time, etc.
ierr = fRunListCollection->GetAsymmetryParams(i, par, dksp);
// set alpha and beta
ierr += fDKS.callSetConsts(dksp.fAlpha, dksp.fBeta);
// 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
chisq = 0.0;
ierr += fDKS.callLaunchChiSquare(FITTYPE_ASYMMETRY, fMemDataAsymmetry[i], fMemDataAsymmetryErr[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 asymmetry failed!" << endl;
exit (EXIT_FAILURE);
}
}
if (dksp.fScaleN0AndBkg)
norm = dksp.fPackedTimeResolution*1.0e3;
return norm*value; return norm*value;
} }
@ -220,6 +247,7 @@ void PFitterFcnDKS::InitDKS(const UInt_t dksTag)
} }
// now ready to init the chisq buffer on the GPU // now ready to init the chisq buffer on the GPU
cout << "debug> minSize=" << minSize << ", parSize=" << parSize << ", funSize=" << funSize << ", mapSize=" << mapSize << endl;
ierr = fDKS.initChiSquare(minSize, parSize, funSize, mapSize); ierr = fDKS.initChiSquare(minSize, parSize, funSize, mapSize);
if (ierr != 0) { if (ierr != 0) {
cerr << ">> PFitterFcnDKS::InitDKS: **ERROR** failed to allocate the necessary chisq buffer on the GPU." << endl; cerr << ">> PFitterFcnDKS::InitDKS: **ERROR** failed to allocate the necessary chisq buffer on the GPU." << endl;
@ -228,28 +256,58 @@ void PFitterFcnDKS::InitDKS(const UInt_t dksTag)
} }
// allocate memory for the data on the GPU/CPU and transfer the data sets // allocate memory for the data on the GPU/CPU and transfer the data sets
fMemData.resize(fRunListCollection->GetNoOfSingleHisto());
PRunData *runData=0; PRunData *runData=0;
// single histos
fMemDataSingleHisto.resize(fRunListCollection->GetNoOfSingleHisto());
for (UInt_t i=0; i<fRunListCollection->GetNoOfSingleHisto(); i++) { for (UInt_t i=0; i<fRunListCollection->GetNoOfSingleHisto(); i++) {
runData = fRunListCollection->GetSingleHisto(i); runData = fRunListCollection->GetSingleHisto(i);
if (runData == 0) { if (runData == 0) {
cerr << ">> PFitterFcnDKS::InitDKS: **ERROR** failed to get data set (i=" << i << ") from fRunListCollection." << endl; cerr << ">> PFitterFcnDKS::InitDKS: **ERROR** failed to get single histo data set (i=" << i << ") from fRunListCollection." << endl;
fValid = false; fValid = false;
return; return;
} }
fMemData[i] = fDKS.allocateMemory<Double_t>(minSize, ierr); fMemDataSingleHisto[i] = fDKS.allocateMemory<Double_t>(minSize, ierr);
if (ierr != 0) { if (ierr != 0) {
cerr << ">> PFitterFcnDKS::InitDKS: **ERROR** failed to allocated data set memory (i=" << i << ") on the GPU" << endl; cerr << ">> PFitterFcnDKS::InitDKS: **ERROR** failed to allocated single histo data set memory (i=" << i << ") on the GPU" << endl;
fValid = false; fValid = false;
return; return;
} }
Int_t startTimeBin = fRunListCollection->GetStartTimeBin(i); Int_t startTimeBin = fRunListCollection->GetStartTimeBin(MSR_FITTYPE_SINGLE_HISTO, i);
if (startTimeBin < 0) { if (startTimeBin < 0) {
cerr << ">> PFitterFcnDKS::InitDKS: **ERROR** startTimeBin undefind." << endl; cerr << ">> PFitterFcnDKS::InitDKS: **ERROR** startTimeBin undefind (single histo fit)." << endl;
fValid = false; fValid = false;
return; return;
} }
fDKS.writeData<double>(fMemData[i], &runData->GetValue()->at(startTimeBin), minSize); fDKS.writeData<double>(fMemDataSingleHisto[i], &runData->GetValue()->at(startTimeBin), minSize);
}
// asymmetry
fMemDataAsymmetry.resize(fRunListCollection->GetNoOfAsymmetry());
fMemDataAsymmetryErr.resize(fRunListCollection->GetNoOfAsymmetry());
for (UInt_t i=0; i<fRunListCollection->GetNoOfAsymmetry(); i++) {
runData = fRunListCollection->GetAsymmetry(i);
if (runData == 0) {
cerr << ">> PFitterFcnDKS::InitDKS: **ERROR** failed to get asymmetry data set (i=" << i << ") from fRunListCollection." << endl;
fValid = false;
return;
}
fMemDataAsymmetry[i] = fDKS.allocateMemory<Double_t>(minSize, ierr);
fMemDataAsymmetryErr[i] = fDKS.allocateMemory<Double_t>(minSize, ierr);
if (ierr != 0) {
cerr << ">> PFitterFcnDKS::InitDKS: **ERROR** failed to allocated asymmetry data set memory (i=" << i << ") on the GPU" << endl;
fValid = false;
return;
}
Int_t startTimeBin = fRunListCollection->GetStartTimeBin(MSR_FITTYPE_ASYM, i);
if (startTimeBin < 0) {
cerr << ">> PFitterFcnDKS::InitDKS: **ERROR** startTimeBin undefind (asymmetry fit)." << endl;
fValid = false;
return;
}
fDKS.writeData<double>(fMemDataAsymmetry[i], &runData->GetValue()->at(startTimeBin), minSize);
fDKS.writeData<double>(fMemDataAsymmetryErr[i], &runData->GetError()->at(startTimeBin), minSize);
} }
// set the function string and compile the program // set the function string and compile the program
@ -281,9 +339,17 @@ void PFitterFcnDKS::InitDKS(const UInt_t dksTag)
void PFitterFcnDKS::FreeDKS() void PFitterFcnDKS::FreeDKS()
{ {
PRunData *runData=0; PRunData *runData=0;
for (UInt_t i=0; i<fMemData.size(); i++) { // single histo
for (UInt_t i=0; i<fMemDataSingleHisto.size(); i++) {
runData = fRunListCollection->GetSingleHisto(i); runData = fRunListCollection->GetSingleHisto(i);
fDKS.freeMemory<Double_t>(fMemData[i], runData->GetValue()->size()); fDKS.freeMemory<Double_t>(fMemDataSingleHisto[i], runData->GetValue()->size());
}
fMemDataSingleHisto.clear();
// asymmetry
for (UInt_t i=0; i<fMemDataAsymmetry.size(); i++) {
runData = fRunListCollection->GetAsymmetry(i);
fDKS.freeMemory<Double_t>(fMemDataAsymmetry[i], runData->GetValue()->size());
fDKS.freeMemory<Double_t>(fMemDataAsymmetryErr[i], runData->GetValue()->size());
} }
fMemData.clear();
} }

5
src/classes/PMsrHandler.cpp
View File

@ -5672,8 +5672,10 @@ Bool_t PMsrHandler::CheckMaps()
for (UInt_t i=0; i<map->size(); i++) for (UInt_t i=0; i<map->size(); i++)
if (map->at(i) != 0) if (map->at(i) != 0)
fNoOfMaps++; fNoOfMaps++;
/*as
if (fNoOfMaps == 0) if (fNoOfMaps == 0)
fNoOfMaps = -1; fNoOfMaps = -1;
*/
} }
// clean up // clean up
@ -6238,11 +6240,10 @@ std::string PMsrHandler::GetDKSTheoryString()
} }
} }
/*
cout << "debug> ++++" << endl; cout << "debug> ++++" << endl;
cout << "debug> " << result << endl; cout << "debug> " << result << endl;
cout << "debug> ++++" << endl; cout << "debug> ++++" << endl;
*/
return result; return result;
} }

36
src/classes/PRunAsymmetry.cpp
View File

@ -62,6 +62,9 @@ PRunAsymmetry::PRunAsymmetry() : PRunBase()
// the fit range can be changed in the command block, these variables need to be accessible // the fit range can be changed in the command block, these variables need to be accessible
fGoodBins[0] = -1; fGoodBins[0] = -1;
fGoodBins[1] = -1; fGoodBins[1] = -1;
fStartTimeBin = -1;
fEndTimeBin = -1;
} }
//-------------------------------------------------------------------------- //--------------------------------------------------------------------------
@ -82,6 +85,9 @@ PRunAsymmetry::PRunAsymmetry(PMsrHandler *msrInfo, PRunDataHandler *rawData, UIn
fGoodBins[0] = -1; fGoodBins[0] = -1;
fGoodBins[1] = -1; fGoodBins[1] = -1;
fStartTimeBin = -1;
fEndTimeBin = -1;
fPacking = fRunInfo->GetPacking(); fPacking = fRunInfo->GetPacking();
if (fPacking == -1) { // i.e. packing is NOT given in the RUN-block, it must be given in the GLOBAL-block if (fPacking == -1) { // i.e. packing is NOT given in the RUN-block, it must be given in the GLOBAL-block
fPacking = fMsrInfo->GetMsrGlobal()->GetPacking(); fPacking = fMsrInfo->GetMsrGlobal()->GetPacking();
@ -189,15 +195,7 @@ Double_t PRunAsymmetry::CalcChiSquare(const std::vector<Double_t>& par)
// calculate chi square // calculate chi square
Double_t time(1.0); Double_t time(1.0);
Int_t i, N(static_cast<Int_t>(fData.GetValue()->size())); Int_t i;
// 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<Int_t>(ceil((fFitStartTime - fData.GetDataTimeStart())/fData.GetDataTimeStep()));
if (startTimeBin < 0)
startTimeBin = 0;
Int_t endTimeBin = static_cast<Int_t>(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. // 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 // This is needed for the LF and user functions where some non-thread-save calculations only need to be calculated once
@ -206,12 +204,12 @@ Double_t PRunAsymmetry::CalcChiSquare(const std::vector<Double_t>& par)
asymFcnValue = fTheory->Func(time, par, fFuncValues); asymFcnValue = fTheory->Func(time, par, fFuncValues);
#ifdef HAVE_GOMP #ifdef HAVE_GOMP
Int_t chunk = (endTimeBin - startTimeBin)/omp_get_num_procs(); Int_t chunk = (fEndTimeBin - fStartTimeBin)/omp_get_num_procs();
if (chunk < 10) if (chunk < 10)
chunk = 10; chunk = 10;
#pragma omp parallel for default(shared) private(i,time,diff,asymFcnValue,a,b,f) schedule(dynamic,chunk) reduction(+:chisq) #pragma omp parallel for default(shared) private(i,time,diff,asymFcnValue,a,b,f) schedule(dynamic,chunk) reduction(+:chisq)
#endif #endif
for (i=startTimeBin; i < endTimeBin; ++i) { for (i=fStartTimeBin; i<fEndTimeBin; ++i) {
time = fData.GetDataTimeStart() + (Double_t)i*fData.GetDataTimeStep(); time = fData.GetDataTimeStart() + (Double_t)i*fData.GetDataTimeStep();
switch (fAlphaBetaTag) { switch (fAlphaBetaTag) {
case 1: // alpha == 1, beta == 1 case 1: // alpha == 1, beta == 1
@ -389,15 +387,15 @@ void PRunAsymmetry::SetFitRangeBin(const TString fitRange)
void PRunAsymmetry::CalcNoOfFitBins() void PRunAsymmetry::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 // 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<Int_t>(ceil((fFitStartTime - fData.GetDataTimeStart())/fData.GetDataTimeStep())); fStartTimeBin = static_cast<Int_t>(ceil((fFitStartTime - fData.GetDataTimeStart())/fData.GetDataTimeStep()));
if (startTimeBin < 0) if (fStartTimeBin < 0)
startTimeBin = 0; fStartTimeBin = 0;
Int_t endTimeBin = static_cast<Int_t>(floor((fFitEndTime - fData.GetDataTimeStart())/fData.GetDataTimeStep())) + 1; fEndTimeBin = static_cast<Int_t>(floor((fFitEndTime - fData.GetDataTimeStart())/fData.GetDataTimeStep())) + 1;
if (endTimeBin > static_cast<Int_t>(fData.GetValue()->size())) if (fEndTimeBin > static_cast<Int_t>(fData.GetValue()->size()))
endTimeBin = fData.GetValue()->size(); fEndTimeBin = fData.GetValue()->size();
if (endTimeBin > startTimeBin) if (fEndTimeBin > fStartTimeBin)
fNoOfFitBins = endTimeBin - startTimeBin; fNoOfFitBins = fEndTimeBin - fStartTimeBin;
else else
fNoOfFitBins = 0; fNoOfFitBins = 0;
} }

149
src/classes/PRunListCollection.cpp
View File

@ -1067,16 +1067,28 @@ const Char_t* PRunListCollection::GetYAxisTitle(const TString &runName, const UI
//-------------------------------------------------------------------------- //--------------------------------------------------------------------------
/** /**
* @brief PRunListCollection::GetStartTimeBin * @brief PRunListCollection::GetStartTimeBin
* @param fitType
* @param idx * @param idx
* @return * @return
*/ */
Int_t PRunListCollection::GetStartTimeBin(UInt_t idx) Int_t PRunListCollection::GetStartTimeBin(Int_t fitType, UInt_t idx)
{ {
// make sure idx is within proper bounds Int_t result = -1;
if (idx >= fRunSingleHistoList.size())
return -1;
return fRunSingleHistoList[idx]->GetStartTimeBin(); switch (fitType) {
case MSR_FITTYPE_SINGLE_HISTO:
if (idx < fRunSingleHistoList.size())
result = fRunSingleHistoList[idx]->GetStartTimeBin();
break;
case MSR_FITTYPE_ASYM:
if (idx < fRunAsymmetryList.size())
result = fRunAsymmetryList[idx]->GetStartTimeBin();
break;
default:
break;
}
return result;
} }
//-------------------------------------------------------------------------- //--------------------------------------------------------------------------
@ -1084,16 +1096,28 @@ Int_t PRunListCollection::GetStartTimeBin(UInt_t idx)
//-------------------------------------------------------------------------- //--------------------------------------------------------------------------
/** /**
* @brief PRunListCollection::GetEndTimeBin * @brief PRunListCollection::GetEndTimeBin
* @param fitType
* @param idx * @param idx
* @return * @return
*/ */
Int_t PRunListCollection::GetEndTimeBin(UInt_t idx) Int_t PRunListCollection::GetEndTimeBin(Int_t fitType, UInt_t idx)
{ {
// make sure idx is within proper bounds Int_t result = -1;
if (idx >= fRunSingleHistoList.size())
return -1;
return fRunSingleHistoList[idx]->GetEndTimeBin(); switch (fitType) {
case MSR_FITTYPE_SINGLE_HISTO:
if (idx < fRunSingleHistoList.size())
result = fRunSingleHistoList[idx]->GetEndTimeBin();
break;
case MSR_FITTYPE_ASYM:
if (idx < fRunAsymmetryList.size())
result = fRunAsymmetryList[idx]->GetEndTimeBin();
break;
default:
break;
}
return result;
} }
//-------------------------------------------------------------------------- //--------------------------------------------------------------------------
@ -1106,7 +1130,7 @@ Int_t PRunListCollection::GetEndTimeBin(UInt_t idx)
* @param shp * @param shp
* @return * @return
*/ */
Int_t PRunListCollection::GetSingleHistoParams(UInt_t idx, const std::vector<Double_t>& par, PSingleHistoParams &shp) Int_t PRunListCollection::GetSingleHistoParams(UInt_t idx, const std::vector<Double_t>& par, PDKSParams &dksp)
{ {
Int_t ierr = 0; Int_t ierr = 0;
// make sure idx is within proper bounds // make sure idx is within proper bounds
@ -1114,80 +1138,141 @@ Int_t PRunListCollection::GetSingleHistoParams(UInt_t idx, const std::vector<Dou
return 1; return 1;
// init param // init param
InitSingleHistoParams(shp); InitDKSParams(dksp);
// get flag if scaling of N0 and Nbkg is wished // get flag if scaling of N0 and Nbkg is wished
shp.fScaleN0AndBkg = fRunSingleHistoList[idx]->GetScaleN0AndBkg(); dksp.fScaleN0AndBkg = fRunSingleHistoList[idx]->GetScaleN0AndBkg();
// check if norm is a parameter or a function // check if norm is a parameter or a function
PMsrRunBlock runInfo = fMsrInfo->GetMsrRunList()->at(idx); PMsrRunBlock runInfo = fMsrInfo->GetMsrRunList()->at(idx);
if (runInfo.GetNormParamNo() < MSR_PARAM_FUN_OFFSET) { // norm is a parameter if (runInfo.GetNormParamNo() < MSR_PARAM_FUN_OFFSET) { // norm is a parameter
shp.fN0 = par[runInfo.GetNormParamNo()-1]; dksp.fN0 = par[runInfo.GetNormParamNo()-1];
} else { // norm is a function } else { // norm is a function
// get function number // get function number
UInt_t funNo = runInfo.GetNormParamNo()-MSR_PARAM_FUN_OFFSET; UInt_t funNo = runInfo.GetNormParamNo()-MSR_PARAM_FUN_OFFSET;
// evaluate function // evaluate function
shp.fN0 = fMsrInfo->EvalFunc(funNo, *runInfo.GetMap(), par); dksp.fN0 = fMsrInfo->EvalFunc(funNo, *runInfo.GetMap(), par);
} }
// get tau // get tau
if (runInfo.GetLifetimeParamNo() != -1) if (runInfo.GetLifetimeParamNo() != -1)
shp.fTau = par[runInfo.GetLifetimeParamNo()-1]; dksp.fTau = par[runInfo.GetLifetimeParamNo()-1];
else else
shp.fTau = PMUON_LIFETIME; dksp.fTau = PMUON_LIFETIME;
// get background // get background
if (runInfo.GetBkgFitParamNo() == -1) { // bkg not fitted if (runInfo.GetBkgFitParamNo() == -1) { // bkg not fitted
if (runInfo.GetBkgFix(0) == PMUSR_UNDEFINED) { // no fixed background given (background interval) if (runInfo.GetBkgFix(0) == PMUSR_UNDEFINED) { // no fixed background given (background interval)
shp.fNbkg = GetBackground(idx); dksp.fNbkg = GetBackground(idx);
} else { // fixed bkg given } else { // fixed bkg given
shp.fNbkg = runInfo.GetBkgFix(0); dksp.fNbkg = runInfo.GetBkgFix(0);
} }
} else { // bkg fitted } else { // bkg fitted
shp.fNbkg = par[runInfo.GetBkgFitParamNo()-1]; dksp.fNbkg = par[runInfo.GetBkgFitParamNo()-1];
} }
// get packed time resolution // get packed time resolution
shp.fPackedTimeResolution = fRunSingleHistoList[idx]->GetData()->GetDataTimeStep(); dksp.fPackedTimeResolution = fRunSingleHistoList[idx]->GetData()->GetDataTimeStep();
// get start time // get start time
// fRunSingleHistoList[idx]->GetData()->GetDataTimeStart() : time of fgb, which is 0-bin of the fit-data-set // 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 // fRunSingleHistoList[idx]->GetStartTimeBin() * dksp.fPackedTimeResolution : time-offset from fgb-time to fit start time
shp.fStartTime = fRunSingleHistoList[idx]->GetData()->GetDataTimeStart() + fRunSingleHistoList[idx]->GetStartTimeBin() * shp.fPackedTimeResolution; dksp.fStartTime = fRunSingleHistoList[idx]->GetData()->GetDataTimeStart() + fRunSingleHistoList[idx]->GetStartTimeBin() * dksp.fPackedTimeResolution;
// get number of bins fitted // get number of bins fitted
shp.fNoOfFitBins = fRunSingleHistoList[idx]->GetNoOfFitBins(); dksp.fNoOfFitBins = fRunSingleHistoList[idx]->GetNoOfFitBins();
// calculate functions // calculate functions
Int_t funcNo = 0; Int_t funcNo = 0;
for (Int_t i=0; i<fMsrInfo->GetNoOfFuncs(); i++) { for (Int_t i=0; i<fMsrInfo->GetNoOfFuncs(); i++) {
funcNo = fMsrInfo->GetFuncNo(i); funcNo = fMsrInfo->GetFuncNo(i);
shp.fFun.push_back(fMsrInfo->EvalFunc(funcNo, *runInfo.GetMap(), par)); dksp.fFun.push_back(fMsrInfo->EvalFunc(funcNo, *runInfo.GetMap(), par));
} }
// get map vector // get map vector
shp.fMap = *runInfo.GetMap(); dksp.fMap = *runInfo.GetMap();
shp.fMap.erase(shp.fMap.begin()+GetNoOfMaps(), shp.fMap.end()); dksp.fMap.erase(dksp.fMap.begin()+GetNoOfMaps(), dksp.fMap.end());
// need to reduce map indexes by 1 since in C/C++ arrays start at 0 // need to reduce map indexes by 1 since in C/C++ arrays start at 0
for (UInt_t i=0; i<shp.fMap.size(); i++) for (UInt_t i=0; i<dksp.fMap.size(); i++)
shp.fMap[i] -= 1; dksp.fMap[i] -= 1;
return ierr; return ierr;
} }
//-------------------------------------------------------------------------- //--------------------------------------------------------------------------
// InitSingleHistoParams (private) // GetAsymmetryParams (public)
//-------------------------------------------------------------------------- //--------------------------------------------------------------------------
/** /**
* \brief PRunListCollection::InitSingleHistoParams * @brief PRunListCollection::GetAsymmetryParams
* @param idx
* @param par
* @param shp
* @return
*/
Int_t PRunListCollection::GetAsymmetryParams(UInt_t idx, const std::vector<Double_t>& par, PDKSParams &dksp)
{
Int_t ierr=0, ival=0;
// make sure idx is within proper bounds
if (idx >= fRunAsymmetryList.size())
return 1;
// init param
InitDKSParams(dksp);
// get alpha
PMsrRunBlock runInfo = fMsrInfo->GetMsrRunList()->at(idx);
ival = runInfo.GetAlphaParamNo();
if (ival > 0)
dksp.fAlpha = par[ival-1];
// get beta
ival = runInfo.GetBetaParamNo();
if (ival > 0)
dksp.fBeta = par[ival-1];
// get packed time resolution
dksp.fPackedTimeResolution = fRunAsymmetryList[idx]->GetData()->GetDataTimeStep();
// get start time
// fRunAsymmetryList[idx]->GetData()->GetDataTimeStart() : time of fgb, which is 0-bin of the fit-data-set
// fRunAsymmetryList[idx]->GetStartTimeBin() * dksp.fPackedTimeResolution : time-offset from fgb-time to fit start time
dksp.fStartTime = fRunAsymmetryList[idx]->GetData()->GetDataTimeStart() + fRunAsymmetryList[idx]->GetStartTimeBin() * dksp.fPackedTimeResolution;
// get number of bins fitted
dksp.fNoOfFitBins = fRunAsymmetryList[idx]->GetNoOfFitBins();
// calculate functions
Int_t funcNo = 0;
for (Int_t i=0; i<fMsrInfo->GetNoOfFuncs(); i++) {
funcNo = fMsrInfo->GetFuncNo(i);
dksp.fFun.push_back(fMsrInfo->EvalFunc(funcNo, *runInfo.GetMap(), par));
}
// get map vector
dksp.fMap = *runInfo.GetMap();
dksp.fMap.erase(dksp.fMap.begin()+GetNoOfMaps(), dksp.fMap.end());
// need to reduce map indexes by 1 since in C/C++ arrays start at 0
for (UInt_t i=0; i<dksp.fMap.size(); i++)
dksp.fMap[i] -= 1;
return ierr;
}
//--------------------------------------------------------------------------
// InitDKSParams (private)
//--------------------------------------------------------------------------
/**
* \brief PRunListCollection::InitDKSParams
* \param param * \param param
*/ */
void PRunListCollection::InitSingleHistoParams(PSingleHistoParams &param) void PRunListCollection::InitDKSParams(PDKSParams &param)
{ {
param.fScaleN0AndBkg = false; param.fScaleN0AndBkg = false;
param.fN0 = -1.0; param.fN0 = -1.0;
param.fNbkg = -1.0; param.fNbkg = -1.0;
param.fTau = -1.0; param.fTau = -1.0;
param.fAlpha = 1.0;
param.fBeta = 1.0;
param.fPackedTimeResolution = -1.0; param.fPackedTimeResolution = -1.0;
param.fStartTime = -1.0; param.fStartTime = -1.0;
param.fNoOfFitBins = -1; param.fNoOfFitBins = -1;

9
src/include/PFitterFcnDKS.h
View File

@ -36,6 +36,11 @@
#include "DKSBaseMuSR.h" #include "DKSBaseMuSR.h"
#include "PRunListCollection.h" #include "PRunListCollection.h"
#define FITTYPE_UNDEFINED 0
#define FITTYPE_SINGLE_HISTO 1
#define FITTYPE_ASYMMETRY 2
#define FITTYPE_MU_MINUS 3
typedef struct { typedef struct {
UInt_t fN0; ///< N0 parameter index UInt_t fN0; ///< N0 parameter index
UInt_t fNbkg; ///< Nbkg parameter index UInt_t fNbkg; ///< Nbkg parameter index
@ -68,7 +73,9 @@ class PFitterFcnDKS : public ROOT::Minuit2::FCNBase
mutable DKSBaseMuSR fDKS; mutable DKSBaseMuSR fDKS;
vector<void *> fMemData; ///< vector holding the initial addresses of the data sets on the GPU vector<void *> fMemDataSingleHisto; ///< vector holding the initial addresses of the single histo data sets on the GPU
vector<void *> fMemDataAsymmetry; ///< vector holding the initial addresses of the asymmetry data sets on the GPU
vector<void *> fMemDataAsymmetryErr; ///< vector holding the initial addresses of the asymmetry error sets on the GPU
vector<PNidx> fNidx; ///< N0 / Nbkg parameter index vector vector<PNidx> fNidx; ///< N0 / Nbkg parameter index vector
virtual void InitDKS(const UInt_t dksTag); virtual void InitDKS(const UInt_t dksTag);

7
src/include/PRunAsymmetry.h
View File

@ -52,6 +52,10 @@ class PRunAsymmetry : public PRunBase
virtual void SetFitRangeBin(const TString fitRange); virtual void SetFitRangeBin(const TString fitRange);
virtual Int_t GetStartTimeBin() { return fStartTimeBin; }
virtual Int_t GetEndTimeBin() { return fEndTimeBin; }
virtual Int_t GetPacking() { return fPacking; }
protected: protected:
virtual void CalcNoOfFitBins(); virtual void CalcNoOfFitBins();
virtual Bool_t PrepareData(); virtual Bool_t PrepareData();
@ -71,6 +75,9 @@ class PRunAsymmetry : public PRunBase
Int_t fGoodBins[4]; ///< keep first/last good bins. 0=fgb, 1=lgb (forward); 2=fgb, 3=lgb (backward) Int_t fGoodBins[4]; ///< keep first/last good bins. 0=fgb, 1=lgb (forward); 2=fgb, 3=lgb (backward)
Int_t fStartTimeBin; ///< bin at which the fit starts
Int_t fEndTimeBin; ///< bin at which the fit ends
Bool_t SubtractFixBkg(); Bool_t SubtractFixBkg();
Bool_t SubtractEstimatedBkg(); Bool_t SubtractEstimatedBkg();

13
src/include/PRunListCollection.h
View File

@ -52,12 +52,14 @@ typedef struct {
Double_t fN0; Double_t fN0;
Double_t fNbkg; Double_t fNbkg;
Double_t fTau; Double_t fTau;
Double_t fAlpha;
Double_t fBeta;
Double_t fPackedTimeResolution; Double_t fPackedTimeResolution;
Double_t fStartTime; Double_t fStartTime;
Int_t fNoOfFitBins; Int_t fNoOfFitBins;
PDoubleVector fFun; PDoubleVector fFun;
PIntVector fMap; PIntVector fMap;
} PSingleHistoParams; } PDKSParams;
//---------------------------------------------------------------------------------------------------------------- //----------------------------------------------------------------------------------------------------------------
/** /**
@ -120,9 +122,10 @@ class PRunListCollection
virtual Int_t GetNoOfParameters() { return fMsrInfo->GetNoOfParams(); } virtual Int_t GetNoOfParameters() { return fMsrInfo->GetNoOfParams(); }
virtual Int_t GetNoOfFunctions() { return fMsrInfo->GetNoOfFuncs(); } virtual Int_t GetNoOfFunctions() { return fMsrInfo->GetNoOfFuncs(); }
virtual Int_t GetNoOfMaps() { return fMsrInfo->GetNoOfMaps(); } virtual Int_t GetNoOfMaps() { return fMsrInfo->GetNoOfMaps(); }
virtual Int_t GetStartTimeBin(UInt_t idx); virtual Int_t GetStartTimeBin(Int_t fitType, UInt_t idx);
virtual Int_t GetEndTimeBin(UInt_t idx); virtual Int_t GetEndTimeBin(Int_t fitType, UInt_t idx);
virtual Int_t GetSingleHistoParams(UInt_t idx, const std::vector<Double_t>& par, PSingleHistoParams &shp); virtual Int_t GetSingleHistoParams(UInt_t idx, const std::vector<Double_t>& par, PDKSParams &dksp);
virtual Int_t GetAsymmetryParams(UInt_t idx, const std::vector<Double_t>& par, PDKSParams &dksp);
private: private:
PMsrHandler *fMsrInfo; ///< pointer to the msr-file handler PMsrHandler *fMsrInfo; ///< pointer to the msr-file handler
@ -135,7 +138,7 @@ class PRunListCollection
vector<PRunMuMinus*> fRunMuMinusList; ///< stores all processed mu-minus data vector<PRunMuMinus*> fRunMuMinusList; ///< stores all processed mu-minus data
vector<PRunNonMusr*> fRunNonMusrList; ///< stores all processed non-muSR data vector<PRunNonMusr*> fRunNonMusrList; ///< stores all processed non-muSR data
virtual void InitSingleHistoParams(PSingleHistoParams &param); virtual void InitDKSParams(PDKSParams &param);
virtual Double_t GetBackground(Int_t idx); virtual Double_t GetBackground(Int_t idx);
}; };