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

Merge branch 'master' into root6

Browse Source
This commit is contained in:
suter_a 2016-02-23 17:37:52 +01:00
commit f3542337fc
7 changed files with 1108 additions and 1118 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
ChangeLog
View File

@ -4,7 +4,7 @@
changes since 0.17.0 changes since 0.17.0
=================================== ===================================
NEW 2016-02-23 It is now possible to export the averaged data/Fourier
changes since 0.16.0 changes since 0.16.0
=================================== ===================================

754
src/classes/PMusrCanvas.cpp
View File

@ -1516,7 +1516,6 @@ void PMusrCanvas::ExportData(const Char_t *fileName)
Int_t xmaxBin; Int_t xmaxBin;
Double_t xmin; Double_t xmin;
Double_t xmax; Double_t xmax;
Double_t time, freq;
Double_t xval, yval; Double_t xval, yval;
switch (fPlotType) { switch (fPlotType) {
@ -1533,29 +1532,12 @@ void PMusrCanvas::ExportData(const Char_t *fileName)
xmax = fHistoFrame->GetXaxis()->GetBinCenter(xmaxBin); xmax = fHistoFrame->GetXaxis()->GetBinCenter(xmaxBin);
// fill ascii dump data // fill ascii dump data
for (UInt_t i=0; i<fData.size(); i++) { // go through all the histogramms if (fAveragedView) {
// clean up dump GetExportDataSet(fDataAvg.diff, xmin, xmax, dumpVector);
dump.dataX.clear(); } else { // go through all the histogramms
dump.data.clear(); for (UInt_t i=0; i<fData.size(); i++) { // go through all the histogramms
dump.dataErr.clear(); GetExportDataSet(fData[i].diff, xmin, xmax, dumpVector);
dump.theoryX.clear();
dump.theory.clear();
// go through all difference data bins
for (Int_t j=1; j<fData[i].diff->GetNbinsX(); j++) {
// get time
time = fData[i].diff->GetBinCenter(j);
// check if time is in the current range
if ((time >= xmin) && (time <= xmax)) {
dump.dataX.push_back(time);
dump.data.push_back(fData[i].diff->GetBinContent(j));
dump.dataErr.push_back(fData[i].diff->GetBinError(j));
}
} }
// if anything found keep it
if (dump.dataX.size() > 0)
dumpVector.push_back(dump);
} }
break; break;
case PV_FOURIER_REAL: case PV_FOURIER_REAL:
@ -1566,28 +1548,12 @@ void PMusrCanvas::ExportData(const Char_t *fileName)
xmax = fData[0].diffFourierRe->GetXaxis()->GetBinCenter(xmaxBin); xmax = fData[0].diffFourierRe->GetXaxis()->GetBinCenter(xmaxBin);
// fill ascii dump data // fill ascii dump data
for (UInt_t i=0; i<fData.size(); i++) { // go through all the histogramms if (fAveragedView) {
// clean up dump GetExportDataSet(fDataAvg.diffFourierRe, xmin, xmax, dumpVector, false);
dump.dataX.clear(); } else { // go through all the histogramms
dump.data.clear(); for (UInt_t i=0; i<fData.size(); i++) { // go through all the histogramms
dump.dataErr.clear(); GetExportDataSet(fData[i].diffFourierRe, xmin, xmax, dumpVector, false);
dump.theoryX.clear();
dump.theory.clear();
// go through all data bins
for (Int_t j=1; j<fData[i].diffFourierRe->GetNbinsX(); j++) {
// get frequency
freq = fData[i].diffFourierRe->GetBinCenter(j);
// check if time is in the current range
if ((freq >= xmin) && (freq <= xmax)) {
dump.dataX.push_back(freq);
dump.data.push_back(fData[i].diffFourierRe->GetBinContent(j));
}
} }
// if anything found keep it
if (dump.dataX.size() > 0)
dumpVector.push_back(dump);
} }
break; break;
case PV_FOURIER_IMAG: case PV_FOURIER_IMAG:
@ -1598,28 +1564,12 @@ void PMusrCanvas::ExportData(const Char_t *fileName)
xmax = fData[0].diffFourierIm->GetXaxis()->GetBinCenter(xmaxBin); xmax = fData[0].diffFourierIm->GetXaxis()->GetBinCenter(xmaxBin);
// fill ascii dump data // fill ascii dump data
for (UInt_t i=0; i<fData.size(); i++) { // go through all the histogramms if (fAveragedView) {
// clean up dump GetExportDataSet(fDataAvg.diffFourierIm, xmin, xmax, dumpVector, false);
dump.dataX.clear(); } else { // go through all the histogramms
dump.data.clear(); for (UInt_t i=0; i<fData.size(); i++) { // go through all the histogramms
dump.dataErr.clear(); GetExportDataSet(fData[i].diffFourierIm, xmin, xmax, dumpVector, false);
dump.theoryX.clear();
dump.theory.clear();
// go through all data bins
for (Int_t j=1; j<fData[i].diffFourierIm->GetNbinsX(); j++) {
// get frequency
freq = fData[i].diffFourierIm->GetBinCenter(j);
// check if time is in the current range
if ((freq >= xmin) && (freq <= xmax)) {
dump.dataX.push_back(freq);
dump.data.push_back(fData[i].diffFourierIm->GetBinContent(j));
}
} }
// if anything found keep it
if (dump.dataX.size() > 0)
dumpVector.push_back(dump);
} }
break; break;
case PV_FOURIER_REAL_AND_IMAG: case PV_FOURIER_REAL_AND_IMAG:
@ -1630,37 +1580,14 @@ void PMusrCanvas::ExportData(const Char_t *fileName)
xmax = fData[0].diffFourierRe->GetXaxis()->GetBinCenter(xmaxBin); xmax = fData[0].diffFourierRe->GetXaxis()->GetBinCenter(xmaxBin);
// fill ascii dump data // fill ascii dump data
for (UInt_t i=0; i<fData.size(); i++) { // go through all the histogramms if (fAveragedView) {
// clean up dump GetExportDataSet(fDataAvg.diffFourierRe, xmin, xmax, dumpVector, false);
dump.dataX.clear(); GetExportDataSet(fDataAvg.diffFourierIm, xmin, xmax, dumpVector, false);
dump.data.clear(); } else { // go through all the histogramms
dump.dataErr.clear(); for (UInt_t i=0; i<fData.size(); i++) { // go through all the histogramms
dump.theoryX.clear(); GetExportDataSet(fData[i].diffFourierRe, xmin, xmax, dumpVector, false);
dump.theory.clear(); GetExportDataSet(fData[i].diffFourierIm, xmin, xmax, dumpVector, false);
// go through all data bins
for (Int_t j=1; j<fData[i].diffFourierRe->GetNbinsX(); j++) {
// get frequency
freq = fData[i].diffFourierRe->GetBinCenter(j);
// check if time is in the current range
if ((freq >= xmin) && (freq <= xmax)) {
dump.dataX.push_back(freq);
dump.data.push_back(fData[i].diffFourierRe->GetBinContent(j));
}
} }
for (Int_t j=1; j<fData[i].diffFourierIm->GetNbinsX(); j++) {
// get frequency
freq = fData[i].diffFourierIm->GetBinCenter(j);
// check if time is in the current range
if ((freq >= xmin) && (freq <= xmax)) {
dump.dataX.push_back(freq);
dump.data.push_back(fData[i].diffFourierIm->GetBinContent(j));
}
}
// if anything found keep it
if (dump.dataX.size() > 0)
dumpVector.push_back(dump);
} }
break; break;
case PV_FOURIER_PWR: case PV_FOURIER_PWR:
@ -1671,28 +1598,12 @@ void PMusrCanvas::ExportData(const Char_t *fileName)
xmax = fData[0].diffFourierPwr->GetXaxis()->GetBinCenter(xmaxBin); xmax = fData[0].diffFourierPwr->GetXaxis()->GetBinCenter(xmaxBin);
// fill ascii dump data // fill ascii dump data
for (UInt_t i=0; i<fData.size(); i++) { // go through all the histogramms if (fAveragedView) {
// clean up dump GetExportDataSet(fDataAvg.diffFourierPwr, xmin, xmax, dumpVector, false);
dump.dataX.clear(); } else { // go through all the histogramms
dump.data.clear(); for (UInt_t i=0; i<fData.size(); i++) {
dump.dataErr.clear(); GetExportDataSet(fData[i].diffFourierPwr, xmin, xmax, dumpVector, false);
dump.theoryX.clear();
dump.theory.clear();
// go through all data bins
for (Int_t j=1; j<fData[i].diffFourierPwr->GetNbinsX(); j++) {
// get frequency
freq = fData[i].diffFourierPwr->GetBinCenter(j);
// check if time is in the current range
if ((freq >= xmin) && (freq <= xmax)) {
dump.dataX.push_back(freq);
dump.data.push_back(fData[i].diffFourierPwr->GetBinContent(j));
}
} }
// if anything found keep it
if (dump.dataX.size() > 0)
dumpVector.push_back(dump);
} }
break; break;
case PV_FOURIER_PHASE: case PV_FOURIER_PHASE:
@ -1703,28 +1614,12 @@ void PMusrCanvas::ExportData(const Char_t *fileName)
xmax = fData[0].diffFourierPhase->GetXaxis()->GetBinCenter(xmaxBin); xmax = fData[0].diffFourierPhase->GetXaxis()->GetBinCenter(xmaxBin);
// fill ascii dump data // fill ascii dump data
for (UInt_t i=0; i<fData.size(); i++) { // go through all the histogramms if (fAveragedView) {
// clean up dump GetExportDataSet(fDataAvg.diffFourierPhase, xmin, xmax, dumpVector, false);
dump.dataX.clear(); } else { // go through all the histogramms
dump.data.clear(); for (UInt_t i=0; i<fData.size(); i++) {
dump.dataErr.clear(); GetExportDataSet(fData[i].diffFourierPhase, xmin, xmax, dumpVector, false);
dump.theoryX.clear();
dump.theory.clear();
// go through all data bins
for (Int_t j=1; j<fData[i].diffFourierPhase->GetNbinsX(); j++) {
// get frequency
freq = fData[i].diffFourierPhase->GetBinCenter(j);
// check if time is in the current range
if ((freq >= xmin) && (freq <= xmax)) {
dump.dataX.push_back(freq);
dump.data.push_back(fData[i].diffFourierPhase->GetBinContent(j));
}
} }
// if anything found keep it
if (dump.dataX.size() > 0)
dumpVector.push_back(dump);
} }
break; break;
default: default:
@ -1740,40 +1635,14 @@ void PMusrCanvas::ExportData(const Char_t *fileName)
xmax = fHistoFrame->GetXaxis()->GetBinCenter(xmaxBin); xmax = fHistoFrame->GetXaxis()->GetBinCenter(xmaxBin);
// fill ascii dump data // fill ascii dump data
for (UInt_t i=0; i<fData.size(); i++) { // go through all the histogramms if (fAveragedView) {
// clean up dump GetExportDataSet(fDataAvg.data, xmin, xmax, dumpVector);
dump.dataX.clear(); GetExportDataSet(fDataAvg.theory, xmin, xmax, dumpVector, false);
dump.data.clear(); } else { // go through all the histogramms
dump.dataErr.clear(); for (UInt_t i=0; i<fData.size(); i++) {
dump.theoryX.clear(); GetExportDataSet(fData[i].data, xmin, xmax, dumpVector);
dump.theory.clear(); GetExportDataSet(fData[i].theory, xmin, xmax, dumpVector, false);
// go through all data bins
for (Int_t j=1; j<fData[i].data->GetNbinsX(); j++) {
// get time
time = fData[i].data->GetBinCenter(j);
// check if time is in the current range
if ((time >= xmin) && (time <= xmax)) {
dump.dataX.push_back(time);
dump.data.push_back(fData[i].data->GetBinContent(j));
dump.dataErr.push_back(fData[i].data->GetBinError(j));
}
} }
// go through all theory bins
for (Int_t j=1; j<fData[i].theory->GetNbinsX(); j++) {
// get time
time = fData[i].theory->GetBinCenter(j);
// check if time is in the current range
if ((time >= xmin) && (time <= xmax)) {
dump.theoryX.push_back(time);
dump.theory.push_back(fData[i].theory->GetBinContent(j));
}
}
// if anything found keep it
if (dump.dataX.size() > 0)
dumpVector.push_back(dump);
} }
break; break;
@ -1785,39 +1654,14 @@ void PMusrCanvas::ExportData(const Char_t *fileName)
xmax = fData[0].dataFourierRe->GetXaxis()->GetBinCenter(xmaxBin); xmax = fData[0].dataFourierRe->GetXaxis()->GetBinCenter(xmaxBin);
// fill ascii dump data // fill ascii dump data
for (UInt_t i=0; i<fData.size(); i++) { // go through all the histogramms if (fAveragedView) {
// clean up dump GetExportDataSet(fDataAvg.dataFourierRe, xmin, xmax, dumpVector, false);
dump.dataX.clear(); GetExportDataSet(fDataAvg.theoryFourierRe, xmin, xmax, dumpVector, false);
dump.data.clear(); } else { // go through all the histogramms
dump.dataErr.clear(); for (UInt_t i=0; i<fData.size(); i++) {
dump.theoryX.clear(); GetExportDataSet(fData[i].dataFourierRe, xmin, xmax, dumpVector, false);
dump.theory.clear(); GetExportDataSet(fData[i].theoryFourierRe, xmin, xmax, dumpVector, false);
// go through all data bins
for (Int_t j=1; j<fData[i].dataFourierRe->GetNbinsX(); j++) {
// get frequency
freq = fData[i].dataFourierRe->GetBinCenter(j);
// check if time is in the current range
if ((freq >= xmin) && (freq <= xmax)) {
dump.dataX.push_back(freq);
dump.data.push_back(fData[i].dataFourierRe->GetBinContent(j));
}
} }
// go through all theory bins
for (Int_t j=1; j<fData[i].theoryFourierRe->GetNbinsX(); j++) {
// get frequency
freq = fData[i].theoryFourierRe->GetBinCenter(j);
// check if time is in the current range
if ((freq >= xmin) && (freq <= xmax)) {
dump.theoryX.push_back(freq);
dump.theory.push_back(fData[i].theoryFourierRe->GetBinContent(j));
}
}
// if anything found keep it
if (dump.dataX.size() > 0)
dumpVector.push_back(dump);
} }
break; break;
case PV_FOURIER_IMAG: case PV_FOURIER_IMAG:
@ -1828,39 +1672,14 @@ void PMusrCanvas::ExportData(const Char_t *fileName)
xmax = fData[0].dataFourierIm->GetXaxis()->GetBinCenter(xmaxBin); xmax = fData[0].dataFourierIm->GetXaxis()->GetBinCenter(xmaxBin);
// fill ascii dump data // fill ascii dump data
for (UInt_t i=0; i<fData.size(); i++) { // go through all the histogramms if (fAveragedView) {
// clean up dump GetExportDataSet(fDataAvg.dataFourierIm, xmin, xmax, dumpVector, false);
dump.dataX.clear(); GetExportDataSet(fDataAvg.theoryFourierIm, xmin, xmax, dumpVector, false);
dump.data.clear(); } else { // go through all the histogramms
dump.dataErr.clear(); for (UInt_t i=0; i<fData.size(); i++) {
dump.theoryX.clear(); GetExportDataSet(fData[i].dataFourierIm, xmin, xmax, dumpVector, false);
dump.theory.clear(); GetExportDataSet(fData[i].theoryFourierIm, xmin, xmax, dumpVector, false);
// go through all data bins
for (Int_t j=1; j<fData[i].dataFourierIm->GetNbinsX(); j++) {
// get frequency
freq = fData[i].dataFourierIm->GetBinCenter(j);
// check if time is in the current range
if ((freq >= xmin) && (freq <= xmax)) {
dump.dataX.push_back(freq);
dump.data.push_back(fData[i].dataFourierIm->GetBinContent(j));
}
} }
// go through all theory bins
for (Int_t j=1; j<fData[i].theoryFourierIm->GetNbinsX(); j++) {
// get frequency
freq = fData[i].theoryFourierIm->GetBinCenter(j);
// check if time is in the current range
if ((freq >= xmin) && (freq <= xmax)) {
dump.theoryX.push_back(freq);
dump.theory.push_back(fData[i].theoryFourierIm->GetBinContent(j));
}
}
// if anything found keep it
if (dump.dataX.size() > 0)
dumpVector.push_back(dump);
} }
break; break;
case PV_FOURIER_REAL_AND_IMAG: case PV_FOURIER_REAL_AND_IMAG:
@ -1870,80 +1689,18 @@ void PMusrCanvas::ExportData(const Char_t *fileName)
xmin = fData[0].dataFourierRe->GetXaxis()->GetBinCenter(xminBin); xmin = fData[0].dataFourierRe->GetXaxis()->GetBinCenter(xminBin);
xmax = fData[0].dataFourierRe->GetXaxis()->GetBinCenter(xmaxBin); xmax = fData[0].dataFourierRe->GetXaxis()->GetBinCenter(xmaxBin);
// fill ascii dump data if (fAveragedView) {
for (UInt_t i=0; i<fData.size(); i++) { // go through all the histogramms GetExportDataSet(fDataAvg.dataFourierRe, xmin, xmax, dumpVector, false);
//----------------------------- GetExportDataSet(fDataAvg.theoryFourierRe, xmin, xmax, dumpVector, false);
// Re GetExportDataSet(fDataAvg.dataFourierIm, xmin, xmax, dumpVector, false);
//----------------------------- GetExportDataSet(fDataAvg.theoryFourierIm, xmin, xmax, dumpVector, false);
// clean up dump } else { // go through all the histogramms
dump.dataX.clear(); for (UInt_t i=0; i<fData.size(); i++) {
dump.data.clear(); GetExportDataSet(fData[i].dataFourierRe, xmin, xmax, dumpVector, false);
dump.dataErr.clear(); GetExportDataSet(fData[i].theoryFourierRe, xmin, xmax, dumpVector, false);
dump.theoryX.clear(); GetExportDataSet(fData[i].dataFourierIm, xmin, xmax, dumpVector, false);
dump.theory.clear(); GetExportDataSet(fData[i].theoryFourierIm, xmin, xmax, dumpVector, false);
// go through all data bins
for (Int_t j=1; j<fData[i].dataFourierRe->GetNbinsX(); j++) {
// get frequency
freq = fData[i].dataFourierRe->GetBinCenter(j);
// check if time is in the current range
if ((freq >= xmin) && (freq <= xmax)) {
dump.dataX.push_back(freq);
dump.data.push_back(fData[i].dataFourierRe->GetBinContent(j));
}
} }
// go through all theory bins
for (Int_t j=1; j<fData[i].theoryFourierRe->GetNbinsX(); j++) {
// get frequency
freq = fData[i].theoryFourierRe->GetBinCenter(j);
// check if time is in the current range
if ((freq >= xmin) && (freq <= xmax)) {
dump.theoryX.push_back(freq);
dump.theory.push_back(fData[i].theoryFourierRe->GetBinContent(j));
}
}
// if anything found keep it
if (dump.dataX.size() > 0)
dumpVector.push_back(dump);
//-----------------------------
// Im
//-----------------------------
// clean up dump
dump.dataX.clear();
dump.data.clear();
dump.dataErr.clear();
dump.theoryX.clear();
dump.theory.clear();
// go through all data bins
for (Int_t j=1; j<fData[i].dataFourierIm->GetNbinsX(); j++) {
// get frequency
freq = fData[i].dataFourierIm->GetBinCenter(j);
// check if time is in the current range
if ((freq >= xmin) && (freq <= xmax)) {
dump.dataX.push_back(freq);
dump.data.push_back(fData[i].dataFourierIm->GetBinContent(j));
}
}
// go through all theory bins
for (Int_t j=1; j<fData[i].theoryFourierIm->GetNbinsX(); j++) {
// get frequency
freq = fData[i].theoryFourierIm->GetBinCenter(j);
// check if time is in the current range
if ((freq >= xmin) && (freq <= xmax)) {
dump.theoryX.push_back(freq);
dump.theory.push_back(fData[i].theoryFourierIm->GetBinContent(j));
}
}
// if anything found keep it
if (dump.dataX.size() > 0)
dumpVector.push_back(dump);
} }
break; break;
case PV_FOURIER_PWR: case PV_FOURIER_PWR:
@ -1954,39 +1711,14 @@ void PMusrCanvas::ExportData(const Char_t *fileName)
xmax = fData[0].dataFourierPwr->GetXaxis()->GetBinCenter(xmaxBin); xmax = fData[0].dataFourierPwr->GetXaxis()->GetBinCenter(xmaxBin);
// fill ascii dump data // fill ascii dump data
for (UInt_t i=0; i<fData.size(); i++) { // go through all the histogramms if (fAveragedView) {
// clean up dump GetExportDataSet(fDataAvg.dataFourierPwr, xmin, xmax, dumpVector, false);
dump.dataX.clear(); GetExportDataSet(fDataAvg.theoryFourierPwr, xmin, xmax, dumpVector, false);
dump.data.clear(); } else { // go through all the histogramms
dump.dataErr.clear(); for (UInt_t i=0; i<fData.size(); i++) {
dump.theoryX.clear(); GetExportDataSet(fData[i].dataFourierPwr, xmin, xmax, dumpVector, false);
dump.theory.clear(); GetExportDataSet(fData[i].theoryFourierPwr, xmin, xmax, dumpVector, false);
// go through all data bins
for (Int_t j=1; j<fData[i].dataFourierPwr->GetNbinsX(); j++) {
// get frequency
freq = fData[i].dataFourierPwr->GetBinCenter(j);
// check if time is in the current range
if ((freq >= xmin) && (freq <= xmax)) {
dump.dataX.push_back(freq);
dump.data.push_back(fData[i].dataFourierPwr->GetBinContent(j));
}
} }
// go through all theory bins
for (Int_t j=1; j<fData[i].theoryFourierPwr->GetNbinsX(); j++) {
// get frequency
freq = fData[i].theoryFourierPwr->GetBinCenter(j);
// check if time is in the current range
if ((freq >= xmin) && (freq <= xmax)) {
dump.theoryX.push_back(freq);
dump.theory.push_back(fData[i].theoryFourierPwr->GetBinContent(j));
}
}
// if anything found keep it
if (dump.dataX.size() > 0)
dumpVector.push_back(dump);
} }
break; break;
case PV_FOURIER_PHASE: case PV_FOURIER_PHASE:
@ -1997,39 +1729,14 @@ void PMusrCanvas::ExportData(const Char_t *fileName)
xmax = fData[0].dataFourierPhase->GetXaxis()->GetBinCenter(xmaxBin); xmax = fData[0].dataFourierPhase->GetXaxis()->GetBinCenter(xmaxBin);
// fill ascii dump data // fill ascii dump data
for (UInt_t i=0; i<fData.size(); i++) { // go through all the histogramms if (fAveragedView) {
// clean up dump GetExportDataSet(fDataAvg.dataFourierPhase, xmin, xmax, dumpVector, false);
dump.dataX.clear(); GetExportDataSet(fDataAvg.theoryFourierPhase, xmin, xmax, dumpVector, false);
dump.data.clear(); } else { // go through all the histogramms
dump.dataErr.clear(); for (UInt_t i=0; i<fData.size(); i++) {
dump.theoryX.clear(); GetExportDataSet(fData[i].dataFourierPhase, xmin, xmax, dumpVector, false);
dump.theory.clear(); GetExportDataSet(fData[i].theoryFourierPhase, xmin, xmax, dumpVector, false);
// go through all data bins
for (Int_t j=1; j<fData[i].dataFourierPhase->GetNbinsX(); j++) {
// get frequency
freq = fData[i].dataFourierPhase->GetBinCenter(j);
// check if time is in the current range
if ((freq >= xmin) && (freq <= xmax)) {
dump.dataX.push_back(freq);
dump.data.push_back(fData[i].dataFourierPhase->GetBinContent(j));
}
} }
// go through all theory bins
for (Int_t j=1; j<fData[i].theoryFourierPhase->GetNbinsX(); j++) {
// get frequency
freq = fData[i].theoryFourierPhase->GetBinCenter(j);
// check if time is in the current range
if ((freq >= xmin) && (freq <= xmax)) {
dump.theoryX.push_back(freq);
dump.theory.push_back(fData[i].theoryFourierPhase->GetBinContent(j));
}
}
// if anything found keep it
if (dump.dataX.size() > 0)
dumpVector.push_back(dump);
} }
break; break;
default: default:
@ -2053,8 +1760,6 @@ void PMusrCanvas::ExportData(const Char_t *fileName)
dump.dataX.clear(); dump.dataX.clear();
dump.data.clear(); dump.data.clear();
dump.dataErr.clear(); dump.dataErr.clear();
dump.theoryX.clear();
dump.theory.clear();
// go through all data bins // go through all data bins
for (Int_t j=0; j<fNonMusrData[i].diff->GetN(); j++) { for (Int_t j=0; j<fNonMusrData[i].diff->GetN(); j++) {
@ -2102,8 +1807,6 @@ void PMusrCanvas::ExportData(const Char_t *fileName)
dump.dataX.clear(); dump.dataX.clear();
dump.data.clear(); dump.data.clear();
dump.dataErr.clear(); dump.dataErr.clear();
dump.theoryX.clear();
dump.theory.clear();
// go through all data bins // go through all data bins
for (Int_t j=0; j<fNonMusrData[i].data->GetN(); j++) { for (Int_t j=0; j<fNonMusrData[i].data->GetN(); j++) {
@ -2123,8 +1826,8 @@ void PMusrCanvas::ExportData(const Char_t *fileName)
fNonMusrData[i].theory->GetPoint(j,xval,yval); fNonMusrData[i].theory->GetPoint(j,xval,yval);
// check if time is in the current range // check if time is in the current range
if ((xval >= xmin) && (xval <= xmax)) { if ((xval >= xmin) && (xval <= xmax)) {
dump.theoryX.push_back(xval); dump.dataX.push_back(xval);
dump.theory.push_back(yval); dump.data.push_back(yval);
} }
} }
@ -2164,72 +1867,95 @@ void PMusrCanvas::ExportData(const Char_t *fileName)
} }
// find out what is the longest data/theory vector // find out what is the longest data/theory vector
UInt_t maxDataLength = 0; UInt_t maxLength = 0;
UInt_t maxTheoryLength = 0;
for (UInt_t i=0; i<dumpVector.size(); i++) { for (UInt_t i=0; i<dumpVector.size(); i++) {
if (maxDataLength < dumpVector[i].dataX.size()) if (maxLength < dumpVector[i].dataX.size())
maxDataLength = dumpVector[i].dataX.size(); maxLength = dumpVector[i].dataX.size();
if (maxTheoryLength < dumpVector[i].theoryX.size())
maxTheoryLength = dumpVector[i].theoryX.size();
} }
// write data to file // write data to file
UInt_t maxLength = 0;
if (fDifferenceView) { // difference view if (fDifferenceView) { // difference view
// write header // write header
switch (fCurrentPlotView) { switch (fCurrentPlotView) {
case PV_DATA: case PV_DATA:
fout << "% "; if (fAveragedView) {
for (UInt_t i=0; i<dumpVector.size()-1; i++) { fout << "% from averaged view" << endl;
fout << "x" << i << " , diff" << i << ", errDiff" << i << ", "; fout << "% x, diff, errDiff" << endl;
} else {
fout << "% ";
for (UInt_t i=0; i<dumpVector.size()-1; i++) {
fout << "x" << i << " , diff" << i << ", errDiff" << i << ", ";
}
fout << "x" << dumpVector.size()-1 << " , diff" << dumpVector.size()-1 << ", errDiff" << dumpVector.size()-1 << endl;
} }
fout << "x" << dumpVector.size()-1 << " , diff" << dumpVector.size()-1 << ", errDiff" << dumpVector.size()-1 << endl;
break; break;
case PV_FOURIER_REAL: case PV_FOURIER_REAL:
fout << "% "; if (fAveragedView) {
for (UInt_t i=0; i<dumpVector.size()-1; i++) { fout << "% from averaged view" << endl;
fout << "freq" << i << ", F_diffRe" << i << ", "; fout << "% x, F_diffRe" << endl;
} else {
fout << "% ";
for (UInt_t i=0; i<dumpVector.size()-1; i++) {
fout << "freq" << i << ", F_diffRe" << i << ", ";
}
fout << "freq" << dumpVector.size()-1 << ", F_diffRe" << dumpVector.size()-1 << endl;
} }
fout << "freq" << dumpVector.size()-1 << ", F_diffRe" << dumpVector.size()-1 << endl;
break; break;
case PV_FOURIER_IMAG: case PV_FOURIER_IMAG:
fout << "% "; if (fAveragedView) {
for (UInt_t i=0; i<dumpVector.size()-1; i++) { fout << "% from averaged view" << endl;
fout << "freq" << i << ", F_diffIm" << i << ", "; fout << "% x, F_diffIm" << endl;
} else {
fout << "% ";
for (UInt_t i=0; i<dumpVector.size()-1; i++) {
fout << "freq" << i << ", F_diffIm" << i << ", ";
}
fout << "freq" << dumpVector.size()-1 << ", F_diffIm" << dumpVector.size()-1 << endl;
} }
fout << "freq" << dumpVector.size()-1 << ", F_diffIm" << dumpVector.size()-1 << endl;
break; break;
case PV_FOURIER_REAL_AND_IMAG: case PV_FOURIER_REAL_AND_IMAG:
fout << "% "; if (fAveragedView) {
for (UInt_t i=0; i<dumpVector.size()/2; i++) { fout << "% from averaged view" << endl;
fout << "freq" << i << ", F_diffRe" << i << ", "; fout << "% x, F_diffRe, F_diffIm" << endl;
} else {
fout << "% ";
for (UInt_t i=0; i<dumpVector.size()/2; i++) {
fout << "freq" << i << ", F_diffRe" << i << ", ";
}
for (UInt_t i=0; i<dumpVector.size()/2-1; i++) {
fout << "freq" << i << ", F_diffIm" << i << ", ";
}
fout << "freq" << dumpVector.size()/2-1 << ", F_diffIm" << dumpVector.size()/2-1 << endl;
} }
for (UInt_t i=0; i<dumpVector.size()/2-1; i++) {
fout << "freq" << i << ", F_diffIm" << i << ", ";
}
fout << "freq" << dumpVector.size()/2-1 << ", F_diffIm" << dumpVector.size()/2-1 << endl;
break; break;
case PV_FOURIER_PWR: case PV_FOURIER_PWR:
fout << "% "; if (fAveragedView) {
for (UInt_t i=0; i<dumpVector.size()-1; i++) { fout << "% from averaged view" << endl;
fout << "freq" << i << ", F_diffPwr" << i << ", "; fout << "% x, F_diffPwr" << endl;
} else {
fout << "% ";
for (UInt_t i=0; i<dumpVector.size()-1; i++) {
fout << "freq" << i << ", F_diffPwr" << i << ", ";
}
fout << "freq" << dumpVector.size()-1 << ", F_diffPwr" << dumpVector.size()-1 << endl;
} }
fout << "freq" << dumpVector.size()-1 << ", F_diffPwr" << dumpVector.size()-1 << endl;
break; break;
case PV_FOURIER_PHASE: case PV_FOURIER_PHASE:
fout << "% "; if (fAveragedView) {
for (UInt_t i=0; i<dumpVector.size()-1; i++) { fout << "% from averaged view" << endl;
fout << "freq" << i << ", F_diffPhase" << i << ", "; fout << "% x, F_diffPhase" << endl;
} else {
fout << "% ";
for (UInt_t i=0; i<dumpVector.size()-1; i++) {
fout << "freq" << i << ", F_diffPhase" << i << ", ";
}
fout << "freq" << dumpVector.size()-1 << ", F_diffPhase" << dumpVector.size()-1 << endl;
} }
fout << "freq" << dumpVector.size()-1 << ", F_diffPhase" << dumpVector.size()-1 << endl;
break; break;
default: default:
break; break;
} }
maxLength = maxDataLength;
// write difference data // write difference data
for (UInt_t i=0; i<maxLength; i++) { for (UInt_t i=0; i<maxLength; i++) {
// write difference data // write difference data
@ -2264,84 +1990,119 @@ void PMusrCanvas::ExportData(const Char_t *fileName)
// write header // write header
switch (fCurrentPlotView) { switch (fCurrentPlotView) {
case PV_DATA: case PV_DATA:
fout << "% "; if (fAveragedView) {
for (UInt_t i=0; i<dumpVector.size(); i++) { fout << "% from averaged view" << endl;
fout << "xData" << i << " , data" << i << ", errData" << i << ", "; fout << "% xData, data, errData, xTheory, theory" << endl;
} else {
fout << "% ";
for (UInt_t i=0; i<dumpVector.size(); i++) {
if (i % 2 == 0)
fout << "xData" << i/2 << " , data" << i/2 << ", errData" << i/2 << ", ";
else
if (i == dumpVector.size()-1)
fout << "xTheory" << (i-1)/2 << " , theory" << (i-1)/2 << endl;
else
fout << "xTheory" << (i-1)/2 << " , theory" << (i-1)/2 << ", ";
}
} }
for (UInt_t i=0; i<dumpVector.size()-1; i++) {
fout << "xTheory" << i << " , theory" << i << ", ";
}
fout << "xTheory" << dumpVector.size()-1 << " , theory" << dumpVector.size()-1 << endl;
break; break;
case PV_FOURIER_REAL: case PV_FOURIER_REAL:
fout << "% "; if (fAveragedView) {
for (UInt_t i=0; i<dumpVector.size(); i++) { fout << "% from averaged view" << endl;
fout << "freq" << i << ", F_Re" << i << ", "; fout << "% freq, F_Re, freqTheo, F_theoRe" << endl;
} else {
fout << "% ";
for (UInt_t i=0; i<dumpVector.size(); i++) {
if (i % 2 == 0)
fout << "freq" << i/2 << ", F_Re" << i/2 << ", ";
else
if (i == dumpVector.size()-1)
fout << "freqTheo" << (i-1)/2 << ", F_theoRe" << (i-1)/2 << endl;
else
fout << "freqTheo" << (i-1)/2 << ", F_theoRe" << (i-1)/2 << ", ";
}
} }
for (UInt_t i=0; i<dumpVector.size()-1; i++) {
fout << "freqTheo" << i << ", F_theo" << i << ", ";
}
fout << "freqTheo" << dumpVector.size()-1 << ", F_theo" << dumpVector.size()-1 << endl;
break; break;
case PV_FOURIER_IMAG: case PV_FOURIER_IMAG:
fout << "% "; if (fAveragedView) {
for (UInt_t i=0; i<dumpVector.size(); i++) { fout << "% from averaged view" << endl;
fout << "freq" << i << ", F_Im" << i << ", "; fout << "% freq, F_Im, freqTheo, F_theoIm" << endl;
} else {
fout << "% ";
for (UInt_t i=0; i<dumpVector.size(); i++) {
if (i % 2 == 0)
fout << "freq" << i/2 << ", F_Im" << i/2 << ", ";
else
if (i == dumpVector.size()-1)
fout << "freqTheo" << (i-1)/2 << ", F_theoIm" << (i-1)/2 << endl;
else
fout << "freqTheo" << (i-1)/2 << ", F_theoIm" << (i-1)/2 << ", ";
}
} }
for (UInt_t i=0; i<dumpVector.size()-1; i++) {
fout << "freqTheo" << i << ", F_theo" << i << ", ";
}
fout << "freqTheo" << dumpVector.size()-1 << ", F_theo" << dumpVector.size()-1 << endl;
break; break;
case PV_FOURIER_REAL_AND_IMAG: case PV_FOURIER_REAL_AND_IMAG:
fout << "% "; if (fAveragedView) {
for (UInt_t i=0; i<dumpVector.size()/2; i++) { fout << "% from averaged view" << endl;
fout << "freq" << i << ", F_Re" << i << ", "; fout << "% freq, F_Re, freqTheo, F_theoRe, freq, F_Im, freqTheo, F_theoIm" << endl;
} else {
fout << "% ";
for (UInt_t i=0; i<dumpVector.size(); i++) {
if (i % 4 == 0)
fout << "freq" << i/4 << ", F_Re" << i/4 << ", ";
else if (i % 4 == 1)
fout << "freqTheo" << (i-1)/4 << ", F_theoRe" << (i-1)/4 << ", ";
else if (i % 4 == 2)
fout << "freq" << (i-2)/4 << ", F_Im" << (i-2)/4 << ", ";
else
if (i == dumpVector.size()-1)
fout << "freqTheo" << (i-3)/4 << ", F_theoIm" << (i-3)/4 << endl;
else
fout << "freqTheo" << (i-3)/4 << ", F_theoIm" << (i-3)/4 << ", ";
}
} }
for (UInt_t i=0; i<dumpVector.size()/2; i++) {
fout << "freq" << i << ", F_Im" << i << ", ";
}
for (UInt_t i=0; i<dumpVector.size()/2; i++) {
fout << "freqTheo" << i << ", F_theoRe" << i << ", ";
}
for (UInt_t i=0; i<(dumpVector.size()-1)/2; i++) {
fout << "freqTheo" << i << ", F_theoIm" << i << ", ";
}
fout << "freqTheo" << (dumpVector.size()-1)/2 << ", F_theoIm" << (dumpVector.size()-1)/2 << endl;
break; break;
case PV_FOURIER_PWR: case PV_FOURIER_PWR:
fout << "% "; if (fAveragedView) {
for (UInt_t i=0; i<dumpVector.size(); i++) { fout << "% from averaged view" << endl;
fout << "freq" << i << ", F_Pwr" << i << ", "; fout << "% freq, F_Pwr, freqTheo, F_theoPwr" << endl;
} else {
fout << "% ";
for (UInt_t i=0; i<dumpVector.size(); i++) {
if (i % 2 == 0)
fout << "freq" << i/2 << ", F_Pwr" << i/2 << ", ";
else
if (i == dumpVector.size()-1)
fout << "freqTheo" << (i-1)/2 << ", F_theoPwr" << (i-1)/2 << endl;
else
fout << "freqTheo" << (i-1)/2 << ", F_theoPwr" << (i-1)/2 << ", ";
}
} }
for (UInt_t i=0; i<dumpVector.size()-1; i++) {
fout << "freqTheo" << i << ", F_theo" << i << ", ";
}
fout << "freqTheo" << dumpVector.size()-1 << ", F_theo" << dumpVector.size()-1 << endl;
break; break;
case PV_FOURIER_PHASE: case PV_FOURIER_PHASE:
fout << "% "; if (fAveragedView) {
for (UInt_t i=0; i<dumpVector.size(); i++) { fout << "% from averaged view" << endl;
fout << "freq" << i << ", F_Phase" << i << ", "; fout << "% freq, F_Phase, freqTheo, F_theoPhase" << endl;
} else {
fout << "% ";
for (UInt_t i=0; i<dumpVector.size(); i++) {
if (i % 2 == 0)
fout << "freq" << i/2 << ", F_Phase" << i/2 << ", ";
else
if (i == dumpVector.size()-1)
fout << "freqTheo" << (i-1)/2 << ", F_theoPhase" << (i-1)/2 << endl;
else
fout << "freqTheo" << (i-1)/2 << ", F_theoPhase" << (i-1)/2 << ", ";
}
} }
for (UInt_t i=0; i<dumpVector.size()-1; i++) {
fout << "freqTheo" << i << ", F_theo" << i << ", ";
}
fout << "freqTheo" << dumpVector.size()-1 << ", F_theo" << dumpVector.size()-1 << endl;
break; break;
default: default:
break; break;
} }
if (maxDataLength > maxTheoryLength)
maxLength = maxDataLength;
else
maxLength = maxTheoryLength;
// write data and theory // write data and theory
for (UInt_t i=0; i<maxLength; i++) { for (UInt_t i=0; i<maxLength; i++) {
// write data // write data/theory
for (UInt_t j=0; j<dumpVector.size(); j++) { for (UInt_t j=0; j<dumpVector.size()-1; j++) {
if (i<dumpVector[j].dataX.size()) { if (i<dumpVector[j].dataX.size()) {
fout << dumpVector[j].dataX[i] << ", "; fout << dumpVector[j].dataX[i] << ", ";
fout << dumpVector[j].data[i] << ", "; fout << dumpVector[j].data[i] << ", ";
@ -2354,19 +2115,10 @@ void PMusrCanvas::ExportData(const Char_t *fileName)
fout << " , , "; fout << " , , ";
} }
} }
// write theory // write last data/theory entry
for (UInt_t j=0; j<dumpVector.size()-1; j++) { if (i<dumpVector[dumpVector.size()-1].dataX.size()) {
if (i<dumpVector[j].theoryX.size()) { fout << dumpVector[dumpVector.size()-1].dataX[i] << ", ";
fout << dumpVector[j].theoryX[i] << ", "; fout << dumpVector[dumpVector.size()-1].data[i];
fout << dumpVector[j].theory[i] << ", ";
} else {
fout << " , , ";
}
}
// write last theory entry
if (i<dumpVector[dumpVector.size()-1].theoryX.size()) {
fout << dumpVector[dumpVector.size()-1].theoryX[i] << ", ";
fout << dumpVector[dumpVector.size()-1].theory[i];
} else { } else {
fout << " , "; fout << " , ";
} }
@ -2382,8 +2134,6 @@ void PMusrCanvas::ExportData(const Char_t *fileName)
dumpVector[i].dataX.clear(); dumpVector[i].dataX.clear();
dumpVector[i].data.clear(); dumpVector[i].data.clear();
dumpVector[i].dataErr.clear(); dumpVector[i].dataErr.clear();
dumpVector[i].theoryX.clear();
dumpVector[i].theory.clear();
} }
dumpVector.clear(); dumpVector.clear();
@ -2394,6 +2144,42 @@ void PMusrCanvas::ExportData(const Char_t *fileName)
// } // }
} }
//--------------------------------------------------------------------------
// GetExportDataSet (private)
//--------------------------------------------------------------------------
/**
* <p> extract data for export.
*
* \param data
* \param xmin
* \param xmax
* \param dumpData
* \param hasError
*/
void PMusrCanvas::GetExportDataSet(const TH1F *data, const Double_t xmin, const Double_t xmax,
PMusrCanvasAsciiDumpVector &dumpData, const Bool_t hasError)
{
PMusrCanvasAsciiDump dump;
Double_t x=0.0;
// go through all difference data bins
for (Int_t j=1; j<data->GetNbinsX(); j++) {
// get time/freq
x = data->GetBinCenter(j);
// check if x is in the current range
if ((x >= xmin) && (x <= xmax)) {
dump.dataX.push_back(x);
dump.data.push_back(data->GetBinContent(j));
if (hasError)
dump.dataErr.push_back(data->GetBinError(j));
}
}
// if anything found keep it
if (dump.dataX.size() > 0)
dumpData.push_back(dump);
}
//-------------------------------------------------------------------------- //--------------------------------------------------------------------------
// CreateStyle (private) // CreateStyle (private)
//-------------------------------------------------------------------------- //--------------------------------------------------------------------------

5
src/include/PMusrCanvas.h
View File

@ -181,8 +181,6 @@ typedef struct {
PDoubleVector dataX; ///< x-axis data set PDoubleVector dataX; ///< x-axis data set
PDoubleVector data; ///< y-axis data set PDoubleVector data; ///< y-axis data set
PDoubleVector dataErr; ///< error of the y-axis data set PDoubleVector dataErr; ///< error of the y-axis data set
PDoubleVector theoryX; ///< x-axis theory set
PDoubleVector theory; ///< y-axis theory set
} PMusrCanvasAsciiDump; } PMusrCanvasAsciiDump;
//------------------------------------------------------------------------ //------------------------------------------------------------------------
@ -333,6 +331,9 @@ class PMusrCanvas : public TObject, public TQObject
virtual Double_t GetInterpolatedValue(TH1F* histo, Double_t xVal); virtual Double_t GetInterpolatedValue(TH1F* histo, Double_t xVal);
virtual void GetExportDataSet(const TH1F *data, const Double_t xmin, const Double_t xmax,
PMusrCanvasAsciiDumpVector &dumpData, const Bool_t hasError=true);
ClassDef(PMusrCanvas, 1) ClassDef(PMusrCanvas, 1)
}; };

67
src/tests/MuTransition/9900.msr Normal file
View File

@ -0,0 +1,67 @@
09900- Mu-frac 1.00, Mu12 134.86 MHz(0.27), Mu23 143.71 MHz(0.23), ionRate 608086.30 MHz, capRate 1.00 MHz, SF rate 0.00, 100 G
###############################################################
FITPARAMETER
# Nr. Name Value Step Pos_Error Boundaries
1 alpha 1.0008 -0.0021 0.0021 0 none
2 asy 0.2717 -0.0014 0.0014 0 0.33
3 phase 1.78 -0.46 0.46
4 field 100.418 -0.035 0.035 0 none
5 rate 0.0000000072 -0.0000000072 0.0013386264 0 100
6 asyMu 0 0 none
7 phaseMu 0 0 none
8 freqMu 35 0 none
9 rateMu 0 0 none
###############################################################
THEORY
asymmetry 2
TFieldCos 3 fun1 (phase frequency)
simplExpo 5 (rate)
+
asymmetry 6
TFieldCos 7 8 (phase frequency)
simplExpo 9 (rate)
###############################################################
FUNCTIONS
fun1 = par4 * gamma_mu
###############################################################
RUN 09900 MUE4 PSI MUSR-ROOT (name beamline institute data-file-format)
fittype 2 (asymmetry fit)
alpha 1
map 0 0 0 0 0 0 0 0 0 0
forward 1
backward 2
backgr.fix 0 0
data 1 12000 1 12000
t0 0.0 0.0
fit 0 8
packing 5
###############################################################
COMMANDS
SET BATCH
MINIMIZE
MINOS
SAVE
END RETURN
###############################################################
FOURIER
units MHz # units either 'Gauss', 'Tesla', 'MHz', or 'Mc/s'
fourier_power 12
apodization STRONG # NONE, WEAK, MEDIUM, STRONG
plot POWER # REAL, IMAG, REAL_AND_IMAG, POWER, PHASE
phase 8
#range_for_phase_correction 50.0 70.0
range 0 200
###############################################################
PLOT 2 (asymmetry plot)
runs 1
range 0 8 -0.35 0.35
###############################################################
STATISTIC --- 2016-02-17 20:38:53
chisq = 1457.5, NDF = 1595, chisq/NDF = 0.913780

841
src/tests/MuTransition/PSimulateMuTransition.cpp
View File

@ -1,382 +1,459 @@
/*************************************************************************** /***************************************************************************
PSimulateMuTransition.cpp PSimulateMuTransition.cpp
Author: Thomas Prokscha Author: Thomas Prokscha
Date: 25-Feb-2010 Date: 25-Feb-2010
$Id$ $Id$
Use root macros runMuSimulation.C and testAnalysis.C to run the simulation Use root macros runMuSimulation.C and testAnalysis.C to run the simulation
and to get a quick look on the data. Data are saved to a root histogram file and to get a quick look on the data. Data are saved to a root histogram file
with a structure similar to LEM histogram files; musrfit can be used to with a structure similar to LEM histogram files; musrfit can be used to
analyze the simulated data. analyze the simulated data.
Description: Description:
Root class to simulate muon spin phase under successive Mu+/Mu0 charge-exchange Root class to simulate muon spin phase under successive Mu+/Mu0 charge-exchange
processes by a Monte-Carlo method. Consider transverse field geometry, and assume processes by a Monte-Carlo method. Consider transverse field geometry, and assume
initial muon spin direction in x, and field applied along z. For PxMu(t) in initial muon spin direction in x, and field applied along z. For PxMu(t) in
muonium use the equation 8.22 of the muSR book of Yaounc and Dalmas de Réotier, in muonium use the equation 8.22 of the muSR book of Yaounc and Dalmas de Réotier, in
slightly modified form (see Senba, J. Phys. B 23, 1545 (1990)); note that PxMu(t) slightly modified form (see Senba, J. Phys. B 23, 1545 (1990)); note that PxMu(t)
is given by a superposition of the four frequencies "nu_12", "nu_34", "nu_23", "nu_14". is given by a superposition of the four frequencies "nu_12", "nu_34", "nu_23", "nu_14".
These frequencies and the corresponding probabilities ("SetMuFractionState12" for These frequencies and the corresponding probabilities ("SetMuFractionState12" for
transitions 12 and 34, "SetMuFractionState23" for states 23 and 14) can be calculated transitions 12 and 34, "SetMuFractionState23" for states 23 and 14) can be calculated
for a given field with the root macro AnisotropicMu.C for a given field with the root macro AnisotropicMu.C
Parameters: Parameters:
1) Precession frequencies of "nu_12", "nu_34", "nu_23", "nu_14" 1) Precession frequencies of "nu_12", "nu_34", "nu_23", "nu_14"
2) fractions of nu_12, nu_34; and nu_23 and nu_14 2) fractions of nu_12, nu_34; and nu_23 and nu_14
3) total Mu0 fraction 3) total Mu0 fraction
4) electron-capture rate 4) electron-capture rate
5) Mu ionization rate 5) Mu ionization rate
6) initial muon spin phase 6) initial muon spin phase
7) total muon decay asymmetry 7) total muon decay asymmetry
8) number of muon decays to be generated. 8) number of muon decays to be generated.
9) debug flag: if TRUE print capture/ionization events on screen 9) debug flag: if TRUE print capture/ionization events on screen
Output: Output:
Two histograms ("forward" and "backward") are written to a root file. Two histograms ("forward" and "backward") are written to a root file.
The muon event simulation with a sequence of charge-changing processes is The muon event simulation with a sequence of charge-changing processes is
done in Event(): done in Event():
simulate muon spin phase under charge-exchange with "4 Mu transitions" simulate muon spin phase under charge-exchange with "4 Mu transitions"
1) according to Mu+/Mu0 fraction begin either with a Mu+ state or Mu state 1) according to Mu+/Mu0 fraction begin either with a Mu+ state or Mu state
2) Mu+: determine next electron-capture time t_c. If t_c is larger than decay time t_d 2) Mu+: determine next electron-capture time t_c. If t_c is larger than decay time t_d
calculate muon spin precession for t_d; else calculate spin precession for t_c. calculate muon spin precession for t_d; else calculate spin precession for t_c.
3) Determine next ionization time t_i; calculate Px(t_i) in Muonium; calculate the 3) Determine next ionization time t_i; calculate Px(t_i) in Muonium; calculate the
muon spin phase by acos(Px(t_i)). muon spin phase by acos(Px(t_i)).
4) get the next electron capture time, continue until t_d is reached; accumulate muon spin 4) get the next electron capture time, continue until t_d is reached; accumulate muon spin
phase. phase.
***************************************************************************/ ***************************************************************************/
/*************************************************************************** /***************************************************************************
* Copyright (C) 2010 by Thomas Prokscha, Paul Scherrer Institut * * Copyright (C) 2010 by Thomas Prokscha, Paul Scherrer Institut *
* * * *
* This program is free software; you can redistribute it and/or modify * * This program is free software; you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by * * it under the terms of the GNU General Public License as published by *
* the Free Software Foundation; either version 2 of the License, or * * the Free Software Foundation; either version 2 of the License, or *
* (at your option) any later version. * * (at your option) any later version. *
* * * *
* This program is distributed in the hope that it will be useful, * * This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of * * but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. * * GNU General Public License for more details. *
* * * *
* You should have received a copy of the GNU General Public License * * You should have received a copy of the GNU General Public License *
* along with this program; if not, write to the * * along with this program; if not, write to the *
* Free Software Foundation, Inc., * * Free Software Foundation, Inc., *
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. * * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
***************************************************************************/ ***************************************************************************/
#include <iostream> #include <iostream>
using namespace std; using namespace std;
#include <TMath.h> #include <TMath.h>
#include <TComplex.h>
#include "PSimulateMuTransition.h"
#include "PSimulateMuTransition.h"
ClassImp(PSimulateMuTransition)
ClassImp(PSimulateMuTransition)
//--------------------------------------------------------------------------
// Constructor //--------------------------------------------------------------------------
//-------------------------------------------------------------------------- // Constructor
/** //--------------------------------------------------------------------------
* <p> Constructor. /**
* * <p> Constructor.
* \param seed for the random number generator *
*/ * \param seed for the random number generator
PSimulateMuTransition::PSimulateMuTransition(UInt_t seed) */
{ PSimulateMuTransition::PSimulateMuTransition(UInt_t seed)
fValid = true; {
fValid = true;
fRandom = new TRandom2(seed);
if (fRandom == 0) { fRandom = new TRandom2(seed);
fValid = false; if (fRandom == 0) {
} fValid = false;
}
fNmuons = 100; // number of muons to simulate
fMuPrecFreq34 = 4463.; // vacuum Mu hyperfine coupling constant fNmuons = 100; // number of muons to simulate
fMuPrecFreq12 = 0.; // Mu precession frequency of a 12 transition fMuPrecFreq34 = 4463.; // vacuum Mu hyperfine coupling constant
fMuPrecFreq23 = 0.; // Mu precession frequency of a 23 transition fMuPrecFreq12 = 0.; // Mu precession frequency of a 12 transition
fMuPrecFreq14 = 0.; // Mu precession frequency of a 14 transition fMuPrecFreq23 = 0.; // Mu precession frequency of a 23 transition
fMuonPrecFreq = 0.; // muon precession frequency fMuPrecFreq14 = 0.; // Mu precession frequency of a 14 transition
fBfield = 0.01; // magnetic field (T) fMuonPrecFreq = 0.; // muon precession frequency
fCaptureRate = 0.01; // Mu+ capture rate (MHz) fBfield = 0.01; // magnetic field (T)
fIonizationRate = 10.; // Mu0 ionization rate (MHz) fCaptureRate = 0.01; // Mu+ capture rate (MHz)
fInitialPhase = 0.; fIonizationRate = 10.; // Mu0 ionization rate (MHz)
fMuonPhase = fInitialPhase; fSpinFlipRate = 0.001; // Mu0 spin flip rate (MHz)
fMuonDecayTime = 0.; fInitialPhase = 0.;
fAsymmetry = 0.27; fMuonPhase = fInitialPhase;
fMuFraction = 0.; fMuonDecayTime = 0.;
fMuFractionState12 = 0.; fAsymmetry = 0.27;
fMuFractionState23 = 0.; fMuFraction = 0.;
fDebugFlag = kFALSE; fMuFractionState12 = 0.;
} fMuFractionState23 = 0.;
fDebugFlag = kFALSE;
//-------------------------------------------------------------------------- }
// Destructor
//-------------------------------------------------------------------------- //--------------------------------------------------------------------------
/** // Destructor
* <p> Destructor. //--------------------------------------------------------------------------
* /**
*/ * <p> Destructor.
PSimulateMuTransition::~PSimulateMuTransition() *
{ */
if (fRandom) { PSimulateMuTransition::~PSimulateMuTransition()
delete fRandom; {
fRandom = 0; if (fRandom) {
} delete fRandom;
} fRandom = 0;
//-------------------------------------------------------------------------- }
// Output of current settings }
//-------------------------------------------------------------------------- //--------------------------------------------------------------------------
/*! // Output of current settings
* <p>Prints the current settings onto std output. //--------------------------------------------------------------------------
*/ /*!
void PSimulateMuTransition::PrintSettings() const * <p>Prints the current settings onto std output.
{ */
cout << endl << "Mu precession frequency 12 (MHz) = " << fMuPrecFreq12; void PSimulateMuTransition::PrintSettings() const
cout << endl << "Mu precession frequency 34 (MHz) = " << fMuPrecFreq34; {
cout << endl << "Mu precession frequency 23 (MHz) = " << fMuPrecFreq23; cout << endl << "Mu precession frequency 12 (MHz) = " << fMuPrecFreq12;
cout << endl << "Mu precession frequency 14 (MHz) = " << fMuPrecFreq14; cout << endl << "Mu precession frequency 34 (MHz) = " << fMuPrecFreq34;
cout << endl << "B field (T) = " << fBfield; cout << endl << "Mu precession frequency 23 (MHz) = " << fMuPrecFreq23;
cout << endl << "Mu+ electron capture rate (MHz) = " << fCaptureRate; cout << endl << "Mu precession frequency 14 (MHz) = " << fMuPrecFreq14;
cout << endl << "Mu ionizatioan rate (MHz) = " << fIonizationRate; cout << endl << "B field (T) = " << fBfield;
cout << endl << "Decay asymmetry = " << fAsymmetry; cout << endl << "Mu+ electron capture rate (MHz) = " << fCaptureRate;
cout << endl << "Muonium fraction = " << fMuFraction; cout << endl << "Mu0 ionizatioan rate (MHz) = " << fIonizationRate;
cout << endl << "Muonium fraction state12 = " << fMuFractionState12; cout << endl << "Mu0 spin-flip rate (MHz) = " << fSpinFlipRate;
cout << endl << "Muonium fraction state23 = " << fMuFractionState23; cout << endl << "!!! Note: if spin-flip rate > 0.001 only spin-flip process is considered!!!";
cout << endl << "Number of particles to simulate = " << fNmuons; cout << endl << "Decay asymmetry = " << fAsymmetry;
cout << endl << "Initial muon spin phase (degree) = " << fInitialPhase; cout << endl << "Muonium fraction = " << fMuFraction;
cout << endl << "Debug flag = " << fDebugFlag; cout << endl << "Muonium fraction state12 = " << fMuFractionState12;
cout << endl << endl; cout << endl << "Muonium fraction state23 = " << fMuFractionState23;
} cout << endl << "Number of particles to simulate = " << fNmuons;
cout << endl << "Initial muon spin phase (degree) = " << fInitialPhase;
//-------------------------------------------------------------------------- cout << endl << "Debug flag = " << fDebugFlag;
// SetSeed (public) cout << endl << endl;
//-------------------------------------------------------------------------- }
/**
* <p>Sets the seed of the random generator. //--------------------------------------------------------------------------
* // SetSeed (public)
* \param seed for the random number generator //--------------------------------------------------------------------------
*/ /**
void PSimulateMuTransition::SetSeed(UInt_t seed) * <p>Sets the seed of the random generator.
{ *
if (!fValid) * \param seed for the random number generator
return; */
void PSimulateMuTransition::SetSeed(UInt_t seed)
fRandom->SetSeed(seed); {
} if (!fValid)
return;
//--------------------------------------------------------------------------
// Run (public) fRandom->SetSeed(seed);
//-------------------------------------------------------------------------- }
/**
* \param histoForward //--------------------------------------------------------------------------
*/ // Run (public)
void PSimulateMuTransition::Run(TH1F *histoForward, TH1F *histoBackward) //--------------------------------------------------------------------------
{ /**
Int_t i; * \param histoForward
if (histoForward == 0 || histoBackward == 0) */
return; void PSimulateMuTransition::Run(TH1F *histoForward, TH1F *histoBackward)
{
for (i = 0; i<fNmuons; i++){ Int_t i;
fMuonPhase = TMath::TwoPi() * fInitialPhase/360.; // transform to radians if (histoForward == 0 || histoBackward == 0)
fMuonDecayTime = NextEventTime(fMuonDecayRate); return;
// initial muon state Mu+ or Mu0?
if (fRandom->Rndm() <= 1.-fMuFraction) for (i = 0; i<fNmuons; i++){
Event("Mu+"); fMuonPhase = TMath::TwoPi() * fInitialPhase/360.; // transform to radians
else fMuonDecayTime = NextEventTime(fMuonDecayRate);
Event("");
if (fSpinFlipRate > 0.001){// consider only Mu0 spin-flip in this case
// fill 50% in "forward", and 50% in "backward" detector to get independent fMuonPhase = TMath::ACos(GTSpinFlip(fMuonDecayTime));
// events in "forward" and "backward" histograms. This allows "normal" uSR }
// analysis of the data else{
// change muon decay time to ns // initial muon state Mu+ or Mu0?
if (fRandom->Rndm() <= 0.5) if (fRandom->Rndm() <= 1.-fMuFraction)
histoForward->Fill(fMuonDecayTime*1000., 1. + fAsymmetry*TMath::Cos(fMuonPhase)); Event("Mu+");
else else
histoBackward->Fill(fMuonDecayTime*1000., 1. - fAsymmetry*TMath::Cos(fMuonPhase)); Event("");
}
if ( (i%100000) == 0) cout << "number of events processed: " << i << endl; // fill 50% in "forward", and 50% in "backward" detector to get independent
} // events in "forward" and "backward" histograms. This allows "normal" uSR
cout << "number of events processed: " << i << endl; // analysis of the data
return; // change muon decay time to ns
} if (fRandom->Rndm() <= 0.5)
histoForward->Fill(fMuonDecayTime*1000., 1. + fAsymmetry*TMath::Cos(fMuonPhase));
//-------------------------------------------------------------------------- else
// NextEventTime (private) histoBackward->Fill(fMuonDecayTime*1000., 1. - fAsymmetry*TMath::Cos(fMuonPhase));
//--------------------------------------------------------------------------
/** if ( (i%100000) == 0) cout << "number of events processed: " << i << endl;
* <p>Determine time of next event, assuming "Poisson" distribution in time }
* cout << "number of events processed: " << i << endl;
* \param EventRate event rate in MHz; returns next event time in micro-seconds return;
*/ }
Double_t PSimulateMuTransition::NextEventTime(const Double_t &EventRate)
{ //--------------------------------------------------------------------------
if (EventRate <= 0.) // NextEventTime (private)
return -1.; // signal error //--------------------------------------------------------------------------
/**
return -1./EventRate * TMath::Log(fRandom->Rndm()); * <p>Determine time of next event, assuming "Poisson" distribution in time
} *
* \param EventRate event rate in MHz; returns next event time in micro-seconds
//-------------------------------------------------------------------------- */
// Phase (private) Double_t PSimulateMuTransition::NextEventTime(const Double_t &EventRate)
//-------------------------------------------------------------------------- {
/** if (EventRate <= 0.)
* <p>Determines phase of the muon spin return -1.; // signal error
*
* \param time duration of precession (us); return -1./EventRate * TMath::Log(fRandom->Rndm());
* \param frequency muon spin precession frequency (MHz); }
*/
Double_t PSimulateMuTransition::PrecessionPhase(const Double_t &time, const TString chargeState) //--------------------------------------------------------------------------
{ // Phase (private)
Double_t muonPhaseX; //--------------------------------------------------------------------------
Double_t muoniumPolX = 0; /**
* <p>Determines phase of the muon spin
if (chargeState == "Mu+") *
muonPhaseX = TMath::TwoPi()*fMuonPrecFreq*time; * \param time duration of precession (us);
else if (chargeState == "Mu0"){ * \param chargeState charge state of Mu ("Mu+" or "Mu0")
muoniumPolX = 0.5 * */
(fMuFractionState12 * (TMath::Cos(TMath::TwoPi()*fMuPrecFreq12*time) + TMath::Cos(TMath::TwoPi()*fMuPrecFreq34*time)) + Double_t PSimulateMuTransition::PrecessionPhase(const Double_t &time, const TString chargeState)
fMuFractionState23 * (TMath::Cos(TMath::TwoPi()*fMuPrecFreq23*time) + TMath::Cos(TMath::TwoPi()*fMuPrecFreq14*time))); {
muonPhaseX = TMath::ACos(muoniumPolX); Double_t muonPhaseX;
} Double_t muoniumPolX = 0;
else
muonPhaseX = 0.; if (chargeState == "Mu+")
muonPhaseX = TMath::TwoPi()*fMuonPrecFreq*time;
return muonPhaseX; else if (chargeState == "Mu0"){
} muoniumPolX = GTFunction(time).Re();
muonPhaseX = TMath::ACos(muoniumPolX);
//-------------------------------------------------------------------------- }
// Event (private) else
//-------------------------------------------------------------------------- muonPhaseX = 0.;
/**
* <p> Generates "muon event": simulate muon spin phase under charge-exchange with return muonPhaseX;
* a neutral muonium state in transverse field, where the polarization evolution }
* PxMu(t) of the muon spin in muonium is determined by a superposition of the
* four "Mu transitions" nu_12, nu_34, nu_23, and nu_14. //--------------------------------------------------------------------------
* 1) according to Mu+/Mu0 fraction begin either with a Mu+ state or Mu state // Mu0 transverse field polarization function (private)
* 2) Mu+: determine next electron-capture time t_c. If t_c is larger than decay time t_d //--------------------------------------------------------------------------
* calculate muon spin precession for t_d; else calculate spin precession for t_c. /**
* 3) Determine next ionization time t_i; calculate Px(t_i) in Muonium; calculate the * <p>Calculates Mu0 polarization in x direction by superposition of four Mu0 frequencies
* muon spin phase by acos(Px(t_i)). *
* 4) get the next electron capture time, continue until t_d is reached. * \param time (us);
* */
* <p> For isotropic muonium, TF: TComplex PSimulateMuTransition::GTFunction(const Double_t &time)
* nu_12 and nu_34 with equal probabilities, probability for both states fMuFractionState12 {
* ni_23 and nu_14 with equal probabilities, probability for both states fMuFractionState23 Double_t twoPi = TMath::TwoPi();
*
* \param muonString if eq. "Mu+" begin with Mu+ precession TComplex complexPol = 0;
*/ complexPol =
void PSimulateMuTransition::Event(const TString muonString) 0.5 * fMuFractionState12 *
{ (TComplex::Exp(TComplex::I()*twoPi*fMuPrecFreq12*time) +
Double_t eventTime, eventDiffTime, captureTime, ionizationTime; TComplex::Exp(-TComplex::I()*twoPi*fMuPrecFreq34*time))
// Double_t muonPrecessionFreq, muoniumPrecessionFreq; // MHz +
// Double_t rndm, frac1, frac2; 0.5 * fMuFractionState23 *
(TComplex::Exp(TComplex::I()*twoPi*fMuPrecFreq23*time) +
fMuonPrecFreq = fMuonGyroRatio * fBfield; TComplex::Exp(TComplex::I()*twoPi*fMuPrecFreq14*time));
// charge-exchange loop until muon decay return complexPol;
eventTime = 0.;
eventDiffTime = 0.; // Double_t muoniumPolX = 0;
// muoniumPolX = 0.5 *
if (fDebugFlag) cout << "Decay time = " << fMuonDecayTime << endl; // (fMuFractionState12 * (TMath::Cos(twoPi*fMuPrecFreq12*time) + TMath::Cos(twoPi*fMuPrecFreq34*time)) +
//cout << muonString << endl; // fMuFractionState23 * (TMath::Cos(twoPi*fMuPrecFreq23*time) + TMath::Cos(twoPi*fMuPrecFreq14*time)));
while (1) { //
if (muonString == "Mu+"){ // return muoniumPolX;
// Mu+ initial state; get next electron capture time
captureTime = NextEventTime(fCaptureRate); }
eventTime += captureTime;
if (fDebugFlag) cout << "Capture time = " << captureTime << " Phase = " << fMuonPhase << endl; //--------------------------------------------------------------------------
if (eventTime < fMuonDecayTime) // Mu0 transverse field polarization function after n spin-flip collisions (private)
fMuonPhase += PrecessionPhase(captureTime, "Mu+"); //--------------------------------------------------------------------------
else{ //muon decays; handle precession prior to muon decay /**
eventDiffTime = fMuonDecayTime - (eventTime - captureTime); * <p>Calculates Mu0 polarization in x direction after n spin flip collisions.
fMuonPhase += PrecessionPhase(eventDiffTime, "Mu+"); * See M. Senba, J.Phys. B24, 3531 (1991), equation (17)
break; *
} * \param time (us);
*/
// now, we have Mu0; get next ionization time Double_t PSimulateMuTransition::GTSpinFlip(const Double_t &time)
ionizationTime = NextEventTime(fIonizationRate); {
eventTime += ionizationTime; TComplex complexPolX = 1.0;
// determine Mu state Double_t muoniumPolX = 1.0; //initial polarization in x direction
// rndm = fRandom->Rndm(); Double_t eventTime = 0;
// frac1 = 1. - fMuFractionState1 - fMuFractionState2; // non-precessing Mu states Double_t eventDiffTime = 0;
// frac2 = 1. - fMuFractionState2; Double_t lastEventTime = 0;
// if ( rndm < frac1 )
// muoniumPrecessionFreq = 0.; eventTime += NextEventTime(fSpinFlipRate);
// else if (rndm >= frac1 && rndm <= frac2){ if (eventTime >= time){
// if (fRandom->Rndm() <= 0.5) muoniumPolX = GTFunction(time).Re();
// muoniumPrecessionFreq = fMuPrecFreq12; }
// else else{
// muoniumPrecessionFreq = fMuPrecFreq34; while (eventTime < time){
// } eventDiffTime = eventTime - lastEventTime;
// else{ complexPolX = complexPolX * GTFunction(eventDiffTime);
// if (fRandom->Rndm() <= 0.5) lastEventTime = eventTime;
// muoniumPrecessionFreq = fMuPrecFreq23; eventTime += NextEventTime(fSpinFlipRate);
// else }
// muoniumPrecessionFreq = fMuPrecFreq14; // calculate for the last collision
// } eventDiffTime = time - lastEventTime;
complexPolX = complexPolX * GTFunction(eventDiffTime);
if (fDebugFlag) cout << "Ioniza. time = " << ionizationTime << " Phase = " << fMuonPhase << endl; muoniumPolX = complexPolX.Re();
if (eventTime < fMuonDecayTime) }
fMuonPhase += PrecessionPhase(ionizationTime, "Mu0");
else{ //muon decays; handle precession prior to muon decay return muoniumPolX;
eventDiffTime = fMuonDecayTime - (eventTime - ionizationTime); }
fMuonPhase += PrecessionPhase(eventDiffTime, "Mu0");
break; //--------------------------------------------------------------------------
} // Event (private)
} //--------------------------------------------------------------------------
else{ /**
// Mu0 as initial state; get next ionization time * <p> Generates "muon event": simulate muon spin phase under charge-exchange with
ionizationTime = NextEventTime(fIonizationRate); * a neutral muonium state in transverse field, where the polarization evolution
eventTime += ionizationTime; * PxMu(t) of the muon spin in muonium is determined by a superposition of the
// determine Mu state * four "Mu transitions" nu_12, nu_34, nu_23, and nu_14.
// rndm = fRandom->Rndm(); * 1) according to Mu+/Mu0 fraction begin either with a Mu+ state or Mu state
// frac1 = 1. - fMuFractionState1 - fMuFractionState2; // non-precessing Mu states * 2) Mu+: determine next electron-capture time t_c. If t_c is larger than decay time t_d
// frac2 = 1. - fMuFractionState2; * calculate muon spin precession for t_d; else calculate spin precession for t_c.
// if ( rndm < frac1 ) * 3) Determine next ionization time t_i; calculate Px(t_i) in Muonium; calculate the
// muoniumPrecessionFreq = 0.; * muon spin phase by acos(Px(t_i)).
// else if (rndm >= frac1 && rndm <= frac2){ * 4) get the next electron capture time, continue until t_d is reached.
// if (fRandom->Rndm() <= 0.5) *
// muoniumPrecessionFreq = fMuPrecFreq12; * <p> For isotropic muonium, TF:
// else * nu_12 and nu_34 with equal probabilities, probability for both states fMuFractionState12
// muoniumPrecessionFreq = fMuPrecFreq34; * ni_23 and nu_14 with equal probabilities, probability for both states fMuFractionState23
// } *
// else{ * \param muonString if eq. "Mu+" begin with Mu+ precession
// if (fRandom->Rndm() <= 0.5) */
// muoniumPrecessionFreq = fMuPrecFreq23; void PSimulateMuTransition::Event(const TString muonString)
// else {
// muoniumPrecessionFreq = fMuPrecFreq14; Double_t eventTime, eventDiffTime, captureTime, ionizationTime;
// } // Double_t muonPrecessionFreq, muoniumPrecessionFreq; // MHz
// Double_t rndm, frac1, frac2;
if (fDebugFlag)
cout << "Mu Ioniza. time = " << ionizationTime << " Phase = " << fMuonPhase << endl; fMuonPrecFreq = fMuonGyroRatio * fBfield;
if (eventTime < fMuonDecayTime)
fMuonPhase += PrecessionPhase(ionizationTime, "Mu0"); // charge-exchange loop until muon decay
else{ //muon decays; handle precession prior to muon decay eventTime = 0.;
eventDiffTime = fMuonDecayTime - (eventTime - ionizationTime); eventDiffTime = 0.;
fMuonPhase += PrecessionPhase(eventDiffTime, "Mu0");
break; if (fDebugFlag) cout << "Decay time = " << fMuonDecayTime << endl;
} //cout << muonString << endl;
while (1) {
// Mu+ state; get next electron capture time if (muonString == "Mu+"){
captureTime = NextEventTime(fCaptureRate); // Mu+ initial state; get next electron capture time
eventTime += captureTime; captureTime = NextEventTime(fCaptureRate);
if (fDebugFlag) cout << "Capture time = " << captureTime << " Phase = " << fMuonPhase << endl; eventTime += captureTime;
if (eventTime < fMuonDecayTime) if (fDebugFlag) cout << "Capture time = " << captureTime << " Phase = " << fMuonPhase << endl;
fMuonPhase += PrecessionPhase(captureTime, "Mu+"); if (eventTime < fMuonDecayTime)
else{ //muon decays; handle precession prior to muon decay fMuonPhase += PrecessionPhase(captureTime, "Mu+");
eventDiffTime = fMuonDecayTime - (eventTime - captureTime); else{ //muon decays; handle precession prior to muon decay
fMuonPhase += PrecessionPhase(eventDiffTime, "Mu+"); eventDiffTime = fMuonDecayTime - (eventTime - captureTime);
break; fMuonPhase += PrecessionPhase(eventDiffTime, "Mu+");
} break;
} }
}
// now, we have Mu0; get next ionization time
if (fDebugFlag) cout << " Final Phase = " << fMuonPhase << endl; ionizationTime = NextEventTime(fIonizationRate);
//fMuonPhase = TMath::ACos(TMath::Cos(fMuonPhase))*360./TMath::TwoPi(); //transform back to [0, 180] degree interval eventTime += ionizationTime;
return; // determine Mu state
} // rndm = fRandom->Rndm();
// frac1 = 1. - fMuFractionState1 - fMuFractionState2; // non-precessing Mu states
// frac2 = 1. - fMuFractionState2;
// if ( rndm < frac1 )
// muoniumPrecessionFreq = 0.;
// else if (rndm >= frac1 && rndm <= frac2){
// if (fRandom->Rndm() <= 0.5)
// muoniumPrecessionFreq = fMuPrecFreq12;
// else
// muoniumPrecessionFreq = fMuPrecFreq34;
// }
// else{
// if (fRandom->Rndm() <= 0.5)
// muoniumPrecessionFreq = fMuPrecFreq23;
// else
// muoniumPrecessionFreq = fMuPrecFreq14;
// }
if (fDebugFlag) cout << "Ioniza. time = " << ionizationTime << " Phase = " << fMuonPhase << endl;
if (eventTime < fMuonDecayTime)
fMuonPhase += PrecessionPhase(ionizationTime, "Mu0");
else{ //muon decays; handle precession prior to muon decay
eventDiffTime = fMuonDecayTime - (eventTime - ionizationTime);
fMuonPhase += PrecessionPhase(eventDiffTime, "Mu0");
break;
}
}
else{
// Mu0 as initial state; get next ionization time
ionizationTime = NextEventTime(fIonizationRate);
eventTime += ionizationTime;
// determine Mu state
// rndm = fRandom->Rndm();
// frac1 = 1. - fMuFractionState1 - fMuFractionState2; // non-precessing Mu states
// frac2 = 1. - fMuFractionState2;
// if ( rndm < frac1 )
// muoniumPrecessionFreq = 0.;
// else if (rndm >= frac1 && rndm <= frac2){
// if (fRandom->Rndm() <= 0.5)
// muoniumPrecessionFreq = fMuPrecFreq12;
// else
// muoniumPrecessionFreq = fMuPrecFreq34;
// }
// else{
// if (fRandom->Rndm() <= 0.5)
// muoniumPrecessionFreq = fMuPrecFreq23;
// else
// muoniumPrecessionFreq = fMuPrecFreq14;
// }
if (fDebugFlag)
cout << "Mu Ioniza. time = " << ionizationTime << " Phase = " << fMuonPhase << endl;
if (eventTime < fMuonDecayTime)
fMuonPhase += PrecessionPhase(ionizationTime, "Mu0");
else{ //muon decays; handle precession prior to muon decay
eventDiffTime = fMuonDecayTime - (eventTime - ionizationTime);
fMuonPhase += PrecessionPhase(eventDiffTime, "Mu0");
break;
}
// Mu+ state; get next electron capture time
captureTime = NextEventTime(fCaptureRate);
eventTime += captureTime;
if (fDebugFlag) cout << "Capture time = " << captureTime << " Phase = " << fMuonPhase << endl;
if (eventTime < fMuonDecayTime)
fMuonPhase += PrecessionPhase(captureTime, "Mu+");
else{ //muon decays; handle precession prior to muon decay
eventDiffTime = fMuonDecayTime - (eventTime - captureTime);
fMuonPhase += PrecessionPhase(eventDiffTime, "Mu+");
break;
}
}
}
if (fDebugFlag) cout << " Final Phase = " << fMuonPhase << endl;
//fMuonPhase = TMath::ACos(TMath::Cos(fMuonPhase))*360./TMath::TwoPi(); //transform back to [0, 180] degree interval
return;
}

207
src/tests/MuTransition/PSimulateMuTransition.h
View File

@ -1,101 +1,106 @@
/*************************************************************************** /***************************************************************************
PSimulateMuTransition.h PSimulateMuTransition.h
Author: Thomas Prokscha Author: Thomas Prokscha
Date: 25-Feb-2010 Date: 25-Feb-2010
$Id$ $Id$
***************************************************************************/ ***************************************************************************/
/*************************************************************************** /***************************************************************************
* Copyright (C) 2010 by Thomas Prokscha, Paul Scherrer Institut * * Copyright (C) 2010 by Thomas Prokscha, Paul Scherrer Institut *
* * * *
* This program is free software; you can redistribute it and/or modify * * This program is free software; you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by * * it under the terms of the GNU General Public License as published by *
* the Free Software Foundation; either version 2 of the License, or * * the Free Software Foundation; either version 2 of the License, or *
* (at your option) any later version. * * (at your option) any later version. *
* * * *
* This program is distributed in the hope that it will be useful, * * This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of * * but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. * * GNU General Public License for more details. *
* * * *
* You should have received a copy of the GNU General Public License * * You should have received a copy of the GNU General Public License *
* along with this program; if not, write to the * * along with this program; if not, write to the *
* Free Software Foundation, Inc., * * Free Software Foundation, Inc., *
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. * * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
***************************************************************************/ ***************************************************************************/
#ifndef _PSIMULATEMUTRANSITION_H_ #ifndef _PSIMULATEMUTRANSITION_H_
#define _PSIMULATEMUTRANSITION_H_ #define _PSIMULATEMUTRANSITION_H_
#include <TObject.h> #include <TObject.h>
#include <TH1F.h> #include <TH1F.h>
#include <TRandom2.h> #include <TRandom2.h>
#include <TComplex.h>
// global constants
const Double_t fMuonGyroRatio = 135.54; //!< muon gyromagnetic ratio (MHz/T) // global constants
const Double_t fMuonDecayRate = 0.4551; //!< muon decay rate (1/tau_mu, MHz) const Double_t fMuonGyroRatio = 135.54; //!< muon gyromagnetic ratio (MHz/T)
const Double_t fMuonDecayRate = 0.4551; //!< muon decay rate (1/tau_mu, MHz)
class PSimulateMuTransition : public TObject
{ class PSimulateMuTransition : public TObject
public: {
PSimulateMuTransition(UInt_t seed = 0); public:
virtual ~PSimulateMuTransition(); PSimulateMuTransition(UInt_t seed = 0);
virtual ~PSimulateMuTransition();
virtual void PrintSettings() const;
virtual void SetNmuons(Int_t value) { fNmuons = value; } //!< number of muons virtual void PrintSettings() const;
virtual void SetDebugFlag(Bool_t value) { fDebugFlag = value; } //!< debug flag virtual void SetNmuons(Int_t value) { fNmuons = value; } //!< number of muons
virtual void SetBfield(Double_t value) { fBfield = value; } //!< sets magnetic field (T) virtual void SetDebugFlag(Bool_t value) { fDebugFlag = value; } //!< debug flag
virtual void SetMuPrecFreq12(Double_t value) { fMuPrecFreq12 = value; } //!< sets Mu transition frequency (MHz) virtual void SetBfield(Double_t value) { fBfield = value; } //!< sets magnetic field (T)
virtual void SetMuPrecFreq34(Double_t value) { fMuPrecFreq34 = value; } //!< sets Mu transition frequency (MHz) virtual void SetMuPrecFreq12(Double_t value) { fMuPrecFreq12 = value; } //!< sets Mu transition frequency (MHz)
virtual void SetMuPrecFreq23(Double_t value) { fMuPrecFreq23 = value; } //!< sets Mu transition frequency (MHz) virtual void SetMuPrecFreq34(Double_t value) { fMuPrecFreq34 = value; } //!< sets Mu transition frequency (MHz)
virtual void SetMuPrecFreq14(Double_t value) { fMuPrecFreq14 = value; } //!< sets Mu transition frequency (MHz) virtual void SetMuPrecFreq23(Double_t value) { fMuPrecFreq23 = value; } //!< sets Mu transition frequency (MHz)
virtual void SetCaptureRate(Double_t value){ fCaptureRate = value; } //!< sets Mu+ electron capture rate (MHz) virtual void SetMuPrecFreq14(Double_t value) { fMuPrecFreq14 = value; } //!< sets Mu transition frequency (MHz)
virtual void SetIonizationRate(Double_t value){ fIonizationRate = value; } //!< sets Mu0 ionization rate (MHz) virtual void SetCaptureRate(Double_t value){ fCaptureRate = value; } //!< sets Mu+ electron capture rate (MHz)
virtual void SetDecayAsymmetry(Double_t value){ fAsymmetry = value; } //!< muon decay asymmetry virtual void SetIonizationRate(Double_t value){ fIonizationRate = value; } //!< sets Mu0 ionization rate (MHz)
virtual void SetMuFraction(Double_t value){ fMuFraction = value; } //!< Muonium fraction virtual void SetSpinFlipRate(Double_t value){ fSpinFlipRate = value; } //!< sets Mu0 spin flip rate (MHz)
virtual void SetMuFractionState12(Double_t value){ fMuFractionState12 = value; } virtual void SetDecayAsymmetry(Double_t value){ fAsymmetry = value; } //!< muon decay asymmetry
virtual void SetMuFractionState23(Double_t value){ fMuFractionState23 = value; } virtual void SetMuFraction(Double_t value){ fMuFraction = value; } //!< Muonium fraction
virtual void SetMuFractionState12(Double_t value){ fMuFractionState12 = value; }
virtual Bool_t IsValid() { return fValid; } virtual void SetMuFractionState23(Double_t value){ fMuFractionState23 = value; }
virtual void SetSeed(UInt_t seed);
virtual Bool_t IsValid() { return fValid; }
virtual Double_t GetBfield() { return fBfield; } //!< returns the magnetic field (T) virtual void SetSeed(UInt_t seed);
virtual Double_t GetCaptureRate() { return fCaptureRate; } //!< returns Mu+ electron capture rate (MHz)
virtual Double_t GetIonizationRate() { return fIonizationRate; } //!< returns Mu0 ionization rate (MHz) virtual Double_t GetBfield() { return fBfield; } //!< returns the magnetic field (T)
virtual void Run(TH1F *histoForward, TH1F *histoBackward); virtual Double_t GetCaptureRate() { return fCaptureRate; } //!< returns Mu+ electron capture rate (MHz)
virtual Double_t GetIonizationRate() { return fIonizationRate; } //!< returns Mu0 ionization rate (MHz)
private: virtual void Run(TH1F *histoForward, TH1F *histoBackward);
Bool_t fValid;
TRandom2 *fRandom; private:
Bool_t fValid;
Double_t fBfield; //!< magnetic field (T) TRandom2 *fRandom;
Double_t fMuPrecFreq12; //!< Mu transition frequency 12 (MHz)
Double_t fMuPrecFreq34; //!< Mu transition frequency 34 (MHz) Double_t fBfield; //!< magnetic field (T)
Double_t fMuPrecFreq23; //!< Mu transition frequency 23 (MHz) Double_t fMuPrecFreq12; //!< Mu transition frequency 12 (MHz)
Double_t fMuPrecFreq14; //!< Mu transition frequency 14 (MHz) Double_t fMuPrecFreq34; //!< Mu transition frequency 34 (MHz)
Double_t fMuonPrecFreq; //!< muon precession frequency (MHz) Double_t fMuPrecFreq23; //!< Mu transition frequency 23 (MHz)
Double_t fCaptureRate; //!< Mu+ electron capture rate (MHz) Double_t fMuPrecFreq14; //!< Mu transition frequency 14 (MHz)
Double_t fIonizationRate; //!< Mu0 ionization rate (MHz) Double_t fMuonPrecFreq; //!< muon precession frequency (MHz)
Double_t fInitialPhase; //!< initial muon spin phase Double_t fCaptureRate; //!< Mu+ electron capture rate (MHz)
Double_t fMuonDecayTime; //!< muon decay time (us) Double_t fIonizationRate; //!< Mu0 ionization rate (MHz)
Double_t fMuonPhase; //!< phase of muon spin Double_t fSpinFlipRate; //!< Mu0 spin-flip rate (MHz)
Double_t fAsymmetry; //!< muon decay asymmetry Double_t fInitialPhase; //!< initial muon spin phase
Double_t fMuFraction; //!< total Mu fraction [0,1] Double_t fMuonDecayTime; //!< muon decay time (us)
Double_t fMuFractionState12; //!< fraction of Mu in state 12, 34 Double_t fMuonPhase; //!< phase of muon spin
Double_t fMuFractionState23; //!< fraction of Mu in state 23, 14 Double_t fAsymmetry; //!< muon decay asymmetry
Int_t fNmuons; //!< number of muons to simulate Double_t fMuFraction; //!< total Mu fraction [0,1]
Bool_t fDebugFlag; //!< debug flag Double_t fMuFractionState12; //!< fraction of Mu in state 12, 34
Double_t fMuFractionState23; //!< fraction of Mu in state 23, 14
virtual Double_t NextEventTime(const Double_t &EventRate); Int_t fNmuons; //!< number of muons to simulate
// virtual Double_t PrecessionPhase(const Double_t &time, const Double_t &frequency); Bool_t fDebugFlag; //!< debug flag
virtual Double_t PrecessionPhase(const Double_t &time, const TString chargeState);
virtual void Event(const TString muonString); virtual Double_t NextEventTime(const Double_t &EventRate);
// virtual Double_t PrecessionPhase(const Double_t &time, const Double_t &frequency);
ClassDef(PSimulateMuTransition, 0) virtual Double_t PrecessionPhase(const Double_t &time, const TString chargeState);
}; virtual TComplex GTFunction(const Double_t &time); //!< transverse field polarization function of Mu0
virtual Double_t GTSpinFlip(const Double_t &time); //!< transverse field polarization function after spin-flip collisions
#endif // _PSIMULATEMUTRANSITION_H_ virtual void Event(const TString muonString);
ClassDef(PSimulateMuTransition, 0)
};
#endif // _PSIMULATEMUTRANSITION_H_

350
src/tests/MuTransition/runMuSimulation.C
View File

@ -1,148 +1,202 @@
/*************************************************************************** /***************************************************************************
runMuSimulation.C runMuSimulation.C
Author: Thomas Prokscha Author: Thomas Prokscha
Date: 25-Feb-2010 Date: 25-Feb-2010
$Id$ $Id$
***************************************************************************/ ***************************************************************************/
/*************************************************************************** /***************************************************************************
* Copyright (C) 2010 by Thomas Prokscha, Paul Scherrer Institut * * Copyright (C) 2010 by Thomas Prokscha, Paul Scherrer Institut *
* * * *
* This program is free software; you can redistribute it and/or modify * * This program is free software; you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by * * it under the terms of the GNU General Public License as published by *
* the Free Software Foundation; either version 2 of the License, or * * the Free Software Foundation; either version 2 of the License, or *
* (at your option) any later version. * * (at your option) any later version. *
* * * *
* This program is distributed in the hope that it will be useful, * * This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of * * but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. * * GNU General Public License for more details. *
* * * *
* You should have received a copy of the GNU General Public License * * You should have received a copy of the GNU General Public License *
* along with this program; if not, write to the * * along with this program; if not, write to the *
* Free Software Foundation, Inc., * * Free Software Foundation, Inc., *
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. * * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
***************************************************************************/ ***************************************************************************/
void runMuSimulation() #include "/apps/cern/root-git/include/TMusrRunHeader.h"
{ #define NDECAYHISTS 2
// load library
gSystem->Load("$ROOTSYS/lib/libPSimulateMuTransition"); void runMuSimulation()
{
// generate data // load library
TFolder *histosFolder; gSystem->Load("$ROOTSYS/lib/libPSimulateMuTransition");
TFolder *decayAnaModule;
TFolder *runInfo; // load TMusrRunHeader class if not already done during root startup
if (!TClass::GetDict("TMusrRunHeader")) {
histosFolder = gROOT->GetRootFolder()->AddFolder("histos", "Histograms"); gROOT->LoadMacro("$(ROOTSYS)/lib/libTMusrRunHeader.so");
gROOT->GetListOfBrowsables()->Add(histosFolder, "histos"); }
decayAnaModule = histosFolder->AddFolder("DecayAnaModule", "muSR decay histograms");
char titleStr[256];
//prepare to run simulation; here: isotropic Mu in Germanium TFolder *histosFolder;
UInt_t runNo = 9903; TFolder *decayAnaModule;
Double_t T = 300.; //temperature TFolder *gRunHeader;
Double_t capRate = 1.0;//*sqrt(T/200.); TString runTitle;
//assume that capture rate varies as sqrt(T), capRate = sigma*v*p , v ~ sqrt(T) TString histogramFileName;
Double_t ionRate; //assume Arrhenius behaviour ionRate = preFac*exp(-EA/kT) TObjArray Slist(0);
Double_t EA = 100.; //activation energy (meV) TMusrRunPhysicalQuantity prop;
ionRate = 2.9e7 * exp(-EA/(0.08625*T)); // Ge: 2.9*10^7MHz "attempt" frequency; 1K = 0.08625 meV
Double_t B = 100.; //field in G //prepare to run simulation; here: isotropic Mu in Germanium
Double_t Freq12 = 4463; //Mu freq of the 12 transition UInt_t runNo = 9900;
Double_t Freq34 = 4463; //Mu freq of the 34 transition Double_t T = 300.; //temperature
Double_t Freq23 = 4463; //Mu freq of the 23 transition Double_t spinFlipRate = 0.001; //if spinFlipRate > 0.001 only spin-flip processes will be simulated
Double_t Freq14 = 4463; //Mu freq of the 14 transition Double_t capRate = 1.0;//*sqrt(T/200.); //assume that capture rate varies as sqrt(T), capRate = sigma*v*p , v ~ sqrt(T)
Double_t MuFrac = 1.0; //total Mu fraction Double_t ionRate; //assume Arrhenius behaviour ionRate = preFac*exp(-EA/kT)
Double_t MuFrac12 = 0.5; //Mu in states 12 and 34 Double_t EA = 100.; //activation energy (meV)
Double_t MuFrac23 = 0.5; //Mu in states 23 and 14 ionRate = 2.9e7 * exp(-EA/(0.08625*T)); // Ge: 2.9*10^7MHz "attempt" frequency; 1K = 0.08625 meV
Int_t Nmuons = 1e7; //number of muons Double_t B = 100.; //field in G
Double_t Asym = 0.27; //muon decay asymmetry Double_t Bvar = 0.; //field variance
Double_t Freq12 = 134.858; //Mu freq of the 12 transition
// feed run info header Double_t Freq34 = 4328.142; //Mu freq of the 34 transition
TString tstr; Double_t Freq23 = 143.713; //Mu freq of the 23 transition
runInfo = gROOT->GetRootFolder()->AddFolder("RunInfo", "LEM RunInfo"); Double_t Freq14 = 4606.713; //Mu freq of the 14 transition
gROOT->GetListOfBrowsables()->Add(runInfo, "RunInfo"); Double_t MuFrac = 1.0; //total Mu fraction
header = new TLemRunHeader(); Double_t MuFrac12 = 2*0.266; //Mu in states 12 and 34
tstr = TString("0"); Double_t MuFrac23 = 2*0.234; //Mu in states 23 and 14
tstr += runNo; Int_t Nmuons = 1e6; //number of muons
tstr += TString(" - Mu-frac 1.0, Mu12 -4463MHz (0.5), Mu34 -4463MHz(0.5), T=300K/EA=100meV, Cap. 1.0MHz, 10mT"); Double_t Asym = 0.27; //muon decay asymmetry
header->SetRunTitle(tstr.Data()); histogramFileName = TString("0");
header->SetLemSetup("trivial"); histogramFileName += runNo;
header->SetRunNumber(runNo); histogramFileName += TString(".root");
header->SetStartTime(0);
header->SetStopTime(1); sprintf(titleStr,"- Mu-frac %3.2f, Mu12 %6.2f MHz(%3.2f), Mu23 %6.2f MHz(%3.2f), ionRate %8.2f MHz, capRate %6.2f MHz, SF rate %6.2f MHz, %5.0f G", MuFrac, Freq12, MuFrac12/2, Freq23, MuFrac23/2, ionRate, capRate, spinFlipRate, B);
header->SetModeratorHV(32.0, 0.01); runTitle = TString("0");
header->SetSampleHV(0.0, 0.01); runTitle += runNo;
header->SetImpEnergy(31.8); runTitle += TString(titleStr);
header->SetSampleTemperature(T, 0.001);
header->SetSampleBField(B, 0.1); PSimulateMuTransition *simulateMuTransition = new PSimulateMuTransition();
header->SetTimeResolution(1.); if (!simulateMuTransition->IsValid()){
header->SetNChannels(12001); cerr << endl << "**ERROR** while invoking PSimulateTransition" << endl;
header->SetNHist(2); return;
header->SetOffsetPPCHistograms(20); }
header->SetCuts("none"); simulateMuTransition->SetMuPrecFreq12(Freq12); // MHz
header->SetModerator("none"); simulateMuTransition->SetMuPrecFreq34(Freq34); // MHz
Double_t tt0[2] = {0., 0.}; simulateMuTransition->SetMuPrecFreq23(Freq23); // MHz
header->SetTimeZero(tt0); simulateMuTransition->SetMuPrecFreq14(Freq14); // MHz
runInfo->Add(header); //add header to RunInfo folder simulateMuTransition->SetMuFraction(MuFrac); // initial Mu fraction
simulateMuTransition->SetMuFractionState12(MuFrac12); // Mu in states 12, 34
TH1F *histo[4]; simulateMuTransition->SetMuFractionState23(MuFrac23); // Mu in states 23, 14
char str[128]; simulateMuTransition->SetBfield(B/10000.); // Tesla
for (UInt_t i=0; i<2; i++) { simulateMuTransition->SetCaptureRate(capRate); // MHz
sprintf(str, "hDecay0%d", (Int_t)i); simulateMuTransition->SetIonizationRate(ionRate); // MHz
histo[i] = new TH1F(str, str, 12001, -0.5, 12000.5); simulateMuTransition->SetSpinFlipRate(spinFlipRate); // MHz
sprintf(str, "hDecay2%d", (Int_t)i); simulateMuTransition->SetNmuons(Nmuons);
histo[i+2] = new TH1F(str, str, 12001, -0.5, 12000.5); simulateMuTransition->SetDecayAsymmetry(Asym);
} simulateMuTransition->SetDebugFlag(kFALSE); // to print time and phase during charge-changing cycle
PSimulateMuTransition *simulateMuTransition = new PSimulateMuTransition(); // feed run info header
if (!simulateMuTransition->IsValid()) { gRunHeader = gROOT->GetRootFolder()->AddFolder("RunHeader", "MuTransition Simulation Header Info");
cerr << endl << "**ERROR** while invoking PSimulateTransition" << endl; gROOT->GetListOfBrowsables()->Add(gRunHeader, "RunHeader");
return; // header = new TLemRunHeader();
} header = new TMusrRunHeader(true);
header->FillFolder(gRunHeader);
simulateMuTransition->SetMuPrecFreq12(Freq12); // MHz gRunHeader->Add(&Slist);
simulateMuTransition->SetMuPrecFreq34(Freq34); // MHz Slist.SetName("RunSummary");
simulateMuTransition->SetMuPrecFreq23(Freq23); // MHz
simulateMuTransition->SetMuPrecFreq14(Freq14); // MHz header->Set("RunInfo/Generic Validator URL", "http://lmu.web.psi.ch/facilities/software/MusrRoot/validation/MusrRoot.xsd");
simulateMuTransition->SetMuFraction(MuFrac); // initial Mu fraction header->Set("RunInfo/Specific Validator URL", "http://lmu.web.psi.ch/facilities/software/MusrRoot/validation/MusrRootLEM.xsd");
simulateMuTransition->SetMuFractionState12(MuFrac12); // Mu in states 12, 34 header->Set("RunInfo/Generator", "runMuSimulation");
simulateMuTransition->SetMuFractionState23(MuFrac23); // Mu in states 23, 14
simulateMuTransition->SetBfield(B/10000.); // Tesla
simulateMuTransition->SetCaptureRate(capRate); // MHz header->Set("RunInfo/File Name", histogramFileName.Data());
simulateMuTransition->SetIonizationRate(ionRate); // MHz header->Set("RunInfo/Run Title", runTitle.Data());
simulateMuTransition->SetNmuons(Nmuons); header->Set("RunInfo/Run Number", (Int_t) runNo);
simulateMuTransition->SetDecayAsymmetry(Asym); header->Set("RunInfo/Run Start Time", "0");
simulateMuTransition->SetDebugFlag(kFALSE); // to print time and phase during charge-changing cycle header->Set("RunInfo/Run Stop Time", "1");
prop.Set("Run Duration", 1.0, "sec");
simulateMuTransition->PrintSettings(); header->Set("RunInfo/Run Duration", prop);
header->Set("RunInfo/Laboratory", "PSI");
simulateMuTransition->Run(histo[0], histo[1]); header->Set("RunInfo/Instrument", "MC-Simulation");
prop.Set("Muon Beam Momentum", 0.0, "MeV/c");
for (UInt_t i=0; i<4; i++) header->Set("RunInfo/Muon Beam Momentum", prop);
decayAnaModule->Add(histo[i]); header->Set("RunInfo/Muon Species", "positive muon and muonium");
header->Set("RunInfo/Muon Source", "Simulation");
// write file header->Set("RunInfo/Setup", "Monte-Carlo setup");
tstr = TString("0"); header->Set("RunInfo/Comment", "Testing effect of charge-exchange or Mu0 spin flip processes on uSR signal");
tstr += runNo; header->Set("RunInfo/Sample Name", "Monte-Carlo");
tstr += TString(".root"); header->Set("RunInfo/Sample Temperature", 300);
TFile *fout = new TFile(tstr.Data(), "RECREATE", "Midas Fake Histograms"); prop.Set("Sample Magnetic Field", MRH_UNDEFINED, B, Bvar, "G");
if (fout == 0) { header->Set("RunInfo/Sample Magnetic Field", prop);
cout << endl << "**ERROR** Couldn't create ROOT file"; header->Set("RunInfo/No of Histos", 2);
cout << endl << endl; prop.Set("Time Resolution", 1.0, "ns", "Simulation");
exit(0); header->Set("RunInfo/Time Resolution", prop);
}
header->Set("DetectorInfo/Detector001/Name", "hDecay001");
fout->cd(); header->Set("DetectorInfo/Detector001/Histo Number", 1);
runInfo->Write(); header->Set("DetectorInfo/Detector001/Histo Length", 12001);
histosFolder->Write(); header->Set("DetectorInfo/Detector001/Time Zero Bin", 0);
fout->Close(); header->Set("DetectorInfo/Detector001/First Good Bin", 1);
cout << "Histograms written to " << tstr.Data() << endl; header->Set("DetectorInfo/Detector001/Last Good Bin", 12001);
delete fout;
header->Set("DetectorInfo/Detector002/Name", "hDecay002");
delete [] histo; header->Set("DetectorInfo/Detector002/Histo Number", 1);
} header->Set("DetectorInfo/Detector002/Histo Length", 12001);
header->Set("DetectorInfo/Detector002/Time Zero Bin", 0);
header->Set("DetectorInfo/Detector002/First Good Bin", 1);
header->Set("DetectorInfo/Detector002/Last Good Bin", 12001);
// simulation parameters
header->Set("Simulation/Mu Precession frequency 12", Freq12);
header->Set("Simulation/Mu Precession frequency 34", Freq34);
header->Set("Simulation/Mu Precession frequency 23", Freq23);
header->Set("Simulation/Mu Precession frequency 14", Freq14);
header->Set("Simulation/Mu Fraction", MuFrac);
header->Set("Simulation/Mu Fraction 12", MuFrac12);
header->Set("Simulation/Mu Fraction 23", MuFrac23);
header->Set("Simulation/Mu+ Capture Rate", capRate);
header->Set("Simulation/Mu0 Ionization Rate", ionRate);
header->Set("Simulation/Mu0 Spin Flip Rate", spinFlipRate);
header->Set("Simulation/Number of Muons", Nmuons);
header->Set("Simulation/Decay Asymmetry", Asym);
histosFolder = gROOT->GetRootFolder()->AddFolder("histos", "Histograms");
gROOT->GetListOfBrowsables()->Add(histosFolder, "histos");
decayAnaModule = histosFolder->AddFolder("DecayAnaModule", "muSR decay histograms");
TH1F *histo[NDECAYHISTS];
char str[128];
for (UInt_t i=0; i<NDECAYHISTS; i++) {
sprintf(str, "hDecay00%d", (Int_t)i+1);
histo[i] = new TH1F(str, str, 12001, -0.5, 12000.5);
}
for (i=0; i<NDECAYHISTS; i++)
decayAnaModule->Add(histo[i]);
// run simulation
simulateMuTransition->PrintSettings();
simulateMuTransition->Run(histo[0], histo[1]);
// write file
TFile *fout = new TFile(histogramFileName.Data(), "RECREATE", "Midas MC Histograms");
if (fout == 0) {
cout << endl << "**ERROR** Couldn't create ROOT file";
cout << endl << endl;
exit(0);
}
fout->cd();
header->FillFolder(gRunHeader);
gRunHeader->Write();
histosFolder->Write();
fout->Close();
cout << "Histograms written to " << histogramFileName.Data() << endl;
delete fout;
delete [] histo;
}