Merged muonspin/musrfit into master
This commit is contained in:
Zaher Salman
commit
3b08cad70b
@ -4,6 +4,7 @@
|
||||
|
||||
changes since 0.17.0
|
||||
===================================
|
||||
NEW 2016-04-22 Added the theory function muMinusExpTF for mu minus fits
|
||||
NEW 2016-02-23 It is now possible to export the averaged data/Fourier
|
||||
|
||||
changes since 0.16.0
|
||||
|
||||
@ -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;
|
||||
|
||||
@ -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
|
||||
//--------------------------------------------------------------------------
|
||||
|
||||
@ -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;
|
||||
|
||||
|
||||
@ -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;
|
||||
}
|
||||
|
||||
@ -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)
|
||||
};
|
||||
|
||||
@ -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);
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user