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

Merge branch 'master' into root6

Browse Source
This commit is contained in:
nemu
2016-04-26 12:58:50 +02:00
parent 0eb017c946 af7b729a5a
commit 95ac103605
42 changed files with 886 additions and 617 deletions
Download Patch File Download Diff File Expand all files Collapse all files

5
src/classes/PMsrHandler.cpp
View File

@ -5866,9 +5866,10 @@ void PMsrHandler::CheckMaxLikelihood()
{
if (!fStatistic.fChisq) {
for (UInt_t i=0; i<fRuns.size(); i++) {
if ((fRuns[i].GetFitType() != MSR_FITTYPE_SINGLE_HISTO) && (fGlobal.GetFitType() != MSR_FITTYPE_SINGLE_HISTO)) {
if ((fRuns[i].GetFitType() != MSR_FITTYPE_SINGLE_HISTO) && (fGlobal.GetFitType() != MSR_FITTYPE_SINGLE_HISTO) &&
(fRuns[i].GetFitType() != MSR_FITTYPE_MU_MINUS) && (fGlobal.GetFitType() != MSR_FITTYPE_MU_MINUS)) {
cerr << endl << ">> PMsrHandler::CheckMaxLikelihood: **WARNING**: Maximum Log Likelihood Fit is only implemented";
cerr << endl << ">> for Single Histogram Fit. Will fall back to Chi Square Fit.";
cerr << endl << ">> for Single Histogram and Mu Minus Fits. Will fall back to Chi Square Fit.";
cerr << endl << endl;
fStatistic.fChisq = true;
break;

36
src/classes/PRunAsymmetry.cpp
View File

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

64
src/classes/PRunMuMinus.cpp
View File

@ -52,6 +52,9 @@ PRunMuMinus::PRunMuMinus() : PRunBase()
fGoodBins[0] = -1;
fGoodBins[1] = -1;
fStartTimeBin = -1;
fEndTimeBin = -1;
fHandleTag = kEmpty;
}
@ -87,6 +90,9 @@ PRunMuMinus::PRunMuMinus(PMsrHandler *msrInfo, PRunDataHandler *rawData, UInt_t
fGoodBins[0] = -1;
fGoodBins[1] = -1;
fStartTimeBin = -1;
fEndTimeBin = -1;
if (!PrepareData()) {
cerr << endl << ">> PRunMuMinus::PRunMuMinus: **SEVERE ERROR**: Couldn't prepare data for fitting!";
cerr << endl << ">> This is very bad :-(, will quit ...";
@ -130,15 +136,7 @@ Double_t PRunMuMinus::CalcChiSquare(const std::vector<Double_t>& par)
// calculate chi square
Double_t time(1.0);
Int_t i, N(static_cast<Int_t>(fData.GetValue()->size()));
// In order not to have an IF in the next loop, determine the start and end bins for the fit range now
Int_t startTimeBin = static_cast<Int_t>(ceil((fFitStartTime - fData.GetDataTimeStart())/fData.GetDataTimeStep()));
if (startTimeBin < 0)
startTimeBin = 0;
Int_t endTimeBin = static_cast<Int_t>(floor((fFitEndTime - fData.GetDataTimeStart())/fData.GetDataTimeStep())) + 1;
if (endTimeBin > N)
endTimeBin = N;
Int_t i;
// Calculate the theory function once to ensure one function evaluation for the current set of parameters.
// This is needed for the LF and user functions where some non-thread-save calculations only need to be calculated once
@ -147,12 +145,12 @@ Double_t PRunMuMinus::CalcChiSquare(const std::vector<Double_t>& par)
time = fTheory->Func(time, par, fFuncValues);
#ifdef HAVE_GOMP
Int_t chunk = (endTimeBin - startTimeBin)/omp_get_num_procs();
Int_t chunk = (fEndTimeBin - fStartTimeBin)/omp_get_num_procs();
if (chunk < 10)
chunk = 10;
#pragma omp parallel for default(shared) private(i,time,diff) schedule(dynamic,chunk) reduction(+:chisq)
#endif
for (i=startTimeBin; i < endTimeBin; ++i) {
for (i=fStartTimeBin; i < fEndTimeBin; ++i) {
time = fData.GetDataTimeStart() + (Double_t)i*fData.GetDataTimeStep();
diff = fData.GetValue()->at(i) - fTheory->Func(time, par, fFuncValues);
chisq += diff*diff / (fData.GetError()->at(i)*fData.GetError()->at(i));
@ -186,15 +184,7 @@ Double_t PRunMuMinus::CalcChiSquareExpected(const std::vector<Double_t>& par)
// calculate chi square
Double_t time(1.0);
Int_t i, N(static_cast<Int_t>(fData.GetValue()->size()));
// In order not to have an IF in the next loop, determine the start and end bins for the fit range now
Int_t startTimeBin = static_cast<Int_t>(ceil((fFitStartTime - fData.GetDataTimeStart())/fData.GetDataTimeStep()));
if (startTimeBin < 0)
startTimeBin = 0;
Int_t endTimeBin = static_cast<Int_t>(floor((fFitEndTime - fData.GetDataTimeStart())/fData.GetDataTimeStep())) + 1;
if (endTimeBin > N)
endTimeBin = N;
Int_t i;
// Calculate the theory function once to ensure one function evaluation for the current set of parameters.
// This is needed for the LF and user functions where some non-thread-save calculations only need to be calculated once
@ -203,12 +193,12 @@ Double_t PRunMuMinus::CalcChiSquareExpected(const std::vector<Double_t>& par)
time = fTheory->Func(time, par, fFuncValues);
#ifdef HAVE_GOMP
Int_t chunk = (endTimeBin - startTimeBin)/omp_get_num_procs();
Int_t chunk = (fEndTimeBin - fStartTimeBin)/omp_get_num_procs();
if (chunk < 10)
chunk = 10;
#pragma omp parallel for default(shared) private(i,time,diff) schedule(dynamic,chunk) reduction(+:chisq)
#endif
for (i=startTimeBin; i < endTimeBin; ++i) {
for (i=fStartTimeBin; i < fEndTimeBin; ++i) {
time = fData.GetDataTimeStart() + (Double_t)i*fData.GetDataTimeStep();
theo = fTheory->Func(time, par, fFuncValues);
diff = fData.GetValue()->at(i) - theo;
@ -243,15 +233,7 @@ Double_t PRunMuMinus::CalcMaxLikelihood(const std::vector<Double_t>& par)
Double_t theo;
Double_t data;
Double_t time(1.0);
Int_t i, N(static_cast<Int_t>(fData.GetValue()->size()));
// In order not to have an IF in the next loop, determine the start and end bins for the fit range now
Int_t startTimeBin = static_cast<Int_t>(ceil((fFitStartTime - fData.GetDataTimeStart())/fData.GetDataTimeStep()));
if (startTimeBin < 0)
startTimeBin = 0;
Int_t endTimeBin = static_cast<Int_t>(floor((fFitEndTime - fData.GetDataTimeStart())/fData.GetDataTimeStep())) + 1;
if (endTimeBin > N)
endTimeBin = N;
Int_t i;
// Calculate the theory function once to ensure one function evaluation for the current set of parameters.
// This is needed for the LF and user functions where some non-thread-save calculations only need to be calculated once
@ -260,12 +242,12 @@ Double_t PRunMuMinus::CalcMaxLikelihood(const std::vector<Double_t>& par)
time = fTheory->Func(time, par, fFuncValues);
#ifdef HAVE_GOMP
Int_t chunk = (endTimeBin - startTimeBin)/omp_get_num_procs();
Int_t chunk = (fEndTimeBin - fStartTimeBin)/omp_get_num_procs();
if (chunk < 10)
chunk = 10;
#pragma omp parallel for default(shared) private(i,time,theo,data) schedule(dynamic,chunk) reduction(-:mllh)
#endif
for (i=startTimeBin; i < endTimeBin; ++i) {
for (i=fStartTimeBin; i < fEndTimeBin; ++i) {
time = fData.GetDataTimeStart() + (Double_t)i*fData.GetDataTimeStep();
// calculate theory for the given parameter set
theo = fTheory->Func(time, par, fFuncValues);
@ -401,15 +383,15 @@ void PRunMuMinus::SetFitRangeBin(const TString fitRange)
void PRunMuMinus::CalcNoOfFitBins()
{
// In order not having to loop over all bins and to stay consistent with the chisq method, calculate the start and end bins explicitly
Int_t startTimeBin = static_cast<Int_t>(ceil((fFitStartTime - fData.GetDataTimeStart())/fData.GetDataTimeStep()));
if (startTimeBin < 0)
startTimeBin = 0;
Int_t endTimeBin = static_cast<Int_t>(floor((fFitEndTime - fData.GetDataTimeStart())/fData.GetDataTimeStep())) + 1;
if (endTimeBin > static_cast<Int_t>(fData.GetValue()->size()))
endTimeBin = fData.GetValue()->size();
fStartTimeBin = static_cast<Int_t>(ceil((fFitStartTime - fData.GetDataTimeStart())/fData.GetDataTimeStep()));
if (fStartTimeBin < 0)
fStartTimeBin = 0;
fEndTimeBin = static_cast<Int_t>(floor((fFitEndTime - fData.GetDataTimeStart())/fData.GetDataTimeStep())) + 1;
if (fEndTimeBin > static_cast<Int_t>(fData.GetValue()->size()))
fEndTimeBin = fData.GetValue()->size();
if (endTimeBin > startTimeBin)
fNoOfFitBins = endTimeBin - startTimeBin;
if (fEndTimeBin > fStartTimeBin)
fNoOfFitBins = fEndTimeBin - fStartTimeBin;
else
fNoOfFitBins = 0;
}

72
src/classes/PRunSingleHisto.cpp
View File

@ -64,6 +64,9 @@ PRunSingleHisto::PRunSingleHisto() : PRunBase()
// the fit range can be changed in the command block, these variables need to be accessible
fGoodBins[0] = -1;
fGoodBins[1] = -1;
fStartTimeBin = -1;
fEndTimeBin = -1;
}
//--------------------------------------------------------------------------
@ -81,6 +84,7 @@ PRunSingleHisto::PRunSingleHisto(PMsrHandler *msrInfo, PRunDataHandler *rawData,
{
fScaleN0AndBkg = IsScaleN0AndBkg();
fNoOfFitBins = 0;
fBackground = 0;
fPacking = fRunInfo->GetPacking();
if (fPacking == -1) { // i.e. packing is NOT given in the RUN-block, it must be given in the GLOBAL-block
@ -99,6 +103,9 @@ PRunSingleHisto::PRunSingleHisto(PMsrHandler *msrInfo, PRunDataHandler *rawData,
fGoodBins[0] = -1;
fGoodBins[1] = -1;
fStartTimeBin = -1;
fEndTimeBin = -1;
if (!PrepareData()) {
cerr << endl << ">> PRunSingleHisto::PRunSingleHisto: **SEVERE ERROR**: Couldn't prepare data for fitting!";
cerr << endl << ">> This is very bad :-(, will quit ...";
@ -173,15 +180,7 @@ Double_t PRunSingleHisto::CalcChiSquare(const std::vector<Double_t>& par)
// calculate chi square
Double_t time(1.0);
Int_t i, N(static_cast<Int_t>(fData.GetValue()->size()));
// In order not to have an IF in the next loop, determine the start and end bins for the fit range now
Int_t startTimeBin = static_cast<Int_t>(ceil((fFitStartTime - fData.GetDataTimeStart())/fData.GetDataTimeStep()));
if (startTimeBin < 0)
startTimeBin = 0;
Int_t endTimeBin = static_cast<Int_t>(floor((fFitEndTime - fData.GetDataTimeStart())/fData.GetDataTimeStep())) + 1;
if (endTimeBin > N)
endTimeBin = N;
Int_t i;
// Calculate the theory function once to ensure one function evaluation for the current set of parameters.
// This is needed for the LF and user functions where some non-thread-save calculations only need to be calculated once
@ -190,12 +189,12 @@ Double_t PRunSingleHisto::CalcChiSquare(const std::vector<Double_t>& par)
time = fTheory->Func(time, par, fFuncValues);
#ifdef HAVE_GOMP
Int_t chunk = (endTimeBin - startTimeBin)/omp_get_num_procs();
Int_t chunk = (fEndTimeBin - fStartTimeBin)/omp_get_num_procs();
if (chunk < 10)
chunk = 10;
#pragma omp parallel for default(shared) private(i,time,diff) schedule(dynamic,chunk) reduction(+:chisq)
#endif
for (i=startTimeBin; i < endTimeBin; ++i) {
for (i=fStartTimeBin; i<fEndTimeBin; ++i) {
time = fData.GetDataTimeStart() + (Double_t)i*fData.GetDataTimeStep();
diff = fData.GetValue()->at(i) -
(N0*TMath::Exp(-time/tau)*(1.0+fTheory->Func(time, par, fFuncValues))+bkg);
@ -266,15 +265,7 @@ Double_t PRunSingleHisto::CalcChiSquareExpected(const std::vector<Double_t>& par
// calculate chi square
Double_t time(1.0);
Int_t i, N(static_cast<Int_t>(fData.GetValue()->size()));
// In order not to have an IF in the next loop, determine the start and end bins for the fit range now
Int_t startTimeBin = static_cast<Int_t>(ceil((fFitStartTime - fData.GetDataTimeStart())/fData.GetDataTimeStep()));
if (startTimeBin < 0)
startTimeBin = 0;
Int_t endTimeBin = static_cast<Int_t>(floor((fFitEndTime - fData.GetDataTimeStart())/fData.GetDataTimeStep())) + 1;
if (endTimeBin > N)
endTimeBin = N;
Int_t i;
// Calculate the theory function once to ensure one function evaluation for the current set of parameters.
// This is needed for the LF and user functions where some non-thread-save calculations only need to be calculated once
@ -283,12 +274,12 @@ Double_t PRunSingleHisto::CalcChiSquareExpected(const std::vector<Double_t>& par
time = fTheory->Func(time, par, fFuncValues);
#ifdef HAVE_GOMP
Int_t chunk = (endTimeBin - startTimeBin)/omp_get_num_procs();
Int_t chunk = (fEndTimeBin - fStartTimeBin)/omp_get_num_procs();
if (chunk < 10)
chunk = 10;
#pragma omp parallel for default(shared) private(i,time,diff) schedule(dynamic,chunk) reduction(+:chisq)
#endif
for (i=startTimeBin; i < endTimeBin; ++i) {
for (i=fStartTimeBin; i < fEndTimeBin; ++i) {
time = fData.GetDataTimeStart() + (Double_t)i*fData.GetDataTimeStep();
theo = N0*TMath::Exp(-time/tau)*(1.0+fTheory->Func(time, par, fFuncValues))+bkg;
diff = fData.GetValue()->at(i) - theo;
@ -359,7 +350,7 @@ Double_t PRunSingleHisto::CalcMaxLikelihood(const std::vector<Double_t>& par)
Double_t theo;
Double_t data;
Double_t time(1.0);
Int_t i, N(static_cast<Int_t>(fData.GetValue()->size()));
Int_t i;
// norm is needed since there is no simple scaling like in chisq case to get the correct Max.Log.Likelihood value when normlizing N(t) to 1/ns
Double_t normalizer = 1.0;
@ -367,14 +358,6 @@ Double_t PRunSingleHisto::CalcMaxLikelihood(const std::vector<Double_t>& par)
if (fScaleN0AndBkg)
normalizer = fPacking * (fTimeResolution * 1.0e3);
// In order not to have an IF in the next loop, determine the start and end bins for the fit range now
Int_t startTimeBin = static_cast<Int_t>(ceil((fFitStartTime - fData.GetDataTimeStart())/fData.GetDataTimeStep()));
if (startTimeBin < 0)
startTimeBin = 0;
Int_t endTimeBin = static_cast<Int_t>(floor((fFitEndTime - fData.GetDataTimeStart())/fData.GetDataTimeStep())) + 1;
if (endTimeBin > N)
endTimeBin = N;
// Calculate the theory function once to ensure one function evaluation for the current set of parameters.
// This is needed for the LF and user functions where some non-thread-save calculations only need to be calculated once
// for a given set of parameters---which should be done outside of the parallelized loop.
@ -382,12 +365,12 @@ Double_t PRunSingleHisto::CalcMaxLikelihood(const std::vector<Double_t>& par)
time = fTheory->Func(time, par, fFuncValues);
#ifdef HAVE_GOMP
Int_t chunk = (endTimeBin - startTimeBin)/omp_get_num_procs();
Int_t chunk = (fEndTimeBin - fStartTimeBin)/omp_get_num_procs();
if (chunk < 10)
chunk = 10;
#pragma omp parallel for default(shared) private(i,time,theo,data) schedule(dynamic,chunk) reduction(-:mllh)
#endif
for (i=startTimeBin; i<endTimeBin; ++i) {
for (i=fStartTimeBin; i<fEndTimeBin; ++i) {
time = fData.GetDataTimeStart() + (Double_t)i*fData.GetDataTimeStep();
// calculate theory for the given parameter set
theo = N0*TMath::Exp(-time/tau)*(1.0+fTheory->Func(time, par, fFuncValues))+bkg;
@ -588,15 +571,15 @@ void PRunSingleHisto::SetFitRangeBin(const TString fitRange)
void PRunSingleHisto::CalcNoOfFitBins()
{
// In order not having to loop over all bins and to stay consistent with the chisq method, calculate the start and end bins explicitly
Int_t startTimeBin = static_cast<Int_t>(ceil((fFitStartTime - fData.GetDataTimeStart())/fData.GetDataTimeStep()));
if (startTimeBin < 0)
startTimeBin = 0;
Int_t endTimeBin = static_cast<Int_t>(floor((fFitEndTime - fData.GetDataTimeStart())/fData.GetDataTimeStep())) + 1;
if (endTimeBin > static_cast<Int_t>(fData.GetValue()->size()))
endTimeBin = fData.GetValue()->size();
fStartTimeBin = static_cast<Int_t>(ceil((fFitStartTime - fData.GetDataTimeStart())/fData.GetDataTimeStep()));
if (fStartTimeBin < 0)
fStartTimeBin = 0;
fEndTimeBin = static_cast<Int_t>(floor((fFitEndTime - fData.GetDataTimeStart())/fData.GetDataTimeStep())) + 1;
if (fEndTimeBin > static_cast<Int_t>(fData.GetValue()->size()))
fEndTimeBin = fData.GetValue()->size();
if (endTimeBin > startTimeBin)
fNoOfFitBins = endTimeBin - startTimeBin;
if (fEndTimeBin > fStartTimeBin)
fNoOfFitBins = fEndTimeBin - fStartTimeBin;
else
fNoOfFitBins = 0;
}
@ -1555,7 +1538,7 @@ void PRunSingleHisto::EstimateN0()
if (paramNo > 10000) // i.e. fun or map
return;
// still missing: set this value in the parameters
// get the parameters
PMsrParamList *param = fMsrInfo->GetMsrParamList();
assert(param);
@ -1564,6 +1547,11 @@ void PRunSingleHisto::EstimateN0()
return;
}
// check if N0 is fixed. If this is the case, do NOT estimate N0
if (param->at(paramNo-1).fStep == 0.0) // N0 parameter fixed
return;
// check that 'backgr.fit' in the msr-file run block is indeed a parameter number.
// in case it is a function, nothing will be done.
Int_t paramNoBkg = fRunInfo->GetBkgFitParamNo();

53
src/classes/PTheory.cpp
View File

@ -507,6 +507,10 @@ Double_t PTheory::Func(register Double_t t, const PDoubleVector& paramValues, co
return Polynom(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
fAdd->Func(t, paramValues, funcValues);
break;
case THEORY_MU_MINUS_EXP:
return MuMinusExpTF(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
fAdd->Func(t, paramValues, funcValues);
break;
case THEORY_USER_FCN:
return UserFcn(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
fAdd->Func(t, paramValues, funcValues);
@ -599,6 +603,9 @@ Double_t PTheory::Func(register Double_t t, const PDoubleVector& paramValues, co
case THEORY_DYNAMIC_TF_NK:
return DynamicNKTF (t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
break;
case THEORY_MU_MINUS_EXP:
return MuMinusExpTF(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
break;
case THEORY_POLYNOM:
return Polynom(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
break;
@ -695,6 +702,9 @@ Double_t PTheory::Func(register Double_t t, const PDoubleVector& paramValues, co
case THEORY_DYNAMIC_TF_NK:
return DynamicNKTF (t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
break;
case THEORY_MU_MINUS_EXP:
return MuMinusExpTF(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
break;
case THEORY_POLYNOM:
return Polynom(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
break;
@ -789,6 +799,9 @@ Double_t PTheory::Func(register Double_t t, const PDoubleVector& paramValues, co
case THEORY_DYNAMIC_TF_NK:
return DynamicNKTF(t, paramValues, funcValues);
break;
case THEORY_MU_MINUS_EXP:
return MuMinusExpTF(t, paramValues, funcValues);
break;
case THEORY_POLYNOM:
return Polynom(t, paramValues, funcValues);
break;
@ -2942,6 +2955,46 @@ Double_t PTheory::GetDynKTLFValue(const Double_t t) const
return fDynLFFuncValue[idx]+df;
}
//--------------------------------------------------------------------------
/**
* <p> theory function: MuMinusExpTF
*
* \f[ = N_0 \exp(-t/tau) [1 + A \exp(-\lambda t) \cos(2\pi\nu t + \phi)] \f]
*
* <b>meaning of paramValues:</b> \f$t_{\rm shift}\f$, \f$N_0\f$, \f$\tau\f$, \f$A\f$, \f$\lambda\f$, \f$\phi\f$, \f$\nu\f$
*
* <b>return:</b> function value
*
* \param t time in \f$(\mu\mathrm{s})\f$, or x-axis value for non-muSR fit
* \param paramValues parameter values
* \param funcValues vector with the functions (i.e. functions of the parameters)
*/
Double_t PTheory::MuMinusExpTF(register Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const
{
// expected parameters: N0 tau A lambda phase frequency [tshift]
Double_t val[7];
assert(fParamNo.size() <= 7);
// check if FUNCTIONS are used
for (UInt_t i=0; i<fParamNo.size(); i++) {
if (fParamNo[i] < MSR_PARAM_FUN_OFFSET) { // parameter or resolved map
val[i] = paramValues[fParamNo[i]];
} else { // function
val[i] = funcValues[fParamNo[i]-MSR_PARAM_FUN_OFFSET];
}
}
Double_t tt;
if (fParamNo.size() == 6) // no tshift
tt = t;
else // tshift present
tt = t-val[6];
return val[0]*exp(-tt/val[1])*(1.0+val[2]*exp(-val[3]*tt)*cos(TWO_PI*val[5]*tt+DEG_TO_RAD*val[4]));
}
//--------------------------------------------------------------------------
// END
//--------------------------------------------------------------------------

7
src/include/PRunAsymmetry.h
View File

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

7
src/include/PRunMuMinus.h
View File

@ -51,6 +51,10 @@ class PRunMuMinus : public PRunBase
virtual void SetFitRangeBin(const TString fitRange);
virtual Int_t GetStartTimeBin() { return fStartTimeBin; }
virtual Int_t GetEndTimeBin() { return fEndTimeBin; }
virtual Int_t GetPacking() { return fPacking; }
protected:
virtual void CalcNoOfFitBins();
virtual Bool_t PrepareData();
@ -65,6 +69,9 @@ class PRunMuMinus : public PRunBase
PDoubleVector fForward; ///< forward histo data
Int_t fStartTimeBin; ///< bin at which the fit starts
Int_t fEndTimeBin; ///< bin at which the fit ends
virtual Bool_t GetProperT0(PRawRunData* runData, PMsrGlobalBlock *globalBlock, PUIntVector &histoNo);
virtual Bool_t GetProperDataRange();
virtual void GetProperFitRange(PMsrGlobalBlock *globalBlock);

14
src/include/PRunSingleHisto.h
View File

@ -51,6 +51,13 @@ class PRunSingleHisto : public PRunBase
virtual void SetFitRangeBin(const TString fitRange);
virtual Double_t GetBackground() { return fBackground; }
virtual Int_t GetStartTimeBin() { return fStartTimeBin; }
virtual Int_t GetEndTimeBin() { return fEndTimeBin; }
virtual Int_t GetPacking() { return fPacking; }
virtual Bool_t GetScaleN0AndBkg() { return fScaleN0AndBkg; }
protected:
virtual void CalcNoOfFitBins();
virtual Bool_t PrepareData();
@ -64,9 +71,12 @@ class PRunSingleHisto : public PRunBase
Double_t fBackground; ///< needed if background range is given (units: 1/bin)
Int_t fPacking; ///< packing for this particular run. Either given in the RUN- or GLOBAL-block.
Int_t fGoodBins[2]; ///< keep first/last good bins. 0=fgb, 1=lgb
Int_t fGoodBins[2]; ///< keep first/last good bins. 0=fgb, 1=lgb
PDoubleVector fForward; ///< forward histo data
PDoubleVector fForward; ///< forward histo data
Int_t fStartTimeBin; ///< bin at which the fit starts
Int_t fEndTimeBin; ///< bin at which the fit ends
virtual Bool_t GetProperT0(PRawRunData* runData, PMsrGlobalBlock *globalBlock, PUIntVector &histoNo);
virtual Bool_t GetProperDataRange();

12
src/include/PTheory.h
View File

@ -70,8 +70,9 @@
#define THEORY_STATIC_TF_NK 24
#define THEORY_DYNAMIC_ZF_NK 25
#define THEORY_DYNAMIC_TF_NK 26
#define THEORY_POLYNOM 27
#define THEORY_USER_FCN 28
#define THEORY_MU_MINUS_EXP 27
#define THEORY_POLYNOM 28
#define THEORY_USER_FCN 29
// function parameter tags, i.e. how many parameters has a specific function
// if there is a comment with a (tshift), the number of parameters is increased by one
@ -102,9 +103,10 @@
#define THEORY_PARAM_STATIC_TF_NK 4 // phase, frequency, damping D0, R_b=DGbG/D0 (tshift)
#define THEORY_PARAM_DYNAMIC_ZF_NK 3 // damping D0, R_b=DGbG/D0, nu_c (tshift)
#define THEORY_PARAM_DYNAMIC_TF_NK 5 // phase, frequency, damping D0, R_b=DGbG/D0, nu_c (tshift)
#define THEORY_PARAM_MU_MINUS_EXP 6 // N0, tau, A, damping, phase, frequency (tshift)
// number of available user functions
#define THEORY_MAX 29
#define THEORY_MAX 30
// maximal number of parameters. Needed in the contents of LF
#define THEORY_MAX_PARAM 10
@ -217,6 +219,9 @@ static PTheoDataBase fgTheoDataBase[THEORY_MAX] = {
{THEORY_DYNAMIC_TF_NK, THEORY_PARAM_DYNAMIC_TF_NK, false,
"dynamicNKTF", "dnktf", "(phase frequency damping_D0 R_b nu_c)", "(phase frequency damping_D0 R_b nu_c tshift)"},
{THEORY_MU_MINUS_EXP, THEORY_PARAM_MU_MINUS_EXP, false,
"muMinusExpTF", "mmsetf", "(N0 tau A lambda phase nu)", "(N0 tau A lambda phase nu tshift)"},
{THEORY_POLYNOM, 0, false,
"polynom", "p", "(tshift p0 p1 ... pn)", "(tshift p0 p1 ... pn)"},
@ -272,6 +277,7 @@ class PTheory
virtual Double_t StaticNKTF(register Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const;
virtual Double_t DynamicNKZF(register Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const;
virtual Double_t DynamicNKTF(register Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const;
virtual Double_t MuMinusExpTF(register Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const;
virtual Double_t Polynom(register Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const;
virtual Double_t UserFcn(register Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const;

7
src/musredit/PTextEdit.cpp
View File

@ -2186,6 +2186,9 @@ void PTextEdit::musrView()
cmd = str + " \"";
str = *fFilenames.find( currentEditor() );
int pos = str.lastIndexOf("/");
if (pos != -1)
str.remove(0, pos+1);
QString numStr;
numStr.setNum(fAdmin->getTimeout());
cmd += str + "\" --timeout " + numStr;
@ -2194,6 +2197,10 @@ void PTextEdit::musrView()
cmd += " &";
int status=system(cmd.toLatin1());
if (status != 0) {
cerr << "**WARNING** musrView: something went wrong ..." << endl;
}
}
//----------------------------------------------------------------------------------------------------

7
src/musredit_qt5/PTextEdit.cpp
View File

@ -2184,6 +2184,9 @@ void PTextEdit::musrView()
cmd = str + " \"";
str = *fFilenames.find( currentEditor() );
int pos = str.lastIndexOf("/");
if (pos != -1)
str.remove(0, pos+1);
QString numStr;
numStr.setNum(fAdmin->getTimeout());
cmd += str + "\" --timeout " + numStr;
@ -2192,6 +2195,10 @@ void PTextEdit::musrView()
cmd += " &";
int status=system(cmd.toLatin1());
if (status != 0) {
cerr << "**WARNING** musrView: something went wrong ..." << endl;
}
}
//----------------------------------------------------------------------------------------------------

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

@ -3,9 +3,7 @@
PSimulateMuTransition.cpp
Author: Thomas Prokscha
Date: 25-Feb-2010
$Id$
Date: 25-Feb-2010, 14-Apr-2016
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
@ -13,26 +11,28 @@
analyze the simulated data.
Description:
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
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
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".
These frequencies and the corresponding probabilities ("SetMuFractionState12" for
transitions 12 and 34, "SetMuFractionState23" for states 23 and 14) can be calculated
Root class to simulate muon spin polarization under successive Mu+/Mu0 charge-exchange
or Mu0 spin-flip 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
muonium use the complex expression of
equation (4) in the paper of M. Senba, J. Phys. B 23, 1545 (1990), or
equation (7) in the paper of M. Senba, J. Phys. B 24, 3531 (1991);
note that PxMu(t) 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 transitions 12 and 34, "SetMuFractionState23" for states 23 and 14) can be calculated
for a given field with the root macro AnisotropicMu.C
Parameters:
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
3) total Mu0 fraction
4) electron-capture rate
5) Mu ionization rate
6) initial muon spin phase
7) total muon decay asymmetry
8) number of muon decays to be generated.
9) debug flag: if TRUE print capture/ionization events on screen
4) Mu+ electron-capture rate
5) Mu0 ionization rate
6) Mu0 spin-flip rate
7) initial muon spin phase
9) total muon decay asymmetry
9) number of muon decays to be generated.
10) debug flag: if TRUE print capture/ionization events on screen
Output:
Two histograms ("forward" and "backward") are written to a root file.
@ -43,10 +43,13 @@
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
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
muon spin phase by acos(Px(t_i)).
4) get the next electron capture time, continue until t_d is reached; accumulate muon spin
phase.
3) Determine next ionization time t_i; calculate Px(t_i) in Muonium; calculate the total
muon spin polarization Px(t_i)*Px(t_c).
4) get the next electron capture time, continue until t_d is reached, and calculate
the resulting polarization.
The Mu0 spin-flip processes are calculated in GTSpinFlip(), using eq. (17) of
M. Senba, J. Phys. B 24, 3531 (1991).
***************************************************************************/
@ -111,8 +114,10 @@ PSimulateMuTransition::PSimulateMuTransition(UInt_t seed)
fMuonDecayTime = 0.;
fAsymmetry = 0.27;
fMuFraction = 0.;
fMuFractionState12 = 0.;
fMuFractionState23 = 0.;
fMuFractionState12 = 0.25;
fMuFractionState34 = 0.25;
fMuFractionState23 = 0.25;
fMuFractionState14 = 0.25;
fDebugFlag = kFALSE;
}
@ -138,19 +143,27 @@ PSimulateMuTransition::~PSimulateMuTransition()
*/
void PSimulateMuTransition::PrintSettings() const
{
cout << endl << "Mu precession frequency 12 (MHz) = " << fMuPrecFreq12;
cout << endl << "Mu precession frequency 34 (MHz) = " << fMuPrecFreq34;
cout << endl << "Mu precession frequency 23 (MHz) = " << fMuPrecFreq23;
cout << endl << "Mu precession frequency 14 (MHz) = " << fMuPrecFreq14;
cout << endl << "Mu0 precession frequency 12 (MHz) = " << fMuPrecFreq12;
cout << endl << "Mu0 precession frequency 34 (MHz) = " << fMuPrecFreq34;
cout << endl << "Mu0 precession frequency 23 (MHz) = " << fMuPrecFreq23;
cout << endl << "Mu0 precession frequency 14 (MHz) = " << fMuPrecFreq14;
cout << endl << "Mu+ precession frequency (MHz) = " << fMuonGyroRatio * fBfield;
cout << endl << "B field (T) = " << fBfield;
cout << endl << "Mu+ electron capture rate (MHz) = " << fCaptureRate;
cout << endl << "Mu0 ionizatioan rate (MHz) = " << fIonizationRate;
cout << endl << "Mu0 spin-flip rate (MHz) = " << fSpinFlipRate;
cout << endl << "!!! Note: if spin-flip rate > 0.001 only spin-flip process is considered!!!";
if (fSpinFlipRate > 0.001)
cout << endl << "!!! Note: spin-flip rate > 0.001 only spin-flip processes are considered!!!";
else{
cout << endl << "!!! spin-flip rate <= 0.001: only charge-exchange cycles are considered!!!";
cout << endl << "!!! if spin-flip rate > 0.001, only spin-flip processes are considered!!!";
}
cout << endl << "Decay asymmetry = " << fAsymmetry;
cout << endl << "Muonium fraction = " << fMuFraction;
cout << endl << "Muonium fraction state12 = " << fMuFractionState12;
cout << endl << "Muonium fraction state34 = " << fMuFractionState34;
cout << endl << "Muonium fraction state23 = " << fMuFractionState23;
cout << endl << "Muonium fraction state14 = " << fMuFractionState14;
cout << endl << "Number of particles to simulate = " << fNmuons;
cout << endl << "Initial muon spin phase (degree) = " << fInitialPhase;
cout << endl << "Debug flag = " << fDebugFlag;
@ -181,23 +194,26 @@ void PSimulateMuTransition::SetSeed(UInt_t seed)
*/
void PSimulateMuTransition::Run(TH1F *histoForward, TH1F *histoBackward)
{
// Double_t muoniumPolX = 1.0; //polarization in x direction
Int_t i;
if (histoForward == 0 || histoBackward == 0)
return;
fMuonPrecFreq = fMuonGyroRatio * fBfield;
for (i = 0; i<fNmuons; i++){
fMuonPhase = TMath::TwoPi() * fInitialPhase/360.; // transform to radians
fMuonDecayTime = NextEventTime(fMuonDecayRate);
if (fSpinFlipRate > 0.001){// consider only Mu0 spin-flip in this case
fMuonPhase = TMath::ACos(GTSpinFlip(fMuonDecayTime));
fMuonPhase += TMath::ACos(GTSpinFlip(fMuonDecayTime));
}
else{
// initial muon state Mu+ or Mu0?
if (fRandom->Rndm() <= 1.-fMuFraction)
Event("Mu+");
fMuonPhase += TMath::ACos(Event("Mu+"));
else
Event("");
fMuonPhase += TMath::ACos(Event("Mu0"));
}
// fill 50% in "forward", and 50% in "backward" detector to get independent
// events in "forward" and "backward" histograms. This allows "normal" uSR
@ -233,28 +249,29 @@ Double_t PSimulateMuTransition::NextEventTime(const Double_t &EventRate)
//--------------------------------------------------------------------------
// Phase (private)
//--------------------------------------------------------------------------
/**
* <p>Determines phase of the muon spin
// /**
/* * <p>Determines phase of the muon spin
*
* \param time duration of precession (us);
* \param chargeState charge state of Mu ("Mu+" or "Mu0")
*/
Double_t PSimulateMuTransition::PrecessionPhase(const Double_t &time, const TString chargeState)
{
Double_t muonPhaseX;
Double_t muoniumPolX = 0;
if (chargeState == "Mu+")
muonPhaseX = TMath::TwoPi()*fMuonPrecFreq*time;
else if (chargeState == "Mu0"){
muoniumPolX = GTFunction(time).Re();
muonPhaseX = TMath::ACos(muoniumPolX);
}
else
muonPhaseX = 0.;
return muonPhaseX;
}
// Double_t PSimulateMuTransition::PrecessionPhase(const Double_t &time, const TString chargeState)
// {
// Double_t muonPhaseX;
// Double_t muoniumPolX = 0;
//
// if (chargeState == "Mu+")
// muonPhaseX = TMath::TwoPi()*fMuonPrecFreq*time;
// else if (chargeState == "Mu0"){
// muoniumPolX = GTFunction(time).Re();
// if (fDebugFlag) cout << "muoniumPolX = " << muoniumPolX << endl;
// muonPhaseX = TMath::ACos(muoniumPolX);
// }
// else
// muonPhaseX = 0.;
//
// return muonPhaseX;
// }
//--------------------------------------------------------------------------
// Mu0 transverse field polarization function (private)
@ -264,29 +281,24 @@ Double_t PSimulateMuTransition::PrecessionPhase(const Double_t &time, const TStr
*
* \param time (us);
*/
TComplex PSimulateMuTransition::GTFunction(const Double_t &time)
TComplex PSimulateMuTransition::GTFunction(const Double_t &time, const TString chargeState)
{
Double_t twoPi = TMath::TwoPi();
TComplex complexPol = 0;
complexPol =
0.5 * fMuFractionState12 *
(TComplex::Exp(TComplex::I()*twoPi*fMuPrecFreq12*time) +
TComplex::Exp(-TComplex::I()*twoPi*fMuPrecFreq34*time))
+
0.5 * fMuFractionState23 *
(TComplex::Exp(TComplex::I()*twoPi*fMuPrecFreq23*time) +
TComplex::Exp(TComplex::I()*twoPi*fMuPrecFreq14*time));
if (chargeState == "Mu+")
complexPol = TComplex::Exp(-TComplex::I()*twoPi*fMuonPrecFreq*time);
else{
complexPol =
(fMuFractionState12 * TComplex::Exp(TComplex::I()*twoPi*fMuPrecFreq12*time) +
fMuFractionState34 * TComplex::Exp(-TComplex::I()*twoPi*fMuPrecFreq34*time))
+
(fMuFractionState23 * TComplex::Exp(TComplex::I()*twoPi*fMuPrecFreq23*time) +
fMuFractionState14 * TComplex::Exp(TComplex::I()*twoPi*fMuPrecFreq14*time));
}
return complexPol;
// Double_t muoniumPolX = 0;
// muoniumPolX = 0.5 *
// (fMuFractionState12 * (TMath::Cos(twoPi*fMuPrecFreq12*time) + TMath::Cos(twoPi*fMuPrecFreq34*time)) +
// fMuFractionState23 * (TMath::Cos(twoPi*fMuPrecFreq23*time) + TMath::Cos(twoPi*fMuPrecFreq14*time)));
//
// return muoniumPolX;
}
//--------------------------------------------------------------------------
@ -308,18 +320,18 @@ Double_t PSimulateMuTransition::GTSpinFlip(const Double_t &time)
eventTime += NextEventTime(fSpinFlipRate);
if (eventTime >= time){
muoniumPolX = GTFunction(time).Re();
muoniumPolX = GTFunction(time, "Mu0").Re();
}
else{
while (eventTime < time){
eventDiffTime = eventTime - lastEventTime;
complexPolX = complexPolX * GTFunction(eventDiffTime);
complexPolX = complexPolX * GTFunction(eventDiffTime, "Mu0");
lastEventTime = eventTime;
eventTime += NextEventTime(fSpinFlipRate);
}
// calculate for the last collision
eventDiffTime = time - lastEventTime;
complexPolX = complexPolX * GTFunction(eventDiffTime);
complexPolX = complexPolX * GTFunction(eventDiffTime, "Mu0");
muoniumPolX = complexPolX.Re();
}
@ -330,130 +342,100 @@ Double_t PSimulateMuTransition::GTSpinFlip(const Double_t &time)
// Event (private)
//--------------------------------------------------------------------------
/**
* <p> Generates "muon event": simulate muon spin phase under charge-exchange with
* <p> Generates "muon event": simulate muon spin polarization under charge-exchange with
* 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.
* four "Mu transitions" nu_12, nu_34, nu_23, and nu_14. Use complex polarization
* functions.
* 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
* 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
* muon spin phase by acos(Px(t_i)).
* calculate muon spin precession for t_d, Px(t_i); else calculate spin precession for t_c.
* 3) Determine next ionization time t_i+1; calculate Px(t_i+1) in Muonium. Polarization
* after ionization process is given by Px(t_i+1)*Px(t_i).
* 4) get the next electron capture time, continue until t_d is reached.
*
* <p> For isotropic muonium, TF:
* 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
*
* <p>Calculates Mu0 polarization in x direction during cyclic charge exchange.
* See M. Senba, J.Phys. B23, 1545 (1990), equations (9), (11)
* \param muonString if eq. "Mu+" begin with Mu+ precession
*/
void PSimulateMuTransition::Event(const TString muonString)
Double_t PSimulateMuTransition::Event(const TString muonString)
{
TComplex complexPolX = 1.0;
Double_t muoniumPolX = 1.0; //initial polarization in x direction
Double_t eventTime, eventDiffTime, captureTime, ionizationTime;
// Double_t muonPrecessionFreq, muoniumPrecessionFreq; // MHz
// Double_t rndm, frac1, frac2;
fMuonPrecFreq = fMuonGyroRatio * fBfield;
// charge-exchange loop until muon decay
eventTime = 0.;
eventDiffTime = 0.;
if (fDebugFlag) cout << "Decay time = " << fMuonDecayTime << endl;
//cout << muonString << endl;
// charge-exchange loop until muon decays
while (1) {
if (muonString == "Mu+"){
// Mu+ initial state; get next electron capture time
if (muonString == "Mu+"){// Mu+ initial state; get next electron capture time
captureTime = NextEventTime(fCaptureRate);
eventTime += captureTime;
if (fDebugFlag) cout << "Capture time = " << captureTime << " Phase = " << fMuonPhase << endl;
if (fDebugFlag) cout << "Capture time = " << captureTime << " PolX = " << complexPolX.Re() << endl;
if (eventTime < fMuonDecayTime)
fMuonPhase += PrecessionPhase(captureTime, "Mu+");
complexPolX *= GTFunction(captureTime, "Mu+");
else{ //muon decays; handle precession prior to muon decay
eventDiffTime = fMuonDecayTime - (eventTime - captureTime);
fMuonPhase += PrecessionPhase(eventDiffTime, "Mu+");
complexPolX *= GTFunction(eventDiffTime, "Mu+");
break;
}
// now, we have Mu0; 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 << "Ioniza. time = " << ionizationTime << " Phase = " << fMuonPhase << endl;
if (fDebugFlag) cout << "Ioniza. time = " << ionizationTime << " PolX = " << complexPolX.Re() << endl;
if (eventTime < fMuonDecayTime)
fMuonPhase += PrecessionPhase(ionizationTime, "Mu0");
complexPolX *= GTFunction(ionizationTime, "Mu0");
else{ //muon decays; handle precession prior to muon decay
eventDiffTime = fMuonDecayTime - (eventTime - ionizationTime);
fMuonPhase += PrecessionPhase(eventDiffTime, "Mu0");
complexPolX *= GTFunction(eventDiffTime, "Mu0");
break;
}
}
else{
// Mu0 as initial state; get next ionization time
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;
cout << "Mu Ioniza. time = " << ionizationTime << " PolX = " << complexPolX.Re() << endl;
if (eventTime < fMuonDecayTime)
fMuonPhase += PrecessionPhase(ionizationTime, "Mu0");
complexPolX *= GTFunction(ionizationTime, "Mu0");
else{ //muon decays; handle precession prior to muon decay
eventDiffTime = fMuonDecayTime - (eventTime - ionizationTime);
fMuonPhase += PrecessionPhase(eventDiffTime, "Mu0");
complexPolX *= GTFunction(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 (fDebugFlag) cout << "Capture time = " << captureTime << " PolX = " << complexPolX.Re() << endl;
if (eventTime < fMuonDecayTime)
fMuonPhase += PrecessionPhase(captureTime, "Mu+");
complexPolX *= GTFunction(captureTime, "Mu+");
else{ //muon decays; handle precession prior to muon decay
eventDiffTime = fMuonDecayTime - (eventTime - captureTime);
fMuonPhase += PrecessionPhase(eventDiffTime, "Mu+");
complexPolX *= GTFunction(eventDiffTime, "Mu+");
break;
}
}
}
muoniumPolX = complexPolX.Re();
if (fDebugFlag) cout << " Final PolX = " << muoniumPolX << endl;
if (fDebugFlag) cout << " Final Phase = " << fMuonPhase << endl;
//fMuonPhase = TMath::ACos(TMath::Cos(fMuonPhase))*360./TMath::TwoPi(); //transform back to [0, 180] degree interval
return;
return muoniumPolX;
}

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

@ -60,7 +60,9 @@ class PSimulateMuTransition : public TObject
virtual void SetDecayAsymmetry(Double_t value){ fAsymmetry = value; } //!< muon decay asymmetry
virtual void SetMuFraction(Double_t value){ fMuFraction = value; } //!< Muonium fraction
virtual void SetMuFractionState12(Double_t value){ fMuFractionState12 = value; }
virtual void SetMuFractionState34(Double_t value){ fMuFractionState34 = value; }
virtual void SetMuFractionState23(Double_t value){ fMuFractionState23 = value; }
virtual void SetMuFractionState14(Double_t value){ fMuFractionState14 = value; }
virtual Bool_t IsValid() { return fValid; }
virtual void SetSeed(UInt_t seed);
@ -88,17 +90,18 @@ class PSimulateMuTransition : public TObject
Double_t fMuonPhase; //!< phase of muon spin
Double_t fAsymmetry; //!< muon decay asymmetry
Double_t fMuFraction; //!< total Mu fraction [0,1]
Double_t fMuFractionState12; //!< fraction of Mu in state 12, 34
Double_t fMuFractionState23; //!< fraction of Mu in state 23, 14
Double_t fMuFractionState12; //!< fraction of Mu in state 12
Double_t fMuFractionState34; //!< fraction of Mu in state 34
Double_t fMuFractionState23; //!< fraction of Mu in state 23
Double_t fMuFractionState14; //!< fraction of Mu in state 14
Int_t fNmuons; //!< number of muons to simulate
Bool_t fDebugFlag; //!< debug flag
virtual Double_t NextEventTime(const Double_t &EventRate);
// virtual Double_t PrecessionPhase(const Double_t &time, const Double_t &frequency);
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 PrecessionPhase(const Double_t &time, const TString chargeState);
virtual TComplex GTFunction(const Double_t &time, const TString chargeState); //!< transverse field polarization function of Mu0 or Mu+
virtual Double_t GTSpinFlip(const Double_t &time); //!< transverse field polarization function after spin-flip collisions
virtual void Event(const TString muonString);
virtual Double_t Event(const TString muonString);
ClassDef(PSimulateMuTransition, 0)
};

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

@ -65,10 +65,13 @@ void runMuSimulation()
Double_t Freq23 = 256.245; //Mu freq of the 23 transition
Double_t Freq14 = 356.245; //Mu freq of the 14 transition
Double_t MuFrac = 1.0; //total Mu fraction
Double_t MuFrac12 = 2*0.487; //Mu in states 12 and 34
Double_t MuFrac23 = 2*0.013; //Mu in states 23 and 14
Double_t MuFrac12 = 0.487; //weight of transition 12
Double_t MuFrac34 = 0.487; //weight of transition 34
Double_t MuFrac23 = 0.013; //weight of transition 23
Double_t MuFrac14 = 0.013; //weight of transition 14
Int_t Nmuons = 5e6; //number of muons
Double_t Asym = 0.27; //muon decay asymmetry
Int_t debugFlag = 0; //print debug information on screen
histogramFileName = TString("0");
histogramFileName += runNo;
@ -89,15 +92,17 @@ void runMuSimulation()
simulateMuTransition->SetMuPrecFreq23(Freq23); // MHz
simulateMuTransition->SetMuPrecFreq14(Freq14); // MHz
simulateMuTransition->SetMuFraction(MuFrac); // initial Mu fraction
simulateMuTransition->SetMuFractionState12(MuFrac12); // Mu in states 12, 34
simulateMuTransition->SetMuFractionState23(MuFrac23); // Mu in states 23, 14
simulateMuTransition->SetMuFractionState12(MuFrac12);
simulateMuTransition->SetMuFractionState34(MuFrac34);
simulateMuTransition->SetMuFractionState23(MuFrac23);
simulateMuTransition->SetMuFractionState14(MuFrac14);
simulateMuTransition->SetBfield(B/10000.); // Tesla
simulateMuTransition->SetCaptureRate(capRate); // MHz
simulateMuTransition->SetIonizationRate(ionRate); // MHz
simulateMuTransition->SetSpinFlipRate(spinFlipRate); // MHz
simulateMuTransition->SetNmuons(Nmuons);
simulateMuTransition->SetDecayAsymmetry(Asym);
simulateMuTransition->SetDebugFlag(kFALSE); // to print time and phase during charge-changing cycle
simulateMuTransition->SetDebugFlag(debugFlag); // to print time and phase during charge-changing cycle
// feed run info header
gRunHeader = gROOT->GetRootFolder()->AddFolder("RunHeader", "MuTransition Simulation Header Info");
@ -168,7 +173,9 @@ void runMuSimulation()
header->Set("Simulation/Mu0 Precession frequency 14", Freq14);
header->Set("Simulation/Mu0 Fraction", MuFrac);
header->Set("Simulation/Mu0 Fraction 12", MuFrac12);
header->Set("Simulation/Mu0 Fraction 34", MuFrac34);
header->Set("Simulation/Mu0 Fraction 23", MuFrac23);
header->Set("Simulation/Mu0 Fraction 14", MuFrac14);
header->Set("Simulation/muon Capture Rate", capRate);
header->Set("Simulation/Mu0 Ionization Rate", ionRate);
header->Set("Simulation/Mu0 Spin Flip Rate", spinFlipRate);