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

first implementation of a RRF. So far, for single histo only.

Browse Source
This commit is contained in:
nemu 2009-12-23 13:43:23 +00:00
parent 19dee34c22
commit e0c893503f
10 changed files with 444 additions and 240 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

283
src/classes/PMsrHandler.cpp
View File

@ -599,16 +599,6 @@ Int_t PMsrHandler::WriteMsrLogFile(const Bool_t messages)
default:
break;
}
} else if (sstr.BeginsWith("rrffrequency")) {
fout.width(16);
fout << left << "rrffrequency";
fout.precision(prec);
fout << fRuns[runNo].GetRRFFreq() << endl;
} else if (sstr.BeginsWith("rrfpacking")) {
fout.width(16);
fout << left << "rrfpacking";
fout.precision(prec);
fout << fRuns[runNo].GetRRFPacking() << endl;
} else if (sstr.BeginsWith("alpha ")) {
fout.width(16);
fout << left << "alpha";
@ -617,14 +607,6 @@ Int_t PMsrHandler::WriteMsrLogFile(const Bool_t messages)
fout.width(16);
fout << left << "beta";
fout << fRuns[runNo].GetBetaParamNo() << endl;
} else if (sstr.BeginsWith("alpha2")) {
fout.width(16);
fout << left << "alpha2";
fout << fRuns[runNo].GetAlpha2ParamNo() << endl;
} else if (sstr.BeginsWith("beta2")) {
fout.width(16);
fout << left << "beta2";
fout << fRuns[runNo].GetBeta2ParamNo() << endl;
} else if (sstr.BeginsWith("norm")) {
fout.width(16);
fout << left << "norm";
@ -670,14 +652,6 @@ Int_t PMsrHandler::WriteMsrLogFile(const Bool_t messages)
fout.width(16);
fout << left << "backward";
fout << fRuns[runNo].GetBackwardHistoNo() << endl;
} else if (sstr.BeginsWith("right")) {
fout.width(16);
fout << left << "right";
fout << fRuns[runNo].GetRightHistoNo() << endl;
} else if (sstr.BeginsWith("left")) {
fout.width(16);
fout << left << "left";
fout << fRuns[runNo].GetLeftHistoNo() << endl;
} else if (sstr.BeginsWith("backgr.fix")) {
fout.width(15);
fout << left << "backgr.fix";
@ -1731,28 +1705,6 @@ Bool_t PMsrHandler::HandleRunEntry(PMsrLines &lines)
}
}
// rphase ------------------------------------------------
if (iter->fLine.BeginsWith("rphase", TString::kIgnoreCase)) {
runLinePresent = false; // this is needed to make sure that a run line is present before and ADDRUN is following
if (tokens->GetEntries() < 2) {
error = true;
} else {
ostr = dynamic_cast<TObjString*>(tokens->At(1));
str = ostr->GetString();
if (str.IsDigit()) {
dval = str.Atoi();
if (dval > 0)
param.SetPhaseParamNo(dval);
else
error = true;
} else {
error = true;
}
}
}
// lifetime ------------------------------------------------
if (iter->fLine.BeginsWith("lifetime ", TString::kIgnoreCase)) {
@ -1987,133 +1939,6 @@ Bool_t PMsrHandler::HandleRunEntry(PMsrLines &lines)
}
}
// rrffrequency --------------------------------------------------
if (iter->fLine.BeginsWith("rrffrequency", TString::kIgnoreCase)) {
runLinePresent = false; // this is needed to make sure that a run line is present before and ADDRUN is following
if (tokens->GetEntries() != 2) {
error = true;
} else {
ostr = dynamic_cast<TObjString*>(tokens->At(1));
str = ostr->GetString();
if (str.IsFloat())
param.SetRRFFreq(str.Atof());
else
error = true;
}
}
// rrfpacking --------------------------------------------------
if (iter->fLine.BeginsWith("rrfpacking", TString::kIgnoreCase)) {
runLinePresent = false; // this is needed to make sure that a run line is present before and ADDRUN is following
if (tokens->GetEntries() != 2) {
error = true;
} else {
ostr = dynamic_cast<TObjString*>(tokens->At(1));
str = ostr->GetString();
if (str.IsDigit()) {
dval = str.Atoi();
if (dval > 0)
param.SetRRFPacking(dval);
else
error = true;
} else {
error = true;
}
}
}
// alpha2 --------------------------------------------------
if (iter->fLine.BeginsWith("alpha2", TString::kIgnoreCase)) {
runLinePresent = false; // this is needed to make sure that a run line is present before and ADDRUN is following
if (tokens->GetEntries() != 2) {
error = true;
} else {
ostr = dynamic_cast<TObjString*>(tokens->At(1));
str = ostr->GetString();
if (str.IsDigit()) {
dval = str.Atoi();
if (dval > 0)
param.SetAlpha2ParamNo(dval);
else
error = true;
} else {
error = true;
}
}
}
// beta2 --------------------------------------------------
if (iter->fLine.BeginsWith("beta2", TString::kIgnoreCase)) {
runLinePresent = false; // this is needed to make sure that a run line is present before and ADDRUN is following
if (tokens->GetEntries() != 2) {
error = true;
} else {
ostr = dynamic_cast<TObjString*>(tokens->At(1));
str = ostr->GetString();
if (str.IsDigit()) {
dval = str.Atoi();
if (dval > 0)
param.SetBeta2ParamNo(dval);
else
error = true;
} else {
error = true;
}
}
}
// right --------------------------------------------------
if (iter->fLine.BeginsWith("right", TString::kIgnoreCase)) {
runLinePresent = false; // this is needed to make sure that a run line is present before and ADDRUN is following
if (tokens->GetEntries() != 2) {
error = true;
} else {
ostr = dynamic_cast<TObjString*>(tokens->At(1));
str = ostr->GetString();
if (str.IsDigit()) {
dval = str.Atoi();
if (dval > 0)
param.SetRightHistoNo(dval);
else
error = true;
} else {
error = true;
}
}
}
// left --------------------------------------------------
if (iter->fLine.BeginsWith("left", TString::kIgnoreCase)) {
runLinePresent = false; // this is needed to make sure that a run line is present before and ADDRUN is following
if (tokens->GetEntries() != 2) {
error = true;
} else {
ostr = dynamic_cast<TObjString*>(tokens->At(1));
str = ostr->GetString();
if (str.IsDigit()) {
dval = str.Atoi();
if (dval > 0)
param.SetLeftHistoNo(str.Atoi());
else
error = true;
} else {
error = true;
}
}
}
// xy-data -----------------------------------------------
if (iter->fLine.BeginsWith("xy-data", TString::kIgnoreCase)) {
@ -2536,6 +2361,10 @@ Bool_t PMsrHandler::HandlePlotEntry(PMsrLines &lines)
param.fLogX = false; // i.e. if not overwritten use linear x-axis
param.fLogY = false; // i.e. if not overwritten use linear y-axis
param.fViewPacking = -1; // i.e. if not overwritten use the packing of the run blocks
param.fRRFPacking = 0; // i.e. if not overwritten it will not be a valid RRF
param.fRRFFreq = 0.0; // i.e. no RRF whished
param.fRRFUnit = RRF_UNIT_MHz;
param.fRRFPhase = 0.0;
// find next plot if any is present
iter2 = iter1;
@ -2777,7 +2606,7 @@ Bool_t PMsrHandler::HandlePlotEntry(PMsrLines &lines)
} else if (iter1->fLine.Contains("view_packing", TString::kIgnoreCase)) {
tokens = iter1->fLine.Tokenize(" \t");
if (!tokens) {
cerr << endl << ">> PMsrHandler::HandlePlotEntry: **SEVERE ERROR** Couldn't tokenize PLOT in line " << iter1->fLineNo;
cerr << endl << ">> PMsrHandler::HandlePlotEntry: **SEVERE ERROR** Couldn't tokenize view_packing in line " << iter1->fLineNo;
cerr << endl << endl;
return false;
}
@ -2797,6 +2626,97 @@ Bool_t PMsrHandler::HandlePlotEntry(PMsrLines &lines)
}
}
// clean up
if (tokens) {
delete tokens;
tokens = 0;
}
} else if (iter1->fLine.Contains("rrf_freq", TString::kIgnoreCase)) {
// expected entry: rrf_freq value unit
// allowed units: kHz, MHz, Mc/s, G, T
tokens = iter1->fLine.Tokenize(" \t");
if (!tokens) {
cerr << endl << ">> PMsrHandler::HandlePlotEntry: **SEVERE ERROR** Couldn't tokenize rrf_freq in line " << iter1->fLineNo;
cerr << endl << endl;
return false;
}
if (tokens->GetEntries() != 3) {
error = true;
} else {
// get rrf frequency
ostr = dynamic_cast<TObjString*>(tokens->At(1));
str = ostr->GetString();
if (str.IsFloat()) {
param.fRRFFreq = str.Atof();
} else {
error = true;
}
// get unit
ostr = dynamic_cast<TObjString*>(tokens->At(2));
str = ostr->GetString();
if (str.Contains("kHz", TString::kIgnoreCase))
param.fRRFUnit = RRF_UNIT_kHz;
else if (str.Contains("MHz", TString::kIgnoreCase))
param.fRRFUnit = RRF_UNIT_MHz;
else if (str.Contains("Mc/s", TString::kIgnoreCase))
param.fRRFUnit = RRF_UNIT_Mcs;
else if (str.Contains("G", TString::kIgnoreCase))
param.fRRFUnit = RRF_UNIT_G;
else if (str.Contains("T", TString::kIgnoreCase))
param.fRRFUnit = RRF_UNIT_T;
else
error = true;
}
// clean up
if (tokens) {
delete tokens;
tokens = 0;
}
} else if (iter1->fLine.Contains("rrf_phase", TString::kIgnoreCase)) {
// expected entry: rrf_phase value. value given in units of degree
tokens = iter1->fLine.Tokenize(" \t");
if (!tokens) {
cerr << endl << ">> PMsrHandler::HandlePlotEntry: **SEVERE ERROR** Couldn't tokenize rrf_phase in line " << iter1->fLineNo;
cerr << endl << endl;
return false;
}
if (tokens->GetEntries() != 2) {
error = true;
} else {
// get rrf phase
ostr = dynamic_cast<TObjString*>(tokens->At(1));
str = ostr->GetString();
if (str.IsFloat()) {
param.fRRFPhase = str.Atof();
} else {
error = true;
}
}
// clean up
if (tokens) {
delete tokens;
tokens = 0;
}
} else if (iter1->fLine.Contains("rrf_packing", TString::kIgnoreCase)) {
// expected entry: rrf_phase value. value given in units of degree
tokens = iter1->fLine.Tokenize(" \t");
if (!tokens) {
cerr << endl << ">> PMsrHandler::HandlePlotEntry: **SEVERE ERROR** Couldn't tokenize rrf_packing in line " << iter1->fLineNo;
cerr << endl << endl;
return false;
}
if (tokens->GetEntries() != 2) {
error = true;
} else {
// get rrf packing
ostr = dynamic_cast<TObjString*>(tokens->At(1));
str = ostr->GetString();
if (str.IsDigit()) {
param.fRRFPacking = str.Atoi();
} else {
error = true;
}
}
// clean up
if (tokens) {
delete tokens;
@ -2835,6 +2755,17 @@ Bool_t PMsrHandler::HandlePlotEntry(PMsrLines &lines)
}
}
}
// check RRF entries
if (param.fRRFFreq != 0.0) {
if (param.fRRFPacking == 0) {
cerr << endl << ">> PMsrHandler::HandlePlotEntry(): **ERROR** found RRF frequency but no required RRF packing.";
cerr << endl << ">> Will ignore the RRF option.";
cerr << endl;
param.fRRFFreq = 0.0;
}
}
fPlots.push_back(param);
}
}

41
src/classes/PMusr.cpp
View File

@ -74,6 +74,35 @@ PRunData::~PRunData()
fTheory.clear();
}
//--------------------------------------------------------------------------
// SetTheoryValue
//--------------------------------------------------------------------------
/**
* <p>Sets a value of the theory vector
*
* \param i index of the theory vector
*/
void PRunData::SetTheoryValue(UInt_t i, Double_t dval)
{
if (i > fTheory.size())
fTheory.resize(i+1);
fTheory[i] = dval;
}
//--------------------------------------------------------------------------
// ReplaceTheory
//--------------------------------------------------------------------------
/**
* <p>Replaces the theory vector.
*
* \param theo vector which is replacing the current theory vector
*/
void PRunData::ReplaceTheory(const PDoubleVector &theo)
{
fTheory = theo;
}
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// implementation PNonMusrRawRunData
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
@ -504,12 +533,6 @@ PMsrRunBlock::PMsrRunBlock()
fFitRange[0] = PMUSR_UNDEFINED; // undefined start fit range
fFitRange[1] = PMUSR_UNDEFINED; // undefined end fit range
fPacking = -1; // undefined packing
fRRFFreq = PMUSR_UNDEFINED; // undefined RRF frequency
fRRFPacking = -1; // undefined RRF packing
fAlpha2ParamNo = -1; // undefined alpha2 parameter number
fBeta2ParamNo = -1; // undefined beta2 parameter number
fRightHistoNo = -1; // undefined right histogram number
fLeftHistoNo = -1; // undefined left histogram number
fXYDataIndex[0] = -1; // undefined x data index (NonMusr)
fXYDataIndex[1] = -1; // undefined y data index (NonMusr)
fXYDataLabel[0] = TString(""); // undefined x data label (NonMusr)
@ -556,12 +579,6 @@ void PMsrRunBlock::CleanUp()
fFitRange[0] = PMUSR_UNDEFINED; // undefined start fit range
fFitRange[1] = PMUSR_UNDEFINED; // undefined end fit range
fPacking = -1; // undefined packing
fRRFFreq = PMUSR_UNDEFINED; // undefined RRF frequency
fRRFPacking = -1; // undefined RRF packing
fAlpha2ParamNo = -1; // undefined alpha2 parameter number
fBeta2ParamNo = -1; // undefined beta2 parameter number
fRightHistoNo = -1; // undefined right histogram number
fLeftHistoNo = -1; // undefined left histogram number
fXYDataIndex[0] = -1; // undefined x data index (NonMusr)
fXYDataIndex[1] = -1; // undefined y data index (NonMusr)
fXYDataLabel[0] = TString("??"); // undefined x data label (NonMusr)

75
src/classes/PMusrCanvas.cpp
View File

@ -81,6 +81,9 @@ PMusrCanvas::PMusrCanvas()
fCurrentFourierPhase = fFourier.fPhaseIncrement;
fCurrentFourierPhaseText = 0;
fRRFText = 0;
fRRFLatexText = 0;
}
//--------------------------------------------------------------------------
@ -105,6 +108,9 @@ PMusrCanvas::PMusrCanvas(const Int_t number, const Char_t* title,
fCurrentFourierPhase = 0.0;
fCurrentFourierPhaseText = 0;
fRRFText = 0;
fRRFLatexText = 0;
}
//--------------------------------------------------------------------------
@ -132,6 +138,9 @@ PMusrCanvas::PMusrCanvas(const Int_t number, const Char_t* title,
fCurrentFourierPhase = 0.0;
fCurrentFourierPhaseText = 0;
fRRFText = 0;
fRRFLatexText = 0;
}
//--------------------------------------------------------------------------
@ -144,6 +153,18 @@ PMusrCanvas::~PMusrCanvas()
{
//cout << "~PMusrCanvas() called. fMainCanvas name=" << fMainCanvas->GetName() << endl;
// cleanup
if (fCurrentFourierPhaseText) {
delete fCurrentFourierPhaseText;
fCurrentFourierPhaseText = 0;
}
if (fRRFLatexText) {
delete fRRFLatexText;
fRRFLatexText = 0;
}
if (fRRFText) {
delete fRRFText;
fRRFText = 0;
}
if (fStyle) {
delete fStyle;
fStyle = 0;
@ -244,6 +265,40 @@ void PMusrCanvas::SetMsrHandler(PMsrHandler *msrHandler)
fFourier.fPlotRange[1] = fMsrHandler->GetMsrFourierList().fPlotRange[1];
}
}
// check if RRF data are present
if ((fMsrHandler->GetMsrPlotList()->at(0).fRRFPacking > 0) &&
(fMsrHandler->GetMsrPlotList()->at(0).fRRFFreq != 0.0)) {
fRRFLatexText = new TLatex();
fRRFLatexText->SetNDC(kTRUE);
fRRFLatexText->SetTextFont(62);
fRRFLatexText->SetTextSize(0.03);
fRRFText = new TString("RRF: ");
if (fMsrHandler->GetMsrPlotList()->at(0).fRRFUnit == RRF_UNIT_kHz) {
*fRRFText += TString("#nu_{RRF} = ");
*fRRFText += fMsrHandler->GetMsrPlotList()->at(0).fRRFFreq;
*fRRFText += TString(" (kHz)");
} else if (fMsrHandler->GetMsrPlotList()->at(0).fRRFUnit == RRF_UNIT_MHz) {
*fRRFText += TString("#nu_{RRF} = ");
*fRRFText += fMsrHandler->GetMsrPlotList()->at(0).fRRFFreq;
*fRRFText += TString(" (MHz)");
} else if (fMsrHandler->GetMsrPlotList()->at(0).fRRFUnit == RRF_UNIT_Mcs) {
*fRRFText += TString("#omega_{RRF} = ");
*fRRFText += fMsrHandler->GetMsrPlotList()->at(0).fRRFFreq;
*fRRFText += TString(" (Mc/s)");
} else if (fMsrHandler->GetMsrPlotList()->at(0).fRRFUnit == RRF_UNIT_G) {
*fRRFText += TString("B_{RRF} = ");
*fRRFText += fMsrHandler->GetMsrPlotList()->at(0).fRRFFreq;
*fRRFText += TString(" (G)");
} else if (fMsrHandler->GetMsrPlotList()->at(0).fRRFUnit == RRF_UNIT_T) {
*fRRFText += TString("B_{RRF} = ");
*fRRFText += fMsrHandler->GetMsrPlotList()->at(0).fRRFFreq;
*fRRFText += TString(" (T)");
}
*fRRFText += TString(", RRF packing = ");
*fRRFText += fMsrHandler->GetMsrPlotList()->at(0).fRRFPacking;
}
}
//--------------------------------------------------------------------------
@ -2531,6 +2586,11 @@ void PMusrCanvas::PlotData()
fData[i].theory->Draw("lsame");
}
}
// check if RRF and if yes show a label
if ((fRRFText != 0) && (fRRFLatexText != 0)) {
fRRFLatexText->DrawLatex(0.1, 0.92, fRRFText->Data());
}
} else { // fPlotType == MSR_PLOT_NO_MUSR
// ugly workaround since multigraphs axis are not going away when switching TMultiGraphs
delete fDataTheoryPad;
@ -2649,6 +2709,11 @@ void PMusrCanvas::PlotDifference()
for (UInt_t i=0; i<fData.size(); i++) {
fData[i].diff->Draw("pesame");
}
// check if RRF and if yes show a label
if ((fRRFText != 0) && (fRRFLatexText != 0)) {
fRRFLatexText->DrawLatex(0.1, 0.92, fRRFText->Data());
}
} else { // fPlotType == MSR_PLOT_NON_MUSR
// ugly workaround since multigraphs axis are not going away when switching TMultiGraphs
delete fDataTheoryPad;
@ -2980,6 +3045,11 @@ void PMusrCanvas::PlotFourier()
break;
}
// check if RRF and if yes show a label
if ((fRRFText != 0) && (fRRFLatexText != 0)) {
fRRFLatexText->DrawLatex(0.1, 0.92, fRRFText->Data());
}
fDataTheoryPad->Update();
fMainCanvas->cd();
@ -3249,6 +3319,11 @@ void PMusrCanvas::PlotFourierDifference()
break;
}
// check if RRF and if yes show a label
if ((fRRFText != 0) && (fRRFLatexText != 0)) {
fRRFLatexText->DrawLatex(0.1, 0.92, fRRFText->Data());
}
fDataTheoryPad->Update();
fMainCanvas->cd();

85
src/classes/PRunBase.cpp
View File

@ -145,3 +145,88 @@ void PRunBase::CleanUp()
fTheory = 0;
}
}
//--------------------------------------------------------------------------
// CalculateKaiserFilterCoeff (protected)
//--------------------------------------------------------------------------
/**
* <p> Calculates the Kaiser filter coefficients for a low pass filter with
* a cut off frequency wc.
* For details see "Zeitdiskrete Signalverarbeitung", A.V. Oppenheim, R.W. Schafer, J.R. Buck. Pearson 2004.
*
* \param wc cut off frequency
* \param A defined as \f[ A = -\log_{10}(\delta) \f], where \f[\delta\f] is the tolerance band.
* \param dw defined as \f[ \Delta\omega = \omega_{\rm S} - \omega_{\rm P} \f], where \f[ \omega_{\rm S} \f] is the
* stop band frequency, and \f[ \omega_{\rm P} \f] is the pass band frequency.
*/
void PRunBase::CalculateKaiserFilterCoeff(Double_t wc, Double_t A, Double_t dw)
{
Double_t beta;
Double_t dt = fData.GetTheoryTimeStep();
UInt_t m;
// estimate beta (see reference above, p.574ff)
if (A > 50.0) {
beta = 0.1102*(A-8.7);
} else if ((A >= 21.0) && (A <= 50.0)) {
beta = 0.5842*TMath::Power(A-21.0, 0.4) + 0.07886*(A-21.0);
} else {
beta = 0.0;
}
m = TMath::FloorNint((A-8.0)/(2.285*dw*TMath::Pi()));
// make sure m is odd
if (m % 2 == 0)
m++;
Double_t alpha = static_cast<Double_t>(m)/2.0;
Double_t dval;
Double_t dsum = 0.0;
for (UInt_t i=0; i<=m; i++) {
dval = TMath::Sin(wc*(i-alpha)*dt)/(TMath::Pi()*(i-alpha)*dt);
dval *= TMath::BesselI0(beta*TMath::Sqrt(1.0-TMath::Power((i-alpha)*dt/alpha, 2.0)))/TMath::BesselI0(beta);
dsum += dval;
fKaiserFilter.push_back(dval);
}
for (UInt_t i=0; i<=m; i++) {
fKaiserFilter[i] /= dsum;
}
/*
for (UInt_t i=0; i<=m; i++)
cout << endl << ">> " << i << ", fKaiserFilter[" << i << "]=" << fKaiserFilter[i];
cout << endl;
*/
}
//--------------------------------------------------------------------------
// FilterData (protected)
//--------------------------------------------------------------------------
/**
* <p> Filters the data with a Kaiser FIR filter.
*/
void PRunBase::FilterData()
{
PDoubleVector theoFiltered;
Double_t dval = 0.0;
const PDoubleVector *theo = fData.GetTheory();
for (UInt_t i=0; i<theo->size(); i++) {
for (UInt_t j=0; j<fKaiserFilter.size(); j++) {
if (i<j)
dval = 0.0;
else
dval += fKaiserFilter[j]*theo->at(i-j);
}
theoFiltered.push_back(dval);
dval = 0.0;
}
fData.ReplaceTheory(theoFiltered);
// shift time start by half the filter length
dval = fData.GetTheoryTimeStart() - 0.5*static_cast<Double_t>(fKaiserFilter.size())*fData.GetTheoryTimeStep();
fData.SetTheoryTimeStart(dval);
theoFiltered.clear();
}

136
src/classes/PRunSingleHisto.cpp
View File

@ -757,11 +757,15 @@ cout << endl << ">> data start time = " << fData.GetDataTimeStart();
*/
Bool_t PRunSingleHisto::PrepareViewData(PRawRunData* runData, const UInt_t histoNo)
{
// check if view_packing is wished
// check if view_packing is wished. This is a global option for all PLOT blocks!
Int_t packing = fRunInfo->GetPacking();
if (fMsrInfo->GetMsrPlotList()->at(0).fViewPacking > 0) {
packing = fMsrInfo->GetMsrPlotList()->at(0).fViewPacking;
}
// check if rrf_packing is present. This is a global option for all PLOT blocks, since operated on a single set of data.
if (fMsrInfo->GetMsrPlotList()->at(0).fRRFPacking > 0) {
packing = fMsrInfo->GetMsrPlotList()->at(0).fRRFPacking;
}
// transform raw histo data. This is done the following way (for details see the manual):
// for the single histo fit, just the rebinned raw data are copied
@ -855,9 +859,10 @@ Bool_t PRunSingleHisto::PrepareViewData(PRawRunData* runData, const UInt_t histo
bkg = par[fRunInfo->GetBkgFitParamNo()-1];
}
Double_t value = 0.0;
Double_t expval;
Double_t time;
Double_t value = 0.0;
Double_t expval = 0.0;
Double_t rrf_val = 0.0;
Double_t time = 0.0;
// data start at data_start-t0 shifted by (pack-1)/2
fData.SetDataTimeStart(fTimeResolution*((Double_t)start-(Double_t)t0+(Double_t)(packing-1)/2.0));
@ -870,21 +875,67 @@ cout << endl << "--------------------------------" << endl;
*/
Double_t normalizer = 1.0;
for (Int_t i=start; i<end; i++) {
if (((i-start) % packing == 0) && (i != start)) { // fill data
// in order that after rebinning the fit does not need to be redone (important for plots)
// the value is normalize to per 1 nsec
normalizer = packing * (fTimeResolution * 1.0e3); // fTimeResolution us->ns
value /= normalizer;
time = (((Double_t)i-(Double_t)(packing-1)/2.0)-t0)*fTimeResolution;
expval = TMath::Exp(+time/tau)/N0;
fData.AppendValue(-1.0+expval*(value-bkg));
//cout << endl << ">> i=" << i << ",t0=" << t0 << ",time=" << time << ",expval=" << expval << ",value=" << value << ",bkg=" << bkg << ",expval*(value-bkg)-1=" << expval*(value-bkg)-1.0;
fData.AppendErrorValue(expval*TMath::Sqrt(value/normalizer));
//cout << endl << ">> " << time << ", " << expval << ", " << -1.0+expval*(value-bkg) << ", " << expval*TMath::Sqrt(value/packing);
value = 0.0;
Double_t gammaRRF = 0.0, wRRF = 0.0, phaseRRF = 0.0;
if (fMsrInfo->GetMsrPlotList()->at(0).fRRFFreq == 0.0) { // normal Data representation
for (Int_t i=start; i<end; i++) {
if (((i-start) % packing == 0) && (i != start)) { // fill data
// in order that after rebinning the fit does not need to be redone (important for plots)
// the value is normalize to per 1 nsec
normalizer = packing * (fTimeResolution * 1.0e3); // fTimeResolution us->ns
value /= normalizer;
time = (((Double_t)i-(Double_t)(packing-1)/2.0)-t0)*fTimeResolution;
expval = TMath::Exp(+time/tau)/N0;
fData.AppendValue(-1.0+expval*(value-bkg));
//cout << endl << ">> i=" << i << ",t0=" << t0 << ",time=" << time << ",expval=" << expval << ",value=" << value << ",bkg=" << bkg << ",expval*(value-bkg)-1=" << expval*(value-bkg)-1.0;
fData.AppendErrorValue(expval*TMath::Sqrt(value/normalizer));
//cout << endl << ">> " << time << ", " << expval << ", " << -1.0+expval*(value-bkg) << ", " << expval*TMath::Sqrt(value/packing);
value = 0.0;
}
value += runData->GetDataBin(histoNo)->at(i);
}
} else { // RRF representation
// check which units shall be used
switch (fMsrInfo->GetMsrPlotList()->at(0).fRRFUnit) {
case RRF_UNIT_kHz:
gammaRRF = TMath::TwoPi()*1.0e-3;
break;
case RRF_UNIT_MHz:
gammaRRF = TMath::TwoPi();
break;
case RRF_UNIT_Mcs:
gammaRRF = 1.0;
break;
case RRF_UNIT_G:
gammaRRF = 0.0135538817*TMath::TwoPi();
break;
case RRF_UNIT_T:
gammaRRF = 0.0135538817*TMath::TwoPi()*1.0e4;
break;
default:
gammaRRF = TMath::TwoPi();
break;
}
wRRF = gammaRRF * fMsrInfo->GetMsrPlotList()->at(0).fRRFFreq;
phaseRRF = fMsrInfo->GetMsrPlotList()->at(0).fRRFPhase / 180.0 * TMath::Pi();
// in order that after rebinning the fit does not need to be redone (important for plots)
// the value is normalize to per 1 nsec
normalizer = (fTimeResolution * 1.0e3); // fTimeResolution us->ns
Double_t error = 0.0;
for (Int_t i=start; i<end; i++) {
if (((i-start) % packing == 0) && (i != start)) { // fill data
fData.AppendValue(2.0*value/packing); // factor 2 needed because cos(a)cos(b) = 1/2(cos(a+b)+cos(a-b))
fData.AppendErrorValue(expval*TMath::Sqrt(error)/normalizer/packing);
value = 0.0;
error = 0.0;
}
time = ((Double_t)i-t0)*fTimeResolution;
expval = TMath::Exp(+time/tau)/N0;
rrf_val = (-1.0+expval*(runData->GetDataBin(histoNo)->at(i)/normalizer-bkg))*TMath::Cos(wRRF * time + phaseRRF);
value += rrf_val;
error += runData->GetDataBin(histoNo)->at(i);
}
value += runData->GetDataBin(histoNo)->at(i);
}
// count the number of bins to be fitted
@ -904,23 +955,58 @@ cout << endl << "--------------------------------" << endl;
Double_t theoryValue;
UInt_t size = runData->GetDataBin(histoNo)->size();
Double_t factor = 1.0;
if (fData.GetValue()->size() * 10 > runData->GetDataBin(histoNo)->size()) {
size = fData.GetValue()->size() * 10;
factor = (Double_t)runData->GetDataBin(histoNo)->size() / (Double_t)size;
UInt_t rebinRRF = 0;
if (wRRF == 0) { // no RRF
// check if a finer binning for the theory is needed
if (fData.GetValue()->size() * 10 > runData->GetDataBin(histoNo)->size()) {
size = fData.GetValue()->size() * 10;
factor = (Double_t)runData->GetDataBin(histoNo)->size() / (Double_t)size;
}
fData.SetTheoryTimeStart(fData.GetDataTimeStart());
fData.SetTheoryTimeStep(fTimeResolution*factor);
} else { // RRF
rebinRRF = static_cast<UInt_t>((TMath::Pi()/2.0/wRRF)/fTimeResolution); // RRF time resolution / data time resolution
fData.SetTheoryTimeStart(fData.GetDataTimeStart());
fData.SetTheoryTimeStep(TMath::Pi()/2.0/wRRF/rebinRRF); // = theory time resolution as close as possible to the data time resolution compatible with wRRF
}
//cout << endl << ">> runData->GetDataBin(histoNo).size() = " << runData->GetDataBin(histoNo).size() << ", fData.GetValue()->size() * 10 = " << fData.GetValue()->size() * 10 << ", size = " << size << ", factor = " << factor << endl;
fData.SetTheoryTimeStart(fData.GetDataTimeStart());
fData.SetTheoryTimeStep(fTimeResolution*factor);
//cout << endl << ">> size=" << size << ", startTime=" << startTime << ", fTimeResolution=" << fTimeResolution;
for (UInt_t i=0; i<size; i++) {
time = fData.GetTheoryTimeStart() + (Double_t)i*fData.GetTheoryTimeStep();
theoryValue = fTheory->Func(time, par, fFuncValues);
if (wRRF != 0.0) {
theoryValue *= 2.0*TMath::Cos(wRRF * time + phaseRRF);
}
if (fabs(theoryValue) > 10.0) { // dirty hack needs to be fixed!!
theoryValue = 0.0;
}
fData.AppendTheoryValue(theoryValue);
}
// if RRF filter the theory with a FIR Kaiser low pass filter
if (wRRF != 0.0) {
// rebin theory to the RRF frequency
if (rebinRRF != 0) {
Double_t dval = 0.0;
PDoubleVector theo;
for (UInt_t i=0; i<fData.GetTheory()->size(); i++) {
if ((i % rebinRRF == 0) && (i != 0)) {
theo.push_back(dval/rebinRRF);
dval = 0.0;
}
dval += fData.GetTheory()->at(i);
}
fData.SetTheoryTimeStart(fData.GetTheoryTimeStart()+static_cast<Double_t>(rebinRRF-1)*fData.GetTheoryTimeStep()/2.0);
fData.SetTheoryTimeStep(rebinRRF*fData.GetTheoryTimeStep());
fData.ReplaceTheory(theo);
theo.clear();
}
// filter data
CalculateKaiserFilterCoeff(wRRF, 60.0, 0.2); // w_c = wRRF, A = -20 log_10(delta), Delta w / w_c = (w_s - w_p) / (2 w_c)
FilterData();
}
// clean up
par.clear();

33
src/include/PMusr.h
View File

@ -116,6 +116,14 @@ using namespace std;
#define FOURIER_PLOT_POWER 4
#define FOURIER_PLOT_PHASE 5
//-------------------------------------------------------------
// RRF related tags
#define RRF_UNIT_kHz 0
#define RRF_UNIT_MHz 1
#define RRF_UNIT_Mcs 2
#define RRF_UNIT_G 3
#define RRF_UNIT_T 4
//-------------------------------------------------------------
/**
* <p> typedef to make to code more readable.
@ -201,6 +209,9 @@ class PRunData {
virtual void AppendXTheoryValue(Double_t dval) { fXTheory.push_back(dval); }
virtual void AppendTheoryValue(Double_t dval) { fTheory.push_back(dval); }
virtual void SetTheoryValue(UInt_t i, Double_t dval);
virtual void ReplaceTheory(const PDoubleVector &theo);
private:
// data related info
Double_t fDataTimeStart;
@ -403,12 +414,6 @@ class PMsrRunBlock {
virtual Int_t GetT0(UInt_t i=0);
virtual Double_t GetFitRange(UInt_t i);
virtual Int_t GetPacking() { return fPacking; }
virtual Double_t GetRRFFreq() { return fRRFFreq; }
virtual Int_t GetRRFPacking() { return fRRFPacking; }
virtual Int_t GetAlpha2ParamNo() { return fAlpha2ParamNo; }
virtual Int_t GetBeta2ParamNo() { return fBeta2ParamNo; }
virtual Int_t GetRightHistoNo() { return fRightHistoNo; }
virtual Int_t GetLeftHistoNo() { return fLeftHistoNo; }
virtual Int_t GetXDataIndex() { return fXYDataIndex[0]; }
virtual Int_t GetYDataIndex() { return fXYDataIndex[1]; }
virtual TString* GetXDataLabel() { return &fXYDataLabel[0]; }
@ -444,12 +449,6 @@ class PMsrRunBlock {
virtual void SetT0(Int_t ival, UInt_t idx);
virtual void SetFitRange(Double_t dval, UInt_t idx);
virtual void SetPacking(Int_t ival) { fPacking = ival; }
virtual void SetRRFFreq(Double_t dval) { fRRFFreq = dval; }
virtual void SetRRFPacking(Int_t ival) { fRRFPacking = ival; }
virtual void SetAlpha2ParamNo(Int_t ival) { fAlpha2ParamNo = ival; }
virtual void SetBeta2ParamNo(Int_t ival) { fBeta2ParamNo = ival; }
virtual void SetRightHistoNo(Int_t ival) { fRightHistoNo = ival; }
virtual void SetLeftHistoNo(Int_t ival) { fLeftHistoNo = ival; }
virtual void SetXDataIndex(Int_t ival) { fXYDataIndex[0] = ival; }
virtual void SetYDataIndex(Int_t ival) { fXYDataIndex[1] = ival; }
virtual void SetXDataLabel(TString& str) { fXYDataLabel[0] = str; }
@ -477,12 +476,6 @@ class PMsrRunBlock {
PIntVector fT0; ///< t0 bins (fit type 0, 2, 4). if fit type 0 -> f0, f1, f2, ...; if fit type 2, 4 -> f0, b0, f1, b1, ...
Double_t fFitRange[2]; ///< fit range in (us)
Int_t fPacking; ///< packing/rebinning
Double_t fRRFFreq; ///< rotating reference frequency (fit type 4)
Int_t fRRFPacking; ///< rotating reference packing (fit type 4)
Int_t fAlpha2ParamNo; ///< rotating reference alpha2 (fit type 4)
Int_t fBeta2ParamNo; ///< rotating reference beta2 (fit type 4)
Int_t fRightHistoNo; ///< rotating reference right histogram number (fit type 4)
Int_t fLeftHistoNo; ///< rotating reference left histogram number (fit type 4)
Int_t fXYDataIndex[2]; ///< used to get the data indices when using db-files (fit type 8)
TString fXYDataLabel[2]; ///< used to get the indices via labels when using db-files (fit type 8)
};
@ -525,6 +518,10 @@ typedef struct {
PDoubleVector fTmax; ///< time maximum
PDoubleVector fYmin; ///< asymmetry/counts minimum
PDoubleVector fYmax; ///< asymmetry/counts maximum
UInt_t fRRFPacking; ///< rotating reference frame (RRF) packing
Double_t fRRFFreq; ///< RRF frequency
UInt_t fRRFUnit; ///< RRF frequency unit. 0=kHz, 1=MHz, 2=Mc/s, 3=Gauss, 4=Tesla
Double_t fRRFPhase; ///< RRF phase
} PMsrPlotStructure;
//-------------------------------------------------------------

2
src/include/PMusrCanvas.h
View File

@ -204,6 +204,8 @@ class PMusrCanvas : public TObject, public TQObject
Double_t fCurrentFourierPhase; /// holds the current Fourier phase
TLatex *fCurrentFourierPhaseText; /// used in Re/Im Fourier to show the current phase in the pad
TString *fRRFText;
TLatex *fRRFLatexText; /// used to display RRF info
TStyle *fStyle;

9
src/include/PRunBase.h
View File

@ -79,10 +79,15 @@ class PRunBase
Double_t fTimeResolution; ///< time resolution in (us)
PIntVector fT0s; ///< all t0's of a run! The derived classes will handle it
virtual Bool_t PrepareData() = 0; // pure virtual, i.e. needs to be implemented by the deriving class!!
PDoubleVector fFuncValues; ///< is keeping the values of the functions from the FUNCTIONS block
PTheory *fTheory; ///< theory needed to calculate chi-square
PDoubleVector fKaiserFilter;
virtual Bool_t PrepareData() = 0; // pure virtual, i.e. needs to be implemented by the deriving class!!
virtual void CalculateKaiserFilterCoeff(Double_t wc, Double_t A, Double_t dw);
virtual void FilterData();
};
#endif // _PRUNBASE_H_

2
src/musrt0.cpp
View File

@ -265,6 +265,7 @@ int main(int argc, char *argv[])
}
break;
case MSR_FITTYPE_ASYM_RRF:
/*
for (unsigned int j=0; j<runList->at(i).GetRunNameSize(); j++) { // necessary in case of ADDRUN
if (!musrt0_item(app, msrHandler, dataHandler->GetRunData(*(runList->at(i).GetRunName(j))), i, runList->at(i).GetForwardHistoNo(), 0, j)) {
musrt0_cleanup(saxParser, startupHandler, msrHandler, dataHandler);
@ -283,6 +284,7 @@ int main(int argc, char *argv[])
exit(0);
}
}
*/
break;
default:
break;

18
src/tests/t0NotEqFirstGoodData/t0NotEqFirstGoodData.C
View File

@ -68,17 +68,20 @@ void t0NotEqFirstGoodData()
// create histos
UInt_t t0[4] = {3419, 3419, 3419, 3419};
UInt_t n0[4] = {200.0, 205.0, 203.0, 198.0};
UInt_t n0[4] = {2000.0, 2005.0, 2003.0, 1998.0};
UInt_t bkg[4] = {11.0, 11.0, 5.0, 8.0};
const Double_t tau = 2197.019; // ns
// asymmetry related stuff
const Double_t timeResolution = 0.1953125; // ns
Double_t a[4] = {0.16, 0.161, 0.16, 0.159};
Double_t phase[4] = {13.0*TMath::Pi()/180.0, 103.0*TMath::Pi()/180.0, 193.0*TMath::Pi()/180.0, 283.0*TMath::Pi()/180.0};
Double_t a0[4] = {0.16, 0.16, 0.16, 0.16};
Double_t a1[4] = {0.0, 0.0, 0.0, 0.0};
Double_t phase[4] = {0.0*TMath::Pi()/180.0, 90.0*TMath::Pi()/180.0, 180.0*TMath::Pi()/180.0, 270.0*TMath::Pi()/180.0};
const Double_t gamma = 0.00001355; // gamma/(2pi)
Double_t bb = 15000.0; // field in Gauss
Double_t bb0 = 15000.0; // field in Gauss
Double_t bb1 = 15702.0; // field in Gauss
// Double_t bb = 100.0; // field in Gauss
Double_t sigma = 1.05/1000.0; // Gaussian sigma in 1/ns
Double_t sigma0 = 0.05/1000.0; // Gaussian sigma in 1/ns
Double_t sigma1 = 0.10/1000.0; // Gaussian sigma in 1/ns
TH1F *histo[8];
char str[128];
@ -96,7 +99,8 @@ void t0NotEqFirstGoodData()
histo[i]->SetBinContent(j, bkg[i]);
} else {
time = (Double_t)(j-t0[i])*timeResolution;
dval = (Double_t)n0[i]*TMath::Exp(-time/tau)*(1.0+a[i]*TMath::Exp(-0.5*TMath::Power(time*sigma,2.0))*TMath::Cos(TMath::TwoPi()*gamma*bb*time+phase[i]))+(Double_t)bkg[i];
dval = (Double_t)n0[i]*TMath::Exp(-time/tau)*(1.0+a0[i]*TMath::Exp(-0.5*TMath::Power(time*sigma0,2.0))*TMath::Cos(TMath::TwoPi()*gamma*bb0*time+phase[i])
+a1[i]*TMath::Exp(-0.5*TMath::Power(time*sigma1,2.0))*TMath::Cos(TMath::TwoPi()*gamma*bb1*time+phase[i]))+(Double_t)bkg[i];
histo[i]->SetBinContent(j, (UInt_t)dval);
}
}
@ -134,7 +138,7 @@ void t0NotEqFirstGoodData()
decayAnaModule->Add(histo[i]);
// write file
TFile *fout = new TFile("010105.root", "RECREATE", "Midas Fake Histograms");
TFile *fout = new TFile("010108.root", "RECREATE", "Midas Fake Histograms");
if (fout == 0) {
cout << endl << "**ERROR** Couldn't create ROOT file";
cout << endl << endl;